Numerical Solution of Two-Point Boundary Value Problems Citation White, Andrew Benjamin, Jr. (1974) Numerical Solution of Two-Point Boundary Value Problems. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/R9VN-9C49. https://resolver.caltech.edu/CaltechETD:etd-01312007-163410 Abstract The approximation of two-point boundary-value problems by general finite difference schemes is treated. A necessary and sufficient condition for the stability of the linear discrete boundary-value problem is derived in terms of the associated discrete initial-value problem. Parallel shooting methods are shown to be equivalent to the discrete boundary-value problem. One-step difference schemes are considered in detail and a class of computationally efficient schemes of arbitrarily high order of accuracy is exhibited. Sufficient conditions are found to insure the convergence of discrete finite difference approximations to nonlinear boundary-value problems with isolated solutions. Newton's method is considered as a procedure for solving the resulting nonlinear algebraic equations. A new, efficient factorization scheme for block tridiagonal matrices is derived. The theory developed is applied to the numerical solution of plane Couette flow. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: (Applied Mathematics) Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Applied Mathematics Thesis Availability: Public (worldwide access) Research Advisor(s): Keller, Herbert Bishop Thesis Committee: Unknown, Unknown Defense Date: 13 March 1974 Record Number: CaltechETD:etd-01312007-163410 Persistent URL: https://resolver.caltech.edu/CaltechETD:etd-01312007-163410 DOI: 10.7907/R9VN-9C49 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 421 Collection: CaltechTHESIS Deposited By: Imported from ETD-db Deposited On: 31 Jan 2007 Last Modified: 29 Jul 2024 20:32 Thesis Files Preview PDF (White_AB_1974.pdf) - Final Version See Usage Policy. 6MB Repository Staff Only: item control page

