Hopf Bifurcation Plot

Question

I am attempting to code a bifurcation diagram to illustrate the values of f for which the Oregonator model yields oscillatory behaviour. I get the "setting an array element with a sequence" error at the solve_ivp line. I suspect it has something to do the time span but I am not sure. I should get a Hopf bifurcation, i.e. a bullet-like cone in the region of oscillations.

Here is the code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Dimensionless constant parameters
eps = 0.04
a = 0.0008

# Dimensionless varying parameter - will reveal limit cycle region
f = np.linspace(-5,5,250)

# Oregonator model
def Oregonator(t, Y):
    x,z = Y;
    return [(x * (1 - x) + ((a - x) * f * z) / (a + x)) / eps, x - z]

# Time span, initial conditions

ts = np.linspace(-5, 5, 250)
Y0 = [1, 0.5]


# Numerical algorithm/method
NumSol = solve_ivp(Oregonator, [0, 30], Y0, method="Radau")
t = NumSol.t
x,z = NumSol.y

# Plot
fig = plt.figure()
ax = fig.gca(projection='3d')


ax.plot(f, x, z, 'm')
ax.set_xlabel('$f$', fontsize=10)
ax.set_ylabel('$x^*$', fontsize=10)
ax.set_zlabel('$z^*$', fontsize=10)

ax.axis([-5, 5, -5, 5])
plt.grid()
plt.show()

score 0 · Answer 1 · answered Feb 20 '20 at 22:08

You can not do this kind of simultaneous computation. (That is, you can, but it requires explicit coding and is still ill-advised as the step size selection may greatly vary over the range of f values.)

You should compute the solution for each of the f values separately, and then plot them. To build a list of all solutions first, it is useful to assemble the construction for a single solution and its last point into a separate function so that the list construction can be done via list processing.

# Dimensionless constant parameters
eps = 0.04
a = 0.0008

def limit(f):
    # Oregonator model
    def Oregonator(t, Y):
        x,z = Y;
        return [(x * (1 - x) + ((a - x) * f * z) / (a + x)) / eps, x - z]

    # Time span, initial conditions

    ts = np.linspace(-5, 5, 250)
    Y0 = [1, 0.5]


    # Numerical algorithm/method
    NumSol = solve_ivp(Oregonator, [0, 30], Y0, method="Radau")
    t = NumSol.t
    x,z = NumSol.y
    return x[-1],z[-1]

# Dimensionless varying parameter - will reveal limit cycle region
f = np.linspace(-5,5,250)

x,z = np.array([ limit(ff) for ff in f ]).T

# Plot
fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot(f, x, z, 'mo', ms=2)

ax.set_xlabel('$f$', fontsize=10)
ax.set_ylabel('$x^*$', fontsize=10)
ax.set_zlabel('$z^*$', fontsize=10)

ax.axis([-5, 5, -5, 5])
plt.grid()
plt.show()

This results in a plot

Hopf Bifurcation Plot

1 Answers1