ode_drv.m

ode routines in matlab

% 
% 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

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


rhs_fctn   = 'ChemReac';

h = 0.2;

% Exact solution
mx    = ceil((x_end-x0)/ h);
x     = (x0:h:x0+mx*h);


[x, y_exp] = ExpEuler( x0, x_end, y0, h, rhs_fctn );

% Plot the computed solution
subplot(2,2,1)
plot(x,y_exp(1,:),'g',x,y_exp(2,:),'r',x,y_exp(3,:),'b');
legend('y_1', 'y_2', 'y_3') 
xlabel(' x ')
title([' Explicit Euler, h = ', num2str(h) ])


[x, y_imp] = ImpEuler( x0, x_end, y0, h, rhs_fctn );
% Plot the computed solution
subplot(2,2,3)
plot(x,y_imp(1,:),'g',x,y_imp(2,:),'r',x,y_imp(3,:),'b');
legend('y_1', 'y_2', 'y_3') 
xlabel(' x ')
title([' Implicit Euler, h = ', num2str(h) ])

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

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

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

% Plot the error in the computed solution
subplot(2,2,2)
plot(x,abs(y_exp(1,:)-y_ex(1,:)),'g'); hold on
plot(x,abs(y_exp(2,:)-y_ex(2,:)),'r'); hold on
plot(x,abs(y_exp(3,:)-y_ex(3,:)),'b'); hold off
legend('y_1', 'y_2', 'y_3') 
xlabel(' x ')
title([' Error, Explicit Euler, h = ', num2str(h) ])

% Plot the error in the computed solution
subplot(2,2,4)
plot(x,abs(y_imp(1,:)-y_ex(1,:)),'g'); hold on
plot(x,abs(y_imp(2,:)-y_ex(2,:)),'r'); hold on
plot(x,abs(y_imp(3,:)-y_ex(3,:)),'b'); hold off
legend('y_1', 'y_2', 'y_3') 
xlabel(' x ')
title([' Error, Implicit Euler, h = ', num2str(h) ])


return
disp(' Hit key to continue .....')
pause



%  Solution of the 
%                         
%    y_1'(x) =  y_2(x)
%               
%    y_2'(x) =  4*y_1(x) + 3*y_2(x) - 8*exp(x) * cos(2x)
%                
%  with initial conditions
%
%    y_1(0) =  2 + 10/13             
%    y_2(0) = -3 + 14/13
%
%  This system is equivalent to the second order ODE
%
%    z''(x) - 3 z'(x) - 4 z(x) = - 8 exp(x) * cos(2x)
%
%  The exact solution of the system is given by
%    y_1(x) =  exp(-x) + exp(4*x) 
%               + 10/13 * exp(x) * cos(2x) + 2/13 * exp(x) * sin(2x)
%               
%    y_2(x) = -exp(-x) + 4 * exp(4*x) 
%               + 14/13 * exp(x) * cos(2x) - 18/13 * exp(x) * sin(2x)
%

clf;

x0    = 0;
x_end = 1;
y0    = [ 2 + 10/13; 3 + 14/13];

rhs_fctn   = 'Ex1';

h = 0.02;

% Exact solution
mx    = ceil(x_end-x0)/ h;
x     = (x0:h:x0+mx*h);
y_ex  = [ exp(-x) + exp(4*x) + (10/13)*exp(x).*cos(2*x) ...
              + (2/13)*exp(x).*sin(2*x);
         -exp(-x) + 4*exp(4*x) + (14/13)*exp(x).*cos(2*x) ...
              - (18/13).*exp(x).*sin(2*x) ];


[x, y_exp] = ExpEuler( x0, x_end, y0, h, rhs_fctn );

% Plot the computed solution
subplot(2,2,1)
plot(x,y_exp(1,:),'g'); hold on
plot(x,y_exp(2,:),'r'); hold off
legend('y_1', 'y_2') 
xlabel(' x ')
title([' Explicit Euler, h = ', num2str(h) ])

% Plot the error in the computed solution
subplot(2,2,2)
plot(x,abs(y_exp(1,:)-y_ex(1,:)),'g'); hold on
plot(x,abs(y_exp(2,:)-y_ex(2,:)),'r'); hold off
legend('y_1', 'y_2') 
xlabel(' x ')
title([' Error, Explicit Euler, h = ', num2str(h) ])

