clear
reset(0)
echo on
%In this program we demonstrate the instability of
%explicit techniques for integrating a stiff
%equation.  We shall look at the behavior of the 
%function dy/dt = f(y) where f = -a*y+b.  This is
%the differential equation which would result from
%a steady filling of a CSTR with salt water if it
%initially contained fresh.  The variable y would
%thus represent the amount of salt present.  The 
%stability depends on the Jacobian J = df/dy, which
%in this case is negative (at least for positive
%values of the parameter a).  Let's look at the 
%stability of the Euler Method.  To begin, we take
%a=5, b=5 and y(0)=0.  The solution to this problem
%is just y=(b/a)*(1-exp(-a*t)).  Thus the function
%looks like:
a=5;
b=5;
x=[0:.01:3];
y=(b/a)*(1-exp(-a*x));
v=[0,3,-.5,2.5];
axis(v)
plot(x,y)
hold on
pause

%Now let's look at what the Euler Method integration
%returns for a result.  Remember that the EM rule
%just locally fits the function with a line tangent
%to the function at a point t(k),y(k) and uses the line
%to predict the value at t(k+1).  The expression
%is thus yem(k+1)=yem(k)+f(yem(k),t(k))*dt.  We thus
%get the series:
dt=.1;
yem(1)=0;
t(1)=0;
plot(t(1),yem(1),'o')
pause(1)
k=1;
while t(k)<3,
  yem(k+1)=yem(k)+(b-yem(k)*a)*dt;
  t(k+1)=t(k)+dt;
  plot(t(k+1),yem(k+1),'o')
  pause(1)
  k=k+1;
  echo off
end;
echo on
%As you can see, the function integration procedure
%does pretty well in tracking the function.  This is
%because the function itself is stable, and the
%step size is small enough that the product hJ 
%satisfies the stability bound |1+hJ| < 1.  We can
%plot the error:
pause
hold off
plot(t,yem-(b/a)*(1-exp(-a*t)))
%From this plot we see that the error starts off at
%zero, increases in magnitude, and then decreases
%again.  The later decrease is because the equation
%is stable - the family of solution curves is
%converging - and hence the earlier errors are
%being wiped out.  Let's contrast this behavior
%with what would happen if we used a larger step
%size.  Now we pick dt = .3, which is still not quite
%enough that we are unstable.
pause
axis(v)
dt=.3;
clear t;
clear yem;
plot(x,y)
hold on
yem(1)=0;
t(1)=0;
plot(t(1),yem(1),'o')
pause(1)
k=1;
while t(k)<3,
  yem(k+1)=yem(k)+(b-yem(k)*a)*dt;
  t(k+1)=t(k)+dt;
  plot(t(k+1),yem(k+1),'o')
  pause(1)
  k=k+1;
  echo off
end;
echo on
%This time the error oscillates before it settles
%down.  The method is almost unstable, but here
%the product hJ is on the order of -1.5, still less
%than the critical value of -2.  If we increase dt
%just a bit more, bad things happen:
pause
hold off
dt=.4;
clear t;
clear yem;
axis(v)
plot(x,y)
hold on
yem(1)=0;
t(1)=0;
plot(t(1),yem(1),'o')
pause(1)
k=1;
while t(k)<3,
  yem(k+1)=yem(k)+(b-yem(k)*a)*dt;
  t(k+1)=t(k)+dt;
  plot(t(k+1),yem(k+1),'o')
  pause(1)
  k=k+1;
  echo off
end;
echo on
hold off
pause
%Now let's look at how we can improve things using
%an implicit technique.  The simplest technique is
%the Backward Euler method, in which we approximate
%the function with a line, this time tangent at
%t(k+1),y(k+1) and use the line drawn from the point
%t(k),y(k) to predict y(k+1).  This leads to the
%BE equation: ybe(k+1)=ybe(k)+f(ybe(k+1),t(k+1))*dt.
%Normally for non-linear functions f(y,t) we have to
%use a non-linear equation solver to solve the
%implicit relationship in the BE algorithm.  For a
%linear system, however, this is not the case.  We
%can solve for ybe(k+1) explicitly.  Thus we get:
dt=.1;
clear t;
clear ybe;
axis(v)
plot(x,y)
hold on
ybe(1)=0;
t(1)=0;
plot(t(1),ybe(1),'o')
pause(1)
k=1;
while t(k)<3,
  ybe(k+1)=(ybe(k)+b*dt)/(1+a*dt);
  t(k+1)=t(k)+dt;
  plot(t(k+1),ybe(k+1),'o')
  pause(1)
  k=k+1;
  echo off
end;
echo on
hold off
pause
%The error in this case is just about the same as
%the normal EM case.  This is because we picked dt
%to be small enough (0.1) so that hJ wasn't large
%enough to cause instability.  Now let's use dt = .4,
%the value that caused the EM to die:
dt=.4;
clear t;
clear ybe;
axis(v)
plot(x,y)
hold on
ybe(1)=0;
t(1)=0;
plot(t(1),ybe(1),'o')
pause(1)
k=1;
while t(k)<3,
  ybe(k+1)=(ybe(k)+b*dt)/(1+a*dt);
  t(k+1)=t(k)+dt;
  plot(t(k+1),ybe(k+1),'o')
  pause(1)
  k=k+1;
  echo off
end;
echo on
hold off
pause
%In this case the error is larger than the dt=0.1
%example (as would be expected), but the oscillations
%due to the EM algorithm are completely wiped out.
%This reflects the fact that the BE implicit method
%is absolutely stable for all stable differential
%equations, no matter how large the step size.  This
%is why most commercial integration routines have
%at least some implicit character to them despite
%the difficulty of solving the non-linear equation
%posed by the implicit problem.  Note that this
%non-linear equation is usually very easy to solve
%because you have an -excellent- initial guess for
%the value of y(k+1) from explicit formulas.  Thus,
%rapid schemes such as Newton's method have a high
%rate of convergence and don't usually run into 
%trouble.
pause
hold off
echo off