function [t,C_a,C_b] = odesintegrate(solverfun,f0,tinterval,k1,k2); %[t,C_a,C_b] = odesintegrate(solvername,f0,tinterval,k1,k2); %GOALS: %demonstrate the use of the various ode solvers to integrate a pair of %interdependent odes % %show how to pass additional problem parameters to the function that %calculates the derivatives % %DESCRIPTION %this example is a well stirred isothermal batch reactor with three species %and two first order reactions: % k1 k2 % A -> B -> C % %only species A and B are followed as C can be gotten from a mass balance %(C is not important and is not followed) % %VARIABLES %solverfun: time integrator used, i.e. @ode45(explicit) or @ode15s(implicit) %f0: initial values of function pair [C_a_0,C_b_0] %tinterval: [0,tmax] integrates the functions from 0 to tmax % if more than 2 times are given, e.g. [0 t1 t2 t3 ...] then % only the results of specified timesteps are returned, however, % MATLAB may make extra "unseen" steps to maintain accuracy %options: a list of options that set, e.g. the tolerance of the solver. % If [] is given, then the default options are used %k1: scalar that affects rate of reaction 1 %k2: scalar that affects rate of reaction 2 % %EXAMPLE sequence, ramping up k1 with constant k2 %[t,C_a,C_b] = odesintegrate(@ode45,[1,0],[0,6],1,1); length(t) %[t,C_a,C_b] = odesintegrate(@ode45,[1,0],[0,6],10,1); length(t) %[t,C_a,C_b] = odesintegrate(@ode45,[1,0],[0,6],100,1); length(t) %[t,C_a,C_b] = odesintegrate(@ode45,[1,0],[0,6],1000,1); length(t) %[t,C_a,C_b] = odesintegrate(@ode15s,[1,0],[0,6],1000,1); length(t) %BEGIN CODE %define param struct of k1 and k2 param.k1 = k1; param.k2 = k2; %Set tolerances (these are the default tolerances) options = odeset('RelTol',10^-3,'AbsTol',10^-6); %call solver, get back list of t values and f values [t,fs] = solverfun(@odefun,tinterval,f0,options,param); %specific calling sequences, without allowing for any solver %%%%%%[t,fs] = ode45(@odefun,tinterval,f0,options,param); %%%%%%[t,fs] = ode15s(@odefun,tinterval,f0,options,param); %results include a column vector of t values and a matrix of f values %the f matrix is really two horizontally concatenated f column vectors %they can be separated using colon operations C_a = fs(:,1); % <= C_a is the first column vector C_b = fs(:,2); % <= C_b is the second column vector %plot results figure(1); plot(t,C_a,t,C_b); legend('C_a','C_b') xlabel('t') ylabel('f') xlim([0,t(end)]); ylim([0,max([max(C_a),max(C_b)])]); title(['k1 = ',num2str(param.k1),' k2 = ',num2str(param.k2)]) return; function derivs = odefun(t,fs,param); %note that in the context of this function, %a scalar t and a column vector [C_a;C_b] are passed to dfdt %also, everything in param has been passed %switch from list of f's to actual names for ease of formulation of %derivative equations C_a = fs(1); C_b = fs(2); k1 = param.k1; k2 = param.k2; %express derivatives in terms of C_a, C_b, and t dC_adt = -k1*C_a; dC_bdt = k1*C_a - k2*C_b; %put derivative results back in column vector format for MATLAB derivs = [dC_adt;dC_bdt]; return;