The Azimuth Project
Stochastic Hopf bifurcation in Python

Contents

Idea

Sometimes, it’s convenient to have a self-contained implementation of an idea which one can then carry around. Here, we present a stochastic Hopf bifurcation model in the Python (also see Python) language, using the Scipy and matplotlib/pylab libraries, which are useful for scientific computations and graphical displays.

Code

#!/usr/bin/python
# Implements the stochastic Hopf-bifurcation model illustrated in
# Week 308 of John Baez's This Week's Finds
# http://johncarlosbaez.wordpress.com/2010/12/24/this-weeks-finds-week-308/

import sys, random             # Python ships with these
import scipy, pylab     # these are extra

def hopf(x, y, beta, dt, lamb):
    """
    Compute the change in coordinates given the current position,
    the parameters which govern the stochastic Hopf dynamics and the 
    Euler integration timestep.
    """
    rsquared = x * x + y * y
    dx = (-y + beta * x - x * rsquared) * dt
    dy = (x + beta * y - y * rsquared) * dt
    if lamb != 0.0:
        sigma = sqrt(dt)
        dx += lamb * gauss(mu = 0.0, sigma = sigma)
        dy += lamb * gauss(mu = 0.0, sigma = sigma)
    return dx, dy

def split_param(param_string):
    """
    Split an item taken from the command line into a variable name and
    a value.
    """
    split = param_string.split("=")
    try:
        param_name = split[0].lower()
        param_value = float(split[1])
    except:
        param_name = param_string.lower()
        param_value = 0.0
    return (param_name, param_value)

def parse_command_line(command_line):
    """
    Use the command line arguments to specify the values of the
    simulation parameters; assign reasonable default values to
    all parameters left unspecified.
    """
    pairs = dict([split_param(item)
                  for item in command_line])
    try:
        beta = pairs["beta"]        # x-y coupling coefficient
    except:
        beta = 0.1
    try:
        lamb = pairs["lambda"]      # noise strength
    except:
        lamb = 0.0
    try:
        dt = pairs["dt"]            # timestep
    except:
        dt = 0.001
    try:
        x_init = pairs["x"]         # initial value for x
    except:
        x_init = 0.0
    try:
        y_init = pairs["y"]         # initial value for y
    except:
        y_init = 0.0
    try:
        T = pairs["T"]              # duration of the simulation
    except:
        T = int(1e4)
    return beta, lamb, dt, T, x_init, y_init

# pull functions out of libraries for our later convenience
gauss = random.gauss
sqrt = scipy.sqrt

# get command-line parameters
beta, lamb, dt, T, x_init, y_init  = parse_command_line(sys.argv[1:])

# initialize the arrays of coordinate values
aX = scipy.zeros(T)
aY = scipy.zeros(T)
aX[0] = x_init
aY[0] = y_init

# MAIN LOOP
for time_step in xrange(1, T):
    x = aX[time_step - 1]
    y = aY[time_step -1]
    dx, dy = hopf(x, y, beta, dt, lamb)
    aX[time_step] = x + dx
    aY[time_step] = y + dy

# display output
# Fig. 1: X vs. Y
pylab.plot(aX, aY)
pylab.xlabel("X")
pylab.ylabel("Y")
pylab.title("beta = " + str(beta)
            + ", lambda = " + str(lamb))
# Fig. 2: X as a function of time
pylab.figure()
pylab.plot(aX)
pylab.xlabel("Time")
pylab.ylabel("X")
pylab.title("beta = " + str(beta)
            + ", lambda = " + str(lamb))
pylab.show()

Operation

This program allows input from the command line, of the form variablename=value. For example, one might type ./stochastic-hopf.py beta=-0.25 lambda=0.1 to get plots resembling those shown below. (Python actually provides a library module for fairly sophisticated handling of command-line arguments; however, the do-it-yourself approach taken here may be more illustrative for those getting started with the language.) A typical Figure 1, showing YY plotted against XX, with β=0.25\beta = -0.25 and λ=0.1\lambda = 0.1, would be the following:

And this is Figure 2, showing XX as a function of time tt, for the same simulation run:

As the Python pseudo-random-number generator is by default seeded with the computer’s clock, one should not expect the plots made by this code to appear exactly the same twice.

Exercises

When learning to program as a novice, picking up a new language as an experienced programmer or studying a new scientific concepts, exercises are useful. (Plus, the Azimuth Project wiki looks more helpful for teachers if it provides ready-made homework problems they can assign.) The point values here are an attempt to indicate the relative difficulty of the various challenges.