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

I am trying to learn how to integrate python SymPy with SciPy for numerically solving ordinary differential equations. However, I was kinda lost about how to actually get the SymPy form of a first order system of ODEs into a format that I can process with scipy.integrate.odeint()

我试图学习如何将python SymPy与SciPy集成以数值求解常微分方程。但是，我有点失去了关于如何将一阶ODE系统的SymPy形式转换为我可以用scipy.integrate.odeint（）处理的格式

Note, some people suggested that this is similar to another post, but that is not the case. The other post is here.
Convert sympy expressions to function of numpy arrays

请注意，有些人认为这与另一篇文章类似，但事实并非如此。另一篇文章在这里。将sympy表达式转换为numpy数组的函数

So this other post is a much more complicated case where the user wants to accelerate computation of an ODE using Theano or some other libraries. I am just trying to understand the basic interface between SymPy and SciPy, so this other post is not at all helpful.

所以这篇文章是一个更复杂的案例，用户希望使用Theano或其他一些库来加速ODE的计算。我只是想了解SymPy和SciPy之间的基本接口，所以这篇文章根本没什么用处。

As a toy example, I was using the Lotka-Volterra equation to test the use of SymPy. The equations are:

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

I can solve this in the conventional way with Scipy and it works. Here is the working code.

我可以用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()

Now I was kind lost about how to do this with SymPy. I have a sense of what needs to be done, but was not sure how to proceed. The only example I have found is way too complicated to learn from. Here is what I have.

现在，我对如何使用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)

I am supposed to put these equations into some sort of system form and then use the sympy.lambdify function to return a new function that I can pass to odeint

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

So can anyone fill in the blanks here on how I set up this ode1,ode2 system to process with SymPy.

所以任何人都可以在这里填写关于如何设置这个ode1，ode2系统来处理SymPy的空白。

2 个解决方案

#1

There is rarely need to use Eq objects in SymPy; equations are best represented by expressions, which are understood to be equated to zero in the context of solvers. So, ode1 and ode2 should just be expressions for the right-hand side of the ODE system. Then they can be lambdified together, as follows:

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

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)

The additional step after lambdification is needed because SciPy's solver passes some arrays, which the lambdified function will not know how to unpack. Example:

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

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

returns [-5, -2] which is a Python list rather than a NumPy array, but the ODE solver should handle that. (Or you can return np.array(F_sympy(x, y, a, b, c, d))).)

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

#2

I wrote a module named JiTCODE, which creates ODE integrator objects (with a handling similar to scipy.integrate.ode) from symbolic expressions (SymPy or SymEngine) describing the right-hand side. Under the hood it uses scipy.integrate.ode and scipy.integrate.solve_ivp for the integration.

我写了一个名为JiTCODE的模块，它从描述右侧的符号表达式（SymPy或SymEngine）创建ODE集成器对象（具有类似于scipy.integrate.ode的处理）。在引擎盖下，它使用scipy.integrate.ode和scipy.integrate.solve_ivp进行集成。

The only drawback in your case is that the symbol for the dynamical variables and time is prescribed, so you might have to do a symbolic substitution – which should not be a big issue, however. Below, is your Lotka–Volterra equation as an example, using y(0) instead of x and y(1) instead of y.

在你的情况下唯一的缺点是动态变量和时间的符号是规定的，所以你可能必须做一个象征性的替换 - 但这不应该是一个大问题。下面，以你的Lotka-Volterra方程为例，使用y（0）代替x和y（1）而不是y。

import numpy as np
from sympy.abc import a,b,c,d
from jitcode import y, jitcode

xy0 = [600,400]
argst = [0.4,0.002,0.001,0.7]

lotka_volterra = [
         a*y(0) - b*y(0)*y(1),
        -d*y(1) + c*y(0)*y(1)
    ]

ODE = jitcode( lotka_volterra, control_pars=[a,b,c,d] )
ODE.set_integrator("dopri5")
ODE.set_initial_value(xy0)
ODE.set_parameters(*argst)

times = np.linspace(0, 50, 250)
xy_t = np.vstack(ODE.integrate(time) for time in times)

Note that the main feature of this module is that the right-hand side is compiled for efficiency. Depending on what you need to do, this may be overkill, but it doesn’t harm if it works (you can also disable this and use lambdification as detailed in the other answer).

请注意，此模块的主要功能是编译右侧以提高效率。根据你需要做的事情，这可能是过度的，但如果它有效则不会造成伤害（你也可以禁用它并使用lambdification，详见另一个答案）。