View source for Monte-Carlo methods for partial differential equations

Deterministic partial differential equations can be solved numerically by probabilistic algorithms such as Monte-Carlo methods, stochastic particle methods, ergodic algorithms, etc. (cf. also [[Differential equation, partial|Differential equation, partial]]).

The simplest example is the [[Heat equation|heat equation]] in <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101901.png" />:

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101902.png" /></td> </tr></table>

For a large class of functions <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101903.png" />, the unique solution is <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101904.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101905.png" /> is the heat kernel. Thus, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101906.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101907.png" /> is any [[Brownian motion|Brownian motion]] such that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101908.png" />, almost surely. The strong [[Law of large numbers|law of large numbers]] applies: if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m1101909.png" /> are independent copies of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019010.png" />, one has, almost surely,

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019011.png" /></td> </tr></table>

A Monte-Carlo procedure to compute <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019012.png" /> at a single point <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019013.png" /> consists in the simulation of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019014.png" /> independent copies of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019015.png" /> and the computation of

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019016.png" /></td> </tr></table>

This procedure requires the simulation of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019017.png" /> independent real-valued Gaussian random variables of mean <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019018.png" /> and variance <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019019.png" /> (cf. [[Random variable|Random variable]]; [[Gaussian process|Gaussian process]]). The convergence rate of the procedure is described by fluctuation results, large-deviation estimates, etc. (cf. also [[Probability of large deviations|Probability of large deviations]]). The most elementary estimate concerns the standard deviation of the random error:

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019020.png" /></td> </tr></table>

More generally, let <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019021.png" /> be an integro-differential operator (cf. also [[Integro-differential equation|Integro-differential equation]]), and consider the evolution problem

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019022.png" /></td> </tr></table>

From now on, it is assumed that there is a <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019023.png" />-valued [[Markov process|Markov process]] <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019024.png" /> such that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019025.png" /> (for a discussion on the assumptions that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019026.png" /> must satisfy, see [[#References|[a3]]], Chap. XIII, Sect. 3). Then, if it is possible to simulate independent copies of the process <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019027.png" />, one can approximate <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019028.png" /> by

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019029.png" /></td> </tr></table>

For example, in neutronics, the process <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019030.png" /> is the pair (position,velocity), and the velocity is a pure [[Jump process|jump process]] (in this case, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019031.png" /> is a transport operator) whose simulation is simple: one simulates the times at which the velocity changes and the new velocity after each change time independently (the change times have an exponential law).

Most often, the law of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019032.png" /> cannot be simulated exactly, but it may be possible to approximate it suitably. For example, when <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019033.png" /> is the solution of a [[Stochastic differential equation|stochastic differential equation]] driven by a Brownian motion (in this case, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019034.png" /> is a second-order differential operator), one can discretize the stochastic differential equation in time. The simplest discretization scheme (the Euler scheme) defines an approximate process <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019035.png" /> depending on the discretization step <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019036.png" />, whose simulation requires the simulation of independent Gaussian random variables only. The approximation formula is

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019037.png" /></td> </tr></table>

The global error of this procedure can be decomposed as the sum of

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019038.png" /></td> </tr></table>

(the discretization error) and of

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019039.png" /></td> </tr></table>

(the statistical error). Under appropriate assumptions, the discretization error can be expanded in terms of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019040.png" />.

The numerical procedure can be adapted to the solution of problems in bounded domains. Dirichlet boundary conditions lead to the simulation of processes stopped at the boundary. The simulation of reflected diffusion processes (cf. [[Diffusion process|Diffusion process]]) solves parabolic or elliptic partial differential equations with suitable [[Neumann boundary conditions|Neumann boundary conditions]].

Let <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019041.png" /> be the formal adjoint of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019042.png" />. Note that the probability distribution <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019043.png" /> of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019044.png" /> is a solution to

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019045.png" /></td> </tr></table>

in the sense of distributions (cf. [[Generalized function|Generalized function]]). The solution <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019046.png" /> can be approximated by the [[Empirical distribution|empirical distribution]]

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019047.png" /></td> </tr></table>

This procedure is used in neutronics.

To certain non-linear partial differential equations one can associate processes <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019048.png" /> which are not Markovian, but are such that the pairs <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019049.png" /> are Markovian, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019050.png" /> is the probability law of the random variable <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019051.png" />. Besides, in some cases, the law <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019052.png" /> of the process <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019053.png" /> can be described by means of the theory of propagation of chaos (see [[#References|[a4]]] and the references therein). That means that the law <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019054.png" /> is the limit in probability, when <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019055.png" /> goes to infinity, of the empirical distributions of the paths of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019056.png" /> particles in a mean field interaction. The interaction between the particles is related to the non-linear differential operator of the partial differential equation under consideration. If one can simulate the stochastic particle system, then an algorithm can be built from the formula

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019057.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019058.png" /> are the locations at time <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019059.png" /> of the interacting particles. For example, consider an interaction kernel <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019060.png" />, a family of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019061.png" /> independent Brownian motions, and the system

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019062.png" /></td> </tr></table>

The time discretization of such systems permits one to approximate equations such as Vlasov–Poisson–Fokker–Planck equations (cf. [[Vlasov–Poisson–Fokker–Planck system|Vlasov–Poisson–Fokker–Planck system]]), the viscous Burgers equation (cf. [[Diffusion equation|Diffusion equation]]; [[Turbulence, mathematical problems in|Turbulence, mathematical problems in]]), the incompressible Navier–Stokes equation for the vorticity in <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019063.png" /> (cf. [[Navier–Stokes equations|Navier–Stokes equations]]), etc. In the latter case, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019064.png" /> is the Biot–Savart kernel, and the algorithm is the Chorin random vortex method (see [[#References|[a2]]] for numerical considerations). Similarly, the simulation of interacting branching processes (cf. [[Branching process|Branching process]]) permits one to approximate the solutions of convection-reaction-diffusion equations. One difficulty of the mathematical study of these methods is to show that the statistical errors (on the distribution function of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019065.png" />, e.g.) have standard deviations of order <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m110/m110190/m11019066.png" /> although the particles are dependent. For theoretical aspects related to probabilistic algorithms for partial differential equations, including the [[Boltzmann equation|Boltzmann equation]], see, for instance, [[#References|[a1]]] and the references therein.

====References====
<table><TR><TD valign="top">[a1]</TD> <TD valign="top">  C. Graham,   T. Kurtz,   S. Méléard,   P. Protter,   M. Pulvirenti,   D. Talay,   "Probabilistic models for nonlinear partial differential equations" , ''Lecture Notes in Mathematics'' , '''1627''' , Springer  (1996)</TD></TR><TR><TD valign="top">[a2]</TD> <TD valign="top">  "Vortex methods and vortex motions"  K.E. Gustavson (ed.)  J.A. Sethian (ed.) , SIAM  (1991)</TD></TR><TR><TD valign="top">[a3]</TD> <TD valign="top">  J. Jacod,   "Calcul stochastique et problèmes de martingales" , ''Lecture Notes in Mathematics'' , '''714''' , Springer  (1979)</TD></TR><TR><TD valign="top">[a4]</TD> <TD valign="top">  A.S. Sznitman,   "Topics in propagation of chaos"  P.L. Hennequin (ed.) , ''Ecole d'Eté de Probabilités de Saint-Flour XI (1989)'' , ''Lecture Notes in Mathematics'' , '''1464''' , Springer  (1991)  pp. 165–251</TD></TR></table>