During my final year at Warwick University I did a project on Numerical Weather Forecasting and one of the methods that came up was the method of finite differences. I recently came across finite differences again over the last few weeks and how they can be used to solve some partial differential equations (PDEs). Here is an example that I found in one of my books that I decided to play around with

Solve, using the method of finite differences, the PDE $$x\dfrac{\partial{f}}{\partial{x}}+(y+1)\dfrac{\partial{f}}{\partial{y}}=0$$ for $0\leq x,y \leq 1$ with $$f(x,0)=x-1$$ $$f(x,1)=\dfrac{x-2}{2}$$ $$f(0,y)=-1$$ $$f(1,y)=-\dfrac{y}{y+1}$$

To solve this problem I used the central difference equations given by

$$\dfrac{\partial{f}}{\partial{x}}\approx \dfrac{f(x+h,y)-f(x-h,y)}{2h}$$

$$\dfrac{\partial{f}}{\partial{y}}\approx \dfrac{f(x,y+k)-f(x,y-k)}{2k}$$

with a step length in both the $x$ and $y$ directions of $\frac{1}{3}$. The domain can be overlaid with a grid as shown in the diagram

Domain with overlaid gridThis makes the whole problem easier to deal with because now we can see visually what is happening rather than just dealing with everything purely algebraically – there’s no need to make things more difficult than they need to be after all.

The idea is to form a system of four simultaneous equations with $A$, $B$, $C$, and $D$ as the unknown quantities. When the system is solved then the values obtained correspond to the approximate values of $f$ at each of the grid points.┬áHere is my full solution to this problem – PDE Solution Using Finite Differences. This fills in some of the gaps in our knowledge about the function $f$ but there is still a lot of information missing – in particular the points between the grid points.

To resolve this we can make the step-lengths smaller and create more interior grid points. This gives us more information but at a price – for example with a step length of $\frac{1}{5}$ in both directions we would have a system of $16$ simultaneous equations in $16$ unknowns – I don’t really fancy sorting that mess out without the help of a computer. And what’s more – as the step-length gets even smaller we get more and more complicated systems with more unknowns so eventually we would reach a point where even a computer would struggle to do anything with the information.

Nevertheless – finite difference methods are used and are useful for all kinds of problems and can provide a very visual way of numerically solving some otherwise very difficult equations.


There are many examples of differential equations that cannot be solved analytically – in fact, it is very rare for a differential equation to have an explicit solution. Euler’s Method is a way of numerically solving differential equations that are difficult or that can’t be solved analytically. Euler’s method, however, still has its limitations.

For a differential equation $y^{\prime}=f(x,y(x))$ with initial condition $y(x_{0})=y_{0}$ we can choose a step-length $h$ and approximate the solution to the differential equation by defining $x_{n}=x_{0}+nh$ and then for each $x_{n}$ finding a corresponding $y_{n}$ where $y_{n}=x_{n-1}+hf(x_{n-1},y_{n-1})$. This method works quite well in many cases and gives good approxiamtions to the actual solution to a differential equation, but there are some differential equations that are very sensitive to the choice of step-length $h$ as the following demonstrates.

Let’s look at the differential equation $y^{\prime}+110y=100$ with initial condition $y(0)=2$.

This differential equation has an exact solution given by $y=1+\mathrm{e}^{-100t}$ but this example is a very good example which demonstrates that Euler’s method cannot be used blindly. Let’s look at what happens for a few different step-lengths.

For the step-length $h=0.019$ step-length we get the following behaviour

Behaviour of numerical solution of an ODE

The red curve is the actual solution and the blue curve represents the behaviour of the numerical solution given by the Euler method – it is clear that the numerical solution converges to the actual solution so we should be very happy. However, look what happens when the step-length $h=0.021$ is chosen

Behaviour of numerical solution of an ODE

Again the actual solution is represented by the red line which on this diagram looks like a flat line because the blue curve gets bigger and bigger as you move along the $x$-axis. So a change of just $0.002$ in the step-length has completely changed the behaviour of the numerical solution. For a step-length $h=0.03$ the graph would look as follows

Behaviour of numerical solution of an ODEThe actual solution can barely be seen and the numerical solution gets out of control very quickly – this solution is completely useless – the scales on the $y$-axis are enormous and increasing the step-length only makes this worse. What has happened?

