Mapped continuation methods for computing all solutions to general systems of nonlinear equations

Computers them . Engng, vol . 14, No . 1, pp. 71-85, 1990 Printed in Great Britain. All rights reserved

0098-1354190 53 .00+0 .00 Copyright © 1990 Pergamon Press pie

MAPPED CONTINUATION METHODS FOR COMPUTING ALL SOLUTIONS TO GENERAL SYSTEMS OF NONLINEAR EQUATIONS J . D . SEADER,' M . KUNO,'t W .-J . LIN,' S . A . JOHNSON," K . UNSWORTH 1 t and J . W . WISKIN 2 'Department of Chemical Engineering and 'Advanced Imaging Methods Laboratory, Department of Bioengineering, University of Utah, Salt Lake City, UT 84112, U .S .A . (Received 30 January 1989 ; final revision received

15

May 1989 ; received for publication 8 June 1989)

Abstract-In many fields of engineering and science, use of mathematical models leads to systems of linear algebraic and nonlinear algebraic and transcendental equations . When all equations are linear, software packages such as LINPACK and YSMP can be applied to obtain the single solution . When the system contains nonlinear equations, more than one solution may exist, but until recently software packages were designed to obtain at best just one solution from a specified starting guess . If the nonlinear equations arc all of polynomial form, recent software packages such as HOMPACK and CONSOL can systematically locate all solutions . The study reported here addresses the general case where the system may contain nonlinear equations with transcendental terms. By funning a fixed-point global homotopy and applying differential arclength continuation in finite mapped space, the two methods described (toroidal mapping and boomerang mapping) have located all solutions from a single starting guess for all cases studied . a two methods are illustrated for the case of an adiabatic continuous stirred-tank reactor operating in a steady-state mode with two consecutive reactions taking place, one of which is catalytic and irreversible with the other nuncatalytic and reversible . All five steady-state solutions are found by each method. classified by solution procedure . When all equations in the system to be solved are real linear algebraic equations, any of the number of readily-available software packages, e.g. UNPACK (Dongarra et al., 1979) for systems having a dense coefficient matrix, and YSMP (Eisenstat et al., 1977) for systems having sparse coefficient matrices, can be applied to directly or iteratively compute the sole solution or declare a singularity . When the system contains N-I linear algebraic equations and one nonlinear equation, the system can be reduced to just one nonlinear equation and a set of linear algebraic equations . This can be done symbolically with a program such as REDUCE (Hearn, 1987) although the resulting nonlinear equation may be exceedingly complicated . For this type of system, more than one solution may be possible. If the single nonlinear equation is a polynomial equation or can be transformed into such, all roots, real and/or complex, can be found by any of a number of iterative search methods, including that of Muller (1956), which begins from a specified initial guess for one of the roots . If the single nonlinear equation contains transcendental term(s), systematic search methods such as the hybrid method of Brent (1971) can be used to isolate all roots . When the system consists of N - M linear algebraic equations and M nonlinear equations, it may still be possible, by clever symbolic manipulation with a program such as REDUCE, to reduce the system to one nonlinear equation and a set of linear algebraic equations . If so, the previous procedure can again be employed to find all roots .

INTRODUCTION

Deterministic sets of equations for mass conservation, energy conservation, phase equilibria, transport and kinetics used to model problems in chemical engineering may include linear algebraic, nonlinear algebraic and transcendental, ordinary differential, partial differential and integral equations, and combinations thereof. Furthermore, the form of these equations is representative of nonlinear problems in other fields of science and engineering . Prior to the availability of digital computers, much of the fundamental theory for constructing such models was already developed, but only linear or relatively simple and often oversimplified models were actually solved, usually in closed analytical form . With the ready availability of high-speed digital computers and advanced numerical methods, nonlinear and much more complex models can now be solved efficiently . Indeed, when phenomena are relatively well understood, computer modeling is confidently used in lieu of experimental measurements to predict or determine the outcome of process events . Systems of linear algebraic and/or nonlinear algebraic and transcendental equations can be

tOn leave from Daitec Co . Ltd, 4-85 Chikara-machi Higashi-ku Nagoya-shi, Japan during the course of the study reported here . tOn leave from the Department of Mathematical Sciences, University of Dundee, Dundee, Scotland DDI 4HN, U .K_ during the course of the study reported here . 71 CAGE 14,1-F



72

J . D .SEADER ci

01 .

Equilibrium off gas S I

CSTR rY 1 Feed mottract BucH map H2 SO,

0.888 0.098 0.014 1 .000

Equilibrium oft gas u 2

T-150°C 1 .5% acid 1000 gallons

Liquid effluent

CSTR 31A2

Make up BuCH

T° 150°C 4 % acid Liquid product

1000 gallons Make up H2504 MBP

mot tactior 001

Fig . 1 . Process flow diagram for production of dibutyl phthalate .

When a system containing M nonlinear equations and N - M linear equations cannot be reduced to a system containing just one nonlinear equation, the problem can become immensely more complicated to solve unless the nonlinear equations are all polynomial equations- If so, the maximum number of roots can be predetermined from the degrees of the polynomial equations, and software packages such as HOMPACK (Watson et at . 1987) and CONSOL (Morgan, 1987) which use robust differential arclength homotopy continuation methods . can be applied to determine all roots . HOMPACK and CONSOL use multiple continuation path techniques amenable to parallel processing . An example of this case, developed from the 1985 AIChE Student Contest Problem, is shown in Fig . 1, where the acid-catalyzed esterification reaction : monobutyl phthalate+n-butanol- . dibutyl phthalate +water is carried out continuously in two CSTRs in series . Vaporjliquid equilibrium is established in each reactor, with vapor being vented . A known amount of acid enters the first CSTR and additional acid enters the second CSTR to achieve a specified acid concentration in the effluent from the second reactor . The modeling equations, consisting of four phase-equilibria relations and three material balance equations are listed in Table 1 . The denominators of these equations are readily cleared to give the seven polynomial equations in Table 2, where two of the equations are linear . Of the other five, one is thirdtSolutiens at infinity arc nonzero solutions, which arc simply the nonexisting solutions, whose paths have nowhere to go, except to infinity .v as discussed in detail by Morgan (1987) .

degree and the other four are second-degree . The ntaximuin number of roots is readily predetermined by the theorem of' Bezout (see Morgan, 1987) to he F - 2 4 3 = 48 . Application of the HOMPACK program, which sets up 48 paths, gives 44 solutions- with the other four paths going (ill to infinity .! Only one real positive solution, as given in Table 3- is found . Of the remaining solutions, eight contain real finite negative values for some variables and 35 contain complex finite values for some variables_ For the general case of N -AN linear equations and A nonlinear equations, where the latter include one or more equations with transcendental terms that arc not readily eliminated by substitution or approximated by polynomial terms, the number and nature of the roots cannot be predetermined, unless the number of transcendental equations are few In number, in which case see Hoenders tmd Slump (1983) . Furthermore, for this general case . n o reported software package is available for systematically finding all roots and no algorithm with a proof of convergence to all toots is known, although the work of Diener ( I9 7) represents a major step in that direction . Chavez et al. (1986) used continuation from multiple starting points to find all real positive solutions to three examples of interlinked distillation systems . The number of simultaneous equations solved numbered, in one case, more than 350, with moic than 50°% of the equations being nonlinear with transcendental terms . Lin et at (1987) were able to extend that work and generate all the real positive solutions from Just one starting point by continuing the path with a global homotopy until all roots were found . For a global fixed-point homotopy . Kuno and Seader (1988) developed a starting-point criterion, which made it possible . for a number of applications, to compute all real roots along a single homotopy



