A couple of months ago, we wrote a **post** on how to use finite difference methods to numerically solve partial differential equations in Mathematica. Several readers have asked for more details about the method. To help everyone out, we are posting a Mathematica notebook that contains explanations and code. Download the notebook by clickingÂ **here**.

The notebook will implement a finite difference method on elliptic boundary value problems of the form:

The comments in the notebook will walk you through how to get a numerical solution. An example boundary value problem is solved, yielding a solution that looks like this:

Last week, Mathematica 9 was released. It has some awesome new features, including enhanced NDSolve capabilities. Sadly, there is no built-in finite difference method solver yet. We will have to continue to use workarounds like the notebook posted here for a while longer.

#### Ben Nolting

#### Latest posts by Ben Nolting (see all)

- Limitations of the negative binomial distribution in spatial models - August 2, 2013
- Volume Rendering and Large Data Sets - March 29, 2013
- Custom interfaces in Mathematica using Dynamic Module - February 4, 2013

Charles HagwoodSmall error in code:

ygrid = Range[ymin, ymax, dx];

Change dx to dy

Ben NoltingPost authorGreat catch! Thanks.

Charles HagwoodAlso,

change dx to dy in

dataPoints =

Table[{xmin + i*dx, ymin + j*dx, solutionArray[[i + 1, j + 1]]}, {i,

0, Length[xgrid] – 1}, {j, 0, Length[ygrid] – 1}];

Did you check your code?

Ben NoltingPost authorOk, I think the version posted now has all of the dx and dy errors corrected. The code was originally used on problems with square grid spacing, so the issue didn’t come up.

Thanks for finding the mistakes in here… this was just a quick afternoon blogging exercise, so it doesn’t surprise me that it isn’t perfect. I appreciate your help.

Dimitris LeivaditisExcellent code. I am trying to imply also some condition around the origin, e.g. that u=5*(x+y) {x,-0.9,-0.7}, {y,-0.9,-0.7} but I cannot solve it. Do you have any idea how can this be done with your code? Thanks in advance.

Ben NoltingPost authorThank you!

For Poisson’s equation to be a well-posed problem, you can only specify boundary conditions. In your example, you could specify that u=5*(-0.9+y) as the boundary condition along the line x= -0.9, and specify similar boundary conditions along x= -0.7, y= -0.9, and y= -0.7. I hope that helps a little.

Dimitris LeivaditisI think I found it. I create another boundary condition

origin = Table[

U[[k + 1, k + 1]] – qorigin[xmin + k*dx, ymin + k*dy], {k, 1,

Length[ygrid] – 18}];

where qorigin=5*(x+y)

and proceed accordingly.

Thanks.

AdokiyeHow do I write the code to solve 2D Diffusion, with t component? I would grateful if this can be done

AdokiyeI will be gratelful.

benHello, really nice program!! I have a small problem, when you make the array of values for the solution, you use

U = Array[u, {xdivisions + 1, ydivisions + 1}, {{xmin, xmax}, {ymin, ymax}}];

This seems a quite complicated form to me. Is there any special reason you make it like this??

I change it to

U = Array[u, {xdivisions + 1, ydivisions + 1}, {xmin, ymin}];

and it works.

Ben NoltingPost authorThanks for your comment, ben! The purpose of specifying the min, max, and number of divisions separately is so that you can control both the range of the mesh and the size of the grid cells. But you are right, if you have one of those parameters specified beforehand, than it is better to streamline things like you suggest.

FranzDoes a nonlinear system of PDEs work the same way?

What has to be changed?

Thanks.