NumMethodsPDEs

Numerical Methods for PDEs: Solving PDEs on a computer

=Introduction=

In this tutorial, we'll take a look at how we might model aspects of the world around us on a computer. When we use mathematics to describe many of the phenomena that we see we end up using Partial Differential Equations (PDEs), and so we seek ways to solve these numerically.

At the outset, let's list some take home-messages from the material below:
 * Finding a numerical solution to a PDE often boils down to solving a system of linear equations, of the (matrix) form $$Ax=b$$.
 * 3rd party linear algebra libraries will be much faster and less error prone than writing your own code for this task. Many of these libraries will also execute in a parallel fashion.
 * The matrix A is often sparse and so iterative methods will often be much faster than performing a direct method for solution (such as an LU decomposition).

=Looking at Nature (& a quick philosophical aside)=

When we observe nature, certain patterns crop up again and again. For example, in the the field of electrostatics, if a function $$f$$ describes a distribution of electric charge, then Poisson's equation:


 * $${\nabla}^2 \varphi = f.$$

gives the electric potential $$\varphi$$.

However, Poisson's equation also describes the steady state temperature of a material when subjected to some heating and also crops up when considering gravitational potentials..

What's going on here? Why the recycling of the same equation? Is Poisson's equation fundamental in some way? Well, yes I suppose it is. However, we can look at it another way and say that Poisson's equation is the way we describe--using the language of maths--steady state phenomena which involve potentials. OK, this sounds a bit out there, but consider the following very general equation for a moment (it's a second-order linear equation in two variables):

$$Au_{xx} + Bu_{xy} + Cu_{yy} + Du_x + Eu_y + Fu = G$$

It turns out that we can categorise certain instances of this equation and then relate these categories to the kind of phenomena that they describe:


 * Parabolic: $$B^2 - 4AC = 0$$.  This family of equations describe heat flow and diffusion processes.
 * Hyperbolic: $$B^2 - 4AC > 0$$.  Describe vibrating systems and wave motion.
 * Elliptic: $$B^2 - 4AC < 0$$.  Steady-state phenomena.

That's pretty handy!

=The Finite Difference Method=

In this section we'll look at solving PDEs using the Finite Difference Method. Finite difference equations have a long and illustrious past. Indeed, they were used by Newton and the gang as a logical route in to talk about differential calculus. So, in fact, we're just going to work back in the other direction..

Discritisation
The functions and variables that we see in, for instance, Poisson's equation are continuous; meaning that they could have an infinite number of values. This is fine when solving equations analytically, but we are limited to considering only a finite number of values when solving PDEs with digital computer. We need a way to discritise our problems. Happily, help is at hand from the Taylor series. Recall that we can approximate the value of a function at a point distant (displaced by an amount h) from x using a linear summation of derivatives of the function at x:


 * $$f(x+h) \approx f(x) + \frac{f'(x)}{1!}h + \frac{f''(x)}{2!}h^2 + \cdots$$

Let's truncate this series after the first differential and rearrange, to get ourselves the forward-difference approximation:


 * $$f'(x) \approx \frac{f(x+h) - f(x)}{h} $$

Thus, following the finite difference approach, we can discritise our equations by replacing our differential operators by their difference approximations, as helpfully provided by a truncated Taylor series. If, like me, you like thinking visually, the schematic below shows the relationship between a differential operator--the gradient at a point $$x_1$$ and the finite difference approximation to the gradient:



If we were to consider a displacement of -h instead, we would get the backward-difference approximation to $$f'(x)$$:


 * $$f'(x) \approx \frac{f(x) - f(x-h)}{h} $$,

Subtracting the backward from the forward we can generate the central-difference approximation:


 * $$f'(x) \approx \frac{1}{2h}(f(x+h) - f(x-h))$$.



By retaining an extra term from the Taylor series, and some similar rearrangements, we can also formulate a finite difference approximation to the 2nd differential. For example the central-difference formula for this is:


 * $$f''(x) \approx \frac{1}{h^2} [ f(x+h) - 2f(x) + f(x-h) ] $$.

In a similar way we can derive formulae for approximations to partial derivatives. Armed with all these formulae, we can convert differential equations into difference equations, which we can solve on a computer. Since we are dealing with approximations to the original equation(s), we must be vigilant that errors in our input and intermediate calculations do not accumulate such that our final solution is meaningless.

Eliptical Equations: Discrete Laplacian Operator
Consider Laplace's equation in 2-dimensions:


 * $$\frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} = 0$$