Mapped continuation methods

73

Table I . Equations for CSTR system 1 . Phase equilibrium for water in CSTRI : 5 .5(13- B) (A-B-C)

B (B+C)

2- Phase equilibrium for butanol in CSTRI : C 3 .5(0.888A -- C - D) (A - H - C) = B + C)' 3 . Phase equilibrium for water in CSTR2 : 5 .5[0.0986A-0 .01(1 .098A-B-9D-E-F+G)-(B+E)l E 1 .098A-B-9D-E-F+G (E+ F)' 4 . Phase equilibrium for butanol in CSTR2 : 3 .5(0.986A-913-F-0 .0986A +0 .01(1 .0 98A - B-913-E-F+ G)] _ F 1 .098A-B-9D-E-F+G (E+F)' 5 . Kinetics of dibutyl phthalate formation in CSTRI : (150 .2)(133 .7 ) (0 .32) D(A - B - C) (A-B-C) = (0 .0986A-D)2' 6. Kinetics of dibutyl phthalate formation in CSTR2 : A(0 .0986)-D (446 .9) (133,7) _ 2. Cr (0 .0032) 1 .098A-B- 9D--E-F+G 1 .098A-B-913-E7. Sulfuric acid material balance for CSTR2 : G- 0 .04 [74 .12(0 .986A-IOD)+222 .24(0 .0986A-D)+18(D-B)+278 .840+98 .09(00136A)J-00136A 98 .01 Where, in Ibmol h - ' : A= inter flow rate, H = Water in off-gas from first reactor . C=BuOH in off-gas from first reactor. D = MBP reacted in first reactor . E- Water in off-gas from second reactor* =BuOH in off-gas from second reactor . G=Makeup H r SO, to second reactor .

path . However, at times, continuation calculations in the complex domain were required . When the starting point criterion was difficult to apply, Kuno and Seader recommended the addition of an auxiliary function to the system of equations to be solved . In the study reported here, the idea of finding all roots along a path from a single starting point is utilized in the development of mapped continuation methods which generate from a single, arbitrary starting guess, a single path . Based upon the successful application of the method to a number of problems, it is conjectured that the path may contain every root to the system of equations . If the conjecture is true, the

significance of the method to modeling problems in chemical engineering or any other discipline involved in the solution of nonlinear equations is as follows : Given a general set of N equations, consisting of N - M linear algebraic equations and M (M > 1) nonlinear equations, one or more of which contain transcendental terms, the proposed method finds all roots. This is in contrast to existing methods, whether they be Newton, steepest descent, Marquardt, Powell or previous types of continuation methods, in which at best a single root is found from a given starting guess . Furthermore, when Newton and/or other methods fail, it is usually not known whether failure

Table 2 . Reduction of CSTR equations to polynomial form Degree of polynomial equation Polynomial or linear equation 2nd 1.5(13-B)(B+C)-B(A-B-C)=0 2nd 2nd 2nd

3 .5(0 .888A-C-D)(B-fC)-C(A-B-C)=0 5 .510 .0986A-B-E-0 .0l(1 .098A-G-90-E-F+G)I(E+F) -E(1 .098A-B-913-E-F+G)=0 3 .5[0 .986A-913-F-0 .0986A+0 .OHI .098A-B-913-F+G)J(E+F)-F(1 .098A-B-90-E-F+G)-0 D(A-B-C) 2 -2056.4(0 .0986A- U)2 =0

3rd Linear

0 .61177-(0 .0986A-D)+0.0l(1 .098A-B-9D-E-F+G)-0

Linear

G-980 ( [74 .12(0 .986A-1013)+222 .24(0 .0986A-D)+18(D-B)+278 .8413+98.09(o

.Ol35A)J-0 .0136A



1 . D . SEAUEK et al .

14

Table 3 . The single real finite solution to CSTR problem (all values in Ihmol h - ) = 80 675 .A -=inlet flow rate B=water in off-gas from first reactor - 6,9841 C = HuOH in off-gas from first reactor - 61 .227 7 .2161 D=DBP formed in first rector 0 .5475 E-- Water in off-gas from second reactor F=BUOH in off-gas from second reactor 3 .6554 O- Make-up H-.SO, to second reactor 0.0591

Three choices for g(x) as discussed by Wayburn and Seader (1987) lead to the following respective homotopies : Newton homotopy :

H(x,t)=f(x)-(I-i)/'t °) ;

(3)

fixed-point homotopy :

(4)

H(.,t)=rf(x)+(I-r)(x-x°) ; scale-invariant affine homotopy : is caused by absence of a solution or inability of the method to converge . The method proposed here, when applied properly, has to date never failed . At worst, it simply completes a path with no solution found because no solution exists . In some cases, the user can specify domain constraints, as discussed by Kuno and Seader (1988), in which case only certain roots are sought and found, if they exist, For example, the domain might be constrained between U and I or in the real positive domain for all or some variables . Although the idea of the conjecture is simple, its implementation can be quite complex because of the nature of the paths encountered . However, once written into a computer code, the proposed method is very reliable . Its only drawback is its potential inefficiency in obtaining a rout when compared to less general, current state-of-the-art methods . Thus, the proposed continuation method is only recommended when other methods fail or when multiple solutions are suspected and desired . Typically, to obtainn a single root, the proposed method, which can supply its own arbitrary starting guess .. takes from three to ten times as much CPU time as Newton-type methods, when the latter converge without difficulty . In the following sections, two implementations of the method arc described and applied to the classical problem of multiple steady-state solutions of a CSTR . DEVELOPMENT OF'THE MAPPING FLNCTIONS

H(x,t)=tf(x)+(1-1)f(x°)(x- .t°,) .

(5)

Equation (3) has an advantage that closed paths for 2 on t are possible within a finite domain . However, because H(x, t) =0 may have multiple roots for f(.x)=f(x°) at t=0, more than one path may exist . Equations (4) and (5) have the advantage that H(x, 1) -= 0 is satisfied by only one root at t = 0 and only one path for a on t can exist . However, as will be shown, the path can consist of branches that are only connected at infinities or by branches in the complex domain . As further shown, branches connected at the infinities traverse x and or t from - x, to +cc. . However, in the algorithm developed here, the fixed-point homotopy is applied within a finite domain by utilizing mapping functions on x and/or t to obtain a single path . We note that the conjecture that all roots to equation (I) may lie on the path is equivalent to the nonexistence of so called isolas in the solution to the system of ordinary differential equations that results from the homotopy method . To illustrate the advantage of using mapping functions, consider the simple fixed-point homotopy example, H(v,t)=t('x'-3x +2)+(I-t ) (x-x°)= 0, given by Lin ei al . (1987), where it is desired to compute both roots (I, 2) of1"(x)=x z -3x -1 2 -U . If we let .v° = 1 .5, the homotopy path shown in Fig . 2 is obtained . Notice that the path has three branches . which, however. are connected at the infinities_ The points of connection are : (I) t=0 at v=±~xr 12) x = 2 .7071 at t-= + a.. ; and (3) t = . 1 2929 at

Let

f(x)=0

(1)

,o r

be a system of nonlinear equations with solutions x e R" . where 5

f: R"-R" .

Branch 2

Embed equation (1) into a convex linear global homotopy vector function to give : H(xrr)=(f(x)+(I-r)g(x)=0,

I Rn^t I

sla>rinc j po~~.t

(2)

/Branch

Rnnt ?

3

where t = scalar homotopy parameter and t e R', e

,

2

T

g(x)

= vector function selected such that the solution to g(x)=0 is known or easily determined .

Fig . 2 . Homotopy path for fixed-point homotopy of xr-3s r2 with a starting point of z°= 1 .5 .



Mapped continuation methods t = ± co . To numerically track such a path, it is desirable to transform the infinite path to a finite domain . Desirable properties of a function defining such a mapping are the existence of a continuous, finite first derivative and the existence of a function defining the inverse mapping . Toroidal hypersurface mapping One such mapping function that we have used with success is : Ay Y (6) (B+y°) where y may be x and/or r in equations (4) or (5) . As y approaches ± co, the asymptotic expansion of y' is given by: y'- Ay - l"(1 - B qyv + - - -) . The requirement that 0 < lim y _ ly'I < oc dictates that n =p/q and, thus, equation (6) is reduced to : yv " Y -A B+y" (7) Because the derivative of equation (7) is : Y °- i dy'/dy =AB p q (B+yn)a+ i the requirement that 0 < dy ;dyly.=a < ac, further dictates that p=q, giving the final form of the generalized toroidal mapping : (B Ay + y v)"v

(8)

Note that, for p even and p = q, (y°)"" has two possible computational interpretations: (y°)"9 = +(yf''a, which is conventional or (y")""=(y)e19=y . The first interpretation produces a mapping which is nonnegative, has a discontinuous first derivative at the origin, and does not have a one-to-one mapping characteristic . The second interpretation has no such problems in the form used in equation (8) . We also require p to be even since (B+y °)=0 is possible if y < 0 for odd p . Note that p > q yields a general toroidal map that does indeed have : d = 0. dy r =a This property allows the use of the absolute value . The inverse can still be obtained in closed form, although it is, of course, double valued . The advantage of using the form in equation (8) is that it obviates the need to jump from ± co to -I - cc . We refer to equation (8) as the generalized toroidal mapping . It is restricted to the domain - A C y' C + A- The inverse of equation (8) exists in the following closed form : B''Dy'

Y = (9) [A"-(-' )"] '°

75

For initial studies, we set A = B = I and p = 2 to give : (10) y'=y/(l+y2)" : This function, which we call the toroidal mapping, is restricted to the domain - 1 C y' C + I and causes the following transformations : Y -co -1 0 +I +cc

Y' -1 -1/f=-0 .7071 0 + 1/v' _ + 0 .7071 +1

The inverse mapping is (11) y =y'/[ 1 - (y') 2] U2 When the toroidal mapping is applied to both x and t for the homotopy in Fig . 2, the result is Fig. 3, where the same three branches are clearly discernible and can be considered to be connected at ± 1, which are the mapped locations of the infinities . However, the complete path is now of finite length and is more rapidly tracked . The starting point is located at t = 0 and the roots are located at t' = +0 .7071 . The choice of the term toroidal to describe this mapping derives from the following observations . Let y=[yl, . .-,y_--,YjT and Y '_ (y„ - - , y;, . . . , y ;, ]T be vectors whose components are related by equation (8), where y'= ±A in the limit or y -+ ± m, respectively . If the boundaries y,' = ± A and t' = ± A (using the same mapping for t' as for y') are joined, then the resulting manifold of dimension n + 1 is a hypersurface with toroidal topology when viewed in a space of n + 2 dimensions . On the surface of the hypertoroid, the mapped homotopy path is continuous, while in the n + I dimensional Euclidean space (y', r'), the homotopy path is continuous only within the domain and terminates at points on

root Branch 2 Brunch 1

_ Starting point Root 2

(X 2-3X+21+ (1-T) (X-X° ) -o X°- 1 .5 X'- X/t1+XCY'/ 2 T'- T / 11+T21 a2 -0 ..5 Branch 3

-1 .

0

-oS

0.5

Fig. 3 . Application of toroidal mapping to homotopy of Fig . 2 (mapped coordinates) .



J . L) . SLAVER et rr! . the domain boundaries, the boundaries defined by i : = + A and I'=+A .

10

Boomerang mapping The toroidal mapping requires one, in Euclidean space, to switch from one branch to another whenever the domain boundary at y' - ± 1 is reached . A continuous path in Euclidean space can be achieved by using a mapping that transforms -x and +x into the same finite value . One such function that we have used with success is :

0.e

Y' =21"41 I-y'), (12) where again y may be x and,/or t in equations (4) or (5) . This function, which we call the boomerang mapping, restricts the function to the domain -- I < y .' _ . I and causes the following transformations : V

