This post is part of a series of Finite Difference Method Articles. Other posts in the series concentrate on Derivative Approximation, Solving the Diffusion Equation Explicitly and the Tridiagonal Matrix Solver/Thomas Algorithm:
In the previous tutorial on Finite Difference Methods it was shown that the explicit method of numerically solving the heat equation lead to an extremely restrictive time step. This motivates another scheme which allows for larger time steps, but with the trade off of more computational work per step. This method is known as the Crank-Nicolson scheme.
The explicit method for the heat-equation involved a forward difference term for the time derivative and a centred second derivative for the second space derivative:
\[ \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} \]
An Implicit Scheme
The Crank-Nicolson scheme modifies this to incorporate a weighted average of the second spatial step at time \(n\) and time \(n+1\). An obvious response is that the values of \(f\) are not known at \(n+1\) and questions arise over how they are calculated. These questions will be answered below. One final question occurs over how to split the weighting of the two second derivatives. The Crank-Nicolson scheme uses a 50-50 split, but others are possible. Stability is a concern here with \(\frac{1}{2} \leq \theta \le 1\) where \(\theta\) is the weighting factor.
The Crank-Nicolson scheme for the 1D heat equation is given below by:
\[ \frac{f^{n+1}_i - f^n_i}{\Delta t} = \frac{f^n_{i+1}- 2f^n_i + f^n_{i-1}}{2(\Delta x)^2} + \frac{f^{n+1}_{i+1}- 2f^{n+1}_i + f^{n+1}_{i-1}}{2(\Delta x)^2}\]
Letting \(r=\frac{\Delta t}{(\Delta x)^2}\), this equation can be rearranged to group the known and unknown terms separately:
\[ -\frac{r}{2} f^{n+1}_{i+1} + (1+r)f^{n+1}_i -\frac{r}{2} f^{n+1}_{i-1} = \frac{r}{2} f^{n}_{i+1} + (1-r)f^{n}_i + \frac{r}{2} f^{n}_{i-1} \]
Since there three unkown terms on the left hand side, it indicates that for each time step a set of simultaneous linear equations must be solved. For any particular time step \(n\) there will be \(I-2\) equations to solve. These correspond to each of the interior spatial grid points. Since we are using Dirichlet boundary conditions (i.e. fixed), the end points do not require such a linear system to be solved.
The linear system can be represented by a matrix equation \(AF=D\), where \(a=-\frac{r}{2}\), \(b=(1+r)\), \(c=-\frac{r}{2}\), \(d_i=\frac{r}{2} f^{n}_{i+1} + (1-r)f^{n}_i + \frac{r}{2} f^{n}_{i-1}\) and \(k\) runs from \(2\) until \(I-1\):
\[ \begin{bmatrix}
b & c & 0 & 0 & ... & 0 \\
a & b & c & 0 & ... & 0 \\
0 & a & b & c & 0 & 0 \\
. & . & & & & . \\
. & . & & & & . \\
. & . & & & & c \\
0 & 0 & 0 & 0 & a & b \\
\end{bmatrix} \begin{bmatrix}
f_1 \\
f_2 \\
f_3 \\
.\\
.\\
.\\
f_k \\
\end{bmatrix} = \begin{bmatrix}
d_1 \\
d_2 \\
d_3 \\
.\\
.\\
.\\
d_k \\
\end{bmatrix} \]
Solving the Linear System
A first attempt at solving this system might involve inverting the matrix \(A\) to give the values of the column matrix \(F\) as \(F=A^{-1}D\). However this is far from the most optimal method and does not take into account the sparsity or the banded structure of the matrix \(A\). In the next tutorial it will be shown how this matrix can be solved in an efficient manner so that the work done at each time step is minimised.
