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.

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.

1 2 3 |
xgrid = Range[0, 10]; ygrid = Range[0, 10]; grid = Outer[{#1, #2} &;, xgrid, ygrid]; |

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

1 |
values = Outer[f[#1, #2] &, xgrid, ygrid]; |

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.

1 |
NDSolve`FiniteDifferenceDerivative[{2, 0}, {xgrid, ygrid}, 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:

1 |
NDSolve`FiniteDifferenceDerivative[{2, 0}, {xgrid, ygrid}, values][[4, 5]] |

1 |
-(1/12) f[1, 4] + 4/3 f[2, 4] - 5/2 f[3, 4] + 4/3 f[4, 4] - 1/12 f[5, 4] |

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!

#### 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

pratipCan I get the code please? I am not able to incorporate the IC and BC.

BenPost authorI’ll send it to you in a Mathematica notebook 🙂

DanThank you for the post. Could I please have the code as well?

Ben NoltingPost authorSure! I’ll send it now.

shaahinhi 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!

Ben NoltingPost authorSure! Feel free to send the problem to me, and I will see what I can do.

Ben NoltingPost authorHere is the latest, fastest version (and includes mixed partials):

Newest Finite Difference Code

Mateus BezerraCan you send me the code, please? =) shewlong@gmail.com

Ben NoltingPost authorDone 🙂

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

Ben NoltingPost authorI’ll send it to you, Nick! Thanks for reading 🙂

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

TaniaThanks 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?

Ben NoltingPost authorOk! I’ll send you a notebook.

JevExcellent and very very helpful post – but I am still learning Mathematica! May I have the notebook as well?

Ben NoltingPost authorI will send it! Let me know if you have any questions.

JosephHi, 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

Ben NoltingPost authorHi 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.

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

Ben NoltingPost authorThere is now a much nicer notebook, with comments and explanations, available for download at: http://datavoreconsulting.com/blog/uncategorized/finite-difference-method-now-free-code/

Charles HagwoodBen 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.

Ben NoltingPost authorThe short code (available here, does not have an adaptive grid. You might be able to modify it for an adaptive grid, though.

MikeHi – please could you send me the notebook? I’d really appreciate it.

Ben NoltingPost authorHi Mike,

I will send it to you. You can also download it from our more recent post, http://datavoreconsulting.com/programming-tips/finite-difference-method-now-free-code/

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

cheers.

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

sol=NDSolve[{x*(5+y)*D[u[x,y],{x,2}]+y/2*(x+5)*D[u[x,y],{y,2}]-x*(5-y)*D[u[x,y],x]-y/2*(x-5)*D[u[x,y],y]==-1,u[0,y]==0,u[x,0]==0,(u^(1,0))[10,y]==0,(u^(0,1))[x,10]==0},u,{x,0,10},{y,0,10}]

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

Ben NoltingPost authorYes, Mathematica 10 has new capabilities for PDE’s. In addition to the built-in stuff with NDSolve, you can also check out the package “NDSolve

`FEM`

” for finite element methods.Fakhar AbbasCan I get the code of this particular equation please, so that I can get help from this to my problem.

MajidHi 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?

Ben NoltingPost authorYes, I think this should work for a Laplace equation with Neumann BC. Here is a link to a more recent blog post, where you can download the code.

http://datavoreconsulting.com/programming-tips/finite-difference-method-now-free-code/

Additionally, Mathematica has added some features in version 10 that could be useful to you. The pakcage “NDSolve

`FEM`

“contains some tools for finite element methods.MajidThank you Ben. You’re a life saver.

sajid aliHi Ben..

I need to solve hyperbolic equation for axially moving string can you help me in that???

Ben NoltingPost authorHi 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.

Yousef AlqawasmehHi 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.

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

Ben NoltingPost authorSure! I will e-mail it to you. Thanks for your interest.

Samir A. El-TantawyCould 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

?(x,t)=?(P/Q)*[-1+4(1+2iPt)/(1+4x²+4P²t²))]exp(iPt).

The following Dirichlet boundary conditions are used

?(-L,?)=?(L,?).

Samir A. El-TantawyCould 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(du/dt)+(1/2)P(d²u/dx²)+Q|u|²u=0, x:[-L,L]

we can put P=0.137 and Q=4.54.

The initial solution is

u(x,t)=Sqrt(P/Q)*[-1+4(1+2iPt)/(1+4x²+4P²t²))]exp(iPt).

The following Dirichlet boundary conditions are used

u(-L,t)=u(L,t).

Detoshcould you send me the code as well. I have an ODE to solve. and i have been trying to write a code and i am running into issues

DetoshI have a PDE rather