This describes purely diffusive phenomena and can give us the steady state temperature over a 2D plate, three sides of which are held at 100°C, while the fourth is held at 0°C. We have divided the plate into a number of grid cells:



Using the central difference approximation to a 2nd differential that we derived in the last section, we can rewrite the differential equation as a difference equation:


 * $$\frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} \approx \frac{1}{h^2}(u^j_{i+1} -2u^j_i +u^j_{i-1}) + \frac{1}{k^2}(u^{j+1}_i -2u^j_i +u^{j-1}_i)$$

where h is the grid spacing in the i-directions and k is the spacing the j-direction. If we set $$h=k$$, then:


 * $$u^j_{i+1} + u^j_{i-1} + u^{j+1}_i + u^{j-1}_i - 4u^j_i = 0$$

solving for $$u^j_i$$:


 * $$u^j_i = \frac{1}{4}(u^j_{i+1} + u^j_{i-1} + u^{j+1}_i + u^{j-1}_i)$$

Note that this last equation says--quite intuitively--that the value at the central grid point is equal to the average of its four neighbours. We can draw a picture of the stencil needed to compute the value:



Placing the stencil on the 4 grid cells on the inner region (orange) in turn, we generate the following 4 difference equations:



\begin{alignat}{5} -4u_{22} & + u_{32} & + u_{23} & + 100 & + 100 & = 0 \\ -4u_{32} & + 0 & + u_{33} & + u_{22} & + 100 & = 0 \\ -4u_{23} & + u_{33} & + 100 & + 100 & + u_{22} & = 0 \\ -4u_{33} & + 0 & + 100 & + u_{23} & + u_{32} & = 0 \\ \end{alignat} $$

This system of equations can be written in matrix form:

\begin{bmatrix}-4 & 1 & 1 & 0\\1 & -4 & 0 & 1\\1 & 0 & -4 & 1\\0 & 1 & 1 & -4\end{bmatrix} \begin{bmatrix}u_{22}\\u_{32}\\u_{23}\\u_{33}\end{bmatrix} = \begin{bmatrix}-200\\-100\\-200\\-100\end{bmatrix} $$

Which may be summarised as the vector equation, Ax=b, which can be efficiently & conveniently solved using a library linear algebra routines, such as LAPACK, PETSc or PLASMA.



Parabolic Equations: Forward (Explicit) Euler
One example of a parabolic equation is the 1D heat equation, which describes the evolution of temperature at various points in a thin, insulated bar over time. Imagine such a bar where the ends are kept at 0°C:



Since this scenario involves the passage of time, just stating the boundary conditions is not enough, and we must also specify some initial condition for the temperature along the bar's length. This is shown below for a bar divided up into grid cells, as before:



The 1D heat equation is:


 * $$\frac{\partial u}{\partial t}=\alpha \frac{\partial^2u}{\partial x^2} $$

where $$\alpha$$--the thermal diffusivity--is set to 1 for simplicity. One way to rewrite this as a difference equation (using a forward difference for the LHS and a central difference for the RHS) is:


 * $$\frac{u^{j+1}_i - u^j_i }{k} = \frac{u^j_{i+1} - 2u^j_i + u^j_{i-1}}{h^2}$$

and rearrange to solve for $$u^{j+1}_i$$:


 * $$u^{j+1}_i = u^j_i + r(u^j_{i+1} - 2u^j_i + u^j_{i-1})$$


 * $$u^{j+1}_i = ru^j_{i+1} + (1- 2r)u^j_i + ru^j_{i-1}$$

where:

$$r=\frac{k}{h^2}$$

This time, we see that the appropriate stencil is:



If we set $$h=\frac{1}{10}$$ and $$k=\frac{1}{1000}$$, then $$r=\frac{1}{10}$$, and the equation becomes:


 * $$u^{j+1}_i = u^j_i + \frac{1}{10}(u^j_{i+1} - 2u^j_i + u^j_{i-1})$$


 * $$u^{j+1}_i = \frac{1}{10}(u^j_{i+1} + 10u^j_i - 2u^j_i + u^j_{i-1})$$


 * $$u^{j+1}_i = \frac{1}{10}(u^j_{i+1} + 8u^j_i + u^j_{i-1})$$

Which we can illustrate using weightings on the stencil nodes:



Forward Euler is as an explicit method, as one unknown value ($$u^{j+1}_i$$) is expressed directly in terms of known values (all values of u at time j are known).

We could implement a simple marching algorithm, where we centre the stencil over a set of known values at time j, and may calculate the unknown temperature for grid cell i at time j+1. However, it would be better to group together all the calculations for a particular time step into a system of simultaneous equations (for our example 1d heat problem):


 * $$u_{22} = \frac{1}{10}(0) + \frac{8}{10}(0.2) + \frac{1}{10}(0.4)$$


 * $$u_{23} = \frac{1}{10}(0.2) + \frac{8}{10}(0.4) + \frac{1}{10}(0.6)$$


 * $$u_{24} = \frac{1}{10}(0.4) + \frac{8}{10}(0.6) + \frac{1}{10}(0.8)$$


 * $$\cdots$$

and then re-write in matrix form:



\begin{bmatrix}u_{22} \\ u_{23} \\ u_{24} \\\vdots \end{bmatrix} = \begin{bmatrix} 0.1   & 0.8 & 0.1    & 0   & 0   & \cdots & 0 \\ 0     & 0.1 & 0.8    & 0.1 & 0   & \cdots & 0 \\ \vdots &    &        &     &     &        & \vdots \\ &    & \cdots & 0   & 0.1 & 0.8    & 0.1 \\ \end{bmatrix} \times \begin{bmatrix}0 \\ 0.2 \\ 0.4 \\ \vdots \\ 0.2 \\ 0 \end{bmatrix} $$

as the computational efficiency of a matrix vector multiplication will be hard to better. Note that we can write the above piece of linear algebra in a general form:



\begin{bmatrix}u_{22} \\ u_{23} \\ u_{24} \\ \vdots \end{bmatrix} = \begin{bmatrix} r    & 1-2r & r      & 0   & 0   & \cdots & 0 \\ 0     & r   & 1-2r   & r   & 0   & \cdots & 0 \\ \vdots &    &        &     &     &        & \vdots \\ &    & \cdots & 0   & r   & 1-2r   & r \\ \end{bmatrix} \times \begin{bmatrix}u_{11} \\ u_{12} \\ u_{13} \\ u_{14} \\ \vdots \end{bmatrix} $$

Solution Stability and the CFL Condition
The Forward Euler scheme has the benefit of simplicity, but is unstable for large time steps. The CFL (Courant–Friedrichs–Lewy) condition states that the following must hold:


 * $$r=\frac{k}{h^2} \le 0.5 $$

for stability.

Parabolic Equations: Backward (Implicit) Euler
An alternative way to formulate our finite difference is to use a backward-difference. The backward Euler scheme requires more computation per time-step than the forward scheme, but has the distinct advantage that it is unconditionally stable with regard to the relative size between the spatial and time step sizes; i.e. we are not limited to $$r \le 0.5$$ when we use the backward Euler scheme.