Note that the maximum and minimum values of v' occur at y = ± 1, respectively . The mapping does not have a one-to-one characteristic (y'=0 maps onto both y = ±oo and onto y=0 ; finite, nonzero values of y' map onto two t= values) . The mapped homotopy path is similar to the flight of a boomerang with a maximum range of - I for each spatial coordinate . The boomerang mapping also has a continuous first derivative everywhere and an inverse mapping, which is given by : Y=jl±[I-(1'}~)'''}, .'y (13) When the boomerang mapping is applied to both .x and r for the homotopy in Fig . 2, the result is Fig . 4 . where a continuous closed path of finite length is now evident . Starting from t'= 0 (t = 0) and x'- 0 .923077 (x = 1 .5), the single path includes both roots at I' = 1 . FOLLOWING THE PATHS OF THE MAPPING FUNCTIONS Methods for following the path of a homotopy function are discussed by Wayburn and Sender (1987) and Lin et at. (1987) . The method used in this development, as summarized next, employs a predictor method applied to a differential arclength form of the homotopy and a corrector method applied to the homotopy itself. The honiofopf predictor equations Let x = x(s) and t = i(s), where s is the arclength along the homotopy path . Then, the derivative of

Start mg Point

• Root 1 ~• Root 2 I

X' 0

I Tixz-

-os

'l1X XOA~ x°= + .5 X, - 2x/(X2 t 1) T' - 2T?( T' t U

t Fig . 4 . Application of boomerang mapping to homotopy of Fig . 2 (mapped coordinates) .

equation (2) with respect to arclength is : BH dc, FUI or --=0 . et ds , ',(IX, ds Note that dH=t cl +(l -t) ng 0,H and - =f (x) - .4 (x ) % .x, The, CA, dt Since f and g are vectors of tr components, the meaning of df/cx,,is : c.r,l - rr and we interpret r19/0v, in a like fashion . The arclength relation may he written : /dx, ~dx;') (dt)(dt I\ ds ds .,, ds d The above two equations can be solved simultaneously for the unit vector s, = [dx,ds ;;dt1s) however, this is itself a nonlinear problem since the second equation is quadratic . A more practical approach is to note that for small step sizes [(dvids)" '', (dr ds)'r-"][(dx/ds)"-r, (dt ;ds)"s]T~ 1 where superscripts (k - 1) and (k) are the previous and present step positions on the llontutopy path . Using this notation and upon combining the ahove two equations, we obtain an approximate, linear, composite equation in matrix notation : dH kl Cc a `-" r ds

J_,711 ntlJ(k)

dx, dsdt ds J

dt '"-'' [0 ds JI This square matrix is of size (n + 1) by (n + l) . Solving this linear system for the unit vector utkl=[dx1ds)'-1, (dt/ds)"')r is of particular significance because a is the unit tangent vector to the homotopy path in R"-' . The predicted point



Mapped continuation methods (x, t)t'`+o is given by (x, t)(k I n= (x t)(k) + hu ( k) where h is the step size chosen . The homotopy corrector equations The above (Euler) predictor method is of first order . Higher-order predictor methods (e .g . Adams-Bashforth), using the values of f and g and their derivatives at multiple previous points, are also known . All predictor methods are approximate and must be followed by a correction step to return more closely to the homotopy path . The classical method is a Newton correction in the t = x'k+ ' ) plane which, however, may fail near turning points where u is or is nearly orthogonal to the t unit vector. A more robust corrector is a Gauss-Newton, least squares, unconstrained correction to the honiotopy path or a Newton-Raphson corrector constrained to lie in the hyperplane perpendicular to the homotopy path and passing through (x,t)'k+° . Finding this hyperplane requires knowing the corrector solution . However, a practical approach is to approximate the correction in the hyperplane by a correction in a practical hyperplane passing through (x,t)'k+ " and perpendicular to the predictor vector u~" ). For small step sizes, these planes are nearly identical . Choosing the latter path we let [OX, 3t] approximate the step from (x. t)tk' 0 to the homotopy path . Then, on applying the Newton-Raphson method, we obtain : OH cH

ax'

a /r

ax,

,

ar tJk,

Cxx l ,k. --

Kubicek (1976), Abbott (1980), Watson and Fenner (1980), PITCON of Rheinboldt (1986), Shacham (1985) and HOMPACK of Watson et al. (1987). In this study we have used HOMPACK and a modified version of Kubicek's algorithm, which can be applied to differential arclength continuation using an Adams-Bashforth predictor and a Newton corrector, as described in detail by Seadcr (1985) . Kubicek's algorithm was amended to perform the predictor step in the mapped space and the corrector step in the original space . To solve the predictor system of ordinary differential equations in transformed space, expressions are needed for dx ;/ds', (i = I, . . . , n + 1), where s' is the transformed arclength parameter, and x,:, the transformed variables (with t' as x ;,, ,) . But dx, ldx, dx,\/ ds (14) I ds' \dx,Ads)(ds' and dx(/dx, is easily obtained from equation (6) as : dx ; 1 (15) it -Y, (1 + x')'rz Kubicek (1976) gives expressions for dx ( /ds, and since : z (16) we obtain, from equation (14) :

at' -

rat, , k+ ., -

Oak+o,

L

[d

The new value of the corrected point given by :

(x,t)tkA

2) is

(i t)tk+2) = (x,t)('N1+(dx,6,)(k+1),

This correction step is repeated until the norm of H is smaller than some selected small number . Clearly, the predictor or corrector steps can be made in the original space (x, t) or in a new mapped space (x', t') . This option of choosing the space in which to make a prediction or correction step is exploited in the next sections which describe application of the two previously-defined mappings .

dx; ' -ta

ds

1]dx '&]T -H_

The constrained condition that the path lies in the chosen practical hyperplane is ]45x, dt] u (k) = 0 . On combining these equations, we obtain a linear matrix equation : aH Ik+n off fk+q

77

Therefore, dx ; ds'

1 1 dxr z (1+x?)'a{ ;=, [(1+r)3'2ds, (i=l, . . .,n+l) .

A number of computer algorithms are available for implementing the path following procedures just described . These algorithms include DERPAR of

n2

dx ds (18)

Equations (18) are used to solve for the next point on the predicted path of the homotopy function using the Adams-Bashforth method, which reduces, for the first-order case, to a Euler predictor. If a value Ix ; I > 1, (1 ~ t <-n + 1) is predicted, a jump to another branch of the solution is signaled . Therefore, before correcting the solution, the jump is completed by simply adding (subtracting) 2 to x ; according to x ; < - I (x' > 1) . The Newton corrector may require more iterations to converge immediately following a jump. Suppose that the predictor-corrector continuation method has produced the approximations x' (j =0, 1, . . . , k) to equation (4) at successive points along the continuation curve . Then, one more application of the Euler predictor gives : x k+I_ xk+hk+l

Toroidal mapping

(17)

I 32 ds L (1 +4)

dxk

(19)

where the E subscript denotes the Euler prediction of xkI' and superscripted indices are not enclosed by parentheses . The ' (pruned) notation has been



1 . D . SEADER ci ¢I .

78

dropped, but as explained above, the predictor is applied in the transformed (primed) space . Hence, at each Enter step a suitable step length, h' ' (k -_ 0), must be chosen . Kubicek (1976) does not include an efficient method for continuously adjusting h . Therefore, an attempt was made to implement a step size algorithm which produces a value for h'=" by adapting itself to the behavior of the path . Such a method has been suggested by Rheinboldt (1980) and den Heijer and Rheinboldt (1981) . Rheinboldt (1980) proves the existence of hk" > 0 such that xp*' lies in a hall, guaranteeing that with this value as a starting point the Newton iteration is well-defined and will converge to a solution which lies in the same hall . lie then describes a technique for actually obtaining such a value for h''', i .e . one such that :

