A wide variety of chemical reactions can be modeled with coupled (often nonlinear) differential equations. These equations describe the time evolution of the concentrations of the various chemical species: reactants, intermediaries, catalysts, and products. Such problems are quite simple to set up and solve with Mathematica. The function NDSolve numerically integrates the differential equations that arise. The resulting concentrations can be plotted as a function of time and can also be used to accurately compute the expected concentration of the molecular species.

Reaction with an intermediate

In the reaction A + 2  B⟶R + S, the mechanism proceeds through an intermediate species X.

({{A+B(Underoverscript[({{⟵}, {⟶}}), 2, 1])X}, {X+B(Overscript[⟶, 3])R+S}})}⟹A+2  B⟶R+S

The rate equations corresponding to this reaction are

({{(∂_t) a[t] ==  -(Subscript[k, 1]) a[t] b[t] + 
              (Subscript[k, 2]) x[t]}, {(∂_t) b[t] ==  -(Subscript[k, 1]) a[t] b[t] + 
              (Subscript[k, 2]) x[t] - (Subscript[k, 3]) b[t] x[t]}, {(∂_t) x[t] == 
            (Subscript[k, 1]) a[t] b[t] - (Subscript[k, 2]) x[t] - (Subscript[k, 3]) b[t] x[t]}})

where a[t], b[t], and x[t] represent the concentrations of the species A, B, and X at a time t at a constant temperature. The parameters (Subscript[k, 1]), (Subscript[k, 2]), and (Subscript[k, 3]), the rate constants, are temperature dependent. If the temperature of the mixture of reactants is allowed to change, then an additional set of equations is needed to model the temperature variation of (Subscript[k, 1]), (Subscript[k, 2]), and (Subscript[k, 3]) with time.

Although these equations do not have an explicit analytic solution, they can be numerically integrated using the built in numerical differential equation solver, NDSolve.  Here is the solution for a particular choice of the rate constants (Subscript[k, 1]), (Subscript[k, 2]), and (Subscript[k, 3]). In this example the initial concentrations of A and B are equal and that of the intermediary X is zero.

(Subscript[k, 1]) = 1/10;
(Subscript[k, 2]) = 1/10;
(Subscript[k, 3]) = 10;

ndsolution = 
NDSolve[{
(∂_t) a[t] == -(Subscript[k, 1]) a[t] b[t] + 
                (Subscript[k, 2]) x[t], (∂_t) b[
                  t] == -(Subscript[k, 1]) a[t] b[t] + (Subscript[k, 2]) x[t] - 
                (Subscript[k, 3]) b[t] x[t], (∂_t) x[t] == 
              (Subscript[k, 1]) a[t] b[t] - (Subscript[k, 2]) x[t] - 
                (Subscript[k, 3]) b[t] x[t], 
a[0] == 
              1, 
b[0] == 
              1, 
x[0] == 0}, 
{a, 
            b, x}, {t, 0, 
            10}][[1]]

The result is expressed in terms of a set of replacement rules that give the functions as InterpolatingFunction objects. The notation InterpolatingFunction[{{0.`, 10. }}, "<>"], for example, is a shorthand for the numerical information needed by the interpolation algorithm to reproduce the solution.

Here is plot of the result. The concentration of X was multiplied by a factor of 50 so that it can be seen in the graph on the same scale as A and B). The colors of the curves for A, B, and X, are black, red, and green respectively.

Plot[Evaluate[{a[t], b[t], 50  x[t]} /. ndsolution], {t, 0, 10}, 
        PlotStyle -> {{AbsoluteThickness[2], 
              RGBColor[0,  0,  0]}, {AbsoluteThickness[2], 
              RGBColor[ .7,  0,  0]}, {AbsoluteThickness[2], 
              RGBColor[0,   .7,  0]}}, Axes -> False, 
        Frame -> True, 
        PlotLabel -> 
          StyleForm[
            A StyleForm[" B", 
                FontColor -> RGBColor[ .7, 0, 0]] StyleForm[" X", 
                FontColor -> RGBColor[0,  .7, 0]] , 
            FontFamily -> "Helvetica", 
            FontWeight -> "Bold"] ];

For this example, as the reaction proceeds, the intermediary is created and then settles down to an approximately steady state value.

Reaction with Catalysts

({{A + X(Overscript[⟶, 1])R + Y}, {B+Y(Overscript[⟶, 2])S+X}})}⟹A+B⟶R+S

The rate equations corresponding to this reaction are

({{(∂_t) a[t] ==  -(Subscript[k, 1]) a[t] x[t]}, {(∂_t) b[t] ==   -(Subscript[k, 2]) b[t] y[t]}, {(∂_t) x[t] ==  -(Subscript[k, 1]) a[t] x[t] + 
              (Subscript[k, 2]) b[t] y[t]}, {(∂_t) y[t] == 
            (Subscript[k, 1]) a[t] x[t] - (Subscript[k, 2]) b[t] y[t]}})

Here is a particular case of the integration of these differential equations. The initial concentrations of A and B are equal and the rate constants (Subscript[k, 1]) and (Subscript[k, 2])chosen as shown.

(Subscript[k, 1]) = (1/10);
(Subscript[k, 2]) = (1/20);
ndsolution = 
NDSolve[{
(∂_t) a[t] == -(Subscript[k, 1]) a[t] x[
                  t], 
(∂_t) b[
                  t] == -(Subscript[k, 2]) b[t] y[
                  t], 
(∂_t) x[
                  t] == -(Subscript[k, 1]) a[t] x[t] + 
                (Subscript[k, 2]) b[t] y[t], 
(∂_t) y[
                  t] == 
              (Subscript[k, 1]) a[t] x[t] - (Subscript[k, 2]) b[t] y[t], 
a[0] == 
              1, 
b[0] == 
              1, 
x[0] == 
              1, 
y[0] == 0}, 
{a, 
            b, x, y}, {t, 0, 
            250}][[1]]

Here is plot of the result. The colors of the curves for A, B, X, and Y are black, red, green, and blue respectively.

Plot[Evaluate[{a[t], b[t], x[t], y[t]} /. ndsolution], {t, 0, 250}, 
        PlotStyle -> {{AbsoluteThickness[2], 
              RGBColor[0,  0,  0]}, {AbsoluteThickness[2], 
              RGBColor[ .7,  0,  0]}, {AbsoluteThickness[2], 
              RGBColor[0,   .7,  0]}, {AbsoluteThickness[2], 
              RGBColor[0,  0,   .7]}}, Axes -> False, 
        Frame -> True, 
        PlotLabel -> 
          StyleForm[
            A StyleForm[" B", 
                FontColor -> RGBColor[ .7, 0, 0]] StyleForm[" X", 
                FontColor -> RGBColor[0,  .7, 0]] 
              StyleForm[" Y", FontColor -> RGBColor[0, 0,  .7]], 
            FontFamily -> "Helvetica", 
            FontWeight -> "Bold"]];

It can be seen that, as the reaction proceeds towards completion, the concentrations of the catalysts X and Y eventually return to their initial values.