6 The Crank-Nicolson method

In the notation established for the explicit method, the so-called Crank-Nicolson scheme can be written

$u_j^{n+1}=u_j^n+\frac{1}{2}r\left(\right.\underset{\text{†}}{\underset{\⏟}{u_{j-1}^n-2u_j^n+u_{j+1}^n}}+\underset{\text{‡}}{\underset{\⏟}{u_{j-1}^{n+1}-2u_j^{n+1}+u_{j+1}^{n+1}}}\left.\right)$

which might, at first glance, look off-puttingly complicated. To aid clarity, certain groups of terms have been gathered together in the above:

$†$

these are the terms that appeared on the right hand side of the explicit method and are involved with approximating $u_{xx}$ at time $t=n\delta t$

$‡$

these are very similar to the $†$ terms, but all the superscripts are $n+1$ instead of $n$ , that is the terms $‡$ approximate $u_{xx}$ at time $t=\left(n+1\right)\delta t$

(the factor of $\frac{1}{2}$ outside the large bracket shows that we take the average of $†$ and $‡$ )

Figure 7 shows another way of thinking of this numerical method. As in the earlier diagram of this type, arrows point away from positions relating to terms on the right-hand side of the numerical scheme.

Figure 7

The new terms in the Crank-Nicolson method, as compared with the explicit method, give rise to two new unfilled circles on the diagram and the horizontal arrows.

The implementation of this method is similar to that used for the explicit method, but there is a key difference. The Crank-Nicolson scheme is implicit , for consider its use in the first time-step when finding $u_j^1$ ,

$u_j^1=\underset{\mathrm{\✓}}{\underset{\⏟}{u_j^0}}+\frac{1}{2}r\left(\right.\underset{\mathrm{\✓}}{\underset{\⏟}{u_{j-1}^0}}-2\underset{\mathrm{\✓}}{\underset{\⏟}{u_j^0}}+\underset{\mathrm{\✓}}{\underset{\⏟}{u_{j+1}^0}}+\underset{\text{?}}{\underset{\⏟}{u_{j-1}^1-2u_j^1+u_{j+1}^1}}\left.\right)$

The terms labelled $\mathrm{\✓}$ are known from the initial condition. But there are other unknown terms on the right-hand side. We cannot simply “read off" the values at the new time-step as we did using the explicit scheme. Instead we have to store all of the equations given by the stencil at a particular time-step and then solve them as a system of simultaneous equations. The following Example illustrates this point.

Example 17

The temperature $u\left(x,t\right)$ of a metal bar of length $\mathrm{\ℓ}=1.2$ at a distance $x$ from one end and at time $t$ is modelled by the partial differential equation

$u_t=\alpha u_{xx}\left(0<x<\mathrm{\&ell;},t>0\right).$

It is given that the metal has diffusivity $\alpha =1$ , that the two ends of the bar are kept at temperature $u=0$ and that the initial temperature distribution is

$u\left(x,0\right)=f\left(x\right)=x\sqrt{{\left(\mathrm{\&ell;}-x\right)}^3}$

Use the Crank-Nicolson difference scheme with $\delta x=0.4$ and $\delta t=0.1$ to approximate $u\left(x,t\right)$ at $t=\delta t$ and $t=2\delta t$ .

Solution

In this case $r=\alpha \delta t∕{\left(\delta x\right)}^2=0.62500$ so that the numerical scheme can be written

$u_j^{n+1}=u_j^n+\frac{0.62500}{2}\left(u_{j-1}^n-2u_j^n+u_{j+1}^n+u_{j-1}^{n+1}-2u_j^{n+1}+u_{j+1}^{n+1}\right)$

Moving the unknowns to the left of the equation we obtain

$-0.31250u_{j-1}^{n+1}+1.62500u_j^{n+1}-0.31250u_{j+1}^{n+1}=0.37500u_j^n+0.31250\left(u_{j-1}^n+u_{j+1}^n\right)$

The first stage is to use the given data to find $u_j^0$

$\begin{array}{ccccc}\hfill u_0^0& \hfill =\hfill & 0\hfill & \hfill & \text{from the boundary condition}\hfill \\ \hfill u_1^0& \hfill =\hfill & f\left(\delta x\right)=f\left(0.4\right)=0.28622\hfill & \hfill & \text{from the initial condition}\hfill \\ \hfill u_2^0& \hfill =\hfill & f\left(2\delta x\right)=f\left(0.8\right)=0.20239\hfill & \hfill & \text{from the initial condition}\hfill \\ \hfill u_3^0& \hfill =\hfill & 0\hfill & \hfill & \text{from the boundary condition}\hfill \end{array}$

