Using Grid Searches to Optimize a Fit

  We were able to estimate the reaction rate constants from functions of three unknown parameters and to improve the fit with a simple search for k₀. For a better fit we need to search for all three unknown parameters. We can choose X₀, X_eq and k₀ as the unknown parameters with X = A - A₀ being the amount of product at some particular time. X simplifies the formulas which is the reason why it was chosen. We have estimates for X_eq and k₀ determined by the data points and can search the neighborhood of the values for the best fit. So a three dimensional grid search will allow us to improve on our initial estimates. Here is a sample of how one might start out. The first three lines determine the bounds of the search and the next three compute the grid points. X₀ was initially assumed to be a small positive number.

I wrote a Mathcad program to do the search for the integers λ, μ and ν which give the lowest value for the variance, V, of A and A'.

Three such searches zoomed in on values very close to the original values.

They give a relative rms error for the fit of about three tenth of one percent. The fit of the "lower left corner" is also much better.

An error propagation formula will help us determine how errors in the grid search parameters will affect the final results. We can estimate how two independent variables, x and y, affect a third dependent variable, z, starting with the chain rule. If dx and dy represent small errors then the chain rule gives us a good estimate of dz. The formula can be modified to give the rms error of dz but we need to use an inequality to express it in terms of the rms errors of dx and dy. One can choose the right hand side of the last equation to represent a bound for the rms error of dz. One usually sees this propagation of error formula written as an equality.

The formula for the propagation of error can be used to give the expected relative errors of the computed parameters for the fit based on the relative errors of the fit variables. The matrix Λ was found by applying the propagation of error formula to each of the computed variables and evaluating it using estimated values. The differences between the fit values and original values can be used to get approximations for the relative errors of the fit variables. The bottom left shows the expected bounds and can be thought of as standard deviations. The bottom right shows the z-scores for the the actual relative errors.

The grid search we used isn't optimal but it allows us to improve on the original estimate which a plot showed to be slightly off. The final result is what can be expected from a grid search.

Supplemental (Nov 12): The program p_min(V) puts the rows of a table containing the pointer p and corresponding V values in the order of the V values. The resulting location for the minimum value of V is the first p value on the new table.

Fit For When Two Reactants are Consumed

  The fit for problem worked in the blog Increasing the Rate of Dissociation by Adding Acid is similar to that in the last blog. It's what one would have to do for the case where the two reactants are consumed in the reaction. There was a minor complication with the complex κ.

The estimate value for κ turned out to be close to the complex conjugate of the original κ. For z = x + i y the complex conjugate is z* = x - i y. It turns out that tanh(κ*) = tanh(κ) when the result is real. There is another formula for the tanh of a complex number involving ordinary and hyperbolic sines and cosines that helps explain the situation. Either the cosine or the sine of the imaginary part in the result has to be equal to zero.

Some data sets are more difficult to fit than others. One has to play with the numbers of terms in the empirical fit for best results. The empirical fit for the initial rate of change might be good but that doesn't guarantee that the fit for the formula will be good. The rounded corner at the lower left of the curve below seems to be hardest to get right. The procedure is sensitive to random errors in the first few data points so some method of adjusting errors might prove useful.

Supplemental (Nov 8): A weakness of the fit procedure is that an error in k'₀ propagates through as one calculates all the other fit constants so if one searches the neighborhood of k'₀ for the value which gives a minimum for the variance of the difference between A' and A then one will get consistently good fits to the data points. The resulting fit is a least squares fit of the data.

Fitting a Simple Catalyzed Reaction

  In a catalyzed reaction the catalyst is not used up in the reaction. With the dissociation of sucrose in an acid solution water takes part in the reaction but the quantity of acid and its ions remain constant. The quantity of hydrogen ions in solution is maintained by the dissociation of water. We can take this into account and it simplifies the equations for the constants in the formula for the amount of sucrose as a function of time. We proceed as we did before determining the coefficients for the quadratic equation for the reaction rate but let B remain unchanged. Integrating the reaction rate equation again produces a formula with a hyperbolic tangent function but the constants in it are slightly different due to changes in relative magnitude. We can assume some initial conditions for the sake of argument.

