An Introduction to the Use of Numerical Methods in Chemical Kinetics

Michelle M. Francl
Bryn Mawr College
101 N. Merion Ave
Bryn Mawr, PA 19010 USA

mfrancl@brynmawr.edu

Copyright M.M. Francl, 2000.

Once a reaction mechanism reaches more than a few steps, it becomes difficult to integrate the differential rate equations analytically.  Most studies of complex mechanisms rely either on simplyfying the reaction mechanism until it looks like one of the simple motifs we have already explored, or by numerically solving the differential equations.  In the lab, we are using the first approach.  The complex mechanism given for the iodination of the ketones is simplified by assuming that we need only consider the kinetics of the rate determining step.  When one has a large collection of reactions in which no steady state intermediates or rate determining steps can be identified, one is driven to using numerical methods.

Differential equations of the sort that commonly arise in chemical kinetics can be solved numerically using the Runge-Kutta method (see Programming in Mathematica, R. Maeder, 2nd edition, pp. 172-173 for a brief description; additional information can be found in texts on differential equations or numerical methods).  I have used Maeder's package from MathSource at Wolfram's WWW site to do the numerical solution (http://www.mathsource.com).

I've included some of the more complex commands (especially for plotting) that are asked for on the accompanying worksheet in this notebook.

Sequential reactions

Consider the following two step mechanism:
A -> B
B -> C

We can easily write down the differential rate equations:

d[A]/dt = -k1[A]
d[B]/dt = k1[A] - k2[B]
d[C]/dt = k2[B]

These can be solved to yield the following expressions: (in Mathematica®  form)

A[t_]:=Ao Exp[-k1 t]
B[t_]:=(k1/(k2-k1)) (Exp[-k1 t] - Exp[-k2 t]) Ao
C[t_]:=Ao-B[t]-A[t]

Using the initial conditions [A]o=3.7, k1=1.2 and k2=1.5 the following graph was obtained:

The analytical solution

A[t_]:=Ao Exp[-k1 t]
B[t_]:=(k1/(k2-k1)) (Exp[-k1 t] - Exp[-k2 t]) Ao
Pr[t_]:=Ao-B[t]-A[t]

Setting up the rate constants:

Ao=3.7;
k1=1.2;
k2=1.5;

Plotting the results:

analsoln=Plot[{A[t],B[t],Pr[t]},{t,0,10}]

Plotting results from a simple minded numerical integrator

❏✀✀First, plot the results from the analytical solutions given above for the initial conditions [A]o=3.7M, k1=1.2/sec and k2=1.5/sec.
❏✀✀Now, try my simple minded solution of the differential equations.  To solve 
dA/dt = k1A
use a recursive solution of the form
A(ti+1) = A(ti) + k1A(ti)Δt
Try a time step of 0.5 seconds for a total reaction time of 5 seconds.  
	✎✀✀Write the formulas for the recursion relation for B and C.
	✎✀✀Are the plots very different?  In what ways are they different?  Does changing the time step affect anything?

The following command should be done only once!!

<< Graphics`MultipleListPlot`

Here is the integration  routine.  To watch the  reaction for a longer period, change Total Reaction Time.  To change the time step (basically the quality control of the integration- the smaller the time step, the smoother the resulting plot and the longer the integration will take)  change DeltaTime.  The rest is automatic!

TotalReactionTime = 5;
DeltaTime = 0.1;
TotalSteps = TotalReactionTime/DeltaTime;
			a = Table[0, {TotalSteps + 1. }];
       b = Table[0, {TotalSteps + 1. }];
       c = Table[0, {TotalSteps + 1. }];
a[[1]] = 3.7;b[[1]] = 0;c[[1]] = 0;
k1 = 1.2;k2 = 1.5;
Do[{a[[i + 1]] = a[[i]] - k1 a[[i]] DeltaTime,
          b[[i + 1]] =
            b[[i]] + (k1 a[[i]]  - k2 b[[i]]) DeltaTime,
          c[[i + 1]] = c[[i]] + k2 b[[i]] DeltaTime}, {i,
          TotalSteps}]\(time = Table[DeltaTime\ i, {i, 0, TotalSteps}]; 

    ConcA = Transpose[{time, a}]; 
ConcB = Transpose[{time, b}]; 

    ConcC = Transpose[{time, c}]; 

    numplot =
      MultipleListPlot[ConcA, \ ConcB, ConcC, 
\ PlotJoined\  -> \ True]\)Show[analsoln, numplot]

Solving the differential rate equations using the Runge-Kutta method

❏✀✀Now use the more sophisticated Runge-Kutta solution.  The commands to do the solution and the plotting are in the section titled "Solving the differential rate equations using the Runge-Kutta method".  I’ve installed a function in the notebook to do the actual Runge-Kutta solution  for you (if you’d like to look at it, it’s under the heading "Runge-Kutta package").  In order to use the package, you must have had Mathematica do the initialization for you – see note at the beginning of the worksheet!!

Use the same initial conditions as you used for the analytical solution, set the time interval (incr) to be 0.1 sec and the total time (total)to be 5 seconds.
	✎✀✀Are the plots very different?  In what ways are they different?
	✎✀✀What happens if you choose  a different  (shorter or longer) time step?  Do the graphs change noticeably?
	✎✀✀Are the numerical methods adequate for this problem?  Does the Runge-Kutta or the simple-minded solver work better?

solve the differential equations and put the concetration vs. time data into a list of points:
{[A],[B],[C]}, one for each time

total = 5 ;
incr=.1 ;
Ao=3.7;Bo=0;Co=0;
k1=1.2;
k2=1.5;

points = Block[{a, b, c},
RK[{-k1 a,k1 a -k2 b,k2 b},
	{a, b, c}, {Ao, Bo,Co}, {total, incr}]];

prepare the data for plotting (this time I've used dotted lines....)

points=Transpose[points];
reactant=points[[1]];
intermediate=points[[2]];
pr=points[[3]];
time=Table[i,{i,0,total,incr}];
rline=Graphics[{Dashing[{.02,.02}],Line[Transpose[{time,reactant}]]}];
iline=Graphics[{Dashing[{.02,.02}],Line[Transpose[{time,intermediate}]]}];
pline=Graphics[{Dashing[{.02,.02}],Line[Transpose[{time,pr}]]}];

Plotting

rkplot = Show[{rline, iline, pline}, Axes -> Automatic]

Now compare the three plots!

Show[numplot, rkplot]

Show[analsoln, rkplot]

Some more complex reaction mechanisms

Now consider a more complex reaction, an 'equilibrium' followed by a single step.

A -> B	kf
B -> A	kr
B -> C	k3

❏  Start by assuming that a rapid pre-equilibrium is established (that is kf and kr>> k3).  
	✎✀✀Derive an expression for the rate of production of product C.  
❏  Using the expression you derived above, find the integrated rate for [C].  Plot the resulting function assuming the following starting conditions:  [A]o=3.7M, [B]o=0, [C]o=0,kf=120/sec kr=150/sec; k3=0.8.  
❏  Now find the numerical solution for [C] which does not rely on the approximation you made above.  (As a bonus you’ll get the values for [A] and [B] as well.)  The Mathematica® commands you need to do that are in the notebook under "Some more complex reaction mechanisms", "Solving the differential rate equations".  Use the same initial conditions and rate constants given in the previous step.  I’d suggest a total time of 3 seconds and an increment of .01 seconds (if the increment gets too large, the numerical results tend to blow up.  Plot the results.  
A hint:  if the numerical solution blows up (you will get a message about overflow), try using a shorter time step combined with a shorter reaction period
	✎✀✀Does the plot look very different from the approximate solution you generated in the last step? (Compare both graphs on the same plot using Show.)Does it look as you expect?
❏   Repeat the Runge-Kutta solution for conditions in which the last step is very fast compared to the equilibrium (k3>>kf,kr).  Use the same initial values of [A], [B] and [C].  
A hint:  if the numerical solution blow up (you will get a message about overflow), try using a shorter time step combined with a shorter reaction period; here you might also try varying the rate constants.  I’d suggest trying k3 150/sec and kf,kr around 0.8 (don’t make them both the same).
	✎✀✀Does it look as you expect?

Approximate analytical soln

int=Integrate[((k3 kf)/(kr+k3)) Ao Exp[((kr kf/(kr+k3))-kf) t],{t,0,t}]

Pr[t_]:= Ao - Ao*E^((-kf + (kf*kr)/(kr + k3))*t)
A[t_]:=Ao E^((kr kf/(kr+k3)-kf) t)
B[t_]:=kf A[t] /(kr+k3)

Ao=3.7;kf=.9;kr=.8;k3=120;
apprx=Plot[{A[t],B[t],Pr[t]},{t,0,10}]

Solving the differential rate equations numerically

set up the constants

I've chosen constants here such that the reverse reaction is very rapid, and that both forward and
reverse reactions are rapid compared to the third step.  This is sometimes called a rapid pre-equilibrium type system.  It's important to keep incr small here (0.01 seems to work well) otherwise the numerical solution blows up.  I've run the reaction for 3.2 seconds.

Ao=3.7;Bo=0;Co=0;
kf=.9;kr=.8;k3=120;
total=6;incr=.01;

apply Runge-Kutta to the solution of the differential equations

rkpoints = Block[{a, b, c},
RK[{-kf a+kr b,kf a - kr b-k3 b, k3 b},
	{a, b,c }, {Ao, Bo, Co}, {total, incr}]];

preparing the data to plot [A], [B] and [C] as a function of time

points=Transpose[rkpoints];
Aconc=points[[1]];
Bconc=points[[2]];
Cconc=points[[3]];
time=Table[i,{i,0,total,incr}];
aline=Graphics[{Thickness[.008],GrayLevel[.5],Line[Transpose[{time,Aconc}]]}];
bline=Graphics[{Thickness[.008],GrayLevel[.5],Line[Transpose[{time,Bconc}]]}];
cline=Graphics[{Thickness[.008],GrayLevel[.5],Line[Transpose[{time,Cconc}]]}];

plot

rkplot=Show[{aline,bline,cline},Axes->Automatic]

Show[{apprx,rkplot}]

The Lotka-Volterra mechanism

The Lotka-Volterra mechanism describes the kinetics of an autocatalytic reaction.  Depending on the conditions under which the reaction is run, the results can be dramatically different.  In some cases, the concentrations of the intermediates are constant, suggesting steady state kinetics.  Under other conditions the autocatalysis is an oscillating reaction, one in which the concentrations of reactants, intermediates and products oscillate in both time and space as the reaction progresses.  For more information see Physical Chemistry, P.W. Atkins, 5th edition, pp. 918-919.

The proposed mechanism is:

A + X -> 2X  k1
X + Y -> 2Y   k2
Y -> B            k3

The associated differential equations for [X] and [Y] are:

d[X]/dt = k1[A][X] - k2[X][Y]
d[Y]/dt = k2[X][Y] - k3[Y]

Typical reaction conditions are such that [A] is constant.  B is typically removed from the
 system, but since the differential rates don't depend on B, this is immaterial to our work.

❏  Apply the steady state approximation for X and Y to the equations above and find expressions for the steady state values of [X] and [Y].
❏  Using the commands given below, integrate the Lotka-Volterra rate equations for X and Y for the following conditions:  [A]0=3.0,[ X]o=0.8, [Y]o=2.0, k1=1.0, k2=1.5 and k3=1.2.  I’d suggest a time step of 0.1 and a total time of 5.0 seconds.
	✎✀✀Does the plot make sense?  Is it what you expect based on your answer to the previous question?
❏  Integrate the Lotka-Volterra rate equations for X and Y for the following conditions:  [A]0=3.0,[ X]o=0.9, [Y]o=2.1, k1=1.0, k2=1.5 and k3=1.2.
	✎✀What do you observe?  Is this unusual?
❏  Try varying the rate constants and intial concentrations (don’t change everything at once! ).
	✎✀✀What happens?
❏  Try checking the behavior of the reaction in concentration space!  Plot [X] vs. [Y] for various conditions.
✓✀✀A reminder:  ListPlot[rkpoints, AxesLabel->{"[X]","[Y]"}]  will plot the [X] and [Y] values
	✎✀Does the plot look like you expect?  Did you have any expectations?  (If not, try the next step!)
❏  The graphs of concentration versus time appear to be oscillatory.  Fourier transform the data to find the frequencies of the oscillation for both X and Y.  What do you notice?  Why do think you observe this?

Solving the differential rate equations

set up the constants

Ao=3.0;Xo=.85;Yo=3.2;
k1=1.0;k2=1.5;k3=1.2;
total=50;incr=.5;

solve the differential equations

rkpoints = Block[{x, y},
RK[{k1 Ao x - k2 x y,k2 x y - k3 y},
	{x, y}, {Xo, Yo}, {total, incr}]];

prepare the data in the form of [X] vs. t, etc.
Ignore this error!!!  General::spell1: 
   Possible spelling error: new symbol name "iline"
     is similar to existing symbol "rline".

points=Transpose[rkpoints];
reactant=points[[1]];
intermediate=points[[2]];
time=Table[i,{i,0,total,incr}];
rline=Graphics[Line[Transpose[{time,reactant}]]];
iline=Graphics[Line[Transpose[{time,intermediate}]]];

Show[{rline,iline},Axes->Automatic]

Variations in space!

ListPlot[rkpoints, AxesLabel->{"[X]","[Y]"}]

try it for sequential reactions!

A[t_]:=Ao Exp[-k1 t]
B[t_]:=(k1/(k2-k1)) (Exp[-k1 t] - Exp[-k2 t]) Ao
Pr[t_]:=Ao-B[t]-A[t]

Setting up the rate constants:

Ao=3.7;
k1=1.2;
k2=1.5;

create the data and plot it!

seqpts=Table[{A[t],B[t]},{t,0,total,incr}];
ListPlot[seqpts,AxesLabel->{"[A]","[B]"}]

periodic behavior (Fourier transform)

FT of the host data (H)

ftH=Abs[Fourier[reactant]];

hostplot=ListPlot[ftH, PlotJoined->True,PlotRange->{0,5}]

ftP=Abs[Fourier[intermediate]];

parasiteplot=ListPlot[ftP, PlotJoined->True,PlotRange->{0,5},
PlotStyle->{Thickness[.008 ],GrayLevel[.5]}]

Show[hostplot,parasiteplot]

Runge-Kutta package

Copyright Notice

Copyright 1989 by Roman Maeder.
	
	Adapted from	Roman E. Maeder: Programming in Mathematica, Addison-Wesley, 1989.
	
	Permission is hereby granted to make copies of	this file for any purpose other than direct profit,	or as part of a commercial product, provided this copyright notice is left intact.	Sale, other than for the cost of media, is prohibited.
	
	Permission is hereby granted to reproduce part or all of this file provided that the source is acknowledged.

Implementation

BeginPackage["RMPackages`RK`"]

RK::usage = "RK[{e1,e2,..}, {y1,y2,..}, {a1,a2,..}, {t1, dt}]
	numerically integrates the ei as functions of the yi with inital values ai.
	The integration proceeds in steps of dt from 0 to t1.
	RungeKutta[{e1,e2,..}, {y1,y2,..}, {a1,a2,..}, {t, t0, t1, dt}] integrates
	a time-dependent system from t0 to t1."

Begin["`Private`"]

RKStep[f_, y_, y0_, dt_] :=
	Block[{ k1, k2, k3, k4 },
		k1 = dt N[ f /. Thread[y -> y0] ];
		k2 = dt N[ f /. Thread[y -> y0 + k1/2] ];
		k3 = dt N[ f /. Thread[y -> y0 + k2/2] ];
		k4 = dt N[ f /. Thread[y -> y0 + k3] ];
		y0 + (k1 + 2 k2 + 2 k3 + k4)/6
	]

RK[f_List, y_List, y0_List, {t1_, dt_}] :=
	NestList[ RKStep[f, y, #, N[dt]]&, N[y0], Round[N[t1/dt]] ] /;
		Length[f] == Length[y] == Length[y0]

RK[f_List, y_List, y0_List, {t_, t0_, t1_, dt_}] :=
	Block[{res},
		res = RK[ Append[f, 1], Append[y, t], Append[y0, t0], {t1 - t0, dt} ];
		Drop[#, -1]& /@ res
	]  /;  Length[f] == Length[y] == Length[y0]

End[]

Protect[RK]

EndPackage[]