无法导入 X 问题。 Oregonator 模型的刚性 ODE 求解器
Cannot import X problem. Stiff ODE solver for Oregonator Model
错误来自于尝试从 scipy.integrate 导入 Radau 方法(需要,因为 Oregonator 模型是一个刚性系统)。
我正在尝试对 Oregonator 模型进行数值积分,以表明参数 f 的 0 和 3 之间必须存在某个过渡点,以便在此区间的特定子集中发生振荡。
请原谅我的经验不足,我是 Python 的新手。
错误:导入错误:无法从 '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 是由于我使用了稳态近似。
使用 solve_ivp
作为求解器 class 的 solve_ivp
界面,例如 RK45
或 Radau
。使用 class 的正确大写。在 ODE 函数中使用正确的参数顺序(您可以在 odeint
中使用 tfirst=True
以在所有地方使用相同的函数)。避免在打算使用浮点除法的地方进行整数除法。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# Dimensionless parameters
eps = 4e-2
q = 8e-4
f = 2.0/3
# Oregonator model
def Oregonator(t,Y):
x,z = Y;
return [(x * (1 - x) + (f*(q-x)*z) / (q + x)) / eps, x - z]
# Time span and inital conditions
ts = np.linspace(0, 10, 100)
Y0 = [1, 0.5]
# Numerical algorithm/method
NumSol = solve_ivp(Oregonator, [0, 30], Y0, method="Radau")
x, z = NumSol.y
y = (f*z)/(q+x)
t = NumSol.t
plt.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
选项或 "dense output" 在 user-defined 个采样点获取解决方案样本。为获得可靠的结果,请手动设置误差容限。
f=0.51262
接近从收敛到振荡行为的过渡点。
错误来自于尝试从 scipy.integrate 导入 Radau 方法(需要,因为 Oregonator 模型是一个刚性系统)。
我正在尝试对 Oregonator 模型进行数值积分,以表明参数 f 的 0 和 3 之间必须存在某个过渡点,以便在此区间的特定子集中发生振荡。
请原谅我的经验不足,我是 Python 的新手。
错误:导入错误:无法从 '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 是由于我使用了稳态近似。
使用 solve_ivp
作为求解器 class 的 solve_ivp
界面,例如 RK45
或 Radau
。使用 class 的正确大写。在 ODE 函数中使用正确的参数顺序(您可以在 odeint
中使用 tfirst=True
以在所有地方使用相同的函数)。避免在打算使用浮点除法的地方进行整数除法。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# Dimensionless parameters
eps = 4e-2
q = 8e-4
f = 2.0/3
# Oregonator model
def Oregonator(t,Y):
x,z = Y;
return [(x * (1 - x) + (f*(q-x)*z) / (q + x)) / eps, x - z]
# Time span and inital conditions
ts = np.linspace(0, 10, 100)
Y0 = [1, 0.5]
# Numerical algorithm/method
NumSol = solve_ivp(Oregonator, [0, 30], Y0, method="Radau")
x, z = NumSol.y
y = (f*z)/(q+x)
t = NumSol.t
plt.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
选项或 "dense output" 在 user-defined 个采样点获取解决方案样本。为获得可靠的结果,请手动设置误差容限。
f=0.51262
接近从收敛到振荡行为的过渡点。