ROOT-FINDING !
IMPLEMENTING MS-EXCEL
IN
Chapter Goals
Learn about implementing root finding methods in MS-EXCEL Learn about GOAL SEEK function in MS-EXCEL
+D=FJAH
INTRODUCTION
In the last chapter we studied how to set up mathematical equations in MS-EXCEL and make straightforward algebraic calculations. While studying environmental systems, it also becomes necessary to obtain roots of the equations and solve a system of linear and nonlinear equations. Algorithms for performing these computations can be coded fairly easily in most cases. MS-EXCEL has a suite of tools for facilitating algebraic calculations. We will look into some of these tools in this lesson and identify any potential pitfalls associated with using them.
Finding Roots for a Single Equation
Implementing Direct Methods
There are several approaches to finding the roots of a single equation. In some instances, an analytical formula is available and can be directly implemented in MS-EXCEL. In particular, for quadratic equations there are standard formulas available: ax2 + bx + c = 0 has the following roots
-b
b 2 - 4ac 2a
(13-1)
MS-EXCEL provides SQRT(X) function that can be used to implement the above formula.
148 Applied Environmental Systems Modeling Example 13-1 An industrial wastewater containing 50 mol/L (M) dimethyl phthalate undergoes reaction in a well-mixed 400 L pilot-reactor with a flowrate of 2000 L/day and operated under steady-state conditions. Determine the steady-state outflow concentration of the compound. Assume, the second order reaction rate constant for the reaction is equal to 6.9 10-2 M-1s-1 at 30oC. The second-order rate constant is converted to Lmol-1day-1 units as follows to be consistent with the units of flowrate (L/day) l
86400 s L L L = 6.9 10 -2 = 5961.6 mol - d mol - d mol - s 1d
(13-2)
The steady-state mass-balance expression for the pollutant under steadystate flow of wastewater (Q) can be written as:
V dC = 0 = QCin QCi lC2V dt
(13-3)
The steady-state concentration C can be solved using Equation 13-1 as: C =
Q
Q 2 + 4 (QC in lV ) 2 lV
(13-4)
Cells containing values of the variables can be named and the following formulas can bed coded in two cells (say A13 and A14) of MS-EXCEL as: A13 = (Q + SQRT(Q^2 + 4(Q*Cin*lambda*V)))/(-2*lambda*V) A14 = (Q - SQRT(Q^2 + 4(Q*Cin*lambda*V)))/(-2*lambda*V) (13-5) (13-6)
Coding the formulas in the spreadsheet the roots were determined to be 0.20 mol/L and -0.21 mol/L. Only the positive root is valid as the concentration cannot be negative. The spreadsheet is depicted in Figure 13-1.
Implementing Bracketing Schemes Bisection Method
When direct formulas are not available to solve for the roots of the equations, indirect methods have to be resorted to. There are several numerical methods for solving the roots of the equations and these can be classified into two categories The closed or bracketing schemes and the open non-bracketing schemes. As their name suggests, the bracketing schemes search for a root within a closed interval. At the root, the value of the function is equal to zero. As such, the value of the function will switch sign on either side of the root. This characteristic is used to identify an approximate range or bracket in which the root is searched. To implement the bisection method, we start by guessing two roots for the function F(x) that bracket the actual root. Let these two roots be xL and xU. The
Implementing Root-Finding in MS-EXCEL
149
Figure 13-1
Spreadsheet setup for finding the roots of the quadratic equations.
values xL and xU are chosen such that the product of their functions is less than zero: F(xL) F(xu) < 0 The first guess for the root is computed as: xr = (13-6)
(x L + x U ) 2
(13-7)
If the value of the function is equal to zero, the program stops. However, if such isnt the case a new bracket is defined using xr and either of xL or xU. Again, the sign change of F(x) is used to identify the new bracket as follows: IF F(xL) F(xr) < 0 IF F(xr) F(xU) < 0 The guess is between xL and xr Then the guess is between xr and xU (13-8) (13-9)
Using the new guesses, a new bracket is obtained and the process is repeated until the root is obtained with desired accuracy. MS-EXCEL provides an IF Statement that can be used to code the conditional expressions in Equations 13-8 and 13-9. Example 13-2 The manufacturer of a proprietary pre-fabricated water disinfection treatment unit claims 99% reductions in coliform concentration after passing secondary treated wastewater through their unit. The removal reaction kinetics was empirically found to vary as the square root of coliform concentration in this unit with a reaction rate constant of 7500 (#/L)0.5d-1. If
150 Applied Environmental Systems Modeling a small treatment plant has purchased a 500 L unit with a maximal allowable flowrate of 1200 L/d, what is the maximum allowable inflow concentration of the coliform into the disinfection unit, if it functions as a completely mixed flow reactor and as per the manufacturers specification. Assume steady-state operating conditions. The steady-state mass-balance equation for coliform transport in the disinfection unit can be written as 0 = Q Cin QC lV C Per manufacturers claims
C in - C = 0.99 C = 0.01 Cin C in
(13-10)
(13-11)
Hence, Equation 13-10 can be written as: 0 = QCin 0.01Q Cin lV 0.01 C in Re-arranging, we obtain 0 = 0.99 Q Cin 0.1 lV C in (13.13) (13-12)
Equation 13-13 can be further re-arranged to obtain the maximal inflow concentration as:
C in =
0.101 lV 0.101 lV Cin = Q Q
FG H
IJ K
2
(13-14)
Substituting the values given we obtain a maximal inflow concentration of 99639.07 #/L (or 99.63 #/100 mL). Now let us solve Equation 13-13 using the Bisection method and compare how this method fares against the direct solution. A spreadsheet set-up for implementing bisection method is shown in Figure 13-2. As a first step let us plot the function over a range of values of Cin which is assumed to be between 106 to zero. As can be seen from Figure 13-2, there is a sign change for the function when Cin is between 105 and 104. This range can be further refined by inserting intermediate values, but is used directly starting points for the Bisection method and provides a test for the method to find correct solution from a rather large bracket. The mid-point, namely 5500 is used as the first guess. Clearly, this is not the solution so we need to refine the bracket. The following formulas are entered in MS-EXCEL in various cells: Cell G7: Cell H7: =(F7+E7)/2 =0.99*Qin*E7-lambd*SQRT(0.01*E7)*Vol (13-15) (13-16)
Implementing Root-Finding in MS-EXCEL
151
Figure 13-2 root finding.
A spreadsheet method for implementing the bisection method for
Cell I7: Cell J7:
=0.99*Qin*F7-lambd*SQRT(0.01*F7)*Vol =0.99*Qin*G7-lambd*SQRT(0.01*G7)*Vol
(13-17) (13-18)
The calculations in cells H7, I7 and I7 evaluate the function (the RHS of Equation 13-13) for lower and upper guess and the mid-point value computed in Cell G7. The values of these functions are used to obtain new guesses in Cells E8 and F8 as follows: Cell E8: =IF(H7*J7<0,E7,G7) Cell F8: =IF(I7*J7<0,F7,G7) (13-19) (13-20)
The Formula in Cell E8 states that IF (F(xL) x F(xr) < 0 ) then the solution is more closer to the lower limit xL which is used the lower limit in the next iteration as well, else the mid-point is used as the lower-limit. Similarly, the formula in Cell F8 states that IF (F(xL) x F(xr) < 0 ) then the solution is more closer to the upper-limit which is used as the upper limit in the next iteration as well, else the solution is close to the computed mid-point which is used as the upper limit.
152 Applied Environmental Systems Modeling If the guess for the root is exact, the true value is equal to zero. The deviation of the function F(xr) from zero is computed in column K and the improvement over previous iteration is computed in Column L to check how the bisection method performs at each iteration. Once the formulas are entered for the first time, they can simply be copied for subsequent iterations by selecting the appropriate cells and using the EDIT FILL DOWN from the menu bar (or pressing CTRL+D). As can be seen from Figure 13-1, the concentration result of 99639.07 #/L (or 99.63 #/100 mL) obtained using the bisection method compares well with the direct solution that we obtained from direct algebraic manipulations. Further more, the solution converges in about 25 iterations even with a very large bracket. However, the method requires two guesses to work, as the root has to be bracketed. One needs to compute the function at several values to identify where the change of sign occurs. Although quite simple, the bisection method may not be the most efficient method in some cases. There are other approaches that work with a single guess and when two guesses are required they need not bracket the root. The Newton-Raphson and Secent method are two such methods which we shall explore next.
Implementing Open Schemes NewtonRaphson/Secant Method
The Newton-Raphson method is based on the Taylor series expansion of a function about a point. If the function F(x) is known at the point x and we intended to compute the value of F(x) at some other point defined as x + Dx, then the Taylor series expansion can be written as: F(x + Dx) = F(x) + DxF(x) +
( Dx ) 2 F(x) + 2!
(13-21)
Where F(x) is the first-derivative of F(x) with respect to x and F(x) the second derivative. Neglecting (Dx)2 and higher terms we can simplify the Taylor series expansion as: F(x + Dx) = F(x) + DxF(x) (13-22)
Now, if x + Dx is indeed the root of the function F(x) then F(x + Dx) would be equal to zero. If such is the case then Equation 13-22 can be further simplified as: Dx =
F ( x) F (x )
(13-23)
If the value of x after the nth iteration is xn (i.e, x = xn) and value of x at the n+1 iteration is x + Dx (i.e., x + Dx = xn+1) then we can write Equation 13-23 as:
Implementing Root-Finding in MS-EXCEL
153
Dx = xn+1 xn =
F ( x) F (x n ) xn+1 = xn F (x ) F (x n )
(13-24)
The computational procedure starts with an initial guess for xn which is used to compute the new value xn+1 using Equation 13-24. The process is repeated until the F(xn+1) becomes reasonably close to zero. The Newton-Raphson method is very sensitive to the initial guess that is supplied to start the iterations. The method may not converge to a solution if the guess is far away from the actual root. Also the derivative F(x) must exist at x and we need to compute it before-hand. Example 13-3 Repeat Example 13-2 using the Newton Raphson Method. Solve, starting with an initial guess of 106 and 105 #/L for Cin. The required function is given by Equation 13-12 which is re-written below: F(x) = 0.99Q Cin 0.1 lV
C implement in
(13-25)
To Newton-Raphson method we need the derivative of the function F(x) which is given as (Equation 13-26): F(x) = 0.99Q
0.1 lV 2 C in
(13-26)
The spreadsheet for implementing Newton-Raphson Method is shown in Figure 13-3.
Figure 13-3
Implementation of Newton-raphson method on a spreadsheet.
The required data is entered in Cells H12, H13, H14 and H15 for Q, l, volume and initial guess respectively. The initial guess is copied into Cell F5 and the following formulas are filled in:
154 Applied Environmental Systems Modeling Cell G5: Cell H5: Cell F6: Cell I5: Cell J6: =0.99*$B$12*F5-0.1*$B$13*$B$14*SQRT(F5) =0.99*$B$12-(0.1*$B$13*$B$14/(2*SQRT(F5))) =F5-(G5/H5) =0-F6 =I5-I6 (13-27) (13-28) (13-29) (13-30) (13-31)
These formulas in these cells are copied down for subsequent iterations by selecting the cells and pressing CTRL+D (EDIT FILL DOWN).
Figure 13-4 guesses.
Convergence of Newton-raphson method with different starting
The number of iterations required to converge to a solution with two different starting conditions one order of magnitude apart are presented in Figure 13-4. As can be seen, the convergence is pretty quick even when the starting condition is one-order of magnitude higher. Clearly, the convergence is further improved when the initial guess is close to the true value. The bisection method took about 20 iterations to find the root against 7 and 3 by Newton-Raphson with different starting conditions. While Newton-Raphson is computationally efficient, there is no guarantee that a solution will be found. This is particularly true if the initial guess is not a good one and the series is divergent.
Implementing Root-Finding in MS-EXCEL
155
The Secant method is a modification of the Newton-Raphson method and is better suited when the derivative is difficult to get. The derivative is computed using the first-order finite-difference approximation as: F(xn)
F ( x n ) - F( x n - 1 ) x n - x n -1
(13-32)
The algorithm starts with initial guesses xo and x1 which are used to compute the derivative F(x1) using Equation 13-32 which in turn is used to compute x2 using Equation 13-24. Note that while the secant method requires two initial guesses, they need not bracket the root as is required with the bisection method. Hence, the method can be modified to require a single initial guess by somewhat arbitrarily assuming that xo equals x1 0.001 (or some other small number). While this assumption introduces arbitrariness, it also reduces the data requirements. This trade-off must be evaluated prior to using the specified heuristic. Example 13-3 Repeat Example 13-2 using the Secant method. Solve, starting with initial guess for Cin,1 as 106 and 105 #/L. Assume Cin,o = Cin,1 0.1 for simplicity. To implement the secant method we need not know the derivative, as it is computed using the Equation 13-32. A spreadsheet to solve the problem using the secant method is presented in Figure 13-5. The required inputs Q, l, V and Cin,1 are specified in Cells B12:B15 respectively and the initial guess is copied into cell F6. The following formulas are entered in the cells to carryout necessary computations: Cell F5: Cell G5: Cell H6: Cell I6: Cell J7: =F6-0.1 (compute the required second guess xo) =0.99*$B$12*F5-0.1*$B$13*$B$14*SQRT(F5) =(G6-G5)/(F6-F5) =0-G6 =I6-I7 (13-33) (13-34) (13-35) (13-36) (13-37)
The results presented in Figure 13-5 indicate that the efficiency of the secant method is good for this problem and comparable to that of the Newton Raphson procedure. The method converges to a solution in five to seven iterations. Again, the selection of the starting condition has a lot to do with the number of iterations required and the ability of the method to converge to a solution. The use of a simple heuristic wherein the guess for xo is calculated from x1 worked well in this case and helped reduce the need to provide two guesses.
156 Applied Environmental Systems Modeling
Figure 13-5
Spreadsheet implementation of secant method.
Implementing Root-Finding in MS-EXCEL
157
Root Finding using GOAL SEEK
The above examples illustrate that root finding algorithms can be easily coded into spreadsheets. While the examples provide useful illustrations of various numerical methods for finding roots, it is often not necessary to implement them as MS-EXCEL provides an in-built tool called GOAL SEEK that can be accessed from the TOOLS menu (TOOLS GOAL SEEK). The GOAL SEEK is an implementation of the modified secant method and requires a single initial guess (xo) for the root. The cell containing the guess must not have a formula in it. The function F(x) is coded into a cell and makes reference to the cell containing the initial guess (xo). The GOAL SEEK menu (see Figure 13-6) asks for the value the cell containing F(x) is to take and which cell to change. For finding root of F(x) we say change F(x) to a value of zero by changing the cell containing xo. Clearly, the cell containing F(x) must contain a formula with x for GOAL SEEK to work.
Figure 13-6
Implementing root finding using GOAL SEEK functionality.
158 Applied Environmental Systems Modeling Example 13-4 Repeat Example 13.3 using GOAL SEEK. Solve, starting with initial guess for Cin,1 as 106 and 105 #/L. The spreadsheet for implementing GOAL SEEK root finding procedure is shown in Figure 13-6 and consists of four steps. In step A - the required inputs Q, l, V and Cin,1 are specified in Cells B12:B15 respectively and the formula F(x) (Equation 13-25) is entered in Cell B16 as: Cell B16: =0.99*B12*B15-0.1*B13*B14*SQRT(B15) (13-35)
In Step B the GOAL SEEK is selected from TOOLS Menu (TOOLS GOAL SEEK); In Step C the cell B16 (containing the formula) is specified as the cell whose value has to be set to zero by changing cell B15 (containing the value for initial guess for Cin). Step D depicts the GOAL SEEK STATUS message box indicating that a solution has been found. Also note that the value in the cell B15 containing the initial guess has been modified from 1E+06 #/L to 9.96E+04 #/L which is indeed the final steady-state value. Remember that GOAL SEEK employs a modified secant method and has some restrictions on the number of iterations it can perform. Hence, sometimes it may not return a solution, especially when the initial guess is too far from the original solution. Also, as the derivatives are being computed numerically, there is a possibility of converging to a wrong solution as well. ...