next up previous contents
Next: How to implement a Up: Example: Feynman Kac plugin Previous: Seed management   Contents

Feynman-Kac plugin description

This plugin solves Exercise 2.2 of the "Introduction to Parallel Computing" book described at page 82 [1].

We solve a partial differential equation inside an elliptical region by Monte Carlo simulations of the Feynman-Kac formula. The following partial differential equation is defined in a three-dimensional ellipsoid E:

$\displaystyle \frac{1}{2}\triangle u(x,y,z)-v(x,y,z)u(x,y,z)=0,
$

where the ellipsoidal domain is

$\displaystyle D= \Bigg\{(x,y,z) \in \mathbf{R}^3;
\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2} < 1 \Bigg\}.
$

The potential is

$\displaystyle v(x,y,z) =
2(\frac{x^2}{a^4}+\frac{y^2}{b^4}+\frac{z^2}{c^4}+\frac{1}{a^2}+\frac{1}{b^2}+\frac{1}{c^2}).
$

Our Dirichlet boundary condition is $ u=g(\mathbf{x})=1$ when $ \mathbf{x} \in \partial{D}.$

The goal is to solve this boundary value problem at $ \mathbf{x}=(x,y,z)$ by a Monte Carlo simulation. At the heart of the simulation lies the Feynman-Kac formula, which in our case ($ g=1$) is
$\displaystyle u(x,y,z)$ $\displaystyle =$ $\displaystyle \mathbf{E} \; g(\mathbf{X}(\tau)) \; exp\Bigg(-\int_{0}^{\tau}v(\mathbf{X}(s))ds \Bigg)$  
$\displaystyle u(x,y,z)$ $\displaystyle =$ $\displaystyle \mathbf{E} \; exp\Bigg(-\int_{0}^{\tau}v(\mathbf{X}(s))ds \Bigg)$  
$\displaystyle u(x,y,z)$ $\displaystyle =$ $\displaystyle \mathbf{E} \; Y(\tau).$  

which describes the solution $ u$ in terms of an expectation value of a stochastic process $ Y$ whose initial value $ Y(0)=1$. Here $ \mathbf{X}(s)$ is a Brownian motion starting from $ \mathbf{X}(0)=\mathbf{x}$ and $ \tau$ is its exit time from $ D$.

It is important to point out that the $ \mathbf{E}$ operator in simulation means $ \mathbf{E} f = \frac{1}{N} \sum_{i=1}^{N} f_{i}$ for samples of size $ N$.


The plugin pde3d.dll can be called in two different ways:

0,[a],[b],[c],feynmankac3d
1,[a],[b],[c],[initialx],[initialy],[initialz],feynmankac3d

The first call passes the ellipses axes as parameter. The function
feynmankac3d stored in pde3d.dll generates first randomly a point $ \mathbf{x}$ inside the ellipse. In the second call, we can specify the initial

$\displaystyle \mathbf{x}=([initialx],[initialy],[initialz]).
$

0 or 1 inform the plugin on how many parameters are loaded on the stack.

From this initial interior point $ \mathbf{X}(0)=\mathbf{x} \in D$ we integrate the following system of stochastic differential equations ( $ \mathbf{W}$ is a Brownian motion). We use $ N$ realizations of $ \mathbf{X}(t)$ to integrate.

$\displaystyle d\mathbf{X}$ $\displaystyle =$ $\displaystyle d\mathbf{W}$  
$\displaystyle dY$ $\displaystyle =$ $\displaystyle -v(\mathbf{X}(t))Ydt$  

For each realization, we integrate this set of equations until $ \mathbf{X}(t)$ exits at time $ t=\tau$ using the trapezoidal rule [1]. Finally, we compute $ u(\mathbf{x})=\mathbf{E}Y(\mathbf{X}(\tau))$ to get the solution.

To plot the error, we compare it to the exact analytical solution of the partial differential equation, computed by two differentiations:

$\displaystyle u(x,y,z)=exp(\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}-1).
$

The error in the plugin's graphical output is plotted as an oscillating green line around the y axis. After oscillating for a while, the error tends to converge to the y axis. Only one out of many random walks ( $ \mathbf{X}(t)$) is plotted with red or white colors.

Roughly speaking, we free a horde of random walkers from an initial point. These walkers diffuse while $ Y$ is integrated with rate $ -v(\mathbf{X})$ on the potential $ v$ until they reach the boundary of the ellipse. All walks are averaged in a similar way as it is done with the pi plugin to get the function value at the chosen initial point. The diffusion process of the walkers, distributed as a trivariate normal density, weights points near the initial point more than points far away.

Figure 13: GPU computing the Feynman-Kac problem and corresponding frontend.


next up previous contents
Next: How to implement a Up: Example: Feynman Kac plugin Previous: Seed management   Contents
Tiziano Mengotti 2004-03-27