Chemical Reaction

 

 

 

 

 

 

Followings are the code that I wrote in Octave to creates all the plots shown in this page. You may copy these code and play with these codes. Change variables and try yourself until you get your own intuitive understanding.

 

< Code 1 >

 

clear;

 

t = 0:0.01:4;

global k1 = 1.0;

global k2 = 1.0;

global k3 = 4.0;

  

function dy_dt = f(y, t)

  global k1;

  global k2;

  global k3;

  dy_dt = zeros (3,1);

  dy_dt(1) = -k1 .* y(1) + k2 .* y(2);

  dy_dt(2) = k1 .* y(1) - (k2 + k3) .* y(2);

  dy_dt(3) = k3 * y(2);

endfunction

 

y3_init = 0.0;

y2_init = 0.0;

y1_init = 1.0;

 

y0 = [y1_init y2_init y3_init];

y = lsode ("f", y0, t);

 

hFig = figure(1);

plot(t,y(:,1),'r-','LineWidth',2,t,y(:,2),'g-','LineWidth',2,t,y(:,3),'b-','LineWidth',2);

xlabel('time');

xlim([0 t(end)]);

ylim([0 1]);

legend('A(t)','B(t)','C(t)');

grid on;

 

tStr = sprintf("A(0)=%0.2f,B(0)=%0.2f,C(0)=%0.2f \nk1=%0.2f,k2=%0.2f,k3=%0.2f",

              y1_init,y2_init,y3_init,k1,k2,k3);

title(tStr);

 

set(hFig,'Position',[300 300 780 300]);