代码之家 › 专栏 › 技术社区 › Tbone Willsone

Python中同时生成和求解ODE

sympy scipy python

Tbone Willsone · 技术社区 · 6 年前

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

我的方法是创建一组变量和系数, (x0,x1,…,xn) 和 (c0,c1,…,cn) 代表性地,在具有函数的列表中 var() . 然后在中构造方程 EOM1() . 初始条件以及方程组都被放在 EOM2() 并使用 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正确处理Symphy生成的变量。我用绳子绕过去了

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

在里面 EOM2() . 不幸的是,我不知道如何将这一行推广到 x0 到 xn公司 ,和来自 c0 到 中国大陆 . 这同样适用于中的前一行 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 然后用 奥丁 ?

2 回复 | 直到 6 年前

Bjoern Dahlgren 6 年前

而不是使用 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)

我们提供了一个关于如何使用 lambdify 具有 奥丁 在2017年SciPy会议上,可在以下位置获取材料: http://www.sympy.org/scipy-2017-codegen-tutorial/

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

user6655984 user6655984 6 年前

如果我理解正确,您希望在一个symphy表达式中进行任意数量的替换。这是如何做到的:

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)})

在这种情况下,SUB的结果是一个浮动(不需要evalf)。

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