#!/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, " nX=[number]" sys.exit(1) # ---------------------------------------------------------------- # SDE: # dX_t = -1/(1+t) X_t dt + 1/(1+t) db_t. # Solution: # X_t = b_t / (1+t). T = 10 nt = 2000 nX = 5 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'^nX=', sys.argv[argi]): nX = int (sys.argv[argi][3:]) else: usage() dt = T/nt sqrtdt = sqrt(dt) t = 0.0 Bt = [0.0] * nX Xexact = [0.0] * nX Xapprox = [0.0] * nX for k in range(0, nX): Xapprox[k] = Bt[k] / (1+t) while (t <= T): for k in range(0, nX): Xexact[k] = Bt[k] / (1+t) print "%11.7f" % (t), for k in range(0, nX): print "%11.7f %11.7f %9.3e" % \ (Xexact[k], Xapprox[k], Xexact[k]-Xapprox[k]), print for k in range(0, nX): dB = random.gauss(0, sqrtdt) Xapprox[k] += -1/(1+t) * Xapprox[k] * dt + 1/(1+t) * dB Bt[k] += dB t += dt