clear
reset(0)
echo on
%Example 31a
%
%Numerical stability problems often occur in the
%solution of partial differential equations.  Here
%we illustrate this with the example of the quenching
%of a finite slab using a very simple numerical method.
%
%Consider a slab of width 2b initially at a temperature T1,
%and then quench the temperature at the surface to a value
%T0.  The thermal diffusivity alpha is assumed to be a
%constant.  We render x dimensionless with respect to the
%half-width b, and time dimensionless with respect to the
%diffusion time b^2/alpha.  Thus we have the equation:
%
% dT/dt = d2T/dx2; T(t=0)=1, T(x=1)=0, dT/dx(x=0)=0
%
%The last is the symmetry condition at the centerline.
pause
%OK, now we solve this as a system of first order ODE's.
%We discretize x:
n=10;
dx=1/n;
x=[0:dx:1];

%and we also discretize T:
T=ones(n+1,1);
pause
%Now we solve this as a function of time.  We shall use
%the Euler method for simplicity:
%
%We construct the tri-diagonal matrix governing the
%discretized equation:
a=-2*eye(n+1);
a=a+diag(ones(1,n),-1);
a=a+diag(ones(1,n),+1);
a=a/dx/dx;
%And we modify the first and last rows to account for the
%boundary conditions:
a(1,:)=zeros(1,n+1);
a(n+1,:)=zeros(1,n+1);
pause
%We can now integrate forward in time:
t=0;
dt=0.005;
T=ones(n+1,1);
T(n+1)=0;
figure(1)
plot(x,T)
drawnow
for k=1:1/dt
 T=T+dt*a*T;
 T(n+1)=0;             %This is the zero Temperature BC at x=1
 T(1)=(4*T(2)-T(3))/3; %This is the zero derivative BC at x=0
 plot(x,T)
 axis([0 1 0 1])
 xlabel('x')
 ylabel('T')
 t=t+dt;
 drawnow
 echo off
end
echo on
pause
%The integration procedure behaved very well.  Now let's do
%it again with a -very- slightly larger time step:
dt=0.0052;
pause
t=0;
T=ones(n+1,1);
T(n+1)=0;
figure(1)
plot(x,T)
drawnow
for k=1:1/dt
 T=T+dt*a*T;
 T(n+1)=0;
 T(1)=(4*T(2)-T(3))/3;
 plot(x,T)
 axis([0 1 0 1])
 xlabel('x')
 ylabel('T')
 t=t+dt;
 drawnow
 echo off
end
echo on
pause
%The dramatic difference between these two results illustrates
%the problem of numerical instability!  We require dt < 0.5 dx^2
%for the method to be stable.  Any explicit integration method
%would have the same stability limitations.  We can get around 
%this limitation by using implicit integration techniques.
pause
%We finish this example by using the implicit Backward Euler Rule
%to integrate this problem forward.  First let's double dt:
dt=dt*2
%This would be *really* unstable for an explicit method!  Next
%we construct the matrices:
al=eye(n+1)-dt*a;
ar=eye(n+1);
pause
%We also have to incorporate the B.C.'s into the matrices:
al(1,:)=zeros(1,n+1);
al(1,1)=-3;
al(1,2)=4;
al(1,3)=-1;
ar(1,:)=zeros(1,n+1);
%and for the other B.C.
ar(n+1,:)=zeros(1,n+1);
%and we get the solution:
pause
t=0;
T=ones(n+1,1);
T(n+1)=0;
figure(1)
plot(x,T)
drawnow
for k=1:1/dt
 T=al\(ar*T);
 plot(x,T)
 axis([0 1 0 1])
 xlabel('x')
 ylabel('T')
 t=t+dt;
 drawnow
 echo off
end
echo on
pause
%Note that this is now stable even for much larger values of dt!
echo off