SOLVING GALERKIN APPROXIMATIONS TO NONLINEAR TWO-POINT BOUNDARY VALUE PROBLEMS BY A GLOBALLY CONVERGENT HOMOTOPY METHOD Layne T. Watson, University of Michigan L. Ridgway Scott, University of Michigan Technical Report 86-17 April 28, 1986

SOLVING GALERKIN APPROXIMATIONS TO NONLINEAR TWO-POINT BOUNDARY VALUE PROBLEMS BY A GLOBALLY CONVERGENT HOMOTOPY METHOD LAYNE T. WATSON+ and L. RIDGWAY SCOTTS Abstract. The Chow-Yorke algorithm is a homotopy method that has been proved globally convergent for Brouwer fixed point problems, certain classes of zero finding and nonlinear programming problems, and two-point boundary value approximations based on shooting, finite differences, and spline collocation. The method is numerically stable and has been successfully applied to a wide range of practical engineering problems. Here the Chow-Yorke algorithm is proved globally convergent for a class of Galerkin approximations to nonlinear two-point boundary value problems. Several numerical implementations of the algorithm are briefly described, and computational results are presented for a fairly difficult magneto-hydrodynamics boundary value problem. Key words. homotopy method, Chow-Yorke algorithm, globally convergent, two-point boundary value problem, Galerkin method, finite element method, nonlinear equations 1. Introduction. The Chow-Yorke algorithm (actually a family of algorithms) for nonlinear systems of equations was proposed in 1976, and since that time both the theory and the scope of its practical applicability have been greatly extended. This homotopy algorithm is accurately described as a globally convergent, probability-one algorithm. It is truly globally convergent in the sense that it will converge to a solution of the problem from an arbitrary starting point. The phrase "probability one refers to the rigorous theoretical results which guarantee convergence for almost all choices of some parameter vector, i.e., with probability one. Homotopy methods (both continuous [lj and simplicial [10, 24 ) were originally very inefficient, and dismissed by some as inherently inferior to quasi-Newton algorithms. A common misconception is that homotopy algorithms are just continuation or an obvious extension of classical continuation. While this is superficially true, there are fundamental philosophical differences between these probability-one algorithms and standard continuation. These differences result in an elegant and rigorous theoretical foundation for globally convergent probability-one homotopy algorithms, and have subtle vet important implications for computer implementations [56;. Furthermore, there has been a series of practical engineering problems, successfully solved by homotopy methods, on *which continuation and quasi-Newton methods either totally failed or experienced great difficulty'28-46-. Current implementations (as in HOMNPACK 155]) of these globally convergent probability-one homotopy algorithms are reasonably efficient. and their robustness, stability, and accuracy have never been in doubt. There are three distinct, but interrelated, aspects of homotopy methods: 1) construction of the right homotopy map, 2) theoretical proof of global convergence for this homotopy map. and 3) tracking the zero curve of this homotopy map. The first aspect is currently still an art, although this is much better understood now due to the accumulation of computational experience 28-46, 49. Although much remains to be done, significant progress has been made on the second aspect. Global convergence has been proved for Brouwer fixed point problems'4. 48`. certain classes of zero finding 48, and nonlinear programming (both unconstrained and constrained) problems 49], and tNvo-point boundary value approximations based on shooting j511. finite differences 152], and spline collocation i54,. Recently Morgan [21-221 obtained some elegant results for polynomial systems. and the present work considers Galerkin approximations to nonlinear two-point boundary value problems. Various curve tracking algorithms have been around for a long time (e.g., r14-16 ). but t Departments of Electrical Engineering and Computer Science, Industrial and Operations Engineering, and Mathematics, University of Michigan, Ann Arbor, hMI 48104. The work of this author was supported in part by AFOSR Grant 85-0250. X Department of Mathematics, University of Michigan, Ann Arbor, MII 48104. 1

