使用 SymPy 表达式和 SciPy 求解器求解一阶 ODE 系统
Posted
技术标签:
【中文标题】使用 SymPy 表达式和 SciPy 求解器求解一阶 ODE 系统【英文标题】:Solving a first order system of ODEs using SymPy expressions and SciPy solver 【发布时间】:2018-09-20 18:06:35 【问题描述】:我正在尝试学习如何将 python SymPy 与 SciPy 集成以数值求解常微分方程。但是,对于如何将一阶 ODE 系统的 SymPy 形式实际转换为我可以使用 scipy.integrate.odeint()
处理的格式,我有点迷茫。
注意,有些人认为这与另一篇文章类似,但事实并非如此。另一个帖子在这里。Convert sympy expressions to function of numpy arrays
所以这篇文章是一个更复杂的案例,用户希望使用 Theano 或其他一些库来加速 ODE 的计算。我只是想了解 SymPy 和 SciPy 之间的基本接口,所以这篇其他帖子根本没有帮助。
作为一个玩具示例,我使用 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
的新函数
那么任何人都可以在这里填写有关我如何设置此ode1,ode2
系统以使用 SymPy 进行处理的空白。
【问题讨论】:
this question 是您的副本吗? @Warren 感谢您找到这个额外的问题。我认为您发布的问题比我正在做的要高级一些。我试图解决一个相对简单的 ODE 来学习 API。在您发布的问题中,用户正在尝试使用 GPU 等来加速数值解的计算。所以我可以使用那个更复杂的答案中的东西,直到我明白如何先做这个更简单的答案。 Convert sympy expressions to function of numpy arrays的可能重复 【参考方案1】:在 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)
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)))
。)
【讨论】:
太好了,我会试一试。是的,也感谢scipy
的提示。我还没有找到关于sympy
到scipy
连接的大量示例。我发现的许多示例都太复杂了,无法开始学习——尽管在我更习惯了这些示例之后它们会很有用。感谢您的帮助。【参考方案2】:
我写了a module named JiTCODE,它从描述右侧的符号表达式(SymPy 或 SymEngine)创建 ODE 积分器对象(处理类似于 scipy.integrate.ode
)。在后台,它使用scipy.integrate.ode
和scipy.integrate.solve_ivp
进行集成。
在您的情况下,唯一的缺点是动态变量和时间的符号是规定的,因此您可能必须进行符号替换——但这应该不是一个大问题。
下面以您的 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)
请注意,该模块的主要特点是右侧为提高效率而编译。根据您需要做的事情,这可能有点过头了,但如果它有效也无妨(您也可以禁用它并使用 the other answer 中详述的lambdification)。
【讨论】:
哦@Wrzlprmft 非常感谢您的回答。太好了,我可以试试 JiTCODE。实际上,我还必须对延迟微分方程做一些工作,所以你的 JiTCDDE 模块看起来也很棒。如果我可以为标准 ODES 和延迟微分方程提供一致的接口,那么这将是一个巨大的帮助。再次感谢。 嘿@Wrzlprmft 只是一个我不清楚的问题。用户是否需要为集成例程指定雅可比矩阵,还是 JiTCODE 会自动生成?我看到类上有一个wants_jacobian
参数和generate_jac_sym
方法。但我不清楚是否需要指定雅可比行列式或者是否自动处理。谢谢。
@krishnab:雅可比将自动计算——如果积分器可以使用它。你提到的东西只是为了微调而存在的。您通常不需要为此烦恼。
哇,太好了。非常感谢这个包:)。以上是关于使用 SymPy 表达式和 SciPy 求解器求解一阶 ODE 系统的主要内容,如果未能解决你的问题,请参考以下文章