Figure 9.16 (page 527):

Model fit to all adsorption experiments.

Code for Figure 9.16

Text of the GNU GPL.

main.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

%%
%% jbr,  11/28/01
%%
global X Y grad Ypred dYdtheta lb ub

Hgas1 = [2.659, 27.590, 141.038, ...
	332.613, 546.857, 608.583, 875.353, 991.072, 1080.134, ...
	1142.782, 1203.178, 1256.244]';

Hads1 = [32.31, 40.02, 45.40, 48.49 ...
	50.19, 53.86, 52.77, 53.19, 53.67, 54.15, 54.51, 54.73]';


Hgas2 = [3.755, 36.331, 189.904, ...
	 413.881, 624.665, 817.592]';
	
Hads2 = [33.17, 40.90, 45.59, 47.72, ...
	 48.74, 49.43]';

Hgas3 = [12.794, 119.331, 339.917, ...
	 545.227, 735.746, 890.642]';

Hads3 = [36.91, 42.90, 45.65, 47.00, ...
	 47.82, 48.60]';

Hgas4 = [7.424, 73.312, 251.177, ...
	 491.691, 707.806, 873.926]';

Hads4 = [39.11, 46.06, 49.87, 51.57, ...
	 52.48, 53.18]';

Hgas5 = [9.176, 72.065, 259.002, ...
	 475.250, 670.034, 818.669]';

Hads5 = [39.52, 46.06, 49.92, 51.86, ...
	 52.96, 53.88]';
	
Hgas6 = [2.962, 34.350, 176.860, 405.233, ...
	 607.476, 799.723, 946.765]';

Hads6 = [32.70, 40.11, 44.92, 47.40, 48.62, ...
	 49.50, 50.12]';

Hgas7 = [3.159, 38.528, 200.877, 433.924, ...
	 647.630, 820.673]';

Hads7 = [32.42, 39.81, 44.51, 46.79, 48.14, ...
	 48.95]';

Hgas8 = [34.127, 172.316, ...
	 385.633, 616.149, 812.932, 961.191, 1061.967, 1155.829, 1219.299]';

Hads8 = [43.79, 48.97, 51.57, ...
	 53.09, 54.05, 54.67, 55.20, 55.60, 55.93]';
%%
%% choose a dataset or datasets to fit
%%
%%X = [Hgas4];
%%Y = [Hads4];
X = [Hgas1; Hgas2; Hgas3; Hgas4; Hgas5; Hgas6; Hgas7; Hgas8];
Y = [Hads1; Hads2; Hads3; Hads4; Hads5; Hads6; Hads7; Hads8];

Ks0 = 2;
cms0 = 40;
cmx0 = 0;
theta0 = [Ks0; cms0; cmx0];



small = 1e-5;
large = 1000;
lb = [small; small; small];
ub = [large; large; large];


options = optimset('GradObj','on');
[theta, phi, info] = fmincon ({@model,@gradient}, theta0, ...
        		[], [], [], [], lb, ub, [], options);

info;
obj = model(theta);
grad = gradient(theta);
nplot = 100;
Xplot = linspace(0,max(X),nplot)';
Ks = theta(1);
cms = theta(2);
cmx = theta(3);
Yplot = cms*Ks*sqrt(Xplot)./(1+Ks*sqrt(Xplot)) + cmx;
%plot(Hgas1, Hads1, '@', Hgas2, Hads2, '@',  Hgas3, Hads3, '@', ...
%     Hgas4, Hads4, '@',  Hgas5, Hads5, '@', Hgas6, Hads6, '@', ...
%     Hgas7, Hads7, '@', Hgas8, Hads8, '@', Xplot,Yplot)
Hgn = 2*dYdtheta'*dYdtheta;
alpha = 0.95;
p     = length(theta);
ndata = length(X);
Fstat = finv(alpha,p,ndata-p);
samplevar = obj/(ndata-p);
b     = 2*p*Fstat*samplevar;
inva  = inv(Hgn);

theta
bbox  = sqrt(b*diag(inva))

table1 = [Hgas1 Hads1];
table2 = [Hgas2 Hads2];
table3 = [Hgas3 Hads3];
table4 = [Hgas4 Hads4];
table5 = [Hgas5 Hads5];
table6 = [Hgas6 Hads6];
table7 = [Hgas7 Hads7];
table8 = [Hgas8 Hads8];
table_all = [Xplot Yplot];
save -ascii adsall.dat table_all;