there is an important distinction between general curve tracking algorithms and homotopy curve tracking algorithms. The object of the former is the zero curve itself, whereas the object of the latter is a point at A = 1. This difference was emphasized in [47], 48], and [56], and the algorithms described here are more in the spirit of homotopy methods than a general tracker. Section 2 outlines the theoretical foundation of globally convergent, probability-one homotopy methods, and proves some convergence theorems for Galerkin approximations to nonlinear twopoint boundary value problems. Section 3 discusses several algorithms for tracking the homotopy curves, as well as some pertinent details of the software package HOMPACK used to obtain the numerical results. Section 4 considers a nontrivial two-point boundary value problem from fluid dynamics, presents numerical results, and illustrates the dependence of the solution on the problem parameters. 2. Theory. Consider the two-point boundary value problem y"(x) = f(x, y(x), y'(x)), O < x 1, (1) y(O) = y(1) = 0, (2) where y(x) is an n-dimensional vector function and f(z, u, v) is an h-dimensional C2 vector function. This is naturally posed in a weak form via a variational formulation using the space V- {f E H1[0, 1] f(O) = f(1) 0}, where HlO, 1] is the Hilbert space of all absolutely continuous functions on [0, 1] whose derivative lies in L"2O, 11. Let Sn C V be a finite-dimensional vector space with basis 61,...,fn Using the finite-dimensional approximations ym(x) V x Am(x) = cmiOi(Z), m =1,..., n, (3) = 1 the standard Galerkin approximation to equations (1-2) is the nonlinear system of equations n E <mi, }ok { (( fm(, AA'), k) =0, i=1 k=,...,n, m=1,...,h. (4) Here (u, v = v u(t)v(t) dt (5) JO is the inner product. Example 1. Continuous piecewise linear functions. Let 0 = xo < xl < "' < < Xn, = 1 be a partition of [0,1', and let Sn be the space of continuous, piecewise linear functions f(x) with breakpoints X,..., x\+l and f(O) = f(1) = 0. The 2

basis functions pi, i= 1,.... n, are defined by 2(zxj) = ijy. Observe that in this case the equations (4) correspond to the difference equations 2 Am(Xk-l) - Am(Xk) Am(Xk) - Am(Xk-i) hk - hk.- 1 hk- 1 hk hk2 h (fm(, A, A'), k= fm(xk, A(Zk), A'(Xk)) +:(h) hk - hk-1 where hi = xi- zx^l, h = max3 hj. This is a standard finite difference approximation to (1) (cf. 752] for the case of a uniform mesh). Example 2. Chebyshev series (spectral methods). Using the Chebyshev polynomials {T,(x)}~o for the interval [-1, 1], take f Tk(2x - 1)- 1, k even dtk-l(z) =- Tk(2 r- 1)- (2x- 1), kodd for k = 2,3.... These bk(x) satisfy the zero boundary conditions, are dense in V, and are relatively well conditioned. The equations (4) in this case are not sparse. See Gottlieb and Orszag I12' for details. A third example, generalizing Example 1, would be a space Sn C Ck-2 0,lj consisting of piecewise polynomials of order k > 2 with B-splines for the basis functions oi. This spline space is developed in Section 4. The accuracy and convergence of the Galerkin approximation (3) to the exact solution o(x) of (1-2) can be shown under various assumptions on the nonlinear right hand side f(x, u. v) and generalizations of y" to elliptic operators L w(x) = Z;o(-l)- ( D(pj(x)liw(z)). To illustrate, some typical theorems from Ciarlet. Schultz, and Varga [5-91 will be quoted here. Following [5., assume that h = 1, a classical solution 6(x) of (1-2) exists, the nonlinear right hand side has the form f(x. y) (i.e., does not involve y'), and f(Z_ U) (w^(t))2dt f?- >'> -A - - inf (6) au we v;o ( w( t(t)) ~dt Jo for all x z 0, 1 and all real u. The appropriate norm (derived from an inner product based on a linearization of (4)) in this context is A. i 1 = [/ (w'(x))2 - (w()) dx], V, where the constant, is the same as that in (6). The accuracy of the Galerkin approximation.4(x) is given by Theorem. Let c(z) be the solution of (1-2) subject to the assumptions above, let S5 be any finite dimensional subspace of 1 of dimension n with basis 9i,....On. and let An(x) = f^l? Qix(r), 3

where an is the (unique under the above assumptions) solution to (4). Then there exist constants C and K independent of n and S, such that 11An- O|oo < KllAn - 4>1 < C inf ||w - 4\1. wESn Theorem. Let 0(x) be the solution of (1-2) subject to the assumptions above, let {Sn}n=1 be any sequence of finite dimensional subspaces of V such that U 1 Sn is dense in V in the norm 1[ * [|, and let {An(x)}~~ be the sequence of functions obtained by solving (4) over the subspaces Sn respectively. Then {A,(x)})l converges uniformly to {(x). Generalizations of these theorems [5-9, 26] tend toward generalizing the differential operator on the left hand side in (1), weakening the smoothness assumptions on the right hand side in (1), and generalizing the boundary conditions (2), rather than the form of and growth conditions on f(x, u, v). The homotopy theorems derived here require C2 smoothness (piecewise C2 is the absolute minimal requirement), and have rather mild growth restrictions on f(x, u, v) compared to (6). Furthermore the f considered here can involve y'(x) in a highly nonlinear way, whereas the f = f(x, u) in (6) cannot involve y'(x) at all. Roughly speaking, with the exception of smoothness, the limitations of the theoretical foundations for applying globally convergent homotopy methods to Galerkin approximations to nonlinear two-point boundary value problems derive more from the applicability of the Galerkin method itself than from the convergence conditions of the homotopy methods. Let Y = (tllCl2,... ln, a21i,...2n,... Cl... aCn) Then the system of equations (4) has the form F(Y) = M Y + N(Y) =0, (7) where /3M 0... 0 0 0... M, N(Y}) ((.fl(xA(z),A'(z)), ~I(X)),... (fl(XA(x),At(z)), >n(x)), (f2(XZ A(x) A'(x))} <^i(z)). (ffi2(x, A(x)} A(x))} on(z)); N( Y (x A(x), A'(x)), 4i(x)),..., (fn(x, A(x), A'(x)), n(x))). Thus the two-point boundary value problem (1-2) is approximated by the nonlinear system of equations (7), which has dimension p = hn. A homotopy method is used to solve (7). Since (utv) -* (U,; v1) is an inner product on V and dl,...,>nd are linearly independent, the Gram matrix Al is symmetric and positive definite. Therefore M is also symmetric, positive definite. 4