It can be shown by induction that for $n \in \mathbb{N}$ that $y_{n}=1+(1-100h)^{n}$. This converges only for $h<0.02$ and diverges for $h>0.02$. $h=0.02$ is a limiting case and gives an oscillating numerical solution that looks as follows

Behaviour of numerical solution of an ODE

For this particular example for $h<0.02$ and as the step-length gets closer to $0$ the solution will converge faster and for $h>0.02$ as the step-length increases the solution will diverge more rapidly.

So even though we have Euler’s method at our disposal for differential equations this example shows that care must be taken when dealing with numerical solutions because they may not always behave as you want them to. This differential equation is an example of a stiff equation – in other words, one that is very sensitive to the choice of step length. In general as the step-length increases the accuracy of the solution decreases but not all differential equations will be as sensitive to the step-length as this differential equation – but they do exist.

Some differential equations are relatively easy to solve – even though at A-Level this may seem to be the norm, in reality it is very rare for a differential equation to have an explicit solution. The types of differential equations at A-Level can be solved by direct integration or by using a range of nifty techniques like making a substitution to change the differential equation into one that can be solved directly, separation of the variables can be used to solve some equations of the form $f(y)y’=g(x)$; sometimes the integrating factor $\mathrm{e}^{\int P(x) \mathrm{d}x}$ can be used to solve equations of the form $y^{\prime}+P(x)y=Q(x)$ and second-order equations of the form $ay^{\prime \prime}+by’+cy=f(x)$ where $a,b$ and $c$ are constants, can be solved by finding the complementary function and particular integral and combining the two.

However, these are very special cases. Many differential equations that come up in real-life situations are very non-linear in nature. Non-linear differential equations are generally very difficult to solve, but quite often, it is the behaviour of the solutions that interests us more than the actual solution. We may not be really interested in a particular solution but a family of solutions so that we can easily compare their behaviours particularly for equations that evolve over time. In some cases – even though it may be very difficult or even impossible to find the solution to a differential equation analytically, we can still analyse the behaviour through use of a phase-diagram.

A phase-diagram is a vector field that we can use to visually present the solutions to a differential equation. For example here is a second-order differential equation – (this is an example that I did that appears in the book by D. W. Jordan and P. Smith titled Nonlinear Ordinary Differential Equations – An Introduction for Scientists and Engineers Fourth Edition)
$$ \ddot{x} = x-x^{2}$$

This second order-differential equation can be separated into a system of first-order differential equations given by $$\dot{x}=y$$ $$\dot{y}=x-x^{2}$$

This is a very common technique for solving differential equations since first-order ODEs are generally easier to solve than second-order ODEs. Note that $$\frac{\mathrm{d}y}{\mathrm{d}x}= \frac{\dot{y}}{\dot{x}}=\frac{x-x^{2}}{y}$$

The phase-portrait of this looks as follows

Phase Portrait for a second order differential equation

A phase portrait drawn in SAGE Math for a second-order ODE

At a point on the phase-diagram there is attached a vector which are represented here by the arrows – the vector tells you which direction to move in. So given a starting point it is simply a case of moving from point to point in the direction that the vector tells you to move in. So different starting points will give different paths through the phase-portrait called trajectories, some of which are shown in the following diagram

Phase Portrait for a second order differential equation

Phase portrait of a second-order ODE with trajectories

Each trajectory is described by the equation $3y^{2}=-2x^{3}+3x^{2}+C$ for different values of $C$. The trajectories above are for $C=-8$, $-4$, $-0.5$, $0$, $1$, $2$, $4$ and $6$. Here the horizontal $x$-axis represents the position of a particle and the vertical $y$-axis represents the velocity of the particle – remember that $y=\dot{x}$.

So now if we start in the top left hand corner of the diagram on the curve corresponding to when $C=8$ (the right-most blue line) the arrows tell us to move to the right and down – so our displacement increases as we move to the right but our velocity decreases as we move down. Eventually we cross the $x$-axis and here our velocity is instantaneously zero; then we start moving backwards since our velocity becomes negative as we continue moving downwards, and our displacement decreases as we move to the left. So I have not even calculated the actual solution to the ODE but I still know how the solution behaves.