plot (table_all(:,1), table_all(:,2), ...
       table1(:,1), table1(:,2), '+', ...
       table2(:,1), table2(:,2), 'x', ...
       table3(:,1), table3(:,2), '*', ...
       table4(:,1), table4(:,2), 'o', ...
       table5(:,1), table5(:,2), '+', ...
       table6(:,1), table6(:,2), 'x', ...
       table7(:,1), table7(:,2), '*', ...
       table8(:,1), table8(:,2), 'o');
title ('Figure 9.16')

%  [x, y, major, minor, box] = ellipse (Hgn, b, 150);
%%
%% check a finite difference Hessian
%%
%  function Hfd = findif(theta, del)
%  p=length(theta);
%  for i = 1:p
%    thetap = theta;
%    thetam = theta;
%    thetap(i) = theta(i)*(1+del);
%    thetam(i) = theta(i)*(1-del);
%    obj = model(thetap);
%    gradp = gradient(thetap);
%    obj = model(thetam);
%    gradm = gradient(thetam);
%    Hfd(1:p,i) = (gradp - gradm)/(2*del*theta(i));
%  end
%  endfunction
%  del= 1e-4
%  Hfd = findif(theta, del)

model.m

function obj = model(theta)
  global X Y grad Ypred dYdtheta
  Ks = theta(1);
  cms = theta(2);
  cmx = theta(3);
  %%
  %% choose associative or dissociative adsorption
  %%
%  Ypred = cm*K*X./(1+K*X);
  Ypred = cms*Ks*sqrt(X)./(1+Ks*sqrt(X)) + cmx;
  resid = Y - Ypred;
  obj = resid'*resid;
  %%
  %% while here, compute the gradient and store in global
  %%
%  dYdtheta = [Ypred.*(1/Ks - X./(1+Ks*X)), Ypred/cms];
  dYdtheta =  ...
[ cms*Ks*sqrt(X)./(1+Ks*sqrt(X)).*(1/Ks - sqrt(X)./(1+Ks*sqrt(X))), ...
	Ks*sqrt(X)./(1+Ks*sqrt(X)), ones(size(Ypred)) ];
  grad = -2*dYdtheta'*resid;

gradient.m

function retval = gradient(theta)
global X Y grad Ypred dYdtheta lb ub

  %%
  %% grad was just computed in model so take it from 
  %% global and return
  %%
  obj = model(theta);
  retval = grad;

hineq.m

function retval = hineq(theta)
%% provides h_ineq(x) for inequality constraint h_ineq(x) >= 0
%% to ensure that parameters are between lower and upper bound
%% parlb <= theta <= parub
  global lb ub
  retval = [theta - lb; -theta + ub];

../../util/common/ellipse.m

%% Copyright (C) 2001, James B. Rawlings and John W. Eaton
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

%% [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)
%%
%% Given a 2x2 matrix, generate ellipse data for plotting.  The
%% arguments N and SHIFT are optional.  If N is an empty matrix, a
%% default value of 100 is used.

function [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)

  if (nargin < 3)
    n = 100;
  end

  if (isempty (n))
    n = 100;
  end

  if (nargin < 4)
    shift = [0, 0];
  end

  ss = size (shift);

  if (any (ss ~= [1, 2]))
    if (ss == [2, 1])
      shift = shift';
    else
      error ('shift must be a 2-element row vector');
    end
  end

  if (nargin > 1)

    [v, l] = eig (amat / level);

    dl = diag(l);
    if (any (imag (dl)) || any (dl <= 0))
      error ('ellipse: amat must be positive definite');
    end

    %% Generate contour data.

    a = 1 / sqrt (l(1,1));
    b = 1 / sqrt (l(2,2));

    t = linspace (0, 2*pi, n)';

    xt = a * cos (t);
    yt = b * sin (t);

    %% Rotate the contours.

    ra = atan2 (v(2,1), v(1,1));

    cos_ra = cos (ra);
    sin_ra = sin (ra);

    x = xt * cos_ra - yt * sin_ra + shift(1);
    y = xt * sin_ra + yt * cos_ra + shift(2);

    %% Endpoints of the major and minor axes.

    minor = (v * diag ([a, b]))';
    major = minor;

    major(2,:) = -major(1,:);
    minor(1,:) = -minor(2,:);

    t = [1; 1] * shift;

    major = major + t;
    minor = minor + t;

    %% Bounding box for the ellipse using magic formula.

    ainv = inv (amat);
    xbox = sqrt (level * ainv(1,1));
    ybox = sqrt (level * ainv(2,2));

    bbox = [xbox ybox; xbox -ybox; -xbox -ybox; -xbox ybox; xbox ybox];

    t = [1; 1; 1; 1; 1] * shift;
    bbox = bbox + t;

  else
    error ('usage: ellipse (amat, level, n, shift)');
  end

%% endfunction