Matlab Consider the chemical reaction X rightarrowk1C Y The

Matlab

Consider the chemical reaction X rightarrow^k_1[C] Y. The molecule X is converted to Y. There is a catalyst [C] that controls the rate of the reaction, such that the rate of conversion of X to Y is k_1[C] molecules/second. The catalyst itself, however, is affected by the amount of X present. Namely, [C] = [X]/k_2 + [X]. Finally, suppose we have an external machine that is adding or removing X periodically, where the time-dependent rate is given by k_3 sin(pi t). Then the differential equation governing the concentration of [X] is: d[X]/dt = -k_1 [C][X] + k_3 sin (pi t). Suppose k_1 = 2, k_2 = 0.25, and k_3 = 1.5. Suppose also that we begin with [X] = 1 at t = 0. We wish to compute [X] as a function of t from t = 0 to t = 1. Integrate the differential equation to t = 1 using the forward Euler\'s method with a time step of h = 0.1. Save the resulting values ([X](t = 0), [X](t = 0.1), ..., [X](t = 1))^T as a column vector in A1.dat. Integrate the differential equation to t = 1 using the fourth order Runge-Kutta method with a time step of h = 0.1. Save the resulting values ([X](t = 0), [X](t = 0.1), ..., [X](t = 1))^T as a column vector in A2.dat.

Solution

(a)

function file

function [ z ] = fun( x,t )
z=-2*x^2/(0.25+x)+1.5*sin(pi*t);

end

script file

t0=0;%input(\'enter the initial value of t = \');
x0=1;%input(\'enter the initial value of y = \');

h=0.1;
d=1/h;
t=zeros(1,d+1);
x=zeros(1,d+1);
t(1,1)=t0;
x(1,1)=x0;
fprintf(\'the value of x at %f is %f\ \',t(1,1),x(1,1));
for n=2:d+1
t(1,n)=t(1,n-1)+h;
p=-h*fun(x(1,n-1),t(1,n-1));
x(1,n)=x(1,n-1)+0.5*(p-h*fun(x(1,n-1),t(1,n-1))+p);

fprintf(\'the value of x at %f is %f\ \',t(1,n),x(1,n));
end

answer

6.316966

(b)

function file

function [ z ] = fun( x,t )
z=-2*x^2/(0.25+x)+1.5*sin(pi*t);

end

script file

%%rk4 method%%
t0=0;%input(\'enter the initial value of t = \');
x0=1;%input(\'enter the initial value of y = \');
h=0.1;
d=1/h;
trk=zeros(1,d+1);
xrk=zeros(1,d+1);
trk(1,1)=t0;
xrk(1,1)=x0;
fprintf(\'the value of x at %f is %f\ \',trk(1,1),xrk(1,1));
for n=2:d+1
trk(1,n)=trk(1,n-1)+h;
end
for n=2:d+1
  
k1=-h*fun(x(1,n-1),t(1,n-1));
k2=-h*fun(x(1,n-1)+k1/2,t(1,n-1)+h/2);
k3=-h*fun(x(1,n-1)+k2/2,t(1,n-1)+h/2);
k4=-h*h*fun(x(1,n-1)+k3,t(1,n-1)+h);
xrk(1,n)=xrk(1,n-1)+(k1+2*k2+2*k3+k4)/6;
fprintf(\'the value of x at %f is %f\ \',trk(1,n),xrk(1,n));
end

answer

1
1.24
1.480055
1.727658
1.998407
2.317281
2.719768
3.253025
3.977256
4.967617
6.316966 (b) function file function [ z ] = fun( x,t ) z=-2*x^2/(0.25+x)+1.5*sin(pi*t); end script file %%rk4 method%% t0=0;%input(\'enter the initial value of t = \'); x0=1;%input(\'enter the initial value of y = \'); h=0.1; d=1/h; trk=zeros(1,d+1); xrk=zeros(1,d+1); trk(1,1)=t0; xrk(1,1)=x0; fprintf(\'the value of x at %f is %f\ \',trk(1,1),xrk(1,1)); for n=2:d+1 trk(1,n)=trk(1,n-1)+h; end for n=2:d+1    k1=-h*fun(x(1,n-1),t(1,n-1)); k2=-h*fun(x(1,n-1)+k1/2,t(1,n-1)+h/2); k3=-h*fun(x(1,n-1)+k2/2,t(1,n-1)+h/2); k4=-h*h*fun(x(1,n-1)+k3,t(1,n-1)+h); xrk(1,n)=xrk(1,n-1)+(k1+2*k2+2*k3+k4)/6; fprintf(\'the value of x at %f is %f\ \',trk(1,n),xrk(1,n)); end answer 1 1.130076 1.261595 1.40068 1.558193 1.750375 1.999462 2.33435 2.791462 3.415978 4.263637