#!/usr/bin/python -Wall # ================================================================ # John Kerl # kerl.john.r@gmail.com # 2008-01-30 # ================================================================ from __future__ import division # 1/2 = 0.5, not 0. import random import kerlutil # ---------------------------------------------------------------- # 1D Brownian motion # As is standard in Python, use an optional argument to implement what would be # called a "static auto" in C or a "save parameter" in Fortran -- a local # variable which is saved between calls. def B(t, dt, uprev=[0.0]): du = random.gauss(0, dt) uprev[0] += du return uprev[0] # ---------------------------------------------------------------- # 2D Brownian motion def B2(t, dt, uprev=[0.0, 0.0]): du = random.gauss(0, dt) dv = random.gauss(0, dt) uprev[0] += du uprev[1] += dv return [uprev[0], uprev[1]] # ---------------------------------------------------------------- def do_1D(): T = 10.0 dt = 0.01 ts = kerlutil.mfrange(0, dt, T) nt = len(ts) # Brownian motion. # Run this first to find out what B(T) turns out to be. Bts = [] for t in ts: Bt = B(t, dt) Bts.append(Bt) BT = Bts[-1] # Brownian bridge: R(t) = B(t) - (t/T) B(T). Rts = [] for i in range(0, nt): t = ts[i] Bt = Bts[i] Rt = Bt - (t/T) * BT Rts.append(Rt) for i in range(0, nt): t = ts[i] Bt = Bts[i] Rt = Rts[i] print "%11.6f %11.6f %11.6f" % (t, Bt, Rt) # ---------------------------------------------------------------- def do_2D(): T = 10.0 dt = 0.01 ts = kerlutil.mfrange(0, dt, T) nt = len(ts) # Brownian motion. # Run this first to find out what B(T) turns out to be. Bts = [] for t in ts: Bt = B2(t, dt) Bts.append(Bt) BT = Bts[-1] # Brownian bridge: R(t) = B(t) - (t/T) B(T). Rts = [] for i in range(0, nt): t = ts[i] Bt = Bts[i] Rt = [ Bt[0] - (t/T) * BT[0], Bt[1] - (t/T) * BT[1] ] Rts.append(Rt) for i in range(0, nt): t = ts[i] Bt = Bts[i] Rt = Rts[i] #print "%11.6f %11.6f %11.6f" % (t, Bt[0], Bt[1]) print "%11.6f %11.6f %11.6f" % (t, Rt[0], Rt[1]) # ---------------------------------------------------------------- #do_1D() do_2D()