Question

Problem 4 (Stiff ODE) We consider a chemical reaction system due to H. Robertson that has been extensively tested as a test problem for stiff solvers yí = -ayi + By2y3, y = ayı – By2y3 – xy2, y = vyž where a = 0.04, B = 1.E +4 and y = 3. E + 7 are slow, fast, and very fast reaction rates. The starting point is y(0) = (1,0,0). (a) It is known that this system reaches a steady state, i.e., where y' = 0. Show that [i-1 yi(t) = 1 for 0 <t<b for any b > 0. Use it to find the steady state. (b) Solve the system using ODE45 solver to obtain y(3). How far is y(3) from the steady state and how many time step is used to achieve it ? (c) The steady state is reached very slowly for this problem. Use ODE15s to integrate the problem for b = 1.E + 6 and plot the solution in t. How far is y(b) from the steady state and how many time step size is used ?

Answer 1

title('t vs. yl(t) plot for ode45') xlabel('t') ylabel('y(t)') subplot(3,1,2) plot (t1, soli(:,2)) title('t vs. y2(t) plot for ode45'). xlabel('t') ylabel('y2(t)') subplot(3,1,3) plot(t1, soli(:,3)) title('t vs. y2(t) plot for ode45') xlabel('t') ylabel('y3(t)') fprintf('At t=1.E+6, yl(1. E+6)=%f, y2(1. E+6)=%f and y3 (1.E+6)= %f\n',... soll (end, 1), soll(end, 2), soll (end, 3)) fprintf('Using ode 15s it required %d time steps to achive it. \n', length(t1)) %Function for evaluating the ODE function duldt = odefcn(t, y, alpha, beta, gamma) eql = -alpha*y (1) +beta*y (2)*y(3); eq2 = alpha*y(1)-beta*y(2) *y(3)-gamma* (y(2)).^2; eq3 = gamma* (y (2)).^2; %Evaluate the ODE for our present problem duldt = [eql;eq2;eq3]; end %%%%%%%%%%%%%%% End of Code % % % % % % % % % % % % % % % % At t=3, y1 (3)=0.921887, y2 (3)=0.000024 and y3(3)=0.078089 Using ode45 it required 8205 time steps to achive it. At t=1.E+6, yl(1.E+6)=0.002032, y2 (1.E+6)=0.000000 and y3(1. E +6)=0.997968 Using ode15s it required 146 time steps to achive it.

title('t vs. yl(t) plot for ode45') xlabel('t') ylabel('y(t)') subplot(3,1,2) plot (t1, soli(:,2)) title('t vs. y2(t) plot for ode45'). xlabel('t') ylabel('y2(t)') subplot(3,1,3) plot(t1, soli(:,3)) title('t vs. y2(t) plot for ode45') xlabel('t') ylabel('y3(t)') fprintf('At t=1.E+6, yl(1. E+6)=%f, y2(1. E+6)=%f and y3 (1.E+6)= %f\n',... soll (end, 1), soll(end, 2), soll (end, 3)) fprintf('Using ode 15s it required %d time steps to achive it. \n', length(t1)) %Function for evaluating the ODE function duldt = odefcn(t, y, alpha, beta, gamma) eql = -alpha*y (1) +beta*y (2)*y(3); eq2 = alpha*y(1)-beta*y(2) *y(3)-gamma* (y(2)).^2; eq3 = gamma* (y (2)).^2; %Evaluate the ODE for our present problem duldt = [eql;eq2;eq3]; end %%%%%%%%%%%%%%% End of Code % % % % % % % % % % % % % % % % At t=3, y1 (3)=0.921887, y2 (3)=0.000024 and y3(3)=0.078089 Using ode45 it required 8205 time steps to achive it. At t=1.E+6, yl(1.E+6)=0.002032, y2 (1.E+6)=0.000000 and y3(1. E +6)=0.997968 Using ode15s it required 146 time steps to achive it.

tvs. y(t) plot for ode45 0.98 0.96 0.94 0 0. 5 1 1.5 tvs. y2(t) plot for ode45 211 y2(t) 0.5 15 2.5 tvs. y2(t) plot for ode45 y3(t) 0. 51 1.5 2 2.5 tvs. y(t) plot for ode45 0.5 0 1 2 3 4 5 6 7 8 9 10 X105 tvs. y2(t) plot for ode45 MEK O 1 2 3 4 5 6 7 8 9 10 *105 tvs. y2(t) plot for ode45 y3(t) 0 1 2 3 4 5 6 7 8 9 10 x 10°

%%Matlab code for solving ode
clear all
close all
%Answering question
%Initial conditions for ode
    u0=[1;0;0];
    alpha=0.04; beta=1*10^4; gamma=3*10^7;
%Solution for equation 1. using ode45
        %minimum and maximum time span
        tspan=[0 3];
        %Solution of ODEs using ode45 matlab function
        [t1,sol1]= ode45(@(t,y) odefcn(t,y,alpha,beta,gamma), tspan, u0);

        %Plotting the solution
        figure(1)
        subplot(3,1,1)
        plot(t1,sol1(:,1))
        title('t vs. y1(t) plot for ode45')
        xlabel('t')
        ylabel('y(t)')
      
        subplot(3,1,2)
        plot(t1,sol1(:,2))
        title('t vs. y2(t) plot for ode45')
        xlabel('t')
        ylabel('y2(t)')
     
        subplot(3,1,3)
        plot(t1,sol1(:,3))
        title('t vs. y2(t) plot for ode45')
        xlabel('t')
        ylabel('y3(t)')
      
        fprintf('At t=3, y1(3)=%f, y2(3)=%f and y3(3)=%f\n',...
            sol1(end,1),sol1(end,2),sol1(end,3))
        fprintf('Using ode45 it required %d time steps to achive it.\n',length(t1))
%Initial conditions for ode
    u0=[1;0;0];
    alpha=0.04; beta=1*10^4; gamma=3*10^7;
%Solution for equation 1. using ode45
        %minimum and maximum time span
        tspan=[0 10^6];
        %Solution of ODEs using ode45 matlab function
        [t1,sol1]= ode15s(@(t,y) odefcn(t,y,alpha,beta,gamma), tspan, u0);

        %Plotting the solution
        figure(2)
        subplot(3,1,1)
        plot(t1,sol1(:,1))
        title('t vs. y1(t) plot for ode45')
        xlabel('t')
        ylabel('y(t)')
      
        subplot(3,1,2)
        plot(t1,sol1(:,2))
        title('t vs. y2(t) plot for ode45')
        xlabel('t')
        ylabel('y2(t)')
     
        subplot(3,1,3)
        plot(t1,sol1(:,3))
        title('t vs. y2(t) plot for ode45')
        xlabel('t')
        ylabel('y3(t)')
      
        fprintf('At t=1.E+6, y1(1.E+6)=%f, y2(1.E+6)=%f and y3(1.E+6)=%f\n',...
            sol1(end,1),sol1(end,2),sol1(end,3))
        fprintf('Using ode15s it required %d time steps to achive it.\n',length(t1))      
      
      

%Function for evaluating the ODE
function du1dt = odefcn(t,y,alpha,beta,gamma)

    eq1 = -alpha*y(1)+beta*y(2)*y(3);
    eq2 = alpha*y(1)-beta*y(2)*y(3)-gamma*(y(2)).^2;
    eq3 = gamma*(y(2)).^2;
  
  
    %Evaluate the ODE for our present problem
    du1dt = [eq1;eq2;eq3];
end

%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%