如何从文本文件中读取微分方程组并使用 scipy.odeint 求解该方程组？

我有一个大型（>2000 个方程）ODE 系统，我想用 python scipy 的 odeint 来求解。

我有三个问题想要解决（也许我需要问 3 个不同的问题？）。 为了简单起见，我将在这里用一个玩具模型来解释它们，但请记住我的系统很大。 假设我有以下 ODE 系统：

dS/dt = -beta*S
dI/dt = beta*S - gamma*I
dR/dt = gamma*I

β = cpI

其中 c、p 和 gamma 是我要传递给 odeint 的参数。

odeint 期待这样的文件：

def myODEs(y, t, params):
    c,p, gamma = params
    beta = c*p
    S = y[0]
    I = y[1]
    R = y[2]
    dydt = [-beta*S*I,
           beta*S*I - gamma*I,
           - gamma*I]  
    return dydt

然后可以像这样传递给 odeint ：

myoutput = odeint(myODEs, [1000, 1, 0], np.linspace(0, 100, 50), args = ([c,p,gamma], ))

我在 Mathematica 中生成了一个文本文件，例如 myOdes.txt，其中文件的每一行对应于我的 ODE 系统的 RHS，所以它看起来像这样

#myODEs.txt

-beta*S*I
beta*S*I - gamma*I
- gamma*I

我的文本文件看起来与 odeint 所期望的类似，但我还没有完全实现。 我有三个主要问题：

1. 如何传递我的文本文件以便 odeint 理解这是我系统的 RHS？
2. 如何以智能的方式（即系统的方式）定义我的变量？由于它们的数量超过 2000 个，我无法手动定义它们。理想情况下，我会将它们定义在一个单独的文件中并读取它。
3. 我如何也将参数（有很多）作为文本文件传递？

I read 这个问题这与我的问题 1 和 2 很接近，并尝试复制它（我直接输入参数值，这样我就不必担心上面的第 3 点）：

    systemOfEquations = []
    with open("myODEs.txt", "r") as fp :
        for line in fp :
            systemOfEquations.append(line)

    def dX_dt(X, t):
        vals = dict(S=X[0], I=X[1], R=X[2], t=t)
        return [eq for eq in systemOfEquations]

    out = odeint(dX_dt, [1000,1,0], np.linspace(0, 1, 5))

但我收到错误：

    odepack.error: Result from function call is not a proper array of          floats.
    ValueError: could not convert string to float: -((12*0.01/1000)*I*S),

编辑：我将代码修改为：

    systemOfEquations = []
    with open("SIREquationsMathematica2.txt", "r") as fp :
        for line in fp :
               pattern = regex.compile(r'.+?\s+=\s+(.+?)$')
               expressionString = regex.search(pattern, line) 
               systemOfEquations.append( sympy.sympify( expressionString) )


    def dX_dt(X, t):
        vals = dict(S=X[0], I=X[1], R=X[2], t=t)
        return [eq for eq in systemOfEquations]

    out = odeint(dX_dt, [1000,1,0], np.linspace(0, 100, 50), )

这是可行的（我不太明白 for 循环的前两行在做什么）。但是，我希望更自动地定义变量的过程，并且我仍然不知道如何使用此解决方案并在文本文件中传递参数。同样，如何在 dX_dt 函数内定义参数（取决于变量）？

提前致谢！

这不是完整的答案，而是一些观察/问题，但它们太长，无法发表评论。

dX_dt被多次调用odeint与一维数组y和元组t。你提供t通过args范围。y是由生成的odeint并且随着每一步的变化而变化。dX_dt应该精简，以便运行速度快。

通常是这样的表达[eq for eq in systemOfEquations]可以简化为systemOfEquations. [eq for eq...]没有做任何有意义的事情。但可能有一些关于systemOfEquations这需要它。

我建议你打印出来systemOfEquations（对于这个小型 3 行案例），为了您和我们的利益。您正在使用sympy将文件中的字符串转换为方程。我们需要看看它会产生什么。

注意myODEs是一个函数，而不是一个文件。它可以从模块导入，模块当然是文件。

重点是vals = dict(S=X[0], I=X[1], R=X[2], t=t)是产生一个字典sympy表达式可以一起使用。更直接（而且我认为更快）dX_dt函数看起来像：

def myODEs(y, t, params):
    c,p, gamma = params
    beta = c*p
    dydt = [-beta*y[0]*y[1],
           beta*y[0]*y[1] - gamma*y[1],
           - gamma*y[1]]  
    return dydt

我怀疑dX_dt运行 sympy 生成的表达式将比像这样的“硬编码”慢很多。

我要添加sympy标签，因为正如所写，这是将文本文件转换为函数的关键odeint可以使用。

我倾向于将方程的变异性放在t参数，而不是 sympy 表达式的列表。

即替换：

    dydt = [-beta*y[0]*y[1],
           beta*y[0]*y[1] - gamma*y[1],
           - gamma*y[1]]  

与类似的东西

    arg12=np.array([-beta, beta, 0])
    arg1 = np.array([0, -gamma, -gamma])
    arg0 = np.array([0,0,0])
    dydt = arg12*y[0]*y[1] + arg1*y[1] + arg0*y[0]

一旦这是正确的，那么argxx定义可以移到外部dX_dt，并通过args. Now dX_dt只是一个简单而快速的计算。

这整个sympy这种方法可能效果很好，但恐怕在实践中会很慢。但有更多的人sympy经历也许会有其他的感悟。