Stopping the integration of an ODE at some condition
In Post 968 we learned how to get the numerical solution to an ODE, and then to use the deval function to solve the solution for a particular value. The deval function uses interpolation to evaluate the solution at other valuse. An alternative approach would be to stop the ODE integration when the solution has the value you want. That can be done in Matlab by using an "event". You setup an event function and tell the ode solver to use it by setting an option.
given that the concentration of a species A in a constant volume, batch reactor obeys this differential equation with the initial condition mol/L and L/mol/s, compute the time it takes for to be reduced to 1 mol/L.
clear all; close all; clc k = 0.23; Cao = 2.3; tspan = [0 10]; % create an options variable options = odeset('Events',@event_function); % note the extra outputs [t,Ca,TE,VE] = ode45(@myode,tspan,Cao,options);
TE % this is the time value where the event occurred VE % this is the value of Ca where the event occurred sprintf('At t = %1.2f seconds the concentration of A is %1.2f mol/L.', TE,VE)
TE = 2.4598 VE = 1.0000 ans = At t = 2.46 seconds the concentration of A is 1.00 mol/L.
plot(t,Ca) xlabel('Time (sec)') ylabel('C_A (mol/L)') % you can see the integration terminated at Ca = 1
ans = done
function dCadt = myode(t,Ca) k = 0.23; dCadt = -k*Ca^2; function [value,isterminal,direction] = event_function(t,Ca)
% when value is equal to zero, an event is triggered. % set isterminal to 1 to stop the solver at the first event, or 0 to % get all the events. % direction=0 if all zeros are to be computed (the default), +1 if % only zeros where the event function is increasing, and -1 if only % zeros where the event function is decreasing. value = Ca - 1.0; % when value = 0, an event is triggered isterminal = 1; % terminate after the first event direction = 0; % get all the zeros
% categories: ODEs % tags: reaction engineering