-

where p"" is an estimate of a "safe" radius of the hall, given by h

p Il w ll . l

A" C

( 20)

k

Clearly, small values of Ia`J need to be avoided . Therefore, the denominator under the square root in equation (21) is replaced by : w"= 21sin

a, rr ;2)I .

(23)

where xif s'Ca .rifa'Cs'Cfl VseR' //ifs >(

2(s_a,/I)=

and a a (0,n ;2) is a lower threshold of a' J . Now let rr k +' = p "'iAv k . Then we may obtain the following bounds for rr" - ' from equations (22) and (23) :

x

2 2x sinx!2_ a - ~sin>< :4=v-

2

a_°Hence

is given by o rr

r

if o,,,a
k

o'/dr ot I V `2 'C sin a'1 d •A x°otherwise, (24)

a +'=

bk

where = 14 - x"a i . Finally, we obtain :

where Wk 1! ,= A

h"" = A '-' Ax -'

a l:

7

stn

,

k

where

k

dx i r/X -x _ '1 ds Ax` }. and as before Ax'=fix'-x' 'hi . A careful analysis of the Rheinboldt method reveals that the heuristic justification for the particular choice of ak as above is that "m"-I . approximates the curvature of the homotopy path . A different (more conservative) estimate of the curvature results from setting : a' = cos - '

