Finite Difference Method (now with free code!)

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.

The following two tabs change content below.

Ben Nolting

Ph.D. Candidate at University of Nebraska

14 thoughts on “Finite Difference Method (now with free code!)

  1. Charles Hagwood

    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?

    1. Ben Nolting Post author

      Ok, 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.

  2. Dimitris Leivaditis

    Excellent 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.

  3. Ben Nolting Post author

    Thank 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.

  4. Dimitris Leivaditis

    I 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.

  5. Adokiye

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

  6. ben

    Hello, 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.

    1. Ben Nolting Post author

      Thanks 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.

    1. Ben Nolting Post author

      It will work, but you will need to modify the code in a few places. In particular, the second order terms (dudx2 and dudy2) will have to be changed to fourth-order FD approximations. This is easily done by changing 2 to 4’s in the code, thanks to Mathematica’s built-in FiniteDifferenceDerivative function. I hope that helps!


Leave a Reply