# Stability of the Implicit Solution of 1D Heat Equation

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.