NUMERICAL SOLUTION OF TWO-POINT BOUNDARY-VALUE PROBLEMS Thesis By Andrew Benjamin White, Jr. In Partial Fulfillment of the Requirements For the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 91109 1974 (Submitted March 13, 1974) Acknowledgements The a u t h o r wishes t o thank P r o f e s s o r Herbert B , K e l l e r f o r h i s p a t i e n t guidence and i n v a l u a b l e a s s i s t a n c e i n t h e p r e p a r a t i o n of t h i s t h e s i s . H i s encouragement and enthusiasm were e s s e n t i a l t o t h e completion of t h i s work. The o p p o r t u n i t y t o c a r r y o u t t h i s r e s e a r c h was provided by C a l i f o r n i a I n s t i t u t e of Technology Graduate Teaching A s s i s t a n t s h i p s f o r which t h e a u t h o r i s g r a t e f u l . H e a l s o wishes t o e x p r e s s h i s a p p r e c i a t i o n t o t h e e n t i r e Department of Applied Mathematics, b o t h f a c u l t y and s t u d e n t s , f o r p r o v i d i n g an i n v i g o r a t i n g c l i m a t e i n which t o pursue t h i s work. The a u t h o r wishes t o thank M r s . J a n e t Smith and Mrs. Mary Beth Briggs f o r t h e i r e x p e r t t y p i n g of a d e t a i l e d and t e d i o u s manuscript. ABSTRACT The approximation of two-point boundary-value problenls by g e n e r a l f i n i t e d i f f e r e n c e schemes i s t r e a t e d . A n e c e s s a r y and s u f f i c i e n t c o n d i t i o n f o r t h e s t a b i l i t y of t h e l i n e a r d i s c r e t e boundary-value problem i s d e r i v e d i n terms of t h e a s s o c i a t e d d i s c r e t e i n i t i a l - v a l u e problem. P a r a l l e l s h o o t i n g methods a r e shown t o be e q u i v a l e n t t o t h e d i s c r e t e boundary-value problem. One-step d i f f e r e n c e schemes a r e c o n s i d e r e d i n d e t a i l and a c l a s s of c o m p u t a t i o n a l l y e f f i c i e n t schemes of a r b i t r a r i l y h i g h o r d e r of accuracy i s e x h i b i t e d . S u f f i c i e n t c o n d i t i o n s a r e found t o i n s u r e t h e convergence of d i s c r e t e f i n i t e d i f f e r e n c e approximations t o n o n l i n e a r boundary-value problems w i t h i s o l a t e d s o l u t i o n s . ~ewton's method i s c o n s i d e r e d a s a procedure f o r s o l v i n g t h e r e s u l t i n g nonlinear a l g e b r a i c equations. A new, e f f i c i e n t f a c t o r i z a t i o n scheme f o r b l o c k t r i d i a g o n a l m a t r i c e s i s d e r i v e d . The t h e o r y developed i s a p p l i e d t o t h e numerical s o l u t i o n of p l a n e C o u e t t e flow. Table of Contents TITLE PAGE Introduction Chapter 1. Linear Two-Point Boundary-Value Problems 1. Existence Theory 2. Numerical Methods 3. Convergence for Linear BoundaryValue Problems Chapter 2. Difference Schemes 1. Taylor Series 2. Quadrature 3. Gap Schemes 4. Other Schemes 5. Stability of Triangular Schemes 6. Asymptotic Error Expansions 7. Increased Accuracy Chapter 3. Parallel Shooting 1. Parallel Shooting 2. Method of Complementary Functions 3. Method of Adjoints Chapter 4. Nonlinear Two-Point Boundary-Value Problems 1. Finite Difference Schemes 2. Existence of Solutions 3. Newton's Method 4. Equivalence of Shooting and Implicit Schemes Chapter 5. Block Tridiagonal Matrices 1. Block LU Decomposition 2. A Split Decomposition Chapter 6 . A Numerical Example: Flow 1. I s o l a t e d S o l u t i o n 2 . Numerical Procedure 3 . Numerical R e s u l t s Appendix A References Plane Couette Introduction T h i s t h e s i s d e a l s w i t h t h e a p p l i c a t i o n of f i n i t e d i f f e r e n c e schemes t o two-point boundary-value problems. The assumption i s made throughout t h a t t h e s e boundary-value problems have i s o l a t e d s o l u t i o n s ; t h a t i s , t h e homogeneous, l i n e a r i z e d problem h a s o n l y t h e t r i v i a l solution, The g e n e r a l t h e o r y developed p l a c e s no r e s t r i c t i o n s on t h e form of t h e d i f f e r e n c e e q u a t i o n s , I n Chapter 1, t h e a p p l i c a t i o n of an a r b i t r a r y , c o n s i s t e n t d i f f e r e n c e scheme t o a l i n e a r boundary-value problem i s t r e a t e d . t h e main theorem of t h i s c h a p t e r , Theorem 1.16, In t h e s t a b i l i t y of t h e d i s c r e t e boundary-value problem i s shown t o be e q u i v a l e n t t o t h e s t a b i l i t y of t h e a s s o c i a t e d d i s c r e t e i n i t i a l - v a l u e problem. This a s s o c i a t e d i n i t i a l - v a l u e problem employs t h e same d i f f e r e n c e e q u a t i o n s t o approximate t h e d i f f e r e n t i a l e q u a t i o n , b u t i n i t i a l c o n d i t i o n s r e p l a c e t h e boundary c o n d i t i o n s . From t h i s r e s u l t , i t i s c l e a r t h a t a simple s h o o t i n g method i s , i n f a c t , a s p e c i f i c procedure f o r s o l v i n g t h e d i s c r e t e boundary-value problem. S p e c i a l emphasis i s placed on one-step d i f f e r e n c e schemes i n Chapter 2 , High o r d e r a c c u r a t e d i f f e r e n c e approximations are developed using b o t h Taylor series and t h e i n t e g r a l form of t h e d i f f e r e n t i a l equation, I n p a r t i c u l a r , one-step schemes of a r b i t r a r y o r d e r a r e d e r i v e d which r e q u i r e t h e e v a l u a t i o n of a minimum number of new f u n c t i o n s (e.g d e r i v a t i v e s of A ( t ) , f ( t ) ) . The e q u i v a l e n c e shown i n Chapter 1 i s used t o examine t h e s t a b i l i t y of t r i a n g u l a r d i f f e r e n c e schemes. In Chapter 3 , t h e e q u i v a l e n c e r e s u l t of Theorem 1.16 i s g e n e r a l i z e d t o i n c l u d e a l l p a r a l l e l s h o o t i n g methods. Theorem 3.22 shows t h a t t h e s e methods a r e each a p a r t i c u l a r procedure f o r s o l v i n g t h e e q u a t i o n s d e r i v e d from approximating l i n e a r boundary-value problems. The Method of Complementary F u n c t i o n s i s examined i n d e t a i l a s an example of methods f o r s o l v i n g problems w i t h s e p a r a t e d boundary conditions. The Method of A d j o i n t s i s a l s o c o n s i d e r e d and i t i s shown t h a t t h i s method i s n o t i n g e n e r a l e q u i v a l e n t t o t h e d i s c r e t e boundary-value problem, Nonlinear boundary-value problems a r e d e a l t w i t h i n Chapter 4. The d i f f e r e n c e schemes examined i n Chapter 2 a r e g e n e r a l i z e d t o be a p p l i c a b l e t o nonlinear d i f f e r e n t i a l equations. Following Keller [ 6 1, e x i s t e n c e and uniqueness of t h e s e d i s c r e t e approximations i s shown. W e n o t e t h a t ~ e w t o n ' smethod converges q u a d r a t i c a l l y . Chapter 5 i s concerned w i t h t h e p r a c t i c a l problem of s o l v i n g t h e systems of a l g e b r a i c e q u a t i o n s a r i s i n g from t h e approximation of boundary-value problems w i t h s e p a r a t e d boundary c o n d i t i o n s . e q u a t i o n s a r e w r i t t e n i n b l o c k t r i d i a g o n a l form, Mx = b. These The s p e c i a l z e r o s t r u c t u r e of t h i s system i s e x p l o i t e d t o show t h a t , w i t h a n a p p r o p r i a t e row s w i t c h i n g s t r a t e g y , such a m a t r i x p o s s e s s e s a simple b l o c k LU decomposition i f and o n l y i f M i s n o n s i n g u l a r . A numerical example i s p r e s e n t e d i n Chapter 6. c o n s i d e r e d model p l a n e C o u e t t e flow. The e q u a t i o n s The Gap4 scheme, a s d e r i v e d i n Chapter 4, i s used t o d i s c r e t i z e t h e n o n l i n e a r boundary-value problem and ~ e w t o n ' smethod i s employed t o s o l v e t h e r e s u l t i n g s e t of nonlinear equations. A c o n s i s t e n t e f f o r t i s made t o u s e o t o r e p r e s e n t z e r o o r a zero v e c t o r and O ' f o r z e r o m a t r i c e s , e x c e p t i n t a b l e s o r e q u a t i o n numbers. The numbering of theorems, e q u a t i o n s , o r t a b l e s i s done c o n s e c u t i v e l y throughout e a c h c h a p t e r . Chapter I L i n e a r Two-Point Boundary-Value Problems 1. E x i s t e n c e Theory W e c o n s i d e r t h e system of n f i r s t - o r d e r , linear ordinary d i f f e r e n t i a l equations: where u , f, 6 a r e n-vectors and A,B , B 0 1 a r e n x n matrices. proceeding t o t h e numerical approximation of ( l . l a , b ) , Before we p r e s e n t a n e x i s t e n c e and uniqueness r e s u l t convenient f o r our purposes. Theorem 1 . 2 . Let A ( t ) E c m [ o , l ] f o r some m 5 o. Define t h e fundamental m a t r i x ~ ( t a) s t h e s o l u t i o n of X' ( t ) Then f o r each f ( t ) - A ( t ) X(t) = o E unique s o l u t i o n y ( t ) Proof: X(O) = I. c m [ o , l ] and 6 E E", E problem ( l . l a , b ) cm+l[0 ,l] i f f [BO+BIX (1) ] i s n o n s i n g u l a r The s o l u t i o n t o t h e i n i t i a l - v a l u e problem is i n general has a . Uniqueness f o r t h e i n i t i a l - v a l u e problem i n s u r e s t h a t X(t) i s nonsingular on [ o , l ] . The boundary-value problem ( l . l a , b ) has a s o l u t i o n i f and o n l y i f we can d e f i n e a n n-vector r such t h a t y ( t ) s a t i s f i e s t h e boundary c o n d i t i o n This r e q u i r e s Thus, ( l . l a , b ) has a unique s o l u t i o n i f and only i f [B nonsingular. That y ( t ) 0 E cm+'[o,l] + B1 X(1) 1 i s i s an o b s e r v a t i o n from t h e form of d i f f e r e n t i a l e q u a t i o n ( 1 . l a ) . 2. Numerical Methods Here w e d i s c u s s some s t a n d a r d concepts of numerical a n a l y s i s and develop some n o t a t i o n . (l.la,b), I n approximating t h e s o l u t i o n of wewillemployanetofpoints It,} i=J i=J d e f i n e d on t h i s n e t . f u n c t i o n {v, 1 Each v i i s a n n-vector and w e d e f i n e V such t h a t on [ o , l ] a n d a n e t where V i s an n ( J + l ) - v e c t o r , hi = ti - W e d e f i n e t h e mesh w i d t h s max t i=l,,..,J,andh = l l i l J h \Je f u r t h e r r e q u i r e i-1' o i' t h a t f o r each hi t h e r e e x i s t s a A i E [ E , ~ ] E, > 0 , such t h a t T h i s c o n d i t i o n merely s t i p u l a t e s t h a t E s min h./max hi 5 1 1 and t h u s t h e mesh becomes dense i n [ o , l ] a s ho + o. The norm we w i l l employ i s 1 IV 1 1 max 0SilJ Ilvill and f o r block m a t r i c e s , t h e u s u a l induced norm max l l ~ l l= Ilvll = 1 I IwI 1 I n t h e n = l c a s e , t h i s induced norm e q u a l s t h e maximum a b s o l u t e row sum, however t h i s r e s u l t does n o t g e n e r a l i z e t o n > 1. We may e a s i l y produce t h e upper bound max I I"I I ' osi<J where M ij J j =o i s a n n x n m a t r i x element of t h e b l o c k m a t r i x M. An approximation t o t h e d i f f e r e n t i a l e q u a t i o n ( l . l a ) may be w r i t t e n as The boundary c o n d i t i o n s ( l . l b ) a r e approximated by We d e f i n e t h e t r u n c a t i o n e r r o r T i t o be where y ( t ) i s any s o l u t i o n of t h e d i f f e r e n t i a l e q u a t i o n ( l o l a ) . S i m i l a r l y , we d e f i n e a t r u n c a t i o n e r r o r T 0 associated with the boundary c o n d i t i o n s where y ( t ) i s any s o l u t i o n of (1. l b ) . Note t h a t zero. I n c o n s i d e r i n g the accuracy of approximations (1.5a), w e a r e concerned w i t h y ( t ) a s o l u t i o n of ( l . l a , b ) ro w i l l always b e and we d e f i n e Combining t h e d i s c r e t e approximations (1.5a) and (1.5b), d i s c r e t e boundary-value problem may be w r i t t e n as (1.5b), the More briefly, we will write where B h is an n(J+l) with ro = f3. x n(J+l) matrix and r is an n(J+l)-vector Employing Euler's method to approximate the differential equation (lala), this formulation of the discrete problem becomes This example will be used throughout to illustrate various points of interest. Consistency. The difference approximation (1.5a) is said to be consistent with the differential equation (l.la) iff where y(t) is a solution of (l.la,b). example. If y(t) E ~uler'smethod provides an cP [o,l], p 2 2, then Thus, ~ u l e r ' smethod i s c o n s i s t e n t and we s a y t h a t i t i s f i r s t - o r d e r a c c u r a t e because t h e l e a d i n g term i n t h e t r u n c a t i o n e r r o r expansion is linear i n h i . Order of accuracy. A numerical scheme has an o r d e r of accuracy p > o i f p i s t h e l a r g e s t i n t e g e r such t h a t f o r a l l n e t s w i t h h o 5 H. Stability. Let V be a s o l u t i o n t o (1.8). The d i f f e r e n c e scheme (1.8) i s s a i d t o be s t a b l e i f f t h e r e e x i s t c o n s t a n t s K 0 2 o, and H > o, such t h a t f o r a l l n e t s w i t h ho 5 H. For t h e l i n e a r c a s e , t h i s c o n d i t i o n i s e q u i v a l e n t t o showing t h a t a KO 2 0 and H > o e x i s t such t h a t f o r a l l n e t s w i t h ho 5 H. Convergence. The numerical scheme (1.8) i s covergent iff max Ily(ti) oSi<J - vill 4 0 as h 0 4 0 . For convenience, we d e f i n e t h e n ( J + l ) - v e c t o r Y t o be A s t a n d a r d r e s u l t may now b e s t a t e d , Lemma 1.10 L e t t h e boundary-value solution y(t). problem ( l . l a , b ) L e t t h e d i s c r e t e boundary-value be c o n s i s t e n t w i t h ( l . l a , b ) and s t a b l e . have t h e e x a c t problem Then t h e method i s convergent, that is, Proof: The d e f i n i t i o n of t h e t r u n c a t i o n e r r o r ~ [ y ( t ) ]combined w i t h (1.5a, b ) g i v e s By h y p o t h e s i s , t h e d i f f e r e n c e scheme i s s t a b l e , t h u s t h e r e e x i s t c o n s t a n t s KO 2 0 , H > o such t h a t f o r a l l n e t s w i t h ho 5 Ha Also by h y p o t h e s i s , t h e method i s consistent, that is Hence, max l ly - 'I l = .,i<J I ly(ti) - vil 1 + 0 a s ho + 0 For any p a r t i c u l a r scheme, we would l i k e t o be a b l e t o show t h a t t h e n e t f u n c t i o n approximates t h e e x a c t s o l u t i o n a t t h e n e t points: t h a t t h e method i s convergent. The t r u n c a t i o n e r r o r and o r d e r of accuracy a r e g e n e r a l l y d e r i v e d v i a T a y l o r ' s Theorem, a s was done i n t h e example of ~ u l e r ' smethod. I n o r d e r t o show t h a t ~ u l e r ' smethod i s convergent f o r t h e d i s c r e t e boundary-value problem, we need t o prove t h a t t h e m a t r i x Bh i n (1.9) h a s a n i n v e r s e w i t h uniformly bounded norm a s h -+ 0 scheme would be convergent. o. Then by Lemma 1.10, t h i s numerical However, even i n t h i s s i m p l e c a s e , t h e n o n s i n g u l a r i t y of Bh i s n o t obvious. I n Theorem 1.2 t h e i n i t i a l - v a l u e problem i s used t o prove well-posedness of t h e boundary-value problem. As Keller [ 5 ]h a s shown f o r t h e c a s e of t h e centered-Euler scheme, t h i s approach i s a l s o u s e f u l i n t h e numerical problem. I f w e approximate t h e s o l u t i o n of t h e i n i t i a l - v a l u e problem ~ ( 0 = ) B b) by E u l e r t s method, t h e f o l l o w i n g e q u a t i o n s r e s u l t S i m i l a r l y , we d e f i n e t h e d i s c r e t e i n i t i a l - v a l u e m a t r i x Th a s s o c i a t e d w i t h any % (1.7) t o be To form I h' w e r e p l a c e Bo, B 1 i n t h e f i r s t block row of i3 w i t h h 1,0 respectively. 3. Convergence f o r l i n e a r boundary-value problems -T h e main r e s u l t of t h i s c h a p t e r i s t h a t a numerical scheme a p p l i e d t o a boundary-value problem w i t h a unique s o l u t i o n i s s t a b l e and c o n s i s t e n t i f and o n l y i f t h e a s s o c i a t e d i n i t i a l - v a l u e problem i s s t a b l e and c o n s i s t e n t . It w i l l be u s e f u l t o prove s e v e r a l lemmas first. Lemma 1.12 ( F a c t o r i z a t i o n ) associated with 'I,. Then Let 7 h be t h e i n i t i a l - v a l u e m a t r i x where L is a (J+l)n x n+l matrix N is an (n+l) x (J+l)n matrix and 6 is an (n+l)-vector Proof: We take (1.7) to be the most general form of the boundary-value matrix Bh. The associated initial-value matrix is It may be d i r e c t l y v e r i f i e d t h a t L,N,C a r e t h e proper f a c t o r s , b u t t h e sequence of o p e r a t i o n s below may c l a r i f y t h e d e r i v a t i o n . The i d e n t i t y ro = 6 i s used h e r e . ' Lemma 1.13 (Reducing) Let L,N be pxq and qxp m a t r i c e s r e s p e c t i v e l y . Let x and b be a p-vector and a q-vector r e s p e c t i v e l y . L have r a n k q 5 p. Further, l e t Then t h e m a t r i x e q u a t i o n has a s o l u t i o n i f and o n l y i f Proof: The s u f f i c i e n c y i s t r i v i a l on s u b s t i t u t i o n of Lw i n t o ( 1 . 1 4 ) . S i n c e L i s a pxq m a t r i x w i t h rank q 5 p , t h e columns of L a r e l i n e a r l y independent. Therefore, t h e r e e x i s t s a p x ( P - ~ ) m a t r i x L such t h a t W e may write any v e c t o r x E EP u n i q u e l y a s Using t h i s r e p r e s e n t a t i o n i n (1.14), The v e c t o r x i s a s o l u t i o n of (1.14) i f and o n l y i f and The n e c e s s i t y f o l l o w s . m Note t h a t t h i s lemma r e d u c e s t h e s o l u t i o n of p simultaneous equations t o s o l u t i o n of q e q u a t i o n s , q 5 p , Now w e s t a t e and prove t h e main r e s u l t of t h i s c h a p t e r . Theorem 1.16 L e t t h e boundary-value problem ( l , l a , b ) have a unique s o l u t i o n y ( t ) E 1 c [0,1]. Then t h e f o l l o w i n g a r e e q u i v a l e n t : a ) The d i s c r e t e boundary-value problem i s s t a b l e , c o n s i s t e n t (and c o n v e r g e n t ) . b ) The a s s o c i a t e d i n i t i a l - v a l u e problem i s s t a b l e , c o n s i s t e n t (and c o n v e r g e n t ) . Proof: (b => a ) The F a c t o r i z a t i o n Lemma (1.12) s t a t e s t h a t f o r t h e boundary-value problem (1.8) The i n i t i a l - v a l u e problem i s s t a b l e by h y p o t h e s i s , t h e r e f o r e Ih i s n o n s i n g u l a r f o r a l l ho 5 H, say. L e f t m u l t i p l y i n g (1.17) by Ih-1 w e o b t a i n Define Z and z by where Z i s an n ( J + l ) x n m a t r i x and z i s an n ( ~ + l ) - v e c t o r , From (1.18), we n o t e t h a t Z s a t i s f i e s and t h u s i s a n approximation t o t h e fundamental m a t r i x s o l u t i o n (1-,3), By Lemma 1.10, t h e d i s c r e t e i n i t i a l - v a l u e problem i s convergent max orilJ IIzi - x ( t i ) I I +- o as h 0 -+ o and i n p a r t i c u l a r I IzJ - x ( ~ ) I I o a s ho -+ -+ o . By v i r t u e of ( 1 , 1 8 ) , we may w r i t e (1.17) a s - v + [ z : , z ] € r w 51 = 0 . The odumns of t h e n ( J + l ) x n + l m a t r i x [Z ,1z ] a r e l i n e a r l y independent provided t h a t r, # o, 1 5 i 5 J. st I f t h e ( n + l ) - - column of L i s t h e I z e r o v e c t o r w e remove i t , r e p l a c e degeneracy i s removed, = k] by €J = [e] and t h e We complete t h e proof assuming r j # o, f o r some j . - 1 By t h e Reducing Lemma, (1.20) h a s a s o l u t i o n i f and o n l y i f , satisfies and t h e s o l u t i o n must be of t h e form W e r e c a l l t h a t Zo = I and t a k e A = 1, s o t h a t (1.21) i s e q u i v a l e n t t o Compare t h i s w i t h the c o n d i t i o n de.ri.vc?(l i n Theorem ( 1 . 2 ) f o r the. s o l u t i o n of t h e boundary-value problem ( l . l a , b ) . E x i s t e n c e and uniqueness of t h e s o l u t i o n , V, of t h e d i s c r e t e boundary-value problem now h i n g e on t h e n o n - s i n g u l a r i t y of (B 0 + BIZJ). We w r i t e By h y p o t h e s i s y ( t ) i s unique, and Theorem 1.2 s t a t e s t h a t [Bo + BIX(l)] must be non-singular. t h u s , by t h e Banach Lemma, B 0 It h a s a l r e a d y been shown t h a t + B1Z J i s non-singular f o r a l l n e t s w i t h ho 5 H, provided H i s s o s m a l l t h a t f o r some p E ( o , l ) It i s c l e a r t h a t t h e d i s c r e t e boundary-value problem i s s t a b l e s i n c e and Z and z a r e s o l u t i o n s of d i s c r e t e i n i t i a l - v a l u e problems, which were assumed s t a b l e . ( a => b) The proof of t h i s p a r t i s e s s e n t i a l l y t h e same. Consider t h e d i s c r e t e i n i t i a l - v a l u e problem I h V - r = o . Following t h e F a c t o r i z a t i o n Lermna (1.12), w e may w r i t e Th v - r = Bhv - L{NV + 6 1 . By h y p o t h e s i s , t h e d i s c r e t e boundary-value problem i s s t a b l e , t h u s A 8 h i s non-singular and we may d e f i n e [ Z ! 21 by The Reducing Lemma (1.13) now i m p l i e s t h a t t h e s o l u t i o n V of (1.23) must be of t h e form and e x i s t s i f and o n l y i f A Therefore, w e wish t o show t h a t Z h 0 5 H, 0 f o r H s u f f i c i e n t l y small. i s non-singular Theorem 1 . 2 d e r i v e s t h e s o l u t i o n of t h e boundary-value problem t o be where r must s a t i s f y Thus, t h e s o l u t i o n t o t h e boundary-value problem A 2' ( t ) - A ( t ) ~ ( t =) o BoX(o) + BIX(l) i s given by f o r a l l n e t s with - I = 0 where R i s a n n x n m a t r i x and must s a t i s f y By h y p o t h e s i s y ( t ) i s u n i q u e , t h u s t h e m a t r i x R i s n o n - s i n g u l a r , A A t h a t i s , X(o) i s n o n - s i n g u l a r . W e n o t e t h a t Z a s d e f i n e d i n (1.24) i s a c o n v e r g e n t approximation t o t h e e q u a t i o n s (1.26), hence A and, a s b e f o r e , Z 0 is non-singular f o r a l l n e t s with h H i s s u f f i c i e n t l y small. 0 5 H, where Now, and t h e d i s c r e t e i n i t i a l - v a l u e problem h a s a s t a b l e s o l u t i o n . m W e r e c a l l t h a t t h e o r d e r of a c c u r a c y of a scheme i s determined by t h e l a r g e s t i n t e g e r p such t h a t I l ~ [ ~ ( t1 )1 ]5 K h: for all h 0 5 H. C o r o l l a r y 1.27. L e t t h e boundary-value problem ( l . l a , b ) solution y(t) C P + ~ [ O , ~ f] o, r some i n t e g e r p 2 1. E a p p r o x i m a t i o n s (1.5a,b) b e p-order accurate, have a u n i q u e Let t h e d i s c r e t e Further, l e t the d i s c r e t e i n i t i a l - v a l u e problem a s s o c i a t e d w i t h (1.5a,b) b e s t a b l e . Then f o r a l l n e t s w i t h ho 5 H , where V i s t h e n ( J + l ) - v e c t o r s o l u t i o n of (1.5a,b) and k i s a constant independent of h o e Proof. The proof f o l l o w s t h a t of Theorem 1.16. These r e s u l t s and t h e m a n i p u l a t i o n s l e a d i n g t o them. s u g g e s t several further topics, Chapter 11 w i l l b e concerned w i t h t h e o p e r a t i o n a l c h a r a c t e r i s t i c s of v a r i o u s s p e c i f i c n u m e r i c a l schemes, Theorem 1.16 w i l l b e e x p l o i t e d t o d e t e r m i n e t h e s t a b i l i t y p r o p e r t i e s of a p a r t i c u l a r c l a s s of schemes. W e w i l l examine a s y m p t o t i c e r r o r expansions and d i s c u s s methods of i n c r e a s i n g t h e a c c u r a c y of d i s c r e t e approximations. I n Chapter 111, we w i l l e n l a r g e upon Theorem 1.16 and t h e e q u i v a l e n c e of d i s c r e t e i n i t i a l - v a l u e and boundary-value problems. D r . H , 0 , Kreiss h a s r e c e n t l y p u b l i s h e d a paper [ 7 ] d e a l i n g w i t h t h e s t a b i l i t y of d i s c r e t e approximations t o a r b i t r a r y o r d e r systems of l i n e a r o r d i n a r y d i f f e r e n t i a l e q u a t i o n s . Thus, i t a p p e a r s t h a t t h i s c h a p t e r i s perhaps a s p e c i a l c a s e of Kreisst r e s u l t s . However, i n Appendix A, we g i v e a n i n d i c a t i o n t h a t Kreisst complementary t o o u r own and n o t i n c l u s i v e . work i s F u r t h e r , Theorem 1.16 g i v e s a n e c e s s a r y and s u f f i c i e n t c o n d i t i o n f o r s t a b i l i t y of t h e most g e n e r a l schemes, whereas D r . K r e i s s d e a l s w i t h k - s t e p methods. Chapter I1 D i f f e r e n c e Schemes I n t h i s c h a p t e r , w e w i l l examine some s p e c i f i c d i f f e r e n c e schemes; of p a r t i c u l a r i n t e r e s t a r e one-step (two-point) methods. One-step methods f o r l i n e a r problems have t h e c h a r a c t e r i s t i c form These methods a r e of s p e c i a l i n t e r e s t f o r s e v e r a l r e a s o n s . nonuniform n e t s and, a s K e l l e r solutions. They admit 1 5 1 h a s shown, p i e c e w i s e smooth I n addition, t h e c a l c u l a t i o n s involved i n solving t h e matrix equation f o r s e p a r a t e d boundary c o n d i t i o n s a r e p a r t i c u l a r l y simple. In t h i s c a s e , t h e m a t r i x 8 may b e p u t i n t o a b l o c k t r i d i a g o n a l form which h p o s s e s s e s a s i m p l e LU-decomposition. One of two methods i s u s u a l l y employed t o g e n e r a t e d i f f e r e n c e approximations t o ( l o l a ) and e v a l u a t e t h e i r o r d e r of accuracy. The f i r s t approach i s v i a Taylor s e r i e s and n e c e s s i t a t e s e v a l u a t i n g t h e t r u n c a t i o n e r r o r f u n c t i o n s a t some r e f e r e n c e p o i n t t ( t ) where R i Quadrature formulae may a l s o be used t o g e n e r a t e u s e f u l approximations. I f we w r i t e (1. l a ) i n t h e e q u i v a l e n t form t h e a p p l i c a t i o n of q u a d r a t u r e i s immediately obvious. 1. Tavlor S e r i e s The centered-Euler o r Box-scheme ( K e l l e r 1 : 5 ] ) i l l u s t r a t e s t h e u s e of Taylor s e r i e s . where t i-112 Applied t o ( l o l a ) t h e c e n t e r e d E u l e r scheme - -hi - ti 2 I f we assume t h a t y ( t ) E 2nrfl [ o , l ] and l e t t R ( t i ) C = ti-1l2, t h e t r u n c a t i o n e r r o r i s shown t o be 2 and the l e a d i n g term, 0 (hi) , i s Here t h e c h o i c e of tR - - ti-1/2 i s i m p o r t a n t due t o t h e symmetry of (2.2). However, t h e f i r s t term of t h e t r u n c a t i o n e r r o r expansion, (2.3b), i s i n v a r i a n t and hence t h e o r d e r of accuracy of any scheme must b e i n v a r i a n t t o t h e r e f e r e n c e p o i n t t tR R. I f we were t o u s e - ti-l' the t r u n c a t i o n e r r o r i s shown t o b e and t h e l e a d i n g term i s However, r e c a l l i n g t h a t y ( t ) i s an e x a c t s o l u t i o n of (1. l a ) , t h i s e x p r e s s i o n i s shown t o b e e q u i v a l e n t t o t h e f o l l o w i n g Thus, t h e l e a d i n g term of (2.3a) t~ = t i-1/2 and t R = i s e x a c t l y t h e same a s (2.4a) w i t h 5-1 r e s p e c t i v e l y . Note t h a t t h e n e x t terms i n each expansion a r e n o t t h e same. The centered-Euler scheme may b e g e n e r a l i z e d t o an a r b i t r a r i l y h i g h o r d e r by employing t h e e x a c t form of t h e t r u n c a t i o n e r r o r . 2 c o r p o r a t i n g t h e f i r s t term, O(hi), In- of t h e t r u n c a t i o n e r r o r expansion Act) and f ( t ) , we must a l s o e v a l u a t e A") ( t ) , A(') ( t ) , )'(f (t ) , f ( l ) ( t ) and v a r i o u s products of t h e s e f u n c t i o n s . Since t h e f u n c t i o n s C(t) and g ( t ) w i l l continue t o appear, we need t o d e f i n e t h e s e q u a n t i t i e s i n more g e n e r a l i t y . The u s u a l procedure w i l l b e t o e v a l u a t e y ( p ) ( t ) i n terms of d e r i v a t i v e s and combinations of A ( t ) , f ( t ) and y ( t ) . For s i m p l i c i t y , w e d e f i n e where C ( t ) i s an n x n matrix and g ( t ) i s an n-vector. v V Thus, Y ( 0 ) ( t ) = y ( t ) => Co = I, go = 0 . The terms of p a r t i c u l a r i n t e r e s t a r e contained i n t h e following t a b l e . Table 2.7 Employing t h i s n o t a t i o n , t h e fourth-order g e n e r a l i z e d centered-Euler scheme (2.5) i s I n t h e n e x t s e c t i o n , t h i s method w i l l b e improved upon i n t h e s e n s e t h a t t h e f u n c t i o n s Co, C1, C 2 , and C 3 w i l l b e used t o d e f i n e a s i x t h - order-method. 2. Quadrature Employing (2.1) we g e t where y ( t ) is t h e e x a c t s o l u t i o n t o ( l . l a , b ) . t r a p e z o i d a l r u l e t o (2.9), I f we a p p l y t h e the result is Thus a second o r d e r a c c u r a t e d i f f e r e n c e scheme ( ~ r a p e z o i d a l~ u l e )may be d e f i n e d by This may b e e q u i v a l e n t l y w r i t t e n a s u s i n g t h e f u n c t i o n s i n Table 2.7. This one-step method and o t h e r s of a r b i t r a r i l y h i g h o r d e r may b e d e r i v e d i n a more g e n e r a l s e t t i n g by means of (2.9) Lemma 2.11. . Let b ( t ) E cdl[a,b]. Also l e t polynomials p i ( t ) be d e f i n e d by wh.er e 5, E [arb]. Then + (-l)&l Proof: (s) b(&') ( s ) ds , Repeated i n t e g r a t i o n by p a r t s y i e l d s ( a ) immediately. I n d u c t i o n on v y i e l d s (b). The i n d u c t i o n h y p o t h e s i s i s s t a r t e d by (2.12a) and, f o r m a l l y , i s IP,-,(t)l 5 Ib-al v-1 Thus, Innumerable schemes p r e s e n t themselves by way of t h i s Lemma. The Euler-Maclarin sum formula may b e d e r i v e d u s i n g polynomials which look l i k e where the s c a l e f o r each polynomial i s n e c e s s a r i l y d i f f e r e n t . Thus, t h e Euler-Maclarin formula may b e used t o d e f i n e a d i f f e r e n c e scheme of a r b i t r a r i l y high o r d e r . A f i f t h - o r d e r method d e f i n e d i n t h i s manner i s One u n d e s i r a b l e a s p e c t of (2.13) i s e v a l u a t i n g c 4 ( t ) and g ( t ) . 4 o r d e r t o employ (2.13) w e must e v a l u a t e (t) , In ( t ) , A(') ( t ) , A(t) , f (3) ( t ) , f C 2 )( t ) , f ( I ) ( t ) , f ( t ) and a s s o r t e d p r o d u c t s and sums of these functions. 3. Gap schemes Thus f a r , w e have c o n s i d e r e d s e v e r a l one-step schemes. The centered-Euler scheme w a s g e n e r a l i z e d i n S e c t i o n 1 t o a f o u r t h o r d e r a c c u r a t e scheme. This i n c r e a s e i n accuracy i s bought a t t h e p r i c e of e v a l u a t i n g AC2)( t ) , A(') ( t ) , f (2) ( t ) , and f ( I ) ( t ) i n a d d i t i o n t o A(t) and f ( t ) . I n S e c t i o n 2, t h e Euler-Maclarin Sum Formula was used t o d e f i n e a f i f t h - o r d e r method, b u t t h i s d i f f e r e n c e scheme r e q u i r e d e v a l u a t i o n of ( t ) , f (3) ( t ) and a l l lower o r d e r d e r i v a t i v e s . S i n c e each new f u n c t i o n t o b e e v a l u a t e d i n c r e a s e s b o t h programming complexity and computation time, we want t o examine h i g h - o r d e r , one- s t e p methods w i t h a minimum of d e r i v a t i v e s r e q u i r e d . The centered-Euler scheme and t h e T r a p e z o i d a l Rule (2.10) possess another u s e f u l property. Each of t h e s e d i f f e r e n c e schemes h a s a t r u n c a t i o n e r r o r expansion which c o n t a i n s o n l y even powers of 1 1 is For r e a s o n s t h a t w i l l become c l e a r , s e e S e c t i o n 6 on asymptotic e r r o r expansions, w e want t o p r e s e r v e this p r o p e r t y i n c o n s i d e r i n g higho r d e r methods. R e c a l l i n g Lemma 2.11, we have the r e s u l t b + (-l)&l f P&l (s) b(&') ( s ) ds. a V=m I f t h e sequence of polynomials ( p V ( t ) lvmo - could be d e f i n e d i n such a way t h a t . v=n+l,. , ,2n 13y(a) = pV(b) = 0 t h e n a d i f f e r e n c e method of o r d e r 2n can b e d e f i n e d which r e q u i r e s t h e e v a l u a t i o n of A - ) ( 1 , f("-') ( t ) , and a l l lower-order d e r i v a t i v e s . For i n s t a n c e , a f o u r t h - o r d e r method s o d e f i n e d r e q u i r e s e v a l u a t i n g ~("(t) and f(') needs ~ ( " ( t ) , Lenrma 2.14 ( t ) where a s t h e g e n e r a l i z e d centered-Euler scheme f(2)(t), etc. L e t pZn(t) = 1 (t-a) n (t-b) polynomials {py ( t ) 3, V 2 0 , such t h a t n . Define t h e sequence of Then ( p v ( t ) ) s a t i s f y conclusion a) Proof: a) of Lemma 2.11. We can r e p l a c e (.2.15a) w i t h and t h e sequence s a t i s f i e s conclusion a ) of Lemma 2.11. 5 Ib-alZn and t h e r e s u l t i s immediate by b. i n d u c t ion. c. Immediate from (2.15a). d. To prove t h i s p a r t of t h e Lemma, i t i s s u f f i c i e n t t o show t h a t p 2 r ( t ) , r 2 n , i s an even f u n c t i o n about t = b+a b+a l o s s of g e n e r a l i t y , we assume t h a t - - 0 , thus tZr= 0 . 2 hypothesis, p 2n ( t ) i s even about 2v-1 BY F,hence we assume t h a t p2 (v-1) i s an even function. Therefore, p Without 2 - ( t ) i s an odd f u n c t i o n . (t> 3 n t-l I w Recalling t h e d e f i n i t i o n of P ( t ) , - f o r L < n, P -(L) ( t i ) = 0. Thus, - Define 8 k - (-1) m! 2n! (m-n) ! (2n-k) ! 2n!(n-k)! (2n-m) n! ! (hi)2n-m* " I and (2.16) becomes k! Table 2.17 A s examples, w e w r i t e out t h e fourth-order and s i x t h - o r d e r schemes. Fourth-order scheme: Sixth-order scheme : The m a t r i c e s C V and t h e v e c t o r s gv a r e d e f i n e d i n Table 2.7. D r . K e l l e r n o t e d t h a t t h e s e gap schemes might b e d e r i v e d u s i n g Hermite i n t e r p o l a t i o n t o e s t i m a t e t h e f u n c t i o n b ( t ) = A(t)y(t) + f(t) . For the f o u r t h - o r d e r gap scheme ( 2 . 1 8 ) , we w r i t e b ( t ) = H(t) + (b(t) - H(t)) where B ( t ) i s t h e lowest o r d e r polynomial t h a t matches b (ti) , b(') b(ti - and b0(tiel). I n this c a s e , (ti) , Recalling t h e i n t e g r a l equation (2.1) , we have ff(t) may be i n t e g r a t e d e x a c t l y t o g e t Since b ( t ) = J1\t), t h i s formulation w i l l y i e l d a d i f f e r e n c e scheme e x a c t l y t h e same a s (2.18) w i t h t h e t r u n c a t i o n e r r o r 4. Other Schemes We a r e p a r t i c u l a r l y i n t e r e s t e d i n one-step d i f f e r e n c e schemes of t h e form (2.2) , (2.5) , and (2.18) , t h u s only a b r i e f mention w i l l b e made of o t h e r methods. The i n t e g r a l e q u a t i o n (2.1) s u g g e s t s k-step methods of t h e form Note t h a t (2.20) r e p r e s e n t s (J-k+l)n e q u a t i o n s i n ( J + l ) n unknowns, and kn more e q u a t i o n must b e found t o determine V uniquely. These e x t r a e q u a t i o n s are given by t h e u s u a l boundary c o n d i t i o n s (1.5b) and by I1 s t a r t i n g " m-step d i f f e r e n c e schemes, m < k , u s u a l l y of lower o r d e r than (2.20). I n many c a s e s t h e " s t a r t i n g " scheme i s a one-s t e p scheme. Combining a l l of t h e s e e q u a t i o n s i n t o t h e form t h e i n i t i a l - v a l u e matrix, I1 Ih, i s lower t r i a n g u l a r . I f t h e one-step s t a r t i n g " scheme is of lower o r d e r than t h e k-step method, t h e mesh s i z e , h ( l ) , a s s o c i a t e d w i t h t h e one-step method should b e s m a l l e r t h a n t h a t of t h e k-s t e p method, h(k) . I f t h e methods a r e of o r d e r p and q r e s p e c t i v e l y (p < q) , then C o r o l l a r y 1.27 i n d i c a t e s t h a t h ( 1 ) and h(k) should b e r e l a t e d by 5. S t a b i l i t y of T r i a n g u l a r D i f f e r e n c e Schemes A l l of t h e methods mentioned h e r e and most of t h o s e commonly used g i v e r i s e t o b l o c k triangular initial-value matrices Ih . That i s , they may be w r i t t e n i n the form The f a c t o r of J appearing w i t h a l l M ij i n v e r s e of t h e mesh s i z e . 's is a n o r m a l i z a t i o n of t h e We would l i k e t o examine t h e s t a b i l i t y of f a m i l i e s of m a t r i c e s 1 of t h i s form. h The s t a b i l i t y theory f o r d i s c r e t e i n i t i a l - v a l u e problems i s well-developed provided t h a t a k-s t e p d i f f e r e n c e scheme i s used and a l l Mij scalar. can be w r i t t e n a s M. . I , M a 1J ij We a r e i n t e r e s t e d i n a r e s u l t which does n o t r e q u i r e t h e s e two r e s t r i c t i o n s on I he Lemma ( 2 . 2 2 ) . L e t Ih be a m a t r i x of t h e form of (2.21) Dv, Nv, by V=0,1, . Define L e t Do b e n o n s i n g u l a r and l e t I IDo-1I I = d1 , I I D ~ I I= dl, 1 1 ~ I ~ = 1 no, 0 and I IN1/ 1 = nl, independent of J. Further, l e t 0 d < 1. Then 0 Proof : T h e matrix 1 i s w r i t t e n a s h or, equivalently where and By h y p o t h e s i s and the Banach Lemma, D J>Jo=- I do , and 0 + J-1DI i s n o n s i n g u l a r f o r Since (Do + J-'D~)-' i s block diagonal and (No + J - h l ) n i l p o t e n t , t h e product (Do Thus, where i s block + J - ~ D ~ )(N- o~ + J m b l ) i s block n i l p o t e n t . Combining Corollary 1.27 and Lemma (2.22) , we s t a t e t h e following Theorem without proof, Theorem 2.24. L e t t h e l i n e a r boundary-value problem ( l . l a , b ) 1 unique s o l u t i o n y ( t ) E C pth-order L e t Bh be of t h e form accurate. [O,l]. have a Also, l e t t h e d i f f e r e n c e scheme be -1 1 . ..,J and I 1 Mvv 11 52 , is nonsingular f o r every i=O , where M ii - I I M I ~5 dl,~ ~independent of J. n 0 and - 5 d 1. Then f o r h 0 0 0 Further, l e t < H, H s u f f i c i e n t l y s m a l l , - This Theorem gives a p r e c i s e e s t i m a t e on the convergence of finite difference solutions. However, we now want t o look a t a more d e t a i l e d accounting of t h e e r r o r i n c u r r e d by approximating ( l . l a , b ) w i t h a d i f f e r e n c e scheme. 6. Asymptotic E r r o r Expansions I n S e c t i o n 3, t h e gap schemes were derived s o t h a t t h e t r u n c a t i o n e r r o r expansion would have only even powers of h i* The advantage of t h i s arrangement w i l l be e v i d e n t i n t h i s d i s c u s s i o n . The g e n e r a l d i f f e r e n c e scheme is defined a s The t r u n c a t i o n e r r o r i s h e r e defined t o be where z(t) e c ~ + ~ [ 0 , 1 ] . For example, E u l e r ' s method has t h e t r u n c a t i o n e r r o r , l e t t i n g - tR - timl' I n general, t h e t r u n c a t i o n e r r o r w i l l be assumed t o be of t h e form = Lz(tR) - ctR) where each f (A,f) i s a l i n e a r d i f f e r e n t i a l o p e r a t o r of a t most o r d e r k k + l which depends, in g e n e r a l n o n l i n e a r l y , upon A ( ~ () t R ) , f ( k ) ( t R ) and a l l lower o r d e r d e r i v a t i v e s . Let t h e d i f f e r e n t i a l equation ( l . l a , b ) have a unique Theorem 2.27. solution y(t) For E u l e r ' s method, E cq + l ~ where ~ , l ]A(t) E cq[O,l] and f ( t ) E cq~0,1~. Let b e a s t a b l e , c o n s i s t e n t d i f f e r e n c e approximation t o (1. l a , b ) . Further, let t h e t r u n c a t i o n e r r o r be given by where z ( t ) E CP+l[O,l] F k ( ~ ( t ), f ( t ) ) z ( t ) E nets with h 0 and tR= ti + O(hi) . Also l e t cv [ 0 , l ] where v = midp-k,q-k} . < H, H s u f f i c i e n t l y small, - Then f o r a l l where yo(t) = y ( t ) and y k ( t ) , k 1, i s defined by Proof: By i n d u c t i o n and t h e c o n t i n u i t y of Fk (A,£), A(t) , and f ( t ) , i t can be shown t h a t Immediately upon s u b s t i t u t i n g i n t o t h e t r u n c a t i o n e r r o r expansion (2.28b), we g e t Since yk(t) s a t i s f y (2.30a), we have immediately and hence The boundary conditions and .ro[w(t) J c l e a r l y y i e l d By h y p o t h e s i s , t h e d i f f e r e n c e scheme i s s t a b l e , thus O r equivalently, The expression (2.29) e r r o r expansion. y(t) E cm[O,l], A(t) = 0 , k=O,l, - y ( t .) i s termed t h e asymptotic 1 I n t h e following d i s c u s s i o n we assume t h a t E cm[O,l], and f ( t ) t r u n c a t i o n e r r o r -ri[z(t)] F2k+l f o r vi E cm[O,l]. has only even powers of h Suppose t h a t t h e i' that is, .... The f u n c t i o n y l ( t ) is defined by Since t h e d i f f e r e n t i a l equation ( l . l a , b ) has a unique s o l u t i o n , Assume t h a t y2k-l(t) defined by Z 0, k=1,2 ,...,v, then y2vtl(t) i s n B + d r-l h) V U h * r-l "4 g , &b? "a z3 r: t g *d U r3 (I] : . " ! g 3 k & 0 k i r l m It a + $ 8*rl h 8 5 - k 0 U k a & 5k " 4 a U a L d 3 a a, rl h U (I] *I4 ) (I] * -rl 2 4 because the h. and h terms a r e 1 i nonvanishing i n (2.32). However, c l o s e r examination of F2 and F4 y i e l d s t h e following i d e n t i t i e s : d - ~ ( A ' ~ ) - ~ A A ( ~ ' - ~ Ad(t~ [Lz ) A ( t R )-f (tR) ] + ~(A(~)A+A~(~)+~A(~)A(~))[Lz( t R )-£( tR) ] 1- Thus, y 2 ( t ) is defined by However, from (2.33) thus, y 2 ( t ) = 0 . i t can b e seen t h a t I n a similar f a s h i o n we can show t h a t y ( t ) = o. 4 Theorem (2.27) and the i d e n t i t i e s ( 2 . 3 3 ) and ( 2 . 3 4 ) show that for t h e Gap6 d i f f e r e n c e scheme where y6(t) i s defined by 7. I n c r e a s e d Accuracy W e w i l l c o n s i d e r two methods of i n c r e a s i n g t h e accuracy of numerical s o l u t i o n s t o t h e d i f f e r e n t i a l equation (1.a , ) . These methods a r e , by name, Richardson h -t o e x t r a p o l a t i o n and FOX'S method of Deferred Corrections (Fox [ 2 1, Pereyra [ 9 1) . Both methods r e l y on t h e e n t i r e t r u n c a t i o n e r r o r s e r i e s rather than t h e o r d e r of accuracy. Richardson e x t r a p o l a t i o n employs a s o l u t i o n of the d i s c r e t e problem on two d i f f e r e n t n e t s t o e l i m i n a t e s u c c e s s i v e e r r o r terms; t h e method of Deferred C o r r e c t i o n s uses a s o l u t i o n t o t h e d i s c r e t e problem t o approximate t h e f i r s t t r u n c a t i o n e r r o r term. To i l l u s t r a t e t h e method of Deferred C o r r e c t i o n s we look a t t h e centered-Euler method (uniform mesh) U [I] k rl '44 Q) k 3 a Q) 0 0 k a [I] hl 6 =rl U I 5 I w n rl I dfd cd Q) U 0 E: 3 .9, cd I-l u I I rld ? f-1 -rl 3 U rlld 3 a [I] 3 C 2 0 d U C) a k k a, n d n rl Q) a U rl I k 0 U rl *rl N 1 a rl d \ rl 0 k k Q) w 8 w 0 U a g d 5 U 3 rl 0 rn 0 8 8 Q) U cd k h ur 50 w k 3 cd Q) a k 0 CV 2 rn disadvantages. At the boundary t=O, the quantity vo must be defined -1 and thus we have to examine the solution outside the domain of interest. Also, the solution procedure for V O is significantly different than (2.38) where the approximation to the first truncation error term is included. Evaluation of higher-order derivatives of y(t) is at best a delicate procedure and must be done with the utmost care. Richardson extrapolation avoids these particular problems. Once two aolutions are calculated, each on a different net, the power structure of the asymptotic error expansion is used to eliminate the leading error term at those points common to both nets. For example, the sixth-order gap scheme has an asymptotic error expansion of the form If we solve the discrete problem twice for {vi (1) 3 and {vi(2) 1 , then Let ti(l) = t.(2), where t (1) is the ith point of the first net and J i t (2) is the jth point of the second net. . J Defh e then This elimination procedure may, of course, be carried on as long as is practical, tkat i s , u n t i l the continuity of y ( t ) or round-off errors make further refinement unproductive. Chapter I11 P a r a l l e l Shooting This c h a p t e r examines t h e r e l a t i o n s h i p between i m p l i c i t f i n i t e d i f f e r e n c e procedures and d i s c r e t e i n i t i a l - v a l u e t e c h n i q u e s . The proof of Theorem 1.16 y i e l d s t h e f o l l o w i n g r e s u l t : I f both t h e boundary-value procedure and t h e i n i t i a l - v a l u e procedure have un'ique s o l u t i o n s , t h e n procedure ( 3 1 ) (3.2a,b,c) vIV = ViBV , i = o , . ..,J. That i s , t h e 1 i s a c t u a l l y a method of s o l v i n g t h e e q u a t i o n s The main r e s u l t of t h i s c h a p t e r s t a t e s t h a t any p a r a l l e l s h o o t i n g t e c h n i q u e i s a l s o e q u i v a l e n t t o (3.1) and, t h u s , t o (3.2a,b,c). P a r a l l e l s h o o t i n g i n c l u d e s such procedures a s simple s h o o t i n g , Method of Complementary F u n c t i o n s , and t h e Godunov-Conte o r t h o g o n a l i z a t i o n procedure ( f o r example, see S e c t i o n 2 on t h e Method of Complementary F u n c t i o n s ) . W e include a s e p a r a t e treatment of t h e Method of Complementary Functions i n o r d e r t o i l l u s t r a t e t h e i m p o r t a n t c a s e of s e p a r a t e d boundary c o n d i t i o n s . For completeness t h e Method of A d j o i n t s w i l l be d i s c u s s e d even though i t h a s no s i m p l e r e l a t i o n s h i p t o (3.1) . 1. P a r a l l e l Shooting The boundary-value problem of i n t e r e s t i s -dudt A(t)u - f ( t ) = o t E [0,11 a) (3. 3 ) The simple s h o o t i n g method (Theorem 1.2) u s e s t h e s o l u t i o n of (n+l) i n i t i a l - v a l u e problems t o g e n e r a t e t h e s o l u t i o n of ( 3 . 3 a , b ) . In p a r a l l e l shooting, t h e i n t e r v a l [ o , l ] i s divided i n t o K n o n o v e r l a p p i n ~ s u b i n t e r v a l s , [Ty, Tv+l], and s e p a r a t e i n i t i a l - v a l u e problems a r e s o l v e d on each s u b i n t e r v a l . he s o l u t i o n of (3.3a,b) i s c o n s t r u c t e d by r e q u i r i n g c o n t i n u i t y a t Tv, v = 1,...,K-1, and s a t i s f a c t i o n of t h e boundary c o n d i t i o n s (3.3b). D i s c r e t e p a r a l l e l s h o o t i n g methods mimic t h i s procedure. we p l a c e a n e t ( t y ~ f :on ~ each s u b i n t e r v a l [ Tv-l, 'fV]. On t h i s n e t , a d i s c r e t e i n i t i a l - v a l u e problem i s s o l v e d w i t h i n i t i a l c o n d i t i o n s V V Zo, zo. These n+l i n i t i a l c o n d i t i o n s a r e a r b i t r a r y except t h a t V we r e q u i r e t h e n x n m a t r i x Z 0 s o l u t i o n can be r e p r e s e n t e d a s be n o n s i n g u l a r . This i n i t i a l - v a l u e o r , e q u i v a l e n t l y we w r i t e z V '1I zV 1 = jrv]. I I n p a r a l l e l s h o o t i n g , i t i s sometimes n e c e s s a r y t o s t a t e an i n i t i a l - v a l u e problem on some i n t e r v a l [T V' Tvcl] backward by f o r m a l l y g i v i n g n + l "boundary" c o n d i t i o n s a t T r a t h e r than v+l " i n i t i a l " conditions a t T v* However, Theoren 1.16 s t a t e s t h a t f o r a l l n e t s w i t h ho 5 H, H s u f f i c i e n t l y s m a l l , t h e s o l u t i o n t o t h i s backward problem i s e q u a l t o t h e s o l u t i o n of (3.4) w i t h zV and z V 0 a p p r o p r i a t e l y chosen. 0 Thus, i n o r d e r t o exaniine t h e most g e n e r a l p a r a l l e l s h o o t i n g methods, i t i s s u f f i c i e n t t o c o n s i d e r d i s c r e t e i n i t i a l - v a l u e problems of t h e form (3.4). V Defining Z , zv by (3.4), v J problem on t h e n e t { t i j i i 0 t h e s o l u t i o n of any i n i t i a l - v a l u e i s g i v e n by Thus, t h e c o n t i n u i t y requirements f o r p a r a l l e l s h o o t i n g a r e I n o r d e r t o meet t h e boundary c o n d i t i o n s , t h e f o l l o w i n g r e l a t i o n s h i p must b e s a t i s f i e d Combining ( 3 . 5 ) , (3.6), (3.,7), and (3.8), t h e p a r a l l e l shooting vPS, i s g i v e n by solution, Stability. For t h e purposes of t h i s d i s c u s s i o n , t h e discrete p a r a l l e l initial value problem is said to be stable if I n t h e p r o o f s t o f o l l o w , d e f i n e t h e nJv x n ( J +1) m a t r i x v M V via the relationship The i n i t i a l - v a l u e m a t r i x I h a s s o c i a t e d w i t h a p a r a l l e l shoot in^ procedure is For convenience, we w i l l w r i t e t h i s m a t r i x a s where t h e (1,o) block element of M' (K = i s t h e (SV+2, 3) block Sv+l) element of 1 h' Theorem 3.13. Let (3.3a,b) have a unique s o l u t i o n y ( t ) Let t h e d i f f e r e n c e scheme be c o n s i s t e n t with (3.3a). E c1 [ o , l ] . Then t h e following a r e equivalent: a) The d i s c r e t e i n i t i a l - v a l u e problem is stable. b) The d i s c r e t e p a r a l l e l initial value problem i s stable. Proof: b) => a) The m a t r i x lh can be f a c t o r e d i n t h e following - By hypothesis , I I(lh)v -1 I I 5 K, v = l , . , . , ~then 111 a) => -K 1 5 max IK , 11- b ) Without l o s s of g e n e r a l i t y , w e w i l l consider t h e c a s e i n which t h e i n t e r v a l [ o , l ] i s divided i n t o 3 s u b i n t e r v a l s , K = 3. The f i r s t s t e p i s t o show t h a t t h e m a t r i x i s nonsingular, 2 M lo' ' ' rh,defined below, By u s i n g t h e F a c t o r i z a t i o n Lemma (1.12), t h e Reducing Lemma (1.13), and t h e convergence of t h e i n i t i a l - v a l u e problem, i t can b e shown that where K 1 i s independent of ho. 1 2 3 L e t V , V ,V b e nJ1, n(J +I), nJ vectors respectively. 2 3- Consider t h e n ( J + l ) e q u a t i o n s where g2 i s any n(J2+1)-vector. g2 t h e r e i s a V 2 S i n c e 7 i s n o n s i n g u l a r , f o r each h such t h a t R e c a l l i n g (3.14), we f i n d t h a t Thus, ..., T h i s argument i s v a l i d f o r any v = 1, v -1 I I (lh) I I ' K1 V = K, 1, t h u s we have a . a , K . Before proceeding w i t h t h e theorem which w i l l t i e a l l t h r e e of t h e s e d i s c r e t e methods t o g e t h e r , i t i s convenient t o p r e s e n t two lemmas. Lemma 3.17. - L e t C1, C 2 b e nonsingular m x m matrices. L e t C1 and C2 be r e l a t e d by where L, N a r e m x q, q x m m a t r i c e s r e s p e c t i v e l y . L have rank q I m. Proof. Further, l e t Then t h e q x q m a t r i x [I + NL] i s n o n s i n g u l a r . By h y p o t h e s i s C1 i s n o n s i n g u l a r , t h e r e f o r e f o r each b E E?, t h e r e e x i s t s a unique m-vector, x, such t h a t R e c a l l i n g (3.18), t h i s e q u a t i o n can b e r e w r i t t e n a s o r equivalently a s The Reducing Lemma (1.13) t h e form where s t a t e s t h a t (3.19) o n l y h a s s o l u t i o n s of S i n c e a unique x, s o l u t i o n of (3.19), matrix ( I e x i s t s f o r every b , t h e + hZ) must b e n o n s i n g u l a r . Lemma 3.20. Let t h e ~ ( J + K )x ~ ( J + K )expanded boundary-value m a t r i x 8' b e d e f i n e d by h L e t t h e boundary-value mat r i x Bh b e n o n s i n g u l a r . Then nonsingular. Proof. V E Suppose t h a t B~ i s s i n g u l a r , t h e n t h e r e i s a ~ ( J + K ) - v e c t o r h such t h a t However, i f we d e l e t e t h e K-1elements v from V E i' ...,K-1, i = SV+v, v = 1, E and c o l l a p s e VE t o form a n(J+l)-vector, V, then By h y p o t h e s i s Rh i s n o n s i n g u l a r , t h u s each element vE d e l e t e d from V i E Note t h a t f o r t o form V and vE a p p e a r s a s a n element of V. i-1 and 1 [ v [1 = O . Thus, E by c o n t r a d i c t i o n , w e have 8 n o n s i n g u l a r . h m Now, w e s t a t e and prove t h e main theorem of t h i s Chapter, which w i l l i n c l u d e Theorem 1.16. Theorem 3.22. Let (3.3a,b) have a unique s o l u t i o n y ( t ) E c1[ o , l ] Let t h e d i f f e r e n c e scheme employed b e c o n s i s t e n t w i t h (3.3a). Denote by B V , IV, and PS the following methods: , BV. D i s c r e t e boundary-value procedure, (3.1) IV. D i s c r e t e i n i t i a l - v a l u e procedure, (3.2a,b,c), PS. Discrete p a r a l l e l shooting procedure, (3.9a,b,c). L e t one of the d i s c r e t e problems (3. I ) , ( 3 . 2 a ) , and ( 3 . 9 a ) be stable. Then, f o r a l l n e t s w i t h ho 5 11, H s u f f i c i e n t l y s m a l l , i) Each method uniquely d e f i n e s a n n ( J + l ) v e c t o r as an approximation t o y ( t ) on t h a t n e t , and ii) This approximation i s t h e same f o r each method. . Proof: Note that stability of (3.1) is equivalent to (3.2a) or (3.9a). Re i). Theorem 1.16 proves that methods SV and IV uniquely define approximations to y(t). To show that method PS uniquely defines an approximation it is sufficient to show that the nK x nK matrix in (3.9b) is nonsingular. Without loss of generality, we consider the case where K = 3. Thus, we want to show that the 3n x 3n matrix P, defined by is nonsingular. From conclusion i), Bh is nonsingular, hence by E Lemma 3.20, 8 is nonsingular. Recall the n(J+l)-vector h E (3.la), and define the ~(J+K)-vector r by Consider the matrix equation E Define the expanded parallel shooting matrix Ih by r in R e c a l l i n g t h e p a r a l l e l s h o o t i n 3 method [3,9a,b,c), eq-uation (3.23) can be w r i t t e n a s where and - I n z e n e r a l , L i s a n(J+K)x nK 1natrix, and N i s a nK + 1 x n(J+K) m a t r i x . < i s a n ( n ~+ 1 ) - v e c t o r , R e c a l l i n g Lemma 3.17, w e n o t e that i s n o n s i n g u l a r , where When t h e proper s u b s t i t u t i o n s a r e made, i t i s c l e a r t h a t (I i s n o n s i n g u l a r i f and o n l y i f P i s n o n s i n g u l a r . R e ii). Theorem 1.16 i s s u f f i c i e n t t o p r o v e t h a t - NL) vBV = vIV* Recalling (3.23), t h e Reducing Lemma 0 . 1 3 ) s t a t e s that where t h e n-vectors, wi, must s a t i s f y Now, we have t h e r e l a t i o n s h i p Recalling (3.23) and t h e d e f i n i t i o n of B:, we n o t e t h a t and thus it is clear that vps = vBv = v IV , Unfortunately, Theorem (3.22) contains no information on the numerical stability of any parallel shooting technique. To say that a unique solution to the difference equations (3.9a,b,c) exists is not to say that it can be accurately approximated numerically. It is well-known that a simple shooting method may produce disastrous numerical solutions if the linearly independent solutions of (3.3a) grow at different rates. For example, consider the problem This problem has a unique solution for all f(t) 2 f3 E E , because the matrix [Bo + BIX(l) 1 = 1 1 + e6o E 1 60 l+e l+e c[o,l] and -..] is nonsingular. However, if we try to invert this matrix numerically, on a machine with less than 85-bit accuracy, the computer will consider (Bo + BIX(l)) singular, because 1 + eg60 will be represented as 1. Thus, no matter how close ZJ gets to X ( 1 ) , this matrix will be numerically singular due to round-off error. For this reason, methods which are equivalent in the sense of Theorem 3.22 may yield greatly different numerical solutions. 2. Method of Complementary Functions The Method of Complementary Functions (Goodman and Lance [ 3 ] ) is used to approximate the solution of (3.3a,b) when the boundary conditions are separated. That is, the boundary conditions (3.3b) can be written as where Bo is a p x n matrix of rank p, B1 is a q x n matrix of rank qy p + q = n, and Bo, B1 are a p-vector, q-vector respectively. In this section, we use the convention that a single matrix in parentheses, (B,), has p rows and a single matrix in brackets, [B 1, has q rows. 1 The original equations (3.3a9b) can now be rewritten as In what follows, we assume that the boundary-value problem (3.27ayb) has a unique solution y(t). If the boundary-value problem has a unique solution, then there exists a vector q 0 such that (Do)qo = (Bo) where qo is orthogonal to the null space of (E 0 ). Let the orthonormal s e t { ~ ~ } span z = ~t h e n u l l s p a c e of (Bo), t h e n (qi, ?'lo) = 0 1 s i S q . With t h e s e n o t i o n s i n hand, t h e s o l u t i o n of (3,27a,b) can b e c h a r a c t e r i z e d i n t h e f o l l o w i n g way. Solve t h e q homogeneous i n i t i a l - v a l u e problems and t h e inhomogeneous problem The unique s o l u t i o n of (3.27a9b) i s g i v e n by where t h e c o n s t a n t s cii, i = . . q , a r e determined from I n t h e i n t r o d u c t i o n t o t h i s c h a p t e r , i t was s t a t e d t h a t t h e Method of Complementary Functions i s a s p e c i a l c a s e of p a r a l l e l shooting. clear. T h i s i s t h e b e s t p o i n t a t which t o make t h i s r e l a t i o n s h i p P a r a l l e l s h o o t i n g n a t u r a l l y c o n t a i n s t h e c a s e of s i m p l e shooting. I n o r d e r t o c h a r a c t e r i z e y ( t ) by means of s i m p l e s h o o t i n g , s o l v e t h e n homogeneous i n i t i a l - v a l u e problems and t h e inhomogeneous problem -dx(t) - A ( t ) x ( t ) + f ( t ) dt where C i s a n o n s i n g u l a r n x n n a t r i x and c i s an n - v e c t o r . where t h e n-vector Now, i s determined by The r e s u l t of s p e c i a l i z i n g t h e s i m p l e s h o o t i n g p r o c e d u r e ( 3 , 3 2 a , b ) , (3.33a,b), (3.34), (3.35) by t a k i n g i s t h e Method of Complementary F u n c t i o n s . The d i s c r e t e Method of Complementary F u n c t i o n s ( h e r e a f t e r r e f e r e n c e t o t h e Method of Conzplementary F u n c t i o n s w i l l imply t h e d i s c r e t e p r o c e d u r e ) i s d e r i v e d by s o l v i n g e q u a t i o n s ( 3 . 2 8 a Y b ) , (3.29a,b), (3.30), and (3.31) numerically. Solve the q homogeneous discrete initial-value problems - c. 'li 0 0 L r and the inhomogeneous discrete initial-value problem For convenience, we write (3.32), (3.33) as where Z is a (J+l)n x q matrix with its ith column being z discrete approximation to y ( t ) is given by where i . The The i n i t i a l - v a l u e e q u a t i o n (3.34) may b e f o r m a l l y w r i t t e n a s (Go) I I I II [ o l I l o I I I I . I I I I l o , The f i r s t n(q+l) e q u a t i o n s of (3.37) s t a t e t h a t t h e columns of Zo a r e orthonormal and belong t o t h e n u l l - s p a c e of Boy and a l s o t h a t 0 Bo Zo = Po and zo0 i s o r t h o g o n a l t o t h e n u l l - s p a c e we want t o show t h a t V? c a l c u l a t e d v i a (3.37), of B 0 (3.35), . Now, and (3.36) i s a s o l u t i o n of t h e e q u a t i o n s The d i s c r e t e boundary-value problem i s R e c a l l i n g t h e F a c t o r i z a t i o n Ler~una (1.12), t h i s m a t r i x can b e w r i t t e n as where T is the n(~+1) x n(J+l) matrix in (3.37). h by Left-multiplying T~~~and recalling (3.37), this equation is shown to be The Reducing Lemma (1.13) states that where w and A are determined from When this is unraveled, we have Recall (3.35) and (3.36) and compare with (3.30) and (3.40) to show that 3. Method of A d i o i n t s The Method of A d j o i n t s completes t h e s t u d y of s o l u t i o n procedures f o r d i s c r e t e boundary-value problems. This method i n t r o d u c e s t h e a d j o i n t e q u a t i o n which can n o t be d e f i n e d v i a l i n e a r t r a n s f o r m a t i o n s of t h e o r i g i n a l problem. For t h i s r e a s o n t h e e q u i v a l e n c e e x h i b i t e d between boundary-value s h o o t i n g procedures i s l a c k i n g h e r e . t e c h n i q u e s and p a r a l l e l - A good e x p l a n a t i o n of t h e Method of A d j o i n t s i s found i n Goodman and Lance [ 3 ] . The problem c o n s i d e r e d i n [ 3 ] i s a v e r y s p e c i a l one: -dudt A(t)u = f ( t ) a) The a d j o i n t d i f f e r e n t i a l e q u a t i o n i s From (3.41), (3.42), t h e following r e l a t i o n s h i p i s derived I n t e g r a t i n g t h i s expression, we o b t a i n By s o l v i n g t h e a d j o i n t problem n-p times w i t h i n i t i a l c o n d i t i o n s ...,n , t h e e q u a t i o n (3.44) r e s u l t s i n a s e t of n-p ( o r q ) e q u a t i o n s f o r t h e q u a n t i t i e s ui(o), i = p + l , ...,n. -V A u (0) = e v' v = p+l, Now, the initial conditions u(o) are explicitly known and the solution of a single initial-value problem determines y(t). If we use the sane discretization for (3.41a,bYc) and (3.42), (3.44), it is not in general true that the solution of the discrete Method of Adjoints, ?, is equal to vBV. However, if we employ the centered-Euler scheme this equivalence does hold. The discrete analogue to (3.44) is naturally -* MA v~ "J -* it 0 0 i-I ) i=l f(ti_ l) = 0- (3.45) 2 From the difference equations and the following identity can be derived ...,n, this identity yields Upon summing from i = 1, - - (vi + vi-1) ik BV BV YA also satisfies (3.66) and vi = vi , i = 0,. . . , J . Thus, V Remember t h a t t h i s e q u i v a l e n c e h a s been shown o n l y i f t h e c e n t e r e d E u l e r scheme i s used t o d i s c r e t i z e t h e Method of A d j o i n t s . Identities s i m i l a r t o (3.47) can be d e r i v e d f o r any number o f o t h e r d i f f e r e n c e schemes, b u t t h e summation which y i e l d e d (3.48) does n o t i n g e n e r a l show e q u i v a l e n c e . Chapter IV Nonlinear Two-Point Boundary-Value Problems T h i s c h a p t e r d e a l s w i t h numerical approximations t o t h e s o l u t i o n of n o n l i n e a r two-point boundary-value problems of t h e form where u, f , b a r e n-vectors. i s o l a t e d s o l u t i o n s of ( 4 . l a , b ) . A f t e r Keller [GI, we a r e i n t e r e s t e d i n The s o l u t i o n of ( 4 , l a , b ) i s i s o l a t e d i f t h e l i n e a r i z e d boundary-value problem has only t h e t r i v i a l solution. Note t h a t w e have used t h e n o t a t i o n When c o n s i d e r i n g t h e d i s c r e t e problern, w e w i l l a l s o examine ~ e w t o n ' smethod a s a means of c a l c u l a t i n g an a p p r o x i n a t i o n t o y ( t ) . The n ( J + l ) - v e c t o r Y i s d e f i n e d t o be where y ( t ) i s a n i s o l a t e d s o l u t i o n of ( 4 . l a , b ) . S,,[g(t)l = Ix We a l s o d e f i n e I x E E", 1 lx - g ( t ) l 1 5 p } , 1. F i n i t e D i f f e r e n c e Schemes Each of t h e d i f f e r e n c e schemes d i s c u s s e d i n Chapter 2 can be modified t o make i t a p p l i c a b l e t o n o n l i n e a r d i f f e r e n t i a l e q u a t i o n s of t h e form ( 4 . l a ) . Thus, t h e centered-Euler scheme (2.2) becomes f o r a general nonlinear equation. I f equation (4.la) is linear i n u, t h a t i s , t h e n t h e d i f f e r e n c e e q u a t i o n (4.3) reduces t o (2.2). The n o n l i n e a r analogue of t h e i n t e g r a l e q u a t i o n (2.1) is clearly A s i n Chapter 2, q u a d r a t u r e fornlulae and i n t e g r a t i o n by p a r t s can be used t o g e n e r a t e c o n s i s t e n t numerical schemes. I n t h i s manner t h e g e n e r a l form of t h e t r a p e z o i d a l r u l e i s shown t o be Again n o t e t h a t i f f (u, t ) i s l i n e a r i n u, t h e n (4.5) r e d u c e s t o which i s t h e t r a p e z o i d a l r u l e a p p l i e d t o a l i n e a r problem. (2.10), The Gap schemes can be d e r i v e d i n t h e same manner a s i n Chapter 2. R e c a l l i n g L e m a 2.11 and Lemma 2.14, i t can be shown that where t h e c o n s t a n t s n a r e g i v e n by Table 2.17. Gap 4 , we need t o e v a l u a t e - I n order t o derive . Using t h e d i f f e r e n t i a l e q u a t i o n (4*la), Thus, we d e f i n e t h e n o n l i n e a r Gap4 scheme by + hi [ f 1(Vi' t .1) f ( v i 9 t i ) + f2(vi,ti) - fl(vi-l,ti-l) T h i s i s a f o u r t h - o r d e r a c c u r a t e scheme i f y ( t ) E f(vi-l'ti-l 5 c [o,l], Coupled w i t h any of t h e s e d i f f e r e n c e schemes, t h e boundary 1 conditions (4,lb) are approximated by b(v ,v ) = o. O J Employing the trapezoidal rule, an approximation to the solution of (4.la,b) can be determined from the following set of n(J+l) equations: More concisely we write where N v = o may represent any difference approximation of the h i differential equation (4.la) . The truncation error associated with such a difference scheme is defined to be where y(t) is a solution of (4.1a9b). The truncation error associated with the boundary condition is ~,[~(t)l and is always zero. n(V) b(y(o), ~(1)) Denote the n(J+l)-equations (4.9) = 0 where the n(J+l)-vector (4.8a,b) by V is 2. E x i s t e n c e of S o l u t i o n s W e w i l l now c o n s i d e r t h e q u e s t i o n of e x i s t e n c e of s o l u t i o n s A s a prelude, t o systems of n o n l i n e a r e q u a t i o n s of t h e form ( 4 . 8 a , b ) . l e t us examine t h e r e s u l t of a p p l y i n g t h e s e n o n l i n e a r d i f f e r e n c e schemes t o t h e l i n e a r i z e d problem. Lermna 4.10, L e t t h e boundary-value problem ( 4 . l a , b ) solution, y(t). have a n i s o l a t e d Let be t h e n ( J + l ) e q u a t i o n s which r e s u l t from t h e a p p l i c a t i o n of t h e d i f f e r e n c e scheme (4.8a,b) t o t h e l i n e a r i z e d problem: where g , f3 a r e n-vectors, and g ( t ) E 1 c [0,1]. Let t h e d i f f e r e n c e scheme be c o n s i s t e n t and s t a b l e when a p p l i e d t o a l i n e a r two-point boundary-value problem w i t h a unique s o l u t i o n , z ( t ) [ y ( t ) ] x [ o , l l ] and b(x,y) E c1 [S [ Y ( o ) ] 2 c [o,l]. Let Sp[~(1)]Is P P Then f o r a l l n e t s w i t h ho 2 H, H s u f f i c i e n t l y s m a l l , L i s n o n s i n g u l a r j=(y,t) r and E where K i s independent of ho. Proof. By h y p o t h e s i s , e q u a t i o n (4,12a,b) h a s a unique s o l u t i o n z (t) c2 [o , I ] . 0 E t o (4.12a,b), Also, t h e d i f f e r e n c e scheme i s s t a b l e when a p p l i e d thus there e x i s t constants Ilvll f o r a l l nets with h Theorem 4.13. KO llrll 0 5 H. L e t (4.la'b) * This i s equivalent t o have a n i s o l a t e d s o l u t i o n y ( t ) Let t h e h y p o t h e s i s of Lemma 4.10 hold. (4.8a,b) be c o n s i s t e n t w i t h ( 4 . l a , b ) , q (V) be continuous f o r V E S v % t o , H > o such t h a t Po [y(t)]. E cL[o,l]. Let t h e d i f f e r e n c e scheme and l e t t h e J a c o b i a n m a t r i x F u r t h e r , l e t "(v) be such t h a t and where Then f o r each p , o < p 5 p0 ' po s u f f i c i e n t l y s m a l l and. for a 1 1 n e t s Il(p) s u f f i c i e n t l y s m a l l , t h e n ( ~ + l ) e q u a t i o n s w i t h ho 5 EL@), rl(v) = 0 have a unique s o l u t i o n W E S [y (t ) ] Proof: . W e define t h e vector function I f t h i s h a s a s o l u t i o n , $(W) = W, t h e n n(W) = o. The C o n t r a c t i n g Mapping Theorem w i l l be used t o show t h a t a unique s o l u t i o n e x i s t s . That i s , we need t o show t h a t f o r some h E ( o , l ) i) 1 IY - $(Y) I I 5 (1 - ii) II1/'(v1) - $ ( v 2 ) I I 5 x I I v ~ - v21[, where Vl, E V2 Sp[y(t)l By t h e h y p o t h e s i s t h e numerical scheme i s c o n s i s t e n t , hence Re i i ) Now, w e employ the mean v a l u e theorem and immediately By h y p o t h e s i s q (V) i s c o n t i n u o u s , t h u s v I IJ, ( 1 -q ( 1 I I 5 KoIK1p + I IL - rlV(Y) I I } I Iv1 - v21 I By h y p o t h e s i s , p r o p e r t y ( 4 . 1 4 ) , Thus, l e t po and H(po) be s u f f i c i e n t l y s m a l l such t h a t and c o n d i t i o n i i ) i s s a t i s f i e d . Let H(P) be s u f f i c i e n t l y s n a l l such that Then J,(V) i s a c o n t r a c t i n g map f o r V E S [ y ( t ) ] and t h e r e e x i s t s P a unique s o l u t i o n W E S [ ~ ( t ) ] . P m C o r o l l a r y 4.17 (Convergence). L e t t h e h y p o t h e s i s of Theorem 4.13 w E sp[y(t)], hold. Then, i f q ( W ) Proof. It i s s u f f i c i e n t t o prove t h a t f o r any E > o , t h e r e e x i s t s = o, an H such t h a t f o r any n e t w i t h ho 5 H. However, t h i s i s t h e c o n c l u s i o n of Theorem 4.13 w i t h E r e p l a c i n g p. rn C o r o l l a r y 4.18. Let t h e h y p o t h e s i s of Theorem 4.13 hold. Let t h e d i f f e r e n c e scheme (4.8b) be pth o r d e r a c c u r a t e , Then, f o r a l l n e t s w i t h ho 5 H, li s u f f i c i e n t l y s m a l l , Proof. From t h e proof of Theorem 4.13, w e know t h a t W E S [ y , t ] P i f p 5 po, ho 5 H(po), and where X E (0,1) and f i x e d . By h y p o t h e s i s and t h e r e f o r e Theorem 4.13 s t a t e s t h a t W E S [ y ( t ) ] i s a s o l u t i o n i f P t h u s i t must be t r u e f o r Hence, Idhat d i f f e r e n c e schemes s a t i s f y c o n d i t i o n (4.14)? It would appear t o be f o r t u i t o u s indeed i f any schene had t h e property t h a t however t h i s i s t h e c a s e f o r many common d i f f e r e n c e schemes. examples, t h e centered-Euler (4.3) and t h e Gap 4 As (4.5) schemes will be examined. Gap4. Application of t h e Gap4 scheme t o the l i n e a r i z e d problem P (4.12a, b) y i e l d s t h e f o l l o w i n g n(J+l) e q u a t i o n s where t h e right-hand s i d e of (4.19b) i s w r i t t e n simply a s si because i t i s n o t of p a r t i c u l a r i n t e r e s t h e r e . From (4.19a,b), completely c h a r a c t e r i z e t h e elements ( L . . ) 1J w e can of L. Let t h e ( i , jlth n X n b l o c k element of nV(Y) be Nij Similarly, N and Noo = Loo. 0j = Lo , j = l,...,J. b e shown t h a t Similarly, N.. 11 = Lii, ...,J. i = 1, . Then From 4 . 6 , i t can Thus, t h e Gap4 scheme s a t i s f i e s c o n d i t i o n (4.14) and i n f a c t Centered-Euler. The centered-Euler scheme s a t i s f i e s e q u a t i o n (4.14), however L i s n o t e q u a l t o nV(y). Employing t h e centered-Euler ( 4 3 ) , t h e J a c o b i a n m a t r i x 17 (Y) h a s t h e f o l l o w i n g nonzero n x n v b l o c k components: The m a t r i x L h a s t h e s e n x n block elements Thus, we have Clearly f o r f ( u , t ) E 1 c [Sp[y(t)l [ o , l ] l and y ( t ) E 1 c [o,l], scheme, p r o p e r t y (4.14) holds. - Note t h a t t h e q u a n t i t y need n o t approach z e r o l i k e hP i n o r d e r f o r a pth o r d e r a c c u r a t e 0 scheme t o s a t i s f y 3. S o l u t i o n of q(V) = o by Newton's Method I n t h e p r e v i o u s s e c t i o n i t was shown t h a t under c e r t a i n c o n d i t i o n s a s o l u t i o n e x i s t s t o t h e n(J+1) n o n l i n e a r e q u a t i o n s W e would l i k e t o examine a s p e c i f i c procedure f o r s o l v i n g t h e s e e q u a t i o n s , namely ~ e w t o n ' smethod. That i s , w e want t o show t h a t each i t e r a t e d e f i n e d by V e x i s t s u n i q u e l y and t h a t t h e sequence {V ) converges when V j u d i c i o u s l y chosen. 0 is The proof of t h e f o l l o w i n g theorem i s i d e n t i c a l t o t h a t given i n K e l l e r [ 6 ] , with t h e a p p r o p r i a t e generalizations. Theorem 4.23 Let t h e h y p o t h e s i s of Theorem 4.13 hold and l e t W be t h e n ( J + l ) - v e c t o r of t h a t theorem where L e t V O b e such t h a t VO small. E S [y(t) P Then f o r a l l n e t s w i t h h 0 f o r a l l p, o < p 5 po, p 0 and [ IvO - w [ I i s s u f f i c i e n t l y 5 H, H s u f f i c i e n t l y s m a l l , and s u f f i c i e n t l y small, t h e Newton i t e r a t i o n (4.22) u n i q u e l y d e f i n e s a sequence {V v 1 and t h i s sequence converges q u a d r a t i c a l l y t o W, Proof: A s noted b e f o r e , see Keller [ b ] . Chapters 1 and 3 examined t h e e q u i v a l e n c e between s h o o t i n g methods and i m p l i c i t t e c h n i q u e s f o r l i n e a r boundary-value problems. I n t h e n e x t s e c t i o n , w e want t o d e f i n e a n o n l i n e a r s h o o t i n g procedure and s t u d y i t s r e l a t i o n s h i p t o t h e method d e s c r i b e d h e r e . 4. Equivalence of Shooting and I m p l i c i t Schemes f o r Nonlinear Problems The n o n l i n e a r s h o o t i n g scheme h a s t h e same u n d e r l y i n g i d e a a s i n Chapter 1. The i n i t i a l - v a l u e problem i s d i s c r e t i z e d , s o l v e d , and t h e n t h e i n i t i a l c o n d i t i o n x(o) = c i s wiggled u n t i l t h e d i s c r e t e boundary c o n d i t i o n b ( c , u J ( c ) ) = o i s satisfied. For example, t h e f o l l o w i n g system of e q u a t i o n s r e s u l t i f we u s e ~ u l e r ' smethod t o approximate (4.24a): Note t h a t t h e s e e q u a t i o n s can be solved e x p l i c i t l y , however a n i t e r a t i o n scheme must be employed t o s a t i s f y t h e n o n l i n e a r boundary condition The a p p l i c a t i o n of Newton's Method t o s o l v e (4.26) i s s t r a i g h t forward 'du However, t h e n x n m a t r i x ac must be e v a l u a t e d a t each s t e p . m a t r i x i s determined by s o l v i n g t h e v a r i a t i o n a l e q u a t i o n s A d i s c r e t e s h o o t i n g procedure r e s u l t s from combining ( 4 , 2 5 a , b ) , (4.28a,b), and (4.27): This An a l t e r n a t i v e t o t h i s procedure i s ~ e w t o n ' smethod a s described i n t h e previous section. That i s , Newton's method i s employed t o s o l v e t h e system of e q u a t i o n s This r e s u l t s i n t h e following i t e r a t i o n scheme It has a l r e a d y been shown t h a t {V v 1 e x i s t s and converges. It i s c l e a r t h a t i n each c a s e t h e method d e f i n e s a s o l u t i o n of t h e same s e t of e q u a t i o n s : The i m p l i c L t method (4,31a,b) a p p l i e s Newton's method t o e a c h of t h e n ( J + l ) e q u a t i o n s of (4,30a, b) . The s h o o t i n g method (4.29a, b , c , d , e ) f o r m a l l y s e p a r a t e s t h e f i r s t n e q u a t i o n s from t h e remaining n J , I n l i e u of t h i s set of e q u a t i o n s , t h e s h o o t i n g method s u b s t i t u t e s the equations which can be s o l v e d e x a c t l y when a n e x p l i c i t d i f f e r e n c e scheme i s used t o approximate (4.23a). Then t h e i n i t i a l v a l u e c i s v a r i e d u n t i l a n i n i t i a l v a l u e i s found such t h a t The p r i n c i p a l d i f f e r e n c e i s t h a t t h e s h o o t i n g procedure s a t i s f i e s t h e l a s t n J e q u a t i o n s e x a c t l y and t h e i m p l i c i t Newton procedure does n o t . where i n g e n e r a l Theorem 4.34 That i s , I IV; I I # 0. L e t t h e h y p o t h e s i s of Theorem 4.16 hold. s h o o t i n g Newton i t e r a t e s be {U V Let t h e 1 , Uv E S P [ y ( t ) 1, and l e t t h e i m p l i c i t v Newton i t e r a t e s b e {V 3, Proof: -v , [[;'I Let ?(vv) = r [ > o. Then uV f vV. Define r V = q(uV) - q(vV) A s i n Theorem 4,13, t h e mean v a l u e theorem i s employed i n t h e following way By hypothesis, q v i s continuous on S P [ y ( t ) ] , hence llrVll 5 Q lluV - vVII A s i n t h e proof of Theorem 4 . 1 3 , i t can b e shown t h a t the m a t r i x 1 V lo v (sU + (1-s) vV)ds i s nonsingular, thus Q > o. Immediately from t h e above, Thus, 1 luV - vV[I > o -> uV # vV. This theorem s t a t e s t h a t t h e i t e r a t i o n sequences defined by i m p l i c t Newton and shooting Newton a r e not t h e same, even though t h e y a r e used t o s o l v e t h e same s e t of nonlinear equations. Chapter V Block T r i d i a g o n a l M a t r i c e s T h i s c h a p t e r develops a computational method f o r s o l v i n g t h e systems of l i n e a r e q u a t i o n s t h a t have a r i s e n i n p r e v i o u s chapters. That i s , we w i l l develop a t e c h n i q u e f o r s o l v i n g t h e n (J+l) l i n e a r equations o r , when t h e boundary-value problem i s n o n l i n e a r , t h e e q u a t i o n s %(vv) (vW1 - vV) + T,(vV) g e n e r a t e d by ~ e w t o n ' s method, = 0 More s p e c i f i c a l l y we a r e i n t e r e s t e d i n m a t r i c e s which r e s u l t from t h e approximation of two-point boundary-value problems w i t h s e p a r a t e d boundary c o n d i t i o n s . That i s , problems which may be w r i t t e n a s u' ( t ) = f (u,t) t E [0,1] a) bo(u(o)) = 0 b) bl(u(l)) = 0 c) where t h e boundary c o n d i t i o n s (5. l b ) and (5. l c ) r e p r e s e n t p and q conditions respectively, p + q = n. To i l l u s t r a t e t h e form t h a t t h e s e m a t r i c e s w i l l t a k e , consider the approximation of by ~uler's method. The resulting finite difference equations ( J = 2 ) are where we have used the convention that a single matrix in parentheses, (B0), has p rows and a single matrix in brackets, [B1], has q rows. The matrix of (5.3) can be rewritten in block tridiagonal form where each Ai, Bi, Ci is an n x n matrix. Block matrices of this type have a very special zero structure: the first p rows of Bi and the last q rows of Ci are zero rows. The solution procedure of interest is a block LU decomposition: I d e a l l y , w e would l i k ; t o show t h a t t h i s decomposition e x i s t s i f t h e matrix 8 is nonsingular. h that clear. Unfortunately t h e r e s u l t i s not For example, i f A ( t ) i n (5.2a) i s a 4 x 4 diagonal m a t r i x and (Bo) = (1 o o o ) , t h e n and A 0 is clearly singular. No m a t t e r how s m a l l h nd x o Ao-[!x!!] row, t h e n o - A 0 becomes, t h i s However, i f t h e 5th row of Rh i s matrix w i l l remain s i n g u l a r . interchanged w i t h t h e 2 0 o , i s n o n s i n g u l a r , and t h e LU f a c t o r i z a t i o n can proceed. Thus, a row s w i t c h i n g s t r a t e g y w i l l be i n c l u d e d i n t h e decomposition procedure. 1. Block LU Decompostion The block t r i d i a g o n a l m a t r i x M i s c o n s i d e r e d i n the f o l l o w i n g generality ; where each Ai, Bi, C. i s a n n x n matrix. 1 W e also require that the f i r s t p rows of each Bi b e z e r o and t h e l a s t q rows of Ci be zero. Such a m a t r i x M w i l l b e r e f e r r e d t o a s a b l o c k t r i d i a g o n a l matrix w i t h p / q z e r o s t r u c t u r e . The s t r a i g h t f o r w a r d b l o c k decomposition i s where t h e m a t r i c e s Li, Di are d e f i n e d by Once t h e LU decomposition i s completed, t h e s t a n d a r d procedure i s adopted t o s o l v e MV = r: LW= r , w = W., The n(J+l) v e c t o r s W and V a r e defined by t h e r e c u r s i o n r e l a t i o n s : and I f i t should occur t h a t some Di i s s i n g u l a r , then t h e s e r e c u r s i o n formulae w i l l have t o be a l t e r e d t o allow f o r t h e a p p r o p r i a t e row interchanges. M denote n(m+l) x n(m+l) block t r i d i a g o n a l Let M', Lemma 5.10. m a t r i c e s w i t h p / q zero s t r u c t u r e . MO, 0 0 M by Ai, Bi, nonsingular. of rows (p A 0 0 Ci and Ai, Bi, 0 Then e i t h e r A 0 Denote t h e block elements of Ci r e s p e c t i v e l y . Let MO be i s nonsingular o r by an i n t e r c h a n g e + 1,..., p + n) of MO a m a t r i x M can be formed such t h a t i s nonsingular. Proof. m = o. - 0 Immediate because M 0 = Ao. m > o. Let (S be t h e p n m a t r i x composed of t h e f i r s t P ) 0 p rows of Ao. Let S be t h e n x n m a t r i x whose f i r s t q rows a r e n t h e l a s t q rows of A' 0 0 of C1. x and whose l a s t p rows a r e t h e f i r s t p rows With t h e s e d e f i n i t i o n s i n hand, t h e m a t r i x MO can be written a s where Q i s a n (nm-tq) x nm m a t r i x . By h y p o t h e s i s MO i s n o n s i n g u l a r , t h u s t h e rows of (S ) a r e l i n e a r l y independent and t h e (n+p) x n P m a t r i x S d e f i n e d by These two c o n d i t i o n s imply t h a t t h e r e a r e q rows of h a s r a n k n. Sn which, t o g e t h e r w i t h t h e rows of (S ), form a set of n l i n e a r l y P independent row v e c t o r s . 0 rows of M Thus, by i n t e r c h a n g i n g t h e a p p r o p r i a t e we can form a m a t r i x M such t h a t A S i m i l a r rows. 0 i s nonsingular. Define rows of a block t r i d i a g o n a l m a t r i x w i t h p/q z e r o s t r u c t u r e t o be s i m i l a r i f by i n t e r c h a n g i n g t h e s e rows t h e m a t r i x remains b l o c k t r i d i a g o n a l w i t h p / q zero s t r u c t u r e . Theorem 5.11, L e t MO be a n n ( J + l ) x n ( J + l ) b l o c k t r i d i a g o n a l matrix with p/q zero s t r u c t u r e . 0 Let M be n o n s i n g u l a r . 0 by a n i n t e r c h a n g e s t r a t e g y among s i m i l a r rows of M Then a matrix M can be formed such t h a t M h a s a b l o c k LU decomposition of t h e form where each Dip Proof. o 5 i 5 J, i s nonsingular. R e c a l l i n g Lemma 5.10, among rows (p + l , . . . , p + n ) of M0 t h a t forms a n MI such t h a t where D1 i s n o n s i n g u l a r and M' 0 with p/q zero s t r u c t u r e . t h e following p r o p e r t i e s : i) where there i s an interchange strategy Mk = F~ G~ i s a block t r i d i a g o n a l m a t r i x k Assume t h a t t h e r e e x i s t s a m a t r i x M w i t h and - 1, i s n o n s i n g u l a r , and ii) i+l ,o 5 i 5 k Each Di iii) Mk i s a b l o c k t r i d i a g o n a l m a t r i x w i t h p / q z e r o s t r u c t u r e formed by a n i n t e r c h a n g e s t r a t e g y of s i m i l a r rows of M0 . 1 Since M has p r o p e r t i e s i ) from M thus k , i i ) , and i i i ) i t i s s u f f i c i e n t t o show t h a t w e can proceed t o Mk+'. By h y p o t h e s i s M0 i s n o n s i n g u l a r , Mk must be n o n s i n g u l a r . Define the n(J+l-k) x n(J+l-k) m a t r i x H by n o t i n g t h a t H i s b l o c k t r i d i a g o n a l w i t h p/q zero s t r u c t u r e . matrix Mk i s n o n s i n g u l a r i f and o n l y i f H i s n o n s i n g u l a r . Lemma 5.10, t h e r e i s a n i n t e r c h a n g e s t r a t e g y among rows (p of H which produces a n such t h a t The Recalling + 1,. . . , p + n ) where k i s nonsingular . Note t h a t an interchange of any rows (p+l, . . . , p + n ) of H i s equivalent t o an interchange of rows (nk + p + 1,...,nk + p + n) of M'. k+l k+l k By d e f i n i n g Lk+l, Dk 9 Bk, and k+l i n t h i s manner, we have found M = Fk+l G ~ + ' which has p r o p e r t i e s i ) , i i ) , and i i i ) . m Corollary 5.12. Let MO be an n ( J + l ) x n ( J + l ) block t r i d i a g o n a l matrix with p/q zero s t r u c t u r e . Let P be a permutation m a t r i x of t h e form i s t h e p x p i d e n t i t y , I i s t h e q x q i d e n t i t y , and each q P. i s an n x n permutation matrix. Then t h e following a r e e q u i v a l e n t where I P 1 a) 0 M i.s n o n s i n g u l a r . b) There e x i s t s a p e r m u t a t i o n m a t r i x P such t h a t where each Di, Proof: o 5 i 5 J, i s n o n s i n g u l a r . a ) => b ) . The proof i s e x a c t l y t h a t of Theorem 5.11. b) => a ) . Immediate from h y p o t h e s i s . Using t h i s r e s u l t , a b l o c k t r i d i a g o n a l f a c t o r i z a t i o n can be designed t o s o l v e t h e e q u a t i o n s where M i s a b l o c k t r i d i a g o n a l m a t r i x w i t h p / q z e r o s t r u c t u r e . R e c a l l i n g C o r o l l a r y (5.12), t h e e q u a t i o n s can be w r i t t e n a s which a l l o w s s o l u t i o n v i a (5.6) , (5.8) , and (5.9) . The o n l y problem remaining i s t h e p r a c t i c a l one of determining each of t h e permutation m a t r i c e s P i comprising P. Each Pi can be found a s t h e LU decomposition (5.6) proceeds. Theorem 5.13. structure. Let M be a block t r i d i a g o n a l m a t r i x w i t h p / q z e r o L e t M be n o n s i n g u l a r . L e t r , V , W be n ( J + l ) v e c t o r s d e f i n e d by Then t h e following procedure can be used t o f i n d t h e s o l u t i o n of MV = r : I, Forward s u b s t i t u t i o n -- Ao, wo = r . A. S t a r t i n g values: B, I t e r a t i o n procedure (k = o , l , 1. Do 0 ...,J-1) I f t h e c u r r e n t D i s s i n g u l a r , i n t e r c h a n g e rows k (kn + p + l , . , , , k n + p + n) of M i n o r d e r t o form a n o n s i n g u l a r D and form new w k k' rk+l by i n t e r c h a n g i n g corresponding elements. 2, Compute: - -1 4c+l - Ck+l Dk Dk+l = %+l- Lk+l Bk+l w = r k+l k+l - Lk+l wk 11. Back s u b s t i t u t i o n A. Proof: Compute v k (k = J , J - ~ , . . . , o ) via D k v - k - Wk - Bk+l Vk+l The proof i s t h e same a s t h e proof of Theorem 5.11, Thus f a r o n l y one method of decomposing t h e m a t r i x M has been considered: t o block m a t r i c e s , t h e n a t u r a l e x t e n s i o n of t h e u s u a l LU decomposition Consider a n LU decomposition where U r a t h e r than L has a u n i t diagonal: The g e n e r a l r e c u r s i o n r e l a t i o n s f o r Di, Ui a r e Allowing f o r row i n t e r c h a n g e s , t h i s . d e c o m p o s i t i o n can always be performed, however i t h a s a p r a c t i c a l d i s a d v a n t a g e . I n t h e LU decomposition ( 5 . 6 ) , L . h a s t h e same z e r o s t r u c t u r e a s C 1 t h e l a s t q rows a r e zero. i' namely Thus, t h e m a t r i c e s L, U c a n b e s t o r e d i n t h e same l o c a t i o n s a s t h e o r i g i n a l m a t r i x M. However, from (5.14) i t i s c l e a r t h a t Ui i s i n g e n e r a l a f u l l n x n m a t r i x , and t h i s decomposition r e q u i r e s pnJ more s t q r a g e l o c a t i o n s t h a n M. 2. A S p l i t decomposition Another p o s s i b i l i t y i s t h e f a c t o r i z a t i o n M = UL where U h a s a u n i t d i a g o n a l . That i s , i n t h e c a s e J = 2, we have Let PI0 be an n<J+1) x ncJ+I) b l o c k t r i d i a g o n a l C o r o l l a r y 5.15. m a t r i x w i t h p/q z e r o s t r u c t u r e , Let P b e a permutation m a t r i x of t h e same form a s C o r o l l a r y 5.12, Then t h e f o l l o w i n g a r e equivalent: 0 M b) There e x i s t s a permutation m a t r i x P such t h a t where each Di, Proof, where q i s nonsingular . a) o S i 5 J, i s nonsingular. b) ==> a) Immediate from h y p o t h e s i s . a ) => b ) Let Q b e t h e n ( J + l ) x n ( J + l ) permutation m a t r i x i , n (J+l)-i+l = 1. Then go = QM'Q m a t r i x w i t h q/p zero s t r u c t u r e . e x i s t s a permutation m a t r i x i s a nonsingular t r i d i a g o n a l R e c a l l i n g C o r o l l a r y 5.11, there P such t h a t where % i s a block t r i d i a g o n a l m a t r i x w i t h q/p z e r o s t r u c t u r e . Note t h a t QQ = I , t h e n w e have where P = QFQ, U = QE, and L = 60. Rather t h a n p u r s u i n g this UL decomposition d i r e c t l y , we w i l l examine a h y b r i d o r s p l i t procedure. The m a t r i x M might be f a c t o r e d i n t h e f o l l o w i n g way: This decomposition i s a hybrid of t h e LU and t h e UL decompositions discussed previously. L e t MO b e a n n ( J + l ) x n ( J + l ) b l o c k t r i d i a g o n a l m a t r i x Theorem 5.19. with p/q zero s t r u c t u r e , Let ~ ( s be ) a permutation m a t r i x of t h e same form a s C o r o l l a r y 5.12. 0 Then t h e f o l l o w i n g a r e e q u i v a l e n t : a) M i s nonsingular., b) For each i n t e g e r s , o I s I J , t h e r e e x i s t s a permutation m a t r i x P ( s ) such t h a t 14 = P(s)M of t h e form: where 0 has a s p l i t f a c t o r i z a t i o n and each Di, Proof: o 2 i 5 J, is nonsingular. The proof i s sketched here. The o n l y element of F o r G t h a t i s n o t guaranteed e x i s t e n c e by C o r o l l a r y (5.12) ... , s ; Di, i = o , ...,s-1) o r Corollary(5.15) (Ui, i =s+l,.. . , J ; Di, i = s + l , ...,J ) i s Ds. I n order f o r t h e (Lip i = 1, factorization (5.20) t o e x i s t , Ds must be d e f i n e d t o be S i n c e Ds i s given i n terms of known q u a n t i t i e s , t h e proof i s complete. m Note t h a t t h i s s p l i t procedure i n c l u d e s t h e LU f a c t o r i z a t i o n (s = J ) and t h e UL f a c t o r i z a t i o n ( s = o ) a s s p e c i a l c a s e s . The g e n e r a l a l g o r i t h m f o r s o l v i n g MV = r i s given below where any q u a n t i t y n o t p r e v i o u s l y d e f i n e d i s assumed t o be zero. S o l u t i o n of MV = r (o 5 s I J ) \ I. Decomposition procedure A. Forward ( i = o , l , . . .,s-1) 1. Compute: 2. -- Ai Di - LiBi n o n s i n g u l a r Di, r i' r 1 + p + 1,...,ni + p + n ) to form i n t e r c h a n g e rows ( n i 3. I f D. i s s i n g u l a r Interchange s t r a t e g y : and new B Change i+19'i+l9*i+lq accordingly. i+l Compute: w = r -Lw i i i-1 i D - ~ Li+l 'i+l i B. Backward ( i = J , 1. Compute: 2. ...,s + l ) Di = *i - 'i+lci+l If Di i s s i n g u l a r , i n t e r c h a n g e rows ...,n i - q - n ) t o form n o n s i n g u l a r ( n i - q - 1, Di, and new B , C i A i l Change ri, accordingly. 3. C. wi = r i - 'i+lwi+1 S p l i t point 1. 11. Compute: Compute: Ds = As - LsBs - us+l Cs+l Solution routine A. S p l i t point 1. B. C. Dsvs = w Solve f o r vs: Backward (i = s-1, ...,o ) 1. Solve f o r vi: Divi Forward (i = s+l, ...,J ) 1. Solve f o r vi: Divi s = wi - Ei+lVi+l = wi ,- Ci~i-l f i-l The numerical example cited in Chapter 6 uses s = o exclusively. In some cases however it has been found that the split point CY s = (5+1)/2 yields good numerical results when the s = o, s = J decompositions are unsuitable. Chapter VI A Numerical Example: P l a n e Couette Flow I n t h i s c h a p t e r , e q u a t i o n s modeling p l a n e C o u e t t e f l o w a r e used a s a t e s t c a s e f o r some of t h e t h e o r y p r e s e n t e d p r e v i o u s l y . This p a r t i c u l a r example was chosen because f o r c e r t a i n s p e c i a l c a s e s t h e e x a c t s o l u t i o n i s known and can b e shown t o b e a n i s o l a t e d s o l u t i o n ( s e e Chapter 4 ) . The e q u a t i o n s a r e g i v e n i n F. K. Moore [ $ 1 a s a s i m i l a r i t y s o l u t i o n of t h e Navier-Stokes e q u a t i o n s d e s c r i b i n g t h e motion of a compressible, v i s c o u s f l u i d c o n t a i n e d between two p a r a l l e l surfaces. These w a l l s a r e i n r e l a t i v e motion and each i s k e p t a t a c o n s t a n t temperature. l- OW Fig, 6.1 The u s u a l n o - s l i p boundary c o n d i t i o n s a r e imposed a t t h e wall-s, and t h e v a r i a b l e s of i n t e r e s t a r e T ( t ) , t h e normalized t e m p e r a t u r e d i s t r i b u t i o n , and u ( t ) , t h e normalized x-component of v e l o c i t y . I n accordance w i t h o u r p r e v i o u s n o t a t i o n , t i s used a s t h e independent variable. The e q u a t i o n s w i l l be c o n s i d e r e d i n t h e f o l l o w i n g form: where p(T) is the viscosity-temperature relation. Isolated Solution In the case where K = o, the differential equations (6.2) reduce to with the same boundary conditions. These equations are easily integrated if we restrict the viscosity-temperature relation to be of the form An exact solution for ( 6 . 3 ) is shown to be The homogeneous, linearized problem becomes where and Substituting the exact solution from (6.4) into A ( t ) , it can be shown that the linearized problem has only the trivial solution, and thus, the solution ( 6 . 4 ) is isolated. 2. Numerical Procedure The Gap4 will be employed to approximate the differential equations in (6,2). in (6.2), Substituting @(T) = p -1 (T) and letting this two-point boundary-value problem can be written as: where In order to use the Gap4 scheme, the vector f f must be evaluated. Y Let a = (a b c d l T , then The Gap 4 difference approximation is given by Coupling these difference equations with the discrete boundary conditions, yields 4J+4 nonlinear equations which must be solved in order to evaluate V. ~ewton'smethod is employed to solve these equations. 2 1 Define the 4 x 4 matrices N (a), N (a) by af (a) 2 (a) N1(a) = aa , N (a)=-. aa Thus, again letting a = (a b c d) T and 2 N (a) = 0 , where 4, 4 ' a r e understood t o s t a n d f o r $ ( c ) , $ ' ( c ) respectively. ~ e w t o n ' smethod y i e l d s t h e following l i n e a r d i f f e r e n c e equations: 3 These d i f f e r e n c e equations may be p u t i n t o a form where t h e matrix, qV(vV), i s block t r i d i a g o n a l m a t r i x w i t h 2/2 zero structure. The s p l i t decomposition r o u t i n e (s=o) was used t o s o l v e each of t h e s e equations f o r (vW1 t h e i n i t i a l guess, 3. - vV). I n each c a s e considered, vO, was taken t o be Numerical r e s u l t s These computations were performed i n double p r e c i s i o n o n an I B M 360/65. Three s e p a r a t e c a s e s a r e r e p o r t e d here: i i ) 1. = -1, p(T) = T, i i i ) K = -1, p(T) = T 3/2 . The left-hand boundary c o n d i t i o n , T(o) * A, was taken to be ~ ( 0 = ) 1/2, I n each c a s e , t e n uniformly spaced n e t p o i n t s (J=9) were placed on t h e i n t e r v a l [ o , l ] . e q u a t i o n s (6.5), ~ e w t o n ' smethod w a s used t o s o l v e t h e (6.6) and t h e i t e r a t i o n procedure w a s terminated when t h e r e s i d u a l , I lI1(vV)I I , was a c c e p t a b l y s m a l l . The r e s u l t s of t h e s e c a l c u l a t i o n s a r e summarized i n Table 6.7. (w*) V [~ll(v~)II Time ( s e c .) .- (O,1) (-191) - 1 0 . 1 1 D 00 1 .7142D-02 2 .1169D-04 3 .5071D-10 0 , 1 1 1 1 D 00 1 .1072D 00 2 -54521)-02 3 .1363D-04 4 ,12421)-09 0 . 1 1 1 1 D 00 1 .1004D 00 2 .1157D-01 3 .7149D-04 4 -6233D-08 2.51 3.67 4.64 *Represents t h e v i s c o s i t y r e l a t i o n p(T) = T~ Table 6.7 Note t h a t i n each c a s e t h e Newton i t e r a t e s e x h i b i t q u a d r a t i c convergence. An e x a c t s o l u t i o n i s known for c a s e i ) , ( o , l ) , s o a more d e t a i l e d i n s p e c t i o n of t h e numerical s o l u t i o n i s i n o r d e r . e x a c t s o l u t i o n f o r K = o, p(T) = T The is 3/4 3/8 + [l-t] + [l-t] . , vil , and y 2 ( t 1.) , vi2, a r e c o n s t a n t s The approximations t o y ( t ) 1 1 f o r a l l i: Recalling t h e exact s o l u t i o n (6.8), t h e corresponding e r r o r s a r e .907 D-05 and .453 D-05 r e s p e c t i v e l y . and (vi4 L e t Ivi3 - y4(ti) 1 = ei4COMPARISON OF NUMERICAL SOLUTION WITH EXACT SOLUTION Table ( 6 - 9 ) - y3(ti) ( = ei 3 The maximum e r r o r o c c u r s i n t h e approximation of t h e c o n s t a n t function yl(t) = ;(t), thus For t h i s example, t h e Gap4 scheme y i e l d s a n e x t r e m e l y a c c u r a t e numerical s o l u t i o n on a comparatively s p a r s e n e t . Also, ~ e w t o n ' s method e x h i b i t s q u a d r a t i c convergence once a n "acceptable" i n i t i a l g u e s s h a s been found. Appendix A Ordinary D i f f e r e n t i a l Equations of Higher Order I n a r e c e n t paper H, 0. Kreiss [ 7 ] examined d i f f e r e n c e schemes used t o approximate l i n e a r two-point boundary-value problems: where % are nxn m a t r i c e s and y,F,g a r e m-vectors. a t r e a t m e n t of e q u a t i o n s (A.la,b) the equations (l.la,b) A t f i r s t glance, seems t o be more g e n e r a l t h a n considered i n Chapter 1. an i n d i c a t i o n h e r e t h a t t h e t h e o r y of Kreiss However, we g i v e [ 7 ] and t h e t h e o r y developed i n Chapter 1 a r e a c t u a l l y concerned w i t h t h e same c l a s s of problems. A s an approximation of ( A . l a ) , Kreiss c o n s i d e r s d i f f e r e n c e schemes of t h e form i = r,.. .,J-s where D + i s t h e forward d i f f e r e n c e o p e r a t o r (D+ viWl 1 = ?; (vi - v ~ - ~ ) ) . The d i s c r e t e boundary c o n d i t i o n s a r e w r i t t e n where D- vi = D+ vi-l. We w i l l c o n s i d e r o n l y t h e c a s e i n which ( A . 2 ) i s "as compact a s possible", With t h i s r e s t r i c t i o n , (A.2) and (A.3) t h a t i s , r+s=n. r e p r e s e n t m(J+l) e q u a t i o n s i n m(J+l) unknowns, vi, i = 0, . . .,J . t h i s c a s e , t h e d i f f e r e n c e e q u a t i o n (A.2) can b e w r i t t e n a s n-1 h L~ vi = A ( t i , h ) ~ + vi-r + 1 ~ ~ k ( t +~ v,i-rh )= ~F~~ k=o i = r , . . .,J'S where n (ti,h) i s a n mxm m a t r i x and r e p l a c e s t h e d i f f e r e n c e o p e r a t o r So (h) The boundary-value problem (A.la,b) can b e w r i t t e n a s a f i r s t - o r d e r system by d e f i n i n g z k+l (t) = y(k)(t), Thus, a f i r s t - o r d e r system of e q u a t i o n s e q u i v a l e n t t o (A.la,b) n-1 dzn(t) dt dzk (t) dt + 1 q ( t ) zk+l ( t ) = F ( t ) k=o = z k+l ( t ) k = ~,...,n-1 is For It remains t o show t h a t , i n some s e n s e , (A.2), (A.3) can be c o n s i d e r e d as a n approximation of ( A b 5 a , b , c ) . Analogous t o t h e d e f i n i t i o n s (A.5b), d e f i n e t h e m-vectors, k k = 1, i' w ...,n , by S u b s t i t u t i n g t h e s e q u a n t i t i e s i n t o (A. 2) and (A. 3 ) , we g e t W n W i-r k+l i-r (A. 7a) and (A. 7b) Kreiss d e f i n e s t h e d i f f e r e n c e scheme (A.4) t o be c o n s i s t e n t w i t h n-l Thus, i f t h e d i f f e r e n c e scheme ( A . 4 ) i s c o n s i s t e n t , then c l e a r l y (A.7a), (A.6b), and (A.7b) form a c o n s i s t e n t approximation of (A05a,b,c). (N = m(nJ Thus, w e have N e q u a t i o n s i n N unknowns - $[n2 - n ] ) ) , and t h e o n l y d i f f e r e n c e between (A. 7 a , b ) , (A.6b) and t h e approximations c o n s i d e r e d i n Chapter 1 i s t h a t t h e k v a r i a b l e wi i s o n l y d e f i n e d f o r i = o , ...,J-k+l. 1- 2 In order to mimic Chapter 1, we introduce T ( n +n) new equations. This can be formally done by expanding the range of i in (A.7a) from J-s to J+r-1 and the range of each i in (ADGb) from J-k to J-1, The only restriction we place on these new equations is that be consistent with the differential equation (A.5a). In this way equations (A.ga,b,c) define a consistent difference approximation of (A.5a,b,c) where w1 = V, V being a solution of (A.2), (A.3). Thus, a stable solution, {w1,w2 ,...,W"], of the difference approximations (~.9a,b,c) yields a stable solution, V, of (A.2) and (A.3). In addition, the theory presented in Chapter 1 provides a necessary condition for stability where no restrictions are placed on the form of the difference approximations. References [I] S. D. Conte, he Numerical S o l u t i o n of L i n e a r Boundary-Value ~ r o b l e m s " , SIAM Review, Volume 8, No. 3, (1966), pp. 309-321. [2] L. Fox, "Numerical S o l u t i o n of Ordinary and P a r t i a l D i f f e r e n t i a l Equations", Addison-Wesley, London, 1962. [3] T. R. Goodman and G. N. Lance, "The Numerical I n t e g r a t i o n of Two-Point Boundary-Value Problems", Mathematics of Computations, Volume 1 0 , (1956), pp. 82-86. [4] E. I s a a c s o n and H. B. K e l l e r , "Analysis of Numerical Methodsv', John Wiley, New York, 1966. [5] H. B. K e l l e r , "Accurate D i f f e r e n c e Methods f o r L i n e a r Ordinary D i f f e r e n t i a l Systems S u b j e c t t o L i n e a r C o n s t r a i n t s " , SIAM J o u r n a l on Numerical A n a l y s i s , Volume 6, No. 1, (1969), pp. 8-30. [6] H. B. K e l l e r , "Accurate D i f f e r e n c e Methods f o r Nonlinear TwoP o i n t Boundary Value Problems", t o appear. (71 H. 0. Kreiss, " ~ i f f e r e n c eApproximations f o r Boundary and Eigenvalue Problems f o r Ordinary D i f f e r e n t i a l ~ q u a t i o n s " , Mathematics of Computation, Volume 26, No. 119, (1972), pp. 605-624. [8] F. K. Moore, e d i t o r , " ~ h e o r yof Laminar ~ l o w s " , P r i n c e t o n U n i v e r s i t y P r e s s , P r i n c e t o n , 1964, pp. 171-175. [9] V. P e r e y r a , " I t e r a t e d Deferred C o r r e c t i o n s f o r Nonlinear Boundary Value Problems", Numerische Mathematik, Volume 11, (1968), pp. 111-125.