在 Python 中生成和求解同时 ODE

Question

我对 Python 比较陌生，在编写一段生成并求解微分方程组的代码时遇到了一些问题。

我这样做的方法是创建一组变量和系数，(x0, x1, ..., xn) 和 (c0, c1 ,..., cn ) 分别位于具有函数 var() 的列表中。然后在EOM1() 中构造方程。初始条件以及方程组都放在 EOM2() 中，并使用odeint 求解。

目前下面的代码运行，尽管效率不高，我认为这是因为odeint 会在每次交互时运行所有代码（这是我需要修复的其他问题，但不是主要问题！）。

import sympy as sy
from scipy.integrate import odeint

n=2
cn0list = [0.01, 0.05]
xn0list = [0.01, 0.01]

def var():
    xnlist=[]
    cnlist=[]
    for i in range(n+1):
        xnlist.append('x{0}'.format(i))
        cnlist.append('c{0}'.format(i))
    return xnlist, cnlist

def EOM1():
    drdtlist=[]
    for i in range(n):
        cn1=sy.Symbol(var()[1][i])
        xn0=sy.Symbol(var()[0][i])
        xn1=sy.Symbol(var()[0][i+1])
        eom=cn1*xn0*(1.0-xn1)-cn1*xn1-xn1
        drdtlist.append(eom)
    xi=sy.Symbol(var()[0][0])
    xf=sy.Symbol(var()[0][n])
    drdtlist[n-1]=drdtlist[n-1].subs(xf,xi)
    return drdtlist

def EOM2(xn, t, cn):
    x0, x1 = xn
    c0, c1 = cn
    f = EOM1()
    output = []
    for part in f:
        output.append(part.evalf(subs={'x0':x0, 'x1':x1, 'c0':c0, 'c1':c1}))
    return output

abserr = 1.0e-6
relerr = 1.0e-4
stoptime = 10.0
numpoints = 20

t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

wsol = odeint(EOM2, xn0list, t, args=(cn0list,), atol=abserr, rtol=relerr)

我的问题是我很难让 Python 正确处理 Sympy 生成的变量。我用这条线解决了这个问题

output.append(part.evalf(subs={'x0':x0, 'x1':x1, 'c0':c0, 'c1':c1}))

在 EOM2() 中。不幸的是，我不知道如何将这一行概括为从 x0 到 xn 以及从 c0 到 cn 的变量列表。这同样适用于 EOM2(),

中的前一行

    x0, x1 = xn
    c0, c1 = cn

换句话说，我将 n 设置为任意数字，Python 有没有办法像我在上面手动输入的那样解释每个元素？我已经尝试了以下

output.append(part.evalf(subs={'x{0}'.format(j):var(n)[0][j], 'c{0}'.format(j):var(n)[1][j]}))

然而这产生了导致我首先使用 evalf 的错误，

TypeError: can't convert expression to float

有什么方法可以做我想做的事，生成一组 n 方程，然后用odeint 求解？

Answer 1

如果我理解正确，您想在 SymPy 表达式中进行任意数量的替换。可以这样做：

n = 10
syms = sy.symbols('x0:{}'.format(n))   # an array of n symbols
expr = sum(syms)                       # some expression with those symbols
floats = [1/(j+1) for j in range(n)]   # numbers to put in 
expr.subs({symbol: value for symbol, value in zip(syms, floats)})

在这种情况下，subs 的结果是一个浮点数（不需要 evalf）。

请注意，函数symbols 可以通过冒号符号直接为您创建任意数量的符号。不需要循环。

Answer 2

您想研究使用sympy.lambdify 来生成用于SciPy 的回调，而不是使用evalf。您将需要创建一个预期签名为odeint 的函数，例如：

y, params = sym.symbols('y:3'), sym.symbols('kf kb')
ydot = rhs(y, p=params)
f = sym.lambdify((y, t) + params, ydot)
yout = odeint(f, y0, tout, param_values)

我们在 SciPy 2017 会议上提供了关于（除其他外）如何使用 lambdify 和 odeint 的教程，材料可在此处获得：http://www.sympy.org/scipy-2017-codegen-tutorial/

如果您愿意使用外部库来处理外部求解器的函数签名，您可能会对我编写的库感兴趣：pyodesys