7.5.3 Crank-Nicolson scheme. \( \theta \)-scheme
One of the bad characteristics of the DuFort-Frankel scheme is that one needs a special procedure at the starting time, since the scheme is a 3-level scheme. Therefore, we try now to find a second order approximation for \( \frac{\partial u}{\partial t} \) where only two time levels are required..
Using central differences centered at the half time interval we obtain: $$ \begin{equation} \tag{7.60} \frac{\partial u}{\partial t}\bigg|^{n+\frac{1}{2}}_j=\frac{u_j^{n+1}-u^n_j}{\Delta t}+O(\Delta t)^2 \end{equation} $$
The problem now is to approximate \( \dfrac{\partial u}{\partial t}\Big|^{n+\frac{1}{2}}_j \) without knowing without including any term evaluate at \( n+\frac{1}{2} \) explicitly in the scheme. This was achieved by the Crank-Nicolson approximation (1947): $$ \begin{equation} \tag{7.61} \frac{\partial^2 u}{\partial x^2}\bigg|^{n+\frac{1}{2}}_j=\frac{1}{2}\left[ \frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{(\Delta x)^2} +\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{(\Delta x)^2} \right]+O(\Delta x)^2 \end{equation} $$
The difference equation is now: $$ \begin{equation} \tag{7.62} u_{j-1}^{n+1}-2(1+\frac{1}{D})u_j^{n+1}+u_{j+1}^{n+1} = -u^n_{j-1}+2(1-\frac{1}{D})u^n_j-u^n_{j+1} \end{equation} $$
with \( D= \dfrac{\Delta t}{(\delta x)^2} \).
The stencil of the resulting numerical scheme is shown in Figure 94.
Figure 94: 2-level stencil of the Crank-Nicolson scheme. Note that the scheme does not require any evaluation of the solution at time level \( n+\frac{1}{2} \).
From (7.62) and Figure 94 we see that the scheme is implicit. This in turn implies that the algebraic equations resulting from the discretization of the original differential problem are coupled and, differently from an explicit scheme, we have to solve them as a system. In this case the system is tridiagonal, so we can use the Thomas' algorithm, for which the code tdma was introduce in the section 4.3 Tridiagonal systems of algebraic equations.
Before proceeding with a stability analysis for (7.62), we will re-write it a more general form, by introducing the parameter \( \theta \): $$ \begin{equation} \tag{7.63} u^{n+1}_j=u^n_j+D\left [\theta(u^{n+1}_{j+1}-2u_j^{n+1}+u^{n+1}_{j-1})+(1-\theta)(u^n_{j+1}-2u^n_j+u^n_{j-1}) \right ] \end{equation} $$ with $ 0\leq\theta\leq 1$. This scheme is a generalization of Crank-Nicolson original scheme and is called the \( \theta \)-scheme.
For \( \theta = 0 \) one recovers the explicit FTCS-scheme. On the other hand, for$\theta = \frac{1}{2}$ one has the original Crank-Nicolson scheme. Finally, for \( \theta = 1 \) one gets an implicit scheme some times called the Laasonen-scheme (1949). In fluid mechanics this last scheme is often called BTCS-scheme. (Backward in Time and Central in Space).
Inserting (7.42) in (7.63) and dividing by \( G^n\cdot e^{i\cdot \beta x_j} \) gives: $$ \begin{equation*} \begin{array}{llc} G&=1+D[G\theta(e^{i\delta}+e^{-i\cdot \delta}-2)+(1-\theta)(e^{i\delta}+e^{-i\delta}-2)]&\\ &=1+D(e^{i\delta}+e^{-i\delta}-2)\cdot(G\theta+1-\theta)&\\ \end{array} \end{equation*} $$ with \( \delta = \beta \cdot h \). Moreover, using (7.50) results in: $$ \begin{equation} \tag{7.64} G=\frac{1-4D(1-\theta)\sin^2(\frac{\delta}{2})}{1+4D\theta\sin^2(\frac{\delta}{2})} \end{equation} $$
The stability condition is \( |G|\leq1 \) or, since \( G \) is real: \( -1\leq G\leq1 \). Since we have that \( 0\leq\theta\leq1 \) , the condition of the right side is satisfied for \( D\geq0 \). For the left side of the inequality we have that: $$ \begin{equation} 2D\sin^2\left(\frac{\delta}{2}\right)(1-2\theta)\leq 1 \tag{7.65} \end{equation} $$ or $$ \begin{equation} D(1-2\theta)\leq \frac{1}{2}\text{ da } \sin^2\left(\frac{\delta}{2}\right)\leq 1 \tag{7.66} \end{equation} $$
For \( \frac{1}{2}\leq \theta \leq 1 \), stability is given for all \( D\geq0 \), i.e.: the scheme is unconditionally stable.
For \( 0\leq \theta \leq \frac{1}{2} \), the stability condition requires that: $$ \begin{equation} \tag{7.67} D(1-2\theta)\leq \frac{1}{2},\ 0\leq \theta \leq 1 \end{equation} $$
Use of derivation
Now we write: $$ \begin{equation*} \begin{array}{l} G=1+D(e^{i\delta}+e^{-i\delta}-2)\cdot(G\theta +1 -\theta)=1+2D\big(\cos(\delta)-1\big)\cdot(G\theta+1-\theta)\\ \dfrac{dG}{d\delta}=2D\left[(G\theta+1-\theta)\dfrac{d}{d\delta}\big(\cos(\delta)-1\big)+\big(\cos(\delta)-1\big)\theta\cdot \dfrac{dG}{d\delta}\right] \end{array} \end{equation*} $$
With \( \dfrac{dG}{d\delta}=0 \) we get the max-min for \( \delta =0,\ \delta=\pm\pi \) (\( \delta=0 \) gives \( G=1 \) as expected) \( \delta= \pm \delta \) gives \( G=\dfrac{1-4D(1-\theta)}{1+4D\theta} \) which is identical to (7.64) for \( \frac{\delta}{2}=90^\circ \). The condition \( |G|\leq1 \) becomes \( -1\leq \dfrac{1-4D(1-\theta)}{1+4D\theta}\leq1 \) with the same results as before.
Accuracy
Let us now look at the accuracy of the scheme presented in (7.63).
To do so, we write a simple Maple-program:
> eq1:= u(x+h,t+k) - 2*u(x,t+k) + u(x-h,t+k):
> eq2:= u(x+h,t) - 2*u(x,t) + u(x-h,t):
> eq:= (u(x,t+k) - u(x,t))/k -(theta*eq1 + (1-theta)*eq2)/h^2:
> Tnj:= mtaylor(eq,[h,k]):
> Tnj:= simplify(Tnj):
> Tnj:= convert(Tnj,diff);
Note that we have used \( h=\Delta x \) and \( k=\Delta t \). \( Tnj \) is the truncation error, which is discussed in detail in the section 7.6 Truncation error, consistency and convergence. Suppose now that \( u(x,t) \) is the analytical solution for the partial differential equation \( \dfrac{\partial u}{\partial t}=\dfrac{\partial^2 u}{\partial x^2} \) so that we can say that \( \dfrac{\partial}{\partial t}()=\dfrac{\partial^2}{\partial x^2} \), etc. If we use these relationships in the Maple-program for \( Tnj \), we will obtain the following result: $$ \begin{equation} \tag{7.68} T^n_j=\left[\bigg(\frac{1}{2}-\theta\bigg)\cdot \Delta t-\frac{1}{12}(\Delta x)^2\right]\cdot \frac{\partial^4u}{\partial x^4}+\frac{1}{6}(\Delta t)^2(1-3\theta)\cdot \frac{\partial^6u}{\partial x^6}+\dots \end{equation} $$
We see that the truncation error is \( T^n_j=O(\Delta t)+O(\Delta x)^2 \) for \( \theta=0 \) and 1, i.e. for Euler-scheme and Laasonen-scheme, while it becomes \( T^n_j=O(\Delta t)^2+O(\Delta x)^2 \) for \( \theta=\frac{1}{2} \), which corresponds to the Crank-Nicolson scheme.
If we add the following line
Tnj:= simplify(subs(theta =(1-h^2/(6*k))/2,Tnj));
to the Maple-program, we find that \( T^n_j=O(\Delta t)^2+O(\Delta x)^4 \) if we choose \( D=\dfrac{1}{6(1-2\theta)} \) when \( D\leq \dfrac{1}{2(1-2\theta)} \).
More about stability
We have found that the Crank-Nicolson scheme is unconditionally stable for \( \frac{1}{2}\leq \theta \leq 1 \). In practice, the scheme can cause oscillations around discontinuities for \( \theta=\frac{1}{2} \). Values of \( \theta \) above \( 0.5 \) will dampen such oscillations with strongest attenuation for \( \theta=1 \). Let's now look at an example where we use the heat equation in dimensional form.
Figure 95: Thin aluminum rod for one-dimensional heat equation example.
Figure 95 shows a thin aluminum rod of length equal to 300 mm.
The one-dimensional heat equation is given by: $$ \begin{equation*} \frac{\partial T}{\partial t}=\alpha \frac{\partial^2T}{\partial x^2},\ T=T(x,t) \end{equation*} $$
Boundary conditions are: $$ \begin{equation*} T(0,t)=T(300, t)=20^\circ C, \end{equation*} $$
while initial conditions are: $$ \begin{equation*} \begin{array}{l} T(x,0)=270^\circ C\text{ for } x\in (100,200)\\ T(x,0)=20^\circ C\text{ for } x\in (0,100) \text{ and } x\in(200,300) \end{array} \end{equation*} $$
Thermal diffusivity is: $$ \begin{equation*} \alpha =100 mm^2/s \end{equation*} $$
The numerical Fourier number is: $$ \begin{equation*} D=\alpha \frac{\Delta t}{(\Delta x)^2} \end{equation*} $$
Furthermore, we set \( \Delta t=0.25 \) s so that \( \Delta x=5/\sqrt{D} \). By selecting \( D=4 \), we get that \( \Delta x = 2.5 \) mm. Figures 96 and 97 show computational results with \( \theta= \frac{1}{2} \) and with \( \theta=1 \); the Crank-Nicolson and Laasonen schemes, respectively.
Figure 96: Solution for the one-dimensional heat equation problem using Crank-Nicolson scheme.
Figure 97: Solution for the one-dimensional heat equation problem using Laasonen scheme.
We see that the solution obtained using the C-N-scheme contains strong oscillations at discontinuities present in the initial conditions at \( x=100 \) of \( x=200 \). These oscillation will be dampen by the numerical scheme, in this case they tend to disappear for \( t=4 \) s. The solution delivered by the Laasonen-scheme presents no oscillations, even around \( t=0 \). To explain the behavior of these solutions we have to go back to equation (7.64): $$ \begin{equation*} G=\frac{1-4D(1-\theta)\sin^2(\frac{\delta}{2})}{1+4D\sin^2(\frac{\delta}{2})} \end{equation*} $$
For \( \theta=\frac{1}{2} \) we get: $$ \begin{equation*} G=\frac{1-2D\sin^2(\frac{\delta}{2})}{1+2D\sin^2(\frac{\delta}{2})} \end{equation*} $$
For \( \delta \) near \( \pi \), \( G \) will have a value close to \( -1 \) for large values of \( D \). This can be clearly seen if one sets \( \delta = \pi=\beta \cdot h \), which results in \( G=\dfrac{1-2D}{1+2D} \).
From (7.44) we have: $$ \begin{equation*} E^n_j=G^n\cdot E^0_j,\ E^0_j=e^{i\cdot \beta x_j} \end{equation*} $$
Therefore, we see that we will get oscillations for \( \delta \) close to \( \pi \). These high wave numbers will be dampened as the simulation proceeds because the initial condition discontinuities are smoothed by dissipation. If we set \( D=4 \) as used in this example, after 16 time increments with $\Delta t=0.25$s we get that \( G^{16}=\left(-\dfrac{7}{9}\right)^{16}\approx 0.018 \), whereas with \( D=1 \) we would obtain \( G^{16}=2.3\cdot 10^{-8} \). This means that high wave-numbers oscillations dampening slows down for large values of \( D \). This is very similar to what happens in the case of the Gibbs-phenomenon for Fourier series of discontinuous functions (See Appendix A.3, in the Numeriske Beregninger).
On the other hand, for the Laasonen-scheme, i.e. the \( \theta \)-scheme, we have that: $$ \begin{equation*} G=\dfrac{1}{1+4D\sin^2\left(\dfrac{\delta}{2}\right)} \end{equation*} $$
We do not get oscillations for any wave number.
If we use terms from the stability analysis of ordinary differential equations we can say that the \( \theta \)-scheme is absolute stable for (A-stable) for \( \theta=\frac{1}{2} \) and \( L \)-stable for \( \theta=1 \) . (See section 1.6.1 in the Numeriske Beregninger)
If we use the \( \theta \)-scheme for the heat equation in a problem setting as the one described above with prescribed temperature at both ends (Dirichlet boundary conditions), we know that the maximum temperature at any given time must be between the largest value of the initial condition and the smallest given on the edge of the domain, or \( T_{min}\leq T^n_j\leq T_{max} \). For the example presented above one has \( T_{min}=20^\circ C \) and \( T_{max}=270^\circ C \).
Writing (7.63) for \( u_j^{n+1} \) yields: $$ \begin{equation} \tag{7.69} u^{n+1}_j=\frac{1}{1+2\theta D}\left[\theta D(u^{n+1}_{j-1}+u^{n+1}_{j+1})+(1-\theta) D(u^{n}_{j-1}+u^{n}_{j+1}) +\big(1-(1-\theta) 2D \big)\right] \end{equation} $$
By summing the coefficients on the right-hand side we find that the sum is equal to 1. Next, we check for which conditions the coefficients are positive, which is equivalent to require that \( 0 < \theta < 1 \) and \( 1-(1-\theta)\cdot 2D > 0 \). The last inequality is met for \( D\cdot(1-\theta) < \frac{1}{2} \). By setting \( \theta=0 \) and \( \theta=1 \), we find that the sum is equal to 1 for these values, so that we get the following condition for fulfilling the PC-criterion: $$ \begin{equation} \tag{7.70} D\cdot (1-\theta)\leq \frac{1}{2} \end{equation} $$
By the von Neumann-analysis we found that the following relation should hold (see (7.67)): $$ \begin{equation} \tag{7.71} D\cdot (1-2\theta)\leq \frac{1}{2} \end{equation} $$
We see that (7.70) is significantly stricter than (7.71). While the \( \theta \)-scheme is unconditionally stable for \( \theta=\frac{1}{2} \) according to the von Neumann stability analysis, we must have \( D\leq 1 \) according to the PC-criterion. The PC-criterion provides a sufficient condition which ensures that the physical max-min criterion is also fulfilled by the solution delivered by the numerical scheme. The example shown above confirms the prediction provided by the PC-criterion. Now, is there a criterion that is both necessary and sufficient for this simple model equation? Kraaijevanger found the following necessary and sufficient criterion in 1992: $$ \begin{equation} \tag{7.72} D\cdot (1-\theta)\leq \frac{2-\theta}{4(1-\theta)} \end{equation} $$
We see that for \( \theta = \frac{1}{2} \) this criterion implies the condition \( D\leq \frac{3}{2} \). For \( \theta = \frac{3}{4} \) it gives (7.72) \( D\leq 5 \), while the PC-criterion yields \( D\leq 2 \).
The purpose of allowing for large \( D \) values is related to the fact that one might be interested in simulations where a large time step is set to let reach a steady state. But on must not forget about accuracy. Moreover, we should remember that this is a simple one-dimensional partial differential equation with constant coefficients, and as such it can be solver very quickly with a modern PC with for any reasonable \( \Delta t \) and \( \Delta x \), as long as we stay within the stability area. However, three-dimensional non-linear problems will in general result in very computational demanding simulations.