CHAPTER 9: Partial Differential Equations 205

9.6 Solving the Heat Equation using the Crank-Nicholson

Method
The one-dimensional heat equation was derived on page 165. Let’s generalize it to allow for the
direct application of heat in the form of, say, an electric heater or a flame:

 T  x, t   2 T  x, t 
D  Papplied  x, t  . (9.97)
t  x2

The new term Papplied  x, t  is the power applied (i.e. the rate at which heat energy is applied) at
point x at time t. It is worthwhile to compare this generalized heat equation with the generalized
wave equation, (9.93). The major difference is that the heat equation has a first time derivative
whereas the wave equation has a second time derivative (if we ignore resistance).

We now wish to approximate the derivatives with finite differences. There are several possible
schemes for doing this and we must see which one works best. The first scheme is called
Forward Time Centered Space (FTCS) because it uses a forward difference to approximate the
time derivative and a centered second difference to approximate the second space derivative.
Thus:
T j n 1  T j n T j n1  2 T j n  T j n1
 D  Pj n
t  x
2
, (9.98)
  
forward from time n centered at j at time n

The definitions of n, j,  t and  x are the same as on page 197 and T  x, t   T j n and
Papplied  x, t   Pj n . Solving for T j n 1 gives

T j n 1  T j n   T j n1  2 T j n  T j n1    t  Pj n , (9.99)

t
where   D is a dimensionless constant that measures the speed of the simulation.
( x) 2
Like Eq. (9.81) for waves, Eq. (9.99) for heat flow is an explicit equation – meaning that unknown
grid point values at timestep n+1 can be calculated directly from known values at timestep n. Unfor-
tunately the scheme is unstable unless  t is made so small that the simulation is impossibly slow.

The second scheme is Backward Time Centered Space (BTCS). It uses a backward difference
to approximate the time derivative:

T j n 1  T j n T j n11  2 T j n 1  T j n11
 D  Pj n 1
t  x
2
, (9.100)
   
backward from time n 1 centered at j at time n 1

Notice that this is a relation among 3 future values and one present value. Let us rearrange this so
that all the future values are isolated on the LHS. This gives:

 T j 1n 1  1  2   T j n 1   T j 1n 1  T j n   t  Pj n 1
     , (9.101)
3 unknowns known value
 206 CHAPTER 9: Partial Differential Equations

(a) FTCS method (b) BTCS (fully implicit) method (c) Crank-Nicolson method

n +1 n +1 n +1 n +1 n +1 n +1 n +1
uj u j 1 uj u j +1 u j 1 uj u j +1
n +1 n +1 n +1
n n n n n n n
time t u j 1 uj u j +1 uj u j 1 uj u j +1
n n n
index n

n 1 n 1 n 1

j1 j j+1 j1 j j+1 j1 j j+1

position index j position index j position index j

Figure 9.34 Stencils for three schemes for finite-differencing the heat equation.

We see that this is an implicit equation – to solve it means to solve a set of simultaneous linear
equations at each timestep. Fortunately this is not a big problem since the system is tridiagonal.
Also, as we will see later, this scheme is unconditionally stable.

The third scheme is called the Crank-Nicolson method. It takes the average of (9.98) and
(9.100):

T j n 1  T j n 1  T j n1  2 T j n  T j n1 T j n11  2 T j n 1  T j n11 

 D  P n
 D  P n1
 (9.102)
t 2   x
2 j
 x
2 j

 

Although this is again an implicit equation, it has the advantage of being symmetrical in time and
is thus more accurate than either FTCS or BTCS. It is helpful to compare the stencils for the three
schemes. They are shown in Fig. 9.34. Solving (9.102) for the new temperatures in terms of old
temperatures gives

   
 T j 1n 1  1    T j n 1  T j 1n 1  T j 1n  1    T j n  T j 1n   t  Pj n  Pj n 1  . (9.103)
2 2 2 2 2

or
 
 T j 1n 1  1    T j n 1  T j 1n 1  R j n
2 2 , (9.104)
   
3 unknowns known value

where
 
Rjn  T j 1n  1    T j n  T j 1n   t  Pj n  Pj n 1  . (9.105)
2 2 2
n
The quantity Rj that appears on the RHS is a known quantity at the beginning of each timestep.
To become familiar with Eq. (9.104) suppose that there are 5 grid points numbered 0 to 4.
Suppose also that the boundary conditions are that the temperature equals a at the left boundary
and b at the right boundary (Dirichlet conditions). Then we must solve this matrix equation:

 1 0 0 0 0  T n 1   a 
 

1 

0 0   T n 1   R n 
0

 2 2  1   1 
  
0  1  0  T2 n 1    R2 n  . (9.106)
 2 2 
    n 1   n 
0 0  1   T3   R3 
 2 2 
 n 1   
 0 0 0 0 1  T4   b 

The first and last equations in it are just the boundary conditions. The middle three equations are
(9.104). Only the three shaded bands are non-zero. This is why the system is called tridiagonal.
 CHAPTER 9: Partial Differential Equations 207

