Numerically solving PDEs in Mathematica using finite difference methods

Mathematica’s NDSolve command is great for numerically solving ordinary differential equations, differential algebraic equations, and many partial differential equations. Most of the integration details are handled automatically, out of the user’s sight. NDSolve switches between integration schemes based on the problem at hand, adapting step sizes and monitoring stiffness as it goes. Advanced users can override these options, customizing NDSolve to their needs.

Sadly, some types of PDEs are beyond NDSolve’s capabilities. Confronted with one of these PDEs, a user must resort to a more “manual” procedure to find a numerical solution. In this post, we’ll examine a few tricks that can make this process easier.

Consider the following PDE:

We seek a solution, f(x,y) on the domain [0,10]x[0,10]. If you try to enter this elliptic PDE into NDSolve, Mathematica will vigorously protest. Instead, you can try to implement a finite difference method.

First, we will divide the domain into a grid. To keep things simple, we will use 1×1 squares.

We next make an array of the (as of yet unknown) values that f takes at each point of the grid.

At each grid point, we will replace the derivatives in our PDE with linear combinations of values at neighboring grid points. If you have ever taken Calculus, this idea should sound very familiar– it is essentially like approximating the slope of a line tangent to a graph with the slope of a secant line.

To determine which linear combination of values to use for the different derivatives at the different grid points, we make use of an obscure Mathematica command, NDSolve`FiniteDifferenceDerivative.

The first argument is the type of derivative (for example {2,0} means the second-order derivative with respect to x). The second argument is the grid. The third argument is the array of (unknown) values.

This returns an array. Each entry in the array is the approximation to desired derivative at a particular grid point. For example, one such entry is:

We can see that the second order derivative with respect to x at the point (3,4) is expressed as a linear combination of values of f. Note that it is at the point (3,4) not (4,5) because the part specifications start at index 1, while the grid coordinates start at zero.

With this tool, we can replace the partial differential equation for the unknown function f with a large number of linear algebraic equations for the unknown values f(i,j). It is necessary to replace some of these equations with appropriate boundary conditions (in our example, two of the boundaries are absorbing, and two are reflecting). Once this is done, we can solve for all of the f(i,j) values.

I’ll spare you the code, and just show the resulting picture of the approximation to the solution:

In conclusion, we can see that, even when NDSolve can’t handle a PDE, Mathematica has some hidden gems that can make life easier!

The following two tabs change content below.

Ben Nolting

Ph.D. Candidate at University of Nebraska

37 thoughts on “Numerically solving PDEs in Mathematica using finite difference methods

      1. shaahin

        hi ben
        i have been trying to solve an extremely nonlinear 2nd order hyperbolic pde (2 vars) with mathematica with no success…can i ask for some pointer or possible an elementary sample code to get started?! help is highly appreciated!

  1. Nick

    Nice post! Makes pdes not so scary.

    I am encountering problems with the boundary conditions as well. Could you please send me the notebook too? Thanks.

  2. Pingback: Numerically solving PDEs in Mathematica using finite difference … | Solve Math & Science Problems -

  3. Tania

    Thanks for your post, it looks very clear, but as the others I’d like to check the code to see how you define the boundary conditions. Could you please send it to me too?

  4. Joseph

    Hi, Ben, thanks a lot. I have solved PDE using C++ but don’t know that mathematica can solve it so neatly. Could I have a copy of your code? Many thanks

  5. Ben Nolting Post author

    Hi Joseph! I’ll send you the code now. Mathematica 9 was released this week and it his many new features for solving PDE’s. Sadly, a fully automated implementation of the finite difference method is not among them. So, for now at least, the type of work arounds you will see in the code I’ll send at necessary.

  6. Pingback: Finite Difference Method (now with free code!) | Datavore Consulting | Blog

  7. Charles Hagwood

    Ben does this code use adaptive grid refinement? I need this because my coefficient c(x,y)=c(x) in front of \partial u/\partial x is a step function in x.

  8. Adokiye

    Hi Ben I just hope that you still remember my request concerning the 2D Diffusion with non rectangular boundary condition, which you suggest it will be fine in polar coordinate, and with NDSolve function in mathematica.

  9. bmarchi

    There doesn’t appear to be any problems directly solving the PDE in Mathematica. I used the following implementation:


    yielding an interpolating function that can be evaluated numerically and visualized by:

    Plot3D[Evaluate[u[x, y] /. sol], {x, 0, 10}, {y, 0, 10}, PlotRange -> All]

    Hope this helps

    1. Ben Nolting Post author

      Yes, Mathematica 10 has new capabilities for PDE’s. In addition to the built-in stuff with NDSolve, you can also check out the package “NDSolveFEM” for finite element methods.

    1. Majid

      Hi Ben,
      I want to solve the Laplace equation with Neumann boundary conditions. Can this code help me? would you care to send it to me too?

    1. Ben Nolting Post author

      Hi Sajid! I’d be happy to try to help. Feel free to e-mail me with the details about your problem. It will probably be a week or two before I have time to help, but if you don’t mind the delay, I’ll take a look at it.

  10. Yousef Alqawasmeh

    Hi Ben,

    I am working on a system of reaction diffusion equations with interface conditions (conditions inside the domain). I need your help to solve the system numerically. Can you please send me your email in order to share the details with you.

    Thank you in advance.

  11. Lu

    Can I have the code please? I would like to solve a higher order nonlinear pde that your code will be very helpful. Thank you!

  12. Samir A. El-Tantawy

    Could you help me (OR i want a code) to find the numerical solution of the following nonlinear Schrödinger equation using finite difference method
    i(??/?t)+(1/2)P(?²?/?x²)+Q|?|²?=0, x?[-L,L]
    we can put P=0.137 and Q=4.54.
    The initial solution is
    The following Dirichlet boundary conditions are used


Leave a Reply