Let EP denote p-dimensional real Euclidean space. The following four lemmas from [511,!51], 50], [521 respectively will be useful. Lemma 1. Let F: EP EP be a C2 map, a E EP, and define pa: [O, 1) x El' EP by pa(A, y) = AF(y) + (1 - A)(y - a). Then for almost all a E EP there is a zero curve - of pa emanating from (0, a) along which the Jacobian matrix Dpa(A, y) has full rank. Lemma 2. If the zero curve I in Lemma 1 is bounded, it has an accumulation point (1. y), where F(y) = 0. Furthermore, if DF(y) is nonsingular, then 7 has finite arc length. Lemma 3. Let F: EP -- EP be a C2 map such that for some r > zxF(xz) > 0 whenever ||zxi = r. Then F has a zero in {z EP " xi! < r}, and for almost all a E E ", all < r, there is a zero curve a of pa(A, x) = AF(x) - (1 - A)(x- a), along which the Jacobian matrix Dp,(A, x) has full rank, emanating from (0, a) and reaching a zero? of F at A = 1. Furthermore, - has finite arc length if DF(x) is nonsingular. Lemma 3 is a special case of the following more general lemma. Lemma 4. Let F: EP - EP be a C2 map such that for some r > 0 and r > 0, F(x) and x- a do not point in opposite directions (pF(x) + (x - a) 0 Vp > 0) for xi = r,''all < i. Then F has a zero in {z E' F x' < r}, and for almost all a E EP, Iai <, there is a zero curve' of pa(A, x) = AF(x) (1 - A)(z- a), along which the Jacobian matrix Dp,(A, x) has full rank, emanating from (0, a) and reaching a zero X of F at A = 1. Furthermore,. has finite arc length if DF(x) is nonsingular. Theorem 1. Let ( V) in (7) be a C2 mapping, and suppose there exist constants C and v such that limsup (Y 2 C, < < 1. (8) For 1W EP, define p r: [0, 1) x EP EP by P H(AX, Y) = A F( r) (1 - )(Y - w). Then for almost all W E EP there exists a zero curve - of pwi, along which the Jacobian matrix Dpw(A, Y) has full rank, emanating from (0, W) and reaching a zero Y of F (at A = 1). Furthermore, if DF(Y) is nonsingular, then' has finite arc length. Proof. As observed above, the matrix M in (7) is symmetric and positive definite. Therefore all its eigenvalues are real and positive. Let F > 0 be the smallest eigenvalue of A.. Then a simple argument shows that for Y # O, YtM Y> F i'llYl > o. (9) Choose e > O. Using (8) and (9) yields y'F(y) = Y-tA;Iy - yt,(Y) >r i -r Y (c- ~) r'" > 0 5

for 11 Y \2 sufficiently large. Therefore there exists r > 0 such that YtF(Y) > 0 for 1| Y\\2 = r. The result now follows directly from Lemma 3. Q.E.D Theorem 2. Let N(Y) in (7) be a C2 mapping and F > 0 defined by (9). If limsup tY112 C < r, then the conclusion of Theorem 1 holds. Proof. This follows immediately from the proof of Theorem 1. Q.E.D. Observe that Theorem 2 is stronger than Theorem 1, but since Theorem 2 requires knowledge of r, Theorem 1 is easier to verify. The next two theorems, giving conditions on the right hand side f(x, u, v) in (1), are what one would use in practice. Theorem S. If f(x, u, v) in (1) is C2 and bounded, then the conclusion of Theorem 1 holds. Proof. f(z, u, v) bounded and C2 implies that N(Y) is also bounded and C2, and thus (8) holds with C = 0 and v = 1/2. Therefore the conclusion of Theorem 1 holds. Q.E.D. Theorem 4. Let f(z, u, v) in (1) be a C2 mapping, and suppose there exist constants u, M, and v such that l f(z, xu, ) 112 (< + (|11u|l2 + IlV 2)'), 0 < V < 1, for all 0 < x < 1 and u, v E En. Then the conclusion of Theorem 1 holds. Proof. Let n Yk() )= Y(k-l)n++ii(), k = 1,...,, n 0 < x < 1, (10a) i=l1 and n Zk(x) = Z(kl)n+i i(x), k = 1,...,, 0 < x < 1. (10b) i= 1 Since the basis functions c6 are in V, it follows from the equivalence of norms in finite dimensional spaces that there exist constants K1 > 0 and K2 < oo such that K(oatC )l/2 < < K2(ctc)'12, (11) i i=1 I UW [O,1] where a = (cl,..., n)t is a real n-vector and Wp[O, 1] is the Sobolev space of functions whose (weak) derivatives lie in LP[O, 1], and the Sobolev norms are (cf. [26]) rr 1 -1I/p V,!il] = -.i 1 tv(x)p p dx 1< < 0o, i VI 1!= <sup< v ( ), v) Z - Y n^^ 0~z< y~l I;XZ- yI) For instance, the constant K2 corresponding to Example 1 in Section 2 has the form K2 - (constant)/ /mmj hj. 6

