% listing of main file "kin3.m" % this file calls the file "kin3F.m" % and uses the standard matlab function fminsearch % The input data are the weight of catalyst (W), volumetric flow % rate (v), input concentrations (cAin cBin) and concentrations % in the reactor (cA cB) for a CSTR. % These data are fit to a rate equation that is specified in the % file "kin3F.m" that contains three rate parameters whose % final values are returned in the array "result" below. % In this example file "kin3F.m", the rate equation specified is: % rcalc = (kf * cA - kr * cB) ./ (1 + KA * cA + KB * cB) % declare global variables used by kin3.m and kin3F.m % run separator fprintf('------------------------------------------------\n') global rexp cA cB % input data W = 100; % 100 g of catalyst v = 10; % 10 liter/hr flow rate cAin = [0.559;1.54;1.13;0.755;0.164;0.326;1.36;0.515;0.596;0.907]; cBin = [0.761;0.035;0.435;0.185;0.186;0.794;0.0488;0.385;0.724;0.0633]; cA = [0.42;0.72;0.62;0.37;0.11;0.32;0.64;0.31;0.43;0.41]; cB = [0.9;0.85;0.95;0.57;0.24;0.8;0.77;0.59;0.89;0.56]; % calculate "experimental" rates = -rA(mol/g hr) % dot . NOT needed before the * here since (v/W) is scalar rexp = (v/W) * (cAin - cA); g = [1;1;1;1]; % initial guesses of kf, kr, KA, KB % do nonlinear multi-variable fitting using fminsearch result = fminsearch('kin3F',g) % fminsearch is a standard MATLAB function % kin3F is user-written function in file kin3F.m % g is the vector of initial guesses set above fprintf('row 1 of result is kf value, row 2 is kr, row 3 is KA, row 4 is KB \n\n') sumSquaredError = kin3F(result) fprintf('SEE PLOT \n\n') % plot kf = result(1); kr = result(2); KA = result(3); KB = result(4); rcalc = (kf * cA - kr * cB) ./ (1 + KA * cA + KB * cB); % get indices of max and min values of rexp % to use in drawing diagonal line % (or just plot rexp vs. rexp) mini = find(rexp == min(rexp)); % indices of min value maxi = find(rexp == max(rexp)); % indices of max value rdiag = [rexp(mini),rexp(maxi)]; plot(rexp,rcalc,'bo',rdiag,rdiag,'k') title('kin3 problem') ylabel('rcalc (mol/g/hr)') xlabel('rexp (mol/g/hr)')