Suppose that in a spreadsheet the following formula is entered into cell B1:

Since cell reference B1 is in the formula of the same cell, this is called a “circular reference.” Normally, Excel will complain about this. However, if the Enable Iterative Calculations (see File/Options/Formulas) box is checked, Excel will immediately do the following:

1. It will automatically use an initial guess of zero (there is no control over this).

2. It will iterate up to 100 times or until successive values are within 0.001 (these values can be changed).

3. The final result will be the answer (hopefully).

In-cell iteration is not to be encouraged when robustness is a goal. It often fails to find a solution.

EXERCISES

Exercise 1.1: For the following functions, graph the function and then use Goal Seek to find the root(s).

a. f(x) = x − x^1/3 − 2

b. f(x) = xtanx − 1

c. f(x) = x⁴ − e^x + 1

d. f(x) = x²e^x − 1

Exercise 1.2: Find the roots of the functions given in Exercise 1.1 using the bisection method. Use the graph of each function to choose points that bracket the root of interest.

Exercise 1.3: Set up a spreadsheet that implements the secant method and then solve each of the problems from Exercise 1.1. Use the graph of each function to select an initial guess. Recall the iteration formula for the secant method:

Hint: Set up the first row of the spreadsheet for your problem with headings such as the following:

Put the formula for the function under the headings f(xk-1) and f(xk). In the cell under xk+1, put the secant method iteration formula. In the second row, replace the previous xk-1 with xk and then xk with xk+1. Now copy the two formulas down one row. At this point, one iteration of the secant method is displayed. To see more iterations, just copy the second row down for as many iterations as desired. If too many iterations are copied and the function difference (the denominator of the iteration formula) becomes exactly zero, a “divide by zero” error will appear.

Exercise 1.4: Use Goal Seek to find root(s) of the following functions. Plot the functions first to obtain an approximation of a desired root.

a. f(x) = x³ – 17x + 12 = 0

b. J₁(x) = 0 J₁ is the Bessel function of the first kind of order 1. It can be computed in Excel using the BESSELJ()function. This function has an infinite number of roots; find the root between 2 and 5.

c. Solve for the molar volume of a gas at 400 K and 1200 kPa using the van der Waals equation of state. The critical temperature and pressure are 500 K and 80 atm, respectively. Use the ideal gas solution for your initial guess.

d. Solve the Colebrook equation for the Darcy friction factor, f, for a Reynolds number (N_Re) of 10⁵ and a roughness factor, ε/D, of 10⁻⁴ (this equation holds for Reynolds numbers > 4000):

e. Repeat part d using the same roughness factor, but for a range of Reynolds numbers from 5000 to 30,000 (pick a reasonable increment). Plot the results (friction factor versus Reynolds number). Automate Goal Seek for this.

Exercise 1.5: Use the secant method as described in Exercise 1.3 to find the root(s) of the functions given in Exercise 1.4. Carefully choose the two initial guesses so that the function values have opposite signs. The roots found may or may not correspond to those found using Goal Seek in Exercise 1.3—it depends on the initial guesses.

Exercise 1.6: Solve the problems of Exercise 1.1 using the Newton method.

Exercise 1.7: Repeat Exercise 1.4 parts a and b using the Newton method. The derivative of J₁(x) is given by

Exercise 1.8: Repeat Exercise 1.4 using the bisection method.

Exercise 1.9: An additional method for solving single nonlinear equation is the “Regula–Falsi” or method of “false position.” Like the bisection method, the false position method starts with two points a₀ and b₀ such that f(a₀) and f(b₀) are of opposite signs, which implies that the function f has a root in the interval [a₀, b₀]. The Regula–Falsi method proceeds by producing a sequence of shrinking intervals [a_k, b_k] that always contain a root of f.

    At iteration number k, the value

is computed. As explained below, c_k is the root of the secant line through (a_k, f(a_k)) and (b_k, f(b_k)). If f(a_k) and f(c_k) have the same sign, then set a_k+1 = c_k and b_k+1 = b_k; otherwise set a_k+1 = a_k and b_k+1 = c_k. This process is repeated until the root is approximated sufficiently well.

The above formula is also used in the secant method, but the secant method always retains the last two computed points, while the false position method retains two points that bracket a root. On the other hand, the only difference between the false position method and the bisection method is that the latter uses c_k = (a_k + b_k)/2.

FIGURE 1.7 Regula–Falsi method.

A schematic description of the Regula–Falsi method is shown in Figure 1.7.

