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.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VI of the course APMA0340
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

 

Return to Mathematica page

Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations

Return to the Part 7 Special Funtions