clear
reset(0)
echo on
%In this example we demonstrate the optimization of
%a function using Newton's method, successive parabolic
%interpolation, and the golden search (modified Fibonacci
%search).  Optimization can mean finding either a
%minimum or a maximum.  In this problem we look for
%the maximum.
%
pause
%Suppose we have the sequential batch reaction:
%
%  A --> B --> C
%
%In this problem B is the desired product, and C is
%an undesirable byproduct.  We want to maximize the
%yield of B.  The two reactions are assumed to be first
%order:
pause
%
%  Ra = -k1 * Ca
%
%  Rb = k1 * Ca - k2 * Cb
%
%where k1 and k2 are rate constants.  If the initial
%concentration of B is zero, and the initial concentration
%of A is Ca0, we have the solution:
%
%  Cb = Ca0*k1/(k2-k1)*(exp(-k1*t)-exp(-k2*t))
%
%which initially increases and decreases with time.  We
%have put this function into the function cb.m.
pause
%Let's plot this up:
%We pick values of Ca0, k1 and k2, and let the function
%cb know what they are:
global ca0 k1 k2
ca0=1;
k1=2;
k2=1;
tt=[0:.1:2];
plot(tt,cb(tt))
set(gca,'FontSize',18)
xlabel('time')
ylabel('concentration of B')
hold on
pause
%As you can see, there is a maximum after about
% t = 0.7 or so.  Let's use the three optimization
%routines discussed in class to find this value more
%accurately.
pause
%First, we use Newton's method starting from t = 1:
%We have the formula:
%
%  x(k+1) = x(k) - f'(x(k))/f''(x(k))
%
%For simplicity, we shall take the derivatives numerically
%using second order difference algorithms:
pause
eps=1e-5;
clear t
t(1)=1;
i=1;
cbc=cb(t(i));
plot(t(i),cbc,'oc')
set(gca,'FontSize',18)
fp=1;
dt=1;
while abs(dt)>1e-4
 cbp=cb(t(i)+eps);
 cbm=cb(t(i)-eps);
 fp=(cbp-cbm)/2/eps;
 fpp=(cbp+cbm-2*cbc)/eps^2;
 dt=-fp/fpp;
 t(i+1)=t(i)+dt;
 i=i+1;
 cbc=cb(t(i));
 plot(t(i),cbc,'ow')
 pause(1)
 echo off
end
echo on
hold off
pause
%The exact solution for the maximum is just:
%
texact=log(k2/k1)/(k2-k1)
%
%thus we can calculate the rate of progression to
%the correct solution:
ernm=abs(t-texact);
nnm=length(ernm);
loglog(ernm(1:nnm-1),ernm(2:nnm))
set(gca,'FontSize',18)
xlabel('error(i)')
ylabel('error(i+1)')
pause
%So you can see that the Newton's method converges
%quadratically for this function.  Note that this 
%works only as long as we can determine the first
%and second derivatives with sufficient accuracy.
%Also, we are required to start where the second
%derivative is of the appropriate sign, in this case
%negative.  If we started at t=2 we would have instead
%tried to capture the minimum at t --> infinity.
pause
%Now let's look at the method of successive parabolic
%interpolation.  At each step we fit a parabola to
%three points, and use the parabola to guess where the
%maximum or the minimum is.  Thus we begin with:
clear t
t(1)=1;
t(2)=0;
t(3)=.5;
f=cb(t);
i=3;
plot(tt,cb(tt))
set(gca,'FontSize',18)
xlabel('time')
ylabel('concentration of B')
hold on
plot(t,f,'oc')
pause
while abs(t(i)-t(i-1))>1e-7,
 tm=t(i-2:i)';
 fm=f(i-2:i)';
 a=[ones(3,1),tm,tm.^2];
 x=a\fm;
 i=i+1;
 t(i)=-x(2)/2/x(3);
 f(i)=cb(t(i));
 plot(t(i),f(i),'ow')
 pause(1)
 echo off
end
echo on
hold off
pause
%We can calculate the rate of progression to the exact
%solution:
erspi=abs(t-texact);
nspi=length(erspi);
loglog(ernm(1:nnm-1),ernm(2:nnm))
set(gca,'FontSize',18)
xlabel('error(i)')
ylabel('error(i+1)')
hold on
loglog(erspi(1:nspi-1),erspi(2:nspi),'g')
hold off
pause
%In this case the progression is only slightly super-
%linear, with a slope of 1.324.  The progression is
%also a bit more random.
pause
%To conclude, we shall finish with the Fibonacci or
%golden search.  This uses the golden ratio which is
%the solution of:
%
%  r^2 = 1 - r
%
%or r = 0.6180 or so.  This search only works on unimodal
%functions, those functions for which the first derivative
%is strictly positive for t < t* and negative for t > t*,
%or the reverse if you are looking for a minimum.  Note
%that the second derivative can locally change sign and
%this method will still behave.  Newton's method would
%fail in this circumstance.
pause
%Let's work out this algorithm.  We start with the
%some known interval over which the function is unimodal,
%say [0,1].  We evaluate the function in the interior at
%two additional points, located r and 1-r from the edges.
%We keep as our new interval either [0,r] or [1-r,1]
%depending on which side contains the maximum (or the
%minimum if that's what you are after).  We then
%proceed to the next step.  The key is that only one
%new function evaluation is required at each step.
pause
%Let's do this:
r=(sqrt(5)-1)/2;
clear t
plot(tt,cb(tt))
set(gca,'FontSize',18)
xlabel('time')
ylabel('concentration of B')
hold on
a=0;
b=1;
t(1)=a+(b-a)*(1-r);
tl=t(1);
t(2)=a+(b-a)*r;
tr=t(2);
f=cb(t);
plot(t,f,'oc')
fl=f(1);
fr=f(2);
i=2;
pause
while abs(t(i)-t(i-1))>1e-5,
 i=i+1;
 if fr>fl,
  a=tl;
  tl=tr;
  fl=fr;
  t(i)=a+(b-a)*r;
  tr=t(i);
  fr=cb(tr);
  plot(tr,fr,'ow')
 else,
  b=tr;
  tr=tl;
  fr=fl;
  t(i)=a+(b-a)*(1-r);
  tl=t(i);
  fl=cb(tl);
  plot(tl,fl,'ow')
 end
 pause(1)
 echo off
end
echo on
hold off
pause
%This technique converges even more slowly than the
%successive parabolic interpolation method.  The rate
%of convergence is linear, and the interval is reduced
%by a factor of 0.6180 at each step.  We can plot this:
pause
ergs=abs(t-texact);
ngs=length(ergs);
loglog(ernm(1:nnm-1),ernm(2:nnm))
set(gca,'FontSize',18)
xlabel('error(i)')
ylabel('error(i+1)')
hold on
loglog(erspi(1:nspi-1),erspi(2:nspi),'g')
loglog(ergs(1:ngs-1),ergs(2:ngs),'w')
hold off
pause
%A comparison of the three techniques for this function
%shows that Newton's method converges the fastest, but
%only if you start close to the correct value.  The
%golden search is slower, but more robust - we could
%have chosen any interval containing the maximum and
%we would have gotten there eventually.  The same cannot
%be said of the two superlinear techniques.
echo off

__________________________


function y=cb(t)
%This function takes it an array of times and returns
%the concentration of B.  We need the reaction rates
%given by a global statement in the main script:

global ca0 k1 k2

y=ca0*k1/(k2-k1)*(exp(-k1*t)-exp(-k2*t));