在Python中生成和解决同时的ODE

Posted 2021-04-26

我对Python比较陌生，在编写一段生成然后解决微分方程系统的代码时遇到了一些问题。

我这样做的方法是在函数var（）的列表中重新创建一组变量和系数，（x0，x1，...，xn）和（c0，c1，...，cn）。然后在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求解？

答案

你不想使用evalf，而是想要使用sympy.lambdify来生成一个与SciPy一起使用的回调。您需要创建一个具有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)

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

如果你打开使用外部库来处理外部解算器的函数签名，你可能会对我创作的库感兴趣：pyodesys

另一答案

如果我理解正确，您希望在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的结果是float（不需要evalf）。

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