#!/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, " X0=[number]" print >> sys.stderr, " Y0=[number]" print >> sys.stderr, " alpha=[number]" print >> sys.stderr, " beta=[number]" sys.exit(1) # ---------------------------------------------------------------- # Stochastically forced vibrating string equation: # dX = Y dt + alpha da # dY = -X dt + beta db # where a and b are independent Brownian motions. T = 100 nt = 10000 X0 = 1 Y0 = 0 alpha = 0.200 beta = 0.200 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'^X0=', sys.argv[argi]): X0 = float(sys.argv[argi][3:]) elif re.match(r'^Y0=', sys.argv[argi]): Y0 = float(sys.argv[argi][3:]) elif re.match(r'^alpha=', sys.argv[argi]): alpha = float(sys.argv[argi][6:]) elif re.match(r'^beta=', sys.argv[argi]): beta = float(sys.argv[argi][5:]) else: usage() t = 0 dt = T/nt At = 0.0 Bt = 0.0 X = X0 Y = Y0 sqrtdt = sqrt(dt) while (t <= T): print "%11.7f %11.7f %11.7f" % (t, X, Y) dA = random.gauss(0, sqrtdt) dB = random.gauss(0, sqrtdt) X += Y * dt + alpha * dA Y -= X * dt + beta * dB At += dA Bt += dB t += dt