%% for & while repeat examples % reaction A > B in batch reactor (equil almost all B) % specify well-mixed, isothermal, constant volume % balance on A is dc/dt = -kc % solve using Euler's method of numerical approximation % Euler's is an approximation because it assumes rate of change % (dc/dt here) is constant between steps of independent % variable (t here) % accuracy will improve by taking smaller time step % at cost of longer computation %% case 1 - use "for", specify final time, find conc clear % get rid of any longer arrays from previous runs k = 1; % 1/s, rate coefficient t(1) = 0; % s, initial time c(1) = 1; % mol/m3, initial conc of reactant A tfinal = 1.2; % s, final time nsteps = 12; % number time steps dt = tfinal/nsteps; % s, time step % for repeat automatically initializes and increments % loop variable, which here is used as array subscript for i = 1:nsteps dcdt = -k * c(i); % rate of change in c, dc/dt dc = dcdt * dt; % estimated change in c over time step c(i+1) = c(i) + dc; % estimate next c t(i+1) = t(i) + dt; % get next t end fprintf('--- case 1, use for repeat ------- \n') fprintf('specify final time & number time steps, find final conc \n') fprintf('final time, final = %0.2f vs. specified = %0.2f\n',t(length(t)),tfinal) fprintf('final conc = %0.3f \n',c(length(c))) fprintf('number time steps specified = %i \n',nsteps) fprintf('computed dt = %0.3f \n',dt) figure(1) plot(t,c,'o') % also plot analytical solution hold on ta = 0:0.001:tfinal; ca = c(1)*exp(-k*ta); plot(ta,ca) tt = sprintf('use "for" repeat, blk = analytical, blue circles %i steps',nsteps); title(tt) ylabel('c, concentration of reactant (mol/m^3)') xlabel('t, time (s)') axis([0 1.4 0 1]) hold off %% case 2 - use "while", specify final conc, find time clear % get rid of any longer arrays from previous runs k = 1; % 1/s, rate coefficient t(1) = 0; % s, initial time c(1) = 1; % mol/m3, initial conc of reactant A cfinal = 0.3; % mol/m3, desired final conc of A dt = 0.1; % s, time step % if want to use variable as array subscript with while repeat % then need to initialize before and increment inside repeat i = 1; % initialize array subscript while c(i) > cfinal dcdt = -k * c(i); % rate of change in c, dc/dt dc = dcdt * dt; % estimated change in c over time step c(i+1) = c(i) + dc; % estimate next c t(i+1) = t(i) + dt; % get next t i = i + 1; % increment array subscript end fprintf('--- case 2, use while repeat ------- \n') fprintf('specify final conc and time step size dt, find final time \n') fprintf('final time = %0.2f \n',t(i)) fprintf('final conc, final = %0.3f vs. specified %0.3f \n',c(i),cfinal) fprintf('number time steps taken = %i \n',i-1) fprintf('specified dt = %0.3f \n',dt) fprintf('while will stop closer to cfinal with smaller dt \n') figure(2) plot(t,c,'bo') % also plot analytical solution hold on ta = 0:0.001:t(i); ca = c(1)*exp(-k*ta); plot(ta,ca,'k') tt = sprintf('use "while" repeat, blk = analytical, blue circles %i steps',i-1); title(tt) ylabel('c, concentration of reactant (mol/m^3)') xlabel('t, time (s)') axis([0 1.4 0 1]) hold off fprintf('--------------------- \n')