In homework #6, we take on a fun challenge: find the best-fit reaction orders for a reversible reaction,

H₂SO₄ + (C₂H₅)₂SO₄ ß à 2C₂H₅SO₄H

Given the following data (from Hellin and Jungers (1957), more recently from Levenspiel (1972)):

Time, min	C2H5SO4H, mol/liter	Time, min	C2H5SO4H, mol/liter
0	0	180	4.11
41	1.18	194	4.31
48	1.38	212	4.45
55	1.63	267	4.86
75	2.24	318	5.15
96	2.75	368	5.32
127	3.31	379	5.35
146	3.76	410	5.42
162	3.81	infinity	(5.80)

This data gives us two important pieces of information: the initial rate of reaction, and the chemical equilibrium. We can calculate the initial rate of reaction by fitting a polynomial curve to the data, recognizing that

, where

Using Excel or Polymath or Matlab, we can get the best fit polynomial to the experimental data,

If we take the derivation of this equation, and calculate at t = 0, we have fit an initial rate of reaction, as follows:

at t= 0.

Thus, mol/L.min.

Okay, now the other piece of information we get is the equilibrium. If we assume arbitrary reaction orders

(we don’t know them) then the equilibrium coefficient can be calculated from the rate of reaction. But first,

cbcause mol/L, and the stoichiometry is identical, and the rate expression can

be written in terms of a lumped forward reaction order, as follows:

, where a' = a + g.

Now, the equilibrium coefficient can be obtained given values for alpha (forward reaction order) and beta

(reverse reaction order), as follows:

Because we know the initial rate of reaction, which accounts for ONLY the forward reaction (mol/L),

we can also calculate the rate coefficient, k:

, or

So, now if we guess a forward reaction order and a reverse reaction order, we can calculate from the data values for

K_eq and k. We can then solve for the concentration of C vs. t using our ASSUMED kinetics, and then compare to

the experimental data.

Guess #1: Assume Elementary Kinetics: a = 1, b = 1, g = 2 (a’ = 2, g = 2).

K_eq = 4.976, and k = 0.000577. In Matlab, we plot up the data vs. the resulting assumed model.

We can calculate the total error of this fit, as follows:

, for 17 data points being fit using 2 degrees of freedom. The error is normalized by the experimental data, so that error at large concentration values (tà inf) weigh as much as error at low concentrations (t à0) for a more balanced fit.

For this guess, the total error, s = 0.252.

Guess #2: Assume Non-Elementary Kinetics: a = 1, b = 1, g = 1 (a’ = 2, g = 1).

K_eq = 0.858, and k = 0.000577. In Matlab, we plot up the data vs. the resulting assumed model.

We can see that this fit is worse. Calculation of the error agrees, as s = 0.418.

Guess #3: Assume Non-Elementary Kinetics: a = ½, b = ½, g = 1 (a’ = 1, g = 1).

In Matlab, we plot up the data vs. the resulting assumed model.

This is a little closer, and the error for this guess is s = 0.2147.

Guess #4: Assume Non-Elementary Kinetics: a = ½, b = ½, g = 2 (a’ = 1, g = 2).

In Matlab, we plot up the data vs. the resulting assumed model.

This is getting even better – the error is now s = 0.130.

Guess #5: Assume Non-Elementary Kinetics: a = ½, b = ½, g = 1.5 (a’ = 1, g = 1.5).

In Matlab, we plot up the data vs. the resulting assumed model.

The error is now down to s = 0.077. Not bad.

We have some matlab files that automate this whole process, as follows:

main.m : the main program, which runs the calculations described below, using the following subroutines:

geterror.m: given a guess in the form ‘[alpha beta]’, solves for C_C vs. t and then calc’s error of fit, and

solve_eqns.m, which contains the differential equations governing C_A and C_C (recognizing that C_A = C_B).

The first thing that our automated Matlab file does is calculate the error of fit for a 15x15 grid of values for alpha and beta, and then plots the error surface, using first a contour plot and then a surface plot in 3-D.

Now we can visualize the error surface over a range of reaction orders – we see that there isn’t one blatantly clear answer; there is a ‘troth’ along this two-dimensional space that gives a pretty-good fit. We can get a bit of a better look by getting a contour plot of the log-error,

We can see that there is our best-fit kinetics. So, the forward reaction is half-order w.r.t Sulfuric Acid, and half-order w.r.t. Diethylsulfate, while the reverse reaction is 1.5^th order w.r.t. Ethylsulfate.

This doesn’t seem too intuitive, at least not just yet. What if we hypothesize that reaction is first-order w.r.t. Sulfuric Acid, first-order w.r.t. Diethylsulfate, and 3^rd-order w.r.t. Ethylsulfate?

Guess #N: Assume Non-Elementary Kinetics: a = 1, b = 1, g = 3 (a’ = 2, g = 3).

In Matlab, we plot up the data vs. the resulting assumed model.

Well, this is okay, it looks alright, but it’s not as good a fit as our solution.

If you’d like to manually guess orders, and then see what the fit looks like, you can use the following codes:

one_run.m, which also needs the file

solve_eqns.m, which contains the differential equations for dCc/dt.

Just download these files into your working directory for Matlab (see Tutorial #1), and then edit the one_run.m file to change the reaction orders.