Stability of Explicit Solution of 1D Heat Equation

(CS 267, Mar 14 1995)

We consider the following explicit algorithm for solving the heat equation

            d u(x,t)/dt = C * d^2 u(x,t)/dx^2

on a 1D mesh from 0 to 1, with constant conductivity C, zero boundary conditions u(1,t)=u(0,t)=0, mesh spacing h=1/n, time step k, and z=k/h^2. We let U(i,m) denote our approximation of u(i*h,m*k).

     Explicit solution of the 1D Heat Equation

     for m = 1 to final m
         for i = 1 to n-1
              U(i,m+1) =  z * U(i-1,m) + (1-2*z) * U(i,m) + z * U(i+1,m)
         end for
     end for

Let norm(U(:,m)) = max_i |U(i,m)|. We make two claims:

When 0 < z <= .5, norm(U(:,m+1)) <= norm(U(:,m)). In other words, the maximum absolute temperature can only decrease.
When z > .5, for any eps>0 there exists some n and some initial data U(:,0) with norm(U(:,0)) < eps such that norm(U(:,m)) becomes arbitrarily large as m grows.

In case 1, the algorithm is called stable; stability is necessary to get an accurate, physically meaningful result. In case 2, the algorithm is called unstable; instability yields an inaccurate and physically meaningless result.

Proof of Case 1. When 0 < z <= .5, then 0 <= 1-2*z < 1, so

   |U(i,m+1)| <= z*|U(i-1,m)| + (1-2*z)*|U(i,m)| + z*|U(i+1,m)|
              <= z*norm(U(:,m)) + (1-2*z)*norm(U(:,m)) + z*norm(U(:,m))
              =  norm(U(:,m))

Since this is true for all i, the result follows. This completes the proof of Case 1.

Proof of Case 2. Let U(i,0)=eps*sin((n-1)*i*pi/n), so norm(U(:,0)) <= eps. Then a little bit of trigonometry shows that

    U(i,1) = eps*sin((n-1)*i*pi/n)* ( 1+2*z*(cos((n-1)*pi/n)-1) )
           = U(i,0) * ( 1+2*z*(cos((n-1)*pi/n)-1) )

After m iterations, we get

    U(i,m) = U(i,0) * ( 1-2*z + 2*z*cos((n-1)*pi/n) )^m

    norm(U(:,m)) = norm(U(:,0)) * ( 1-2*z + 2*z*cos((n-1)*pi/n) )^m

The proof will be complete if we show that

      1-2*z + 2*z*cos((n-1)*pi/n) < -1

since at each stage the solution U(:,m) will be multiplied by a number exceeding 1 in absolute value, and so grow without bound. After a little algebra, we get

      1 - cos((n-1)*pi/n) > 1/z

Since z > .5, 1/z < 2. As n grows, cos((n-1)*pi/n) approaches cos(pi) = -1, so 1 - cos((n-1)*pi/n) approaches 2, and so eventually exceeds 1/z. This completes the proof of Case 2.

The background behind the Proof of Case 2 is this. We can express the inner loop of the algorithm in one line as

       U(:,m+1) = T(z) * U(:,m)

where U(:,m+1) and U(:,m) are (n-1)-by-1 vectors, and T(z) is an (n-1)-by-(n-1) tridiagonal matrix:

              [ 1-2*z    z                 ]
              [   z    1-2*z    z          ]
       T(z) = [            ...             ]
              [          z    1-2*z    z   ]
              [                 z    1-2*z ]

One can use trigonometry to confirm that the eigenvalues and eigenvectors of T(z), T(z)*x(i) = lambda(i)*x(i), are given by

      lambda(i) = (1-2*z) + 2*z*cos(i*pi/n)
  and
      x(i,j) = sqrt(2/n) * sin(i*j*pi/n)

The sqrt(2/n) factor guarantees that each eigenvector has unit length. In case 2, the largest eigenvalue in absolute value is lambda(n-1). We chose U(:,0) proportional to the corresponding eigenvector x(n-1), so at each step of the algorithm, U(:,m) just gets multiplied by the eigenvalue lambda(n-1).