Showing changes from revision #6 to #7: Added | Removed | Changed

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, (also using see thePython) 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?