1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 | % 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. % % Revised 7/24/2018 p = struct(); % Create structure to pass parameters to ode15s function p.k1 = 1.0; p.k2 = 1.0; Initial = [1,0,0]'; t = [0:.01:1]'; opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps)); [tsolver, solution] = ode15s(@(t,x) rxrate(t,x,p),t,Initial,opts); answer = [t solution]; p.k2 = 10.0; [tsolver, solution] = ode15s(@(t,x) rxrate(t,x,p),t,Initial,opts); answer = [answer solution]; p.k2 = 50.0; [tsolver, solution] = ode15s(@(t,x) rxrate(t,x,p),t,Initial,opts); answer = [answer solution]; p.k2 = 1000.0; [tsolver, solution] = ode15s(@(t,x) rxrate(t,x,p),t,Initial,opts); answer = [answer solution]; p.k2 = 10000.0; [tsolver, solution] = ode15s(@(t,x) rxrate(t,x,p),t,Initial,opts); answer = [answer solution]; p.k2 = 100000.0; [tsolver, solution] = ode15s(@(t,x) rxrate(t,x,p),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(-p.k1*t_tmp); %solution for k2 = Inf 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]; save -ascii schm1_error.dat answer3; if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING semilogy (answer3(:,1), answer3(:,2:4)); % TITLE end % PLOTTING |

1 2 3 4 5 6 7 8 9 | function dcdt = rxrate(t,x,p) c1 = x(1); c2 = x(2); c3 = x(3); r1 = p.k1*c1; r2 = p.k2*c2; dcdt = [-r1; r1-r2; r2]; |