Cannot import X problem. Stiff ODE solver for Oregonator Model

Question

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 the Oregonator model to show that there must be some transition point between the 0 and 3 for the parameter f such that in a particular subset of this interval, oscillations occur.

Forgive my inexperience, I'm new to Python.

Error: ImportError: cannot import name 'radau' from 'scipy.integrate'

In my ignorance, I constructed a 4th Order Runge-Kutta method from scratch. Having found stock prices instead of chemical oscillations, I transitioned to using odeint. This still failed. It was only after this that I discovered the idea of a stiff system, so I have been working on the Radau method.

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]

Expected results should be oscillations like those found in (page 12): https://pdfs.semanticscholar.org/0959/9106a563e9d88ce6442e3bb5b242d5ccbdad.pdf for x and z only. The absence of y is due to a steady-state approximation I have used.

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

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.