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.
#!/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.lower() param_value = float(split) 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 = x_init aY = 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()
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 plotted against , with and , would be the following:
And this is Figure 2, showing as a function of time , 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.
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.
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.)