#!/usr/bin/python -Wall # ================================================================ # John Kerl # kerl at math dot arizona dot edu # 2008-05-12 # ================================================================ from __future__ import division # 1/2 = 0.5, not 0. import random from math import * import copy import sys import re # ---------------------------------------------------------------- def usage(): print >> sys.stderr, "Usage: %s [options]" % (sys.argv[0]) print >> sys.stderr, "Options:" print >> sys.stderr, " T=[number]" print >> sys.stderr, " nt=[number]" print >> sys.stderr, " A=[number]" print >> sys.stderr, " B=[number]" sys.exit(1) # ---------------------------------------------------------------- # Brownian motion on an ellipse. # SDE: # dX_t = -1/2 X dt - A/B Y_t db_t # dY_t = -1/2 Y dt + B/A X_t db_t # Solution: # X_t = A cos(b_t) # Y_t = B sin(b_t). T = 2.0 nt = 20000 A = 2.0 B = 2.0 argc = len(sys.argv) for argi in range(1, argc): if re.match(r'^T=', sys.argv[argi]): T = float(sys.argv[argi][2:]) elif re.match(r'^nt=', sys.argv[argi]): nt = int (sys.argv[argi][3:]) elif re.match(r'^A=', sys.argv[argi]): A = float(sys.argv[argi][2:]) elif re.match(r'^B=', sys.argv[argi]): B = float(sys.argv[argi][2:]) else: usage() dt = T/nt t = 0.0 Bt = 0.0 sqrtdt = sqrt(dt) Xapprox = A * cos(Bt) Yapprox = B * sin(Bt) while (t <= T): Xexact = A * cos(Bt) Yexact = B * sin(Bt) print "%11.7f %11.7f %11.7f %11.7f %11.7f %9.3e %9.3e" % \ (t, Xexact, Xapprox, Xexact-Xapprox, Yexact, Yapprox, Yexact-Yapprox) dB = random.gauss(0, sqrtdt) Xapprox += -0.5 * Xapprox * dt - A/B * Yapprox * dB Yapprox += -0.5 * Yapprox * dt + B/A * Xapprox * dB Bt += dB t += dt