tHESiS LIBRARIES ax. r‘.v,.g'.s....re,rm STATE UNIVERW EAQF LANsme. MICH. 48824 This is to certify that the dissertation entitled HOMOTOPY CONTINUATION METHOD FOR NONLINEAR EQUATIONS presented by Mahmoud Mohseni Moghadam has been accepted towards fulfillment of the requirements for Ph . D . degree in Mathematics 4(L'fft— L,‘ Major professor MS U is an Affirmative Action/Equal Opportunity Institution 0-12771 MSU LIBRARIES “ RETURNING MATERIALS: Place in book drop to remove this checkout from your record. FINES will be charged if book is returned after the date stamped below. HOMOTOPY CONTINUATION METHOD FOR NONLINEAR EQUATIONS BY Mahmoud Mohseni Moghadam A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 1984 ABSTRACT HOMOTOPY CONTINUATION METHOD FOR NONLINEAR EQUATIONS BY Mahmoud Mohseni Moghadam The essence of the homotopy continuation method is path following. Chapter one of this dissertation is devoted to a detailed discussion concerning certain important aspects of the path following technique. It is shown that the determination of the orientation of the path is a by-product of this computation. In chapter two we used the homotopy continuation method to determine all roots of a symmetric polynomial system. Chapter three and four are contributing to the applications of the polynomial systems. In chapter three an algorithm is derived to approximate all real roots of an analytic function in a bounded domain. In chapter four we developed an algorithm to approximate all eigenvalues of a matrix. Finally in chapter five a special homotopy is constructed to show that there are exactly n distinct smooth curves connecting n trivial solutions to n eigenpairs of an n.xn given matrix. This homotopy may be used to approximate the eigenpairs of a matrix. DEDICATED TO MY PARENTS AND MY WIFE ii ACKNOWLEDGMENTS It is a great pleasure to express my sincere appreciation to my thesis advisor Professor Tien-Yien Li. His invaluable guidance, expert advice, stimulating discussion, and useful insights made this work possible. I would also like to thank Professor Richard Hill for careful reading of this thesis. Finally, I thank my wife Tahereh Ghasemi for her confidence, patience, and understanding displayed throughout the duration of my graduate study. iii TABLE OF CONTENTS Chapter 0. INTRODUCTION HOMOTOPY CONTINUATION METHOD 1.1 Introduction . . 1.2 Homotopy . . . . 1.3 Existance of a path 1.4 Movement along the path FINDING THE ZEROS OF AN ANALYTIC FUNCTION 2. Introduction and notations Method and homotopy Following the curve Generalization . . . . . . . . . NNNN {)1wa?“ Numerical results . . . . . . . REAL ZEROS OF A FUNCTION 3.1 Introduction . . . . . . . . . . . 3.2 Integration around the boundary of a strip . . . . . . . . 3.3- Error analysis Numerical results EIGENVALUE PROBLEM . . . . . . . . 4.1 Introduction . . . . . . . 4.2 Method and Theorems 4.3 Numerical results . iv 48 49 51 64 65 65 66 73 Chapter 5. EIGENPAIR PROBLEM 5.1 Introduction 5.2 The homotopy 5.3 Theorems. . . 5u4 The algorithm BIBLIOGRAPHY Page 81 81 82 86 92 94 CHAPTER ZERO INTRODUCTION In the past decade,considerab1e advances have been made in the problem of obtaining numerical solutions of systems of nonlinear equations. Two closely related methods have become available. (See [24,38] for general references.) The simplicial methods developed by C.E. Lemke and J.T. Howson [30], H. Scarf [41], B.C. Eaves [14.15.16], R. Saigal [40], and others was employed initially for finding the Brouwer fixed point. Their method is based on using a simplicial approximation of the maps as is used in the Sperner's Lemma Proof of the Brouwer fixed point Theorem. R.B. Kellogg, T.Y. Li and J.A. Yorke developed an alternative approach [26,27] for the numerical solution of the Brouwer fixed point Theorem. By a twist of a nonconstructive proof of M. Hirsch [23], they obtained a constructive proof. Given a smooth map F of a ball in 391 into itself they choose a point p on the boundary of the ball. For «in almost any choice of p there is a smooth curve ‘p which leads to a fixed point. The curve can be followed efficiently by computer. 2 The Brouwer fixed point theorem itself is not used extensively in applications, though similar degree theoretic results are. S.N. Chow, J. Mallet-Paret and J.A. Yorke [4] showed that by using elementary homotopy arguments, numerical methods become available for many problems. (While B.C. Eaves introduced the homotopy idea to the simplicial approach, O.H. Merrill [34] and independently H. Kuhn and J.G. Mackinnon [27] pioneered the sandwiCh approach.) Their idea is to take any existance proof based on degree theory and to convert it using elementary homotopy into a constructive, computer implementable method for finding solutions. It is shown [4] that the method can be expected to converge with probability one and even is some situations where elementary degree theory is not applicable. This method is more wildly known as "Homotopy continuation method". The essence of the homotOpy continuation method is path following. One starts a certain path with a solution which is easily solved, and follows the path until the desired solution is reached. The path is given by the integral of some differential equation under certain regularity assumption, the path is well behaved and can be followed successfully. In fact, simplicial methods and continuation methods are not unrelated as they may initially seem. Both are techniques for following a certain path and are both related to Newton's method. In chapter one of this thesis, we shall present a detailed discussion concerning certain important aspects of the path following technique. It is shown that the QR decomposition can be applied to the computation of the vector field of the differential equation such that numerical stability can be achieved. It is also shown that the determination of the "orientation" of the path is a by product of this compuation. The homotopy continuation method is applicable to a considerably large number classes of problems. Among others, its application to locating all the isolated roots of polynomial systems has been most effective and profound. Recently, considerable attention has been given to the problem of finding all solutions of systems of n poly- nomials in n unknowns by homotopy methods. This problem is important in application and no methods other than homotopy methods are available for finding all the isolated solutions of such systems. Let p1(zl,...,zn) = O (l) Pn(al,...,zn] = 0 be n polynomials in n unknowns and write 2 = (21,. ,zn), and p(z) = (p1(z).....pn(z))- (2) F.J3 Drexler[12,13] and C.B. Garcia & W.I. Zangwill [20,21] showed independently and almost simultaneously that it is possible to find all isolated solutions of p(z) = 0. More precisely, it is possible to find a homotopy n H: [0.1] xcn -v C (3) which starts from a trivial set of polynomials 01(21) = O Qn(zn) = 0 That is 0(2) = 0, where Q(z) = (01(2),...,Qn(z)), H(O,z) = 0(2), H(1,z) = p(z) , and has the property that for each isolated zero ‘3 of p there is a smooth curve I of zeros of H in [0,1] an leading from some zero of 'Q at t = O to a zero of H(1[§). In [5] S.N. Chow, J. Mallet-paret and J.A. Yorke basically proved Drexler's results but replaced his algebraic geometry arguments with the results involving the generalized Sard's Theorem on which their general homotopy method is based. Garcia and T.Y. Li [19] used those developments to prove a classical theorem of Noether and Van der Waerden [43] concerning the exact count of the number of solutions of a polynomial system and to give a generalization of the fundamental theorm of algebra. The proof is again free from the algebraic geometry argument used in [37] and is done in such a way that all the isolated solutions can be explicitly calculated. The homotopy used by S.N. Chow, J. Mallet- Paret and J.A. Yorke [5] has been simplified by T.Y. Li [32]. In chapter two of this dissertation, we study a special system of polynomials with strong symmetric properties. The polynomial system is defined by p(z) = (p1(z)....,pn(z)) (4) where z = (21,...,zn) and pi(z) = zi-tz:-+-..-+z; ; i = l,2,...,n. (5) We shall present a special homotopy to exploit the symmetries of this system in such a way that the compu- tations of the curve following can be made much simpler. ~We also discuss certain techniques to overcome the difficulties when multiple roots are presented in (5). Chapter three and four are contributing to the applications of the polynomial system (5). In chapter three an algorithm is derived to approximate all the real roots of an analytic function F(z) in a bounded 6 domain. Numerical results up to seven roots are presented to show the effectiveness of this algorithm. In chapter four, we develop an algorithm to approximate all the eigenvalues of an n.xn matrix. In fact, the polynomial system (5) represents the trace formulas for the ith power of a matrix A when zj's, j = 1,2,...,n, are the eigenvalues of A. In chapter five, we consider the problem of determining all the eigenpairs of a symmetric matrix A. A special homotopy is constructed to show that there are exactly n distinct smooth curves connecting trivial solutions to the desired eigenpairs. Incorporated with sparse matrix techniques, this method may be used to solve eigenvalue prdblem for large Sparse matrices. CHAPTER ONE HOMOTOPY CONTINUATION METHOD §(1.1) Introduction Let F be a smooth function from l‘Rn to IRn . In this chapter we shall consider the problem of finding the solution of F(x) = O, by homotopy continuation method. The term "continuation method" is derived from a class of numerical methods dating at least back to E. Lahaye [28,29], and also known as "embedding method". Detailed discussion of these methods can be found in articles by H. Wacker [44], and E. Allgower & K. Georg [2]. One starts with a trivial equation, one to which the solution is obvious and immediately known. Then the system is deformed continuously to F(x) = O. In general, the solution of the trivial system will prescribe, under this deformation, a smooth curve which is connected to the solution of F(x) = O. This curve can be characterized by a solution to an initial value problem ofaniordinary differential equation. Our discussion here is limited to following the curve. As a main result, we show that the determination of the orientation of the curve is a 8 'byproduct of the computation of the vector field of the ordinary differential equation. §(1.2) Homotopy In early 1960's, Davidenko [8], [9] introduced a method of solving F(x) = 0 where F is a smooth function from JRn to En . Let H :anx [0,1] 41Rn be defined as H(x,t) = (1 -t) (x -a) +tF(x) (1.2.1) with a E HUI given. It is clear that H(x,O) = x-a and H(x,l) = F(x). Suppose (A): the partial derivative of H with respect to x is always nonsingular Then by repeat application of the Implicit Function Theorem, there exists a curve x(t), as a function of t, such that H(x(t) ,t) = O . (1.2.2) We differentiate (1.2.2) with respect to t, to get the differential equation 1:2- -1 dt x t (1.2.3) x(O) ll DJ 9 Thus, finding a zero of F(x) is equivalent to solving the initial value problem (1.2.3), and finding its value at t = 1. The assumption (A) is rather strong, and hence the power of Davidenko's method is restricted. Let us consider the homotopy H : 1Rn len x (0,1) .. an (1.2.4) defined by H(x,a,t) = (1 -t) (x -a) +tF(x) with x 5 1R“ , a e 1Rn and t 6 (0,1). For a fixed a E 1Rn , define n Ha :13“ x (0.1) .. 1R (1.2.5) by Ha(x,t) = H(x,a,t). Given this homotopy system, we define -1 _ . _ Ha (O) — [(x.t) -Ha(x,t) — 0} §(l.3) Existance of a path We shall use the generalized Sard's Theorem to assure _1( a 0) contains a smooth that for almost every a 6 H31, H curve Ta which will lead from a trivial solution to the desired solution. 10 Definition (1.3.1) Let U be an. open set in Ian, and f :U 4 mm be a Cl map, then y 6 Hg“ is called regular value of f provided m -1 range Df(x) = Hi for all x e f (y) . where Df(x) denotes the n.xm. matrix of partial derivative of f. The following theorem may be found in [1]. Theorem (1.3.2) (Generalized Sard's Theorem) Let v g IRn, w 5 112‘“ be open and let G:VxW alRp be smooth. If 0 6 NJ) is a regular value for G. then for almost every a E V (in the sense of either Baire category or Lebesgue measure), 0 is a regular value for Ga(-) 5 G(a,-). For our homotopy defined in (1.2.4) we have the following. Lemma (1.3.3) n . For almost every a 6 n1 , zero is a regular value of Ha : an x (0,1) -D IRn where Ha is given in (1.2.5). 11 Proof: Consider the homotopy H : an len x (0,1) 4 JRn , with H(x,a,t) = (1 -t)(x-a)-ttF(x) Let (£3.33) 6 H-1(O), i.e. H(§,§,'t') = o. It is clear that Dal-I = -(1 -t)I where DaH is the partial derivative of H with respect to a, and I is the n.xn identity matrix. Since Range DH (33.3,?) 3 Range DaH(§,a—,t‘) = IRn , we conclude that DH(§}§,E) has rank n. Hence 0 is a regular value for H. Thus, by applying the generalized Sard's Theorem we have for almost every a E nf‘, O is a regular value of Ha' This completes the proof. The overall idea is to start from a trivial solution of Ha(-,O) at t = O, and follow the path generated in Ha(-,t) as t goes from zero to one. We hope the trivial solution deforms into the solution of the original system, and hence we would be able to follow the connected path from the trivial system to the solution of F(x) = 0. Of course this is quite an idealized process , and there are a number of difficulties. First of all, in general a path need not exist. Second if one exists, it might be very illebehaved. In other words the set 12 ((x.t) :x e m“. t e (0.1). Ha(X.t) = 0) may consist of different solutions, such as isolated points, self-convergings, bifurcations, endless spirals, closed orbit, and smooth paths. But we are interested only in smooth paths (Figure (1.1)) «we Figure (1.1) Let a be choosen so that O is a regular value for (because of Lemma (1.3.3) this can be done with probability one). Ha(x.t) Then repeated use of the Implicit Function Theorem implies that H;l(0) consists of one dimensional manifolds. Detailed discussion for the existance of paths 13 are given by Garcia and Zangwill [22], [20], [21], and Chow, Mallet-Paret and Yorke [4]. Let Fa be the 1(O) with a as one endpoint. Also let component of H; us assume this component is parameterized by s. For notational convenience we refer to Ha(x,t) by H(x,t). Therefore, ' H(x(s),t(s)) = 0 (1.3.6) x(0) = a Differentiation of H with respect to parameter 5 yields HX(X(S),t(S)) GE +Ht(x(5).t(S)) -é = 0 (1.3.7) x(0) = a Here Hx and Ht are respectively the partial derivative of H with respect to x and t. The ordinary differential equations(1.3.7) can be written in the following matrix form :2 [H H ] , = O x t t (1.3.8) x(O) a t(O) o The integral solution of this differential equation, namely (x(s),t(s)) is a simple curve starting from (a,O). In the next section we carefully examine the movement along this curve. 14 §(1.4) Movement along the path 1 We have seen that H; (0) consists of only arcs and closed curves. These curves are the solutions of the ordinary differential equations i [H H] . =o x t t (1.4.1) x(O) O t(O) 0 where Hx is an n.xn matrix, Ht is an n.xl matrix and - = ii for some parameter s. For the remainder of this chapter, we will let s be the arc length. Since 0 is a regular value of H, [Hx Ht] is of full rank. Hence kernel of [Hx Ht]. is one-dimensional, by above the vector [x,t]t lies in this kernel. 0 X Let A = [H H 1 and g = . . Then (1.4.1) X t t simplifies to Ag, = o a (1.4.2) y(0) = O with Hi H2 = 1. The equation Ay = 0 means that y is perpendicular to the row space of A. In order to see how y can be determined by A, we first give the following definition. 15 Definition (1.4.1) A Housholder transformation in HUI is a matrix of the form U = I-2vvt where v 6 n9“, vtv = 1, and I is the identity matrix in LRn len . It is easily verified that Housholder transformations satisfy the following properties: 1. Ut = U (symmetric) 2. U2 = I (involutary) 3. UtU = I (orthogonal) Finally the matrix U has the important property that given any two vectors of equal length, x and y, we can find a matrix U such that Ux = y. To this end we take =x- 0 .o V.U;T-§T|—2- (143) By using these properties, a sequence of Housholder transformations, P1,P2,...,Pn in EJHJ' can be constructed such that t ...p = Pn 1A R where R = (rij) is an (nntl) xn upper triangular matrix, that is, rij = O for i > j. AlSC> Pi's can 16 be choosen in such a way that rii > O for i = 1,2,...,n. For details see [40,41] . Set ) . (1.4.4) Then Q is an orthogonal matrix and At = QR . (1.4.5) Suppose at a certain point y (s) and s are known, hence A is known. The following lemma enables us to find y (s-tAs) for the next point and trace the path by an ordinary differential equation solver. Lemma (1 .4.2) Let q be the last column of the orthogonal n+1 matrix Q, then y = i qn+l O (1.4.5.1) Proof: Since A = OR, and Ay = 0, we have Rtoty = 0 Since matrix R has rank n, we get rii ¢ 0 1 = 142,...,I1 . (1.4.6) I Suppose Qty = (81,62,...,Bn+1)t, then we have l7 r11 51 = 12 51"r22 52 = r1n Bl +r2n BZ'+°'°'+rnn Bn = 0 Because of (1.4.6) this system implies Bl = 82 - = Bn = 0 Hence t° _ t Q y - (0:0: oan+1) So 1 _ t Y "' Q(OIOIo--IBn+1) Therefore y is a scalar multiple of the last column of Q. Since Hy ”2 = l, we get y - i qn+1 In order to determine the orientation of y , we give the following theorem which can be found in [18]. Theorem.(1.4.3) Let H:]Rn+l aRn be a C1 map, and let 2(5) = (21(5),...,zn(s)) be a C1 curve in an+1 satisfying Then either 18 Ii(z(s)) 0 sgn 2],:(5) = sgn det Hi(z(s)) (1.4.7) or sgn z{(s) = -sgn det Hi(z(s)) (1.4.7.1) dz. . for all s, where z£(s) = 7§%' and H1 is the Jacobian of H with 13h column deleted. Applying this theorem to our homotopy, we get either That is sgn t (s) = sgn det Hx(x(s) ,t(s)) (1.4.8) or sgn t (s) = -sgn det Hx(x(s),t(s)) (1.4.8.1) for all 5. However, at s = 0, we have sgn det HX = 1. We may assume t (0) > 0, therefore (1.4.8) holds for all s > o, If we know sgn t (s) the sign in (1.4.5.1) is determined, to determine sgn t (s) for any sh> 0 we prove the following proposition: Proposition (1.4.4) Let Q = (qij) be as in (1.4.4), then ' _ n sgn t (s) - (Al) sgn(qn+l'n+1) (1.4.9) Proof: n+1 . Let en+1 6 Hz be the (nntl)£h unit vector. l9 _ t en+1 — (O'O'oootl) then t t _ t t = [R Q en+1] By property of the Housholder transformations ° det P = (-1)n t _ det Q — det Pn det Pn 1 _1 .. Hence t _ t , t det[R Q en+l] — det Q det[A en+1] t H O = (-1)n det X = (-1)n det H: Ht 1 t On the other hand, since R is an upper triangular matrix w1th rii > 0, we have sgn det[R Qte ). (1.4.10) n+1] = sgn(qn+1,n+l Therefore sgn t = sgn det H: = (-l)n sgn(qn+1 n+1). (1.4.11) From the above discussion we see that in order to follow the curve I computation of the vector field al x(s) t(s) at s .can be summerized as follows: 20 We first compute Hx(x(s),t(s)) and Ht(x(s),t(s)). Then write t H At = x t Ht as a product of an (n-tl) x(nntl) orthogonal matrix Q = (qij) and an (ni-l) xn upper triangular matrix R. x(s) Then , is given by the last column of Q with a t(s) possible sign change, and the sign of this vector is given by (1.4.9). CHAPTER TWO FINDING THE ZEROS OF AN ANALYTIC FUNCTION §(2.1) Introduction and notations Let f be an entire function and R be a bounded domain in the complex plane with a closed boundary C. Let us assume that C does not pass through a zero of f(z). It is well known from the theory of complex variables that 1 kf'(z) n k -—?' z = Z) 2. (2.1.1) ZUI EC f(z) i=1 1 where zi(i = 1,2,...,n) are all the zeros of f(z) which lie in R, (a multiple zero is counted according to its multiplicity in its formula). By using the algorithm explained in chapter three of this dissertation or any other proper algorithm we may find an approximation for n k s = 23 z., k=0,1,2,.... (2.1.2) =1 Therefore in order to determine the zeros of f we may solve the following system of polynomials. 21 22 f = zl+zz+~~+zn sl 2 2 _ z1+zz+ +2 —32 n n n _ KZl--]-22-i----+zn-Sn Delves and Lyness [11,33] developed an algorithm for finding all zeros of f(z) inside R as follows: i. Let f(z) be an analytic function on and within a given contour C. 11. Let 21 2.. of f(z) which have so far been obtained. ,2 ..,zk be a list of the known zeros iii. Let m, k, and e be given constants. iv. Search routine. The main routine attempts to calculate the number of zeros of f(z) within C using trapezoidal role approximation to contour integral I S = -£- ( é—LEL dz 0 2ni C f(z) Since the exact result of S0 is known to be an integer, the accuracy required is low. We only need to determine unambiguously which integer is involved. There are three possible outcomes. These are: 23 a. It finds that f(z) becomes unduly small on the contour and takes this to imply there area zero of f(z) close to the contour. In this case the integration method, if continued, would converge slowly. The routine does not continue the integration, but returns control to the search routine, which in turn chooses a different contour. b. It finds a value so. On checking the list of known zeros it finds that q of these lie within the region R, and hence there are ‘-q unknown zeros within So R. If SO-q > M it returns control to the search routine. c. It proceeds as in (b) but finds SO-q g_M. It then evaluates the SO-q unknown zeros as follows. It evaluates approximations to the sums of power of zeros So _ N - SN - 1:31 zi , (N — o,1,....sO-q) (2.1.4) using trapezoidal approximation to the integral -__1_ Nfla). S — Zvi I 2 dz N C f(z) Since the locations of the known zeros zl,z2,...,zn are available the sums Sfi of the power of unknown zeros are '3' 82(3) N 32 N 1 = z.=S - z. , N=0, ,...,s -q. N i=q+1 1 hi i=1 1 0 Having these numbers, a polynomial of degree SO-q may be constructed as follows: 24 Define 01 = -(zl+-22+- i-zso_q) 02 = (2122+- 4-zs -q-l 28 -q) 0 0 (2.1.5) SO-q o = (-1) z z z . SO-q l 2 SO-q Then by the Newton's identities Sl-t-o1 = 0 Sz-tslol-+202 = 0 (2.1.6) S +S o+---+(s —q)o =0 SO-q SO-q—l 1 0 SO-q One may construct the polynomial S -q S -q-1 _ 0 0 ... (2.1.7) P(z) — 2 4-01 2 + i-os -q 0 which has 2i (i = 1,2,...,SO-q) as roots. This polynomial is solved by using the polynomial root finding subroutine. In the above algorithm M is generally choosen to be 5. Although this is a big restriction, it is necessary due to the highly sensitive dependence of high degree polynomial on its coefficients (cf: Wilkinson's example [45]) . In other words M cannot be choosen too large 25 in order to avoid the ill-conditioning of the polynomial (2.1.7). Professor Li [32] recently solved this problem by homotopy continuation method directly without forming the polynomial (2.1.7). Hence he removed the search routine as well as the step (b) from the procedure. Here we shall present a special homotopy to exploit the symmetries of this system in such a way that the computation of the curve can be made much simpler. That is we will utilize the strong symmetry structure of (2.1.3) to write the corresponding differential equations explicitly in such a way that it will be ready to be solved by any ordinary differential equations solver. We also generalize the problem in the sense that the assumption "f(z) must have single roots" is relaxed. Namely, in section (2.5) we discuss certain techniques to overcome the difficulties when multiple roots are presented in (2.1.3). §(2.2) Method and Homotopy Let Q :Cn 4 Cn be a polynomial system, and Q = (Ql""'Qn)' By a polynomial system, we mean each term of Qi is of the form C O O n azl 22 2 (2.2.1) where a is a complex number, zi a complex variable, and 26 ri a nonnegative integer. For each term (of the form (2.2.1)) 1n 01' con31der the sum rli-r24-o-o4-rn Let qi be the maximum such sum in Qi' We assume qi > 0, for all i. We call qi the degree of Qi' With these notations we state the following theorem. Theorem (2.2.1) n Let Q :Cn 4 C be a polynomial system with degree _ ._ o _ o o Qi - qi (1 - lozo-ooon) 0 Let Q (2) - (01(2) 002(z)loool 02(2)) where 02(2) consists of the terms in 01(2) with degree q.. If 00(2) has only the trivial 1 solution 2 = 0, then 0(2) = O has n n q. (2.2.1.1) solutions (counting the multiplicity). 2222:: See Garcia and Li (on the number of solution to polynomial system of equation, Theorem 3.1, p. 543, [19]). Their proof based on homotopy approaches. Noether and Van der Waerden [37], and Friedland [17] gave a proof using classical Bavout's Theorem [43]. As a result of this theorem we prove the following Lemma. Lemma (2.2.2) 1. The system (2.1.3) has n1 solutions. 27 11. If 2 = (21,22,...,zn) IS a solution of P(z) = S. Then for any permutation u of (1,2,...,n), 2 = (ZH(1)'ZH(2)"'°'ZH(n)) is also a solution. arise: i. A straightforward consequence of the Newton's identities implies that 0 is the only solution of P(z) = 0. Therefore by Theorem (2.2.1), P(z) = S has n1 solutions. ii. It is a trivial consequence of symmetry of P(z) = 3. Let us consider the homotopy H ; (0,1) xU an ... Cn defined by H(t,a,z) = (l-t)P(a) -P(z) +ts (2.2.2) where U = ((21,22,...,zn) 6 CD :zi # 2j for i # j] is an open subset of cm. Lemma (2.2.3) For almost every a E U (in the sense of Bair category or Lebesgue measure) 0 is a regular value of Ha(t,z). where Ha : (0,1) an 4 cn is defined by HEJt,z) = H(t,a,z) 28 Proof: Let H(t,a,2) = o for some (t,a,2) e (0.1) xecn. Then DaH(t,a,z) = (1—t)p’(a)- (2.2.3) Therefore det DaH(t,a,2) = (l-t) det p’(a) l 1 l 2a1 2a2 - - - 2an = (l -t) det na?-1 nag"1 na '1 l 1 1 al a2 an n10.-t) det This involves the determinant of the Vandermond matrix, which is well-known, and we obtain det DaH(t,a,z) = (l.-t)n1 H (aj-a.) # 0 1:10ng 1 29 Thus for t # 1, DaH(t,a,z) is nonsingular. Hence for t ¥ 1 Range DH(t,a,2) 3 Range DaH(t,a,z) = C“. (2.2.4) Therefore DH has rank n. This means 0 is a regular value of H(t,a,2). Thus by the generalized Sard's theorem (1.3.2) for almost every a 6 U, 0 is a regular value for Ha(t,2) = H(t,a,2). As we mentioned in section three of chapter one, each component of H;1(O) = ((t,2) :Ha(t,z) = o] is a smooth curve which is diffeomorphic to either a circle or an open interval. Let a be choosen so that 0 is a regular value of Ha(t,z). Let Ta = [(t(>.).z(>.)) :0 g l < 11. (t(O).z(O)) = (0,a)] 1 be a component of H; (O) which contains (0,a) and it is parameterized by a parameter A. Differentiating Ha(t(1),2(1)) = 0 with respect to the parameter 1 implies Dz Ha. z +Dt Ha' t = 0 (2.2.5) 82 62 az ' _ ._l_ __ __r1t ‘ -511; where 2 — ( 51' 51’ ..., 81 and t — d1 . Theorem (2.2.4) Let 0 be a regular value for H then with the a! above notation t(k) is a monotonic function of l. 30 Moreover t # 0 for any (t,2) 6 Pa. Proof: . n n Let us rewrite Ha :C x[0,1] 4 C and = ‘ +' 2 (XI-tiyl, xn 1yn) 111 terms (xf their real and imaginary components obtaining G:rflx[ou]-+mm and W = (x1.y1.....xn.yn) where m.= 2n. By Theorem (1.4.3) we have either (*0 sgn = sgn det G; for all 1 (2.2.6) or sgn t -sgn det Gé for all x (2.2.6.1) where G; is the mjxm Jacobian of G with respect to x and y. It is shown in [21] (also in [5]) that det a; 2,0. Hence by (2.2.6) and (2.2.6.1), t 2_o or t g.0 for all x. This proves the monotonicity of t. If t = 0 for some (t,z) 6 Ta, then (2.2.5) implies that Dz Ha = -P’(z) is singular. So 2 is a multiple root of Ha(t,z) which contradicts the result of Lemma (2.2.2). This completes the proof of the theorem. Without loss of generality we may assume t > 0. If we start at t = 0 and increase the parameter, then t(l) increases. This means t cannot turn around and return back to t = 0. In other words, the monotonicity of t implies that Fa cannot be diffeomorphic to a 31 circle, and each component of Ta with t = 0 at one end will either go to infinity or to a solution of P(z) = 3. At this stage it is important to know whether the component of Ta will diverge to infinity or not. For this we prove the following: Lemma (2 .2.5) 1 The connected component Ta of H; (0) which contains a is a bounded curve. Proof: Suppose this is not the case. Then there is a X'< 1, such that (t(l).z(A)) diverges to infinity as 1 goes to ‘T. Since for any k, Pk(z) is a homogeneous poly— k nomial of degree k (i.e. Pk(12) = A Pk(z)), P (“F—MW) = ——1-——P (2(1)) = -—-1——[ts +(1-t)P (a)1 k (2(A) ')Z(A)Hk k “2(1)”k k k So Pk( :(§)]) goes to zero as 1 goes to I. for all k. Hence, if E'is a cluster point of :(i) we have, P(2)==0 with HE“ = 1. However, this cannot be true because 0 is the only solution of P(z) = 0. §(2.3) Following the curve Integrating the initial value problem P’(2) ~2' 2(0) S«-P(a) a (2.3.1) 32 is what we mean by following the curve. Since P'(2) is nonsingular (2.3.1) can be written as r d2_ I -l ( a- (P (2)) (S-P(a)) (2.3.2) 2(0) = a Nowadays, there are efficient initial value problem solvers which can be used to solve this problem. In general we form the Jacobian matrix p’(z) and solve for 2 which satisfies in (2.3.2). However this procedure uses 0(n3) operations each time. Also 0(n2) function evaluation are needed for computation of p'(2). Reducing the cost of computation is the purpose of this section. We shall utilize the symmetrical structure of p(z) in such a way that (p'(z))-1 may be obtained by 0(n2) operations. To do this we introduce the following: Suppose we have a function g which is known at a set of disjoint points (x1.x2....,xn]. Let us define n (x -xi) q.(x) = II ——-_—, (2.3.3) i7‘j j = 1,2,...,n Clearly qj(xi) = aij (the so called Kronecker delta). The interpolating polynomial defined by 33 n y(x) = Z q.(x)g(x.) (2.3.4) 3:1 3 . J is called the Lagrangian interpolating polynomial. It is well—known that the polynomials ql.q2,...,qn defined by (2.3.3) are linearly independent and form a basis of vector space V of dimension n. Now let us define the following notation which will be needed in the next proposition. For each i = 1,2,... n define oi,...,o;_l as in (2.1.5) with zi deleted. For example 1 _ ol--(22+23+ +2) 01 = (z z -+ -+2 2 ) 2 2 3 n—l n (2.3.4.1) 1 n-l On-l ‘- (‘1) 22 23 Zn Also we define oj = 1 l = 1 2 n . O I I I Proposition (2.3.1) The differential equations(2.3.2) can be written explicitly in the following form: f . dz. n+1 n o1 . ‘57:": ff” - 2: 411(3 -p (a)) .=1 3 II (zk-zl) 3 ( Lzi 0) = a1 i = 1.2, ,n. 34 Proof: We make use of the Lagrangian interpolation formula 11 gj(z) = Z gj(zi)qi(z) (2.3.6) i=1 for functions gj(2) defined by gj(z) = 23 , j = 0,1,...,n -l . (2.3.7) and qi(2) defined by n (2 -2.) qi(z) = n -—j——J—-r i = 1,2,...,n . (2.3.7.1) l . n . 23 = Z) 2; q.(2), j = 0,1,...,n -1 . (2.3.8) :1 1' 1 1 1 --- 1 q1(2) 2 21 22 "° zn q2(2) n- n-1 n-1 ... n- 2 21 22 Zn qn(z) Multiply both sides of this equation by 35 we thus get the following equation, * 2 = W0 (2.3.9) where 1 1 1 221 222 22 W = n-1 n-1 n-1 n21 112.2 nz 0(2) = (q (2) q (2) q (2))t 1 ' 2 n and 2* = (1,22,3Z2,...,n2n-l) Let a11 aln w"1 — anl 0.0. o o ann we can rewrite (2.3.9) as 36 q1(z a11 ' ' ° a1n l \\ q2(2) . 22 ‘ = : . ' - (2.3.10) Also the polynomials (2.3.7.1) can be written as the following +1 . . - (-1)n 1 1 ... i n-1 qi(z) — n [On-1 + on_224- 4—00 2 ], (2.3.11) II ( -z.) i 1 = 1,2, ,n . Comparing the coefficients on both sides of (2.3.10). we get +1 . _ (-l)n ,‘l 1 ij - n j Gn-j . (2.3.12) H( -z.) k=1 2k 1 k¥i The differential.equations(2.3.2) thus become as in (2.3.5). Based on the recursion formula Jj = oj-t-zi Oj-l (2.3.13) for 1 g,i g,n, l g_j g n-l, The computation of Oh-j involves only 0(n2) operations. The values oj's may be obtained by (2.1.6), using forward substition which needs 37 21%ZLL operations, while the total operations needed for the second term in (2.3.13) is (n-1)2. ((2.4) Generalization In this section we will discuss the case that two zeros of f(z) are equal (multiple roots), or nearly equal. Without loss of generality we may assume that 21 and 22 are such zeros. In this case equationsdé2.3.3) are not suitable to us, because evaluation of 757' and d2 d1 is too expensive. Because 1 n ( -2.) viz-k 1 tends to infinity as 1 tends to 1 for i = 1,2. To avoid this difficulty, we write a new differential equations whose solutions are easily led to the solutions of f(z). We define 01 0’2 A. = “'3' + rn‘j (2 4 1) J H (zk-zl) H (ZR-22) ° ° k¥l k¥2 j = 1,2,...,n Proposition (2.4.1) For 1gjgn, Aj(z) is a rational function of z and contains no factor of the form (221-21) Proof: Let us write Aj(2) as follows 38 1 B.(2) A.(2) = x (2.4.2) 3 n ( )( ) (22-21) H 2 -2 2 -2 1511.2 k 1 k 2 where 1 n 2 n B.(2) = o . II ( -2 ) -0’ . IT ( -2 ).(2.4.3) 3 “‘3 k¥l,2 zk 2 “‘3 195,2 21‘ 1 For i # j let us define i i . . _ Sj(ok) 1,3,k - 1,2,...,n to be the same as o; with zj changed to 2i. With these notations we have 1 l _ 2 2 2 _ l 81(on_j) - Un-j . (2.4.4.1) Let us view Bj(2) as a function of 21, and evaluate this function at 22. Then because of equations (2.4.4) and (2.4.4.1) we have Bj(22) = 0. That is Bj(z) = (22-21)qj(2) , (2.4.5) For some polynomial qj(z). Thus B. (22)-Bj(21) . a 11m 3 = —— B. (2.4.6) 22421 22--21 821 j n - a 2 _ ‘ ‘az 0 E (2k 21) 39 A repeat use of recursion formula (2.3.13) combined with equation (2.4.6) imply n-j n qjm = "5?; “30 GMT '22 ei.2(z*<'zl) . (2.4.7) n-j n = - 213—5- (0 II (zk-z1)) Hence n 1 . A.(Z) =q.(Z) n ’ 3 = 1,2,...,n O J 3 k¥1,2 (zk-zl)(Zk-22) Propositiong(2.4.2) * * For lgjgnr let Aj(2) and gj(z) be defined by l 2 2 o . z o . 213(2) = 3231— + n 1 “‘3 (2.4.8) B ( -Z ) H ( -z ) katl zk 1 k¥2 2k 2 *() a 2 F ( ) (2 4 e 1) q. 2 = -—2 0' . 1 -Z . . . . j le l n.) k#1,2 z~k 1 Then * * (qj(z) .Aj(2) = r1 - (2.4.9) II z)( --22) The proof of this proposition is similar to the proof of proposition (2.4.2) and hence is left out. 40 Let us define v = (211(2) ,A2(2),...,An(2))t , (2.4.10) ))t * * * (181(2) .A2(z) .....An(z (2.4.10.1) < ll Then the initial value problem (2.3.1) can be written as d2. 1 _ n+1 t . = d1 - (-1) vi (s-P(a)) 1 1,2 (2.4.11) d2 n+1 n . j ___ (-l) I J _ d1 n L21 2 On-E (Si PL(a)) E (zk-z.) kn 3 j = 3, ,n . . . . _ t w1th the initial value 2(0) - (al-taz, a1a2,a3,.. ,an) . §(2.5) Numerical results In this section a series of numerical examples are given. We have used the ordinary differential equations solver subroutine "DE" (see "Computer solution of ordinary differential equations: The initial value problem” by L.F. Shampine and M.K. Gordon). Subroutine "DE" is probably the most sophisticated and advanced subroutine of this kind. 41 Briefly subroutine "DE" integrates a system of up to 20 first order ordinary differential equations of the form dy. 1 _ . _ dt —-F(t,y1,y2,...,yn), 1 - l,...,n . (2.5.1) The complete description of this subroutine can be found in [40] (the computer program of "DE" is given on pp. 186- 209). The following procedure employing the subroutine "DE" was used. 1) The computer's random function routine generates the starting vectors 21.22,...,2n. We choose two small positive numbers 61. 62 to list the absolute and relative errors respectively. ii) Designate the exist time, called tout, and call the subroutine "DB" in order to follow the path. iii) We write down the output of (ii) and compare them. a) If all solutions are isolated in the sense that they are far from each other by at least 6 > 0, we use equation (2.3.5) to initialize the differential equation and go to (iv). b) Otherwise we use equations (2.4.11) to initialize the differential equations and go to (iv). iv) If tout = 1 we read the solutions and stop. v) Otherwise go to (ii). 42 In the following tables we have shown some numerical results for four systems of equations with known solutions in order to evaluate the accuracy of our algorithm. n Starting Vectors l (.580113649E-POO, .950512735E-t00) 2 (.786371425E-t00, .297620264E-t00) 3 (.453699900E4-00, .626194160E-02) 4 (.275736426E4-00, .305650944E-t00) 5 (.689100711E-t00, .382662239E-t00) 6 (.132902705E-POO, .831857903E-+00) 7 (.582979796E4-00. .98625338BE-01) 8 (.276548455E-t00, .620446028E-t00) 9 (.835029668E-01, .990377121E-t00) 10 (.979346943E-t00, .693884438E-t00) ll (.934477014E-+00, .212092955E-t00) 12 (.130652748E + 00, .862596980E + 00) 13 (.818909294E-+00, .862596980E-t00) l4 (.18799409lE-01, .314116017E-t00) 15 (.765182100E4-00, .941526108E-t00) Table (2.1) The starting vectors for various problems came from Table (2.1). Its first 5 vectors have been used as the starting vector for examples 1 and 4. Also its first 10 vectors are the starting vectors for example 2. Finally all these vectors are used as the starting vectors for example 3. Example 1. We solve system (2.1.3) for n = 5. 51,...,s5 are chosen so that the exact solutions are 5,9,7,4,2. The computed results are as follows: 43 n input vector constant vector 1 .270002+01 (.2421502+02, -.1942712+01) 2 .17500E4—03 (.l47520E-PO3. -.227252E-t01) 3 .1269OE + 04 (.127004E + 04 - .116076E + 01) 4 .98590E«t04 (.985999E-+04, -.411588E4-00) 5 .8003 7E + 05 (.800367E + 05 , .100286E + 01) Table (2.2) Homotopy solutions at time t = 1 2(1) = (.499999999E4—01, -.644220655E-—08) 2(2) = (.9000OOOOOE:P01, -.194643929E‘-O9) z(3) = (.699999999E+01, .119773958E -09) 2(4) = (.4000000002 +01. .627337758E -08) 2(5) = (.ZOOOOOOOOE4-Ol, -.837647611E-—O9) Table (2.2.1) The relative and absolute errors in this problem are respectively 61 = 10.15,e2 = 10-12. Example 2. This problem solves system (2.1.3) for n = 10. sl""'le are chosen so that the exact solutions are l,2,3,...,lO, (in same order). Here 61 = 10-12, 82 = 10-10. The computed results are as follows: 44 n input vector constant vector 1 .5500000E + 02 (.501597E + 02 , - .517790E + 01) 2 .3850000E4-03 (.385669E-t03, -.447630E4-01) 3 .30250002 + 04 (.3027152 + 04 , -.1344092 + 01) 4 .2533300E4-05 (.253342E-t05, -.l30606E-+00) 5 .2208250E + 06 (.220826E + 06 . - .309642E + 00) 6 .1978405E4-07 (.19784lE-+07, .695913E-t00) 7 .l808043E-t08 (.180804E-t08, .242440E4-01) 8 .1677313E-t09 (.l67731E4-09, .283408E-+01) 9 .1574305E4-10 (.157440E-th, .233227E4—01) 10 .1491434E-tll (.149143E-tll, .214617E4-01) Table (2.3) Homotopy solutions at time t = 1 2(1) = (.5000000042+01. -.157844480E -06) 2(2) = (.900000000E-+01, -.115510367E-07) 2(3) = (.599999998E-+01, .146901334E-06) 2(4) = (.399999996E-tOl, .119657234E-—06) 2(5) = (.999999999E-t01, .116296633E -08) 2(6) = (.300000003E-FOl, -.544245438E -O7) 2(7) = (.799999997E-+01, .708308777E-07) 2(8) = (.199999999E4-01, .136081468E-07) 2(9) = (.100000000E4-01, -.l25990729E-08) z(lO) = (.700000003E + 01, -.128213396E -06) Table (2.3.1) Example 3. Here we solve system (2.1.3) for n = 15. 31.32,...,s15 are chosen so that the exact solutions are 1,2,...,15, (in some order). We choose 61 = 62 = 10-8. The computed results are as follows: 45 :5 input vector constant vector \Omflmmbwwr‘ .1200000E + 03 .1240000E‘. + 04 .1440000E + 05 .1783120E + 06 .2299200E + 07 . .3048292E + 08 .4124208E + 09 .6666482E + 10 .7880094E + 11 .1166533E + 13 .1866217E + 14 .8231603E + 15 .3197504E + 16 .4663402E + 17 .6664785E + 18 AAA/‘AAAAAAAAAAA .112498E + 03 : .124059E + 04, .14403 SE + 05 , .178315E + 06. .229920E + 07. .304829E + 08: .412421E + 09. .566648E + 10. . 788009E + 11 , .119653E + 13 . .156622E + 14 , .2231603 + 15 o .319750E + 16, .46034OE + 17, .665479E + 18: .80489534-01) .74363834-01) .301077E + 01) .259791E + 00) .923576E4-00) .251047E-t01) .283309E + 01) .132393E-t00) .326317E4—01) .162573E4-01) .262221E-t01) .365233E4-01) .3907OSE + 01) .134212E-F02) .172283E4—02) Table (2. 4) Homotopy solutions at time t 1 2(1) 2(2) 2(3) 2(4) 2(5) 2(6) 2(7) 2(8) 2(9) 2(10) 2(11) 2(12) 2(13) 2(14) 2(15) AAAAAAAAAAAAAAA . 90003 903 9E + 01 , .119999114E + 02 , . 799988123E + 01 . .60001 748 7E + 01 , . 700006205E + 01 , .999997814E + 00 . .999971337E + 01 , . 200001713E + 01 . .499975294E + 01 , .130000293E + 02 . .139999773E + 02 , .299992826E + 01 , .150000010E + 02 , .400017191E + 01 , .110001766E + 00. .177185166E - 02) .182960120E -03) .218749510E -02) .144389314E -02) . 205918686E - 02) .1656887023 -O6) .110689513E ~02) .417415352E - 05) . 735860613E - 03) .47507804'7E - 04) .391695467E - 05) .531830915E - 04) ' .363097092E - 06) .256633723E - 03) . 522092661E - 03) Table (2.4. 1) 46 Remark. The above examples along with a lot of other numerical experiments suggest that as long as n, the number of the equations, increases, we would have less accurate results. This is due to the subroutine "DE" and presumably because this subroutine is coded for at most 21 equations. It should be mentioned that in each problem we have converted n complex variables into 2n real variables and with the parameter time t totally we worked with 2n4-l real variables. Other subroutines such as "DeVerk", "D Gear" from IMSL are also used to solve these problems. But overall our numerical experiments show that among all available differential equations solvers "DE" is the best one, even though it is not very accurate for problems with relatively large dimension. Example 4. This problem solves system (2.1.3) for n = 5, and its exact solutions are 2, 5, 5, 7, 8. Here we follow the curves until two of them get close to each other, that is,the distance between them be less than 10-5. Then we follow the curves by integrating of € differential equations (2.4.11). The relative and absolute errors in this problem are respectively 6 = 10-15, 1 -12 62 = 10 47 n input vector constant vector 1 .27000E-t02 (.242150E4-02, -.l9427lE4—01) 2 .167OOE-+03 (.166520E-+03, -.227252E4-01) 3 .11130E-+04 (.111404E-t04, -.ll6076E+—01) 4 .7763OE-t04 (.776399E-FO4, .411588E-+00) 5 .55857E-+05 (.558567E-t05, .100286E-t01) Table (2.5) Homotopy solution at the time n t = .899 where the distance between two curves is less than c. 1 (.500000358E4-01, .226529962E-04) 2 (.799999999E4-01, .459513577E -lO) 3 (.700000000E4-01, -.l34110233E -09) 4 (.499999642E-t01, -.226529034E -04) 5 (.l99999999E-t01, -.834203779E -11) Table (2.5.1) Homotopy solutions at t = 1 2(1) = (.999999999E-+01, .220664352E-09) 2(2) = (.799999999E-+01, .728877445E1-10) 2(3) = (.7000000002+01. -.198322803E-09) 2(4) = (.ZSOOOOOOOE-tOZ, .101079963E-—08) 2(5) = (.199999999E4-01, -.96786l485E -lO) Table (2.5.2) CHAPTER THREE REAL ZEROS OF A FUNCTION §(3.l) Introduction In this chapter, we shall be finding the real roots of f(z) = 0 (3.1.1) in the interval [-R,R], (R > 0), where f(z) is an analytic function on an open set U containing the interval [-R,R]. We shall use Cauchy's formula _ 1 f'(z) sO — 27d (as f(z) dz (3.1.2) Where S is a rectangle with vertices tR:tie and entirely contained in U and BS is the boundary of S, (e is a small positive number). S0 is the number of zeros of f(z) inside the contour as, and hence is an integer number. Thus we need only a crude estimate of S When we find an approximation for SO, say so, 0' then the dimension of our problem will be the integer nearest to so. Having the dimension (number of zeros) we will employ the generalized Cauchy's formula 48 49 _ 1 i f'(z) . _ Si - ———2m._ (as 2 f(z) dz 1 - 1,2,...,n (3.1.3) to approximate the numbers 51' $2,...,Sn, where n 1 Si = Z) 2. , 1 = 1,2, ,n i=1 3 These numbers along with the homotopy continuation method presented in chapter two will determine an approximation for the real zeros of f(z). In order to obtain a more accurate result, the approximated zeros 21,22,...,zn could be a good starting point for using the Newton's iteration method. That is to refine the approximations obtained by the homotopy continuation method, one or two- step method can be used as the final stage in calculation of the real zeros of f(z). §(3.2) Integration around the boundary of a strip Suppose g 6 C[0,l] and define m le'“(g) = fi . ’g(j/m) J=0 (3.2.1) = {13% 9(0) 4.9%) + ~-+g(1m"—1-) +% 9(1)). The prime on the summation indicates that the first and the last terms are assigned a weighting factors -%. Using Romberg's method of integration, RIm'1](g) is evaluated for mesh m = 1,2.4.8,... 50 Let us define k T3 = R[2 '1] k = 0.1.2.3.... . (3.2.2) 0 l 2 . . Then we evaluate TO, To. TO' ..., and With the aid of the following recursion formula Tk = —JL— (4 m Tk+l - Tk ) (m > 0) (3.2.3) “‘ 4m—1 m-l m-l we form the following T-table: 0 To 1 0 T0 T1 T3 Ti TS (3.2.4) 3 2 1 0 T0 T1 T2 T3 4 3 2 1 0 T0 T1 T2 T3 T4 Equation (3.2.3) implies that Tk is a linear combination m k k+l k+2 k+m . of T0, T0 . TO . .... TO . That 18 k m m-i k+i Tm = 12% cm TO . (3.2.5) Since 9 is Riemann integrable on [0,1], as k goes to infinity the first column of the T-table, namely the 1 sequence (TE? converges to ( g(t)dt. O 51 The leading term in the error of TE} is of the order 1 2k+2 OTEFI , (we will see this in the next section) while 2 k 1 2 k that of T0 is (ERII) Thus Tm converges much more rapidly than Tk After all the convergence of the 0' first column of the T-table implies the convergence of all further columns and all diagonal sequences (T:. with k constant) to the same limit. We observe that dividing h successively by 2 is a sound one, since each application of the composite rule (namely, equation (3.2.1)) utilizes all of the function values previously used for smaller values of k and m.= 0. A division of h by some other number would require introducing a complete set of new function values. Hence the above form of the Romberg's method of integration based on trapezoidal rule requires fewer function values than might at first appear. Moreover the recursion formula (3.2.4) is a simple manipulation and requires minimal computer time. The most time consuming part of the algorithm is the calculation of TE. k = 1,2,... , the original composite trapezoidal rule. §(3.3) Error analysis Here we consider the error estimate of numerical approximation to i f’(2) . 1555— d2, 1 1 S. =-—. 1,2,...,n (3.3.1) 1 2W1 yes 52 where S denotes a narrow strip with vertices :tR:tic. (e is a small positive number). a b = -R+—ie a = Ri-ic 7| 1. c = -R-ie d = R‘-ie Figure (3.1) f’ 2) Let us assume M.= sup f(z) , and choose a small 2688 enough (with S < R). Then -1— ~ 2'1 _’(_ _ zi f’ _(_z) 26M |21ri ($2 dz z+27ri i—z f—_(2) dz ) 1“”. mi (3.3.2) 1 = 0,1,2,... Having a small enough, the amount of integration on the two vertical sides of the strip S is contributing a small amount in the whole process of integration. Let us denote this amount by the initial error EI‘ Next we consider the integration along the other two sides. By using a simple change of variable we may consider the 53 integration on the interval [0,1]. In order to analyze the corresponding error we define the following. Definition (3.3.1) A function o is of class E (even) if i) m(X) is of definite sign for 0 < x < 1, ii) o'(x) is of definite sign for O < x < 1/2, iii) cp(0) = ep(l). cp(x) = p(l-x) and cp(x) = cp(1+x). Definition (3.3.2) A function u is of class 0 (odd) if i) ((x) is of definite sign for 0 < x < 1, ii) t(O) ((1/2) = ((1). iii) $(x) -w(1-x), and ((x) = u(l-tx). For example p(x) = [sin wx| and t(x) = sin 2vx are respectively of class E and 0. The following lemma is a straightforward result of the definitions. Lemma (3.3.3) a) If p(x) is of class B, then x N1(X) = ( (0(2X)-o(x))dx , O and x l 42(x) = ( m(x)dx-x ( o(x)dx O O are of class 0. 54 b) If ((x) is of class 0, then x p(x) = I0 t(x)dx is of class E. Let us define the following functions inductively: -% x(l-—x) if 0 g.x g,1 K2(x) = (3.3.3) K2(x-n) if ngxgn+lo Suppose K2m is defined, then we define x K2m+1(x) = (O K2m(x)dx (3.3.3.1) where ‘— ._ 1 K2m(x) — 4m-1 (K2m(2x)-K2m(x)) .(3.3.3.2) Also we define x K2m+2(x) = (O K2m+1(x)dx . (3.3.3.3) It is not hard to show that the functions K2m(x) are all positive in the interval [0,1]: more exactly monotically increasing from x = 0 to x = 1/2 and from there again decreasing to zero [33]. Theorem (3 .3 .4) Let g eczm+2[o,1]. Then 55 1 ‘(11H'1)k .i‘ K (Zm-I'Z) Tfi-Ig = 4 (2kx)g (x)dx . (3.3.4) 2m+2 1 where Ig = ( g(x)dx. 0 Proof: The remainder term for trapezoidal rule is given by ([3]) TO I - [15(1 t) ”(t)dt 0" g I a 2 ‘ g ' 0 and likewise k -k 1 k ,, TO-Ig= 4 F K2(2 t)g (t)dt “0 We prove the theorem by induction with respect to m. Assume 1 k _ -km k (2m) Tm_l -Ig - 4 (‘ K2m(2 t)g (t)dt , (3.3.5) where sz is a function of class B. This assumption is certainly true for m = l. with K2(x). Substitution of' (3.2.3) into (3.3.5) gives Tk m = ———(l m-t)-4'”‘"+1)‘“1<2m(2k"14‘2kmxmk(2 t) “L1 0 .19 g(zm)(t)dt . Hence 56 4km 1 k+1 k ( 4m_1 (O [K2m(2 t) -x2m(2 t)]g T:-Ig = 4 2m)(x)dx. (3.3.6) By the way that kernel functions are defined and lemma (3.3.3) Rém(x) is of class B. Also 1 __ 1 l x2m+1(1) = j K2m(x)dx = Ty (K2m(2x) —K2m(x))dx 0 4 -1 0 1 l (2 (1 = [—, K (x)dx- K (x)dx] = 0 .. (3.3.7) 4m_1 2 "0 2m 0 2m Similarly (1 1/2 K2m+2(1) = JO K2m+1(x)dx = (11/2 K2m+l(x)dx = o . (3.3.8) Clearly Rém(x) is also even and has period one. According to this and equations (3.3.7) and (3.3.8) a twofold integration by parts yields 1 1 - k (2m) _ _k k (2m-l-2) (O_K2m(2 t)g (t)dt — 4 (0 K2 I2(2 t)g (t)dt Therefore 1 k _ -(m+1)k k (2m)-2) Tm-Ig—4 (‘0 K2 I2(2 t)g (t)dt This completes the proof of the theorem. Let us define l __ k Em;k — Tm-(O g(t)dt Then equation (3.3.4) implies 57 1 - m+l k k 2m+2 {IEm'kl g4 ( ) [(0 K2m+2(2 t)g( )(t)dt| l = 4‘(m+1)k|g(2m+2)(§)| If K2m+2(2kt)dt| (3.3.9) 0 for some 4 6 [0,1], (the last equation is true because the kernel function K2 has a fixed sign in [0,1]). m+2 In order to arrive at an exact error estimate, we need to derive a compact form for 1 K = (O K2m+2(t)dt . x2m+2 = (2m+2)'. ' the“ Let us apply equation (3.3.4) for g(x) 0 1 1 Tm-(O g(t)dt = (O K2m+2(t)dt = K Because of the recursion formula (3.2.3), every entry of the T—table is a linear combination of values of the first column,that is m-i Tk-I-l Cm 0 (3 .3 .10) *3 B:# II II 0 MB 1 where C3 are independent of k. In fact with the assumption Cm-1 = C-1 = 0 we have m m 1 l. m i i-l c = (4 c -C ) m 4m_1 m-l m-l X2m+2 If the Euler-MacLaurin sum for g(x) = 755:577 is written in terms of Bernolli coefficients 58 (see [3 J. p. 207), we can see that -m(m&l) B 4 _ m 2m+2 K" (‘1) (2m+2)T tm (4 i ' where ( 4 1)( 4 l) (l 1) 1 -— 1 ___ o g g -— ufil)__ 4 42 4m 1 l 1 (1 --) (1 -—) ... (l -—) 4 42 4m and B2m+2 is the 2m4-2pg_Bernoulli number. Therefore B2m+2 Thus because of (3.3.9) we have 4-(nH-l)k .32 (2 2) m+2 m+ -(m+l)k _ 4 (2m+2) ' "('2m+'2): 5 “3(5) l ' BZm-I-Z where ‘5==;me:I) and goes superlinearly to zero [33]. Let p be the shortest distance of any singularity of g(z) from the interval [0.1]. Then for any 50 > 0 there exists a constant such that (2m+2) (2m+2) '.K(6O) [9 (5) i (p-é )2m+2 O Let us consider the special case where,_g(2) has the form 59 _ i f’ z) 9(21-2 T‘s—((2‘)— That is g(z) may have only simple singularities. Thus a constant K exists such that (9)352) I g -(———)-—2'§+mi2"K (3.3.13) P In order to consider the error of integration on 33 we use transformation 2 = a-+(b-a)t. Thus d = 2pR .is the shortest distance of any singularity of g from line segment dbl 4'41““)k 6 , (2m+2) '.K 2R lEch —<- (2m+2): p2m+2 _ 2m+2 _ (14k)(m+1) . K 2m+4 — 62 2m+2 R d l 2 . Let Em,k and Em,k be the errors corresponding to the horizontal sides 53' and 33 respectively. Then (14k)(m+1) l3; kl g 59.21“” fem”) K -2 = 1.2 (3.3.14) Let = E19 -+E2 then EF m,k m,k ' IEFI g 25R2m+4 d-(2m+2) K .2(1-k)(m+1) g 25R2m+4 ,K 2(1-k) (m+1) ,e-(2m+2) (3.3.15) 60 Now let us consider the corresponding error for SO and use the fact that S is an integer number. In this case 0 equation (3 .3 .2) is of the form 2 (Ell $7.3M Note that the integration along the vertical sides of the strip 8 are ignored. In order to determine SO it is enough to assume the following .1. l lEIl g7, “31?“ g 3 (3.3.16) 2 2 Hence 26R2m-I-4K 2(1-k)(m+l) €-(2m+2) £13 (3.3.17) 2 and 2_€ 1 .T Mia-3;. (3.3.17.1) This leads to the following inequalities 14k (_2_' n 02, gag—4:, (3.3.18) 2 M __l_ 2 where C = R(166R2K)2 Based on the above argument we present the following flow chart. This flow chart gives a general framework for the algorithm. Input 61 F,t Compute Ti Check inequalities (3.3.18) NO Compute T1 i Check Inequalities (3.3.18) 4!” ‘j - j + 1 YES Write the diff. eqs. (2.4.11) YES Is there any two solutions close to each other? N0 Determine dimension n k I k + l Y Write the diff. eqs. (2.6.11) Compute and rec rd S k L Y Subroutine "DE" Flow chart (3.1) Increase time and continue 444-8 integration until t - l YES 62 §(3.4) Numerical results A series of computer experiments were conducted to show the accuracy of the algorithm. As the Flow chart shows, we limit our discussion here to the homotopy type algorithm of chapter two of this thesis. The computed results for the well known Bessel function are shown in Table (3.1). Example: Here we consider the real roots of Jl(z) in the interval [-4w,4n], where Jp(z) is the Bessel function of the first kind of order P defined by to (_1)m(z/2)P+2m me = “EC miI‘(m+P+l) and is the solution of the following differentiable equation w”+ z‘lw’ + (1 -P2z'2)w = 0 First of all by a crude estimate of S we observe 0 that J1(z) has seven zeros in [-4w,4v]. In fact 2 = 0 is a zero of J1(z) (trivial solution). Then we compute S ,S .,S Having these we employ the homotopy 1 2).. 6o continuation method to determine nontrivial zeros of J1(z). First we choose an arbitrary vector A = (A1,A2,...,A6), then we follow the solution curve of the homotopy equation (2.2.2) until t = 1. 63 n Input Vector l (.100893900E + 00, -.158603000E -01) 2 (.334024966E + 03 . -.2304932313 + 00) 3 (.447584441E + 01 . .458962664E + 01) 4 (.266768698E + 05 , —.138646669E + 03) 5 (.113324558E + 03 . .900923 546E + 03) 6 (.245160906E + 07 . -.232377074E + 05) Table (3.1) This table shows the input Cauchy vector se= ($1,...,s6), where si is the sum of the ipthower 6 of zeros of J1(2). Table (3.2) is an arbitrary vector in C . Finally, the computed zeros of Jl(2) is given in Table (3.3). Arbitrary Chosen Vector A l (.580113649E +00 , .950512735E + 00) 2 (.786371143E +00 . .297620264E + 00) 3 (.453366999E +00 , .626194161E - 02) 4 (.275736426E +00 . .305650943E + 00) 5 (.689100711E +00 , .382662239E + 00) 6 (.132902705E +00 , .83185‘7903E + 00) Table (3 .2) 64 Computed Zeros of J1(z) 2(0) 2(1) 2(2) 2(3) 2(4) 2(5) 2(6) .367073133E + 01 . .144951670E - 01) .101561514E + 02 , —.1o43299oos -o1) (-.707456470E + 01 . .536099550E - 02) (-.101568583E + 02 . .257022866E - 01) ( .707868536E + 01 , .450648139E - 02) (-.366324131E + 01. -.5549224o42 —01) 0. (trivial solution) ( ( Table (3.3) CHAPTER FOUR EIGENVALUE PROBLEM §(4.l) Introduction In this chapter we shall describe an algorithm for finding the eigenvalues of a given n.xn matrix. The homotopy continuation method of chapter two is implemented to approximate the eigenvalues of a matrix. We shall present an efficient recursive formula to find trace (Ck), (k = 1,2,...,n), where C is an n‘xn companion matrix. It has long been known that the eigenvalue of a matrix can be found by solving its characteristic equation. Namely the eigenvalue problem reduces to solving a polynomial, P(A) of degree n with 1 — n n- 0.0 P(X) — ). +P 7. + +Pn_11+Pn (4.1.1) 1 A large number of methods are available for the determination of the zeros of P(X). The most well known and simple classes of methods are iterative methods. For example the Newton's iteration method P (in) A -r§TT:;7 (4.1.2) xn+1 = n 65 66 is one of the most popular methods. Such methods have a number of difficulties, among them i) the starting point must be fairly close to desired solution, ii) even if we choose a "good" starting point, the rate of convergence may not be as fast as we wish, iii) some polynomials of high degree are highly sensitive to their coefficients. §(4.2) Method and Theorems Let us start this section with the following theorem which will be needed later. Theorem (4.2.1) Let A = (aij) be an n xn. matrix With akk-l # 0, for some 1 < k gin. Then A is similar to a matrix B = (bij) With hkk-l = 1 and bkj = O for j #'k-l. Proof: Let us define a matrix U = (Uij) as follows: r -1-— if i = k-l. j = k-l. akj U.. = (ii if i = k-l. jello-1. (4.2.1) 13 a'kk-l L aij 1f 1 # k-l, j = 1,2,...,n Then if j¥k-lo u 8 Define C = AU, with C 67 n ij = 1E: akiUij = akk-lUk-lji'aijjj ..ak' - akk-l (__lakk-l) + akj — 0 , (4.2.2) and if j = k-l n ij = Ckk-l = iii akiUij = akk-lUkk-l = 1° (4°2'3) -The construction of U implies U is invertible, and -1 _ . U - (Vij) With akj if i = k‘J-p j = 102: Inl V. . = 13 éij if 1 ¥ k-l, j = 1,2, .,n It is clear that the kph_row of AU is the same as the kph_row of U-lAU. Hence if we define B = U‘lAU, then A and B are similar and B satisfies the stated condition. This completes the proof. Let A = (aij) and ann-l # 0, then A 13 Similar to a matrix B1 with the npp row of the form (0,0,...,l,0), that is, there exists an invertible matrix U1 such that 1 l l b11 b12 °°' b1n -1 B1 — U1 AU1 — t 1 l l bn-11 bn-12 °°' bn-1n 0 0 ... 1 0 (4.2.4) 68 Lets assume bh-ln-Z # 0, then there exists an invertible matrix U such that 2 _ -1 _ - -1 _ -1 132 — U2 131112.. 1121111 AUle — (Ule) A(U1U2) 2 2 2 b11 b12 bln = 2 2 2 bn-21 bn—22 bn-2n o 0 1 o o O O C) 1 0 Suppose the above procedure can be done successfully n-l times (the case which this procedure fails will be discussed later). Then there exists a sequence of matrices Ul'U2’°"'Un-l such that B = (U U U )‘1A(U U U ) = U'lAU n—l l 2 n-l l 2 "' n n-1 n-1 n-1 b11 b12 ' ° ' bln l 0 0 = O 1 0 (4.2.5) 0 C). l 0 with U = U U ...U 69 Remark 1. Suppose after n-k steps of the above procedure we end up with a matrix B = (bij) where bkk-l = 0, but still there exists some j(1 g,j < k-—l) with bkj ¥ 0. In this case we postmultiply and premultiply this matrix by the matrix E (elementary row operation), where .. 3m row E = (4.2.6) ' 4 k-lpp row The effect of this multiplication will not change the similarity, but it will interchange the jpp_and k-lph columns, and simultaneously changes the rows with the same numbers. So we can proceed with the above algorithm. Remark 2. If after n-k steps of the above procedure we end up with a matrix of the form B with 70 11 12 ' ' ' 1n B = , 4 kph row Then B can be viewed as the following A1 A2 B = 0 A3 with “11 ' ' ' “1k-1 akk ° ' ' C1kn A = ' A = I 1 3 “k-11 ' ' ' ak-lk-l 0 Clearly l is an eigenvalue of A if and only if A is an eigenvalue of A1 or A3. A3 is already in companion form, and we continue to find the corresponding companion matrix for A1. Overall a certain matrix A can be transformed to a companion matrix C with the same set of eigen- values. In order to find the eigenvalues of C we first state the following lemma [38]. 71 Lemma (4 .2 .2) Let A be an n.xn matrix with eigenvalues 11.12.....1n. then R n trace(A ) = 23 i. k= 1,2.... =1 In the following we will present a recurrence formula to compute the entries of matrix Ck from matrix Ck.l for a companion matrix C (here C1 is the ipp power of C). Proposition (4.2.3) Let c1 C2 Cn-l Cn 0 C = I O . . k k be a companion matrix. Suppose C = (Cij) (the kph_ power of C) has been calculated. then the entries of Ck+1 = (CEEI) are given by k+l _ k k - _ k+1 _ k Cln — Cncll , (4.2.7.1) k+1 _ k - = Cij - ci-lj , 1 2,3,...,n (4.2.7.2) 72 Proof: Since Ck+l = C ~Ck = C .(Cij)' it is clear that Ckfl = C3 1 = 2,3, ,n 1] 1-1] 3 = 1,2, ,n Also since Ck+l = Ck ~C, for j = 1,2,...,n—l we get k+l ’1 k k 1c Cij g LE1 ciLCLj — C3 ll'Fclj+l ’ and similarly k+l _ k C1n ' CnC11 Remark 1. This proposition enables us to find the entries of Ck+l from Ck. Namely we only need to compute the first row of Ck+1. Because the second up to npp rows of Ck+1 are respectively the first up to n-lpp,rows of Ck by (4.2.11.2). Hence Ck+1 is obtained from Ck at the cost of n multiplication. Therefore n(n-l) multiplication is necessary to find all entries of C2,C3,...,Cn. Having found the entries of C1‘ (k = 1,2,...,n), we compute Sk = trace(Ck) : k = 1,2,...,n . In order to determine the eigenvalues 11,12....,xn of matrix C (and hence the eigenvalues of A) we may 73 solve the following system of equations: r A1 + X2 + "' + An ’ S1 2 2 2__ "1+’\2+ "‘+"n‘52 ) (4.2.8) n n n _ (*1 + 12 + ... + xn - Sn From here on we employ the homotopy continuation method which is described in chapter two of this dissertation. §(4.3) Numerical results In this section we shall consider two examples. Example 1. Let 1 2 3 A = -l -2 l 1 0 —1 Since 332 = 0, we make use of the elementary matrix B where and transfer A to a similar matrix B, with B = EAE Let T L1 then -1 U1 Hence -1 U1 BUl Define U2 (3 O fiflH 74 I hflw H Then the corresponding matrices as follows. I l )4 O hflm V C, C and C are 75 _ -1 c .. (EUIUZ) A(EU1U2) — 1 o 6 4 ~16 02 = —2 2 8 . and c3 Hence we have the following system 11 + 12 + 13 = 2 2 2 A1 + 12 + 13 1.3.43: 3 11 At this stage we use the algorithm 2 8 0 O , 1 O of chapter two to find 11, 12, and 13. The computed results are shown in table (4.1). 76 . t 5(1) (-O.ZE+01, 0.02+o.o) mp: 5(2) (+0.8E+Ol, o.or.+o.0) ““3 °r3 5(3) (+0.4E+Ol 0.01~:+o.0) arbitrary 6(1) (.580113649E-+00. .950512735E-t00) choosen 6(2) (.786371425E4-00, .297620264E-t00) eigenvalues 6(3) (.453699900E-t00, .62619416lE-t00) constant Pa (1) (- .382018497E + 01 . -.125439494E + 01) vector Pa (2) (0 . 783133603E + 01 - .157657305E + 01) Pa Pa (3) (0.500647786E + 01 . -.630500849E + 00) eigenvalues x(l) (-.l87555035E4-01, .102511753E-t01) of the given 1(2) (0 .17511007OE + 01 . .238749298E - 11) matrix A 1(3) (-.l87555035E+01, -.102511753E+ 01) output 5 ’ (1) (-.2000000002 + 01 . -.2295053042 - 11) vector 5 ’ (2) (o . 7999999992 + 01 . o .873 967565E - 11) s ' s ’ (3) (o .399999999E + 01 . 0 .840714165E - 10) Table (4.1) Here 81’ 82. and S are respectively, trace(A), 3 trace(Az), and trace(AB). Likewise Pa(1), Pa(2), Pa(3) are trace(M), trace(Mz), and trace (M3) for an arbitrary choosen matrix M with known eigenvalues 6(1), 6(2), and 6(3). Finally 1(1), 1(2) and 1(3) are the computed eigenvalues of the given matrix A. Example 2. In order toillustrate the accuracy of our algorithm. we will find the eigenvalues of a diagonal matrix D, 77 1 d2 0 D = . 0 . d10 where di (for i = 1,2,...,10) are given in Table (4.2). The corresponding companion matrix C is computed as c1 c2 ... C10 0 C= I 0 where ci (for i = 1,2, .,10) are also given in Table(4.2). 1 l 2 3 4 5 6 7 8 9 10 d1104O-10-3O 5.00020 —40 20 3.0 C 5.0 30.0 -150. -273. 1365. 820. -4100. -576. 2880. 0.0 Table (4.2) 78 Since C10 = 0, the last column of C is identically zero. This implies that zero is an eigenvalue of C. Hence we may delete the last column and last row of C and denote the remaining part by C1. To compute trace(Ci) (for i = 1,2,...,9) we observe the recursion formulas (4.2.7), (4.2.7.1), and (4.2.7.2) and compute all elements below or on the main diagonal of C3. These are shown in Table (4.3). In fact the first 9 rows from the top of this table are the elements of C3, and from the second row up to the tenth rows are the elements of CE, and so on. Overall we form a system of nine equations with unknowns 11,12, . . . ,kg . Then we use the homotopy continuation method of chapter two to solve this system. The computed results are shown in Tables (4.3) to (4.5). 79 1m.co manna OOOOOOOO Ommm OO¢¢H OOOOOOOH mhml O OmNhHI OOOOOOHO OOH¢I @hOHNI OOmmmml Omhvwaal OOOr-(OOO O moma memo ¢NO¢N thm5 0 mammmm Ommwmw Omnmmbm OHmOOHVH O N OCDOOOOOHOO O C) 0 r4 0 C>..t) = ([D+t(A-D) -11]... é—(xtxq'n. (5.2.3) where D is an arbitrary diagonal matrix with distinct elements. Applying this homotopy has the following disadvantages: a) if we follow the n distinct curves suggested by M. Chu we may not get all n eigenpairs, since two of these curves may link into a pair of ~eigenpairs of the form (x,l) and (ex,k). (so they actually represent one eigenpair). b) to get all eigenpairs we actually must follow 2n distinct curves rather than n curves. In order to remedy this problem and solve (5.2.2) at a reasonable cost,a special homotopy is constructed as follows: Let D be an arbitrary diagonal matrix with distinct elements on its diagonal. Construct the homotopy equation‘ Hun“ x112 XJR 41R“ xm defined by H(X.X.t) 1 n t 2' = ([(1-t)D+tA-xI]X,e Z xi-t(x x) -1). (5.2.4) i=1 where e is a small positive number. It is clear that vectors x(0) ei/e = 1 = 1.2. ,n 1(0) di are the eigenpairs of H(x,x,0), where ei is the standard ipp unit vector and di is the ipp element of the diagonal matrix D. We should mention here that the crucial step in applying the homotopy continuation method is the construction of an appropriate homotopy, such as (5.2.4) so that i) the existence of a curve connecting the trivial solution and desired solution is assured and ii) the numerical work in following this curve has a reasonable cost. In the next section we shall show that the homotopy equation (5.2.3) guarantees the existance of n distinct smooth curves. Each of them leads from an obvious starting point to a desired eigenpair. Furthermore if a certain curve links to an eigenpair (:91), then there is no other curve that may link to (-:nk). These curves are characterized by an explicit ordinary differential equation with distinct initial values, and hence they can be easily followed by any ordinary differential equations solver Coupled with the large scale matrix techniques, this method can be used to solve eigenvalue problems for sparse matrices [22]. 86 §(5.3) Theorems In this section we present some theorems which serve as a theoretical basis for our algorithm. Theorem (5.3.1) 0 E IRn le is a regular value for H. In other words, for each (E. I, f) 6 IRn le le with H(;.I,t-) = 0. H The Jacobian D(- .1. t) has rank n4—l. Proof: Let ($55.35) 6 En xIR le and H(x,1-,t_:') = 0. Observe that ”(iraH = r ._ ._ ._ ._ ;_‘ (5.3.1) (l-t)p+tA-).I -x (A-D)x -1 _1 1 _ 2 _ _ 2 _ _ 2 L(e-t(§tx) x1.....e-t(§tx) xn) o-(F‘x) Since H(;hiuf) = 0, we have n __ .1. e: Z". x. -t(x x)21 (5.3.2) .=l i and ((1-E)D+'EA-II) ~35: 0 (5.3.3) We claim that the (nntl) x(ni—l) matrix 87 (1 -‘E)D+'€A-§I ..x (32.33) - 1 (e-thtE) x1....,e-t( 2 NH4 _-t— — x x) xn) O is of full rank. Since otherwise there exists a vector (37,“)t with y em“, and tem such that y (0) ’ H' = o O . (5.3.5) Thus ((1-t)D+tA-3:I) °y-ux This implies Qt . ((1-E)D+EA-3.'I) 'y = Qt}. Since A is symmetric, and E’ is orthogonal to the rowspace of the (1 -t)D+tA -i-I, we have (THEM; = (SIZE: (it - (1 -E)D+EA -fI\ oyxt t ... ... — .— = y '((l-—t)D-+tA-AI) -x = O This implies u = 0. Therefore ((1-E)D+tA-II) -y = o . (5.3.6) Since the matrix (l-—t)D-+tA-—ii has a simple spectrum ([10] , Lemma 6 .l) , we conclude that the matrix B = (l-t)D-+th«-Xl has a set of orthogonal eigenvectors. say, zl,z2,...,zn, with corresponding eigenvalues 88 61.62,...,6 . Let a -+c z -+----+a z 121 22 nn' xl n Y = 3121+3222+ "'+ann Then by (5.3.3). (5.3.6) and orthogonality we get a. = 8.5 i = 1,2,...,n 1 1 Therefore y = 5;, for some 5. Substituting y = 5x and u = O in (5.3.4) we get (1-E)D+tA-KI -x 5x 0 NH4 1 2; .e-t(§t§) E o o o ._ -t— E. t(XX) 1:... n Thus by (5.3.2) we have, 1 n _— o= we 2: xi-t(§tx) 2) = a. (5.3.7) i=1 Hence D(;'X)H is of rank n-tl. This completes the proof. Remark. We have restricted our discussion with the Jacobi structure of the matrix A. This is needed only as a sufficient condition for Theorem (5.3.1). This condition may be rephrased as "choosing D so that the matrix (1 -t)D-ttA has a simple spectrum.for any t E [O.lft Overall, it is only needed that the matrix D( H be x.k.t) 89 of full rank for any (x.)..t) e mnxm xm with H(x,k.t) = O. Apparently the sparse matrix techniques can be incorporated in any of these cases. As is the usual procedure of the homotopy continuation method we start from a trivial solution of H(-,-,O) at t = O, and follow the generated path as t increases from zero to one. We hope the trivial eigenpairs deform into the eigenpairs of the original matrix A. Hence we would be able to follow the n distinct connected paths from the trivial system to the original problem. In order to assure that this process works. we prove the following; Theorem (5.3.2) Let us define .T‘= ((x,>.,t) e Zan xIR le :H(x,>..t) = o) a) F is a one dimensional smooth manifold, b) as t increases, the curve T will never turn back. 3232:: Part (a) is in fact a standard result from the differential topology [35]. That is, a repeat use of the Implicit Function Theorem implies that P consists of one dimensional manifold. In order to prove (b), let T be parameterized with a parameter, 9. Along each component, we may take the 9O derivative with respect to the parameter a. The set P is then characterized by 2122 dt de - —-+ H = . . D(t)H d6 Dow) O (5 3 8’ :11. d9 . dt . . We claim ‘d§'# 0, Since otherW1se d6 D(x.>.)H = 0 El Lde Hence D(x l)H would be singular. This completes the proof. Now consider starting a path at t = 0. Since t' # 0. without loss of generality, we may assume t’ > 0. So as we increase the parameter The following lemma shows that as F remains bounded. the set F O ((x.x) e 1Rn le: H(x.1.t) is bounded. Lemma (5.3.3) The set Tb 9. t(e) cannot reverse. t goes to one the curve In other words for any 0 < tO < 1, O for some t E [to.l]] (5.3.9) is a bounded set. 91 Proof: By equation (5.3.3) lMt) IHXH = HUI -t)D+tAIXH g H(1 -t>D+tAEHIXH 3 (ND!) +))A‘.))HXH Hence l).(t)l g. HDH + HAM . (5.3.10) Also for any (x.l) 6 TO, equation (5.3.2) implies n Il-te .2) xi) touxn s. (x) . i=1 n g.l-+e Z) Ixil == l4—eHle i=1 S.l'*€JSMXH Thus for 9 small I ) 1 (X) g -—7=_ . (5.3.11) t0 6 n Therefore To is boundedo (Here by H H we mean ||||2.) Part (b) of Theorem (5.3.2) implies that P will never turn‘badk.Thus T can be parametrized by the variable t. Then (5.3.8) becomes r- "' -1 r- T 3ng (1 -t)D + tA -).I v] (D-A)X = -l _l .1. L%% Ls--t(XtX) 2 x1.....e-t(xtx) 2 Xn O -(xtx)2 A ‘ (5.3.12) X(O) ei/e = l = 1,2, ,n t l(0) d 92 These differential equations may be solved by any ordinary differential equations solver. Over all for each 1 g_i g_n the numerical solution of (5.3.12) at t = 1 is an approximation for an eigenpair of the matrix A. The following lemma assures that among the computed eigenpairs we will not have a pair of eigenpairs of the form [X.).]t and [-x.l]t. Lemma (5.3.4) There is no pair of eigenpairs of the form [x.l]t. [-x.l]t. Proof: Since otherwise both vectors satisfy equation (5.3.2) with t l. Subtracting the two resulting equations n gives 2) xi = 0. Substituting this in (5.3.2) yields i=1 HXHZ = -1. Which is not true. §(5.4) The Algorithm We have seen that the set I‘ = {(x,).,t) :H(x,).,t) = O} , (5.4.1) consists of exactly n curves, and each of them is the solution of the differential equations (5.3.12). We proceed as follows: I) Choose any diagonal matrix D with distinct elements on its diagonaL II) Set i = O. BI BLIOGRAPHY BIBLIOGRAPHY Abraham, R. and Robbin, J. (1967). Transversal mapping and flows, Benjamin, New York. Allgower, E. and Georg, K. (1980). Simplicial and continuation methods for approximating fixed- points and solutions to systems of equations, SIAM Review, vol. 22, No. 1. Bauer, F.L., Rutishauser, H. and Stiefel, E. (1963). New aspects in numerical quadrature, Proc. Sympos. Appl. Math. Vol. 15. pp. 199-218. Amer. Math. Soc. Providence, R.I., MR3O #4384. Chow, S.N., Mallet-paret, J. and Yorke, J.A. (1978). Finding zeros of maps: Homotopy methods that are constructive with probability one, Math. Comput., 32, pp. 887-899. . . (1979). A homotopy method for locating all zeros of a system polynomials, Functional differential equations and approximation of fixed points, Springer lecture notes in mathematics, No. 730, pp. 77-88. Chu, M.T. (1983). A simple application of the homotopy method to symmetric eigenvalue problems, "to appear". Collatz, L. (1966). Functional analysis and numerical mathematics, Acd. Press. New York. Davidenko, D. (1953a). On a new method of numerically integrating a system of nonlinear equations, Dokl. Akad. Nauk SSSR, 88, pp. 601-604. (in Russian). , (1953b). On the approximation solution of a system of nonlinear equations, Ukrain, Math. 2. 15, pp. 190-206 (in Russian). 94 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 95 Deift, P., Lund, F. and Trubowitz, E. (1980). Nonlinear wave equation and constrained Harmonic motion, Comm. Math. Phys. 74, pp. 141-188. Delves, L.M. and Lyness, T.N. (1967). A numerical method for locating the zeros of an analytic function, Math. Comp., vol. 21, pp. 543-560. Drexler, F.J. (1977). Eine method zur Berechnung II . . Samtlicher l'sunger von polynomgleichunges systemen, Numerische Mathematic 29, pp. 45-58. , (1978). A homotopy method for calculating of zeros of zero-dimensional polynomial ideals, in continuation methods, ed. H. Wacker, New York: Academic Press. Eaves, B.C. (1971a). On the basic theory of complementary, Math. Programming, 1.1, pp. 68-75. , (1972). Homotopy for computation of fixed point, Math Programming, Ibid., 3.1, pp. 1-22. , (1976). A short course in solving equations with PL homotopies, SIAMqAMS Proc., 9, pp. 73-143. Friedland, S. (1977). Inverse eigenvalue problems, Linear algebra and Appl., 17, pp. 15-51. Garcia, C.B. and Gould, F.J. (1978). A theorem on homotopy paths, Mathematics of operational research, vol. 3, No. 4. , and Li, T.Y. (1980). On the number of solutions of polynomial system of equation, SIAM Journal of numerical analysis, No. 17, pp. 540-546. , and Zangwill, W.I. (1978). Global continuation methods for finding all solutions to polynomial systems of equations in n variables, in: Symposium on external methods and systems analysis, Springer Verlag, Berlin. , (1979). Determining all solutions to certain systems of nonlinear equations, Mathematics of operations Research, vol. 4, No. 1. , (1979). An approach to homotopy and degree theory, Mathematics of operations Research, vol. 4, No. 4, pp. 390-405. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 96 Hirsch, M. (1963). A proof of the nonretractability of a cell onto its boundary, Proc. Amer. Math. Soc. 14. pp. 364-365. Karamardian ed. (1977). Fixed points algorithms and applications, Academic Press, New York. Kellogg, R.B., Li, T.Y. and Yorke, J.A. (1976). A constructive proof of the Brouwer fixed- point Theorem and computational results, SIAM J. Numer. Anal., vol. 13, No. 4. , , . (1977). A method of continuation for calculating a Brouwer fixed point, computing fixed points with applications, Kuramardian, ed., Academic Press, New York. Kuhn, H. and Mackinnon, J.G. (1975). Sandewich method for finding fixed points, J. optimization Theory Appl.. 17, pp. 189-204. LahayeviE. (1934). Une methode de resolution d'une categorie d'équations 'trancendantes, C.R. Acad. Sci., Paris, 198, pp. 1840-1842. . (1935). Sur la representation. des racines systemes d' equations transcendantes, Deuxieme congres national de sciences, 1, pp. 141-146. Lemke, C.B. and Howson, J.T. (1964). Equilibrium points of bimatrix games, Siam J. Applied Math., 12, 2. pp. 413-423. Li, T.Y. (1976). Computing the Brouwer fixed point by following the continuation curve, Proc. Sem. Dalhousia Univ., Halifax, Academic Press, New York, pp. 131-135. , (1982). On locating all the zeros of an analytic function within a bounded domain,"to appear". Lyness, J.N. and Delves, L.M. (1967). On numerical contour integration round a close contour, Math. Comp. vol. 21, pp. 561—577. Merrill, O.H. (1972). Applications and extensions of an algorithm that computes fixed points of a certain upper semi-continuous point to set mapping, Ph.D. Thesis, Dept. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 97 Milnor, J. (1965). Topology from the differentiable viewpoint, Univ. Virginia Press, VA. Morris, J. Li. (1983). Computational methods in elementary numerical analysis, John Wiley and Sons. Noether, E. , Van der waerden, B.L. (1928). Die alternative bei nichlinearen Gleichvgen, Nachrichten der Gesselschaft der Wissenschaften zu G6ttingen, Math. Phys. Klass, pp. 77-87. Ortega, J.M. and Rheinboldt, W.C. (1970). Iterative solution of nonlinear equations in several variables, Academic Press, New York - London. , Raleston, A. and Rabinowitz, P. (1978). A first course in numerical analysis, McGraw-Hill, Inc. Saigal, R. (1977b). On the convergence rate of algorithms for solving equations that are based on complementarity pivoting, Math. Oper. Res., 2. pp. 108-124. Scarf, H. (1967). The approximation of fixed points of a continuous mapping, Siam J. Appl. Math., 15,5, pp. 1328-1343. Shampine, L.F. and Gordon, M.K. (1975). Computer solution of ordinary differential equations: the initial value problem,. Freeman and company, San Francisco. Van der Waerden, B.L. (1931). Modern algebra, I and II. New York: F. Ungar Publishing. Walker, H., ed. (1978). Continuation methods, Academic Press, New York. Wilkinson, J.H. (1965). The algebraic eigenvalue problem, Clarendon Press, Oxford.