# 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, NDSolveFiniteDifferenceDerivative.

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 ## 61 thoughts on “Numerically solving PDEs in Mathematica using finite difference methods”

1. pratip

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

1. Ben Post author

I’ll send it to you in a Mathematica notebook 🙂

1. Dan

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

1. Ben Nolting Post author

Sure! I’ll send it now.

2. 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. Ben Nolting Post author

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

1. Ben Nolting Post author

Done 🙂

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

3. Ben Nolting Post author

I’ll send it to you, Nick! Thanks for reading 🙂

4. 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?

1. Ben Nolting Post author

Ok! I’ll send you a notebook.

5. Jev

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

1. Ben Nolting Post author

I will send it! Let me know if you have any questions.

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

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

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

1. Ben Nolting Post author

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

9. Mike

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

10. 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.
cheers.

11. bmarchi

There 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

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.

12. Fakhar Abbas

Can I get the code of this particular equation please, so that I can get help from this to my problem.

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

Yes, 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 “NDSolveFEM`“contains some tools for finite element methods.

1. Majid

Thank you Ben. You’re a life saver.

13. sajid ali

Hi Ben..

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

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.

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

15. 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!

1. Ben Nolting Post author

Sure! I will e-mail it to you. Thanks for your interest.

16. 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
?(x,t)=?(P/Q)*[-1+4(1+2iPt)/(1+4x²+4P²t²))]exp(iPt).
The following Dirichlet boundary conditions are used
?(-L,?)=?(L,?).

1. 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(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).

17. Detosh

could 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

18. Detosh

I have a PDE rather

1. Ben Nolting Post author

I apologize for the delayed reply. I’ve sent the code to you. Thanks!

19. MuhammadIrfan

Hi I need a finite difference scheme code for the Korteweg-de Veries (KdV) equation. Can anyone help me?

1. Ben Nolting Post author

I apologize for the delayed reply. I’ve sent the code to you for an example system, which you can adapt to the KdV equation.

20. Fahim

Can I get the code please

1. Ben Nolting Post author

Sure! I apologize for the delayed response.

1. Ben Nolting Post author

My e-mail got bounced. Feel free to e-mail me if you would like the code.

21. Maha

Is it valid to solve coupled system of PDEs?
Would you please provide me with the code?

1. Ben Nolting Post author

I’ve sent the code for the single equation example. It will require some modification for a system of PDEs, but the same general framework will work.

22. climate503

Would you send me a PDE code with Boundary Conditions Implied?

1. Ben Nolting Post author

Sure!

1. GGG

Could I please have the Mathematica code?

1. Ben Nolting Post author

Sure! I’ll send it now.

1. brian luna

if possible, could i receive the code too ?

23. Ale

Please, can you send me the code? Thanks

1. Ben Nolting Post author

Sure! I’ll send it.

24. tommy

Could you please send me the code ? Thanks.

1. Ben Nolting Post author

I will! 🙂

25. sunitha G K

Please, can you send me the mathematic a code for coupled pde.

1. Ben Nolting Post author

Sure!

26. Zubair

Hi Ben , i appreciate your work .. i am working on nonlinear coupled ODE and PDE.
Kindly send me the finite difference code on mzaqhashmi@gmail.com
Thanks

1. Ben Nolting Post author

Thanks for reading! I’ll send the code to you.