a'=cos -'

-dS

d 4-T .

lL J Cd At this point, however, numerical experiments do not yield a significant increase in efficiency with this second approximation . Three methods for selecting p" - ' are proposed by den Heijer and Rheinboldt (1981), but in the numerical experiments of this study, their Method Ill appeared to perform consistently better than the other two and was the only one implemented . This was done as follows : equation (20) may be written as : hl`- I = rIk `' AF ,

A

where

and we enforce the condition that (x > 1) .

k

with w given by equation (23) and rk - ' by equation (24). den Heijer and Rheinboldt (1981) regard equation (25) as the "tentative step length" ; they indicate that it may subsequently be altered to ensure that . h,,,,,,

(22)

h"`'

h,. .1

where Iz,, and hm „ are user-supplied bounds on In our implementation, we check to ensure that by using equation (25) the difference in magnitude between each component of x, and x' will be no more than some user-specified value . Since the predictor is used in the transformed space . and since each component of 4' will therefore lie in the range ( - 1, t), a value of 0 .1 for this difference seems reasonable- If this condition is not satisfied, the value of h" " is reduced accordingivOnce h +' has been chosen and xt"' obtained . the corrector iteration starts from H The method allows a user-specified maximum number of iterations in which to converge (typically five) at each step, although this may need to be relaxed when considering paths in the complex domain as discussed below . If the method fails to converge, then It' "' is halved and is recalculated . The algorithm described here has been used on a large number of examples, with a = 0 .05 and w = 3 , as suggested by den Heijer and Rheinholdt (1981) . The use of this step-size algorithm significantly improves the efficiency of the Kubicek algorithm .

H

5

4' I

2sin s,

1 S A k " S x-

(25)

4'



Mapped continuation methods Homotopy paths are not always confined to real space . The original system of equations may have complex roots (possibly as well as real roots) . Secondly, the original system of equations may have only real roots, but at some stage of the computations the homotopy equations may have complex roots for certain values of t . For example, the equation x' - 2-v' + 4x - 8 = 0 has a single real root of 2 and complex roots of ±2i . The system of transcendental equations given by Allgower and Georg (1983) : (x, - xl)(x, - sin x,) = 0, (cos x 2 - x, ) (x2 - cos x,) = 0, where x, and x, are in radians, has an infinite number of roots, including four real roots in the region :

79

10

0

ranch

2

4

Branch

I a .-I .__

arena,

3

•

x o Ttx a - 25X . 2)+l1-TI

-5

