如何在 ODE 系统的解决方案中使用 sympy 生成的 Jacobi 矩阵？

Question

我有一个由 3 个差分组成的一阶 ODE 系统。 eqn的。我想用scipy.integrate.solve_ivp的BDF方法来解决它。所以我需要计算系统的 jacobi 矩阵（并在 SymPy 的帮助下完成）。

如果我没有误会的话；根据 scipy.integrate.solve_ivp 文档，您必须以 jac(t,u) 的形式引入 jacobien 矩阵，其中 u 应该是 ODE 系统的状态变量。为此，我正确地 lambdify jacobien 矩阵。

我的问题就出在这里。尽管我能够用一些 (t,u) 计算 jac(t,u)，例如 ((1/800),(150,1E-6,3))，但我无法将数组参数发送到我的 jac。当我引入 jac(t,u) 作为 solve_ivp 的参数时，它会给出一条错误消息。那么我应该如何介绍 jac 矩阵呢？还是我的 lambdify 不合适？

这是我的代码。任何帮助我都很感激。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def cvs(t,u):
    u1,u2,u3 = u
    def Qmi(t):
        return t**2
    u1p = Qmi(t)*u3
    u2p = (u1**2)*np.cos(2*np.pi*200*t)
    u3p = (np.sin(2*np.pi*t))*u2**-1
    return [u1p,u2p,u3p]
def jac_func():
    ######### DEFINE THE ODE SYSTEM #########
    import sympy
    sympy.init_printing()
    t = sympy.symbols("t")
    Q_mi = sympy.Function("Q_mi")(t)
    u1 = sympy.Function("u1")(t)
    u2 = sympy.Function("u2")(t)
    u3 = sympy.Function("u3")(t)
    Q_mi = t**2
    u1p = (u3*Q_mi)
    u2p = (u1**2)*sympy.cos(2*sympy.pi*200*t)
    u3p = sympy.sin(2*sympy.pi*5*t)*u2**-1
    ####### CALCULATE THE JACOBIEN ########
    ode_rhs = sympy.Matrix([u1p,u2p,u3p])
    ode_var = sympy.Matrix([u1,u2,u3])
    jac = sympy.Matrix([[ode.diff(var) for var in ode_var]for ode in ode_rhs])
    u = (u1,u2,u3)
    jac_np = sympy.lambdify((t,u),jac,"numpy")
    return jac_np

jac_np = jac_func()
U_0 = [500,20,20]
t = np.linspace(0,100,10000)

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

错误信息：

ValueError                                Traceback (most recent call last)
<ipython-input-1-8b86ffb3a7cf> in <module>()
41 t = np.linspace(0,100,10000)
42 
---> 43 solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

<lambdifygenerated-1> in _lambdifygenerated(t, _Dummy_188)
  1 def _lambdifygenerated(t, _Dummy_188):
  2     [_Dummy_185, _Dummy_186, _Dummy_187] = _Dummy_188
----> 3     return (array([[0, 0, t**2], [2*_Dummy_185*cos(400*pi*t), 0, 0], [0, -sin(10*pi*t)/_Dummy_186**2, 0]]))

ValueError: setting an array element with a sequence.

Answer 1

您遇到问题是因为您按照错误消息的说明进行操作，您正在传递一个数组，其中过程需要一个数字。在

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

您正在尝试将常量矩阵 jac_np(t,U_0) 应用于雅可比行列式参数。但是，此时 t 包含您想要从中输出样本的所有 t 值。 [ array, scalar, scalar ] 列表与 numpy 数组不兼容。

长话短说，删除参数，将雅可比行列式作为可调用函数传递，很可能如您所愿，

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np, t_eval=t)

如何在 ODE 系统的解决方案中使用 sympy 生成的 Jacobi 矩阵？

How to use a sympy generated Jacobi matrix in the solution of an ODE system ?

sympy

scipy

ode