1

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.

AlphaArgonian
  • 61
  • 1
  • 7
  • If you checked the documentation you would have found that the solver class is called "Radau", with uppercase first letter. // Why do you not use the `solve_ivp` interface? – Lutz Lehmann Oct 19 '19 at 10:23

1 Answers1

0

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

enter image description here

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. enter image description here

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Everything checks out now. It seems I have messed up a lot! I thought you'd have to directly import the Radau method itself, hence the capitalisation issue. Your phrase portraits look better than mine, I'll just add flow arrows for good measure. Thanks for your help. Edit: I have some questions though. 1. With the definition, why CAN you use Y as an input and then use x and z normally in the return line? 2. What do the NumSol. lines mean and why are they needed? – AlphaArgonian Oct 19 '19 at 18:47
  • The Radau class itself contains an initialization and a stepper, similar but not identical to the ode class of the older scipy interface. You can use it directly, but chances are that the solve_ivp framework with event mechanism captures most of the use cases. 1.) The first line unpacks the array Y to the variable list on the left. 2.) solve_ivp returns a "bunch object", a structure that contains the internal step points and their values and additional information on the solution. Note that Matlab's ode45 interpolates 3 additional points per internal step to get rounder curves. – Lutz Lehmann Oct 19 '19 at 22:01