(X-x°l=o

5 0 .1 .5

-o T

Fig .

5.

Application of switching between original and mapped paths for homotopy of Fig. 2 .

0
However, to reach the four roots, complex domain computations somewhere within the interval 0 < t C 1, as shown by Kuno and Seader (1988) tray be necessary, depending upon the selected starting point . To deal with problems in complex space, two possible options seem appropriate . First, the entire complex computation can be done directly in complex arithmetic; or secondly, the original equations can be put into real form (thus generating twice as many equations) and solved for both the real and imaginary parts of the roots, both being treated as real variables . Followin the second option and writing x = u + iv, (i = - ), the equation x 2 - 3x + 2 = 0, solved in Figs 2-4 becomes the two equations : a 2 -t: 2 -3u+2=0, v(2u-3)=0 .

This system is then solved for u and r . When following paths in the complex domain, additional Newton iterations may be required when the solution attempts to jump to another branch, because when t changes sign, the imaginary part o changes from nonzero to (approximately) zero . A suitable way of dealing with this is to provide extra Newton iterations at a jump . As a jump is approached, the number of Newton iterations is allowed to increase . If convergence is not obtained, the predictor step length is halved and the prediction repeated . Once the jump is successfully negotiated, the maximum number of Newton iterations is reset to the original value. Boomerang mapping

With this type of mapping, all solutions are located at turning points on a boundary of the single path at t' = + 1 . Because the defined region is y ; < 1, the variables may tend toward the undefined region if a predictor-corrector step is used . This disadvantage of the continuous mapping can be overcome by switching

back and forth between the original and mapped version of the path while tracking it . First, a set of y,,, greater than one is chosen for t and the set of variables x„ where c are arbitrarily chosen limits for calculating in the original space . The corresponding mapped values, yare from equation (8) : (26) which arc located inside the boundary of the mapped space . If y J is less than y,, r , then the variable y; is tracked in the original coordinates ; otherwise, this variable is transferred into the mapped space . To illustrate the details of this idea, consider the previous example, H(x, t) = t (x' -3x+2)+ (1-t)(x-x°)=0 . If we chose the limits of the original space as x, = 5 and t, = 1 .5, with x " - 1 .5, the homotopy path of the original space is shown in Fig . 5 . The mapped variables x' and t' arc defined as : J

x'= 2x/(x2 +1),

(27)

t'=21/0 2 +1) .

(28)

At the starting point, x° and t are less than x,, and i,. . The path is followed according to the original function H(x, t) = l (X2 - 3x + 2) + homotopy (I - t) (x - x°) = 0 on Branch I of the original space until the point Q1 .5, 2 .1937) is reached, corresponding to t,. and x = 2 .1937, which is less than the limiting x, . Because the value of t is greater than t,.

Table 4 . Application of the boomerang mapping Computation Itomotopy branch Region function I 1 .2 2 23 3 3 .4 4

A-C C-D

D-E E-F F-0 G-H

H-A

H(x,t)=0 H(,'')=0 H(x,t)-0 H(x',t)-0 H(x,r)=0 H(z, r')=0 H(x Q=0



J . D . SEADER et al.

80

The second reaction, which is reversible and noncatalytic, is :

Table 5 . Results for adiabatic CSTR problem Concentrations (kgmvbs) T(K) C C, Cc 2 78732 '_38042 0.126397 0.003901 0.000380

310-213 333 .493 452 .569 594 .027 691 .624

0.212680 02619577 2 .84991 1 .71361 0 698897

6.4686 x 10 ° 1.6427 x 10 ° 0.023694 L28256 2 .11073

B with rate expression r,=k7C E -k_C c ,

beyond the point C . the mapped i-variable homotopy function : I t [l +

(t") 2]°°

H(x,1')=-

(33)

where the rate constants, which are consistent with chemical equilibrium, are given in terms of temperature as : ki =3 x 10 4 exp(- 80,0008 .314T) .

(34)

k==3x 10'exp(-90,000/8 .314T),

(35)

--(x 2 -3x12) as (29)

is traced until point D ( - t .5, 3 .30-5) is reached on Branch 2. The path is then tracked in the original space between points D and E on Branch 2, at which point x = x . = 5 . Now, x will become greater than x,, so the variable x is mapped . The homotopy path becomes the solution of H(.x', x)=0. This process is repeated until a closed loop is reached . The procedure is summarized in Table 4 . For any branch, if either x or t are found to be outside their limiting values . then that limiting value is mapped and computation of the branch is repeated . To follow the path, Kubicck's algorithm was modified so as to swap between the original and mapped coordinates automatically . In Fig . 5, two roots are found at x = 1 and 2 . Although Fig. 5 shows at most one solution (at I = 1) . to the original equations for any one computation branch, more than one root could be located on any one computation branch . The example of Fig . 5 also leads to a closed loop formed by the computation branches of Table 4 . Two other types of paths are possible : one fails to form a closed path and the other requires some computations in the complex domain to bridge certain computation branches .

APPLICA'1'IOr

Consider an adiabatic continuous stirred-tank reactor (CSTR) operating in a steady-state mode with two consecutive reactions taking place . The first reaction, which is catalytic and irreversible is : B -f R,

with the rate expression r,-k :CAr(l-I KAC,),

(30)

where the kinetic and adsorption equilibrium constants are given in terms of temperature as : k 1 =4x 10'exp( - 6000/8 .314T), KA = 17 exp( -- 7000/8 .314'/) .

(31) (32)

In equations (3(f-35) . SI units are used . The mass balance equations for the CSTR reactor in terms of residence time v are (36) n (9k, C',, and

C,--Cc,=Hk_Cs-©k_C6 .

(38)

The enthalpy balance for the CSTR reactor, which accounts for sensible enthalpy and enthalpies of the two reactions in terms of quadratic expressions in temperature, is in SI units : 85(T- Tt,) I- 0 .02(T 2 - 7 ;) (16,000 } 31 -0 .00277(C A" (30,000+4T-1)003 T 1 )(-

(39)

Equations (36-39) constitute a set of four highly-nonlinear equations in the reactor effluent quantities C 5 , C'a , C, and T in terms of parameters consisting of reactor inlet conditions C,,,,, C B ,, . C< ., and T4 and residence time 0 . Depending upon the values selected for the parameters, it is well-known (e.g . Biious and Amundson, 1955) that equations (30-39) can have multiple roots. For the particular case here, as many as five steady-state solutions are possible . As presented in aa number of textbooks (e .g . Smith, 1981 ; Fogler, 1986) . a commonly used method for determining the multiple steady states is to substitute equations (36-38) into the right-hand-side of equation (39) so as to eliminate all unknowns except T. Equation (39) then has the form_

G(7) = F(T) .

1,40)

By plotting G(T) and F(T) each against T on the same graph, the multiple steady states are located graphically at the intersections of the two curves . Alternatively, equation (40) can he treated as a single nonlinear equation with transcendental terms, which



Mapped continuation methods

R2 RI Branch

,

R5

200 Branch 5

B anch 3

Fig . 6 . Homotopy path for CSTR problem using toroidal mapping (original coordinates) .

can be solved for all roots by the method of Brent (1971) . The elimination in equation (39) of all unknowns except effluent temperature T is readily accomplished if the equations utilized to eliminate the effluent concentrations [e .g . equations (36 38)] are linear . If not, it may be difficult or impossible by substitution, for a complex reaction network to reduce analytically the set of CSTR equations to one equation . In the application here to the toroidal mapping, the variable C, was eliminated, reducing the set of equations to three nonlinear equations in CA , Ca and T, consisting of equations (36), (37) and (39) after eliminating C c by using the following rearrangement of equation (38) : C,1,+Ok,C0 CC = (4 1) 1 + Ok ; The set of three equations, or further modification thereof, was solved by using the global fixed-point homotopy continuation with the toroidal mapping for all real roots for the case of_ To =298K, CA , = 3 kgmol m - ', C R, = Cc, = 0 and O = 300 s . The resulting multiple steady-state solutions for both mappings are listed in Table 5 . Solution by toroidal mapping The solution by means of the toroidal mapping was achieved as follows : First, a necessary preliminary step involved the elimination of the denominators . Numerical instabilities appeared to be insurmountable when this step was not taken . There are many ways in which denominators can be eliminated . One result is : AC, +(1 +k,e)(CA-CA,)

=0,

( 42)

Sl

AK, CP+AC5 -k,O(1+k,e)CA =0, (43) 85(T - Ta ) + (0 .02) (T 2 - Ta ) + (- 16,000 - 3 T + 0.002T2) (C,b - CA ) JCAn +(-30,000-4T+0 .003 T 2) a (CA, - CA - C B N CA, = 0, (44) where A = I + O (k,,, + kz ). Equations (42-44), with equations (31-32) and (34-35) substituted for k, . k s , k ; and KA were the equations then solved . To avoid numerical problems for small and negative values of T, the quantity T in all terms of equations being solved, including the fixed-point term (1 - t)(T - Ta ) of the homotopy equation, was replaced by 171 . Figure 6 shows the homotopy path for T as a function of the homotopy parameter t for the arbitrary starting point (CA , C a , T, t) _ (1 .9, 1 .9, 250, 0) . The direction of the path, indicated by the arrows, consists of five branches as numbered in Fig . 6. The first branch begins from the starting point Sand finds three of the five real roots before moving tot = F 00 . The roots RI . R2 and R3 found for T are 462 .569, 594 .027 and 691 .624, respectively, as listed in Table 5 together with the corresponding roots for the other variables, CA , Ca and Cc . Branch 2 begins at t = -eo, for the same value of Tat which the first branch reached t _ +a, and terminates at T=+x . without finding a root . Branch 3 . which is a mirror image of Branch 2 with respect to the T-axis, begins at T = -oo and terminates at t =-rr: without finding a root . Branch 4 is a mirror image of Branch l . Branch 5 begins at t = +oo and finds two roots, corresponding to T = -333 .493 and -310 .213 . Like the three roots of Branch 4, these are not roots of the original set of equations, but arise only because of the above-mentioned transformation of T to DTI- The other two roots, R4 and R5, at T = 310 .213 and T = 333 .493, respectively, are then found before Branch 5 tends to t = +cc . Branch 1 is resumed at r = -oc and returns to the starting point without locating any other roots. An interesting bifurcation phenomenon occurs when the starting point xo =(CA ,CB ,T) is varied . Below critical starting values for the variables, the homotopy path is traversed in the direction of the arrows of Fig . 6 . Above these critical values, the path is traversed in a direction opposite in direction to the arrows shown . Therefore, depending upon where x o is located with respect to the critical starting point all roots would be missed unless the path were allowed to continue and jump from ±o to T c- as in the procedure described here. The critical value of xo when To is fixed at 250 is (1 .99511514633, 1 .99511514633, 250, 0) . Another important point is the utility of employing values for A and B not equal to unity in the generalized toroidal mapping associated with T . That is, AT /B+T'



1- D . Ssnooa e( uL

n2 tss

e Brunch s•R' IN R4

Brand 6 60

p 4 E

S R3

8-Ch 4

Branch

2

Branch 3 (

s

I

Branch t

20

0

Homotopy Fig .

2

parameter fhi

7 . Homotopy path for CSTR problem boomerang mapping (original coordinates)_

using

For example, for A = 10' and B = 10 6 , the path is completed in 823 steps . On the other hand, if A = 1, B = 106 (corresponding to a simple scaling of "T by a factor of 1000) the path is very far from complete after 800 steps . The reason for the greater elllciency of the full generalized toroidal mapping (A -4 1) as opposed to simple scaling (A= 1) is that the curvature of the homotopy path in primed coordinates is increased enormously with simple scaling and concomitant rapid decrease in step size (as per Rheinboldt's algorithm) greatly hampers progress . Solution by boomerang mapping It is possible to reduce the system of equations (42-44) further to just a single, but very complicated nonlinear equation in temperature . This was accomplished readily with the symbolic algebra program REDUCE (Hearn, 1987) . The result was equation (44) with the following expressions for C,c =C„(T), CB = CB (CA , 7') and Cc = C C (CB , T) : CA- - [(T'k=ki+2T'k z k_k 4-T"(k ;) 2ki +2Tk k,+12T'k,k_k,K,,+4T'k2 k_k, +2Tk,ki+12T'(ki) 2 k,K +2Tr(k') 2 k, + 2Tk 2k i + T=kd + 2T 2k z k_ + 12T 2 k,k ; K + 4Tk,k, + T2 (k_ )2 + 4T2k2 k, + 24T2kik, KA + T2 ki + 2Tk-, -I- 2Tk : + 12Tk, KA + 2Tk, + 1)1'`- T2 k,k, -T-k,k, - Tk, - 6Tk _ KA - Tk ; - Tk, -6KA -lJ/[2KA (Tk,'+1)J, (45) Ca =-(CA Tk_+ C,,-3Tk,'-3)/(Tk,,+Tk ;+I), (46) Cr.=(C,Tk,)l(Tk4+1) (47)

The resulting equation obtained by combining equations (31-32), (34-35) and (45-47) was then solved with boomerang mapping as described above from a starting guess for the temperature of 300 K . The resulting homotopy path is shown in Fig . 7, in terms of the unmapped coordinates . with the branches of the path marked from one to six . The starting point is labeled S and the five roots for T, which are identical to those in Table 5, are labeled RI-- R5 . Starting front Sand moving upward and to the right, root RI is found before Branch I moves toward t=+w- where T-TB , . Branch 2 is then tracked from T='TB , starling from t=-x': and moving toward t = 0 . However, before reaching t =0, Branch 2 turns upward and then to the left toward r=-eu, where T=TB_- Branch 3 begins at (t = + A ., T=I er) and finds two roots, R2 and R3 before moving back toward t---+x. at T=T„ . Branches 4 and 5 are similar to Branches 2 and 3 . but' are located at higher levels of T . Finally . Branch 6 is tracked starting at t = -rr, and T-TB , . Before reaching t =0, Branch 6 turns upward toward T = + s, for t just less than zero and equal to t,, . A closed mapped curve is obtained by moving to T = -cc at i n, on Branch I and moving upward to the starting point S . The toroidally mapped homotopy method was also applied to the above single equation in T and gave identical results . It is instructive to plot the modulus of the residual of this function as gray scale values in the (T, t') toroidal domain . Let black be zero, white be 250 or greater and let intermediate gray values lie proportionally between these values . Such a picture can be displayed on the video screen of many workstations or personal computers . A photograph of such a display is shown in Fig . 8 . Solution by other globed fixed-point homofopt merritrids The above described CSTR problem was also solved with a modified version of HOMPACK by applying the mapping concept from 4-x to - ,- or vice versa . but without using a mapping function . This was accomplished by computing separately each branch of the homotopy path until an apparent asymptotic limit was reached as the path vended to + :L: or -x' in v or t . This limit was then mapped manually into the next branch at the opposite point of infinity . For example, in one case- while tracking a branch, nearly-asymptotic values of T =353 .1 and CB = 0 .9686 were achieved at r -- 7 . Computations for this branch were terminated and another branch was initiated at T = 353 .1. Ca = 0.9686 and t = -7, with the path moving in the direction of increasing values of r . Three different forms of the equations were solved in this manner . The first form involved the same three equations [equations (42-44)] used for the solution by the toroidal mapping. as discussed above . However, T was not replaced by TI in the fixed-point term . Accordingly . the paths centered



Mapped continuation methods

83

(T)

d V N La, aEE

-ov n n sm

(150C) 1 (~)

0 (0)

+1 ( c)

t' (t)

Mapped Homotopy Parameter Fig. 8 . Photograph of residuals of homolopy path in the toroidal domain for CSTR reactor problem . around T = 0 were not symmetrical as they are in Fig . 6 . However, the same five positive solutions were obtained . The second and third forms of the equations, obtained by eliminating one more variable, involved two equations in two variables . The set in T and Ca was obtained by solving equation (36) for C A followed by substitution of the resulting equation and the stochiometric relation C c = CA , - CA - Ca into equations (37) and (39) to eliminate CA and Cc . The resulting equations, after clearing denominators were : AKA CB+A(l +k,O)CO -k,OCA „(1 +ki0)=0 (48) and [85(T - T0 ) + 0 .02(T2 - To)) (1 + k, O )C, 0 + (- 16,000 - 3T + 0002T 2 )AC, +(-30,000-4T+0 .003T t)k,OC a =O. (49) A set of equations in T and C A were obtained in a similar manner- The second and third sets were solved successfully from the same starting point used for the three-equation set discussed earlier after replacing T with ITI, except in the fixed-point term . The three-equation form, given by equations (42-44) was also solved by the method given by Kuno and Seader (1988) . The variable T was not replaced by DTI, but Twas scaled by a factor of 100 . Application of the starting-point criterion developed by Kuno and Seader to give a starting point of CA = 10 . Ca = 10 and (T/100) = 5 proved to be successful in

locating all positive roots on a single path . Similar success was achieved by using the single-equation form given by equations (44-46), where T was again scaled and a starting point of (T/100)=2 was selected, to obtain all five roots on a single path . STEREOSCOPIC DISPLAY TOOLS

For a function f(x) of a single independent variable, the homotopy path and the value of the modulus of the residual of the function h(x, t) may be plotted in the 2-D domain (x, t) as seen in Figs 2-7 . A function f(x,, .x2 ) of two independent variables generates a homotopy path in (x, x,, t) space which may be plotted as stereo-pairs as shown in Fig . 9 for the CSTR reactor example discussed above . The solution to this example for just temperature is given in Fig . 6 . Two variables, CA and T, are shown in Fig. 9 . When the C A axis is positive into the plane of the page, the similarity to Fig . 6 is apparent, but the dependence on CA is also clear . For domains of higher dimensions than two, stereo display of all independent variables is impossible, but display of a projected domain (u, v, t) is possible where u=u(x„x x„), v =v(x,,x x,J . The variables u and r could be some subset as in Fig. 9 or a more complicated function . One such function pair we suggest is the radius u = (x'-x)t't , x = (x„ x,,-, . . .,x„), and the inner product v = , where t = r(x, t) is the unit tangent to the homotopy path and t,, is the unit vector along the t-axis . A



J. D.

84

SEADER et a),

(b) Left eye

to) Right eye heft panel)

fright

panel)

Ttet

200

Pig . 9 . Stereoscopic display of hurnotopy path for CSTR reactor problem (viewed by crossing eyes),

further suggestion for r is the local signed curvature with the sign determined by the parity of the curvature vector components (i .e, the product of the signs of the components) . A change in parity may occur if the step size is too large and a new branch is leachedCONCLUSIONS A global fixed-point homotopy has been applied to find front a single, arbitrary starting guess all solutions to several sets of nonlinear equations, even when transcendental terms are present . However, depending upon the starting guess selected, the homotopy path may consist of branches that are only connected at infinities, which traverse the dependent variables and/or the homotopy parameter from --x to -I- . By using mapping functions, the fixed-point homotopy path may be transformed into a finite domain, wherein all solutions lie on the path . Two mapping functions have been developed, namely a toroidal mapping and a boomerang mapping . Both have been applied successfully to a number of applications, including an adiabatic CSTR operating in a steady-state mode with two consecutive complex reactions taking place . The use of mapped continuation should make possible the solution of a number of modeling problems which may not have been previously amenable to solution, particularly if all solutions are desired The method is an alternative to the method of Kuno and Seader (1988) . which requires the application of a starting point criterion . Acknowledgements-We appreciate the initial programming

of the toroidal mapping and the predictor-corrector steps by Mr Ming Tar Tsai . The computations for the two-equation forms of the CSTR problem were carried out by Mr Jui-Jung Chen . Partial financial support for this study was provided by the Phillips Petrnlcum Fomulation, Inc .

NOMENCLATURE A = Constant i n equation (6) ; constant in equations (42), (43), (48) and (49) B-Constant in equation (6) C; = Concentration of species i F= Feed flow rule f= Function vector or mapping J Derivative of function vector with respect to vector of independent variables g=An arbitrary function H - Homotopy function h = Step length k = Iteration index k„ k ., k~ - Reaction rate constants K, =Adsorption equilibrium constant n=Index, exponent in equation (6) p=Constant in equation (6) q=Constant in equation (6) R"= Vector space of n-tuples of real numbers r, =Reaction rate for the nth reaction s =Arclcngth 7 _=Temperature 1=Fixed-point homotopy parameter u = Complex extension vector s-Vector of independent variables ; an independent variable V-Vapor flow rate u=Complex extension vector x = Approximation to x y = Vecto, of independent variables Greek

Parameter defined below equation (20) Ii=Parameter in expressions below equation (23) -= Finite differential s=Element of A=Parameter defined by equation (21) n - Given value in equation (22) .=Parameter below equation (23) 0 = Time ps" =Estimated safe radius of the ball at step k for den Heifer and Rheinbotdt step size algorithm Parameter defined below equation (23) so = Parameter defined below equation (20) V = For all



Mapped continuation methods Subscripts and superscripts E = Eider prediction T = Transpose = Mapped space for I and x n = Index k - Index REFERENCES Abbott J . P ., Algorithm 110 . Comput . J . 23, 85-89 (1980) . Allgower E . L . and K . Gcorg, Relationships between deflation and global methods in the problem of approximating additional zeros of a system of nonlinear equalions . Free . of the NATO Advanced Research Institute on Homotopy Methods and Global Convergence (C . B . Evans et al ., Eds) . Plenum, New York (1983) . Bilous O . and N . R . Amundson, Chemical reactor stability and sensitivity_ AIChE Jt 1, 513-521 (1955) . Brent A . P., An algorithm with guaranteed convergence for finding a zero of a function . The Comput . J 14, 422-425 (1971) . Chavez C . R ., J. D . Seader and T . L . Wayburn, Multiple steady-state solutions for interlinked separation systems . Ind. Engng Chem . Fundam . 25, 566 (1986) . Diener I ., On the global convergence of path-following methods to determine all solutions to a system of nonlinear equations . Math . Program . 39, 181 (1987) . Dongarra J . J ., C- B . Molcr, J . R . Bunch and C . W . Stewart, LINPACK U.eers' Guide. SIAM, Philadelphia (1979) . Eisenstat S ., M . C . Gursky, M . Schutz and T_ Sherman, Yale sparse matrix package, II . Nonsymmetric codes . Report 114, Department of Computer Science, Yale University (1977) . Folger H . S ., Elements of Chemical Reaction Engineering . Prentice-Hall, Englewood Cliffs, New Jersey (1986) . He .. A . C ., REDUCE t.ser' .s Manual The Rand Corporation, Santa Monica (1987) . den Heijer C . and W . C . Rheinboldt, On steplength algorithms for a class of continuation methods . SIAM J. Numer . Analysis 18, 925 (1981) .

85

Hoenders B. J . and C . H . Slump, On the calculation of the exact number of zeros of a set of equations . Computing 30, 137 (1983) . Kubicek M ., Dependence of solution of nonlinear systems on a parameter. ACM Trans . Math . Software 2,98 (1976) . Kuno M . and J . D . Seader, Computing all real solutions to systems of nonlinear equations with a global fixed-point homotopy . Ind. Engng Chem . Res . 27, 1320-1329 (1988) . Lip W . J ., J . D . Seader and T . L . Wayburn, Computing multiple solutions to systems of interlinked separation columns. AIChE Jl 33, 886 (1987) . Morgan A . P ., Soloing Polynomial Systems using Continuation for Engineering and Scientific Problems . PrenticeHall, Englewood Cliffs, New Jersey (1987) . Muller D . E ., A method for solving algebraic equations using an automatic computer . Math Tables Aids Comput . (MTAC) 10, 208-215 (1956) . Rheinboldt W . C ., Solution field of nonlinear equations and continuation methods. SIAM J. Numer. Analysis 17, 221 (1980) . Rheinboldt W . C. . Numerical Analysis of Parametrized Nonlinear Equations . Wiley, New York (1986)_ Seader J . D ., Computer modeling of chemical processes. .47ChE Monograph Ser. 81(15) . AIChE, New York (1985) . Shacham M ., Comparing software for the solution of systems of nonlinear algebraic equations arising in chemical engineering. Comput . chem . Engng 9, 103-112 (1985). Smith J . M ., Chemical Engineering Kinetics, Third Edn, p . 426 . McGraw-Hill, New York (1981) . Watson L . T . and D . Fenner, Algorithm 555, Chow Yorke algorithm for fixed points or zeros of C 2 maps . ACM Trans . Math . Software 6, 252 (1980) . Watson L . T ., S . C . Billups and A . P . Morgan, Algorithm 652, HOMPACK : a suite of codes for globally convergent homotopy algorithms . ACM Trans . Math . Software 13, 281 (1987) . Wayburn T . L . and J . D . Sender, Homotopy continuation methods for computer aided process design . Comput . chem . Engng 11, 7-25 (1987) .