Now using (10), (11), the theorem hypothesis, the Cauchy-Schwarz inequality, Jensen's inequality, and the relation (a + b)2 < 2(a2 + b2) gives i(Y)!2 = sup (y)tZ- SUp (fm(- YY,), Zm()) 1IZli2 =lc Zl 2 = 1 m= 1 n < Sup I lim' Y y)i/2 liZm( )!2. lZill2=l m=l < sup llf(., Y,')l z(.)11|z 11 Z12=1 r1 1/2 < sup f[ fl(x, y(x), y'(x))2 + + ft(x, y(x), y'())2 dx K2 1 Z2'ZI! 2=1lO LJ -'-K [p2 ~2 1/ v I2 d _5 y2 2|1 (Iy(x)\\ iy'(x)2> dz,u<K 2'S:;-I i:,t 1 dx y,"! =y ( )l;* ( / ) 1 ] )/2 -h' (,t (- t )V/2 h- i; \,y 2v) (rf) Y(X) 1 2 Y, i! 2 LI x 2V. Q.E.D. Even though the growth conditions on f(z, u, v) in Theorems 3 or 4 do not hold. it may turn out that for some constant K > O, the solution to (1) is the sanme as the solution to (1) with f(. u, ) replaced by A( ffuv) if Iu2 < h'and t2 <Kh. f(xu) l(x'(u).(v)) if u > >'. C2 2P-",2 -.;12) ( ) a globalrly the convergent homotop algorithm can still be successfull applied to the problem. that for some constant K > 0, the solution to (1) is the same as the solution to (1) with. u, )

3. Algorithm. The general idea of the algorithm is apparent from Theorem 1: just follow the zero curve - emanating from (0, W) until a zero Y of F(Y) is reached (at A = 1). Of course it is nontrivial to develop a viable numerical algorithm based on that idea, but at least conceptually, the algorithm for solving the nonlinear system of equations (7) is clear and simple. The homotopy map is pw(A, Y) = AF(Y) + (1 - A)(Y - W), (13) which has the same form as a standard continuation or embedding mapping. However, there are two crucial differences. In standard continuation, the embedding parameter A increases monotonically from 0 to 1 as the trivial problem Y - W = 0 is continuously deformed to the problem F(Y) = 0. The present homotopy method permits A to both increase and decrease along - with no adverse effect; that is, turning points present no special difficulty. Of course, methods have been developed [14-19, 23] to navigate turning points, but the method considered here does so automatically. The second important difference is that there are never any "singular points" which afflict standard continuation methods. The way in which the zero curve - of pw is followed and the full rank of Dpw along?- guarantee this. Observe that Lemma 1 guarantees that y cannot just "stop" at an interior point of [0, 1) x EP. The zero curve y of the homotopy map pw(A, Y) in (13) can be tracked by many different techniques; refer to the excellent survey [l] and recent work by Rheinboldt and Burkhardt 1231 and MIejia 1181. The numerical results in Section 4 were obtained with the software package HOMPACK, currently under development at Sandia National Laboratories, General Motors Research Laboratories, Virginia Polytechnic Institute and State University, and the University of'Michigan. HO.MPACK is a suite of codes for tracking zero curves of probability one homotopy maps, and provides both high-level and low-level subroutines for three different approaches to tracking. The three algorithmic approaches provided by HOMPACK are: 1) an ODE-based algorithm derived from that in 47l, with several refinements; 2) a predictor-corrector algorithm whose corrector follows the flowA normal to the Davidenko flow (a "normal flow" algorithm); 3) a version of Rheinboldt's linear predictor, quasi-Newton corrector algorithm'23' (an "augmented Jacobian" method). There are qualitatively different algorithms for dense and sparse Jacobian matrices; only algorithms for dense Jacobian matrices are discussed here. See [56] for a discussion of sparsity in relation to homotopy methods. First the ODE-based algorithm will be discussed. Assuming that F(]Y) is C2 and W is such that Theorem 1 holds, the zero curve -7 is C1 and can be parametrized by arc length s. Thus A = A(s). Y = Y(s) along y, and p (X(s), Y (s)) =0 (14) identically in s. Therefore d (( YA a -P M(A(s), Y(s)) = Dp (XA(s), Y(s)) I 0 (15) ds ds d T ds d)s 2)' (6) 8