As usual, if we have Neumann boundary conditions (i.e. insulated boundaries) then we need
special versions of (9.104). For the left boundary we get it by letting j = 0 in (9.104) and using the
fact that T1n 1  T1n 1 . This gives

1    T0 n 1   T1n 1  R 0n , where R 0n  1    T0 n   T1n   t  P0 n  P0 n 1  (9.107)

2

For the right boundary we let j = J in (9.104) and use the fact that TJ 1n 1  TJ 1n 1 . This gives

  T J n11  1    TJ n 1  RJ n , where R J n  1    T J n   T J 1n   t  PJ n  PJ n 1  (9.108)

2

Thus for Neumann boundary conditions we must solve this matrix equation:

 1    0 0 0  T0 n 1   R 0n 
   1   0 0   n 1   n 
 2 2   T1   R1 
 0   
 1  0 T2 n 1    R2 n  . (9.109)
 2 2   n 1   n 
 0   
0  1  T3   R3 
 2 2  T n 1   R n 
 0 0 0  1     4   4 

The first and last equations are the Neumann the boundary conditions, (9.107) and (9.108). The
middle three equations are (9.104). Again the matrix is tridiagonal.

On the next page is Visual Basic code that solves the heat flow problem in Example 9.1 and Fig.
9.1. To load it into Excel follow the instructions given on page 194. To use it run the macro
called CrankNicolson. To create a graph of the temperature profile select columns And B in the
spreadsheet, click on the Insert tab, and choose a chart of the Scatter type. Problem Set 9.6
explores this code and makes modifications to it to solve other heat flow problems.
208 CHAPTER 9: Partial Differential Equations

Public Sub CrankNicolson ( )

Const TLeft = 5, TRight = 15, Alpha = 1 '[1]
Const Iterations = 1, Right = 10
Dim J As Integer, Timestep As Integer, UserResponse 'Index j runs from 0 to Right. [2]
Dim T(0 To Right) 'Various arrays: Temperature distribution T.
Dim A(0 To Right), B(0 To Right), C(0 To Right) 'Tridiagonal matrix's bands, A, B, C.
Dim R(0 To Right), U(0 To Right) 'R and U in equation: [Matrix]*U=R.
For J = 0 To Right 'Set up the tridiagonal matrix. This has to be done just once.
A(J) = -Alpha / 2 'below-diagonal band.
B(J) = 1 + Alpha 'diagonal band.
C(J) = -Alpha / 2 'above-diagonal band.
Next J
B(0) = 1: C(0) = 0: A(Right) = 0: B(Right) = 1 'For Dirichlet conditions. [3]
Timestep = 0
For J = 0 To Right: T(J) = 0: Next J 'Set up initial temperature profile. [4]
Cells(2, 1) = “Position”: Cells(2, 2) = “Temperature” 'Set up spreadsheet columns.
For J = 0 To Right: Cells(J + 3, 1) = J: Cells(J + 3, 2) = T(J): Next J
MsgBox ("Timestep = 0, the initial temperature profile")
Do 'This is the loop. Each loop does 1 timestep.
For J = 1 To Right - 1 'Fill in the matrix R (the RHS matrix).
R(J) = Alpha / 2 * (T(J - 1) + T(J + 1)) + (1 - Alpha) * T(J)
Next J
R(0) = TLeft: R(Right) = TRight 'For Dirichlet conditions. [5]
Call TriDiagonalSolver(Right, A, B, C, R, U)
Timestep = Timestep + 1
For J = 0 To Right: T(J) = U(J): Next J 'Copy output of solver to temperature profile.
For J = 0 To Right: Cells(J + 3, 2) = T(J): Next J 'Update the spreadsheet.
If Timestep Mod Iterations = 0 Then 'Periodically ask the user whether to stop.
UserResponse = MsgBox("Timestep = " & Timestep & " Continue?", vbYesNo)
End If
Loop Until UserResponse = vbNo
End Sub
Sub TriDiagonalSolver (J, A, B, C, R, U)
'Tridiagonal system solver. The input (unchanged) is J, A, B, C, R. The output is array U.
'Integer J, arrays A, B, C, R, U are defined in this picture:
Dim I, Beta, Gamma()
 b (0 ) c (0) 0 0 0   u (0)   r (0 ) 
ReDim Gamma(0 To J)  a (1)   u (1)   r (1) 
 b (1) c (1) 0 0     
Beta = B(0)  0    0        
     
U(0) = R(0) / Beta  0 0 a ( J  1) b ( J  1) c ( J  1)   u ( J  1)   r ( J  1) 
 0 0 0 a (J ) b ( J )   u ( J )   r ( J ) 
For I = 1 To J
Gamma(I) = C(I - 1) / Beta
Beta = B(I) - A(I) * Gamma(I)
U(I) = (R(I) - A(I) * U(I - 1)) / Beta
Next I
For I = J - 1 To 0 Step -1
U(I) = U(I) - Gamma(I + 1) * U(I + 1)
Next I
End Sub