Solving the Diffusion Equation Explicitly

This post is part of a series of **Finite Difference Method Articles**. Other posts in the series concentrate on Derivative Approximation, the Crank-Nicolson Implicit Method and the Tridiagonal Matrix Solver/Thomas Algorithm:

- Derivative Approximation via Finite Difference Methods
*Solving the Diffusion Equation Explicitly*- Crank-Nicolson Implicit Scheme
- Tridiagonal Matrix Solver via Thomas Algorithm

In Part 1 of the series on Finite Difference Methods it was shown that continuous derivatives could be approximated and applied to a discrete domain. The next step is to apply these derivatives to a parabolic PDE. The heat equation is the canonical example of a parabolic PDE and this will now be discretised.

The heat equation with initial condition \( g \) is given below by:

\[ \frac{\partial f}{\partial t} = \frac{\partial^2 f}{\partial x^2}, \qquad f(x, 0) = g(x) \]

This is discretised by applying a *forward difference* to the time derivative and a *centered second difference* for the diffusion term to give:

\[ \frac{f^{n+1}_i - f^n_i}{\Delta t} = \frac{f^n_{i+1}- 2f^n_i + f^n_{i-1}}{(\Delta x)^2} \]

This equation can be rearranged for \( f^{n+1}_i \) to give:

\[ f^{n+1}_i = f^n_i + \frac{\Delta t}{(\Delta x)^2} \left( f^n_{i+1} - 2f^n_i + f^n_{i-1} \right) \]

For \(n=1\) all of the approximations to the solution \(f\) are known on the right hand side of the equation. This makes the equation *explicit*. Thus the only remaining task is to determine \(\Delta t\), which is set to \( (\Delta x)^2 \). Thus the equation reduces to:

\[ f^{n+1}_i = f^n_{i+1} - f^n_i + f^n_{i-1} \]

This formula will allow calculation of \(f^2_i\) for all \(i \in {1,...,I} \) as an approximation to \(g\) provides the discretised initial condition. Let \(g_i\) be such an approximation, with three 'heat bumps' and take \(I=13\), then the values of \(g_i\) are given by:

\[ \begin{matrix}

0 & 0 & 1 & 0 & 0 & 1 & 2 & 1 & 0 & 0 & 1 & 0 & 0

\end{matrix} \]

This will allow a solution to be generated by progressing from the initial condition in a temporal manner. This is referred to as a *time marched* solution:

\[ \begin{matrix}

n = 1 & 0 & 0 & 1 & 0 & 0 & 1 & 2 & 1 & 0 & 0 & 1 & 0 & 0 \\

n = 2 & 0 & 1 & -1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & -1 & 1 & 0 \\

n = 3 & 1 & -2 & 3 & -1 & 1 & 0 & 2 & 0 & 1 & -1 & 3 & -2 & 1 \\

n = 4 & -3 & 6 & -6 & 5 & -2 & 3 & -2 & 3 & -2 & 5 & -6 & 6 & -3

\end{matrix} \]

For \(n=4\) the discretised solution is wildly oscillating and bears no resemblance to the continuous solution, which has the physical interpretation of diffusing the heat bumps. Clearly one of the assumptions made about the formulae is incorrect.

It turns out that the issue is with the choice of time step size \( \Delta t \) relative to the spatial step size \( \Delta x \). The correct stability criterion is given by:

\[ \frac{\Delta t}{(\Delta x)^2} \leq \frac{1}{2} \]

It is beyond the scope of this article to prove this result, but note that it puts a significant restriction on the time step - a problem which will be remedied in the next article.

Applying this new stability criterion, with \( \Delta t = 0.5 \times (\Delta x)^2 \), to the above formulae produces:

\[ f^{n+1}_i = \frac{1}{2} f^n_{i+1} + \frac{1}{2} f^n_{i-1} \]

The calculation for this formulae provides:

\[ \begin{matrix}

n = 1 & 0 & 0 & 1 & 0 & 0 & 1 & 2 & 1 & 0 & 0 & 1 & 0 & 0 \\

n = 2 & 0 & 0.5 & 0 & 0.5 & 0.5 & 1 & 1 & 1 & 0.5 & 0.5 & 0 & 0.5 & 0 \\

n = 3 & 0.25 & 0 & 0.5 & 0.25 & 0.75 & 0.75 & 1 & 0.75 & 0.75 & 0.25 & 0.5 & 0 & 0.25 \\

n = 4 & 0 & 0.375 & 0.125 & 0.625 & 0.5 & 0.875 & 0.75 & 0.875 & 0.5 & 0.625 & 0.125 & 0.375 & 0

\end{matrix} \]

This is an accurate approximation to the true solution. There are no negative values and the physical interpretation of the heat diffusing through a 1D bar fits with the solution.

The choice of time step is very restrictive. In practice, with a refined spatial mesh, the time step would be proihibitively small. The next article will introduce another method to overcome the small timestep issue.

comments powered by DisqusYou'll get instant access to a free 10-part email course packed with hints and tips to help you get started in quantitative trading!

Every week I'll send you a wrap of all activity on QuantStart so you'll never miss a post again.

Real, actionable quant trading tips with no nonsense.