We then generate some random data to work with and allow sufficient time for equilibrium to be reached.

If one checks the initial reaction rate one gets dA/dt = −k₀B₀A which agrees with the formula that Wilhelmy used so the initial solution is approximately an exponential decay function. To make the change easier to fit we can use the natural log of A instead of A and fit a polynomial equation. F(t) below is a vector function containing the powers of t and is needed by linfit to find the coefficients of the polynomial. We can estimate the initial value for A by averaging the values of A divided by powers of e with x values corresponding to the data points. A check of the standard deviation shows that the rms error is about one unit.

Once we have the fit we can draw a smooth curve through the data points using the empirical power law.

By taking the derivative of the fit polynomial at t = 0 we get the initial rate constant and this gives us an estimate of the first rate constant, k₀, on division by B₀. To find the second rate constant, k₁, we need an estimate of the equilibrium value for A and the equilibrium condition k₀B₀A_eq = k₁(A₀ - A_eq)² which states that the amount of product that is produced is canceled by the reverse reaction.

With the estimates for the rate constants we can now estimate the constants needed to fit the formula to the data. Using an auxiliary constant γ makes it easier to determine an estimate for κ. The fit agrees fairly well with the initial assumptions.

Increasing the Rate of Dissociation by Adding Acid

  Sucrose is a disaccharide produced by the hydrolysis of glucose and fructose molecules in which an hydroxyl group (-OH) in each sugar molecule join producing an oxygen bridge between them by releasing an H₂O water molecule. One can see the oxygen atom linking the hexagonal ring (glucose) to the pentagonal ring (fructose) in the .gif below.

animation credit: Wikipedia article Sucrose

Dissociating a disaccharide molecule requires a hydration reaction to replace the missing atoms of the simple sugars. There are always hydrogen and hydroxyl radicals normally present in a sugar solution which can react with the sucrose molecules but adding an acid as Wilhelmy did acts as a catalyst to increase the reaction rate.

  To study the effect of adding the acid to the solution we need to modify the equation for the reaction and deduce the corresponding reaction rate equation. In the reaction below we can let A correspond to the sucrose molecule, B to the hydrogen molecules present in the solution resulting from the addition of acid and C and D to the glucose and fructose products. We can also let X represent the progress of the reaction in order to determine the amount of sucrose and other molecules present at any given time. The acid molecules easily dissociate so the amount of hydrogen present is essentially equal to the concentration of acid. Plugging the modified amounts into the reaction rate formula again reduces it to a first order differential equation involving the sucrose concentration [A].

The reaction again follows a binomial rate law and we can simplify it by "completing the square" and rewriting the equation using the auxiliary parameters u and v.

The result is a similar formula involving the hyperbolic tangent function but here λ is replaced by βc so the net effect is that the reaction rate constants are multiplied by the factor β increasing the reaction rate considerably. We can check our math by working an example. From the positive polynomial coefficients a, b and c in the rate equation we can compute the constant parameters α and β and κ. Here κ is complex but that isn't a problem since we can use tan(iθ) = i tanh(θ) with a complex θ. The complex κ is easily handled by Mathcad. The complex value is the results from ([A]₀ + α)/β being greater than 1.

The results from an iterative evaluation of the reaction rate formula and those from the derived formula produce identical curves for the reaction, justifying the use of the derived formula. Setting the reaction rate polynomial to zero allows us to compute the equilibrium sucrose concentration A_eq which agrees with the result of the numerical computation.

Supplemental (Nov 4): Try k₀=1×10^-4, k₁=1.234567901234567×10<sup>-6

Suppemental (Nov 5): For the hyperbolic tangent of a complex number one can use the following identity which is similar to the identity for tanh(u+v).</sup>

