Stochastic Hopf bifurcation in Python (Rev #7)

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[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()

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.

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?

category: computational methods, software