Reaction Kinetics and Differential Equations

The METRIC Project

Introduction: Kinetic problems and differential equations

Differential equations (DEs) are a universal feature of the modelling of processes in chemistry. For example, in the last Maths Lab, we used Newton's second law of motion to set up three DEs as a model for the interaction of three ions, and then asked Mathematica's numerical solver, NDSolve, to generate solutions for the ion motions from a rather formidable set of coupled, second-order differential equations (which are, by the way, an instance of one of the great unsolved problems of applied mathematics, the three body problem).

As you've seen, it is one thing to write down a DE model for some process, and quite another thing to solve the DEs. Today, we are going to examine these issues in the context of reaction kinetics, considering the solution of DEs in Mathematica by both analytical and numerical methods. The next few paragraphs are given to remind you of the basic terminology of kinetics.

In kinetics, DEs form a link between two aspects of the chemistry: empirical reaction data, and kinetic theory. The DE is an expression for the rate of change of concentrations of reactants and products, and the solutions to the DE give the model's prediction for the time behaviour. A rate law--also known as a rate equation--is a DE that relates reaction rate to the concentrations of chemical species present.  For example, in a simple bimolecular reaction A + B -> C, the rate law is

((d[A])/dt)=-k [A]^α [B]^β     where     ((d[A])/dt)=((d[B])/dt)=-((d[C])/dt)

The number k is the rate constant, the numbers α and β (which need not be integer) are the orders of the reaction. The overall order is the sum of the individual orders. Confusingly, in the mathematics of DEs the term order has a quite different meaning: it refers to the highest derivative appearing in an equation. So for example if ((d^2)/(d (t^2)))is the highest derivative the DE is said to be of second order.

1. Simple kinetic models and DSolve

The rate laws and their solutions ("integrated rate laws") for a reaction A -> Products are normally tabulated in this form, where a(t) = [A] at time t and (Subscript[a, 0]) = the initial concentration:

({{Order, Rate law, Integrated rate law}, {0, ((d a)/(d t))=-k, k t=(Subscript[a, 0]) - a}, {1, ((d a)/(d t))=-k a, k t=ln((((Subscript[a, 0]))/a))}, {2, ((d a)/(d t))=-k (a^2), k t=(1/a)-(1/(Subscript[a, 0]))}, {3, ((d a)/(d t))=-k (a^3), k t=(1/(2(a^2)))-(1/(2((Subscript[a, 0])^2)))}, {n≥2, ((d a)/(d t))=-k (a^n), k t=(1/(n - 1))((1/((a^(n - 1))))-(1/((Subscript[a, 0])^(n - 1))))}})

Notice that this is not how the solutions are conventionally written in mathematics, or are expressed by Mathematica, where they are always written with a(t) as the subject--for example, the first-order solution would be written as a(t) = (Subscript[a, 0]) (e^(-k t)).

This traditional chemist's way of writing makes a lot of sense in the context of doing a kinetic analysis manually, when a quick, approximate answer is required. It says that a certain quantity involving a and (Subscript[a, 0]) is linearly related to the time t, and that the gradient of the line is k. This is useful because linear relationships are the easiest things to spot, and to calculate by simple statistical techniques (such as are found on a pocket calculator).

DSolve: Mathematica's symbolic DE solver

The first-order rate law is written in Mathematica as:

a'[t]  ==  -k a[t]

