+关注
已关注

分类  

暂无分类

标签  

暂无标签

日期归档  

暂无数据

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

发布于2020-06-01 21:56     阅读(740)     评论(0)     点赞(17)     收藏(4)


我正在尝试学习如何将python SymPy与SciPy集成在一起,以数值方式求解常微分方程。但是,我对如何将ODE的一阶系统的SymPy形式实际转换为可以处理的格式感到有些困惑。scipy.integrate.odeint()

请注意,有人建议这类似于另一篇文章,但事实并非如此。另一个帖子在这里。
将sympy表达式转换为numpy数组的函数

因此,这另一篇文章的情况要复杂得多,用户希望使用Theano或其他一些库来加速ODE的计算。我只是想了解SymPy和SciPy之间的基本接口,因此这篇其他文章根本没有帮助。

作为一个玩具示例,我使用Lotka-Volterra方程测试了SymPy的使用。等式为:

\ dot {x} = ax-bxy

\ dot {y} = cxy-dy

我可以使用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做到这一点感到迷lost。我对需要做什么有所了解,但不确定如何进行。我发现的唯一例子太复杂了,无法学习。这就是我所拥有的。

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进行处理的空白


解决方案


很少需要在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)))。)



所属网站分类: 技术文章 > 问答

作者:黑洞官方问答小能手

链接: https://www.pythonheidong.com/blog/article/400996/

来源: python黑洞网

任何形式的转载都请注明出处,如有侵权 一经发现 必将追究其法律责任

17 0
收藏该文
已收藏

评论内容:(最多支持255个字符)