%% 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. % This program "schm1_error.m" generates curves for the error in a % series problem as the rate constant ratios for the first and % second first-order reactions are varied. % Last edited 1/30/97. global k1 k2 k1 = 1; k2 = 1.0; Initial = [1,0,0]'; t = [0:.01:1]'; opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps)); [tsolver, solution] = ode15s(@rxrate,t,Initial,opts); answer = [t solution]; k2 = 10.0; Initial = [1,0,0]'; opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps)); [tsolver, solution] = ode15s(@rxrate,t,Initial,opts); answer = [answer solution]; k2 = 50.0; Initial = [1,0,0]'; opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps)); [tsolver, solution] = ode15s(@rxrate,t,Initial,opts); answer = [answer solution]; k2 = 1000.0; Initial = [1,0,0]'; opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps)); [tsolver, solution] = ode15s(@rxrate,t,Initial,opts); answer = [answer solution]; k2 = 10000.0; Initial = [1,0,0]'; opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps)); [tsolver, solution] = ode15s(@rxrate,t,Initial,opts); answer = [answer solution]; k2 = 100000.0; Initial = [1,0,0]'; opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps)); [tsolver, solution] = ode15s(@rxrate,t,Initial,opts); %% Strip out first row to avoid creating NaNs. answer = [answer solution]; answer = answer (2:end, :); t_tmp = t(2:length(t)); c_Css = 1-exp(-k1*t_tmp); answer2 = [t_tmp c_Css]; c_Cexact = answer(:,13); E=(c_Cexact - c_Css)./c_Cexact; EE=abs(E); error=log10(EE); %answer3 = [t_tmp error]; % JBR, 2/22/98 answer3 = [t_tmp EE]; c_Cexact = answer(:,16); E=(c_Cexact - c_Css)./c_Cexact; EE=abs(E); error=log10(EE); %answer3 = [answer3 error]; % JBR, 2/22/98 answer3 = [answer3 EE]; c_Cexact = answer(:,19); E=(c_Cexact - c_Css)./c_Cexact; EE=abs(E); error=log10(EE); %answer3 = [answer3 error]; % JBR, 2/22/98 answer3 = [answer3 EE]; semilogy (answer3(:,1), answer3(:,2:4)); title ('Figure 5.11')

function dcdt = rxrate(t,x) global k1 k2 c1 = x(1); c2 = x(2); c3 = x(3); r1 = k1*c1; r2 = k2*c2; dcdt = [-r1; r1-r2; r2];