Notice the "prime" notation (a'(t) being the same thing as da/dt), and the need to indicate explicitly the dependence of a upon t. The general solution to this equation is given by the following application of DSolve:

solution1=DSolve[a'[t]  ==  -k a[t],  a[t],  t]

DSolve takes three arguments: the equation(s) to be solved, the dependent variable(s), and the independent variable. In the solution, constants of integration appear as C[1], C[2], etc.

Getting this solution into a usable form is a surprisingly (annoyingly!) intricate task. As with solutions to algebraic equations given by Solve or FindRoot, the solution here is given as a rule, which has to be used in a substitution operation:

a[t]  /.  First[solution1]

The First command is required to remove the redundant outer pair of curly brackets in solution1. Now, we'd quite like a(t) as a proper Mathematica function (but we'll retain "a[t]" itself as an unassigned variable for use in expressing DEs):

asol[t_]  =  a[t]  /.  First[solution1]

(Note: this breaks our general recommendation about always using ":=" to define functions; "=" is preferable here because we want the substitution to happen right away; with := the substitution would happen repeatedly every time the function was called.)

Here we plot a(t) for some typical real values of C[1] (which equals the initial concentration in this case, in mol (dm^(-3))) and k ((s^(-1))):

\(Plot[
      asol[t] /. {C[1] -> 10\^\(-3\), k -> 2.0\ 10\^\(-4\)}, {t, 0, 10000}, 
      AxesLabel -> {"t(s)", \*""[A](mol \!\(dm\^\(-3\)\))""}]\)

We can make the solving process a bit more elegant by specifying an initial concentration in the call to DSolve: when you want such a particular solution, you just specify the initial condition by means of a second equation (n.b. with "==") in the first argument (which thus becomes a list). For example

DSolve[{a'[t]  ==  -k a[t], a[0] == (10^(-3))}, a[t], t]

The initial condition need not be numeric, however:

DSolve[{a'[t]  ==  -k a[t],  a[0]  ==  a0},  a[t],  t]

Exercises and explorations 1

1. Finding integrated rate laws with DSolve

(i) Using DSolve, find solutions to all the rate laws given in the table above. Verify that Mathematica's solutions are equivalent to the integrated rate laws stated; note that you can make t the subject of an equation using Solve--for example, using solution1 from above:

Solve[a == (a[t] /. First[solution1]), t]

(ii) Suppose that you had obtained some experimental data for a reaction A -> Products where \(TraditionalForm
      \`\([A]\)\_0 = \(10\^\(\(-3\)\ \)\) mol\ dm\^\(-3\)\) and the reaction had proceeded to 75% completion after 30 minutes. On the same axes, plot graphs of a(t) for the zeroth, first, second and third-order integrated rate laws which "fit" the time scale of the observed reaction. State the value of k (with units, please!) which you use in each case.

2. Second order reactions

The rate law for a bimolecular, overall second order reaction is

((d [A])/dt)=-k [A] [B]

However, the integrated rate law depends on the stoichiometry of the reaction, thus:

A + B -> Products  .................  k t =(1/((Subscript[b, 0]) - (Subscript[a, 0])))ln((((Subscript[(a ), 0]) b)/(a (Subscript[b, 0]))))

A + 2B -> Products  ...............  k t =(1/((Subscript[b, 0]) - 2 (Subscript[a, 0])))ln((((Subscript[(a ), 0]) b)/(a (Subscript[b, 0]))))

where a and b denote [A] and [B] at time t, and (Subscript[a, 0]) and (Subscript[b, 0]) are the initial concentrations.

(i) Derive these two integrated rate laws using DSolve. (Hint: express b(t) in terms of a(t), (Subscript[a, 0]) and (Subscript[b, 0]).) Note that if you wish to use Solve for the rearrangement to "k t =  ..." you may need to help Mathematica out by manually cancelling some common terms.

(ii) For the A + B -> P reaction, it appears that the integrated rate law "blows up" when (Subscript[b, 0]) and (Subscript[a, 0]) become equal, which we know cannot be physically sensible. Go back to your rate law and solve it again for the case (Subscript[b, 0]) = (Subscript[a, 0]).

(iii) Make a plot on the same axes of graphs of a(t) for the cases (Subscript[b, 0]) = 0.9 (Subscript[a, 0]),  (Subscript[b, 0]) = (Subscript[a, 0]) and (Subscript[b, 0]) = 1.1 (Subscript[a, 0]); take k = 100 (dm^3) (mol^(-1)) (s^(-1)) and (Subscript[a, 0]) = (10^(-3)) mol (dm^(-3)), and use different plot styles to distinguish the three graphs. Produce one or two further plots to show that the graphs for (Subscript[b, 0]) ≠ (Subscript[a, 0]) do in fact converge onto the (Subscript[b, 0]) = (Subscript[a, 0]) graph as (Subscript[b, 0]) approaches (Subscript[a, 0]) from above and below. 

(Note: if you know something about limits in mathematics, you should be able to show, perhaps with Mathematica's Limit command, that in the limit of (Subscript[b, 0]) -> (Subscript[a, 0]) the general ((Subscript[b, 0]) ≠ (Subscript[a, 0])) expression for a(t) is equal to the alternative expression for when (Subscript[b, 0]) = (Subscript[a, 0]).)

2. Consecutive reactions

Let's consider the simplest case of the consecutive reactions:

A  (Overscript[-> , (Subscript[k, 1])]) B  (Overscript[-> , (Subscript[k, 2])]) C

where both reactions are first order and [B] = [C] = 0 at time 0. The coupled rate laws for these reactions can be solved either symbolically or numerically. In the following exercises, we ask you to construct both types of solution and to compare the results.

Exercises and explorations 2

(i) Construct the rate laws for d[A]/dt, d[B]/dt and d[C]/dt as Mathematica DEs for the above reaction. Denote concentrations in Mathematica by a[t] for [A] at time t, etc. Using subscripts in Mathematica DEs is a bit problematic, so it's easiest to use k1 and k2 to stand for the rate constants.

(ii) Using code of the form given in Maths Lab 6 (Newtonian dynamics), type in a call to NDSolve that generates numerical solutions for a(t), b(t) and c(t), for a time range t =0 up to t = 10, a(0) = (10^(-3)) mol (dm^(-3)) , (Subscript[k, 1])= 1.0 (s^(-1)) and (Subscript[k, 2])=1.5 (s^(-1)).

(iii) Now (again following the pattern of Maths Lab 6), use those solutions to create a plot that shows on the same axes the graphs of a(t), b(t) and c(t) against t. Please use different plotting styles so that the different graphs can be distinguished.

(iv) Using DSolve, or working by hand if you prefer, construct the exact, symbolic expressions for a(t), b(t) and c(t) in terms of (Subscript[k, 1]), (Subscript[k, 2]) and (Subscript[a, 0]) = (Subscript[([A]), 0]).  Devise a procedure that plots the graphs of the concentrations against t using these symbolic expressions. Devise a test for the equivalence of the numerical and exact solutions, and apply it to the numerical solution found in part (ii).

(v) Generate some plots showing the characteristic behaviours for reactions where (Subscript[k, 1]) ≫ (Subscript[k, 2]) or (Subscript[k, 1]) ≪ (Subscript[k, 2]).

(vi) Looking at your exact solutions for the DEs, it would appear that the b(t) and c(t) solutions "blow up" for the case of (Subscript[k, 1]) = (Subscript[k, 2]), which we know is not physically sensible, and also you can check that the numerical solutions are perfectly well-behaved in this situation. Go back to your rate laws and re-solve them for when (Subscript[k, 2]) = (Subscript[k, 1]). Make plots which compare these solutions with (I) the numerical solutions, and (II) the general ((Subscript[k, 2]) ≠ (Subscript[k, 1])) solutions for when (Subscript[k, 2]) is just below and just above (Subscript[k, 1]) (cf. the previous exercise on A + B -> P).

(Note: again, if you know something about limits in mathematics, you should be able to show that in the limit of (Subscript[k, 2]) -> (Subscript[k, 1]) the general ((Subscript[k, 2]) ≠ (Subscript[k, 1])) expressions for the solutions are equal to the alternative expressions for when (Subscript[k, 2]) = (Subscript[k, 1]).)