Reaction Kinetics and Differential Equations

Introduction

Differential equations (DEs) are a universal feature of the modelling of processes in chemistry. In this Lab we introduce Mathematica's handling of DEs in the context of reaction kinetics. We'll be mostly concerned here with solving DEs analytically, that is, producing explicit algebraic formulae for the solutions. Unfortunately, though, analytical techniques simply don't work for most of the interesting DE models in physical chemistry (including, for example, almost all practical quantum mechanics calculations). You need to know something about Mathematica's numerical solving capabilities, therefore, which we make a start on in the final section of this Lab, and which you'll make use of extensively in the next Lab.

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).

1.1 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 name for use in setting up DEs), so let's create a function called "asol":

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.)

In kinetics, we're usually not much interested in seeing the general solution to a rate law DE with its arbitrary constant of integration. Instead, we specify an initial condition for the concentration as part of the call to DSolve, inserting it as an extra equation 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:

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

This form of the solution is useful for exploring the effects of varying the initial concentration and rate constant. For example, here is a plot of a(t) for some typical real values of a0 (in mol (dm^(-3))) -- and k ((s^(-1))):

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

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

Exercises and Explorations

1: Solving rate laws with DSolve

(i) Using DSolve, find solutions (in terms of an initial concentration a0) to all the rate laws given in the table above. Satisfy yourself that Mathematica's solutions are equivalent to the integrated rate laws stated -- the easiest way to do this is to start with the "k t = ..." form and rearrange it to make a the subject. It turns out that Solve is the best command for doing this; for example

Solve[k t  ==  (1/a) - (1/a0), a]

This won't work with the first-order rate law, however; for that it is best to start with "a =..." and rearrange to make t the subject:

Solve[a  ==  a0 (E^(-k t)), t]

(ii) Suppose that you had obtained some experimental data for a reaction A -> Products where (Subscript[([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.

Note: if you have assigned values to a0 and k, you should "clear" them before starting Question 2. Simply type

Clear[a0, k]

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]))))    or   \(TraditionalForm\`a = \(a\_0\ \((a\_0 - b\_0)\) e\^\(a\_0\ k\ t\)\)\/\(\(a\_0\) e\^\(a\_0\ k\ t\) - \(b\_0\) e\^\(b\_0\ k\ t\)\)\)

A + 2B -> Products  ........  k t =(1/((Subscript[b, 0]) - 2  (Subscript[a, 0])))ln((((Subscript[a, 0]) b)/(a (Subscript[b, 0]))))  or   \(TraditionalForm\`a = \(a\_0\ \((2  a\_0 - b\_0)\) e\^\(2  a\_0\ k\ t\)\)\/\(2 \( a\_0\) e\^\(2  a\_0\ k\ t\) - \(b\_0\) e\^\(b\_0\ k\ t\)\)\)

where a and b denote [A] and [B] at time t, and (Subscript[a, 0]) and (Subscript[b, 0]) are the initial concentrations. Of course, a and b are related through b =  (Subscript[b, 0]) - ((Subscript[a, 0]) - a) in the first case, and b =  (Subscript[b, 0]) - 2 ((Subscript[a, 0]) - a) in the second.

(i) Derive the integrated rate laws (in the form "a = ...") for the two different stoichiometries using suitable calls to DSolve. Hint: express b[t] in terms of a[t], a0 and b0.

(ii) Plot graphs of a(t) for the two different stoichiometries, taking k = 100 (dm^3) (mol^(-1)) (s^(-1)) and t ranging from 0 to 100 s, and for the following choices of (Subscript[a, 0]) and (Subscript[b, 0]):
                   (a)  (Subscript[a, 0])=(10^(-3))  and  (Subscript[b, 0]) = 0.5×(10^(-3))  mol (dm^(-3))
                   (b)  (Subscript[a, 0])=(10^(-3))  and  (Subscript[b, 0]) = 2.5×(10^(-3))  mol (dm^(-3))
Explain your graphs in terms of the reaction chemistry.

(iii) For the A + B -> P reaction, it appears that when (Subscript[b, 0]) and (Subscript[a, 0]) become equal, a(t) is equal to (0/0), which is sensible neither physically nor mathematically! Go back to the rate law and solve it again for the case (Subscript[b, 0]) = (Subscript[a, 0]). Comment on the result.

(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])) formula for a(t) is equal to the alternative formula 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. In the following exercises, we ask you first to construct a symbolic solution for these coupled rate laws, and then to use Mathematica's numerical solver function NDSolve to construct an approximate numerical solution.

Exercises and Explorations

1: Constructing the rate equations

Construct the rate laws for d[A]/dt, d[B]/dt and d[C]/dt as Mathematica DEs for the above reaction. We've made a start for you ...

a'[t] == -k1 a[t]
b'[
              t]== ...  ...  ... 
c'[t]== ...  ...  ...

Now produce symbolic solutions using DSolve. The output can be simplified slightly, hence we call Simplify.

DSolve[{... INSERT YOUR RATE LAWS HERE ... , a[0] == a0, b[0] == 0, c[0] == 0}, {a[t], 
        b[t], c[t]}, t]

exactsolution = Simplify[%]

2: The exact solution

Use the exactsolution rules to define Mathematica functions for a(t), b(t) and c(t) as before. That is,

aexact[t_]  =  a[t]  /.  First[exactsolution]

and so on.

Hence, create a plot that shows on the same axes the graphs of a(t), b(t) and c(t) against 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)). Use different plotting styles so that the different graphs can be distinguished.

3: The approximate solution

Now for a numerical solution: the inputs to NDSolve are similar to those for DSolve, except that (Subscript[a, 0]), (Subscript[k, 1]) and (Subscript[k, 2]) must be given numeric values, and a specific time range must be given.

numericalsolution  = 
      NDSolve[{... INSERT YOUR RATE LAWS HERE ... , a[0] == (10^(-3)), b[0] == 0, 
            c[0] == 0} /. {k1 -> 1, k2 -> 1.5}, {a[t], b[t], c[t]}, {t, 0, 
          10}]

The solutions are given in the form of three "hidden" functions of t, and are defined only for 0 ≤ t ≤ 10. If you ask for function values outside this range you will get error messages. Define Mathematica functions for these numerical results using

anumer[t_]  =  a[t]  /.  First[numericalsolution]

and so on. Hence, create a plot of these functions as you did for the exact functions in part (ii).

4: Two regimes

(iv) Generate some plots -- using either exact or numerical solutions as you wish -- to show the characteristic behaviours for reactions where (Subscript[k, 1]) ≫ (Subscript[k, 2]) or (Subscript[k, 1]) ≪ (Subscript[k, 2]) -- that is, illustrate the effects due to having a "rate-determining step".

5: Numerical precision

Numerical solutions for DEs are (in general) nowhere near as reliable as symbolic ones, and should always be thoroughly checked. Devise a test for the comparison of the numerical and exact solutions, and apply it first to the numerical solution from part3, and then to a numerical solution for your own choice of values for (Subscript[a, 0]), (Subscript[k, 1]), (Subscript[k, 2]) and time range. 

On the basis of your test, how accurate (in terms of an approximate number of decimal places) do you think NDSolve's solutions are? (What do the Help System entries for NDSolve say about this?)