# Numerical technique to solve 2nd order PDE

Currently I am trying to numerically solve a linear partial differential equation of the form $$\sum_{l=0}^2 \sum_{m=0}^2 \sum_{n=0}^2 C_{lmn} (x,y,\phi) \frac{\partial^l \partial^m \partial^n}{\partial x^l \partial y^m \partial \phi^n } \Psi(x,y,\phi) = 0$$ with the following domain $x \geq 0$, $y \geq 0$, and $0 \leq \phi \leq \pi$. $\Psi$ and $C_{lmn}$ are real functions. The boundary conditions are $$\Psi(0,y,\phi) = \Psi(x,0,\phi) = \Psi_x(0,y,\phi) = \Psi(x,y,0) = \Psi_\phi(x,y,0) = 0$$ and $$\Psi(x_0,y,\phi) = f(y)$$ where $x_0$ is a positive constant. I hope the number of boundary conditions are sufficient to make the solution unique. For 1D ODE I can use 4th order Runge-Kutta for instance, I wonder if a similar technique exists for multivariable PDE.

In the case of ordinary differential equations (ODE), Runge-Kutta methods are appreciated because of their generic nature. However, they are not always a suitable choice. For instance, in the case of stiff systems, it is better to use an adaptative method, such as Rosenbrock methods or backward differentiation formulas (BFD). See e.g. the book [1] for details.

In the case of partial differential equations (PDE), there is no such generic method. The overview given in chapter 20 of [2] states that partial differential equations are classified into three categories, hyperbolic, parabolic, and elliptic, on the basis of their characteristics (curves of information propagation). From a computational point of view, there is also an essential distinction between time-evolution problems and time-independent (static) problems, which gives rise to different numerical methods. Nevertheless, some numerical methods can be adapted to many different PDE problems, in particular the finite-difference method (see e.g. [2, 3]), and the finite element method. For some problems, other methods such as the finite volume method and the discontinuous Galerkin method can be used too.

To go back to the comparison with Runge-Kutta methods, let us consider the case of finite-difference methods (FDM) for PDEs. As we will see with two time-evolution PDEs, the efficiency of a FDM is highly problem-dependent, and such a derivation must be carried out carefully.

• Let us consider the one-dimensional heat equation $u_t = \alpha u_{xx}$. An explicit FD method with centered differentiation in space writes $$\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_{i}^n + u_{i-1}^n}{\Delta x^2} \, ,$$ i.e. $$u_i^{n+1} = \left(1 - 2 \frac{\alpha\Delta t}{\Delta x^2}\right) u_i^n + \frac{\alpha\Delta t}{\Delta x^2} (u_{i+1}^n + u_{i-1}^n) \, .$$ This method is stable and convergent for $\alpha\Delta t/\Delta x^2 \leq 1/2$.

• Now, we consider the first-order linear advection equation $u_t + c u_x = 0$ and we proceed similarly. An explicit FD method with centered differentiation in space writes $$\frac{u_i^{n+1} - u_i^n}{\Delta t} + c \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta x} = 0 \, ,$$ i.e. $$u_i^{n+1} = u_i^n - \frac{c\Delta t}{2\Delta x} (u_{i+1}^n - u_{i-1}^n) \, .$$ Despite its quite natural derivation, this method suffers from severe stability problems and is useless in practice [4].

As shown on these two simple examples, the derivation of convergent FD methods is not generic. Depending on the problem, an approach may be successful, or not.

[1] J.C. Butcher: Numerical Methods for Ordinary Differential Equations, 2nd ed., John Wiley & Sons, 2008.

[2] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery: Numerical Recipes. The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007.

[3] A.R. Mitchell, D.F. Griffiths: The Finite Difference Method in Partial Differential Equations, John Wiley & Sons, 1980.

[4] R.J. LeVeque: Numerical Methods for Conservation Laws, Birkhäuser, 1992.

• Thank you for your insightful answer. It seems like FDM still has too much uncertainty around the stability issue. May be for the first try I will try FEM. – nougako Sep 17 '17 at 14:13
• @nougako In case of a generic approach, I would suggest to do so too. If your first FEM trial is unstable, you can still tune the interpolating functions, the integration algorithm, and the mesh. For info, the Comsol software has a build-in FEM tool for custom PDEs. – Harry49 Sep 17 '17 at 14:20

You can use finite differences, finite elements or finite columns, all of them are based on discretizing the domain, which is the same idea behind RK. But I guess finite differences is intuitively more similar to RK

• Thank you. Can you recommend me some references (other than wikipedia)? The closest I have worked with is the finite difference method for 1D where the ODE is transformed to a matrix equation but this method is for eigenvalue problem. – nougako Sep 17 '17 at 13:11
• @nougako you're right. Finite elements is perhaps the most suitable method for your problem this is a good reference – caverac Sep 17 '17 at 13:17