# Preface

A numerical approximation to the solution of a linear elliptic partial differential equation at a single point. Simulation of the motion of a random particle may be used to approximate the solution to linear elliptic equations.

Introduction to Linear Algebra with Mathematica

# Monte Carlo for Elliptic Equations

Simulation of the motion of a random particle may be used to approximate the solution to linear elliptic equation. In particular, consider the two-dimensional Dirichlet boundary value problem:

$\begin{split} L[u] &= f(x,y), \qquad x,y \in R, \\ u &= \phi (x,y,t) , \qquad x,y \in \partial R, \end{split}$
with the differential operator L[ ] defined by
$L[u] = A\,u_{xx} + 2B\,u_{xy} + C\,u_{yy} + D\,u_x + E\,u_y ,$
where A, B, C, D, E, F are all functions of {x,y,t}. Here R is two dimensional domain on ℝ² with smooth boundary ∂R.

The steps for this method are straightforward. First, approximate the given elliptic partial differential equation by a finite difference method. Rewrite the finite difference formula as a recursive function for the value of the unknown at any given point. Then interpret this recursive formula as a set of transition probabilities that determine the motion of a random particle.

Now, write a computer program that will allow many (say K) particles to wander randomly around the domain of interest, based on the transition probabilities found from the difference formula. Simulate particles one at a time, with every particle starting off at the same point (say the point $\bfz$).

If the boundary data are not necessarily of the Dirichlet type (i.e., Neumann or mixed boundary conditions) then, when the particles reach the boundary, they will be given a finite probability to leave the boundary, and re-enter the domain of the problem. If the particle leaves the boundary, then continue the iteration process. If it does not leave the boundary, then the value at the boundary is stored away, and a new particle is started off at the point~$\bfz$.

The simulation is finished after all $K$ particles have been absorbed into the boundary. If the original elliptic equation was homogeneous, then an approximation to the solution, at the point $\bfz$, will be given by the average of all the values obtained (recall that when the particles stop at the boundary they obtain a value).

step = 0.10;
SUM = 0.0;
disp("Number of particles    Approximation")
for i = 1:10000
X = [2;0];
r = norm(X);
while r <= 3 && r > 1
X = X + step*sign(rand(2,1) - .5);
r = norm(X);
end
if r <= 1 % When a particle hits the boundary, sum the value
SUM = SUM + 4;
end
if  r >= 3
SUM = SUM + 6;
end
if mod(i,1000) == 0
fprintf("%10d %21.3f \n", i, SUM/i)
end
end