#!/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, " c=[number]" print >> sys.stderr, " cos" print >> sys.stderr, " sin" print >> sys.stderr, " nX=[number]" sys.exit(1) # ---------------------------------------------------------------- def get_mean(X): n = len(X) sum = 0.0 for i in range(0, n): sum += X[i] return sum / n # ---------------------------------------------------------------- # X_t = exp(t/2) cos(Bt) T = 2 nt = 1000 c = 0.5 do_cos = 1 nX = 40 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'^c=', sys.argv[argi]): c = float(sys.argv[argi][2:]) elif sys.argv[argi] == 'cos': do_cos = 1 elif sys.argv[argi] == 'sin': do_cos = 0 elif re.match(r'^nX=', sys.argv[argi]): nX = int (sys.argv[argi][3:]) else: usage() dt = T/nt t = 0 sqrtdt = sqrt(dt) Bt = [0.0] * nX X = [0.0] * nX while (t <= T): for k in range(0, nX): if (do_cos): X[k] = exp(c*t) * cos(Bt[k]) else: X[k] = exp(c*t) * sin(Bt[k]) print "%11.7f" % (t), for k in range(0, nX): print "%11.7f" % (X[k]), print "%11.7f %11.7f %11.7f" % (get_mean(X), exp(c*t), -exp(c*t)) for k in range(0, nX): dB = random.gauss(0, sqrtdt) Bt[k] += dB t += dt