Repeat either Exercise 1.1 or 1.4 using the method of false position.

Exercise 1.10: Refer to the diagram and data of Example 1.7. The dew point of a vapor is the temperature at a given pressure at which the first drop of liquid is formed. Thus, at the dew point, the ratio of vapor to feed (V/F) is essentially one. Conversely, the bubble point of a liquid is the temperature at a given pressure at which the first bubble of vapor is formed and at this condition V/F = 0.

Assume, as in Example 1.7, an ideal mixture (where Raoult’s law applies). Also, assume the same mixture and data as shown in Figure 1.4.

a. At the dew point, V/F = 1 [there is an infinitesimal amount of liquid, but it does exist and has a mole fraction for each component (x_i); however, since there is only an infinitesimal amount of liquid, V = F and y_i = z_i]. Applying the fact that the sum of x_i = 1, the following equation results:

Since k_j are functions of temperature, this nonlinear equation can be solved for the temperature (dew point).

b. At the bubble point, V/F = 0 (L/F = 1) and x_i = z_i, and if that fact that the sum of y_i = 1 is applied, the following equation results:

This nonlinear equation can be solved for the temperature (bubble point).

The specific assignment is as follows.

1. Prepare an Excel spreadsheet that calculates for a mixture the bubble point at pressures of 15, 16, …, 25 atm and produces a plot of the bubble point versus pressure with suitable annotations and title. Use Goal Seek to solve the single nonlinear equation at each pressure. Data for the system are given in Figure 1.4.

If it is (somehow) known that a root lies in the interval [a, b], then by simply halving the interval in which the root lies, the interval can be reduced to an acceptable level. This idea is at the heart of the bisection method as shown in Figure 1.2.

The restriction is that f(a) and f(b) must have opposite signs—one of them must be positive, the other negative (it does not matter which). Then, because f is assumed to be continuous, it must be a zero somewhere in [a, b]. Let c be the midpoint of [a, b]. Either c is the root, or the root lies in [a, c] or in [a, b]. If f(c) is close enough to zero (see below regarding tolerance), then the root has been found. Otherwise, one pair of [f(a),f(c)]or [f(c),f(b)] has opposite signs. Keep the half-interval with opposite signs and discard the other. Repeat the process until either (1) f, evaluated at the midpoint of the interval, is sufficiently small or (2) the interval has been shrunk to a suitably small value.

FIGURE 1.2 Bisection: c = (a + b)/2.

The term “sufficiently small” is usually tested using a “tolerance,” which represents a number small enough to be considered zero based on the application. This can be stated more formally as follows:

If |f(x)|<tolerance then (x) is sufficiently small.

Example 1.3: Bisection Applied to the van der Waals EOS

Recall the van der Waals EOS in the form

The table below shows the progression of applying the bisection method. It uses the Excel IF function to choose whether f(a) or f(b) is replaced by f(c) (and correspondingly whether a or b is to be replaced by c). The syntax for the IF function is as follows (see the explanation of the Excel spreadsheet below for a full explanation of how the IF function is used):

IF(test, True Value, False Value)

In this example (again for ammonia at 10 atm and 250°C), the function in the cell below the label V1 (corresponding to the point a in the nomenclature of the figure) is = IF(F2<0, E2, A2), where cell F3 holds f(Vnew) [corresponding to f(c)], E2 contains Vnew (corresponding to c), and A2 holds V1. So, if f(Vnew) is negative (test is TRUE), V1 (corresponding to a) gets the Vnew value; otherwise, it retains the old V1. Likewise, the formula in the cell below the V2 label is = IF(F2<0,B2,E2), which replaces V2 (corresponding to b) with the value Vnew if the test (F2 < 0) is FALSE. After the first row of formulas has been entered, these cells are copied down to repeat the iterations. Iterations should be repeated until f(Vnew) is sufficiently small (in absolute value). The initial guesses for V are 3 and 5, respectively, which give opposite signs for the function.

Did You Know?: That there are two ways to copy cell contents down. One way is to select the cells to be copied and then grab the small box in the lower right corner of the selection. Pulling down will copy the cells. The second method is to select the cells to be copied and while holding down the mouse button pull down as far as desired and finally hitting CTRL/d (hold down the CTRL key and hit the d key).

1.2.4 NEWTON’S METHOD

