python - 使用 SymPy 表达式和 SciPy 求解器求解 ODE 的一阶系统

Question

我正在尝试学习如何将 python SymPy 与 SciPy 集成以数值求解常微分方程。但是，我有点迷失了如何将 ODE 的一阶系统的 SymPy 形式实际转换为我可以使用 scipy.integrate.odeint()

处理的格式

请注意，有些人认为这与另一篇文章相似，但事实并非如此。另一篇文章在这里。
Convert sympy expressions to function of numpy arrays

所以另一篇文章是一个更复杂的案例，用户想要使用 Theano 或其他一些库来加速 ODE 的计算。我只是想了解 SymPy 和 SciPy 之间的基本接口(interface)，所以这篇文章一点帮助也没有。

作为玩具示例，我使用 Lotka-Volterra 方程来测试 SymPy 的使用。方程式是:

我可以用 Scipy 以常规方式解决这个问题，而且它有效。这是工作代码。

import numpy as np
import scipy
from scipy.integrate import odeint, ode, solve_ivp
import sympy
import matplotlib.pyplot as plt
sympy.init_printing()

def F_direct(X, t, args):
    F = np.zeros(2)
    a, b, c, d = args
    x,y = X
    F[0] = a*x - b*x*y
    F[1] = c*x*y- d*y
    return F

argst = [0.4,0.002,0.001,0.7]
xy0 = [600, 400]
t = np.linspace(0, 50, 250)
xy_t, infodict = odeint(F_direct, xy0, t, args=(argst,), full_output=True)

plt.plot(t, xy_t[:,1], 'o', t, xy_t[:,0])
plt.grid(True)
plt.xlabel('x'); plt.ylabel('y')
plt.legend(('Numerical', 'Exact'), loc=0)
plt.show()

现在我对如何使用 SymPy 做到这一点感到迷茫。我知道需要做什么，但不确定如何进行。我发现的唯一示例太复杂而无法学习。这是我的。

x, y, a, b, c, d = sympy.symbols('x y a b c d')
t = np.linspace(0, 50, 250)
ode1 = sympy.Eq(a*x - b*x*y)
ode2 = sympy.Eq(c*x*y - d*y)

我应该将这些等式放入某种系统形式，然后使用 sympy.lambdify 函数返回一个我可以传递给 odeint 的新函数/p>

关于我如何设置此 ode1,ode2 系统以使用 SymPy 进行处理，任何人都可以在这里填写空白。

Answer 1

在 SymPy 中很少需要使用 Eq 对象；方程最好用表达式表示，在求解器的上下文中，表达式被理解为等于零。因此，ode1 和 ode2 应该只是 ODE 系统右侧的表达式。然后它们可以一起被 lambdified，如下所示:

ode1 = a*x - b*x*y
ode2 = c*x*y - d*y
F_sympy = sympy.lambdify((x, y, a, b, c, d), [ode1, ode2])
def F_from_sympy(X, t, args):
    a, b, c, d = args
    x, y = X    
    return F_sympy(x, y, a, b, c, d)

lambdification 之后需要额外的步骤，因为 SciPy 的求解器传递了一些数组，lambdified 函数不知道如何解包。示例:

F_from_sympy([1, 2], np.linspace(0, 1, 100), (3, 4, 5, 6))

返回 [-5, -2]，这是一个 Python 列表而不是 NumPy 数组，但 ODE 求解器应该处理它。 (或者您可以返回 np.array(F_sympy(x, y, a, b, c, d)))。)