The first time-step will find $u_j^1$ . First we note that the boundary condition implies that $u_0^1=u_3^1=0$ . Two uses of the stencil give

$\begin{array}{ccc}\hfill -0.31250u_0^1+1.62500u_1^1-0.31250u_2^1=0.37500u_1^0+0.31250\left(u_0^0+u_2^0\right)& \hfill =\hfill & 0.17058\hfill \\ \hfill -0.31250u_1^1+1.62500u_2^1-0.31250u_3^1=0.37500u_2^0+0.31250\left(u_1^0+u_3^0\right)& \hfill =\hfill & 0.16534\hfill \\ \hfill \end{array}$

The implicit nature of this method means that we have to do some extra work to complete the time-step. We must now solve the simultaneous equations

$\left(\begin{array}{cc}\hfill 1.62500\hfill & \hfill -0.31250\hfill \\ \hfill -0.31250\hfill & \hfill 1.62500\hfill \end{array}\right)\left(\begin{array}{c}\hfill u_1^1\hfill \\ \hfill u_2^1\hfill \end{array}\right)=\left(\begin{array}{c}\hfill 0.17058\hfill \\ \hfill 0.16534\hfill \end{array}\right)$

In this case there are only two unknowns and it is a simple matter to solve the pair of equations to give $u_1^1=0.12932$ and $u_2^1=0.12662$ .

The second time-step will find $u_j^2$ . First we note that the boundary condition implies that $u_0^2=u_3^2=0$ . Two uses of the stencil give

$\begin{array}{ccc}\hfill -0.31250u_0^2+1.62500u_1^2-0.31250u_2^2=0.37500u_1^1+0.31250\left(u_0^1+u_2^1\right)& \hfill =\hfill & 0.08806\hfill \\ \hfill -0.31250u_1^2+1.62500u_2^2-0.31250u_3^2=0.37500u_2^1+0.31250\left(u_1^1+u_3^1\right)& \hfill =\hfill & 0.08789\hfill \\ \hfill \end{array}$

The implicit nature of this method means that we have to do some extra work to complete the time-step. We must now solve the simultaneous equations

$\left(\begin{array}{cc}\hfill 1.62500\hfill & \hfill -0.31250\hfill \\ \hfill -0.31250\hfill & \hfill 1.62500\hfill \end{array}\right)\left(\begin{array}{c}\hfill u_1^2\hfill \\ \hfill u_2^2\hfill \end{array}\right)=\left(\begin{array}{c}\hfill 0.08806\hfill \\ \hfill 0.08789\hfill \end{array}\right)$

In this case there are only two unknowns and it is a simple matter to solve the pair of equations to give $u_1^2=0.06707$ and $u_2^2=0.06699$ .

Figure 8 depicts the numerical solutions found in Example 17 above. (Again, the dotted lines are intended to aid clarity, they are not part of the numerical solution.)

Figure 8

Task!

The temperature $u\left(x,t\right)$ of a metal bar of length $\mathrm{\&ell;}=0.9$ at a distance $x$ from one end and at time $t$ is modelled by the partial differential equation

$u_t=\alpha u_{xx}\left(0<x<\mathrm{\&ell;},t>0\right).$

It is given that the metal has diffusivity $\alpha =0.25$ , that the two ends of the bar are kept at temperature $u=0$ and that the initial temperature distribution is

$u\left(x,0\right)=f\left(x\right)=sin\left(\pi x∕\mathrm{\&ell;}\right)$

Use the Crank-Nicolson difference scheme with $\delta x=0.3$ and $\delta t=0.2$ to approximate $u\left(x,t\right)$ at $t=\delta t$ and $t=2\delta t$ .

Answer

The second time-step will find $u_j^2$ . First we note that the boundary condition implies that $u_0^2=u_3^2=0$ . Two uses of the stencil give

$\begin{array}{ccc}\hfill -0.27778u_0^2+1.55556u_1^2-0.27778u_2^2=0.44444u_1^1+0.27778\left(u_0^1+u_2^1\right)& \hfill =\hfill & 0.35352\hfill \\ \hfill -0.27778u_1^2+1.55556u_2^2-0.27778u_3^2=0.44444u_2^1+0.27778\left(u_1^1+u_3^1\right)& \hfill =\hfill & 0.35352\hfill \\ \hfill \end{array}$

The implicit nature of this method means that we have to do some extra work to complete the time-step. We must now solve the simultaneous equations