The next algorithm to be considered is the Newton’s method for finding roots. Newton’s method does not always converge. But, when it does converge, it usually does so very rapidly (at least once it is “close enough” to the root). Newton’s method also has the advantage of not requiring a bracketing interval with positive and negative values. So Newton’s method allows the solution of equations such as

whereas bisection does not. This parabola just touches the zero axis and has no negative values.

The basic idea of the Newton algorithm is this: given an initial guess, call it x¹ to a root of f(x) = 0, a refined guess, x², is computed based on the x-intercept of the line tangent to f(x) at x¹. That is, consider the equation of the line tangent to f(x) at x¹ (this is just the Taylor series expansion of the function ignoring all but linear terms):

This is the point-slope equation of a line, where x¹ is the base point and f΄(x¹) is the slope [derivative of f(x) evaluated at x¹]. Solving Equation 1.8 for x at which f(x) = 0 gives

or more generally,

The value of x^k+1 is the new guess at the root. The process is repeated, computing successively x², x³, x⁴,… until an x^K is found at which

where tol is a prescribed tolerance.

Example 1.4: Newton’s Method Applied to the van der Waals Equation

Newton’s method is now applied for ammonia at 10 atm and 250°C with an initial guess for V = 1 L/gmol. To begin using Newton’s method, the derivative of Equation 1.2 is required and is as follows:

Here is an Excel spreadsheet that solves this problem:

The root is found in about 3 iterations. Note the rapid rate of convergence.

An examination of Equation 1.9 for the new iterate x^k reveals a potential for failure of Newton’s method, namely,

This may lead to a wildly divergent iterative process. There are other possible reasons why this method might not converge. In general, Newton’s method is prone to failure, but when it does work, it converges rapidly. Good initial guesses are the key to success.

Excel has very good plotting capabilities. Unfortunately, it is not possible in Excel to simply give a command such as plot(f(x)). It is necessary to produce a list or table of x and f(x) values and to graph the resulting data. This is best illustrated by an example.

Example 1.1: Plotting the Equation

A table of data and an Excel graph of the van der Waals equation for ammonia at 250°C and 10 atm are shown in Figure 1.1.

From the graph, it is easy to see that there is one real root between 0 and 5 (the other roots are complex conjugates). Remember that the van der Walls equation of state is an attempt to model the nonideality of the gas. Since the given temperature and pressure are not severe, it is expected that the calculated molar volume from the equation of state should not be greatly different from that predicted by the ideal gas law. From the ideal gas law, the molar volume is 4.29 L/gmol, which is very close to the root shown in the graph.

1.2.2 FIXED-POINT ITERATION (DIRECT SUBSTITUTION)

To apply this method, the equation must be cast into the form

x = g(x)

or, more generally,

FIGURE 1.1 Roots of the van der Walls equation for ammonia at 250°C and 10 atm.

where k is an iteration counter.

Example 1.2: Direct Substitution

Note that this can always be accomplished by adding x to each side of f(x) = 0, if necessary. The van der Walls equation can be cast into the following form:

The pertinent data for ammonia are shown in Figure 1.1. If a value for V = 4 is guessed (based on the graph of Figure 1.1) and is used on the right-hand side, a new (and hopefully better) value is calculated from Equation 1.5. The iterations produce the following sequence:

4.00000 4.21854 4.22946 4.22996 4.22998.

The solution is 4.22998 L/gmol. If more significant digits are required, then more iterations can be carried out.

It should be noted that direct substitution can be a divergent process. That is, the successively calculated values actually get worse rather than closer to the correct value. Without proof, the following statements apply:

Let g΄ be the first derivative of the function g in Equation 1.4. Then,

• If |g΄| < 1, the error will decrease with each iteration.

• If |g΄| > 1, the error grows at each iteration.

• If g΄ > 0, the error will have the same sign at each iteration.

• If g΄ < 0, the error will alternate signs at each iteration.

Clearly, the equation should be arranged so that the magnitude of g΄ is less than 1. This might take some experimentation. Often, a form of g is tried, and if the process does not converge, then other forms are attempted.

In Chapter 4, methods for approximating the derivative of a function using finite differences are presented. The secant method uses the idea of finite differences to approximate the derivative in the Newton method formula. Starting with two initial guesses x⁰ and x¹, which need not bracket the root of interest, the approximation to f΄(x) can be written as follows:

Or, in general, after k steps or iterations,

Substituting this approximation into the Newton formula (Equation 1.9), the following iteration formula results for the secant method:

Example 1.5: The Secant Method for van der Waals Equation

