SymPy/SciPy：使用不同变量解决普通微分方程组

Question

我对SymPy和Python不太了解，目前正在使用Python 2.7和SymPy 0.7.5工作，目标是： a）从文本文件中读取一组微分方程 b）解决这个系统

我已经阅读了这个问题和这个其他问题，它们几乎是我所寻找的，但我有一个额外的问题：我事先不知道方程组的形式，因此不能像在这个例子中一样，在脚本内部使用def创建相应的函数。整个过程必须在运行时管理。

因此，这里有我的一些代码片段。假设我有一个文本文件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;此外，它们以数据列X的相同顺序出现，该数据将被odeint使用。因此，我使用的函数是

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]

在我的特定问题中，我只对变量"time"做了一个例外。再次感谢！:-)

- Alberto

---

“我事先不知道方程组的形式。”：在“system.txt”文件中，变量总是被称为“x”和“y”吗？方程组总是只有两个吗？如果不是，那么其他附加变量将被称为什么？ - Warren Weckesser

进一步澄清，我已经知道所有涉及符号的名称（我将从另一个文件中读取它们），而系统中方程的数量可能会有所变化。此外，从这个问题和下面Warren的回答中，我开始看到了一个模式；我可以从一个包含数据的文件中构建一个字典，使用列名作为符号，然后尝试Warren的解决方案。 - Alberto

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 会话中，我已经运行了创建 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]]

- Warren Weckesser