MSU LIBRARIES .—_—. RETURNING MATERIALS: P1ace in book drop to remove this checkout from your record. FINES w111 be charged if book is returned after the date stamped be10w. NUMERICAL METHODS OF BIFURCATION PROBLEMS VIA SINGULAR VALUE DECOMPOSITIONS AND HOMOTOPY METHODS By Yun-qiu Shen A DISSERTATION Submitted to Michigan State University In partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY ' Department of Mathematics 1988 ABSTRACT NUMERICAL METHODS OF BIFURCATION PROBLEMS VIA SINGULAR VALUE DECOMPOSITIONS AND HOMOTOPY METHODS By Yun-qiu Shen A relation between bifurcation theory and the singular value decomposition, homotopy methods in numerical analysis is studied. Given a nonlinear equation, we give a local analysis in a neighborhood of a solution via the Liapunov-Schmidt method and the singular value decomposition. This analysis is applicable to regular, turning or bifurcation points. In the case of a bifurcation point, homotopy methods are used for solving the bifurcation equation. A numerical method for global bifurcation prob- lems based on the above analysis is presented. ACKNOWLEDGMENTS It is a great pleasure to express my sincere appreciation to Professor Shui-Nee Chow, my thesis advisor, for all his invaluable guidance, expert advice, stimulat- ing discussions and encouragement during the course of my research. I would also like to thank Professors B. Drachman, J. C. Kurtz and E. Palmer of Mathematics Department and Professor L. Ni of Computer Science Department for their patient reading my thesis and attending my defense. A special thank goes to my wife, Guobei, for her support and patience. Finally I wish to thank all those people who have taught and helped me during my stay at MSU. ii TABLE OF CONTENTS Chapter 1. Introduction ...................................................................................... 1 Chapter 2. The Singular Value Decomposition ................................................ 4 2.1. Definitions and Theorems ........................................................................ 4 2.2. Computational Methods ........................................................................... 6 Chapter 3. The Singular Value Decomposition and the Liapunov-Schmidt Method .......................................................... 11 3.1. The Liapunov-Schmidt Method ............................................................... 11 3.2. A Connection of SVD and LSM .............................................................. 12 Chapter 4. A Numerical Bifurcation Equation ................................................ 16 4.1. Theorems .................................................................................................. 16 4.2. Proofs ....................................................................................................... 20 Chapter 5. A Probability One Homotopy Method ........................................... 24 5.1. Introduction .............................................................................................. 24 5.2. A Probability One Homotopy Method for the Bifurcation Equation ....... 25 Chapter 6. The Singular Value Decomposition and Regular Points .............. 30 6.1. A Matrix Result ........................................................................................ 30 6.2. A Newton’s Method via SVD .................................................................. 33 Chapter 7. A Numerical Method of Global Bifurcations ................................ 35 7.1. Rank Detecting, Curve Following and Switching via SVD ..................... 35 7.2. A Numerical Method of Global Bifurcations via SVD ............................ 37 7.3. An Example .............................................................................................. 40 List of References .................................................................................................. 42 iii LIST OF FIGURES Figure 1. An Illustration of the Algorithm ........................................................ 39 Figure 2. Five Computer Graphs for the Example .......................................... 41 iv CHAPTER 1 INTRODUCTION Let F: X xR‘ -> Z be a Fredholm operator, where X, Z are Banach spaces and R‘ is the usual 5 dimensional space in which the parameter vector It sets. The Liapunov - Schmidt method can be used to reduce F(x,h) = O to a finite dimensional system f(y,7L) = O with f: Rme‘ —9 R". The bifurcation behavior of the solution set of f determines the bifurcation behavior of the solution set of F. In the boundary value problem of ODE or PDE, a finite element method or a finite difference method also can be used to approximate the problem by a finite dimensional system. If the original differential equation contains parameters, then the finite dimensional system which is of the same form as f obtained from the Liapunov- Schmidt method also contains parameters. Furthermore, many problems themselves are finite dimensional prob- lems. Therefore studying the behavior of the finite dimensional system with parame- ters is meaningful. In this paper we present a new numerical methods for the bifurcation problem f(x,?») = O with f : Rme —->R”‘. The methods include using reliable ways to distin~ guish the bifurcation and non-bifurcation, to factor out the bifurcation equation, and to solve the problem in either cases. Bifurcation occurs at the points of the solution set when the J acobian matrix of f at these points is rank deficient and a point of the solution set is near a bifurcation point when the Jacobian matrix at that point is nearly rank deficient. We relate this J acobian 1 matrix with the singular value decomposition(SVD). Because an effective way to detect and neat rank deficient problems in numerical analysis is to compute the singu- lar value decomposition. see [6]. If the singular value decomposition of the Jacobian matrix at some point is performed, we certainly want to use the informations from it as much as we can. In this thesis, we connect it with the Liapunov- Schmidt method in the bifurcation theory. The major results which relate this are: Theorem 3.2.2 The Liapunov-Schmidt method via S VD ; Theorem 4.1.1 A numerical bifurcation equation via SVD; Theorem 6.1.1 A matrix result via 5 VD for the Newton's iterates; Algorithm 7.2.1 A numerical method of global bifurcations via S VD. After a numerical bifurcation equation, which is a system of special polynomi- als, is obtained via SVD, we need a reliable method to solve it. The development in 1976 by Chow, Mallet-Paret and Yorke [2] is an advance in the homotopy methods. The methods are called the probability one homotopy methods which provide practical ways for solving nonlinear equations [10]. For the general case of a syStem of polyno- mials, the problem of finding all zeros has already been solved in [3]. The numerical bifurcation equation which appears in this paper is a system of special polynomials. Only some of the solutions are to be solved. The others can be obtained by symmetry. We develop a special homotopy equation to find these required solutions by using some techniques in [2] [3]. We obtain a method which only requires half of the usual number of computations. The result is obtained in: Theorem 5.2.6 A probability one homotopy method with symmetry for solving the numerical bifurcation equation. In most chapters, we assume f satisfies the following two conditions: (a) fis C", 1:22; (b) there are only finite smooth curves passing each bifurcation point of the solu- tion set of f(x,A) = 0. CHAPTER 2 THE SINGULAR VALUE DECOMPOSITION § 2.1. Definitions and Theorem Singular value decomposition (SVD) is one of the most important tools in matrix computations. It is a reliable method for detecting and treating rank deficient prob- lems. In this section we briefly describe the results of SVD which are needed in this thesis. We will present them by restricting to real matrices. For proofs or details, see, for example, Golub and Van Loan [6], Dongan‘a et a1 [4]. Similar results for complex matrices also can be found in Stoer and Burlisch [9]. But for the purpose of our research, we only need the case for real mauices which is presented in the follow- ing definitions and theorems. Let A be a real m by n matrix. It is known that there exist an m by m orthogonal mauix U and an n by n orthogonal mauix V such that A=U£VL (LU where D 2, = 0’ g (1.2) is an m by n matrix, and D, is a r by r diagonal mauix. Furthermore D, = diag(0'1,0'2, ° - - ,o,), where 612022 - - - 20', > 0. VT is the transpose of V. 4 Definition 1.1. (1.1) is called the singular value decomposition of the matrix A. Let p = min(m,n), and define o,+1=o,+2=...=cp = 0. Then 01,62, - - ' ,op are called the singular values of A. Theorem 1.2. (a) The number of nonzero singular values is equal to the rank of A ; (b) oioé, . . . , a? are the positive eigenvalues of ATA and AAT ; (C) the columns of V are corresponding eigenvectors of A TA ; (d) the columns of U are corresponding eigenvectors of AA T . Definition 1.3. The columns of V are called the right singular vectors of A, while the columns of U are called the left singular vectors of A. Definition 1.4. The pseudo-inverse A1” is defined tobe the mauix A+=v2+ UT, (1.3) where D"1 0 2*: 0' 0 (1.4) . 1 1 m D 1=c1 —,—-, , — w1 rag( 1 02 0r ) It’s easy to verify the properties of the pseudo-inverse: Theorem 1.5. (a) AA+A =A; (b) A+AA+ =A+; "I '1 (c) AA+=U (33 UT; I l (d) A+A = V 11. 0 L00 (c) (1111+)2 =AA+; (r) (A+A)2 =A+A. Regarding A as a linear transformation from R" to R’" under certain bases, denoting U = [ u1,u2, ° - - ,um] and V = [ v1,v2, - ' - ,vn ], the following theorem is obtained: Theorem 1.6. (a) NUll(A) = Span{ vr+lavr+2’ . ° . 9vn}; and (b) Range (A) =Span { u1,u2, -°-,u,}. § 2.2. Computational Methods Theoretically Theorem 1.2 already gives a way to compute SVD of A. i.e., ATA and AAT can be used to obtain SVD of A. However it is not a satisfactory way for the computational purpose due to that the round off errors often destroy pertinent infor- 1 O mation. For example, letA = [a 0] ,where (1 satisfies 80 < a < V80 and 80 is the O a l-l-a2 0 machine precision. Then theoretically ATA = 0 a2] gives the singular values of A with 0'1 = Vl+a2 and 02 = a since obviously the eigenvalues of ATA are 1+ a2 and a2. But the computational results due to the round off errors yield ATA = [(1) 8] gives the singular values of A withol = 1 and 02 = 0 . The second singular value is qualitatively incorrect. The basic computational method is the Golub-Reisch SVD algorithm [6] in 1970 which contains Householder bidiagonalization, the decoupling calculation and the Golub-Kahan SVD step of a bidiagonal square matrix having no zeros on its diagonal and superdiagonal. We will briefly illustrate these methods. For more details, see [4] [5] [6] [7] . Let A be an m by n matrix. In this section, we assume rn 2 n, otherwise con- sider AT. Two orthogonal mauices are involved in this algorithm, they are: Definition 2.1. A Householder matrix U is a mauix of the form U = I - 2uu TluTu, where u is a column vector. Definition 2.2. A Givens rotation matrix U is a mauix of the form .1 ] cosO sine -sin6 c056 Then the matrix A can be transformed into a bidiagonal mauix by: Theorem 2.3 ( Householder Bidiagonalization ). There exist products of House- holder mauices U3 = Ule - - - U" and VB = V1V2 - - ~ Vn_2 such that B UBTAVB = ....... (2.1) 0 , where difZ dzfs B = d3 .. . . f" L dud cat11- attest 11:00 #1100 #11100 $11100 tar-00 taste it! not one 111-1- n-c #:1- tottulxgcctxvlgcctuzxgottxvz80¢9U3xgotgu‘x80c9 -> —) —> -> -) —) not: 0:111:11: 0:11:11: 00:12-13 0011-11- 000: 000* Ito-11:1: 0111*: 0:11:31: 00:11:11 0011:11- 000: 0000 U 2C OC ‘1’} .6. (VBVC)T. If B has a zero on its superdiagonal, B can be immediately decoupled into two smaller upper bidiagonal square matrices. If B has a zero on its diagonal, B also can be decoupled into two smaller upper bidiagonal square matrices by multiplying by a Givens mauix. Therefore using the decoupling calculations, upper bidiagonal square submatrices with no zero on their diagonal and superbidiagonal always can be obtained except for these submatrices on the diagonal of B with size 1. An example of a 5 x 5 upper bidiagonal mauix with (3,3)-th entry zero is illus- trated as following( +,0 denote the enuies which are changed in that step): sacoo Now we discuss a method to diagonalize the above mentioned upperbidiagonal submatrices . Without losting generality we just assume B with no zero on its diagonal or superdiagonal since the decoupling calculations always can be used to change to smaller mauices. Choosing Givens rotations 55.3.1.1, Tim. i=1,2,...,n-1, leaving T13 open, B is transformed as following: at 1 ’*#+ 1 1:0 1 + III III 0 at a: a1 :11 XTIJ .1. :1: be a1: :1: xTu + :1- :1: B _.’ * ... _) * _) * . a1: :1: :1- 1: :1: a1: '1: :1: 1 :1- 1: 1 * * + T * * sgsx 0 a1: :11 San—1.11x 111 —> * -> -> :1: ... :1- 1: 0 :1 ~ _ T _ T B -(Sr.252.3 ‘ ' '5n—1.n) 3(T1,2T2,3 "'Tn-1.n) -5 BT (2.2) which is again an upper bidiagonal mauix. Now T13 is chosen as following: 10 ' 1 costtz simttz ‘Slmbiz COS¢12 T12 ___ 1 (2.3) 1 such that cos¢1,2 -Sin¢1,2 di 'IJ- __ [1‘] (2 4) Sin¢1.2 COS¢1.2 d if 2 0 . . where u is the eigenvalue of 2 2 dn-1+fn-l daz-lf; (2.5) dn-lfn dn+fn which is closer to d?, + f3. The algorithm for obtaining E from B by (2.2) - (2.5) is just the Golub—Kahan SVD step. Let B“) = B. Recursively 8““) can be obtained by using B“) instead of B in the above step. If for some B“), there are zeros on its diagonal or superdiagonal, smaller matrices will be considered. It is known that fan—>0 and dam->0" at least qua- dratically when i —9 co. Disregardin g exception, the convergent rate is even cubic. Hence the singular value decomposition of B always can be obtained by itera- tions. The subroutine SSVD in Linpack [4], SVD and MINIF in Eispack [5] are all based on the Golub-Reisch SVD algorithm. A different way of computing the bidiagonalization in the Lawson-Harson algo- rithm in 1974 is to upper triangularize the mauix A first, which is faster when m > n. The subroutine LSVDF in IMSL Library [7] is based on the Lawson-Harson SVD algorithm. CHAPTER 3 THE SINGULAR VALUE DECOMPOSITION AND THE LlAPUNOV-SCHMIDT METHOD § 3.1. The Liapunov-Schmidt Method Many problems in analysis and applied mathematics can be reduced to the deter- mination of the zeros of a function in a Banach space. A bifurcation occurs when a multiple zero exists. A technique which is called the Liapunov-Schmidt Method (LSM) can be used to simplify bifurcation problems. In this section we want to estab- lish the connection between the singular value decomposition and the Liapunov- Schmidt method. Therefore a reliable numerical method can be used to solve the bifurcation problems. We first state the following theorem which gives the Liapunov-Schmidt method. Theorem 1.1 ( Chow-Hale [1]). Suppose X, Z are Banach spaces, A: X —> Z is a continuous linear operator, N: X —9 Z is a continuous nonlinear operator and I is the identity operator in X. Let W and E be continuous projections in X and 2 respectively. Suppose Null(A) = Xw , Range(A) = 25 , where XW is the range of W in X and 25 is the range of E in Z. Then there is a 11 12 bounded linear operator K: Z5 —9 X1_w. called the right inverse of A , such that AK = I on Z3 , KA = I - W on X. Moreover, the equation Ax-Nx=0 (1-1) is equivalent to the following equations z-mm@n)=a um demon)=a am where x = y+z, yeXw, zeX,_w. For the proof of this theorem, see [1]. § 3.2. AConnection of SVD and LSM Consider a C 1 function f: Rme‘ —) R’". Such a function can be viewed as a function from R’" to itself with 5 parameters. In our discussion, we treat Rme‘ as Rm+3 The first derivative of f at some point x0 6 Rm” can be represented as a real m by m+s Jocabian matrix, which is denoted by Df(xo). Assume the rank of this mauix is r, 0 S r S m. The singular value decomposition of Df(x0) is denoted by Df(xo) = UoZng . The following Lemma is obvious. 13 Lemma 2.1. DO 0 U5 is the projection from R" to {Range[ Df(xo)] }; (a) U0 IVrO 0V ois the projection from Rm” to { Null[Df(xo)] }—; (13) V0 (c) [Df(xo) 1+ = V023 US is the right inverse of Drag). Proof: Identify Df(xo) as A in Theorem 2.1.5, then Theorem 2.1.5 /(c),(e) give (a); Theorem 2.1.5/(d),(f) give (b). (Df(xo) )+ is just the pseudo-inverse of Df(xo), therefore Theorem 2.1.5/(c),(d) give I, 0 Df 0:0)th (xoir = U0 [0 OJUE I, 0 lDf (xo)]+Df (X0) = Vo [ OJVT’ then (a),(b) and the definition implies (c). C] Theorem 2.2. Let f: R’"+’—>R"‘ be C1, f (x0)=0. Let Df (x0)=U0>:0V[; has rank r with OSrSm. Then there is a neighborhood 1] of x0 such that f(x)=0 , x e n if and only if z—zo=Vo£3Ug [Df (x0)(x—x0)-f(x)] and (2.1.a) [m ' §u0.i,r+1fi(x) =0 , (2.1.b) §u0.t.M(x) 14 where z - 20 is the projector of x - x0 in Rm” to { Null [ Df( x0) ] }l which is the orthogonal complementary subspace of Null[Df( xo)] in R’"” and f1(x) f(x) = . . . . firm In the case r=m, the second equation disappears. Proof: Consider f(x)=0, which is same as Df (xo)(x -xo) - [Df(onx-xo) -f(x)]=0 - (22) Regarding Df(xo) as the operator A, Df(x0)(x- x0) -f(x) as the continuous operator N(x), x- x0 as x in Theorem 1.1, the right hand side of (1.1) in the theorem gives that (2.2) is equivalent to: (Z ‘Zo)-KE [Df (onx *xo)-f (101:0 (I-E)[Df (xo)(x-Xo)-f (x)]=0 , (2.3) where (z- 20) is the projector of x - x0 in R’"” to { Null[Df(x0)] jl . From Lemma 2.1/(a), we have: (l-E)[Df (xo)(x-xo)-f (10] ’ 1 =Uo g I": UgonZoV€(x-Xo)-f(X)l q , =U0 8 ,0 Hit-foo] . “"1 15 -(uo,,1 1f1(x)+uo.,21f2(x)+-,+u0m.1fm(x)) (“our mfl(x)+u02 mf2(x)+ +110»th (16)) 01m- 0 0 '- 2 uO,i,r+1fi(x) i=1 --Zuo.i,ntfi(x) ~ i=1 From Lemma 2.1/(a),(c), we also have: I [(5 = (vozanTXUo 0' 0 011107) = virivi Combining (2.3) and the above two equalities, the conclusion of the above theorem connecting SVD and LSM is followed. III In Chapter 4 and Chapter 6, we will restrict to the case s = 1. We are going to change (2.1) further in different situations according the rank r. In the deficient rank case, we obtain a bifurcation equation via the singular value decomposition in Chapter 4, and in the full rank case, we obtain a mauix result via the singular value decomposition in Chapter 6. CHAPTER 4 A NUMERICAL BIFURCATION EQUATION § 4.1. Theorems SVD can be used for detecting bifurcations because that it is a way for detecting rank deficiency, and a bifurcation occurs at the point where the rank of the J acobian matrix of the derivative of the map is deficient. Now we are going to derive a bifurca- tion equation by using the Liapunov-Schmidt method via SVD which is a system of polynomials whose coefficients are expressed in terms of the enuies of the orthonormal matrices in SVD. In Chapter 4, 6, 7, we use s=l, i.e. f: R’"+1—>R’" is a function containing one parameter. In Other words, we deal with one-parameter problems or equivalent one-parameter problems ( if several parameters are involved, we give more condi- tions to change them to one-parameter problems). Furthermore we assume that there are only finite smooth curves passing each bifurcation point of the solution set of f(x)=0, and f e C " , k 2 2 for some integer k which is discussed later in this chapter such that (1.1) holds. In this chapter, we assume f(xo)=0 for some x0=(x0,1,x0.2, - ~,xo.,,,.,1)T in Rm“. Df(x0)=Uozovg withrankr, OSrSn—l, where V0.1.1. ' ° ° V0.1.m+1 V0 = O O O O O I V0,m+1,l ° ° ° v0,m+l,m+1 . l6 17 Ziumnifior) Consider the Taylor series of about x=xo. Obviously the constant £u0.i,n1fi(x) i-l L term is zero, Also its linear term is zero, since 1 §u0,i,r+lfi(x) D( .. (xo))=D(Uo g lilugfxxo) .Zluosnfdx) =Uo 3 ,3, U3(Df(xo))=(Uo 3 ,3, levozovg) =0. Denote x=(x1,x2, - - - ,x,,,.,1)T . We assume that f has enough smoothness such that Taylor’s Theorem holds: "" ' ' a ,. "' i§u0,i,r+lfi(x) dr+l! [(xi-xo,1)3;+° - -+(xm+1-xo,m+1) <9me Id l(gr‘iuon',+1fi(xo))l = + m I a d,” m guahmfdx) HUM-10,1 )Xl—i' ’ ' ' +(xm+1-Xo.m+1) axm+1 I (iéuo.i,mfl(xo)) L J . ___.__..l _ __8__ __ ‘9 dr+1+1m . . _. (dr+1+1)! [(Xi JCo,1)axl + +(xm+1 xo'm+l)8x,,.+1 I (Euo,.,,+1fi(xo+e(x xo)) 3 axm+l WK? 1-xo,1 )fii' ' ' ' +(xm+1-xo,m+1) 1d"I+l (gammaxowu-xo» .J (1.1) 18 where O < 9 < 1, d122, j=r+1, r+2, , m are positive integers and the second term of the right hand is continuous. Actually if there exists a positive integer k for the smoothness such that: k> max (dp +1), (1.2) r+l 0, hence the second term in each component goes to zero. For the first terms in all 22 components, the orders of the derivatives are just the same as the orders of the powers of II x - x0 ll , therefore we distribute through to each of xl—on, xz—xo; ,..., a denominator II x - x0 ll . By all the hypotheses of f assumed, the xm+l—x0,m+l limit can be performed in all above mentioned quotients, and the limits of these quo- tients are just the components of the unit tangent vector, by (a), they are m+1 m+1 m+1 z YjVOI. j, Z iijOHZj’ 2 yJ-vanHJ. The limit of the second term of each j=r+l j=r+l j=r+l component is zero. Therefore the limit of (2.3) gives (b), 1 e . m+l "1+1 d In - [( Z ij0.,1,')% +( Z ijOHm+lj)ax I '+1(2uo.i,r+ifi)(xo) j=r+l j=r+l ”1+1 i=1 =0 (2.4) m+1 m+1 [( z ijOHl '93—:- +( 2 Yij, 111+] j)ax Idm(zu0,i,mfi)(x0) j=r+l j=r+l i-I (c) of Theorem 1.1 is obvious. Hence the theorem is proved. Cl Proof of Corollary 1.2: Corollary 1.2 can be immediately followed since the second derivative can be written in the matrix form: m+1 m+1 [( Z ijOHl j)—:—a +( z ij0,m+l ,j)ax—'—T1I2(g(x0» j=r+l j=r+l "1+1 01+] 01+] 32800) =2 ( 2 ijOHij 2 lijOqu)——’ ax Max p,q=l jar-+1 j=r+l "1+1 n+1 = [( Z 1ijO,,l j) ( Z lijOm-i-l j)I j=r+l j=r+l _ 823010) 32 g(xo) q r mil .v . I axnlaxl a151339;.“ j=r+lyl 0'” 328 628 m+l . ___.(_XL)- . . . _—(x(.))_ Z ijO,m+l,j dxml 8x1 8xm+1 met ] j=r+l 23 V0.1,r+l v0.m+l,r+l = (YH-lv ' - - 9ym+l) V0.1.m+1 v0.m+l,m+l P 2 1 a gtxo) 328(10) 8x 131 1 axlaxnwl Vo,i.r+1 V0.1.m+1 Yr+1 3280(0) 328(10) v0.m+l,r+l Vo,m+1.»t+1 Ynt+1 axm+l 8x1 axm-I-laxm-I-l J The bifurcation equation in Theorem 1.1 and Corollary 1.2 are systems of polyno- mials. They can be solved by the probability one homotopy methods which can guarantee finding all the roots. We will describe them in the next chapter. CHAPTER 5 A PROBABILITY ONE HOMOTOPY METHOD § 5.1. Introduction The bifurcation equation reduced from bifurcation problems via SVD in Chapter 4 is a system of special polynomial equations. There are m+1~r equations and vari- ables. For sake of convenience, in this chapter, we use n instead of m+1-r. Denmc ( n+1 .yr+2. "'1)’m+1 ) by Z = ( 21.22. “.21. ) and n polynomials by P1(Z),P2(Z), - ~ ~ ,Pn(Z) which can be regarded as components of a polynomial vec- tor P(Z). Homotopy methods can be used to solve nonlinear equations, i.e., if we want to solve P(Z) = 0 , we first solve a simple equation Q(Z) = O, and then set a homotopy function H(Z,t) = (l-t) Q(Z) + t P(Z). Solve H(Z,t) = 0 by following the homotopy curves ( solution set of H(Z,t) = O ) from t=0 to t=1, hence the zeros of Q(Z) lead to the zeros of P(Z). The nonsingular Jacobian matrices of H(Z,t) are important for tracing the homotopy curves. The development in 1976 by Chow et al can finally avoid singular Jacobian matrices by constructions called the probability one homotopy methods . In the methods, for almost all the choices of the homotopy parameters, the methods are globally convergent. This is an advance over earlier homotopies, since the philosophy and the resulting software are fundamentally new [ 10]. The numerical bifurcation equation is a system of special polynomial equations. 24 25 Observe that if Z = ( 21,22, . . . ,zn ) is a solution of this system, then -Z = ( -zl,-zz, . . . ,—z,, ) is also a solution of it. Therefore we wish to construct a sym- metric homotopy function such that only half of the solutions need to be computed. In the next section, we will use the fundamental idea of "probability one" and the tech- niques in [2][3] to construct such a function. § 5.2. A Probability One Homotopy Method for the Bifurcation Equation The system of polynomials appear in the bifurcation equation above is a system of homogeneous polynomials with degree d;22 for each polynomial P;(Z), i = 1, 2, ..., n-l, except the last one Pn(Z) = 2% +z§+...+z?,-1. Without losing generality, we can assume that the first 5 polynomials are of odd degree, i.e., d,-23, fori = 1, 2, ..., s, and the rest of the polynomials are of even degree, OSsSn —1. Note first that there is at least one even degree polynomial ( the last one ), and secondly that Z = 0 is not a solution of the above system, although it satisfies all the polynomial equations except the last one. Thirdly DP(O) =0, since P(Z) has no linear term. We constrict the following symmetric homotopy function: 1‘11(Z.t) ”(2,1) = ° Hn(Z,t) . d - 211-17121 ,, a, .d. 0 2‘11,ij Z,‘-b,2_, [31(2) 131. =(1-t) 2311‘ -b.+t + t - ~- + NH) . d (2.1) . . . Pn(Z) an-lszn-l 2:11:11 Tbn-l 1:1 2 ’ .2121-bn j zlanrjzj j I1: 26 where Z = ( zl,z2, ° ° ' ,zn)T e C" and H(Z,t) = (H 1(Z,t),I-12(Z,t), - - - ,Hn(Z,t))T e C". The parameters are chosen in random by 2 a=(a11,... ,a1ma21, . . . ,a2,., . . . ,am.)Te C" andb =(b1, . . . ,bn)Te C". The difference of this homotopy function and the known homotopy function is that b 121, . . . , bsz, are used instead of b1, . . . ,bs in the first 3 components. There- fore the transversality condition should be checked. We have: Lemma 2.1. Let W= { (a,b)e anxC" |b1,-~-,b,, #0 }, and A={ (2.0 e n . . 8(H11'HaHn) C xR }. Then the submatrix of the Jacobi of (2.1) 8( Z a b) has full rank on A x W ( i.e., rank is n if one regards the mauix as a n x [n + n2 +n ] complex mauix, or the rank is 2n if one regards the matrix as a 2n x [ 2n + (2n)2 + 2n] real matrix ). Proof: Case 1: t = 0. Consider a submatrix : l' l -1 dIZl' —b! ..zI 1,-1 8” 61,2, "b‘ —— a In -1 312.1») d ' .MIzloI _1 d...‘ “l dl-Izl-l (2.2) . . d--l . . which has full rank Since one of djzj’ — bj, -zj lS nonzero, for 1 s 1 SS. There- fore the assertion is true in this case. Case 2: 16 (0,1) and 21, °--,z,, are not all zero. Consider a submatrix: 27 1(1-1 "11-1 1" a” )2" '( )2 10-021" ..... ((14):? an ...... 1" *4 ..... l :4 1( [)11 K .4): ((14):? ..... “1'02: (2.3) which has rank n as soon as one of 21, - - - ,z,, is nonzero. Case3: 16 (0,1) and zl=--- =z,,=0. In this case 3—2 = 0 and also the partial derivatives to Z of the terms in (2.1), which have degrees of Z greater than one, are zero. iii-obi 0 . aH _ —(1-t)b, ' ' 0 . 3(Zb) - 0 —(1_:) 0 -(1-1) (2.4) which gives rank n. [I] Now define the homogeneous part of (2.1) by: - Putz.» H(Z.t)= ,_ H,(Z,t) ct n d1 .. 2a 1' '2 ' 21" 191(2) ,.t ’ ’ =<1-t) a... H pn'jjz) +t(1-t) . a.-. (2.5) :21 2%4- : ° ° +23 [Earn—1,11} . Xaw'zi i=1 : 28 Since 6_H is just as case 2 of (2.3) , we have: 6a Lemma 2.2. The Jacobi of (2.5) has fullrank w.r. t. (z,t) x (a,b) on { (C” - {0}) x (0,1)] x W. Now we need a transversality theorem (see [ 3 D. Definition 2.3. Let F be a smooth map : open set A c; Rd-eR” , then a point y e R” is called a regular value of F on S g A provided that Range { DF(x) } = RP for all x e S n F ‘1 (y). Those x’s are called regular points. Theorem 2.4 ( Transversality Theorem ). Let A 1; Rd and W ; R4 be open sets, and F: A x W —> R" be C' smooth with r > max {0, d-p}. Suppose for some set S g: A that y e R” is a regular value of F on S x W. Then for almost every w e W ( in the sense of either Baire category or Lebesgue measure ), y is a reg- ular value of F( o ,w) on 8. Lemma 2.1 and Lemma 2.2 give the full rank of the Jacobians of Hand H ( regarding them as real mauices by C = R2), therefore they give two onto linear transformations. This implies that 0 e C" is a regular value of H on [ C” x [0,1)]xW and of H on [(C" - {0} )x (0,1)]xW. Also direct computation gives that 0 is a regular value of H on [(C" - {0} )x 0 ]xW. Hence the transversality theorem gives: 29 Lemma 2.5. For almost every (a,b) e C"2 xC", 0 e C" is a regular value both of H( - , . ,a,b) on C"x[0, 1) and ofH(o,-,a) on (c" -{0})x[0,1) . As soon as Lemma 2.5 is established, the rest of work is same as [3] . i.e., the first part implies the homotopy curves are one-dimensional manifolds, also gs!- > 0 where s is the arclength; the second part guarantees the curves not going to the infinity before t -) l. The degree theory argument guarantees all the solutions of P are the end points of these curves. Therefore we have: Theorem 2.6. For almost every (a,b) e C"2 xC", the solution set of H(Z,t) = 0 forms d = dlx - - - xdn_1 x2 one-dimensional homotopy curves beginning with d dis- tinct roots of H(Z,0) = 0 which are easily obtained and leading to all solutions of H(Z,1) = P(Z) = 0 with each curve reaching one zero of P(Z) or approaching the infinite ( this occurs if the number of zeros of P(Z) is less than d including multipli- city) when t —-> 1. Observe that if Z is a solution of (2.1), so is -Z. Thus only half of the curves are needed to be followed. We can pick the beginning points by choosing 21 coordinate zero or one of (d1 -1)-th roots of b1, , z, coordinate zero or one of (d3 -1)-th roots of b3, zs+1 coordinate one of dHI-th roots of b,“ , , z,,_1 coordinate one of dn-1 -th roots of b,,_1, but choosing 2,, coordinate only the positive square root of b,, ( or only the negative square root of b,. ). As soon as half of the solutions of P(Z)=0 are obtained, the another half are just negative of them. The solutions of the numerical bifurcation should be real numbers, and the homo- topy method gives the complex numbers, hence we only pick those solutions of the homotopy function with imaginary part zero theoretically and near zero numerically. CHAPTER 6 THE SINGULAR VALUE DECOMPOSITION AND REGULAR POINTS § 6.1. A Matrix Result Let f: R’"+1 —-)R”', f(xo) = 0, f 6 C2. Suppose x0 is aregular point of the solution set. i.e., Df(xo) has full rank. The question is how to find some points near x0 in the solution set, on a one-dimensional smooth manifold. Let g be the unit tangent vector along this one-dimensional smooth manifold. Using a limit procedure as in Chapter 4, we have: f(x)=0 => §l=o vl,m+l i.e., g: (1.1) Vm+1,m+l . Let f(x) = o, Df(x) = usz with rank m. '21:. . . . aft ' V1.1 ° ° ° V1,m+1 8x1 31ml V= Df(x)= Vm+1,1 ' ' ' Vm+l,m+l ’ 3:11 _ . . afm axl axm+l J 30 31 Let le‘m+l l=1smsilnx1|vi'm+ll it 0. Denote Df (x) as the matrix from Df(x) by 1 + deleting one column : 2a an art aft“ - 3X 1 311-1 31m 316nm Df(x)= (1.2) an. an an. 8f... bit—1 . . . axii-1 311ml axm-I-l and we have: Theorem 1.1. Among all m x m submatrices of Df(x) , [if (x) is the only submatrix which is always nonsingular. Proof: First we prove that it is nonsingular. From f(x) = 0, we have by (1.1) ' 1 ' 'l V 1,111+] Bfl vl,m+l .. vk-l +1 dxk 0 =Df(x) ~~ =Df(X) ,k 1"" l + Vimn --- Vm+1.m+1 ‘2 37* 8f... _Vm+l,n1+ld _ axk J P - r- - V l,m+l 3f] . . . -vk.m+1$ " vk—l,m+l k Df(x) vk+i n+1 = O . .8f (13) ' m . . . ’vk,m+1 — _vm+l,m+l L axk .. vkmm at 0 means as soon as vhmsfl is known, the rest of the components of ( v1,m+19 . . . ,vmr1,m+1)r are known since the null space is one-dimensional. That is equivalent to say (1.3) is uniquely solved, i.e., 15f (x) is nonsingular. L‘— 32 To prove that it is the only m x m submatrix which is always nonsingular, we V 1.m+1 pick the case when § = - Vm +1,m+l Df(x) =U z VT r 0 '1 0 l --- k-th . Therefore 0 . 0 . ' T ”1.1 V1,m OI Vk—IJ vk-IJII 0 o o 1 Vk+1.l Vk+1,m 0 me+1,1 vm+l,m 0 V1.1 viz—1,1 0 Vk+l,l Vm+1.1 V1,». viz-1,111 0 vk+l,m vm+l,m 0 . o . 0 l 0 . . . 0 Thus any in x m submatrix deleting the i-th column of Df(x) with i 3* k must include the column 0 0 0 omO 1 which is the zero column. Thus the theorem is proved. Cl 33 By this theorem, among all m+1 coordinate directions, this is the best direc- tion for keeping the derivative invertible. § 6.2. A Newton’s Method via SVD We may use the following method to find some point near x0 on the solution set of f. Method 2.1. Compute SVD of f at x0. If k is the coordinate index such that ISiSm+l set predictor: x‘O) =xo =h -—-k-th (2.1) ozo~ogo for some small real number h>0 . Obtain x‘ = lim x(") by keeping the k-th coordinate fixed xiv) =x10), and v—-)oo changing the other coordinates from Pxfiv+l)_xiv) ' ~ \NU;:V) Df (1M) :égl)_:’{{): = ‘f (1M) (22) vH v XIn+1)-x§nli .. .I 34 forv =0,1, 2, . If f( x0 ) = 0, in the neighborhood of x0, f(x) is small. And 15f (x0) is non- singular. f e C2 implies Df is Lipschitz in that subspace. Therefore the following theorem can guarantee the above convergence quadratically. After a few iterates, x' can be approximated by high accuracy. Theorem 2.2 ( Newton-Kantorovich [9] ). Given g : Q t: R’" —>R”‘ and the convex set (21 c: Q , let g be continuously differentiable on 91 and satisfy the con- ditions: (a) ll Dg (y)- Dg (z) II 57 II y - 2 II for all y, z in (21; (1» II ng‘0’)“gR"’, and suppose f is sm00th enough so that the Taylor Theorem in Chapter 4 holds. Df(x)=U 2 V‘ , where I 01 0 D, 0 21— 0 0 — 0, cm 0 by setting O',+1=O',+2="'=O'm=0. We also assume that there are only finite smooth curves passing each bifurcation point of the solution set of f(x)=0. In this chapter, we assume that the solution set is connected, otherwise we just consider a component of the solution set. We present in this section several numerical methods by using the singular value decomposition which are going to be used in the next section for an algorithm. Theoretically 0,, = 0 is a criterion for a bifurcation point. Numerically we can choose a very small number 8 which is machine dependent, hence 35 36 Om <8 (1.1) is a criterion for the bifurcation point. Since SVD is a reliable way for detecting rank deficiency and near deficiency, the above method should be accurate to decide the bifurcation point. If a point x is detected as a regular point, according to the tangent vector which is obtained by SVD, i.e., the m+1-th right singular vector (v1.,,,+1,v2.,,,+1. - - - ,v,,,+1,,,,+1)T , we can find xk the largest coordinate changing the solution. Consequently Newton’s method, i.e., Method 6.2.2., can be used. Namely choose a suitable positive number h which is dependent on 0’," monotonically, denoted it by: h =g(om) (1.2) such that . 0 . 0 x<°>=x +1; (1) —--kth; (1.3) o and P xsv+l)_xiv) I .. at): v) Df(xv) :jéimfi'fi = -f) (1.4) ASP-x1111 for v = 0,1,2, , where 37 Ea an an aft” .. 8x1 311-1 311m 3Xm+1 Df(xv)= 31a at». at». 6f»- 8x1 311-1 8x“, axm-t-l After a few iterates, we get a point x(v°) for some v0 which is a numerical solu- tion of f(x) = 0 different from x. Note that (1.2) also can be used to avoid missing bifurcation points since by our methods the step size h varies and the h is smaller when x is close to the bifurcation point, especially h is smaller when x is near the bifurcation point. Using Newton’s method via SVD, we can follow the branch to the bifurcation point x. When h is very small, we can detect the bifurcation point by detecting the last singular values along the tangent direction. If a bifurcation point x is detected, using Theorem 4.1.1 and 5.2.6, all the unit tangential directions i can be determined. Deleting one direction which is opposite the direction while the bifurcation point is obtained, we have points on the other branchs near the bifurcation point numerically: := x +8§, (1.5) where 8 is a small positive number. Now connecting all methods above, we have a numerical method of global bifur- cations which is described in the next section. § 7.2. A Numerical Method of Global Bifurcations via SVD We consider an open bounded region. So if the solution set is bounded in this region, then the algorithm which we will describe below gives all branchs numerically. 38 If the solution set is unbounded, or the set is too large comparing the open region we set, then the algorithm gives the subset of the solution inside the region numerically. A suitable chosen region will give the major character of the solution set. The method begins at one point of the solution set in the region. Usually if the problem has a trivial solution, we will take it with some value of the parameter as a starting point which is a regular point. Then trace the solution curve to get the global behavior. Algorithm 2.1. Given f(x)=0 with f: R’"+1 —) R’" and the hypotheses in the beginning of this chapter, the following algorithm gives the solution set numerically. Procedure Solution_Via_SVD begin define a bounded open region B; initialize x=(x1, . . . ,xm+1 )T ( * in B and fails to satisfy (1.1) am < e * ) and two directions ( * based on the last right singular vector of Df * ); L=2 ( * L is the label of branchs to be followed * ); while L>0 do L:=L-1; while 0‘," 2 e in (1.1) and x in B do generate x1 ( * based on (1.2),(1.3),(1.4) * ), set x:= x1 ; end while if x in B then case (to judge if the bifurcation point x has been met before) begin x is a new bifurcation point: find points with directions on N branchs to be followed ( * based on (1.5) * ), set L:=L+N; x is a previously found bifurcation point: then the branch to this point x determines a branch with opposite direction emanating from x to the starting point. Delete this branch that emanates from x from the list of branches to be followed when x is used as a starting point; Set L: = L-l; ( * See Figure 1 * ) end case end if end while end procedure ; 39 H 11 x2 I'IN lel For example if path 3 starting at x1 leads to x2, then the path 5 leading from x2 to x1 is deleted when x2 is used as a starting point Figure 1. An illustration of the algorithm. 40 § 7.3. An Example We give an example to show the above algorithm. This example is obtained from the bifurcation problem in Banach space of codimension 4, reducing into the null space by the Liapunov-Schmidt method which has been theoretically done in [1] . Let 1:112 xR4 —)R2 bywx (tt,v,ot1,ot2 )—->0with C(w) + ot1 1.1 w + 0:2 1.2 w = 0 , where 3 2 w1 w1 +11th2 __ 10 __00 W: [W2], C(W)= [vwfivz + w% ]' L1 — [01] and L2 — [01]. By selecting 3 parameters, we get an equivalent one-parameter problem , and so the Algorithm 2.1 can be used to get the global bifurcation numerically. Figure 2 is the computer graph of this example. The value on the horizontal axis represents the value of the parameter a1, and the value on the vertical axis represents the sum of the values of WI and w; . The following table gives five different choices of 11 , v, a2 ( for the reason of such a choice also see [1]). For all such choices beginning with( wl , W2 , 011) = ( 0.0, 0.0, -1.5 ), five different diagrams are obtained. case (a) (b) (C) (d) (c) tt 2.0 0.75 0.5 0.5 1.5 012 1.0 1.0 1.0 1.0 1.0 41 3- 3 j 2 2 I o a}; : -r -1 ~ ; -... -2 § -3 —3 I -2 -1 O 1 2 3 4 -2 -1 0 1 2 3 4 (a) (b) 3_ a p ‘ ’1' ° A -1 -.. § ...3' Figure 2. Five computer graphs for the example. LIST OF REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] LIST OF REFERENCES Chow,S.N. and Hale,J.K., Methods of Bifurcation Theory, Grundlehren 251, Springer-Verlag, New York, 1982. Chow,S.N., Mallet-ParetJ. and Yorke,J.A., Finding zeros of maps: homotopy methods that are constructive with probability one, Math Comp, 32, p 887- 899, 1978. Chow,S.N., Mallet-Paret,J. and Yorke,J.A., A homotopy method for locating all zeros of a system of polynomials, Functional Differential Equations and Approx- imation of Fixed Points, H.-O. Peitgen and H.O.Wa1ther eds., Lecture Notes in Math. , 730, p 77-88, Springer -Verlag, New York, p 77-88, 1979. Dongarra,J.J., Moler,C.B., BunchJR. and Stewert,G.W., Linpack Users’ Guide , SIAM, 1979. Garbow,B.S., Boyle,J.M., Dongarra,J.J. and Moler,C.B., Matrix Eigensystem Routines - Eispack Guide Extension , Springer Verlag, New York, 1983. Golub,G.H. and Van Loan,C.F., Matrix Computations , John Hopkins Univ. Press, Baltimore, 1983. IMSL Library, Fortran Subroutines for Mathematics and Statistics, User’s Manual, Edition 9.2, IMSL, November, 1984. Keller,H.B., Numerical solution of bifurcation and nonlinear eigenvalue prob- lems, Applications of Bifurcation Theory , Robanowitz ed., Academic Press, New York, p 359-384, 1977. Stoer,J. and Bulirsh,R., Introduction to Numerical Analysis , Springer-Verlag, New York, 1980. [10] Walson,L.T., Numerical linear algebra aspects of globally convergent homotopy methods, SIAM Review , V0128 (4), p 529-545, 1986. 42