无法导入 X 问题。Oregonator 模型的刚性 ODE 求解器

该错误来自尝试从 scipy.integrate 导入 Radau 方法(需要,因为 Oregonator 模型是一个刚性系统)。


我试图对俄勒冈州模型进行数值积分,以表明参数 f 在 0 和 3 之间必须存在某个过渡点,以便在该区间的特定子集中发生振荡。


原谅我的经验不足,我是 Python 新手。


错误:ImportError:无法从“scipy.integrate”导入名称“radau”


在我的无知中,我从头开始构建了一个四阶龙格-库塔方法。在找到股票价格而不是化学波动后,我转而使用 odeint。这仍然失败。直到这之后我才发现了刚性系统的想法,所以我一直在研究 Radau 方法。


import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate. import radau


# Dimensionless parameters

e = 0.04

q = 0.0008

f = 1.0


# Oregonator model

def Oregonator(Y, t):

    return [((Y[0] * (1 - Y[0])  - ((Y[0] - q) * f * Y[1]) // (q + Y[0]))) 

    // e, Y[0] - Y[1]]


# Time span and inital conditions

ts = np.linspace(0, 10, 100)

Y0 = [1, 3]


# Numerical algorithm/method

NumSol = radau(Oregonator, 0, Y0, t_bound=30)

x = NumSol[:,0]

z = NumSol[:,1]

预期的结果应该是类似(第 12 页)中的振荡: https ://pdfs.semanticscholar.org/0959/9106a563e9d88ce6442e3bb5b242d5ccbdad.pdf 仅适用于 x 和 z。y 的缺失是由于我使用了稳态近似值。


慕哥9229398
浏览 102回答 1
1回答

慕村225694

用作求解器类(如或solve_ivp)的单行接口。使用类的正确大写。在 ODE 函数中使用正确的参数顺序(您可以使用in来使用相同的函数)。在您打算使用浮点除法的地方避免整数除法。RK45Radautfirst=Trueodeintimport numpy as npimport matplotlib.pyplot as pltfrom scipy.integrate import solve_ivp# Dimensionless parameterseps = 4e-2q = 8e-4f = 2.0/3# Oregonator modeldef Oregonator(t,Y):    x,z = Y;    return [(x * (1 - x) + (f*(q-x)*z) / (q + x)) / eps, x - z]# Time span and inital conditionsts = np.linspace(0, 10, 100)Y0 = [1, 0.5]# Numerical algorithm/methodNumSol = solve_ivp(Oregonator, [0, 30], Y0, method="Radau")x, z = NumSol.yy = (f*z)/(q+x)t = NumSol.tplt.subplot(221);plt.plot(t,x,'b'); plt.xlabel("t"); plt.ylabel("x");plt.subplot(222);plt.plot(t,y,'r'); plt.xlabel("t"); plt.ylabel("y");plt.subplot(223);plt.plot(t,z,'g'); plt.xlabel("t"); plt.ylabel("z");plt.subplot(224);plt.plot(x,z,'k'); plt.xlabel("x"); plt.ylabel("z");plt.tight_layout(); plt.show()然后这会产生情节表现出周期性振荡。进一步的步骤可能是使用tspan选项或“密集输出”在用户定义的采样点处获取解决方案样本。为获得可靠的结果,请手动设置误差容限。f=0.51262接近从收敛到振荡行为的过渡点。
打开App,查看更多内容
随时随地看视频慕课网APP

相关分类

Python