Stability of Implicit Solution of 1D Heat Equation

(CS 267, Mar 14 1995)

We consider the Crank Nicholson 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,j*k).

     Crank-Nicholson solution of the 1D Heat Equation

     for m = 1 to final m
         Solve ( I+(z/2)*T)*U(:,m+1) = ( I - (z/2)*T ) * U(:,m) for U(:,m+1)
     end for

Here T is the (n-1)-by-(n-1) tridiagonal matrix

         [  2 -1          ]
         [ -1  2 -1       ]
     T = [     ...        ]
         [       -1  2 -1 ]
         [          -1  2 ]

Let norm2(U(:,m)) = sum_{i=1}^{n-1} |U(i,m)|. We claim that for all z>0, norm2(U(:,m+1)) <= norm2(U(:,m)), so the norm of the absolute temperature can only decrease, independent of time step k.

Proof. We can write

        U(:,m+1) = inv( I+(z/2)*T) ) * ( I - (z/2)*T ) * U(:,m)
                 = X * U(:,m)

So each step is equivalent to matrix-vector multiplication by the matrix X = inv( I+(z/2)*T) ) * ( I - (z/2)*T ). We will compute the eigenvalues and eigenvectors of X. Since T is a symmetric matrix, it has real eigenvalues and orthogonal eigenvectors. Using trigonometry, on can confirm that they are given by lambda(i) and v(i), where T*v(i) = lambda(i)*v(i) and

      lambda(i) = 2 - 2*cos(i*pi/n)
      v(i,j) = sqrt(2/n) * sin(i*j*pi/n)

Note that all the eigenvalues lie strictly between 0 and 4. Letting Lambda = diag(lambda(1), ... , lambda(n-1)) and V(i,j) = v(i,j), we may write

     T*V = V*Lambda

Since the v(i) are orthgonal (and and unit length), one can verify that V*V'=I. so

     T = V*Lambda*V'

Substituting this in X = inv( I+(z/2)*T) ) * ( I - (z/2)*T ) yields

     X = inv( I + (z/2)*V*Lambda*V' ) * ( I - (z/2)*V*Lambda*V' )
       = inv(V*(I + (z/2)*Lambda)*V') * V*(I - (z/2)*Lambda)*V'
       = V*inv(I+(z/2)*Lambda)*V'*V*(I-(z/2)*Lambda)*V'
       = V*inv(I+(z/2)*Lambda)*(I-(z/2)*Lambda)*V'
       = V*diag( (1-(z/2)*lambda(i))/(1+(z/2)*lambda(i)))*V'

So we see X has the same eigenvectorx v(i) as T and eigenvalues

       x(i) = (1-(z/2)*lambda(i))/(1+(z/2)*lambda(i))

Since 0 < lambda(i), one can confirm that -1 < x(i) < 1. Multiplying a vector by V essentially just rotates the vector, it does not change its length so norm2(s) = norm2(V*s) = norm2(V'*s) for any vector s. Thus

      norm2(U(:,m+1)) = norm2(X*U(:,m))
                      = norm2(V*diag(x(i))*V'*U(:,m))
                      = norm2(diag(x(i))*V'*U(:,m))     
                           ... since V*s and s have the same norm2
                      < norm2(V'*U(:,m))                
                           ... since each |x(i)| < 1
                      = norm2(U(:,m))                   
                           ... since V'*s and S ahve the same norm2

completing the proof.