#!/usr/bin/python -Wall # ================================================================ # John Kerl # kerl at math dot arizona dot edu # 2008-06-24 # ================================================================ from __future__ import division import sys from math import * import random import kerlutil # ================================================================ # Feynman-Kac formula: # # Function: phi(t,x). # # dphi f(x)^2 d^2 phi dphi # PDE: ---- = ------ ------- + e(x) ---- + g(x) phi # dt 2 dx^2 dx # # with initial conditions # phi(0,x) = h(x). # # Solution # # phi(t,x) = E^x[ h(X_t) exp{int_0^t g(X_u) du} ] # # with # # dX_t = e(X_t) dt + f(X_t) db_t # # where E^x[...] means E[... | X_0 = x]. # ================================================================ def h(x): return 1 / (1 + x**2) # ---------------------------------------------------------------- def f(x): return 1.414214 # ---------------------------------------------------------------- def e(x): return 1.0 # ---------------------------------------------------------------- def g(x): #return -0.1 * x**2 return 0.0 # ---------------------------------------------------------------- Exreps = 10000 xlo = -2 xhi = 2 dx = 0.1 xmesh = kerlutil.mfrange(xlo, dx, xhi) thi = 2.0 dt = 0.01 tmesh = kerlutil.mfrange(0.0, dt, thi) sqrtdt = sqrt(dt) for x in xmesh: bts = [0.0] * Exreps Xts = [x] * Exreps Is = [0.0] * Exreps for t in tmesh: Ex = 0.0 for i in range(0, Exreps): bt = bts[i] Xt = Xts[i] Is[i] += g(Xt) * dt Ex += h(Xt) * exp(Is[i]) db = random.gauss(0, sqrtdt) bts[i] += db Xts[i] += e(Xt) * dt + f(Xt) * db Ex /= Exreps print "%11.7f %11.7f %11.7f" % (x, h(x), Ex) sys.stdout.flush()