Crank-Nicolson Scheme for the 1D Heat Equation ME 448/548
1. Combine finite difference approximations for ∂u/∂t at x = xi
uki − uk−1
∂u i
= + O(∆t) (1)
∂t tk ,xi ∆t
and the average of the ∂ 2 u/∂x2 operators at time tk and at tk+1
" # " #
k+1 k+1 k+1
∂ 2 u α ui−1 − 2ui + ui+1 α uki−1 − 2uki + uki+1
= + + O(∆x2 ). (2)
∂x2 tk ,xi 2 ∆x2 2 ∆x2
to get a system of equations having the same structure as the BTCS method
α 1 α α
− uk+1 + + uk+1 − uk+1 =
2∆x2 i−1 ∆t ∆x2 i
2∆x2 i+1
α k 1 α α
u + − uki + uk (3)
2∆x2 i−1 ∆t ∆x2 2∆x2 i+1
Equation (3) is the computational formula for the Crank-Nicolson scheme. It is an implicit
scheme because all uk+1 values are coupled and must be updated simultaneously.
2. Computational Molecule
t .. .. .. .. .. .. ..
. . . . . . . Crank-Nicolson scheme
requires simultaneous
k+1
calculation of u at all
k nodes on the k+1 mesh line
k 1
Solution is known
for these nodes
t=0, k=1
i=1 i 1 i i+1 nx
x=0 x=L
3. Stability:
The Crank-Nicolson method is unconditionally stable for the heat equation.
The benefit of stability comes at a cost of increased complexity of solving a linear system of
equations at each time step. The Crank-Nicolson scheme is not significantly more costly to
implement than the BTCS Scheme
4. The Crank-Nicolson scheme has a truncation error that is O(∆t2 ) + O(∆x2 )
2
5. For the one-dimensional heat equation, the linear system of equations for the Crank-Nicolson
scheme can be organized into a tridiagonal matrix that looks just like the tridiagonal matrix
for the BTCS scheme.
k+1
u1
a1 b1 0 0 0 0 k+1 d1
c2 a2 b2 u2
0 0 0 d2
0 c3 a3 k+1
b 3 0 0
u d3
3 = . (4)
0 0 ... .. ..
. ..
. . 0 ..
0 0 0 cnx −1 anx −1 bnx −1 uk+1 dnx −1
nx −1
0 0 0 0 cnx anx dnx
uk+1
nx
The coefficients of the interior nodes (i = 2, 3, . . . , N − 1) are
ai = 1/∆t + α/∆x2 = 1/∆t − (bi + ci ),
bi = ci = −α/(2∆x2 ),
di = −ci uki−1 + (1/∆t + bi + ci ) uki − bi uki+1 .
As with the BTCS scheme, this system of equations is efficiently solved with a form of LU
factorization. The LU factors need to be computed only once before the first time step.
6. Matlab implementation: code from demoCN
% --- Coefficients of the tridiagonal system
b = (-alfa/2/dx^2)*ones(nx,1); % Super diagonal: coefficients of u(i+1)
c = b; % Subdiagonal: coefficients of u(i-1)
a = (1/dt)*ones(nx,1) - (b+c); % Main Diagonal: coefficients of u(i)
at = (1/dt + b + c); % Coefficient of u_i^k on RHS
a(1) = 1; b(1) = 0; % Fix coefficients of boundary nodes
a(end) = 1; c(end) = 0;
[e,f] = tridiagLU(a,b,c); % Save LU factorization
% --- Assign IC and save BC values in ub. IC creates u vector
x = linspace(0,L,nx)’; u = sin(pi*x/L); ub = [0 0];
% --- Loop over time steps
for k=2:nt
% --- Update RHS for all equations, including those on boundary
d = - [0; c(2:end-1).*u(1:end-2); 0] ...
+ [ub(1); at(2:end-1).*u(2:end-1); ub(2)] ...
- [0; b(2:end-1).*u(3:end); 0];
u = tridiagLUsolve(e,f,b,d); % Solve the system
end
A more general implementation is in heatCN.
7. A comparison of FTCS, BTCS and Crank-Nicolson shows Truncation Errors
that all three have the same spatial truncation error. FTCS Scheme Spatial Temporal
and BTCS have the same temporal truncation error. Crank- FTCS ∆x2 ∆t
Nicolson has superior temporal truncation error. BTCS ∆x2 ∆t
C-N ∆x2 ∆t2
3
8. The compHeatSchemes function shows that our Matlab implementation of all three schemes
demonstrate the correct behavior of truncation error.
>> compHeatSchemes
Reduce both dx and dt within the FTCS stability limit
------------- Errors ------------
nx nt FTCS BTCS CN
4 5 2.903e-02 5.346e-02 1.304e-02
8 21 6.028e-03 1.186e-02 2.929e-03
16 92 1.356e-03 2.716e-03 6.804e-04
32 386 3.262e-04 6.522e-04 1.630e-04
64 1589 7.972e-05 1.594e-04 3.984e-05
128 6453 1.970e-05 3.939e-05 9.847e-06
256 26012 4.895e-06 9.790e-06 2.448e-06
512 104452 1.220e-06 2.440e-06 6.101e-07
Reduce dt while holding dx = 9.775171e-04 (L=1.0, nx=1024) constant
------------- Errors ------------
nx nt FTCS BTCS CN
1024 8 NaN 2.601e-02 1.291e-03
1024 16 NaN 1.246e-02 2.798e-04
1024 32 NaN 6.102e-03 6.534e-05
1024 64 NaN 3.020e-03 1.570e-05
1024 128 NaN 1.502e-03 3.749e-06
1024 256 NaN 7.492e-04 8.154e-07
1024 512 NaN 3.742e-04 8.868e-08
1024 1024 NaN 1.871e-04 9.218e-08
−1 −1
10 10
FTCS
BTCS ∆ x = 9.8e−04 (constant)
−2
−2 CN 10
10
E ∝∆ x2
−3
10
−3
10
−4
10
E(∆ x,∆ t)
E(∆ x,∆ t)
−4
10
−5
10
−5
10
−6
10
FTCS
−6 BTCS
10 −7
CN
10
E∝∆ t
2
E∝∆ t
−7 −8
10 −3 −2 −1 0
10 −3 −2 −1 0
10 10 10 10 10 10 10 10
∆x ∆t
In the plot of truncation error versus ∆t (right hand plot), there is an irregularity at ∆t ∼
3.9 × 10−3 . At that level of ∆t, and for the chosen ∆x, which is held constant, the truncation
error due to ∆x is no longer negligible. Further reductions in ∆t alone will not reduce the total
truncation error.
