d u(x,t)/dt = C * d^2 u(x,t)/dx^2on 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 forHere 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*LambdaSince 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 norm2completing the proof.