python - SymPy/SciPy : solving a system of ordinary differential equations with different variables

Question

总体而言，我是 SymPy 和 Python 的新手，目前我正在使用 Python 2.7 和 SymPy 0.7.5，目标是:a) 从文本文件中读取微分方程组b) 求解系统

我已经阅读了this question和 this other question ，它们几乎是我要找的，但我还有一个问题:我事先不知道方程组的形式，所以我无法使用 创建相应的函数def 在脚本中，如 this example .整个事情必须在运行时进行管理。

所以，这是我的一些代码片段。假设我有一个包含以下内容的文本文件 system.txt:

dx/dt = 0.0387*x - 0.0005*x*y
dy/dt = 0.0036*x*y - 0.1898*y

我做的是:

# imports
import sympy
import scipy
import re as regex

# define all symbols I am going to use
x = sympy.Symbol('x')
y = sympy.Symbol('y')
t = sympy.Symbol('t')

# read the file
systemOfEquations = []
with open("system.txt", "r") as fp :
   for line in fp :
            pattern = regex.compile(r'.+?\s+=\s+(.+?)$')
            expressionString = regex.search(pattern, line) # first match ends in group(1)   
            systemOfEquations.append( sympy.sympify( expressionString.group(1) ) )

此时，我被 systemOfEquation 列表中的两个符号表达式卡住了。如果我可以从另一个文件读取 ODE 系统的初始条件，为了使用 scipy.integrate.odeint，我必须将系统转换为 Python 可读的函数，例如:

def dX_dt(X, t=0):
return array([ 0.0387*X[0] - 0.0005*X[0]*X[1] ,
              -0.1898*X[1] + 0.0036*X[0]*X[1] ])

有没有一种在运行时创建它的好方法？比如将函数写入另一个文件，然后将新建的文件作为函数导入？ (也许我在这里很愚蠢，但请记住我是 Python 的新手:-D)

我已经看到使用 sympy.utilities.lambdify.lambdify 可以将符号表达式转换为 lambda 函数，但我想知道这是否可以帮助我...lambdify 似乎有效当时只有一种表达方式，而不是系统。

提前感谢您的任何建议:-)

编辑:

经过最少的修改，Warren 的回答完美无缺。我在 listOfSymbols 中有一个所有符号的列表；此外，它们的出现顺序与 odeint 将使用的数据 X 列的顺序相同。所以，我用的函数是

def dX_dt(X, t):
    vals = dict()
    for index, s in enumerate(listOfSymbols) :
            if s != time :
                    vals[s] = X[index]
    vals[time] = t
    return [eq.evalf(subs=vals) for eq in systemOfEquations]

我只是在我的特定问题中为变量“时间”做了一个异常(exception)。再次感谢! :-)

Answer 1

如果您要在读取文件的同一脚本中求解系统(因此 systemOfEquations 可用作全局变量)，并且如果仅变量被使用systemOfEquations 中有 x、y 和可能的 t，您可以在中定义 dX_dt同样的文件是这样的:

def dX_dt(X, t):
    vals = dict(x=X[0], y=X[1], t=t)
    return [eq.evalf(subs=vals) for eq in systemOfEquations]

dX_dt 可以在odeint 中使用。在下面的 ipython session 中，我已经运行了创建 systemOfEquations 并定义 dX_dt 的脚本:

In [31]: odeint(dX_dt, [1,2], np.linspace(0, 1, 5))
Out[31]: 
array([[ 1.        ,  2.        ],
       [ 1.00947534,  1.90904183],
       [ 1.01905178,  1.82223595],
       [ 1.02872997,  1.73939226],
       [ 1.03851059,  1.66032942]]