Appendix D — Numerical Methods

Reaction Engineering Basics assumes that most equations will be solved and parameters will be estimated numerically. This appendix briefly describes how numerical methods for accomplishing these tasks work, the input typically required by the numerical methods, and the quantities they return.

D.1 Parameter Estimation

Parameter estimation involves using statistical methods and experimental data to find the best values for unknown quantities appearing in a mathematical model. This is sometimes called fitting the model to the data. In this context, the unknown quantities are referred to as model parameters.

In each experiment, the values of one or more adjusted experimental input variables are set by the person doing the experiments, and then the value of an experimental response is measured. Reaction Engineering Basics only considers parameter estimation using data where the adjusted experimental input variables are the same in every experiment, where only one reaponse is measured in each experiment, and where the measured experimental response variable is the same in every experiment. The values of the experimental input variables and the response variable will change from one experiment to the next, but the variables are the same.

The difference between a measured experimental response and the corresponding model-predicted response is called the experiment residual. Generally, to find the best values for the parameters, the experiment residuals are squared to prevent positive and negative residuals from offsetting each other. The best values for the parameters are then taken to be those that minimize the sum of the squares of all of the experiment residuals. Thus, this is sometimes called least-squares fitting or least-squares parameter estimation.

D.1.1 Linear Model Parameter Estimation

A model where the response is linearly dependent upon the adjusted experimental inputs is called a linear model. The linear models encountered in Reaction Engineering Basics always have either the form shown in Equation D.1 or Equation D.2. In both linear models \(y\) is the response and \(x\) is the adjusted input. Equation D.1 has two parameters: the slope, \(m\), and the intercept, \(b\), while Equation D.2 has only one parameter, the slope.

\[ y=mx+b \tag{D.1}\]

\[ y=mx \tag{D.2}\]

The sum of the squares of the experiment residuals, \(\Psi\), is given in Equation D.3, where \(y_{expt,i}\) is the measured experimental response in experiment \(i\), and \(y_{model,i}\) is the model-predicted response for experiment \(i\). Substituting Equation D.1 results in Equation D.4 where \(x_i\) is the value of the experimental input for experiment \(i\).

\[ \Psi = \sum_i \left( y_{expt,i} - y_{model,i} \right)^2 \tag{D.3}\]

\[ \Psi = \sum_i \left( y_{expt,i} - mx_i - b \right)^2 \tag{D.4}\]

Equations for calculating the best values for \(m\) and \(b\) can be generated by taking the derivatives of \(\Psi\) in Equation D.4 with respect to \(m\) and \(b\) and setting them equal to zero. Analogous approach can be used to generate an expression for the best value for \(m\) in Equation D.2.

Assuming that the experiment residuals conform to a Normal distribution, statistical methods can be used to calculate the upper and lower limits of the 95% confidence intervals for \(m\) and \(b\) and the coefficient of determination, \(R^2\), among other statistical metrics.

D.1.1.1 Linear Model Assessment

The objective of parameter estimation is twofold: to find the best values for the model parameters and to assess the accuracy of the resulting model. Several criteria can be used to assess accuracy. First, if the model has perfect predictive power, the coefficient of determination will equal 1.0, and if it has absolutely no predictive power, \(R^2\) will equal zero. Hence the closer \(R^2\) is to zero, the greater the accuracy of the model.

Second, if lower and upper limits of the 95% confidence intervals for \(m\) and \(b\) are very close to the values of \(m\) and \(b\), that also indicates that the model is accurate. In other words, when the model is accurate, the 95% confidence interval for each parameter will be very small compared to the estimated value of that parameter value.

A graphical assessment is also possible when both the experimental \(y_{expt,i}\) vs. \(x_i\) data points and the \(y_{model}\) vs. \(x\) line are plotted on the same axes. When the model is accurate, the data points will all lie close to the line there will be no systematic trends in deviations.

Ultimately, assessing the accuracy of a fitted model involves a judgement call.

D.1.1.2 Linear Model Parameter Estimation Functions

Many mathematics software packages include a function that performs linear model parameter estimation. Linear least squares analysis can also be performed using many spreadsheet programs or handheld calculators. The required input consists of the set of experimental inputs, \(\underline{x}\), and the corresponding set of measured experimental responses, \(\underline{y}_{expt}\).

Reaction Engineering Basics assumes that whatever tool is used returns the best estimates for the slope, \(m\), and intercept, \(b\), the 95% confidence intervals for \(m\) and \(b\), and the coefficient of determination, \(R^2\). At the minimum, every linear least-squares fitting tool will return the best estimates for \(m\) and \(b\). Often they can be configured to additionally return specific statistics. The more common statistics include the coefficient of determination, \(R^2\), the correlation coefficient, \(r\), the standard errors for each estimated parameter, and the 95% confidence interval for each estimated parameter. If the tool can’t be configured to return a specific statistic, that statistic can usually be computed using the results that are returned.

D.1.2 Non-Linear Model Parameter Estimation

to be added

D.2 Solving Algebraic-Transcendental Equations

Algebraic-transcendental equations (ATEs) are one of the equation types that must be solved when performing reaction engineering tasks. (The acronym ATE is used in this book, but it is not generally used in other books or online.) In the context of Reaction Engineering Basics, ATEs are easily identifiable because they are the only equation type that does not contain derivatives. Slightly more specifically, a set of ATEs is a group of 1 or more mathematical equations that may involve or contain math operations (addition, subtraction, multiplication, and division). Terms in ATEs can be exponents or bases that are raised to some power, and they can appear in transcendental functions. Exponential functions are the most common transcendental functions appearing in the ATEs found in Reaction Engineering Basics. They arise any time an equation includes a rate coefficient that displays Arrhenius temperature dependence, Equation 3.21.

