8.4. MATLAB Program 4#
%Solve a system of three 1st order ODEs, using a stiff solver, ode15s,
% from the MATLAB ODE solver routines.
%%1-> Define the IVP
% declare an ode function (definition given at the end of the code)
f=@odeFunc;
%set the initial values
y0= [1.0 0.0 0.0];
%set specific output points
tspan=[0,1.0e-5,1.0e-3,0.1,0.5,1,2,4,6,8,10,20,30,40,50,60];
%%2-> Set options
options = odeset('RelTol',1e-6,'AbsTol',1e-6); %set tolerances
%%3-> Call the ODE solver routine ode15s
[t,y]=ode15s(@odeFunc,tspan,y0,options);
%%4-> Output the result
disp(' i t y1 y2 y3 ');
fprintf(' -----------------------------------------\n');
for i=1:length(t)
fprintf('%5i %10.2e %10.2e %10.2e %10.2e\n',i,t(i),y(i,1),y(i,2),y(i,3));
end
figure(1)
semilogy(t,y(:,1),t,y(:,2),t,y(:,3),'LineWidth',2);
title(['Solution of the Robertson Problem for t=0:0.60s']);
xlabel('time (Sec)'); ylabel('solution y');
axis([tspan(1) tspan(end) 1.0e-12 10]);
figure(2)
semilogy(t,y(:,1),t,y(:,2),t,y(:,3),'LineWidth',2);
title(['Solution of the Robertson Problem for t=0:0.01s']);
xlabel('time (Sec)'); ylabel('solution y');
axis([tspan(1) 0.01 1.0e-12 10]);
% --------------------------------------------------------------------------
% Defining the ODE problem (or function) to be solved
function f = odeFunc(t,y)
f=zeros(3,1);
f(1)= -0.04*y(1)+1.0e4*y(2)*y(3);
f(2)= 0.04*y(1)-1.0e4*y(2)*y(3)-3.0e7*y(2)*y(2);
f(3)= 3.0e7*y(2)*y(2);
end
Computed result using Matlab stiff ODE solver – ode15s, for t=0:60, RelTol=1.e-6, AbsTol=1.e-6
i t y1 y2 y3
------------------------------------------------
1 0.00e+00 1.00e+00 0.00e+00 0.00e+00
2 1.00e-05 1.00e+00 4.00e-07 4.05e-11
3 1.00e-03 1.00e+00 2.78e-05 1.22e-05
4 1.00e-01 9.96e-01 3.58e-05 3.89e-03
5 5.00e-01 9.82e-01 3.33e-05 1.82e-02
6 1.00e+00 9.66e-01 3.07e-05 3.35e-02
7 2.00e+00 9.42e-01 2.70e-05 5.84e-02
8 4.00e+00 9.06e-01 2.24e-05 9.45e-02
9 6.00e+00 8.79e-01 1.96e-05 1.21e-01
10 8.00e+00 8.59e-01 1.77e-05 1.41e-01
11 1.00e+01 8.41e-01 1.62e-05 1.59e-01
12 2.00e+01 7.82e-01 1.23e-05 2.18e-01
13 3.00e+01 7.44e-01 1.04e-05 2.56e-01
14 4.00e+01 7.16e-01 9.19e-06 2.84e-01
15 5.00e+01 6.93e-01 8.34e-06 3.07e-01
16 6.00e+01 6.74e-01 7.71e-06 3.26e-01
