Cannot Import X Problem. Stiff Ode Solver For Oregonator Model

The error comes from attempting to import the Radau method from scipy.integrate (needed because the Oregonator model is a stiff system). I am attempting to numerically integrate t

Solution 1:

Use solve_ivp as a one-line interface to the solver classes like RK45 or Radau. Use the correct uppercasing of the class. Use the correct argument order in the ODE function (you can use tfirst=True in odeint to use the same function all over). Avoid integer division where you intent to use floating point division.

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 modeldefOregonator(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()

This then produces the plot

exhibiting periodic oscillations.

Further steps could be to use the tspan option or "dense output" to get the solution samples at user-defined sampling points. For reliable results set the error tolerances manually.

f=0.51262 is close to the transition point from convergent to oscillating behavior.