The following shows the van der Waals example solved using the secant method. The two initial guesses for V are 1.50 and 1.51—the second was chosen arbitrarily close to the first one. Note that the rate of convergence is about the same as that of Newton’s method. Beware that the secant method is subject to the same potential shortcomings of Newton’s method.

Many engineering problems require the solution of a single nonlinear equation. Such an equation can always be cast into the form

The objective of this chapter is to study methods and learn of Excel^® tools for finding the root(s) of a nonlinear equation, that is, for finding x such that f(x) = 0.

Simple algebra provides the root for a linear equation. However, for more complex (nonlinear or transcendental) equations, it is often the case that no analytical solution is available, or is difficult to obtain, so that numerical methods must be used.

An example of a nonlinear equation is the van der Waals equation of state, which is given by

where

with subscript c referring to the “critical” values of temperature and pressure for the gas. Note that, if a and b are zero, this reduces to the ideal gas equation of state.

A typical problem is to find the molar volume, V, given the temperature and pressure (and the type of gas). While it is possible to find analytic solutions to the van der Waals equation of state for V (this is simply a cubic equation), a numerical solution is often preferred. The equation of state in the form f(x) = 0 is obtained by simple rearrangement:

Note that this particular form f(x) = 0 is not unique, and other algebraic rearrangements are possible.

1.2 ALGORITHMS FOR SOLVING f(x) = 0

Clearly there is a need to find good methods for determining roots. Four such methods are now presented (though many more have been developed): fixed-point iteration (direct substitution), bisection, Newton’s method, and the secant method. These methods are iterative. That is, given a guess of a root, or the interval in which a root lies, the algorithm refines that guess repeatedly, obtaining (hopefully) better and better guesses, until a value “close enough” to the true root is found. Also, graphing the equation can add insight into the roots of interest and can provide good initial estimates of the roots for the iterative algorithms. It is recommended to always prepare a graph of the function to give insight into possible solutions.

1.2.1 PLOTTING THE EQUATION

Excel has very good plotting capabilities. Unfortunately, it is not possible in Excel to simply give a command such as plot(f(x)). It is necessary to produce a list or table of x and f(x) values and to graph the resulting data. This is best illustrated by an example.

Example 1.1: Plotting the Equation

A table of data and an Excel graph of the van der Waals equation for ammonia at 250°C and 10 atm are shown in Figure 1.1.

From the graph, it is easy to see that there is one real root between 0 and 5 (the other roots are complex conjugates). Remember that the van der Walls equation of state is an attempt to model the nonideality of the gas. Since the given temperature and pressure are not severe, it is expected that the calculated molar volume from the equation of state should not be greatly different from that predicted by the ideal gas law. From the ideal gas law, the molar volume is 4.29 L/gmol, which is very close to the root shown in the graph.

1.2.2 FIXED-POINT ITERATION (DIRECT SUBSTITUTION)

To apply this method, the equation must be cast into the form

x = g(x)

or, more generally,

FIGURE 1.1 Roots of the van der Walls equation for ammonia at 250°C and 10 atm.

where k is an iteration counter.

Example 1.2: Direct Substitution

Note that this can always be accomplished by adding x to each side of f(x) = 0, if necessary. The van der Walls equation can be cast into the following form:

The pertinent data for ammonia are shown in Figure 1.1. If a value for V = 4 is guessed (based on the graph of Figure 1.1) and is used on the right-hand side, a new (and hopefully better) value is calculated from Equation 1.5. The iterations produce the following sequence:

4.00000 4.21854 4.22946 4.22996 4.22998.

The solution is 4.22998 L/gmol. If more significant digits are required, then more iterations can be carried out.

It should be noted that direct substitution can be a divergent process. That is, the successively calculated values actually get worse rather than closer to the correct value. Without proof, the following statements apply:

Let g΄ be the first derivative of the function g in Equation 1.4. Then,

• If |g΄| < 1, the error will decrease with each iteration.

• If |g΄| > 1, the error grows at each iteration.

• If g΄ > 0, the error will have the same sign at each iteration.

• If g΄ < 0, the error will alternate signs at each iteration.

Clearly, the equation should be arranged so that the magnitude of g΄ is less than 1. This might take some experimentation. Often, a form of g is tried, and if the process does not converge, then other forms are attempted.

1.2.3 BISECTION

If it is (somehow) known that a root lies in the interval [a, b], then by simply halving the interval in which the root lies, the interval can be reduced to an acceptable level. This idea is at the heart of the bisection method as shown in Figure 1.2.