If we take A(0) = 0, Y (0)= W. (17) the zero curve 7 is the trajectory of the initial value problem (15-17). When A(s) = 1, the corresponding Y(s) is a zero of F(Y). Thus all the sophisticated ODE techniques currently available can be brought to bear on the problem of tracking - [271, [48]. ODE software requires (dA/d, d dY/'ds) explicitly, and (15). (16) only implicitly define the derivative (dAXids, d Y/lds). This can be calculated by finding the kernel of the p x (p 1) Jacobian matrix DpW (A(s), Y(s)), which has full rank by Theorem 1. It is here that a substantial amount of computation is incurred, and it is imperative that the number of derivative evaluations be kept small. Once the kernel has been calculated, the derivative (dA//ds, dY/'ds) is uniquely determined by (16) and continuity. Complete details for solving the initial value problem (15-17) and obtaining Y(s) are in [47] and 51]. Remember that tracking A- was merely a means to an end, namely a zero Y of F(Y). Since' itself is of no interest (usually), one should not waste computational effort following it too closely. However, since'- is the only sure way to Y, losing - can be disastrous. The tradeoff between computational efficiency and reliability is very delicate, and a fool-proof strategy appears difficult to achieve. This is the reason HOMIPACK provides several algorithms: no single algorithm is superior overall. and each of the three beats the other two (sometimes by an order of magnitude) on particular problems. The normal flow algorithm due to Georg'11i has three phases: prediction. correction, and step size estimation. (13) and (14) are the relevant equations here. For the prediction phase, assume that several points p(1) (A(s1), Y(s1)), p( = (A(s2), Y(s2)) on; with corresponding tangent vectors (dA ds(s1), dY ds(si)). (dA ds(s2), dY'ids(s2)) have been found, and h is an estimate of the optimal step (in arc length) to take along,. The prediction of the next point on 3 is (~) = p(s2 t A), (18) where p(s) is the Hermite cubic interpolating (A(s), Y(s)) at s1 and s2. Precisely, p(si) (A(sl), 1 (si)),'(s1) = (dAX/ds(sl), dY! ds(s8)), p(s) = (A(s-), Y(st)), p'(s2) = (dA> ds(s2), d )'ds(s)), and each component of p(s) is a polynomial in s of degree less than or equal to 3. Starting at the predicted point Z(), the corrector iteration is z( =Z - [Dp-(Z()] pwZ( ), n=O. 1.... (19) where Dp~ -(Z(~))' is the M\loore-Penrose pseudoinverse of the p x (p -1) Jacobian matrix Dpe. Small perturbations of IW produce small changes in the trajectory -. and the family of trajectories - for varying IW is known as the "Davidenko flow". Geometrically, the iterates given by (19) 9