A set of \(N\) ATEs can be solved to find the values of the \(N\) variables appearing in them. A set of ATEs can have more than one solution. An ATE solver is a computer function that solves a set of ATEs numerically. This sub-appendix describes the input required by ATE solvers, how they work, the results they return, and potential issues that can arise when using them.

D.2.1 ATE Solver Input

Typically, an ATE solver must be provided with two inputs. The first is an initial guess for the solution, and the second is a residuals function. Depending upon the mathematics software being used, the solver may accept other input. The order and formatting of the arguments also vary between mathematics software packages.

D.2.1.1 The Initial Guess for the ATE Variables

As just noted, when an ATE solver is first called, it must be provided with an initial guess for a solution. Providing an acceptable guess is not usually difficult. In reaction engineering analyses, the ATE variables have a physical meaning. That, together with a qualitative understanding of reactor performance, provides guidance for making an initial guess. In this book, examples that involve solution of sets of ATEs will have “Click Here to See What an Expert Might be Thinking at this Point” callouts that discuss the choice of the initial guess for that set of equations. Similarly, some of the ATEs that are solved in Reaction Engineering Basics have multiple solutions. The initial guess determines which solution the solver will find, and a callout will describe how the initial guesses leading to each solution were generated.

D.2.1.2 The Residuals Function

In preparation for numerically solving a set of ATEs, each equation must be rewritten as a residual expression. Doing so is trivially simple. If there is a zero on either side of the equals sign, no rearrangement is necessary. If not, everything on one side of the equals sign should be subtracted from both sides of the equation. This will result in an equation with a zero on one side of the equals sign. The nonzero side of that equation is called a residual. A residual expression is created by choosing a variable to represent the residual and setting it equal to the nonzero side of the equation. In Reaction Engineering Basics, \(\epsilon\) is usually used to represent residuals.

As an example, consider the typical ATE mole balance shown in Equation D.5. It does not have a zero on either side of the equals sign, so the right side of the equation, \(\dot{n}_{A,1} + rV\), is subtracted from both sides, leading to Equation D.6. Letting \(\epsilon\) represent the residual, the corresponding residual expression, Equation D.7, is written by setting \(\epsilon\) equal to the non-zero side of Equation D.6

\[ \dot{V}C_{A,0} = \dot{n}_{A,1} + rV \tag{D.5}\]

\[ \dot{V}C_{A,0} - \dot{n}_{A,1} - rV =0 \tag{D.6}\]

\[ \epsilon = \dot{V}C_{A,0} - \dot{n}_{A,1} - rV \tag{D.7}\]

If a set of \(N\) ATEs is being solved, each of them must be converted into a residual expressions. In general each residual can be a function of all of the ATE variables. Substitution of a solution will cause all of the residuals to evaluate to zero.

The purpose of the residuals function is to evaluate the residuals, given a guess for the solution. The engineer solving the ATEs must write the residuals function, but because it will be called by the ATE solver, the arguments passed to it and the values it returns are specified by the ATE solver. The arguments to the residuals function are guesses for the ATE variables, and it always returns the corresponding values of the ATE residuals.

If additional input is needed to evaluate the residuals, that input must be made available to the residuals function by some means other than passing them as arguments. The documentation for the solver may suggest a preferred way to do this, but in Reaction Engineering Basics quantities that cannot be passed to the resdiuals function as arguments are assumed to be made available to it as global variables.

D.2.2 How ATE Solvers Work

Effectively, the ATE solver uses an iterative, trial and error process to find a solution.

The ATE residuals corresponding to the initial guess are calculated and retained as the best solution so far.
The ATE solver generates a new guess and calculates the corresponding ATE residuals.
The residuals corresponding to the new guess and the best solution so far are compared.
Whichever guess gave residuals that are closer to zero is retained as the best solution so far
Steps 2 through 4 are repeated until the solver determines that either
1. the best solution so far is acceptably close to the exact solution, or
2. it cannot find a solution that is acceptably close to the exact solution.

For details on how the solver generates guesses, how it determines which set of residuals is closer to zero, and the criteria for the determinations in step 5, one can consult the documentation for the solver being used and reference works on numerical methods.

D.2.3 Convergence and Solver Issues

Ideally, the ATE residuals should get closer and closer to zero with each iteration. This is called convergence to a solution. The solution returned by the solver will not be exact, but it will be “very close” to the exact solution. Put differently, when a converged solution is found, the difference between the solution returned by the solver and the exact solution is negligible. If the solver is unable to converge to the point where the residuals are “very close” to zero, it will print an error message and/or return a flag variable indicating that it did not converge and why.

ATE solvers can fail to converge for different reasons, and for each reason, there are corrective actions that can be taken. Nonetheless, if the initial guess provided to the solver is “close enough” to a solution, the solver will converge. When the ATE solver fails to converge it is a good idea to check the inaccurate solution it returned. If an ATE variable in that solution is absurdly large or small, it often indicates that these is an error in the computer code or that the units being used in the calculations are not consistent. Otherwise, the inaccurate solution, together with a qualitative understanding of reactor performance, may suggest a way to revise the initial guess so that the solver does converge.

Finally, as already noted, sets of ATEs can have more than one solution. If a set of ATEs has multiple solutions, most ATE solvers will only find one of them. To find other solutions, the solver should be called again using a different initial guess.

D.2.4 ATE Solver Return Values

Even if it fails to converge, an ATE solver typically will return two things. The first is the best solution it was able to find and the second is some indicator of whether it converged. The latter might be a boolean variable or an integer where the value of the integer indicates convergence or the reason why the solver stopped iterating. Some solvers also return a text string that provides a brief explanation of the reason why it stopped trying to improve the solution.

D.3 Solving Initial-Value Ordinary Differential Equations