previous up next ','..','$myPermit') ?> SYLLABUS  Previous: 5.3 Particle orbits  Up: 5 MONTE-CARLO METHOD  Next: 5.5 When should you


5.4 A scheme for the advection diffusion equation

Slide : [ scheme - Box-Mueller || VIDEO 99) echo " modem - ISDN - LAN "; else echo "login";?> ]

Let us use now the Feynman-Kac theorem and evolve a large number of stochastic particle orbits to approximate the advection diffusion equation

$\displaystyle \frac{\partial f}{\partial t}+\frac{\partial }{\partial x}\left( {v}f\right) - \frac{\partial^2}{\partial x^2}\left( \frac{b^2}{2}f \right) = 0$ (1)

The rules of the Itô calculus require the time discretization to be explicit

$\displaystyle {X}_{t+\Delta t}={X}_t + {v}({X}_t,t)\Delta t + b({X}_t,t)({W}_{t+\Delta t}-{W}_t)$ (2)

Note that it is in fact possible to construct implicit discretizations of a stochastic differential equation [22], but the construction it not as simple as it is for ordinary equations (5.3#eq.3). Since $ (W_{t+\Delta t}-W_t)\in\mathcal{N}\left(0,\sqrt{\Delta t}\right)$ , the Wiener process can be written explicitly

$\displaystyle {W}_{t+\Delta t}-{W}_t={\zeta} \sqrt{\Delta t}$ (3)

where $ {\zeta}\in\mathcal{N}(0,1)$ are normally distributed random numbers. The central limit theorem states that any large sum involving $ n$ random numbers equally distributed with zero mean and unit variance will eventually converge to $ \mathcal{N}(0,\sqrt{n})$ . Any such random number could therefore be used for $ {\zeta}$ provided that the number of time steps used to approximate the differential is large enough; choosing $ {\zeta}\in\mathcal{N}(0,1)$ however leads to the fastest convergence.

A numerical Monte Carlo scheme for the advection-diffusion equation is now constructed, using $ N$ particles to sample the possible outcomes of the entire ensemble of orbits

$\displaystyle {X}_{t+\Delta t}={X}_t + {v}({X}_t,t)\Delta t + {\zeta} b({X}_t,t) \sqrt{\Delta t}$ (4)

This is simply be implemented as 
      for(int j = 0; j < numberOfParticles; j++){
        particlePosition[j] += velocity * timeStep +
          random.nextGaussian() *
          Math.sqrt(2 * diffusCo * timeStep);
        // Insert here your periodic boundary conditions (exercise 5.03)
      }
and can be tested in the JBONE applet below
JBONE applet:  press Start/Stop to solve the advection-diffusion equation with 1000 random walkers initialized so as to approximate a box function. Switch to Particle methods* to display the position of each particle (red dots) in addition to the density distribution (black line) obtained from a linear assignment of the position on a grid.
"; if ($user_nbr<100) echo " "; echo "

"; ?>

Numerical experiments: Monte Carlo scheme for advection-diffusion
  1. Compare the solution obtained with 128 steps of duration TimeStep=0.5 with the one achieved with a single step of duration TimeStep=64; which is faster and which is more precise?
  2. Observe how the periodic boundary conditions have been implemented so that particles can take steps that are larger than the domain size (exercise 5.03).

JAVA is one of the few programming languages that has a normal pseudo random numbers generator  $ \mathcal{N}(0,1)$ . Other programming languages generally provide only uniform pseudo random numbers  $ \mathcal{U}(0,1)$ and the Box Müller method can then be used to transfrom them into  $ \mathcal{N}(0,1)$ with the algorithm:

The rest of the implementation is equivalent if one remembers the definition of the command x+=dx as being equivalent to x=x+dx.

SYLLABUS  Previous: 5.3 Particle orbits  Up: 5 MONTE-CARLO METHOD  Next: 5.5 When should you