We will use a central-difference on the RHS, as before, but shifting the time index to $$j+1$$. This will allow us to formulate a backward-difference on the RHS between the time indices $$j$$ and $$j+1$$:


 * $$\frac{u^{j+1}_i - u^j_i }{k} = \frac{u^{j+1}_{i+1} - 2u^{j+1}_i + u^{j+1}_{i-1}}{h^2}$$


 * $$u^{j+1}_i - u^j_i = \frac{k}{h^2}(u^{j+1}_{i+1} - 2u^{j+1}_i + u^{j+1}_{i-1})$$


 * $$u^{j+1}_i - u^j_i = r(u^{j+1}_{i+1} - 2u^{j+1}_i + u^{j+1}_{i-1})$$


 * $$u^{j+1}_i - u^j_i = ru^{j+1}_{i+1} - 2ru^{j+1}_i + ru^{j+1}_{i-1}$$


 * $$u^{j+1}_i - ru^{j+1}_{i+1} + 2ru^{j+1}_i - ru^{j+1}_{i-1} = u^j_i$$


 * $$u^j_i = - ru^{j+1}_{i+1} + (1 + 2r)u^{j+1}_i - ru^{j+1}_{i-1}$$

If we again choose $$h=\frac{1}{10}$$ and $$k=\frac{1}{1000}$$, such that $$r=\frac{k}{h^2}=\frac{1}{10}$$, as before, we end up with the following system of equations when we consider all values of $$i$$ for our first time-step--from $$j$$ to $$j+1$$--where we know values of all the cells at time $$j$$:


 * $$0.2 = -\frac{1}{10}(0) + \frac{12}{10}(u_{22}) -\frac{1}{10}(u_{23})$$


 * $$0.4 = -\frac{1}{10}(u_{22}) + \frac{12}{10}(u_{23}) -\frac{1}{10}(u_{24})$$


 * $$0.6 = -\frac{1}{10}(u_{23}) + \frac{12}{10}(u_{24}) -\frac{1}{10}(u_{25})$$
 * $$\cdots$$

We can see that this is akin to looping over the domain with the following stencil:



Rearranging, so that we have only constants on the RHS:


 * $$\frac{1}{10}(0) + \frac{12}{10}(u_{22}) -\frac{1}{10}(u_{23}) = 0.2 + \frac{1}{10}(0)$$


 * $$\frac{1}{10}(u_{22}) + \frac{12}{10}(u_{23}) -\frac{1}{10}(u_{24}) = 0.4$$


 * $$\frac{1}{10}(u_{23}) + \frac{12}{10}(u_{24}) -\frac{1}{10}(u_{25}) = 0.6$$
 * $$\cdots$$

We can re-write this system of equations in matrix form:



\begin{bmatrix} 1.2   & -0.1 & 0      & 0   & \cdots & 0 \\ -0.1  & 1.2  & -0.1   & 0   & \cdots & 0 \\ \vdots &     &        &     &        & \vdots \\ &     & \cdots & 0   & -0.1   & 1.2 \\ \end{bmatrix} \times \begin{bmatrix} u_{22} \\ u_{23} \\ u_{24} \\\vdots \end{bmatrix} = \begin{bmatrix} 0.2 + 0.1(0) \\ 0.4 \\ \vdots \\ 0.2 +0.1(0) \end{bmatrix} $$

The general form is:



\begin{bmatrix} 1+2r  & -r   & 0      & 0   & \cdots & 0 \\ -r    & 1+2r & -r     & 0   & \cdots & 0 \\ \vdots &     &        &     &        & \vdots \\ &     & \cdots & 0   & -r     & 1+2r \\ \end{bmatrix} \times \begin{bmatrix} u_{22} \\ u_{23} \\ u_{24} \\\vdots \end{bmatrix} = \begin{bmatrix} u_{12} + r(u_{11}) \\ u_{13} \\ \vdots \\ u_{1n-1} + r(u_{1n}) \end{bmatrix} $$

Parabolic Equations: Crank-Nicolson
A third scheme for solving the 1D heat equation is to create a finite difference equation using the Crank-Nicolson stencil:



The stencil amounts to an average of both the forward and backward Euler methods. The resulting difference equation is:


 * $$\frac{u^{j+1}_i - u^j_i}{k} = \frac{1}{2}(\frac{u^{j+1}_{i+1} - 2u^{j+1}_{i} + u^{j+1}_{i-1}}{h^2} + \frac{u^j_{i+1} - 2u^j_i + u^j_{i-1}}{h^2}).$$