[x, y_imp] = ImpEuler( x0, x_end, y0, h, rhs_fctn );
% Plot the computed solution
subplot(2,2,3)
plot(x,y_imp(1,:),'g'); hold on
plot(x,y_imp(2,:),'r'); hold off
legend('y_1', 'y_2') 
xlabel(' x ')
title([' Implicit Euler, h = ', num2str(h) ])

% Plot the error in the computed solution
subplot(2,2,4)
plot(x,abs(y_imp(1,:)-y_ex(1,:)),'g'); hold on
plot(x,abs(y_imp(2,:)-y_ex(2,:)),'r'); hold off
legend('y_1', 'y_2') 
xlabel(' x ')
title([' Error, Implicit Euler, h = ', num2str(h) ])


disp(' Hit key to continue .....')
pause


%  Predator Prey model 1
%
%
%    y_1'(x) = alpha*y_1(x) - beta*y_1(x)*y_2(x)
%    y_2'(x) = gamma*y_2(x)  + delta*y_1(x)*y_2(x)
%
%
%

clf;

x0    = 0;
x_end = 20;
y0    = [ 80; 30];

rhs_fctn   = 'PredPrey';

h = 0.05;

[x, y_exp] = ExpEuler( x0, x_end, y0, h, rhs_fctn );

% Plot the computed solution
subplot(2,2,1)
plot(x,y_exp(1,:),'g'); hold on
plot(x,y_exp(2,:),'r'); hold off
legend('Population 1', 'Population 2') 
xlabel(' Time ')
title([' Explicit Euler, h = ', num2str(h) ])


[x, y_imp] = ImpEuler( x0, x_end, y0, h, rhs_fctn );

% Plot the computed solution
subplot(2,2,2)
plot(x,y_imp(1,:),'g'); hold on
plot(x,y_imp(2,:),'r'); hold off
legend('Population 1', 'Population 2') 
xlabel(' Time ')
title([' Implicit Euler, h = ',  num2str(h) ])


% Plot the difference in computed solutions
subplot(2,2,3)
plot(x,abs(y_exp(1,:)-y_imp(1,:)),'g'); hold on
plot(x,abs(y_exp(2,:)-y_imp(2,:)),'r'); hold off
legend('Population 1', 'Population 2') 
xlabel(' Time ')
title([' Difference in Computed Solutions, h = ',  num2str(h) ])



disp(' Hit key to continue .....')
pause



%  Predator Prey model 2
%
%                                         y_1(x) * y_2(x)
%    y_1'(x) = 1.2 * y_1(x) - y_1^2(x) -  ---------------
%                                          y_1(x) + 0.2
%
%               1.5 * y_1(x) * y_2(x)
%    y_2'(x) =  ---------------------  -  y_2(x)
%                   y_1(x) + 0.2
%
%
%

clf;

x0    = 0;
x_end = 30;
y0    = [ 1; 0.75];

rhs_fctn   = 'PredPrey2';

h = 0.001;

[x, y_exp] = ExpEuler( x0, x_end, y0, h, rhs_fctn );

% Plot the computed solution
subplot(2,2,1)
plot(x,y_exp(1,:),'g'); hold on
plot(x,y_exp(2,:),'r'); hold off
legend('Population 1', 'Population 2') 
xlabel(' Time ')
title([' Explicit Euler, h = ', num2str(h) ])


[x, y_imp] = ImpEuler( x0, x_end, y0, h, rhs_fctn );

% Plot the computed solution
subplot(2,2,2)
plot(x,y_imp(1,:),'g'); hold on
plot(x,y_imp(2,:),'r'); hold off
legend('Population 1', 'Population 2') 
xlabel(' Time ')
title([' Implicit Euler, h = ',  num2str(h) ])


% Plot the difference in computed solutions
subplot(2,2,3)
plot(x,abs(y_exp(1,:)-y_imp(1,:)),'g'); hold on
plot(x,abs(y_exp(2,:)-y_imp(2,:)),'r'); hold off
legend('Population 1', 'Population 2') 
xlabel(' Time ')
title([' Difference in Computed Solutions, h = ',  num2str(h) ])