在 python 中求逆矩阵的问题

Question

我正在尝试计算雅可比矩阵的逆矩阵。我使用 sympy 计算了 Jacobian。但是，现在当我尝试 inv 矩阵时，出现了一个我不明白的错误。如果有人能帮助我就太好了

import numpy as py
from numpy.linalg import inv
from sympy import Matrix
import warnings

def f1(y1, y2, y3, y1_old, dt):
    return y1_old + (-0.04*y1 + 10**4*y2*y3)*dt

def f2(y1, y2, y3, y2_old, dt):
    return y2_old + (0.04*y1 - 10**4*y2*y3 - (3.10**7)*(y2**2))*dt

def f3(y1, y2, y3, y3_old, dt):
    return y3_old + ((3.10**7)*(y2**2))*dt

def j(y1,y2,y3):
    a = Matrix(('-0.04*y1 + 10**4*y2*y3','0.04*y1 - 10**4*y2*y3 - (3.10**7)* 
    (y2**2)','(3.10^7)*(y2**2)')).jacobian(('y1', 'y2', 'y3'))
    a = str(a)
    a = eval(a)
    return a

warnings.filterwarnings("ignore", category=DeprecationWarning)

y_old = py.zeros((3,1))
y_old[0] = 1

#Guess values for the implicit variable
y_guess = 2*py.ones((3,1))

#New values - assumed initially
y_new = py.ones((3,1))

F = py.copy(y_new)

start_time = 0
end_time = 10
dt = 0.01
nt = py.arange(start_time,end_time,dt)

error = 9e9
tol = 1e-10
alpha = 0.8
#for i in range(0,len(nt)):

#while error>tol:

jac = j(y_guess[0],y_guess[1],y_guess[2])

F[0] = f1(y_guess[0],y_guess[1],y_guess[2], y_old[0], dt)
F[1] = f2(y_guess[0],y_guess[1],y_guess[2], y_old[1], dt)
F[2] = f3(y_guess[0],y_guess[1],y_guess[2], y_old[2], dt)

j = inv(jac)#error

这是我遇到的错误：

enter image description here

Answer 1

这里有很多错误。首先不要关闭警告：它们存在是有原因的。您应该看到的警告是：

sympy/matrices/repmatrix.py:102: SymPyDeprecationWarning: 

non-Expr objects in a Matrix has been deprecated since SymPy 1.9. Use
list of lists, TableForm or some other data structure instead. See
https://github.com/sympy/sympy/issues/21497 for more info.

  deprecated_since_version="1.9"

看到该警告是因为您将 SymPy 矩阵的某些元素设置为 NumPy 数组，因此您的矩阵看起来像

In [2]: jac
Out[2]: 
⎡-0.04     [20000.0]      [20000.0] ⎤
⎢                                   ⎥
⎢0.04   [-31005.0456444]  [-20000.0]⎥
⎢                                   ⎥
⎣  0    [11005.0456444]       0     ⎦

将来类似的事情会立即引发异常。该警告是一种礼貌的提醒，您应该更改您的代码，否则它将不适用于 SymPy 的未来版本。

在您更好地了解 SymPy 和 NumPy 各自的功能之前，我建议您只坚持使用其中一个，不要尝试将它们混合在一起。

你有像 3.10^7 这样的数字，我认为它应该是 30000000.0 但实际上应该写成 3e7 或 3*10**7.

在 SymPy 中执行此操作的正常方法是：

In [9]: y1, y2, y3 = symbols('y1, y2, y3')

In [10]: a = Matrix([-0.04*y1 + 10**4*y2*y3, 0.04*y1 - 10**4*y2*y3 - 3e7*y2**2, 3e7*y2**2])

In [11]: a
Out[11]: 
⎡        -0.04⋅y₁ + 10000⋅y₂⋅y₃        ⎤
⎢                                      ⎥
⎢                       2              ⎥
⎢0.04⋅y₁ - 30000000.0⋅y₂  - 10000⋅y₂⋅y₃⎥
⎢                                      ⎥
⎢                         2            ⎥
⎣            30000000.0⋅y₂             ⎦

In [12]: jac = a.jacobian([y1, y2, y3])

In [13]: jac
Out[13]: 
⎡-0.04          10000⋅y₃           10000⋅y₂ ⎤
⎢                                           ⎥
⎢0.04   -60000000.0⋅y₂ - 10000⋅y₃  -10000⋅y₂⎥
⎢                                           ⎥
⎣  0          60000000.0⋅y₂            0    ⎦

然后，如果您想要逆矩阵，可以使用 jac.inv() 方法，以便 SymPy 可以计算矩阵的逆矩阵。或者，您可以替换符号的值，然后转换为 NumPy 数组并使用 NumPy 的 inv 函数。

无论哪种方式，您在这里都会遇到问题，因为矩阵对于符号 y1、y2 和 y3:

的所有可能值都是奇异的

In [14]: jac.det()
Out[14]: 0

In [15]: jac.inv()
---------------------------------------------------------------------------
NonInvertibleMatrixError

很容易看出行之和为零，因此矩阵没有满秩。

Answer 2

显示问题的另一种方式，没有 sympy 漂亮的打印

In [5]: jac = j(y_guess[0],y_guess[1],y_guess[2])
In [6]: jac
Out[6]: 
Matrix([
[-0.04,        [20000.0],  [20000.0]],
[ 0.04, [-31005.0456444], [-20000.0]],
[    0,  [11005.0456444],          0]])

numpy inv 适用于数值数组。由于 jac 中列表和标量的混合，结果是一个对象 dtype 数组。

In [7]: np.array(jac)
Out[7]: 
array([[-0.0400000000000000, [20000.0], [20000.0]],
       [0.0400000000000000, [-31005.0456444], [-20000.0]],
       [0, [11005.0456444], 0]], dtype=object)

In [8]: np.linalg.inv(_)
Traceback (most recent call last):
  File "<ipython-input-8-61bbaec16d5b>", line 1, in <module>
    np.linalg.inv(_)
  File "<__array_function__ internals>", line 5, in inv
  File "/usr/local/lib/python3.8/dist-packages/numpy/linalg/linalg.py", line 545, in inv
    ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
TypeError: No loop matching the specified signature and casting was found for ufunc inv

我不知道您是否可以清理 jac 公式，使其不混合使用标量和列表。无论如何，很明显混合使用 sympy 和 numpy 不适合新手，尤其是当他忽略检查中间变量时。

清理多余的[]:

In [14]: M = np.array([[-0.0400000000000000, 20000.0, 20000.0],
    ...:        [0.0400000000000000, -31005.0456444, -20000.0],
    ...:        [0, 11005.0456444, 0]])
    ...: 
    ...: 
In [15]: M
Out[15]: 
array([[-4.00000000e-02,  2.00000000e+04,  2.00000000e+04],
       [ 4.00000000e-02, -3.10050456e+04, -2.00000000e+04],
       [ 0.00000000e+00,  1.10050456e+04,  0.00000000e+00]])

这样做运行，但有不同的错误：

In [16]: np.linalg.inv(M)
....
LinAlgError: Singular matrix

y_guess 是 (3,1) 数组，因此评估 j 时传递 (1,) 数组：

In [17]: y_guess
Out[17]: 
array([[2.],
       [2.],
       [2.]])

相反，如果我们将标量传递给 j，则生成的矩阵没有额外的括号：

In [18]: j(2,2,2)
Matrix([[-0.04*y1 + 10000*y2*y3], [0.04*y1 - 2751.2614111*y2**2 - 10000*y2*y3], [2751.2614111*y2**2]])
Matrix([[-0.0400000000000000, 10000*y3, 10000*y2], [0.0400000000000000, -5502.5228222*y2 - 10000*y3, -10000*y2], [0, 5502.5228222*y2, 0]])
Out[18]: 
Matrix([
[-0.04,          20000,  20000],
[ 0.04, -31005.0456444, -20000],
[    0,  11005.0456444,      0]])

现在可以创建 np.linalg 可以使用的数值数组

In [20]: np.array(j(2,2,2),float)
Out[20]: 
array([[-4.00000000e-02,  2.00000000e+04,  2.00000000e+04],
       [ 4.00000000e-02, -3.10050456e+04, -2.00000000e+04],
       [ 0.00000000e+00,  1.10050456e+04,  0.00000000e+00]])

在 python 中求逆矩阵的问题

Problem with finding the inverse of a matrix in python

python

sympy