Which we may rearrange to become:


 * $$u^{j+1}_i - u^j_i = \frac{1}{2}\frac{k}{h^2}(u^{j+1}_{i+1} - 2u^{j+1}_i + u^{j+1}_{i-1} + u^j_{i+1} - 2u^j_i + u^j_{i-1})$$


 * $$2(u^{j+1}_i - u^j_i) = r(u^{j+1}_{i+1} - 2u^{j+1}_i + u^{j+1}_{i-1} + u^j_{i+1} - 2u^j_i + u^j_{i-1})$$,

where $$r=\frac{k}{h^2}$$.

Collecting the unknowns (values of u at the next timestep, $$j+1$$) to the LHS and the knowns (values of u at current timestep, $$j$$), we get:


 * $$2u^{j+1}_i + 2ru^{j+1}_i - ru^{j+1}_{i+1} - ru^{j+1}_{i-1} = 2u^j_i +ru^j_{i+1} -2ru^j_i + ru^j_{i-1}$$


 * $$(2+2r)u^{j+1}_i - ru^{j+1}_{i+1} - ru^{j+1}_{i-1} = (2-2r)u^j_i +ru^j_{i+1} + ru^j_{i-1}$$.

If we choose $$h=1/10$$ and $$k=1/100$$ (note that the time-step is 10x larger than than chosen for the forward Euler example), then $$r=1$$, and we get:


 * $$-u_{i-1}^{j+1} +4u_i^{j+1} -u_{i+1}^{j+1} = u_{i-1}^j + u_{i+1}^j$$.

We can again create a system of linear equations for all the cells j along the bar at a given time step i:


 * $$-0 + 4u_{22} - u_{23} = 0 + 0.4$$


 * $$-u_{22} + 4u_{23} - u_{24} = 0.2 + 0.6$$


 * $$-u_{23} + 4u_{24} - u_{25} = 0.4 + 0.8$$


 * $$\cdots$$

We have $$j$$ unknowns and $$j$$ equations. Then we can again translate into matrix notation and solve using a library linear algebra routines:



\begin{bmatrix} -4    & 1      & 0       & 1  & 0  & 0  & \cdots & 0 \\ 1    & -4     & 1       & 0  & 1  & 0  & \cdots & 0 \\ \vdots &       &         &    &    &    &        & \vdots \\ & \cdots &        & 0  & 1  & 0  & 1      & -4\\ \end{bmatrix} \times \begin{bmatrix} u_{22} \\ u_{23} \\ u_{24} \\\vdots \end{bmatrix} = \begin{bmatrix} 0 + 0.4 \\ 0.2 + 0.6 \\ \vdots \\ 0 + 0.4 \end{bmatrix} $$

In general:



\begin{bmatrix} -(2+2r) & r      & 0       & r  & 0  & 0  & \cdots & 0 \\ r     & -(2+2r) & r       & 0  & r  & 0  & \cdots & 0 \\ \vdots &         &         &    &    &    &        & \vdots \\ & \cdots &         & 0  & r  & 0  &      r & -(2+2r)\\ \end{bmatrix} \times \begin{bmatrix} u_{22} \\ u_{23} \\ u_{24} \\\vdots \end{bmatrix} = \begin{bmatrix} ru_{11} + (2-2r)u_{12} + ru_{13} \\ \vdots \end{bmatrix} $$

=Other Numerical Methods=

We have looked at the finite difference method (FDM) above, but other (more mathsy) approaches exist, such as the finite element and finite volume methods. Each approach has its pros and cons. For example, the FDM is limited to regular grids spread over rectangular domains, but is easy to program. In contrast the FEM is applicable to intricately shaped domains and the resolution of the individual elements can be varied, but is requires more involved programming and a deeper level of understand of the relevant mathematics. It should also be noted, however, that under certain circumstances, the FEM and FVM boil down to be exactly the same as the FDM, even though we've approached the task of solving a PDE from a different direction.

One thing which is always common to the approaches is that--the problem boils down to a system of linear equations to solve, and that a library of linear algebra routines will be the workhorse when finding solutions to PDEs on a computer. Repeat after me, "a fast linear algebra library is key!"

