I n d . Eng: Chem. Res. 1988,27, 1320-1329
1320
Mochimaru, Y. J . Non-Newtonian Fluid Mech. 1983,12, 135-152. Nada, I.; Yamada, Y.; Nagasawa, M. J. Phys. Chem. 1968, 72, 2890-2898. Ottinger, H. C., private communication, 1986. Ottinger, H. C. J . Non-Newtonian Fluid Mech. 1987, 26, 201-246. Peterlin, A. J. Polym. Sci., Polym. Lett. 1966, 4B, 287-291. Tanner, R. I. Trans. SOC.Rheol. 1975, 19, 37-65. Thomas, R. H.; Walters, K. Q.J. Mech. Appl. Math. 1964,17,39-53. van Kampen, N. G. Stochastic Processes in Physics and Chemistry; North-Holland: Amsterdam, 1981; Section VIII.6. van Wiechen, P. H.; Booij, H.C. J . Eng. Math. 1971,5, 89-98.
Walters, K.; Savins, J. G. Trans. SOC.Rheol. 1965, 9, 407-416. Walters, K.; Waters, N. D. Rheol. Acta 1964, 3, 312-315. Warner, H. R., Jr. Ind. Eng. Chem. Fundam. 1972, 11, 379-387. Wedgewood, L. E.; Bird, R. B. In Integration of Fundamental Polymer Science and Technology; Kleintjens, L. A., Lemstra, P. J., Eds.; Elsevier,,Applied Science: London, 1986; pp 337-345. Wedgewood, L. E.; Ottinger, H. C. J . Non-Nontonian Fluid Mech. 1988,27, 245-264. Received for review August 3, 1987 Accepted February 22, 1988
Computing All Real Solutions to Systems of Nonlinear Equations with a Global Fixed-point Homotopy Minoru Kuno* and J. D. Seader Department of Chemical Engineering, University of Utah, Salt Lake City, Utah 84112
A global fixed-point homotopy has not found wide application in chemical engineering for solving systems of nonlinear equations by continuation. Rather, the Newton and problem-dependent homotopies have been favored. However, it is conjectured here that the fixed-point homotopy would be expected to always place all real roots of the system on a global homotopy path because the path is forced t o begin from a single starting point. If so, all roots could be computed by following the single path with a suitable continuation method. This would be particularly desirable when the number of real roots to a system cannot be predetermined and one wishes t o compute all solutions. In this study, it was found that such a path does exist provided that a criterion is used for selecting a starting point which minimizes the number of real roots of the global fixed-point homotopy function a t a n infinite value of the homotopy parameter. Several examples, including a n adiabatic reactor with multiple solutions are presented t o illustrate the application of the criterion. While methods have been devised recently t o find all roots of polynomial equations, this is the first method for systematically locating all roots t o general systems of nonlinear equations. The necessity of solving systems of nonlinear equations often arises in simulating and designing a chemical plant or optimizing a process. Newton's method or Powell's method are commonly applied to solve such systems. However, as described by Seader (1985), both methods have well-known disadvantages such as the requirement that the starting point (initial guess of the solution) must be in the neighborhood of a root. That is, these methods are only locally convergent. In addition, both methods are designed to locate, a t best, just one root even though multiple solutions may exist. Continuation methods are rapidly being applied to numerous problems in engineering analysis as summarized in a recent review by Seydel and Hlavacek (1987). A solution method for nonlinear equations that is globally convergent from any starting point is the homotopy continuation method. This method, as described, e.g., by Garcia and Zangwill (1981)) does not solve the nonlinear equations directly, but gradually reaches a root by beginning from a starting point which satisfies a second, simpler system of equations. Both systems of equations are embedded into a so-called homotopy function. Thus, if f ( ~ ) is the system of nonlinear equations to be solved and g ( x ) is a second simpler system of the same number of equations, the homotopy function might be constructed as H(x,t) = tf(x)
+ (1- t ) g ( x ) = 0
(1)
*Author to whom a11 correspondence should be addressed.
Present address: Daitec Co., Ltd., 4-85 Chikara-machi Higashi-ku Nagoya-shi, Japan. Minoru Kuno was on leave from Daitec Co., Ltd., during the course of the study reported here. 0888-5885/88/2627-1320$01.50/0
where t is a scalar homotopy parameter which is gradually varied from 0 to 1as the path is tracked from the starting point to a solution. The homotopy function can be constructed in accordance with the characteristics of each system of nonlinear equations, but this effort can be time consuming. Alternatively, canonical homotopy functions can be used, two of which, as described by Garcia and Zangwill (1981) and Wayburn and Seader (1987), are the Newton homotopy and the fixed-point homotopy, where g ( x ) = f ( x ) - f ( x o ) and x - xo, respectively. Until now, most studies have focused on the use of homotopy continuation to obtain only one root even though other roots may exist. Studies to obtain all roots are recent and have been restricted largely to systems of polynomial equations, as discussed by Morgan (1987))where multiple starting points are used to obtain multiple solutions. Allgower and Georg (1980, 1983) discuss methods for approximating additional, but not necessarily all, roots of general systems of nonlinear equations. The aim of this paper is to present a method for determining a starting point from which one can find all the roots along one branch of the path when using a global fixed-point homotopy. The homotopy parameter is permitted to take on values greater than one, but the homotopy path remains bounded for all values of the homotopy parameter within the region where the real path connects all the zeros of f ( x ) . The real path may form an isola, which is a closed loop that does not contain the starting point. In general, the global Netwon homotopy and global fixed-point homotopy encounter difficulties when an at0 1988 American Chemical Society
Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 1321 tempt is made to locate all roots of f ( x ) = 0. For instance, when we use the global Newton homotopy to obtain all real roots from one starting point, we cannot guarantee that all solutions will be obtained. Part of the cause for failure may lie in the fact that multiple values of the starting point, xo, may satisfy g ( x ) = f ( x ) - f(xo) = 0, as discussed by Lin et al. (1987). Multiple paths may be required to obtain all solutions. From the structure of the global fixed-point homotopy, only one starting point can exist for a selected value of xo because H ( x , t ) = x - xo = 0 at the starting point. Thus, it would seem that all of the roots should lie on a single branched path provided that the complex domain is included. However, it is often found that unbounded multiple branches exist that are only connected by mapping + m into -. Such unbounded multiple branches can be avoided and all real roots can be found in the real domain on a single branch of the path for which t 2 0 by use of the criterion presented here for selecting the starting point. Some applications of the homotopy continuation method for the chemical industry have been discussed by Bhargava and Hlavacek (1984), Byrne and Baird (1985), Chavez et al. (1986), Hlavacek and Rompay (1985),Lin et al. (19871, Salgovic et al. (1981), Vickery and Taylor (1986), Wayburn and Seader (1987), and Kovach and Seider (1987). The homotopy function used in most of these references is the Netwon homotopy. The fixed-point homotopy has been used hardly a t all. The fixed-point homotopy differs significantly from the Newton homotopy in that the vector homotopy function has only one starting point once the set xo is selected. It is conjectured here that if xo is properly selected, according to a criterion presented later, a single branch of a bounded real homotopy path will contain all real roots to f ( x ) . As will be shown, a preferred starting point is one which minimizes the number of real roots of the homotopy function at t = +a. Various examples, including a chemical reactor with multiple steady states, are given to illustrate applications of the criterion. The method is suitable for general systems of nonlinear equations that may contain transcendental terms. Theoretical Development Global Fixed-point Homotopy. Consider a system of nonlinear equations
-
f(x) = 0
(2)
where f = Rn R", the continuously differentiable function vector; and x E R", variables. Two canonical homotopy functions which are used to solve this system are Newton homotopy
H(x,t) = f(x)
- (1- t)f(xO) = 0
t E [O, 11
(3)
fixed-point homotopy
H ( x , t ) = (1- t ) ( x - XO)
-
+ tf(x) = 0
t E [O, 11 (4)
where H = Rn+l R", the homotopy function vector, t = homotopy parameter, and xo = starting point vector. At the starting point, t = 0, eq 3 becomes f ( x ) = f ( x o ) and eq 4 becomes x = xo, respectively. A t t = 1, both eq 3 and 4 become H ( x , l ) = f ( x ) = 0 and the original system of equations is satisfied. The solution of H ( x , t ) = 0 for t E [0, 11 is obtained by numerical calculations, usually by either of the two general methods, one discrete and the other continuous.
In the former (e.g., Shacham (1986)), a sequence of values of t is selected or generated for solving H(x,t) = 0: H(x,t,) = 0 0 = to
< tl < tz < ... < tl-l < t, < ... < t , t, =
=1
+ At
(5) where At is a suitably small increment, so as to obtain convergence for each incremental move along the path. Typically, at each increment, Newton's method would be used to solve H(x,t) = 0, employing as a starting guess the value of x that was previously found to satisfy H(x,t-1) = 0. In the latter (e.g., Klopfenstein (1961)), H ( x , t ) is transformed to ordinary differential equations by differentiating H with respect to some parameter such as the arc length of the path. The ordinary differential equations are then solved as an initial-value problem using a predictor-corrector technique. The latter method is greatly preferred, but the application of either method is normally initiated at t = 0 and terminated when one root is obtained at t = 1. Nonlinear equations, however, have many roots in general, and it is highly desirable in many chemical engineering problems to seek all real solutions. Methods for finding additional roots are discussed by Allgower and Georg (1983) and include use of multiple starting points, deflation, global Newton methods, global homotopy methods, and a d-trick homotopy developed by Allgower and Georg (1980). By far, the most successful of these approaches, which however is limited to systems of polynomial equations and does not apply to general systems of nonlinear equations, is the global homotopy method as discussed, e.g., by Garcia and Zangwill (1979), Chow et al. (1979), Watson (1986), Morgan (19861, Li (1987), and Watson et al. (1987). Because the maximum number of roots for systems of polynomial equations can be predetermined by Bezout's theorem, as discussed by Morgan (1987), the methods described by the above authors mostly utilize a different starting point for each root sought. Thus, multiple homotopy paths are followed, one for each root. A starting point for each desired root is associated with each path, and the root is obtained by following the path from t = Otot=l. Alternatively, especially for general systems of nonlinear equations for which the maximum number of roots cannot be predetermined, it would seem desirable to have a method that would find all roots from a single starting point. For example, the global Newton homotopy function might be applied in the following modification of eq 3: t,-l
H(x,t) = f ( x ) - (1 - t ) f ( x o ) = 0 t ER (6) where values of t are not now confined to the interval [0, 11. Garcia and Gould (1980) used this homotopy to determine multiple roots and compared the method to Newton's method. Allgower and Georg (1980,1983) developed the d-trick method which uses eq 6 and seeks additional roots after the first root is found. Wayburn and Seader (1987) and Lin et al. (1987) found all real roots from one starting point by using eq 6 for several examples. However, other examples can lead to multiple paths, which require multiple starting points of unknown number. A second alternative is the global fixed-point homotopy. The equation is as follows:
H(x,t) = (1- t ) ( x - no) + t f ( x ) = 0
t ER
(7)
1322 Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988
':$
I-
I
-
2-
x (a) X a t t=O i s - 2
0
-.
-
(d)
-2
x
1:-
a t t=O i s 1.5
-5
I
-
-I -7
-1
-9 -IC
1
-
t
f
-2
-I
~
*
,
2
/I' (bl X a t t=O i s 0
X (C)
X a t t=O i s 1
-1 -1
(f) X a t t=O i s 4
-I -I7 -7 7
-I
-9
-
t
Figure 1. Effect of starting point on homotopy path for f ( x ) = 2 - 32 + 2 = 0.
where H i s global in the sense that t is not confined to the interval [0, 11, but can be allowed to take on negative and/or any positive values. Equation 7 is satisfied at the starting point, t = 0, by only one value of xo. This is not generally so for eq 6. The question to be asked is what is the nature of the homotopy path as a function of the seemingly arbitrary choice of xo. This nature is crucial to the development of the method presented here. To illustrate, consider the very simple, single nonlinear equation f(x) =
x2
- 3x
+2 =0
(8)
which has obvious roots of 1 and 2. If this function is embedded into eq 7, the global fixed-point homotopy, we obtain after simplification
H(x,t)=
t(X2
- 4x
+2+
xO)
+ (X - x0) = 0
tER
(9)
Equation 9 can be plotted in terms of x as a function of t for a particular value of xo, as in shown in Figure 1for six different values of xo. The curves in the figure are termed homotopy paths. The filled-in circles at t = 0 and t = 1 represent the starting point xo and roots of f ( x ) = 0, respectively. The crosses (+) on the branches of the paths are computed points, which were obtained by the path-following predictor-corrector continuation algorithm of Kubicek (1976). Although more efficient continuation codes are available [e.g., HOMPACK of Watson et al. (1987) and PITCON of Rheinboldt (1986)], Kubicek's code was more readily adapted to our purposes. As shown in Figure 1,the nature of the homotopy path (locus of points which satisfy eq 9) for this homotopy function depends
strikingly on the assumed value for xo ( t = 0). When xo = -2 (Figure l a ) or 0 (Figure lb), three unbounded real branches exist. In both cases, the two roots (at t = 1)are located on one of the branches, but the roots cannot be reached from the starting point (at t = 0) without making computations in complex arithmetic on the interval in t of [0.0670, 0.93301 for x o = -2 and [0.1464,0.85361 for xo = 0. For values of xo less than -2, the nature of the path is the same. In the limit, for very large negative values of xo, the complex region is almost the entire interval [0, 11 of t. For xo = 1 (which is one of the two roots), Figure IC shows that the homotopy path consists of four unbounded real branches with a bifurcation point, where paths intersect, a t t = 0.5. Both roots cannot be reached from the starting point without locating the bifurcation point and tracking two branches of the path to the right of the bifurcation point. The complex region disappears when xo is increased to 1.0 (one of the two roots) and does not reappear until xo is increased to a value greater than 2.0 (the other root). For x o = 1.5 (which is intermediate to the two roots), Figure Id shows that three unbounded real branches exist. Only the root at x = 2 can be reached from the starting point. The same is almost the case for x o = 2 (one of the two roots), as shown in Figure le. For xo = 2, the two rightmost branches are actually connected at t = +a,as are the two leftmost branches at t = -a. Thus, only one unbounded real branch actually exists for xo = 2. For values of xo larger than 2, as typified by the curve shown in Figure If for xo = 4, a single real branch is ob-
Ind. Eng. Chem. Res., Vol. 27, No. 7 , 1988 1323 Table I branch 1 2 3 1
beginning x t -2 0 +m
4 0
0 +m --OD
ending x -m
comments
t
0
4
-m
0 -2
+m
0
tained. Of most importance is that both real roots are located on this branch and the branch is not unbounded in t. Both roots can be reached from the starting point without entering a complex arithmetic region and without traversing a large distance in t , as would be required for xo = 2. Thus, from the six cases shown in Figure 1, it appears that the proper selection of xo (in this case, x o > 2) can greatly facilitate the computation of all real roots from just one starting point using a global fixed-point homotopy. For other problems involving a single nonlinear equation or a system of nonlinear equations, different types of real branches containing all roots may occur for the proper selection of xo. For the above example, which involves an even number of real roots, the branch beginning from t = 0, and heading initially in the direction of increasing t , never reaches t = +m in the real domain. If the function is continuous and has an odd number of real roots, the branch will terminate at t = +a. In another case, particularly involving systems of equations, all roots may lie on a real isola (closed curve), as will be shown later. Unfortunately, the isola will not contain the starting point, but can be reached from the starting point by making computations in the complex domain until the first root is located. Alternatively, all roots can be found on the isola by starting from one known root. Next, we develop a criterion for the selection of a suitable xo, called x,O and then a method based on an auxiliary function for applying it when the criterion is difficult to compute. The criterion and the auxiliary function apply to all systems of nonlinear equations, including those with transcendental terms, which occur frequently in chemical engineering problems. On the basis of numerous computational experiments, the application of the criterion and the auxiliary function for selecting x o is the global fixedpoint homotopy permits the computation of all real roots from a single starting point. Criterion for Selecting the Starting Point. With a global fixed-point homotopy, it is possible to find all real roots on one finite-length real branch of the path provided that a proper starting point vector is selected. We give the following criterion for selecting x: from the infinite number of values of xo: x,"
E (xo; xo E R", min,o N(xo))
(10)
where N(xo) = the number of real roots of the equation f(x) - x
2
contains both roots
+ xo = 0
(11)
This criterion is based on the observation, as shown in Figure 1 and verified in numerous other examples not included here, that a homotopy path containing all real roots can always be formed from a global fixed-point homotopy when x," is selected as the starting point. However, if a proper x o vector is not selected, it may be necessary, in order to obtain all roots, to (1) resume a path that appears to terminate a t t = fa or x = f m by reentering the path a t the opposite extreme of x or t (that is, map + m into -a) or (2) traverse an infinite distance between successive roots. For example, in Figure l a for xo = -2, a closed real path is formed by the three branches in Table
xo
1
0
-1
-2
X
Figure 2. Application of the criterion t o f ( x ) = x 2 - 3 x
+ 2 = 0.
I beginning from the starting point and initially increasing the homotopy parameter t , with the third branch containing both roots. For xo = 2, the closed path contains just one branch, which, however, is unbounded in t and must be connected a t t = 0 by exiting a t x = -m and reentering a t x = +m. The two roots are separated by an infinite distance in t. For xo = 4, or any value of xo greater than 2, the real path now becomes bounded in t , but remains unbounded in x . However, both roots are in relatively close proximity and are quickly located from the starting point. This value of xo is a desirable selection. Of particular importance is the observation, from the above discussion for the example of Figure 1,that the most desirable path for locating all roots is one where the homotopy function has no roots at t = +a. To determine this, the global fixed-point homotopy of eq 7 is divided by t and allowed to go to the limit as t +m to give f ( x ) x + xo = 0, which is eq 11. The relationship between x o and x in eq 11for the example of eq 8 is shown for the real domain in Figure 2. For this case, eq 11 becomes x 2 - 4x + 2 + xo = 0, and the curve in Figure 2 is the locus of points which satisfy this equation. It is clear that, as shown in Figure 2, when xo > 2, the equations cannot be satisfied in the real domain. That is, the number of real roots, N(xo) in eq 10, is at a minimum of zero for xo > 2. For xo < 2, two real roots exist. For example, for xo = 1,the roots are 1 and 3 as seen in Figure 2 . Thus, x," > 2 = xo should be selected and is consistent with the above discussion accompanying Figure 1. Ideally, we then select x," from all xo such that eq 11, which accompanies the criterion of eq 10, has no real roots. Then all real roots of the equations being solved, f ( x ) = 0, will be forced to lie on a bounded portion of the real branch of the homotopy path starting from xo in the direction of increasing t. However, if more than one equation is involved, a complex domain may need to be traversed to connect the starting point to the real roots, which will lie on an isola. This technique always works when f ( x ) has an even number of real roots. However, when f ( x ) has an odd number of real roots, a modification is needed; viz., we select x o such that eq 11 has one real root. In general, the criterion is given by eq 10 in conjunction with eq 11. That is, we seek an xo that minimizes the number of real roots of the homotopy function at t = +a. For complicated problems, it would appear that finding the number of real roots, N(xo),of eq 11 for different values of xo might be more difficult than finding all solutions of the original problem, f ( x ) = 0. In practice, this is not always the case.
-
1324 Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 5
0
-5
-*?-
0
__
.- 7-1
-
~ -
T
- - ~ r
2
-10
3
0
1
7
X
Figure 3. Application of the criterion to f ( x ) = x 3 - 6 x 2 + llx - 6 = 0.
I
I
4i
Path for Xo
i
0.0
/ -
I
//
1
As a more complicated example, consider the problem of finding all stationary points of Himmelblau's function (Reklaitis et al., 1983). The function is as follows:
i 1
f I'
31
t
Figure 5. Homotopy path for Himmelblau's function.
1
g(x) =
(XI'
+ ~2
- 11)'
+
(XI
+ ~2~
- 7)'
(14)
Accordingly, the stationary points are the roots of the following equations:
4% 4x1
+ 2x1xg - 21x, + x2' - 7 = o
- = f l ( x ) = 2x13
$gg/4x2 = f 2 ( x ) = x 1 2 + 2 x 1 x 2
Path for Xo = 1.5
+ 2xZ3- 13x2- 11 = o (15)
1 i 01
4 I
I,
-1
-5
5
0
'
10
t
Figure 4. Homotopy path for f ( x ) = x 3 - 6x2 + llx - 6 = 0.
xIo = - 2 x I 3 - 2x,x2
One can often quickly find suitable values of x: that satisfy eq 10 and 11 by considering either large positive and/or negative values of xo and examining the magnitudes of the resulting terms of eq 11. An example of an odd number of real roots is the simple cubic equation
f ( x ) = x 3 - 6x2 + l l x - 6 = 0
(12)
which has three roots (all real) at one, two, and three. For this example, eq 11 for selecting x,O becomes x0 = -x3
+ 6x2 - lox + 6
Lin et al. (1987) obtained all of the real roots of eq 15 from one starting point, x: = x: = 0.4, by using the global Newton homotopy. Here we solve eq 15 by the global fixed-point homotopy. For the starting point criterion, eq 11 becomes
(13)
As shown in Figure 3, which is the locus of points satisfying this equation, the number of real roots is a minimum (in this example, one) for xo > 3.0887 and x o < 0.9113. For 0.9113 < x o < 3.0887, Figure 3 shows that three roots of x exist. The successful application of the criterion of eq 10 and 11for this example is shown by the homotopy path of eq 7 in Figure 4, where all three roots (at t = 1)to eq 12 are obtained on one branch for x,O = xo = 0, which satisfies the criterion; and not for xo = 1.5 # x:, which does not satisfy the criterion. In Figure 4, closed circles are starting points ( t = 0) and roots ( t = 1). The dashed path corresponds to the selection of an xo according to the criterion. In this case of an f ( x ) = 0 with an odd number of real roots, the path, for a proper selection of xo = x,O, is not closed and is unbounded in both x o and t. However, all the roots lie on a finite length of that path.
xzo = -2x23 - 2 7 ~ 1 x 2
+ 2 2 2 , - xZ2+ 7 14x2 - XI' -k 11
(16)
where f i ( x ) is associated with xl and fz(x) is associated with x2. This association can be made in an arbitrary manner. Values of xIo and xZowhich satisfy the criterion of eq 10 and 11 are expected to have large absolute values. For instance, if we assign the value of -50 to both x: and x: and solve the global fixed-point homotopy formulation of eq 15, all nine roots are obtained as shown by the path in Figure 5 . Also, because the number of real roots of eq 15 can be shown to be one for these values of xo, one end of the path goes to t = +a. If we try to solve this problem using starting points of even larger absolute value, all of the real roots are still obtained. However, if we use a starting point which does not satisfy the criterion of eq 10 and 11,e.g., xIo = xZo= 0.0, which results in nine real roots of eq 16, only one root is obtained on the path from the starting point before it goes to t = +m. Sometimes the global Newton homotopy can find all real roots from one starting point as in the previous example. However, for other cases, it may fail. For instance, consider the following nonlinear equations:
+ 2 x 1 x 2 - 2 2 x , + xZ2+ 13 = 0 + 2x1x2 + 2xz3 - 14x2 + 9 = 0
f,(x) = 2 x I 3 f 2 ( ~ = )
(17)
The difference between eq 15 and 17 is small. Equation 17 has seven real roots given in Table 11. If we solve this problem by using the global Newton homotopy with the
Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 1325 Table 11. Real Roots of Eauation 17 -4.255 -3.454 -3.369 0.987
-3.845 1.144 2.514 1.801
0.713 0.793 3.252
Embed f(x) into global fixed-point homotopy Eq 7 I
I
0.856 -2.814 -2.712
Path goes t o t =
I I
Apply Criterion o f Eqs 10 and 11 to determine xo
-
Path goes t o x = Minimum Number o f Real Roots of H(x,t) at I = +a
-
t
OD
I /
t
1 (Path goes to t = +-)
Domain Calculations All Real
Roots
Real Domain Calculations All Real Roots
Some Complex Domain Calculations All Real
Roots (C)
Path goes t o x = w i t h isola
2
-
Figure 7. Types of homotopy paths for global fixed-point homotoPY.
1,and t = m have important meanings. The starting point on the t axis is the only point that joins the paths of the domains of positive and negative values of t. Accordingly, the t axis serves as a barrier to other branches of the path, which cannot cross the t axis. All real roots are located at t = 1. If a single branch of the path contains all roots and the starting point, that branch cannot go to t = + m until all roots are located. The number of branches of the path that go to t = + m is equal to the number of real roots of eq 11because those branches continuously approach the roots of eq 11. Therefore, in order to obtain all of the real roots on one branch of the path from one starting point, it is necessary that one selects a value of xo that minimizes the number of the real roots of eq 11. However, this is not a sufficient condition. Even if we choose such a starting point, we have no guarantee that a continuous real path will contain the minimum starting point and all of the real roots unless one real root of eq 11 exists. If a minimum of zero roots of eq 11exists, an isola may be present which is only connected to the starting point through a complex domain portion of the path. Use of an Auxiliary Function To Satisfy the Criterion. The criterion given by eq 10 and 11is very simple, but not always easy to apply. In a large complicated system of nonlinear equations, we cannot, in general, easily find a starting point which satisfies the criterion. Accordingly, an auxiliary function was developed as an alternative to the use of the criterion. The auxiliary function is always easy to apply, but it increases by one the number of equations that must be solved simultaneously. The auxiliary function, defined as follows, is added to the system of equations to be solved:
1326 Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988
where
Intermediate
A(xl,...,x,) = (B(xl,...,x , ) ) ~+~ c
flow
I .o
(19)
/
or
A(xl,...,x,) =
dB(xl,-,xn)
+c
(20) L
where x,+~ = auxiliary variable, B = scalar function (e.g., B(xl,...,x,) = xl), c = constant (e.g., O), d = positive constant (e.g., the exponential 2.718...), and m = positive integer (e.g., 1). Equation 19 is preferred when equations are dominated by polynomial terms, while eq 20 may be useful when exponential terms are present. Of importance is that both functions are designed to assume only positive values. Typical quantities for A are xi2 and exp(xi). The constant c is included, if desired, to assist in shifting A to smaller or larger values. If the original system of equations does not contain, in one of the terms of one of the equations, a factor which depends upon a variable and can only assume positive values, such a factor must be introduced in order to use the auxiliary function. This can be done in a variety of ways. For example, for the function f ( x ) = 2x1x2 + x3x4 = 0, incorporate exle-x* into the first term to obtain f ( x ) = 2exle-x1xlx2 + x3x4. Then define the auxiliary function, using eq 20, as x5 - ex' = 0. The original function then becomes f ( x ) = 2x5e-x1x1x2+ x 3 x 4 = 0. When the auxiliary function is introduced for a general system of equations, the criterion, eq 11,now becomes as follows: f i ( X 1 , ...,X",X,+l) f k..., ,X " , X " O * )
- Xi + X1° = 0 - X2 + X20 = 0
...,Xn,Xn+J - X n + zoo = 0 -A(X~,...,X,,) + x,+: = 0
fn(X1,
(21)
When we substitute eq 19 or eq 20 into the last equation of eq 21, the following equations are obtained: x,+10
-c
= {B(x~,...,xn)]zm
(22)
or x,+10
-c =
dB(xlv-,xn)
(23)
Because the right-hand side of eq 22 or 23 is always positive, the simplest starting condition for x,+: that satisfies the criterion is Xn+lo
0.36, we obtained all roots from the starting point by tracking a real path. However, when xlo was < 0.36, an isola with all real roots was formed, which could only be reached by tracking a complex path between the starting point and the first root.
Applications Multiple Steady States of a CSTR Reactor. The classical example of important multiple solutions to nonlinear equations in chemical engineering is multiple steady states for the operation of a continuous-flow adiabatic stirred-tank reactor (CSTR). The multiple solutions are established by the simultaneous solution of the mass balance, based on a kinetic expression for an exothermic first-order reaction, and the energy balance, based on an enthalpy of reaction. As discussed by Smith (1981) and shown schematically in Figure 8, a low or high feed concentration can result in just one steady-state solution, while an intermediate feed concentration can lead to three possible steady-state solutions. Alternatively, the three cases can be achieved by changing the feed temperature instead of feed concentration. Consider the example cited by Smith, which leads to the following two balances in terms of the fractional conversion, X, and the reactor effluent temperature, T , with the reactor feed temperature, T,, as a parameter: energy balance
X
1 = -150 (T-
T,)
mass balance x=-
R l+R
where (25)
(26)
R = 1.34 X
lo9 exp[-62800/(8.314T)]
Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 1327 .J5
1
'0°
600
500
-
7 400
X1
.65
04
6 5
0
99
1
The solution to this problem is generally achieved by plotting the two equations as X vs T and locating the points of intersection or by substituting eq 28 into eq 29 to eliminate X and solving the resulting transcendental equation
(T - Tf)-
1.34 X
1
lo9 exp(-7554/T)
+ 1.34 X lo9 exp(-7554/T) (150) = 0 (30)
by some numerical method for the roots. Let us apply the algorithm of Figure 6 to eq 30 to obtain all real roots from one starting point. When eq 30 is embedded into the global fixed-point, we obtain H(T,t) = (1 - t ) ( T - T o )+ f ( T ) = 0
(31)
where To is the starting point and not the reactor inlet temperature. Applying the criterion of eq 10 and 11,which gives for eq 11, gives
T - Tf-
lo9 exp(-7554/T) 1 + 1.34 X lo9 exp(-7554/T)
(150)1.34 X
or
To-Tf=
1
lo9 exp(-7554/T)
+ 1.34 X lo9 exp(-7554/T)
(33)
But the right-hand side is always positive. Therefore, the criterion for the starting point that will minimize the number of real roots of eq 33 is that To- Tf < 0 or To < TP
Figure 10. Example of kola formation.
sider the following example taken from Allgower and Georg (1983): fl(xl, x 2 ) = ( x , - x22)(x1- sin x 2 ) = O
f&
= (cos x2 - x J ( x 2 - cos X I ) = 0
(34) where x1 and x 2 are in radians. An infinite number of real roots exist, but we are only interested in those roots (known to be four) which lie in the regions 0 6 x1 6 1.0 and 0 6 x 2 6 1.0. To accommodate the inequality constraints, we can introduce a slack variable, x 3 , with the equation x2)
f3(x1, x 2 , x3)
For Tf = 298 K, homotopy paths are shown in Figure 9 for values of To = 100, 200, and 250, each of which satisfies the criterion. The three roots obtained are 300.2, 348.0, and 445.5. For Tf = 258 K, a single root of 260 is obtained, while for Tf = 318 K, a single root of 486 is obtained. Although this problem can be solved quite readily by graphical or other methods as mentioned above, more realistic cases involving more complex kinetics, multiple reactions, multiple reactors, and temperature-dependent properties cannot be solved by classical methods. The method described here can solve these more complex cases. Constrained Systems. The method described here for finding all real roots can also be applied to the determination of all real roots within a constrained region. Con-
= x 2 ( x 2 - 1) + x3' = 0
(35)
Because x: is positive, x 2 cannot take on values larger than one or less than zero. Then, from eq 34, x1 cannot take on values larger than one or less than zero. Equations 34 and 35 are embedded in a global fixedpoint homotopy and the criterion is applied first to eq 35 to give x30
-T+To=O (32)
(150)1.34 X
1 12
t
Figure 9. Multiple steady-state solution for a CSTR reactor by global fixed-point homotopy continuation method.
f(T) =
101
1
t
=
x3
- X ' ( X 2 - 1) - x32
(36)
which can be rearranged to x30
=
y2 - ( x 2 - 72)2 - ( x , - '/2)2
(37)
Equation 37 will have no real roots if x: > lI2. Furthermore, values of xl0 and x20 can then be selected arbitrarily. The global fixed-point homotopies of eq 34 and 35 were solved for the following starting points, with all four roots [(0.7071,0.7854,0.4105);(0.6792,0.8241,0.3807); (0.6416, 0.8011, 0.3992); (0.6946, 0.7679, 0.4222)] being obtained: starting pt for x I o = x: = x: 0.60 0.65 0.70 0.80 0.85 0.90
results all four roots all four roots all four roots all four roots all four roots isola
In the last case, the path, shown as a projection for x1 in Figure 10, had to include complex number computations from t = 0.91 to 0.99 to reach an isola, which contains all four roots. The isola appears to include a bifurcation point at t = 1.0027. However, homotopy paths for the other two variables, x 2 and x 3 , do not intersect a t t = 1.0027. Therefore, this point is not a bifurcation point. Bifurcation points can be located by methods described by Doedel (1981).
1328 Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 y = (1 - 1 . 2 5 ~ ) / ( 1- X )
(A-3)
Conclusions The ability of a global fixed-point homotopy to find all of real roots of systems of nonlinear equations from one starting point is dependent upon the starting point. In this study, a criterion was developed for choosing a correct starting point. The criterion given is not always sufficient for locating all real roots by using continuation solely in the real domain. However, if failure occurs when using real numbers, success can be achieved by computing in the complex domain until the first root is located and then transferring to the real domain to locate all of the remaining roots. The algorithm presented should find wide application to chemical engineering problems where multiple solutions may exist and it is desirable to compute them. Although no mathematical proof of the validity of the criterion for finding all real solutions is presented, the criterion has been applied successfully and withou failure to numerous examples, some of which are discussed above. Additional examples taken from the literature are given in the Appendix.
f ( x ) = (1/63) In x
Nomenclature A = auxiliary function B = function c = constant d = constant f = function vector or mapping g = function vector H = homotopy function N = number of real roots of eq 11 R" = vector space of n-tuples of real numbers t = homotopy parameter ti = ith t defined in eq 5 to = starting point t At = increment oft T = effluent temperature Tf= feed temperature x = vector of independent variables; an independent variable xo = value of x at starting point of homotopy function x," = value of x at correct starting point of homotopy function X = conversion
+ X 2 + X 3 / 2 - Xg/X7 = 0 (A-7) f z ( x ) = x3 + x4 + 2x5 - 2/xI = 0 (A-8) f&) = x1 + X' + x5 - l/x, = 0 (-4-9) f4(x) = -28837x1 - 139009x2 - 78213x3 + 18927x4 -t
Subscripts c = correct i = ith 0 = starting point Superscripts m = dimension of vector 0 = starting point
Appendix The following three examles, taken from the chemical engineering literature, further illustrate the application of the starting point criterion presented above as eq 10 and 11. All examples have multiple roots, two contain logarithmic terms, and one involves a set of seven equations. For each example, the equation(s) is (are) first reduced, transformed, and/or simplified to obtain a more appropriate form for the application of the continuation method. Example A-1, Shacham (1986). The following equation gives the fractional conversion of a species in a chemical reactor: f(x) = x/(l - x ) 5 In [0.4(1 - x)/(0.4 - 0.5~11+ 4.45977 = 0 (A-1) This equation can be transformed to an exponential form f(y) = y - exp(O.8~- 1.691954) = 0 (A-2) where
The starting point criterion is satisfied for yo < 0, since eq 11 would then have no real roots. If yo = -0.5 is selected, real roots are obtained: 0.757 396 and 1.098984. The latter is infeasible because the maximum conversion is 1.0. Example A-2, Paterson (1986). The following equation gives the mole fraction of a species being separated by batch distillation a t infinite reflux:
+ (64/63) In [ l / ( l- x ) ] + In (0.95 - x ) - In (0.9) = 0 (A-4)
This equation can be transformed to an exponential form
fb)= y - [ ( l + 0 . 0 5 ~ ) / 0 . 9 5 ] ~ ~ ( 0=. 90) ~ ~(A-5) where
y = ~ / ( 0 . 9 5- X )
(A-6)
The starting point criterion is satisfied for yo < 0, since eq 11 would then have no real roots. If yo = -1.0 is selected, all real roots are obtained: 0.5 and 0.036 210 08. Example A-3, Carnahan et al. (1969). The following set of equations represents three atom balances, an energy balance, a mole fraction constraint, and two equilibrium relations for the production of synthesis gas in an adiabatic reactor: = X1/2
fi(X)
8427x5-13492 - 10690-x 6 = 0 (A-10) XI f5(~)
=
f6(x)
XI
+ ~2 + ~3 + ~4 +
XI XS
- 1= 0
= P 2 X 1 X , 3 - 17837oX3X5 = 0
f , ( ~= ) XlXB - 2 . 6 0 5 8 ~ 2 ~=4 0
(A-11) (A-12) (A-13)
where P = 20 atm, x1 - x 5 are mole fractions at chemical equilibrium for five species, and xg and x , are mole ratios. Equations A-7 to A-10 can be linearlized by letting x 8 = X g / + and x9 = 1/x7. The resulting equations together with eq A-11, which is already linear, can be commined with eq A-12 and A-13 to reduce the system to two simultaneous nonlinear polynomial equations in x1 and x4: = -0.001x1X~ 0.16517~4'- 0.046354xlX4 0.20707~4- 0.262098~1'+ 0.17967~1+ 0.043528 = 0 (A-14)
fs(X)
and f7(X)
=
O.lx4'
+ 0.0213875~40.086512~1~ - 0.016415~1= 0 (A-15)
- 0.185963x1X4
To constrain the variables, x1 and x4, which represent mole fractions between 0 and 1.0, it is convenient to introduce two additional equations: f&) = q ( x 1 - 1) + XI02 = 0
(A-16)
f&) = x4(x4 - 1) + X I 1 2 = 0
(A-17)
where xl0 and xI1 can be considered slack variables. Equations A-16 and A-17 can only be satisfied for values of x1 and x 2 between the constraints of 0 and 1. As shown in the subsection of this paper entitled Constrained Sys-
Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 1329 tems, the criterion of eq 10 and 11will be satisfied for xloo and xl: > 0.5. The transformed system of equations was solved for four different sets of starting points: xl0 = x40 = x12 = xllo = 0.6,0.7,0.8, and 0.9. In all cases, computations were made in the complex domain because of the existence of an isola similar to that shown in Figure 7c. The following two sets of real roots were located on the isola for each of the four starting points: x = (0.322 87, 0.0092235, 0.046017,0.61817,0.0037169,0.57672, 2.9779) and x = (0.456 93, -0.000 407 198, -0.002 125, 0.915 17, -0.369 57, 2.6105, 11.5). These are the only two solutions that satisfy the constraints on x1 and x4. However, the second solution is infeasible because of the negative mole fractions for x 2 , x 3 , and x 5 . Literature Cited Allgower, E. L.; Georg, K. “Homotopy Methods for Approximating Several Solutions to Nonlinear Systems of Equations”. In Numerical Solution of Highly Nonlinear Problems; Forester, W., Ed.; North-Holland: Amsterdam, 1980; p 253. Allgower, E. L.; Georg, K. “Relationships Between Deflation and Global Methods in the problem of Approximating Additional Zeros of a System of Nonlinear Equations”. In Proc. NATO Adu. Res. Inst. on Homotopy Method and Global Convergence;Eaves, B. C., et al., Eds.; Plenum: New York, 1983; p 31. Bhargava, R.; Hlavacek, V. “Experience with Adapting One-Pmameter Imbedding Methods toward Calculation of Countercurrent Separation Processes”. Chem. Eng. Commun. 1984, 28, 165. Byme, G. D.; Baird, L. A. “Distillation Calculations Using a Locally Parameterized Continuation Method”. Comput. Chem. Eng. 1985, 9, 593. Carnahan, B.; Luther, H. A.; Wilkes, J. 0. Applied Numerical Methods; Wiley: New York, 1969. Chavez, R.; Seader, J. D.; Wayburn, T. L. “Multiple Steady State Solutions for Interlinked Separation System”. Znd. Eng. Chem. Fundam. 1986,25, 566. Chow, S. N.; Mallet-Paret, J.; York, J. A. “A Homotopy Method for Locating all Zeros of a System of Polynomials, in Functional Differential Equations and Approximation of Fixed Points”. In Lecture Notes in Mathematics; Peitgen, H. O., Walther, H. O., Eds.; Springer: 1979; No. 730, pp 228-237. Doedel, E. J. “AUTO: A Program for the Automatic Bifurcation Analysis of Autonomous Systems”. Cong. Num. 1981,30,265-284. Presented at the Proc. 10th Manitoba Conf. on Num. Math and Comp., Univ. of Manitoba, Winnipeg, Canada, 1980. Garcia, C. B.; Gould, F. J. “Relations Between Several Path Following Algorithms and Local and Global Newton Methods”. SZAM Reu. 1980, 22(3), 263. Garcia, C. B.; Zangwill, W. I. “Determining all Solutions to Certain Systems of Nonlinear Equations”. Math. Oper. Res. 1979, 4(1), 1.
Garcia, C. B.; Zangwill, W. I. “Pathways to Solutions, Fixed Points, and Equilibria”. Prentice Hall: Englewood Cliffs, NJ, 1981. Hlavacek, V.; Rompay, P. V. “Simulation of Counter-Current Separation Processes via Global Approach”. Comput. Chem. Eng. 1985, 9, 343. Klopfenstein, R. W. “Zeros of Nonlinear Functions”. J. Assoc. Comput. Mach. 1961,8, 366. Kovach, J. W., 111; Seider, W. D. “Heterogeneous Azeotropic Distillation: Experimental and simulation Results”. AZChE J. 1987, 33, 1300. Kubicek, M. “Algorithm 502, Dependence of Solution of Nonlinear Systems on a Parameter”. ACM Trans. Math. Software 1976,2, 98-107. Li, T.-Y. “Solving Polynomial Systems”. Math. Intelligence 1987, 9(3), 33. Lin, W. J.; Seader, J. D.; Wayburn, T. L. “Computing Multiple Solutions to Systems of Interlinked Separation Columns”. AZChE J. 1987, 33, 886. Morgan, A. P. “A Homotopy for Solving Polynomial Systems”. Appl. Math. Comput. 1986,18, 87. Morgan, A. P. Soluing Polynomial Systems using Continuation for Engineering and Scientific Problems;Prentice-Hall: Englewood Cliffs, NJ, 1987. Paterson, W. R. “A New Method for Solving a Class of Nonlinear Equations”. Chem Eng. Sci. 1986, 41, 1935-1937. Reklaitis, G. V.; Ravindran, A.; Ragsdell, K. M. Engineering Optimization; Wiley: New York, 1983, p 69. Rheinboldt, W. C. Numerical Analysis of Parametrized Nonlinear Equations,” Wiley-Interscience: New York, 1986. Salgovic, A.; Hlavacek, V.; Ilavsky, J. “Global Simulation of Countercurrent Separation Processes via One-Parameter Imbedding Techniques”. Chem. Eng. Sci. 1981, 36, 1599. Seader, J. D. “Computer Modeling of Chemical Processes”. AZChE Monogr. Ser. 1985,81, 15. Seydel, R.; Hlavacek, V. “Role of Continuation in Engineering Analysis”. Chem. Eng. Sci. 1987, 42, 1281-1295. Shacham, M. “Numerical Solution of Constrained Nonlinear Algebraic Equations”. Znt. J. Numer. Methods Eng. 1986, 23, 1455. Smith, J. M. Chemical Engineering Kinetics, 3rd ed.; McGraw-Hill: New York, 1981; p 246. Vickery, D. J.; Taylor, R. “Path-Following Approaches to the Solution of Multicomponent, Multistage Separation Process Problems”. AZChE J. 1986, 32, 547. Watson, L. T. “Numerical Linear Algebra Aspects of Globally Convergent Homotopy Methods”. SZAM Reu. 1986, 28, 529. Watson, L. T.; Billups, S. C.; Morgan, A. P. “Algorithm 652, HOMPACK: A Suite of Codes for Globally Convergent Homotopy Algorithms”. ACM Trans. Math. Software 1987, 13, 281-310. Wayburn, T. L.; Seader, J. D. ”Homotopy Continuation Methods for Computer-Aided Process Design”. Comput. Chem. Eng. 1987,11, 7. Received for review October 20, 1987 Accepted March 1, 1988
