Differential equations - Chemical Reaction Rates

Hi everyone!

I have been trying to figure out how to solve and graph this set of equations in MATLAB. I have tried different variations of ODE formatting and cannot get it right.

The goal is to be able to change the initial concetrations and see how the concentration vs. time graph changes.

k1, k2, k3 are constant so I have been setting those to equal 1. The time interval is also constant [1,10] seconds. Initially, [E] is higher and [ES] is zero. As the reaction progresses, [ES] increases and [E] decreases.

If you can help, that would be much appreciated, thank you!!

E + S < -- > ES --> E + P;
dsdt = -k1[E][S]+ k2[ES];
dpdt = +k3[ES];
dedt = -k1[E][S]+k2[ES]+k3[ES];
dcdt = +k1[E][S]-k2[ES]-k3[ES];

[James Tursa]

Can you post the code you have written so far? Then we can make suggestions for improvements/corrections.

[Ariel Patterson-Gonzales]

I would also like to view dsdt, dedt, dpdt, AND dcdt in the graph and dont know how to add that.. below is my code

function [ f ] = dXdT( time,x )

%dedt+dcdt = 0 (all rate laws cancel out)
%  SO: (de+dc)/dt = d(e+c)/dt = 0      Eo=e+c and c=Eo-e

s = x(1);
p = x(2);
e = x(3);

   k1 = 5; %1/M*sec
   k2= 0.1;  %1/s
   k3 = 3; % 1/s
   E0 = 0;

      dsdt = -k1*e*s + k2*(E0-e);
      dpdt = k3*(E0-e);
      dedt = -k1*e*s + k2*(E0-e) + k3*(E0-e);

   f = [s; p; e];
  % to Find rates, run in the command line "dXdT([0,10],[ 1 0 1 ])" the numbers are

%GRAPH REPRESENTATION:
% where tspan = [t0 tf]
%for this example: [time,x] = ode15s(@dXdT,[0,10],[1 0 1]);
% type "time" and enter for the time and "x" and enter to see output values
% TO SEE A GRAPH: "plot(time,x)"

Answer 1

[Star Strider]

Try this:

f = [dsdt; dpdt; dedt];

Note that you didn’t code ‘dcdt’. You need to include it in your calculations, and in the ‘f’ vector.

[Ariel Patterson-Gonzales]

Thank you, I had f = [dsdt; dpdt; dedt]; and it resulted in the same graphs. How do I add dcdt into it? I am new with MATLAB and have tried to get this for a month now. eek!

[Star Strider]

My pleasure.

Add ‘dcdt’ by coding this equation:

dcdt = +k1[E][S]-k2[ES]-k3[ES];

as:

dcdt =  k1*e.*s - k2*.(E0-e) - k3*(E0-e);

with ‘f’ changed to:

f = [dsdt; dpdt; dedt; dcdt];

since it seems ‘dcdt’ and ‘dedt’ are the same except for the signs of the constants.

To plot it, use:

figure(1)
plot(time, x)
grid
legend('ds/dt', 'dp/dt', 'de/dt', 'dc/dt')

NOTE — I did not run this, but it should work if your previous code worked.

[Ariel Patterson-Gonzales]

Thank you so much! Do you know how to define time in the script?

[Star Strider]

My pleasure!

The ODE integration function you use (here ode15s) calculates and returns a time vector corresponding to the times it integrated your ODE as the first output. The integrated values of your ODE are the second output.

Quoting your original Question:

[time,x] = ode15s(@dXdT,[0,10],[1 0 1]);

See the documentation for ode15s for a complete description and explanation.

[Ariel Patterson-Gonzales]

So I have to define time in the command window?

[Star Strider]

No. The integration time vector (here ‘time’) is created by ode15s and appears in your workspace as a result of the integration. You don’t have to do anything. To plot the integrated equations as a function of time, you would plot them as:

[time,x] = ode15s(@dXdT,[0,10],[1 0 1]);
plot(time, x)

That’s all there is to it!

[Ariel Patterson-Gonzales]

Thank you so much for the clarification on time! You have been so helpful, sorry if my questions are redundant. After running the script with your additions, dc/dt is not showing on the graph..

[Star Strider]

It is for me:

function [ f ] = dXdT( time,x )

%dedt+dcdt = 0 (all rate laws cancel out)
%  SO: (de+dc)/dt = d(e+c)/dt = 0      Eo=e+c and c=Eo-e

s = x(1);
p = x(2);
e = x(3);

     k1 = 5; %1/M*sec
     k2= 0.1;  %1/s
     k3 = 3; % 1/s
     E0 = 0;

        dsdt = -k1*e*s + k2*(E0-e);
        dpdt = k3*(E0-e);
        dedt = -k1*e*s + k2*(E0-e) + k3*(E0-e);
        dcdt =  k1*e.*s - k2.*(E0-e) - k3.*(E0-e);

f = [dsdt; dpdt; dedt; dcdt];
end

[time,x] = ode15s(@dXdT,[0,10],[1 0 1 0]);

figure(1)
plot(time, x)
grid
legend('ds/dt', 'dp/dt', 'de/dt', 'dc/dt')

%

[Ariel Patterson-Gonzales]

I still only get 3 lines and the errors:

"Error using odearguments (line 92) DXDT returns a vector of length 3, but the length of initial conditions vector is 4. The vector returned by DXDT and the initial conditions vector must have the same number of elements.

Error in ode15s (line 148) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ..."

and

"Warning: Ignoring extra legend entries. > In /Applications/MATLAB_R2014b.app/toolbox/matlab/scribe/private/legendHGUsingMATLABClasses.p>set_children_and_strings at 649 In /Applications/MATLAB_R2014b.app/toolbox/matlab/scribe/private/legendHGUsingMATLABClasses.p>make_legend at 312 In /Applications/MATLAB_R2014b.app/toolbox/matlab/scribe/private/legendHGUsingMATLABClasses.p>legendHGUsingMATLABClasses at 241 In legend at 118 "

[Star Strider]

Did you use my code in my previous Comment? It ran for me without error. I guessed at the initial condition for ‘dc/dt’ of zero, so change that if necessary, but my posted code works.

[Ariel Patterson-Gonzales]

I copy pasted that exact code.. not sure why it isnt working.

[Ariel Patterson-Gonzales]

and yeah, dc/dt is initially zero

[Star Strider]

I cannot understand your using my exact code — all of it including the ‘f’ vector that I changed to add ‘dcdt’ — and getting any errors.

[Ariel Patterson-Gonzales]

I am reinstalling MATLAB.. I do not know how it can work for you and not for me.. I'll keep in touch

[Star Strider]

You may not have to reinstall MATLAB.

First, use:

which dXdT.m -all

to check to be sure you only have one copy of ‘dXdT.m’ and it is the one that matches my code. If you have more than one copy, name them ‘dXdT_1.m’ or something similar so you can still have them for reference but not run them with your current ode15s call.

[Ariel Patterson-Gonzales]

I reinstalled, and deleted all files and it works!!! Thank you so much and for your patience!

[Star Strider]

My pleasure!

If my Answer solved your problem, please Accept it.

[John BG]

Mr Start Strider straight to the point, and improving marketing skills ;)