The Finite Element Method (FEM)
The general approach here is to create a potential function from our PDE and to find a solution to it via minimising the potential function. Key steps include:


 * Discritisation of the original PDE, often by proposing that we can approximate the solution using a linear sum of (polynomial) basis functions.
 * Create the potential function by invoking the method of weighted residuals, and in particular the Galerkin specialisation of this approach.
 * Use whatever algebraic manipulation and integration (possibly numerical) is required to create a system of linear equations.
 * Solve said system and in doing so, solve the original PDE.

For example, consider, for example, Poisson's equation over a 2-dimensional region, $$R$$:


 * $$\nabla^2u=f(x,y)$$

where $$u=0$$ at the boundary of the domain.

Following the method of weighted residuals, we multiply both sides by weighting functions $$\phi_i$$, and take integrals over the domain:


 * $$ \iint\limits_R \nabla^2u\phi_i\, dA = \iint\limits_R f\phi_i\, dA $$

We can apply integration-by-parts to the LHS ($$\phi_i\nabla^2u$$), for example:


 * $$\int \phi_i\nabla^2 u = \phi_i\nabla u - \int \nabla\phi_i \cdot \nabla u$$

Differentiating both sides:


 * $$\phi_i\nabla^2 u = \nabla \cdot (\phi_i\nabla u) - \nabla\phi_i \cdot \nabla u$$

So inserting this into our equation for weighted residuals, above:


 * $$ \iint\limits_R [\nabla \cdot (\phi_i\nabla u) - \nabla\phi_i \cdot \nabla u]\, dA = \iint\limits_R f\phi_i\, dA $$

and rearranging, we get:


 * $$ \iint\limits_R \nabla\phi_i \cdot \nabla u\, dA = \iint\limits_R \nabla \cdot (\phi_i\nabla u)\, dA - \iint\limits_R f\phi_i\, dA $$

We can see that the 1st term on the RHS is a divergence(the divergence of a vector field $$F$$, $$div F = \nabla \cdot F$$), and so we can apply Gauss' Divergence Theorem which states, for the region $$R$$ bounded by the curve s:


 * $$\iint\limits_R \nabla \cdot F\, dA = \oint\limits_s F \cdot n\, ds $$

where $$n$$ is the unit vector normal to the surface.

and so,


 * $$ \iint\limits_R \nabla\phi_i \cdot \nabla u\, dA =  \oint\limits_s \phi_i\nabla u \cdot n\, ds - \iint\limits_R f\phi_i\, dA $$

but our homogeneous boundary conditions also tell us that $$u$$ is equal to zero at the boundary and so the 1st term on the RHS goes to zero:


 * $$ \iint\limits_R \nabla\phi_i \cdot \nabla u\, dA = - \iint\limits_R f\phi_i\, dA $$

So far, all of the operations above have been done on our unaltered, continuous differential equation. At this point, let's introduce the discritisation--and hence the approximation--that we need to make to work on a digital computer. We intend to approximate the solution $$u(x,y)$$ using a linear combination of basis functions:


 * $$U(x,y) = \sum_{j=1}^m U_j\phi_j(x,y)$$

The Galerkin method is being used here, as our basis functions are the same as the weighting functions used above.


 * $$ \sum_{j=1}^m U_j \iint\limits_R \nabla\phi_i \cdot \nabla\phi_j \, dA = - \iint\limits_R f\phi_i\, dA$$

At this point, we can introduce a symmetric 'stiffness matrix' with entries:


 * $$K_{ij} = \iint\limits_R \nabla\phi_i \cdot \nabla\phi_j \, dA$$

and a vector F, with entries:


 * $$F_i = - \iint\limits_R f\phi_i\, dA$$

and thus our $$m$$ equations:


 * $$\sum_{j=1}^m K_{ij}U_j = F_i $$

which, we see, yet again is a linear system of equations of the form $$Ax=b$$. Cue a call to a linear algebra library!



FEM packages include:


 * Elmer: http://www.csc.fi/english/pages/elmer