SymPy：从采用 numpy 数组的对角矩阵创建一个 numpy 函数答案

【问题标题】：SymPy : creating a numpy function from diagonal matrix that takes a numpy arraySymPy：从采用 numpy 数组的对角矩阵创建一个 numpy 函数
【发布时间】：2017-06-07 01:58:21
【问题描述】：

基于我找到的here 示例，我正在尝试从使用sumpy.diag 创建的对角矩阵创建一个函数

myM = Matrix([
[x1, 4, 4],
[4, x2, 4],
[4, 4, x3]])

这是使用此例程创建的，例如：

import sympy as sp
import numpy as np

x1 = sp.Symbol('x1')
x2 = sp.Symbol('x2')
x3 = sp.Symbol('x3')
X = sp.Matrix([x1, x2, x3])

myM = 4 * sp.ones(3, 3)
sp.diag(*X) + myM - sp.diag(*np.diag(myM))

现在我想创建一个函数，使用ufuncify 中的lambdify，它将numpy.array 或长度为3（如np.array([0.1,0.2,0.3])）作为输入，并根据@ 将输出作为矩阵给出987654331@

myM = Matrix([
[0.1, 4, 4],
[4, 0.2, 4],
[4, 4, 0.3]])

最终我需要使用这种方法来象征性地创建一个雅可比矩阵：并且由于函数形式在计算过程中可能会发生变化，因此以符号方式计算雅可比将非常有用。

【问题讨论】：

标签： python numpy matrix sympy symbolic-math

【解决方案1】：

从数字向量创建数字 3 x 3 矩阵并不是真正的 SymPy 事情，因为不涉及任何符号。考虑以下情况，其中参数 d 是一个包含对角线元素的数组。

def mat(d):
    return np.diag(d-4) + 4

上述函数返回一个二维 NumPy 数组。要改为返回 SymPy 矩阵，请使用

def mat(d):
    return sp.Matrix(np.diag(d-4) + 4)

当 d 的值非常小时，减法后加法可能会导致精度损失：例如，(1e-20 - 4) + 4 的计算结果为零。更安全的选择是

def mat(d):
    diagmat = np.diag(d) 
    return diagmat + np.fromfunction(lambda i, j: (i != j)*4, diagmat.shape)

【讨论】：

不完全是我的想法，但我会使用它，所以谢谢

【解决方案2】：

您可以 .subs() 将值浮动到相应的符号中：

import sympy as sp
import numpy as np

x1 = sp.Symbol('x1')
x2 = sp.Symbol('x2')
x3 = sp.Symbol('x3')
X = sp.Matrix([x1, x2, x3])

myM = 4 * sp.ones(3, 3)
smyM=sp.diag(*X) + myM - sp.diag(*np.diag(myM))

fcoefs = [(a, f) for a, f in (zip([x1, x2, x3], np.array([0.1,0.2,0.3])))]

fmyM = smyM.subs(fcoefs)

smyM
Out[105]: 
Matrix([
[x1,  4,  4],
[ 4, x2,  4],
[ 4,  4, x3]])

fmyM
Out[106]: 
Matrix([
[0.1,   4,   4],
[  4, 0.2,   4],
[  4,   4, 0.3]])

似乎是一个很好的sympy.matrices.dense.MutableDenseMatrix矩阵后：

fmyM @ myM
Out[107]: 
Matrix([
[32.4, 32.4, 32.4],
[32.8, 32.8, 32.8],
[33.2, 33.2, 33.2]])

可能需要转换为 np.array 才能完全用于 numpy

下面是我的一些代码，展示了我使用的更多模式：

def ysolv(coeffs):
    x,y,a,b,c,d,e = symbols('x y a b c d e')
    ellipse = a*y**2 + b*x*y + c*x + d*y + e - x**2
    y_sols = solve(ellipse, y)
    print(*y_sols, sep='\n')

    num_coefs = [(a, f) for a, f in (zip([a,b,c,d,e], coeffs))]
    y_solsf0 = y_sols[0].subs(num_coefs)
    y_solsf1 = y_sols[1].subs(num_coefs)

    f0 = lambdify([x], y_solsf0)
    f1 = lambdify([x], y_solsf1)
    return f0, f1

f0, f1 = ysolv(t[0])

y0 = [f0(x) for x in xs]
y1 = [f1(x) for x in xs]
...

来自：https://stackoverflow.com/a/41232062/6876009（是的，我的“feeloutXrange”有一个非常糟糕的 hack，必须显示）

【讨论】：