ode23_drv.m

来自「ode routines in matlab」· M 代码 · 共 61 行

% 
% ode_drv.m: simple driver for ode solvers
%
%


%  Solution of the chemical reaction system
%                         
%    y_1'(x) = -k1*y_1(x) 
%    y_2'(x) =  k1*y_1(x) - k2*y_2(x)
%    y_3'(x) =              k2*y_2(x)
%
%

clf;

% declare reaction rates as global variables

global k1 k2


y0    = [ 1; 0; 0];
k1    = 1;
k2    = 10;


rhs_fctn   = 'ChemReac23';

reltol = 1.e-3;
abstol = 1.e-3;
opts   = odeset('reltol',reltol,'abstol', abstol);

[t, y] = ode23tx( rhs_fctn, [0,10], y0, opts );

% Plot the computed solution
figure(1)
plot(t,y(:,1),'g+',t,y(:,2),'r+',t,y(:,3),'b+');
legend('y_1', 'y_2', 'y_3') 
xlabel(' t ')
title(['ode23tx (abstol = ', num2str(abstol), 'reltol = ', num2str(reltol),')'])

% Compute exact solution
A = [ -k1   0   0;
       k1  -k2  0;
       0    k2  0 ];
[V,D] = eig(A);

y_ex = zeros(length(t),3);

for j = 1:length(t)
    y_ex(j,:)  = (V*(exp(diag(D)*t(j)).*(V\y0)))';
end

% Plot the error in the computed solution
figure(2)
semilogy(t,abs(y(:,1)-y_ex(:,1)),'g'); hold on
semilogy(t,abs(y(:,2)-y_ex(:,2)),'r'); hold on
semilogy(t,abs(y(:,3)-y_ex(:,3)),'b'); hold off
legend('y_1', 'y_2', 'y_3') 
xlabel(' t ')
title(['ode23tx (abstol = ', num2str(abstol), 'reltol = ', num2str(reltol),')'])