# 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 $Y$ plotted against $X$, with $\beta = -0.25$ and $\lambda = 0.1$, would be the following:

And this is Figure 2, showing $X$ as a function of time $t$, 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.

• (5 points) Modify the program so that the appropriate bits of text are nicely formatted; e.g., have it say $\beta$ instead of beta. (Hint: dollar signs make the world go round. But, be forewarned: both TeX and Python use the backslash as the escape character, so you may end up using more of them than you anticipated.)
• (5 points) Use the program to explore parameter space: try it out for the values of $\beta$, $\lambda$, $x_0$ and $y_0$ used to make the figures on the Hopf bifurcation page.
• (10 points) Modify the program so that it reads parameter values from a text or CSV file instead of the command line. You’ll have to decide how that file is formatted! Try to devise a scheme which you can reuse for another simulation.
• (10 points) Add an option to the program so that it can write its results to a file, and write a script which can read in that file and display it using scipy and pylab.
• (10 points) Write a program which runs the simulation twice and displays the results of both runs on the same plot. Try using colours and/or line widths to distinguish them.
• (15 points) Have the program output a third figure: the histogram of $r = \sqrt{x^2 + y^2}$, as shown on the Fokker-Planck equation example page. Experiment: do transient effects show up in your histogram? What effect does sampling different portions of the simulation results have on the histogram?