return to the zero curve along the flow normal to the Davidenko flow, hence the name "normal flow algorithm". A corrector step AZ is the unique minimum norm solution of the equation [Dpw] AZ = -p w. (20) Fortunately AZ can be calculated at the same time as the kernel of [Dpw], and with just a little more work. Normally for dense problems the kernel of [Dpug] is found by computing a QR factorization of [Dpw], and then using back substitution. By applying this QR factorization to -pw and using back substitution again, a particular solution v to (20) can be found. Let u - 0 be any vector in the kernel of [Dpw]. Then the minimum norm solution of (20) is v t AZ= v- t u. (21) utu Since the kernel of [Dpw]l is needed anyway for the tangent vectors, solving (20) only requires another 0(p2) operations beyond those for the kernel. The number of iterations required for convergence of (19) should be kept small (say < 4) since QR factorizations of [Dpwj are expensive. The alternative of using [Dpw (Z(()); for several iterations, which results in linear convergence, is rarely cost effective. When the iteration (19) converges, the final iterate Z(n~1) is accepted as the next point on y and the tangent vector to the integral curve through Z() is used for the tangent-this saves a Jacobian matrix evaluation and factorization at Z(nl). The step-size estimation described next attempts to balance progress along y with the effort expended on the iteration (19). Define a contraction factor!Z(2) _ Z()i L Z(1) _ Z()' (22) a residual factor R iPw (Z(1)) I (3) a distance factor (Z' = lim-oo Z() 1Z(1) - z" D =- (24) and ideal values L, R, D for these three. Let h be the current step-size (the distance from Z' to the previous point found on oy ). and h the "optimal" step-size for the next step. The goal is to achieve L R D hi - " - ~ - (25) L R D hq for some q. This leads to the choice h = (min(L/L,./R, D/D}) 1/ h, (26) 10

a worst case choice. To prevent chattering and unreasonable values, constants hiin (minimum allowed step-size), hmax (maximum allowed step-size), B,,i,, (contraction factor), and Bm,, (expansion factor) are chosen, and h is taken as h min {max { hmin, Bminh, h}, Bmc h hmI }. (27) There are eight parameters in this process: L, R, D, hnin, hnx Bmin, BmL,c, q. HOMPACK permits the user to specify nondefault values for any of these. The choice of h from (27) can be refined further. If (19) converged in one iteration, then h should certainly not be smaller than h, hence set h:=max {h, h} (28) if (19) only required one iteration. To prevent divergence from the iteration (19), if (19) has not converged after K iterations, h is halved and a new prediction is computed. Every time h is halved the old value hIod is saved. Thus if (19) has failed to converge in K iterations sometime during this step, the new h should not be greater than the value hold known to produce failure. Hence in this case i:= min { hold, (29) Finally, if (19) required the maximum K iterations, the step-size should not increase. so in this case set,:= min{h.h}. (30) The logic in (28-30) is rarely invoked, but it does have a stabilizing effect on the algorithm. Rheinboldt's augmented Jacobian algorithm together with step size strategies has been described very well elsewhere [2, 23, and will not be repeated here. 4. Numerical results. As a first example, consider Jeffery-Hamel flow [17X given bS y1 4y -4 6y2= c, y(O) = y(l) = 0, *where c is a constant. Using the piecewise linear basis functions {oi}^1 described in Example 1 in Section 2, the Galerkin approximation for y(x); Y(x) L= i i'i{(x) is the nonlinear system of equations I (s, 6k') - ( C - 4 - 6, ok = 0, k..., n. i=l The remark following Theorem 4 applies to this situation, and the homotopy algorithm can be applied with no difficulty. Rather than presenting numerical results for Jeffery-Hamel flow (which are well known, cf. i17), we instead tackle a problem where y is a vector, the boundary conditions are unbalanced (so integration by parts does not give (-A",Ok) = A, k), the interval is 0. xA). and shooting and finite difference methods are known to have difficulty!25i. 11

Consider the magneto-hydronamics (MHD) nonlinear two-point boundary value problem i25]: H"' = HH"- (H')2 + H 2 G (31)' = H G- - H' G- mG (32) 0" = Pr HO' (33) with the boundary conditions H(O) = -A, H'(0) = 0, G(0) = 1, 0(0) = 1, (34) H'(r]) - 0, G(r7) -O 0, (77) - 0 as rt - oo, (35) where Pr denotes the Prandtl number and m is a magnetic parameter. Note that H(r7) and G(r7) are determined by equations (31) and (32) independent of 0(7t), and thus 0(7r) can be uncoupled from the system. Once H(r) has been found from (31), (32). (34), and (35), 9(r7) can be determined from the one-dimensional two-point boundary value problem constituted by (33) and the 9 parts of (34), (35). The (routine) computation for 0(77) will not be discussed here. In practice the boundary conditions (35) are replaced by H'(I) G(r) = 0(r) = 0 (35a) for some r sufficiently large. Let Sn be the finite dimensional vector space with basis {Bj,k,t (x) n1 where Bjk,t (X) is the j-th B-spline of order k (degree < k - 1) defined on the knot sequence t = (t, t2,..., tnk). When there is no ambiguity Bjkt () is simply written as Bj(x). For this problem the knot sequence t is based on the breakpoint sequence = (0,.25,.50,.75, 1.0, 1.25, 1.50,1.75 2.0, 2.252.50,2.75, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0,8.0,9.0, 11.0, 13.0, 15.0, 18.0, 21.0, 24.0,28.0, 32.0, 36.0, 41.0,46.0, 51.0,60.0, 70.0, 80.0,90.0,100.0), following the convention tl = t2 * = tk and tn+1 = tn-2 =' * = tn+kDepending on the values of n and k. only an initial subsequence of < may be needed. The approximations are N+2 H() = - ajBj(rB), (36) j=l A B1(0) -a__B__(_-_ t =-A, a2 =- B,(0)' a2 B'2- (37) 12

Table 1. A m k Nr- - 2 ) H(r (c) NFE CPU time arc length -1.0 1.0 6 12 -.80700 -.43166 29 i3:28 1.309 -1.0 2.0 6 12 -.88173 -.78156 29 13:28 1.387 -1.0 4.0 6 12 -.94477 -.93015 30 14:37 1.518 0.0 1.0 6 12.11991.25331 25 11:13 1.067 0.0 2.0 6 12.07372.10858 25 11:24 1.103 0.0 4.0 6 12.03558.04078 20 9:07 1.198 1.0 1.0 6 12 1.06079 1.0898 21 9:38 1.141 1.0 2.0 6 12 1.03959 1.0481 20 8:57 1.187 1.0 4.0 6 12 1.02103 1.0225 15 6:42 1.272 2.0 1.0 6 12 2.02744 2.0318 18 8:00 1.257 2.0 2.0 6 12 2.01968 2.0213 24 10:38 1.296 2.0 4.0 6 12 2.01193 2.0123 15 6:37 1.365 4.0 1.0 6 12 4.00625 4.0064 18 7:58 1.471 4.0 2.0 6 12 4.00530 4.0054 18 7:58 1.491 4.0 4.0 6 12 4.00402 4.0041 15 6:38 1.530 -1.0 1.0 6 24 -.43877 -.43166 32 1:09:37 1.916 -1.0 2.0 6 24 -.78196 -.78156 27 58:20 1.604 -1.0 4.0 6 24 -.93019 -.93015 24 51:32 1.953 0.0 1.0 6 24.25286.25331 26 55:45 1.986 0.0 2.0 6 24.10852.10858 21 44:53 1.804 0.0 4.0 6 24.04073.04078 25 53:27 1.903 -1.0 1.0 6 32 -.43165 -.43166 33 2:15:48 2.765 -1.0 2.0 6 32 -.78158 -.78156 34 2:20:28 2.334 -1.0 4.0 6 32 -.93018 -.93015 27 1:52:50 2.752 -1.0 1.0 4 24 -.42854 -.43166 41 55:37 2.196 -1.0 2.0 4 24 -.78043 -,.78156 38 51:12 1.810 -1.0 4.0 4,24 -.93062 -.93015 29 39:06 2.128 N+2 G(r?7)- )=,Bj() (38) j=i 3 = 1, 3N-2 = 0. (39) The boundary conditions (34-35) force the equations (37) and (39). The Galerkin approximation is the nonlinear system of equations. -H"' - H" " H H" (H')2 m H'' G B,\ = O. i = 3,...., N- 1, (40) -G H - G-H - m G, B, =, i= 2.... N i- 1, where (u, t = J u(r7)z}(;) drj. Let V = (a3.a4,... a. 1aA,32,33,...,3N Yl)t and F(Y) = 0 be given by the p = 2N - 1 2n - 5 equations (40). Table 1 shows some numerical results obtained by applying HOMPACK to (40). The value: H(oc) are from'25i, and 7 can be inferred from n, k, and the breakpoint sequence l listed above. The integrals in (40) were computed by 10-point Gaussian quadrature over each subinterval. and are are thus essentially exact. NFE is the number of Jacobian matrix evaluations, and the format of the CPU time (on a VAX 11'780) is hh:mm:ss. The local curve tracking tolerance was 10-4 and the final accuracy (the end game tolerance [55j) was 10-s. Alt-hough the theory in Section 2 is not directly applicable to this MHD problem. these numerical results suggest that the homotopy algorithm is more widely applicable than the theory 13

indicates, and the accuracy is exactly what would be expected given the spline order and knot spacing. References. 1F E. Allgower and K. Georg, Simplicial and continuation methods for approximating fixed points, SIAM Rev., 22 (1980), pp. 28-85. I21 S C. Billups, An augmented Jacobian matrix algorithm for tracking homotopy zero curves, M.S. Thesis, Dept. of Computer Sci., VPI & SU, Blacksburg, VA, Sept., 1985. 3l P. Boggs, The solution of nonlinear systems of equations by A-stable integration techniques, SIAM J. Numer. Anal., 8 (1971), pp. 767-785. [41 S. N. Chow, J. Mallet-Paret, and J. A. Yorke, Finding zeros of maps: homotopy methods that are constructive with probability one, Math. Comp., 32 (1978), pp. 887-899. 51 P. G. Ciarlet, M. H. Schultz, and R. S. Varga, Numerical methods of high-order accuracy for nonlinear boundary value problems. I: One-dimensional problems. Numer. Math., 9 (1967), pp. 394-430. 6, P. G. Ciarlet,'M. H. Schultz, and R. S. Varga, Numerical methods of high-order accuracy for nonlinear boundary value problems. II: Nonlinear boundary conditions, Numer. Math., 11 (1968), pp. 331-345. 71 P. G. Ciarlet, MI. H. Schultz, and R. S. Varga, Numerical methods of high-order accuracy for nonlinear boundary value problems. III: Eigenvalue problems, Numer. Math., 12 (1968), pp. 120-133.,81 P. G. Ciarlet, M. H. Schultz, and R. S. Varga, Numerical methods of high-order accuracy for nonlinear boundary value problems. IV: Periodic boundary conditions, Numer. Math., 12 (1968). pp. 266-279. 9 P. G. Ciarlet, M. H. Schultz, and R. S. Varga, Numerical methods of high-order accuracy for nonlinear boundary value problems. 1: Monotone operators, Numer. Math., 13 (1969), pp. 51-77. 10' B. C. Eaves and R. Saigal, Homnotopies for computation of fixed points on unbounded regions, Math. Programming, 3 (1972), pp. 225-237. 11' K. Georg, private communication, 1981. 12i D. Gottlieb and S. A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, 1977. 13 H. B. Keller. Numerical Solution of Two-point Boundary Value Problems. Society for Industrial and Applied Mathematics, Philadelphia, 1976. 14!, Numerical solution of bifurcation and nonlinear eigenvalue problems, in Applications of Bifurcation Theory, Academic Press, New York, 1977. 151 R. W. Klopfenstein, Zeros of nonlinear functions, J. Assoc. Comput. MIech., 8 (1961), pp. 336-373.'16l M. Kubicek, Dependence of solutions of nonlinear systems on a parameter, ACM-TOMS. 2 (1976), pp. 98-107. 17 L. D. Landau and E. IM. Lifshitz, Fluid Mechanics, Pergamon, Oxford, 1959. 18; R. Mejia, CONKUB: A conversational path-follower for systems of nonlinear equations. J. Comput. Phys., to appear. 14

119] R. Menzel and H. Schwetlick, Zur Losung parameterabhdngiger nichtlinearer Gleichungen mit singularen Jacobi-Matrizen, Numer. Math.. 30 (1978), pp. 65-79. I201 G. Meyer, On solving nonlinear equations with a one-parameter operator embedding, SIAM J. Numer. Anal., 5 (1968), pp. 739-752. 211 A. P. Morgan, A transformation to avoid solutions at infinity for polynomial systems, Appl. Math. Comput., 18 (1986), pp. 77-86. 221], A homotopy for solving polynomial systems, Appl. Math. Comput., 18 (1986), pp. 87-92. 223j TW. C. Rheinboldt and J. V. Burkardt, Algorithm 596: A program for a locally parameterized continuation process, ACM Trans. Math. Software, 9 (1983), pp. 236-241. [24] R. Saigal and M. J. Todd, Efficient acceleration techniques for fixed point algorithms, SIAM. J. Numer. Anal., 15 (1978), pp. 997-1007. [L25 K. K. Sankara, W. I. Thacker, and L. T. Watson,,manuscript. [261 M. H. Schultz, Spline Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1973. {27j L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary Differential Equations: The Initial Value Problem, W. H. Freeman, San Francisco, 1975. 28 C.. Wang and L. T. Watson, Squeezing of a viscous fluid between elliptic plates, Appl. Sci. Res., 35 (1979), pp. 195-207. [29], Viscous flow between rotating discs with injection on the porous disc, Z. Angew. Math. Phys., 30 (1979), pp. 773-787. 30, On the large deformations of C-shaped springs, Internat. J. Mech. Sci.. 22 (1980), pp. 395-400. [31] _, Theory of the constant force spring. J. Appl. Mech., 47 (1980), pp. 956-958. 32] _, The fluid-filled cylindrical membrane container. J. Engrg. MIath., 15 (1981). pp. 81-88. [33, Equilibrium of heavy elastic cylindrical shells, J. Appl. Mlech., 48 (1981), pp. 582-586. L34 __, The elastic catenary, Internat. J. Mech. Sci., 24 (1982), pp. 349-357.:35, Free rotation of a circular ring about a diameter, Z. Angew. Math. Phys.. 34 (1983). pp. 13-24. 36, C. Y. Wang, L. T. W'atson, and M. P. Kamat, Buckling, postbuckling and the flow through a tethered elastic cylinder under external pressure, J. Appl. lech., 50 (1983), pp. 13-18. 37 L. T. Watson, T. Y. Li and C. Y. Wang, Fluid dynamics of the ellilptic porous slider. J. Appl. Mech., 45 (1978), pp. 435-436.;38 L. T. Watson and C. Y. W7ang, A homotopy method applied to elastica problems, Internat. J. Solids Structures, 17 (1981), pp. 29-37. [39; __. Deceleration of a rotating disc in a viscous fluid, Phys. Fluids, 22 (1979), pp. 2267-2269. r40, The circular leaf spring, Acta Mech., 40 (1981), pp. 25-32.:41 ___. Hanging an elastic ring, Internat. J. Mech. Sci., 23 (1981). pp. 161-168. 142], Overhang of a heavy elastic sheet, Z. Angew. Math. Phys., 33 (1982). pp. 17-23. 43 _, Periodically supported heavy elastic sheet, J. Engrg. Mech., 109 (1983). pp. 811-820. [441. The equilibrium states of a heavy rotating column, Internat. J. Solids Structures. 19 (1983), pp. 653-658. 15

UNIVERSITY OF MICHIGAN 1I IIIIIIIUIN 11 II IIII[11111111111111111111111 3 9015 03994 8032 [45] L. T. Watson and W. H. Yang, Optimal design by a homotopy method, Applicable Anal., 10 (1980), pp. 275-284. [46, Methods for optimal engineering design problems based on globally convergent methods, Comput. & Structures, 13 (1981), pp. 115-119. [47j L. T. Watson and D. Fenner, Chow- Yorke algorithm for fixed points or zeros of C2 maps, ACM TOMS, 6 (1980), pp. 252-260. [48] L. T. Watson, A globally convergent algorithm for computing fixed points of C2 maps, Appl. Math. Comput., 5 (1979), pp. 297-311. [49], Computational experience with the Chow- Yorke algorithm, Math. Programming, 19 (1980), pp. 92-101. 50], Solving the nonlinear complementarity problem by a homotopy method, SIAM J. Control Optim., 17 (1979), pp. 36-46. 51i], An algorithm that is globally convergent with probability one for a class of nonlinear two-point boundary value problems, SIAM J. Numer. Anal., 16 (1979), pp. 394-401. [52], Solving finite difference approximations to nonlinear two-point boundary value problems by a homotopy method, SIAM J. Sci. Stat. Comput., 1 (1980), pp. 467-480. [531, Engineering applications of the Chow-Yorke algorithm, Appl. Math. Comput., 9 (1981), pp. 111-133. [54j L. T. Watson and M. R. Scott, Solving spline collocation approximations to nonlinear two-point boundary value problems by a homotopy method, Tech. Rep. CS84015-R, Dept. of Computer Sci., VPI&SU, Blacksburg, VA, 1984. 55] L. T. Watson. S. C. Billups, and A. P. Morgan, HOAIPACK: A suite of codes for globally convergent homotopy algorithms, Tech. Rep. 85-34, Dept. of Industrial and Operations Eng., Univ. of Michigan, Ann Arbor, MI, 1985. [56 L. T. Watson, Numerical linear algebra aspects of globally convergent homotopy methods, Tech. Report TR-85-14, Dept. of Computer Sci., VPI&SU, Blacksburg, VA, 1985. 16