I'm now tired...
------------
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
so
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).
