Chemical Reactions

We will study a standard kind of chemical reaction, model it with a differential equation and use Mathematica to describe the eventual outcome of the reaction. Typically, two or more substances combine to form one or more new substances. Suppose two chemicals combine to form a compound X. Start with A grams of the first chemical and B grams of the second. We will denote by x(t) the amount of the new compound after t minutes. We observe that the reaction between the two chemicals is such that the proportion A:B of chemicals used in forming X is u:1-u, where 0 < u < 1. For example, u = 1/3 means that for every gram of the first chemical, two grams of the second are required. The basic mathematical model of the reaction stipulates that the rate at which the amount x(t) is changing is proportional to the product of the residual amounts of the two chemicals still present. The rate of change of x(t) is its derivative x'(t), and since it takes uX grams of the first, and (1 - u)X grams of the second to produce X grams of the new substance, we can express that mathematically by the differential equation
x'(t) = k (A - ux(t))(B - (1 - u)x(t)).
The constant of proportionality k is not known at the outset, but can be determined by measurement, as we shall see. Let's set Mathematica loose on the equation.

sol = DSolve[{x'[t] == k*(A - u*x[t])(B - (1 - u)*x[t]),
 x[0] == 0}, x[t], t]

We used the initial condition x(0) = 0, since we assume there are no grams of X initially. Next we extract a function:

f[t_, A_, B_, u_, k_] = x[t] /. First[sol]

Now let us look at a couple of concrete examples.

Example 1

Suppose that 50 grams of the first substance and 32 grams of the second are present initially. Suppose also that it takes 4 grams of the latter for each gram of the former in the reaction. Finally, suppose we measure 30 grams of X after 10 minutes. Then we can use this data to solve for the coefficient of proportionality k.

coeff = FindRoot[f[10, 50, 32, 0.2, k] == 30, {k, 0.005},
 MaxIterations -> 20]

(Note: We had to increase the default number of iterations in FindRoot to do the analysis.)
Now we define a specific function corresponding to this data, check that it has the correct value for t = 10, compute its long-term behavior, plot it, and compute the point where the reaction is 99% finished.

f1[t_] = f[t, 50, 32, 0.2, k/.First[coeff]] //FullSimplify

f1[10] //N

Limit[f1[t], t -> Infinity]

Plot[f1[t], {t, 0, 40}, PlotRange -> {0, 40},
 AxesLabel -> {"Time","Quantity"}];

FindRoot[f1[t] == 0.99*40, {t, 35}]

Example 2

Now assume that there are still 50 grams of the first substance, but 80 grams of the second. Assume u = 0.6, that is, it takes 3 grams of the former for each 2 of the latter for the reaction to proceed. Once again, assume that 30 grams of the new substance are present after 30 minutes. We carry out the same analysis.

coeff2 = FindRoot[f[10, 50, 80, 0.6, k] == 30, {k, 0.005},
 MaxIterations -> 20]

f2[t_] = f[t, 50, 80, 0.6, k/.First[coeff2]]//FullSimplify

f2[10] //N

lim2 = Limit[f2[t], t -> Infinity]

FindRoot[f2[t]==0.99*lim2,{t,74}]

The picture looks flat around t = 80 minutes, but this reaction is slower than the first--it takes 143.543 minutes to complete 99% of the process.