Even though κ is complex, tanh(κ) is real and consequentially, A is real.

Note also that since tanh(κ + βct) > 0, the new formula for A subtracts 1 from the tanh instead of subtracting the tanh from 1.

Wilhelmy's 1850 Paper on the Dissociation of Cane Sugar

  Wilhelmy published a paper, Ueber das Gesetz, nach welchem die Einwirkung der Säuren auf den Rohrzucker stattfindet (On the law, which takes place after the action of acids on cane sugar) in Annalen der Physik und Chemie in 1850 in which he gives the amount of dissociation of cane sugar (sucrose) into its simpler components (glucose and fructose) as a function of time. One can observe the change in the optical rotation of the sugar solution and determine the rate at which the amount of sugar decreases by means of an inversion coefficient, μ. Wilhelmy uses the optical rotation angle corresponding to the amount of sugar present to also specify its quantity, Z, present in the solution. When a portion, X, of the amount of sugar in solution is inverted it no longer rotates the light passing through the polarimeter X degrees to the right but instead rotates it μX degrees to the left. Wilhelmy explains as follows,

 "The herein stated task will now be to determine, depending on whether and in what way M depends on the different physical conditions of the process, i.e., whether and in what manner M is a function of time, the amount of sugar, the amount of acid, the amount of Resolution agent, the quality of the acid, the temperature and air pressure. Answering these questions should be attempted in the following order.
  The above however is yet an important point to discuss for the calculation of the experiments. The action of the acid on the sugar not only removes a portion of the dextrorotatory sugar but transforms it into opposite rotation. What therefore was originally the rotation =Z°, is, after a rotating X° amount of sugar, converted to =(Z - X - μX)°. To find out that given by observation quantity X itself, one needs to know μ. Hr. Biot has introduced for μ, i.e., for the amount of reverse rotation to the left, produced by a quantity of sugar which had a clockwise rotation of 1°, the name of the inversion coefficient. X then also gives Z = Z₀ - X. The determination of the inversion coefficients has caused some initial difficulties for me. In addition Hr. Biot has been wavering in his statements about the same. In his work on the sugar content of corn ¹) he gives for μ
     for hydrochloric acid  = 0.38
     for sulfuric acid         = 0.3867,
It can be said however that he himself noticed that he received very different results, which he thinks he can ascribe to different purities of the sugar. In a later essay ²) he then gives:
     for sulfuric acid         μ = 0.417
     for nitric acid             μ = 0.394
     for hydrochloric acid  μ = 0.38.
  Thus it is initially to be noted that Hr. Biot, so far as I know, did not state for which temperature these coefficients apply; but since the rotatory power of galactose is dependent on the temperature, as well as the value of μ must be also differ according to the temperature of the reading. If one therefore wishes to know the value of μ for each temperature, one must first determine the law according to which the rotatory power of the galactose is dependent on the temperature."

¹) Comptes rend. 15, 529.
²) Comptes rend. 17, 757.

  The results of Wilhelmy's observations are presented in Table III of his paper. The column headings are time, T, optical rotation, Drehung, the log of the relative change in the amount of sugar, Z, and the temperature. There is a footnote at the bottom of the page which gives a formula for computing the amount of sucrose, Z, present given the rotation angle, D. One can derive the formula as follows,

If we plot Wilhelmy's data we get a reasonably good fit to a decaying exponential function.

Supplemental (Nov 1): The value of the inversion coefficient used for the experimental data is μ = 0.3966. Wilhelmy refers to Biot's corn sugar in German as Schleimzucker which is literally translated as "slime sugar" but in the German Wikipedia it is another name given for galactose. Wilhelmy's definition of the inversion coefficient appears to be based on the relation Z₀ - D = (1 + μ)(Z₀ - Z) where the quantity or concentration of Z is given in units of 1° rotation. For problems in measuring the optical rotation associated with dissociation of complex molecules see the Wikipedia article on mutarotation.