The restriction is that f(a) and f(b) must have opposite signs—one of them must be positive, the other negative (it does not matter which). Then, because f is assumed to be continuous, it must be a zero somewhere in [a, b]. Let c be the midpoint of [a, b]. Either c is the root, or the root lies in [a, c] or in [a, b]. If f(c) is close enough to zero (see below regarding tolerance), then the root has been found. Otherwise, one pair of [f(a),f(c)]or [f(c),f(b)] has opposite signs. Keep the half-interval with opposite signs and discard the other. Repeat the process until either (1) f, evaluated at the midpoint of the interval, is sufficiently small or (2) the interval has been shrunk to a suitably small value.

FIGURE 1.2 Bisection: c = (a + b)/2.

The term “sufficiently small” is usually tested using a “tolerance,” which represents a number small enough to be considered zero based on the application. This can be stated more formally as follows:

If |f(x)|<tolerance then (x) is sufficiently small.

Example 1.3: Bisection Applied to the van der Waals EOS

Recall the van der Waals EOS in the form

The table below shows the progression of applying the bisection method. It uses the Excel IF function to choose whether f(a) or f(b) is replaced by f(c) (and correspondingly whether a or b is to be replaced by c). The syntax for the IF function is as follows (see the explanation of the Excel spreadsheet below for a full explanation of how the IF function is used):

IF(test, True Value, False Value)

In this example (again for ammonia at 10 atm and 250°C), the function in the cell below the label V1 (corresponding to the point a in the nomenclature of the figure) is = IF(F2<0, E2, A2), where cell F3 holds f(Vnew) [corresponding to f(c)], E2 contains Vnew (corresponding to c), and A2 holds V1. So, if f(Vnew) is negative (test is TRUE), V1 (corresponding to a) gets the Vnew value; otherwise, it retains the old V1. Likewise, the formula in the cell below the V2 label is = IF(F2<0,B2,E2), which replaces V2 (corresponding to b) with the value Vnew if the test (F2 < 0) is FALSE. After the first row of formulas has been entered, these cells are copied down to repeat the iterations. Iterations should be repeated until f(Vnew) is sufficiently small (in absolute value). The initial guesses for V are 3 and 5, respectively, which give opposite signs for the function.

Did You Know?: That there are two ways to copy cell contents down. One way is to select the cells to be copied and then grab the small box in the lower right corner of the selection. Pulling down will copy the cells. The second method is to select the cells to be copied and while holding down the mouse button pull down as far as desired and finally hitting CTRL/d (hold down the CTRL key and hit the d key).

1.2.4 NEWTON’S METHOD

The next algorithm to be considered is the Newton’s method for finding roots. Newton’s method does not always converge. But, when it does converge, it usually does so very rapidly (at least once it is “close enough” to the root). Newton’s method also has the advantage of not requiring a bracketing interval with positive and negative values. So Newton’s method allows the solution of equations such as

whereas bisection does not. This parabola just touches the zero axis and has no negative values.

The basic idea of the Newton algorithm is this: given an initial guess, call it x¹ to a root of f(x) = 0, a refined guess, x², is computed based on the x-intercept of the line tangent to f(x) at x¹. That is, consider the equation of the line tangent to f(x) at x¹ (this is just the Taylor series expansion of the function ignoring all but linear terms):

This is the point-slope equation of a line, where x¹ is the base point and f΄(x¹) is the slope [derivative of f(x) evaluated at x¹]. Solving Equation 1.8 for x at which f(x) = 0 gives

or more generally,

The value of x^k+1 is the new guess at the root. The process is repeated, computing successively x², x³, x⁴,… until an x^K is found at which

where tol is a prescribed tolerance.

Example 1.4: Newton’s Method Applied to the van der Waals Equation

Newton’s method is now applied for ammonia at 10 atm and 250°C with an initial guess for V = 1 L/gmol. To begin using Newton’s method, the derivative of Equation 1.2 is required and is as follows:

Here is an Excel spreadsheet that solves this problem:

The root is found in about 3 iterations. Note the rapid rate of convergence.

An examination of Equation 1.9 for the new iterate x^k reveals a potential for failure of Newton’s method, namely,

This may lead to a wildly divergent iterative process. There are other possible reasons why this method might not converge. In general, Newton’s method is prone to failure, but when it does work, it converges rapidly. Good initial guesses are the key to success.