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 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 here. (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.)

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.

category: methodology, software