当方程中的变量自动生成时，如何用 Scipy 求解 ODE 系统

Question

我正在以方程列表的形式生成 ODE 系统，其中变量为 sympy.Symbol()，例如 [3*sympy.Symbol('x')**3+sympy.Symbol('y') , sympy.Symbol('x')-sympy.Symbol('y')**4]。

所以用这个例子，我可以通过代码解决这个系统

import numpy as np
from sympy import Symbol  
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def kin(t, z):
    x, y = z
    eqs = [3*x**3+y  , x-y**4]
    return eqs

y0 = [0, 10]
t = np.linspace(0, 10000, 10000)
sol = solve_ivp(kin, [0, 10000], y0, dense_output=True)
print(sol.sol(t).T)
plt.plot(t, sol.sol(t).T)
plt.show()

它会给出这个结果https://imgur.com/3KvOkrU

但由于 x 和 y 的变量实际上是由不同的代码生成的，我想我可以做到

import numpy as np
from sympy import Symbol  
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

dic = {'x': Symbol('x'), 'y': Symbol('y')}
vars = ('x', 'y')
#eqs =[3x^3+y , x-y^4]
eqs = [3*dic[vars[0]]**3+dic[vars[1]], dic[vars[0]]-dic[vars[1]]**4]

def kin(t, z):
    dics = {}
    for v in range(len(vars)):
        dics[vars[v]] = z[v]

    for e in range(len(eqs)):
        eqs[e] = eqs[e].subs(dics)
    return eqs

y0 = [0, 10]
t = np.linspace(0, 10000, 10000)
sol = solve_ivp(kin, [0, 10000], y0, dense_output=True)
print(sol.sol(t).T)
plt.plot(t, sol.sol(t).T)
plt.show()

它只是用 z[0] 和 z[1] 替换了第一部分的 x 和 y 的值，但是这段代码的输出是完全不同的 https://imgur.com/ovi3mk8

我看不出问题出在哪里，如果有人代码指出问题出在哪里或帮助我找到一种方法来评估这些自动生成的 ODE 系统，我将不胜感激。

Answer 1

您得到该结果的原因是因为您在每次迭代时都在 kin 中覆盖了 eqs。此外，正如 Lutz 所提到的，您应该使用 lambdify 它将符号表达式转换为数值函数，以便它们可以被 Numpy 计算（更快）。

import numpy as np
from sympy import Symbol  
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

dic = {'x': Symbol('x'), 'y': Symbol('y')}
vars = ('x', 'y')
eqs = [3*dic[vars[0]]**3+dic[vars[1]], dic[vars[0]]-dic[vars[1]]**4]
# convert symbolic expressions to numerical functions
eqs = [lambdify(list(dic.values()), e) for e in eqs]

def kin(t, z):
    # evaluate each numerical function with the current values
    return [e(*z) for e in eqs]

y0 = [0, 10]
t = np.linspace(0, 10000, 10000)
sol = solve_ivp(kin, [0, 10000], y0, dense_output=True)
plt.plot(t, sol.sol(t).T)
plt.show()

这将为您提供您发布的第一张图片。

当方程中的变量自动生成时，如何用 Scipy 求解 ODE 系统

How to solve a system of ODEs with Scipy when the variables in the equations are autogenerated

python

numpy

sympy

scipy

ode