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) \]

Computation of Solution

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.

Stability

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.

Conclusion

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.