$\left(\begin{array}{cc}\hfill 1.55556\hfill & \hfill -0.27778\hfill \\ \hfill -0.27778\hfill & \hfill 1.55556\hfill \end{array}\right)\left(\begin{array}{c}\hfill u_1^2\hfill \\ \hfill u_2^2\hfill \end{array}\right)=\left(\begin{array}{c}\hfill 0.35352\hfill \\ \hfill 0.35352\hfill \end{array}\right)$

In this case there are only two unknowns and it is a simple matter to solve the pair of equations to give $u_1^2=0.27667$ and $u_2^2=0.27667$ .

6.1 In general

Having now seen some instances with a relatively large $\delta x$ , we now look at the general case where the space step may be much smaller. In this case there will be a larger system of equations to solve at each time-step than was the case above.

In general, the procedure of moving the unknowns to the left hand side of the equation leads to

$-\frac{r}{2}{u}_{j-1}^{n+1}+\left(1+r\right)u_j^{n+1}-\frac{r}{2}{u}_{j+1}^{n+1}=\frac{r}{2}{u}_{j-1}^n+\left(1-r\right)u_j^n+\frac{r}{2}{u}_{j+1}^n$

which we apply all the way along the $x$ -axis. That is, we put $j=1,2,3,\dots ,J-1$ in the above expression and hence derive a system of equations for all the $u$ with superscript $n+1$ .

$\left(\begin{array}{cccccc}\hfill 1+r\hfill & \hfill -\frac{r}{2}\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill -\frac{r}{2}\hfill & \hfill 1+r\hfill & \hfill -\frac{r}{2}\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill 0\hfill & \hfill -\frac{r}{2}\hfill & \hfill 1+r\hfill & \hfill -\frac{r}{2}\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill \\ \hfill \&vellip;\hfill & \hfill \hfill & \hfill \&dtdot;\hfill & \hfill \&dtdot;\hfill & \hfill \&dtdot;\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill -\frac{r}{2}\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill -\frac{r}{2}\hfill & \hfill 1+r\hfill \end{array}\right)\left(\begin{array}{c}\hfill u_1^{n+1}\hfill \\ \hfill u_2^{n+1}\hfill \\ \hfill u_3^{n+1}\hfill \\ \hfill \hfill \\ \hfill \&vellip;\hfill \\ \hfill \hfill \\ \hfill u_{J-1}^{n+1}\hfill \end{array}\right)=\left(\begin{array}{c}\frac{r}{2}\underset{̲}{u_0^n}+\left(1-r\right)u_1^n+\frac{r}{2}{u}_2^n+\frac{r}{2}\underset{̲}{\underset{̲}{u_0^{n+1}}}\hfill \\ \frac{r}{2}{u}_1^n+\left(1-r\right)u_2^n+\frac{r}{2}{u}_3^n\hfill \\ \frac{r}{2}{u}_2^n+\left(1-r\right)u_3^n+\frac{r}{2}{u}_4^n\hfill \\ \hfill \\ \&vellip;\hfill \\ \hfill \\ \frac{r}{2}{u}_{J-2}^n+\left(1-r\right)u_{J-1}^n+\frac{r}{2}\underset{̲}{u_J^n}+\frac{r}{2}\underset{̲}{\underset{̲}{u_J^{n+1}}}\hfill \end{array}\right)$

The underlined terms on the right-hand side will be known from the boundary conditions. The doubly underlined quantities are “new" at the current time-step and involve the only appearances of $n+1$ on the right-hand side. All the other $u$ approximations at time level $n+1$ are unknown at this stage and appear on the left.

The matrix on the left-hand side of the system has the following properties

• It is independent of $n$ . In other words, the same matrix appears at each time-step. (We saw this in the example and exercise above in which the same $2\times 2$ matrix appeared at each of the two time-steps carried out).
• It is tridiagonal . That is, the only non-zero entries are either on, or adjacent to, the diagonal. Furthermore, there are only two different values ( $\frac{r}{2}$ and $1+r$ ) which appear. This is good news as far as storage is concerned. Gaussian elimination (seen in HELM booklet  30, for example) works extremely well on tridiagonal matrices.

It is also true that the matrix is strictly diagonally dominant . (That is, the diagonal element on each row is greater in size than the sum of the absolute values of the off-diagonal elements on that row.) This means that methods such as Jacobi and Gauss Seidel (see HELM booklet  30 for details) would work very well.

6.2 Stability of the Crank-Nicolson scheme

This is the big pay-off when using the Crank-Nicolson method.

This is excellent news. It means that there is no hideously restrictive constraint on the size of $\delta t$ .