Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Difference methods for ordinary differential equations with applications to parabolic equations Doedel, Eusebius Jacobus 1976

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1976_A1 D63.pdf [ 5.6MB ]
Metadata
JSON: 831-1.0080139.json
JSON-LD: 831-1.0080139-ld.json
RDF/XML (Pretty): 831-1.0080139-rdf.xml
RDF/JSON: 831-1.0080139-rdf.json
Turtle: 831-1.0080139-turtle.txt
N-Triples: 831-1.0080139-rdf-ntriples.txt
Original Record: 831-1.0080139-source.json
Full Text
831-1.0080139-fulltext.txt
Citation
831-1.0080139.ris

Full Text

D I F F E R E N C E M E T H O D S F O 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 W I T H A P P L I C A T I O N S T O P A R A B O L I C E Q U A T I O N S by Eusebius Jacobus Doedel B . S c . , U n i v e r s i t y of B r i t i s h C o l u m b i a , 1971 M . S c , U n i v e r s i t y of B r i t i s h C o l u m b i a , 1973 A Thes i s Submitted i n P a r t i a l Fu l f i lmen t of the Requi rements for the Degree of Doctor of Phi losophy i n the Institute of A p p l i e d Mathemat ics We accept this thesis as conforming to the r equ i red standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A January, 1976 In p r e s e n t i n g t h i s t h e s i s in p a r t i a l f u l f i l m e n t o f the r e q u i r e m e n t s f o r an advanced degree at the U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by the Head o f my Department o r by h i s r e p r e s e n t a t i v e s . It i s u n d e r s t o o d tha t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Department o f t ^l/^/jC f The U n i v e r s i t y o f B r i t i s h C o l u m b i a 2075 Wesbrook P l a c e Vancouver, Canada V6T 1W5 A B S T R A C T The f i r s t chapter of the thesis i s concerned with the c o n s t r u c t i o n of finite d i f f e r e n c e approximations to boundary value p r o b l e m s i n l i n e a r nth o r d e r o r d i n a r y d i f f e r e n t i a l equations. T h i s c o n s t r u c t i o n i s based upon a l o c a l c o l l o c a t i o n p r o c e d u r e with polynomials, which i s equivalent to a method of undetermined c o e f f i c i e n t s . It is shown that the c o e f f i c i e n t s of these finite d i f f e r e n c e approximations can be e x p r e s s e d as the deter-minants of m a t r i c e s of r e l a t i v e l y s m a l l dimension. A basic theorem states that these approximations a r e consistent, p r o v i d e d only that a c e r t a i n n o r m a l i z a t i o n factor does not vanish. T h i s i s the case for compact di f f e r e n c e equations and for d i f f e r e n c e equations with only one c o l l o c a t i o n point. The o r d e r of c o n s i s t e n c y may be i m p r o v e d by suitable choice of the c o l l o c a t i o n points. S e v e r a l examples of known, as w e l l as new d i f f e r e n c e approximations a r e given. A p p r o x i m a t i o n s to boundary conditions a r e a l s o treated i n d e t a i l . The s t a b i l i t y theory of H. O. K r e i s s i s applied to investigate the s t a b i l i t y of finite d i f f e r e n c e schemes based upon these approximations. A number of n u m e r i c a l examples a r e a l s o given. In the second chapter i t i s shown how the c o n s t r u c t i o n method of the f i r s t chapter can be extended to i n i t i a l value p r o b l e m s for systems of l i n e a r f i r s t o r d e r o r d i n a r y d i f f e r e n t i a l equations. Sp e c i f i c examples a r e 'included and the well-known st a b i l i t y theory for these d i f f e r e n c e equations i s s u mmarized. It is then shown how these d i f f e r e n c e methods may be applied to l i n e a r p a r a b o l i c p a r t i a l d i f f e r e n t i a l equations i n one space v a r i a b l e after f i r s t d i s c r e t i z i n g i n space by a suitable method f r o m the f i r s t chapter. -111-The s t a b i l i t y of such d i f f e r e n c e schemes for p a r a b o l i c equations i s investigated using an eigenvalue-eigenvector a n a l y s i s . In p a r t i c u l a r , the effect of v a r i o u s approximations to the boundary conditions i s c o n s i d e r e d . The r e l a t i o n of this a n a l y s i s to the s t a b i l i t y theory of J. M. V a r a h i s indicated. N u m e r i c a l examples a r e also included. - i v -T A B L E O F C O N T E N T S A B S T R A C T i T A B L E O F C O N T E N T S i v A C K N O W L E D G M E N T S v Chapter I. B O U N D A R Y V A L U E P R O B L E M S 1 1. Introduction 1 2. The Cons t ruc t ion of F in i t e Difference Approx imat ions 6 3. Effect of the Choice of Co l loca t ion Points 19 4. Examples of F in i t e Difference Approx imat ions 21 5. F in i t e Difference Approx imat ions to Boundary Condit ions 29 6. Examples of Approx imat ions to Boundary Condit ions 34 7. The Stabi l i ty of F in i t e Difference Schemes 37 8. N u m e r i c a l Example s 46 II. I N I T I A L V A L U E P R O B L E M S 54 1. Introduction 54 2. In i t i a l Value P rob lems in Ord ina ry Di f fe ren t ia l Equations 55 3. Examples of F in i t e Difference Approx imat ions 58 4. The Stabi l i ty of F in i t e Difference Approx imat ions to In i t i a l Value P rob l ems 65 5. In i t i a l Boundary Value P rob lems 71 6. The Stabi l i ty of F in i t e Difference Approx imat ions to In i t i a l Boundary Value P rob lems 76 7. N u m e r i c a l Examples 91 B I B L I O G R A P H Y 117 -v-A C K N O W L E D G M E N T S The author wishes to thank P r o f e s s o r J . M. V a r a h f o r stimulating and guiding this r e s e a r c h . He i s al s o thankful to his wife A d r i e n n e f o r her patience. F i n a l l y , he is indebted to M r s . V i v i a n Davies f o r her excellent typing of the m a n u s c r i p t and to D r . Fau s t o M i l i n a z z o f o r pro o f r e a d i n g the ma n u s c r i p t . -1-Chapter I B O U N D A R Y V A L U E P R O B L E M S 1. I N T R O D U C T I O N In this chapter the c o n s t r u c t i o n of finite d i f f e r e n c e approximations to the l i n e a r d i f f e r e n t i a l equation (1.1) Ly(x) = y ( n ) ( x ) + n j 1 a k ( x ) y ( k ) ( x ) = f(x), 0 < x < 1 , k=0 i s investigated. Define a m e s h S. = {x. : 0 = x„ < x, < ••• < x T = 1 } with mesh-h L l 0 1 J J spacings h. = x. - x. ^ , (1 < j < J). L e t h = max h.. It w i l l always be 1 j > -J assumed that m i n > h for some constant K. F o r any function w(x) J on S, let w. = w(x.). h J J F i n i t e d i f f e r e n c e approximations to the d i f f e r e n t i a l equation (1.1) at the point x = x^ a r e assumed to have the f o r m s. J (1. 2) L h u . ^ I <L . u j + . = f., k Q < j < J - n + k Q , k Q > 0 , i=-r. J where f. i s some app r o x i m a t i o n to f.. Thus for each m e s h S, there a r e J J Im-p r e c i s e l y J - n + 1 such equations. The width co.. of each fini t e d i f f e r e n c e a p p r o x i m a t i o n i s defined as co. s r. + s. + 1. If co. = n + 1 then the J J J J a p p r o x i m a t i o n i s said to be compact. The assumption i s made that co. < co f o r some integer co that does not depend on S, . j — max b max r h The l o c a l truncation e r r o r i s given by T . = L^y^ - f.. If there i s a n s constant C independent of h such that |T.|<Ch with n > 0, and n as J s s high as possible, then the difference approximation is said to be consistent with the differential equation and the order of consistency is equal to n . s To specify y(x) completely it is necessary to impose boundary conditions, which are assumed to take the form n k(0) (1.1a) B k(0)y(0)= £ b ^ .(0)y ( l )(0) = b k(0), l < k < n 0 , i=0 and n k ( l ) (1.1b) B k ( l ) y ( l ) = b k > . ( l ) y ( l ) ( l ) = b k ( l ) , l < k < n - n 0 , i=0 where n k(0) < n and n_r.(l) n • It will always be assumed that the differential equation (1.1) subject to the boundary conditions (l.la,b) admits a unique solution that is as many times continuously differentiable as needed. Similarly, to define a finite difference scheme that determines the approximations u. to y. completely one has approximations to the boundary conditions (1.1a, b) that look like sk<0> (1.2a) B h j k ( 0 ) u 0 E E YJ d k , i ( 0 ) u i = b k ( 0 ) ' i < k < n 0 ' i=0 and -3-0 (1.2b) B h ( k ( l ) u j S - 2 d k > . ( l ) u J + i = b k ( l ) . l < k < n - n 0 . i= - r k ( l ) Here b k(0) and b k ( l ) are approximations to b k(0) and b k ( l ) respectively. The definitions of truncation e r ro r and consistency are s imi lar to the previous definitions and the equations are compact if s k(0) = n k (0 ) and r k ( l ) = n k ( l ) . A finite difference scheme is said to be (accurate) of  order s_ if the approximate solution Uy (0 <j < J), is uniquely defined by (1.2) and (1.2a,b) and |y^-u. | < K ^ h S for some constant that does not depend on h. Numerous texts and papers are partly or entirely concerned with finite difference approximations to (1.1) and ( l . l a , b ) . In part icular , the books by Collatz (I960) and Kel ler (1968) need to be mentioned here. The general theory of difference approximations to systems of boundary value problems is contained in papers by Grigorieff (1970) and Kre iss (1972). A n extensive survey of finite difference methods for boundary value problems can be found in a paper by Kel ler (1975). Var ious finite difference approximations of the form (1.2) appear in the l i terature. To the second order differential equation (1.3) y"(x) + a^xjy'fx) + a ° (x )y (x ) = f(x) , the simplest difference approximation is perhaps given by (1.4) -4-Here D. u. = (u h J J+1 -u. ,) / 2h and the mesh is uniform, i.e., h. = h. The order of consistency is equal to two, as is the order of accuracy of a difference scheme based upon (1.4), provided that the approximations to the boundary conditions are consistent. Frequently i t is desirable to use equations for which the order of con-sistency is more than two. One way of obtaining such equations is by allowing the width of the approximation to be greater than n+1. For example, a fourth order finite difference approximation to (1.3) is The width of this approximation is equal to five. If a difference approxi-mation is not compact then an entire difference scheme cannot be based upon that equation alone. Near the boundaries other equations need to be used. In (1.5) for example, L^u. is not defined if j = 1 or j = J-1 unless periodicity is assumed. Hence other approximations need to be used here. Examples of these are given by Shoosmith (1975). In this chapter such extra boundary conditions will be discussed i n Sections 7 and 8. It is not necessary to increase the width of the approximation to obtain higher order equations, for one can also require more "local information" to be available. In this case each difference equation involves the values of the coefficient functions and the inhomogeneous term at more than one point. For example, a fourth order approximation to (1.5) L h u . EE (1.6) y"(x) + a°(x)y.(x) = f(x) -5-is given by (1.7) D*u. + K h ( f ) (a°u.) = K h ( f ) f . , with K h(c)w. s i|£. w.^ + cw. + 1 f £ w . + 1 . The construction and analysis of such higher order compact approxi-mations is also included in this chapter. Usually the statement that a certain difference approximation has a given order of consistency can be justified by simply Taylor expanding. Systematic methods of deriving finite difference equations appear infre-quently in the literature however, unless one considers finite element methods such as Galerkin, Ritz and collocation procedures with piecewise polynomials as finite difference methods. (For a discussion of these, see e.g., Varga (1966), Ciarlet Schultz and Varga (1967), Russell and Shampine (1972), deBoor and Swartz (1973), as well as the books of Strang and F i x (1973) and Schultz (1973)). Birkhoff and Gulati: (1974) compare various compact equations, but no attempt is made to methodically construct such compact difference approximations. Swartz (1974) discusses a procedure for deriving difference approximations that resembles the procedure given in this thesis, but no complete analysis is given. The main purpose of the present chapter is then to derive a large class of finite difference approximations to the problem (1.1) and to give their order of consistency. This is done in the next two sections. In Section 4 a number of examples is given, including known as well as new -6-d i f f e r e n c e equations. In S e c t i o n 5 the c o n s t r u c t i o n of f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s to the b o u n d a r y c o n d i t i o n s i s c o n s i d e r e d . E x a m p l e s of these appear i n S e c t i o n 6. In S e c t i o n 7 the t h e o r y of K r e i s s (1972) i s a p p l i e d to i n v e s t i g a t e the s t a b i l i t y of d i f f e r e n c e s chemes b a s e d upon the d i f f e r e n c e a p p r o x i m a t i o n s g i v e n i n this c h a p t e r . F i n a l l y , i n S e c t i o n 8 the r e s u l t s of some n u m e r i c a l e x p e r i m e n t s a r e g i v e n . 2. T H E C O N S T R U C T I O N O F F I N I T E D I F F E R E N C E A P P R O X I M A T I O N S To c o n s t r u c t a f i n i t e d i f f e r e n c e a p p r o x i m a t i o n of the f o r m (1.2) to the d i f f e r e n t i a l equation (1.1) one m a y p r o c e e d as f o l l o w s . F o r g i v e n i n d e x j l e t z , < z_ < • • • < z be p o i n t s i n the i n t e r v a l Tx. ,x.. 1, J 1 2 m r L j - r j + s J s a t i s f y i n g | z . . - z , | >—- h, i # k\ (To s i m p l i f y the n o t a t i o n the s u b s c r i p t 1 K j denoting dependence on the p a r t i c u l a r d i f f e r e n c e a p p r o x i m a t i o n w i l l be o m i t t e d h e r e , so that, e.g., z. . b e c omes z.). J> 1 1 L e t P_j denote the space of a l l p o l y n o m i a l s w i t h degree not exceeding d = r + s + m - l . A s s u m e that f o r some i n t e g e r i, (-r < l< s), SL a p o l y n o m i a l p (x) e P^ i s d e t e r m i n e d u n i q u e l y by the r + s + m equations r < i < s , i + t 1 < i < m . (The equations (2.2) w i l l be r e f e r r e d to as " c o l l o c a t i o n " equations.) T h e n the equation (2. 1) P (x )_ = u j+i (2.2) Lp (z.) f(z.) (2.3) r e l a t e s the quantities u j ^ » ( - r < i < s ) , by an e x p r e s s i o n (difference equation) of the f o r m m (2.4) u. = V d.u.,. = V e.f(z.), l = - r i = l In o r d e r to be able to exhibit e x p l i c i t e x pressions for the coe f f i c i e n t s d^ and e^ i n (2.4), let the polynomials w 1(x), (0 < i < d), f o r m a basis of P^, so that one can w r i t e P*(x) = Yi c iw 1(x), i=0 A l s o , define a m a t r i x by (2.5) D. w ( X j _ r ) 0, w (x. , ) j + s L w ( z ^ Lw^(z^) w (x. ) j _ r ' w (x. , ,) w (x. , , ) x j-r+1 v j-r+1' w (x. , ) L w ^ z ^ ) Lw"(z ) Lw*(z ) m m d. . w (x. ) v J-r' w (x. , ) j + s' u w (x. , ,) u. , , x j-r+1' J-r+1 U j + B f ( z A ) f ( z 2 ) L w u ( z ) f(z ) m' v m If the relations (2.1), (2.2) and (2.3) a r e s a t i s f i e d then det(D^) = 0. Ev a l u a t i n g the determinant by the l a s t column and int r o d u c i n g a n o r m a l i z -ing factor E one finds that the coe f f i c i e n t s d^ and e^ i n (2.4) are given by (2.6) d i = c o f L [ V i ] / E ' -8-and (2.7) e. = - c o f L [ f ( z . ) ] / E where cof-^f • ] i s the cofactor of the given element i n the m a t r i x D ^ , wi th the operator L as i n equation (1.1) . It i s convenient to define the n o r m a l i z i n g factor as m E = " E c o f L ' where i = l L Q y ( x ) EE y ^ x ) . If the bas i s functions w 1(x) are chosen such that (2.8) w ( x j . r + k ) = 6 i k f o r 0 < i < r + s , and (2.9) w ( x j + k ) = 0 f o r r r + s + l < i < d , - r < k < s , then (2. 10) d m i " E j r+i . * , T r+s + l \ L w . ( z ^ . , lLw.u , ( z ^ , ; , , j r+i . > T r+s + l , . L w (z ) L w (z ) m x m Lw^(z , L w d ( z ) m , - r < i < s , and -9-( 2 . 1 1 ) e, Lzll m+i+1 E T r+s+1 , . Lw ( z ^ T r+s+ 1 , . L w (z. ,) l - l j r+s+ 1 , . ( z i + 1 ) L w d ( z Lw (z^ ^ ) L w d ( z . + 1) Lw r^ S"'"^(z ) " ' Lw^(z ) m m , 1 < i < m With the above choice of b a s i s functions, one can al s o express the n o r m a l i z i n g factor as ( 2 . 1 2 ) E = ( -1) m , j r+s + 1, > . . . T r+s+m - 1 , . 1 L Q w ( z t ) L Q w ( Z l) , T r+s + 1, . . T r+s+m - 1 , . 1 L „ w (z ) L n w (z 0 v m/ 0 v rrr If m = 1 then E = -1 and e, = 1. F u r t h e r , the coe f f i c i e n t s d. and e. a r e 1 i i r e l a t e d by m d i = L L , w r + 1 ^ z k ^ e k ' " r - 1 - S ' k=l Remark. The above f o r m u l a s for d^ and e.. a r e e s s e n t i a l l y intended to show the f o r m of these c o e f f i c i e n t s . In p a r t i c u l a r , note the type of products of co e f f i c i e n t functions that w i l l appear. F o r actual a p p l i c a t i o n to n u m e r i c a l problems, one may f i r s t evaluate the determinants s y m b o l i c a l l y , c o l l e c t i n g the v a r i o u s m u l t i p l i c a t i v e and additive constants. T h e r e a f t e r , the n u m e r i c a l evaluation of the d^ and e^ may be or g a n i z e d i n such a way, that these coefficients a r e obtained i n the s m a l l e s t p o s s i b l e 10-number of a r i thmet ic operations. B a s i s functions w 1(x) sat isfying (2.8) and (2.9) are eas i ly constructed. In fact, i t i s convenient to assume that the w 1(x) are given by (2. 13) w 1(x) = r + s rr k=0 k ^ i (x -x . ,. ) j - r + k r+s (x. , . - x . , . ) k = 0 ' J -r+i J-r+k' k#i 0 < i < r + s , and (2. 14) w r + S + 1 ( x ) = m-1 TT(*-tt k=l k^a  m - 1 TT<t_-tk) k=l 1 M i r+s TT (x-x. k=0 j - r + k ' r+s Lo 1 J " r + k 1 < i < m - 1. Here the t., (1 < i < m-1 ) , are dis t inct points i n the i n t e rva l Tx. , x . , 1. I — — r L j - r j+s J It i s not diff icult to check that the points t\ can be chosen i n such a way that the polynomials w 1 (x), (1 < i < d), f o rm a l i n e a r l y independent set. The procedure leading up to the definit ion of the coefficients d^ and e^  i s equivalent to a method of undetermined coeff icients . To i l l u s t r a t e this cons ider the following example. L e t n = 2 , r = s = 1 and m = 2; i . e., we are concerned wi th the cons t ruc t ion of a difference approximat ion of the fo rm -11-(2.15) d _ ! u j _ l + d 0 u j + d l u j + i = e i £ ( z i ) + e 2 f ( z 2 ) ' to the differential equation Ly(x) = y (x) + a (x)y (x) + a (x)y(x) = f(x) . This can be accomplished by requiring that (2.15) be satisfied for all polynomials p(x) e P^. (Replacing f(z^) by Lp(z.)). With w1(x), ( 0 < i < 3), as in (2. 13) and (2. 14), this requirement leads to the following linear system of equations: 0 " 0 CO 0 ' Assume that the rank of the system is equal to four, so that one can rewrite the equations as 1 0 0 . d - l Lw^(z^) e^ 0 1 0 . -Lw'(z ^) d o L w1 ( z 2 ) e 2 0 0 1 . 2 -Lw (z^) d l L w2 ( z 2 ) e 2 0 0 o • 3 -Lw (z^) _ e l _ 3 Lw ( z 2 ) e 2 where the matrix on the left is nonsingular, i.e., Lw (z.)±0. Partition 1 0 0 0 1 0 0 0 1 0 0 0 •Lw^(z -Lw (z^) 2 -Lw (z p 3 -Lw (z ,) L •Lw^(z. -Lw (z 2) -Lw^(z 2) t h i s m a t r i x i n t o I A ' 0 B 12-a s i n d i c a t e d . T h e n I A 0 B -. -1 I - A B 0 B " , a n d -1 1 b y e l e m e n t a r y l i n e a r a l g e b r a B = ^g^g ^ C , w h e r e C i s the m a t r i x o f T -1 -1 c o f a c t o r s o f B . In t h i s s i m p l e e x a m p l e B = L w (z )^ H e n c e = [ L w ° ( z 2 ) - L w ° ( z 1 ) L w 3 ( z 2 ) / L w 3 ( z 1 ) ] e 2 , d Q = [ L w 1 ( z 2 ) - L w 1 ( z 1 ) L w 3 ( z 2 ) / L w 3 ( z 1 ) ] e 2 , a n d [ L w 2 ( z 2 ) - L w 2 ( z 1 ) L w 3 ( z 2 ) / L w 3 ( z e 2 > 3 3 [ - 1 / L w ( z 4 ) ] L w ( z 2 ) e 2 -L e t e 2 = • L w (z^) E T h e n the c o e f f i c i e n t s a r e e q u a l to [ L w " ( Z j ) L w ^ ( z _ ) L w ^ ( z 2 ) L w 3 ( z p ] / E d p = [ L w ^ ( z j ) L w ^ ( z _ ) - L w ^ ( z 2 ) L w 3 ( z ^ ) ] / E , [ L w 2 ( z 1 ) L w 3 ( z 2 ) L w 2 ( z 2 ) L w 3 ( z ^ ) ] / E = L w (z_) / E , = - L w (z )^ / E w h i c h i s i n a g r e e m e n t w i t h the e q u a t i o n s ( 2 . 1 0 ) a n d ( 2 . 1 1 ) . M o r e g e n e r a l l y the a b o v e p r o c e d u r e i s w e l l - d e f i n e d p r o v i d e d tha t the m a t r i x -13-T r + s + 1. > T r+s+m-1/ L w (z^) • • • L w (z j r + s+1, . T r+s+m-1, > L w (z ) . . . L w (z ) v m m has rank equal to m-1. In p a r t i c u l a r i f the n o r m a l i z i n g factor E i s non-z e r o then the m a t r i x has rank m-1 for a l l s m a l l enough h. With the definitions given this far, i t i s p o s s i b l e to state the b a s i c theorem of this chapter. T h e o r e m 2. 1 L e t the co e f f i c i e n t s d. and e. and the n o r m a l i z i n g factor E be given by (2. 10), (2. 11) and (2. 12) r e s p e c t i v e l y . A s s u m e that E ± 0 and that r + s > n. If h is s m a l l enough then at l e a s t n + 1 of the coefficients d. are nonzero and the o r d e r of c o n s i s t e n c y of the finite l 1 d i f f e r e n c e approximation (2.4) i s at l e a s t equal to r + s + m - n . Proo f . Define (2. 16) d. = c o f T [ u . , . l / E, - r < i < s , i L L j + i J ' - -(2. 17) e. = - c o f T [f(z.)] / E, 1 < i < m , J0 m so that Y ei - 1-i = l F i r s t i t w i l l be shown that i f E # 0 then at l e a s t n + 1 of the coefficients d. a r e nonzero. F o r suppose on the c o n t r a r y that d. + 0, 1 _ k l < k < k ^ < n + 1 and d^ = 0 otherwise. L e t , 14-p(x) EE 1 L k=i k-l q(x), where q(x) e P , i s chosen such that p(x) 1 has degree n and leading c o e f f i c i e n t — . B y assumption n < r + s, so c e r t a i n l y n < d. Hence by c o n s t r u c t i o n m I d i p j + i - E ei LoP< zi> = °-i = - r i = l But f r o m the definition of p(x) i t follows that ^ d^ Pj+j_ = 0 and LQP(X)=1, m _ i = - r T h e r e f o r e ^ e^ = 0, which i s a c o n t r a d i c t i o n . So i f E + 0 then n + 1 or _ i = l _ m o r e d. a r e nonzero. Since d. / d. = 1 + 0(h), the same is true for the l l ' l d^, p r o v i d e d h i s s m a l l enough. The truncation e r r o r T. i s defined as J m s m r j - V j ~ I e i f ( z i > = I V j + i - I e i £ ( z i ) • i = l i = - r i = l where y(x) i s the solution to the d i f f e r e n t i a l equation (1.1) subject to the boundary conditions (1. l a , b). T a y l o r expand to get T. = d y<w k=0 m r- d i=-r £=1 y (k) k! k=0 + V d.(x. ,.-x. ,d+l y ( d + 1 ) ( - ) ™ (d+l)! i = - r - - I ^ i = l (z.-x.) . 1 J d+l y( d + 1 ) ( t i : (d+l)! d k=0 -15-^ V A i xd+1 y [Si} v 1 _ i _ r / ^+1 + > d.(x.,.-x.) — . ,, ,, , — - > e. — , , , ... L(z.-x.) , t-i l j+i j (d+1)! Li I (d+1)! v I y i=-r i = l for some s. between x. and x., . and t- between x. and z. . B y c o n s t r u c t i o n i J J+i i J i the quantity between square brackets equals zero for each k. F r o m the definitions of the coefficients d. and e. and the basis functions w 1(x), 1 1 v together with the underlying assumption that ^ < tu < h, (1 <_ j <_ J), i t follows that there a r e constants C^, and that do not depend on h, such that d^[ < C^h n , - r < i < s , e i I — ^ 2 ' 1 < i < m , and T / ^d+1I ^ „ ,d+l-n L(z-x.) < C.h , z. < z < z J 1 — 3 1 — — m With these estimates the second c o n c l u s i o n of the theorem follows imme dia tely. B y means of two examples the n e c e s s i t y of the condition that E be nonzero w i l l be i l l u s t r a t e d . E xample 2. 2 C o n s i d e r the case where n = 2, m = 2, r = 1, and s = 2. Thus we a r e co n c e r n e d with the c o n s t r u c t i o n of a finite d i f f e r e n c e I I 1 ' 0 a p p r o x i m a t i o n to the p r o b l e m L y = y (x) + a (x)y (x) + a (x)y(x) = f(x) at the point x., that involves the approximations u. <, u., u.,., u., 0 to y. J F F J - l J J+1 J+2 y J - l ' y-' y-iJ-1 y-x?- T h i s i s a c c o m p l i s h e d by the p r o c e d u r e of this section, J J' i ' J+^ -16-with two c o l l o c a t i o n points. F o r s i m p l i c i t y i t i s also assumed that the m e s h i s uniform. The basis functions w X(x) a r e defined as i n (2. 13) and (2. 14). In this c ase they can be w r i t t e n as w°(x) = -(x-x.)(x-x. + 1 ) ( x - x . + 2 ) / 6h , w^x) = (x-x j_ 1)(x-x. + 1 ) ( x - x . + 2 ) / 2h , w 2(x) = -(x-x._ 1)(x-x.)(x-x j + 2) / 2h , 3 , , / W W V / / ! 3 w (x) = ( x - X j _ 1 ) ( x - x ^ ) ( x - X j + 1 ) / 6h , 4 4 w (x) = (x-x^._ 1)(x-x j)(x-x. + 1 ) ( x - x . + 2 ) / h . L e t the c o l l o c a t i o n points be defined by 1 = X j - 1 + e i h a n d z 2 = X j - 1 + ^ 2 h The n o r m a l i z i n g factor E i s i n this example given by E = L^w ( z 2 ) -L Q W ^ ( Z ^ ) = ( < _ , 2 - < r ; ^ ) ( £ 2 + 3) / h 2 . The c o l l o c a t i o n points z^ a r e always assumed to be distinct, so £ 2 ± S_^. Hence E - 0 i f and only i f ^  + i_, 2 = 3, i . e . , when the points z^ and z 2 a r e p l a c e d s y m m e t r i c a l l y about the midpoint (x^ + xj_|__) / 2,x, which would seem the m o s t natural placement of these points. One can not a l l e v i a t e the p r o b l e m by redefining the n o r m a l i z i n g factor E . F o r example, with E = l / h , ^ = 1 and = 2 the f o r m u l a d. = h' l L Q w ( Z l) L Q w ( Z l) T \ T 4. . L Q w ( z 2 ) L Q w ( z 2 ) -1 < i < 2 , - 1 7 -2 2 2 2 g i v e s d _ 1 = -2/h , d Q - 6/h , = - 6 / h a n d d 2 = 2/h . F u r t h e r , — 2 4 — 2 4 e^ = h L Q W ( Z ^ ) = -2 a n d e^ = -h L Q W ( Z ^ ) = 2, so t h a t f o r t h e p r o b l e m y (x) = f(x) the d i f f e r e n c e e q u a t i o n b e c o m e s (-2u. j . l + - 6u. + 1 + 2 u j + 2 ) /h* - 2f( x j + 1) - 2 f ( x . ) . D i v i d i n g b y 2 h one o b t a i n s ( " V l + 3 U j - 3 U j + l + U j + 2 ) / h 3 = ( f h + l ) " f ( X j } ) / h ' i n i I I w h i c h i s c o n s i s t e n t w i t h y (x) = f (x), r a t h e r t h a n y (x) = f ( x ) . E x a m p l e 2. 3 A s i n t h e p r e v i o u s e x a m p l e l e t n = 2, r = 1 a n d s = 2, b u t take m = 3. So t h e r e a r e t h r e e c o l l o c a t i o n p o i n t s . I n a d d i t i o n to the b a s i s f u n c t i o n s w 1 ( x ) , ( 0 < i < 4 ) , i n e x a m p l e (2.2) t a k e w (x) = ( x - x . •. ) ( x - x . ) ( x - x . A ) ( x - x . . ) ( x - x . _) / h , J - } J J + 2 J+1 i+C w h e r e x. _i = \ (x. + x. . ) . T h e n the p o l y n o m i a l s w 1 ( x ) , ( 0 < i < 5 ) , f o r m J+2 ^ J J + 1 a b a s i s o f P^. N o w E i s g i v e n b y E = L Q w 4 ( z 2 ) L Q w 5 ( z 2 ) 4 5 L Q w ( z 3 ) L Q w ( z 3 ) 4 * - 5 • L Q W ' ( Z J ) L Q W ( Z ^ ) 4 5 L Q w ( z 3 ) L Q W ( Z 3 ) 4 5 L , Q W ( Z ^ L Q W ( Z J ) L Q w 4 ( z 2 ) L Q w 5 ( z 2 ) I f t h e m e s h i s t a k e n to be u n i f o r m a n d i f t h e c o l l o c a t i o n p o i n t s a r e d i s t r i b u t e d s y m m e t r i c a l l y , i . e . , z ^ = x j _ ^ i . - £h, z 2 = x j X i ' z 3 = x j + - * 3 ? 5 t h e n a f t e r s o m e c o m p u t a t i o n i t i s f o u n d t h a t E = c£ (4£ - 3 ) / h , w h e r e c i s -18-some constant. Since the points z^ are dist inct , § cannot be equal to ze ro . Hence i f the points z.. a re placed s y m m e t r i c a l l y , then E = 0 i f and only i f _ = j ^3' = 0.866. R e m a r k . A necessary condit ion for the n o r m a l i z i n g factor E to be non-zero i s that r + s > n . F o r example, the width of a difference a p p r o x i -mat ion to a second order different ia l equation i s at least three. This should not be s u r p r i s i n g . To prove the a s se r t i on suppose r + s < n . Then d n w 1 (x) / dx 1 1 = 0 , (0 < i <.r + s), because the f i r s t r + s + l basis functions have degree r + s. Hence the f i r s t co lumn i n the determinant defining d.. consis ts of zeros for each i . Thus d.. = 0, (-r < i < s), and f r o m the proof of theorem (2. 1) i t follows that E = 0. Some sufficient conditions for E to be nonzero are given i n the following theorem. Theo rem 2.4 The n o r m a l i z i n g factor E i s nonzero i f m = 1 and r + s > n, or i f r + s = n . ( i . e . , when the difference equation is compact . ) P roo f . If m = 1 then E = - 1 . If- r + s = n, m > l and E = 0, then i t follows f rom (2.12) that there are constants c , (r + s < i < d), not a l l d 1 zero , such that the po lynomia l q(x) = ( __, c^w 1(x)^ e P ^ sat isf ies i=r+s+l q^n\z/) - c r + g > (1 < i < m ) . Thus the po lynomia l q(x) = q(x) - c _ + g x n / n ! satisfies q ^ ( z ^ ) =0. But this i s imposs ib le since q^n^(x) € P _ n _ ^  because r + s = n . R e m a r k . If the approximat ion i s compact then for each n the coefficients dj, are uniquely determined, independent of the number of co l loca t ion points . F o r example, i f n = 2 and r + s = 2 then d ^ V V V i 1 -19-— -2 — 2 d.g = and d^ = ^ (h J . h f ' s o ^h.at t n e approximat ion j j+1 j+l j j + 1 second der ivat ive is always given by u . , . - u . u. - u... , __±i 1 _ J L l i h . . h . J+1 J 3. E F F E C T O F T H E C H O I C E O F C O L L O C A T I O N POINTS It i s w e l l known, (see, e . g . , R u s s e l l and V a r a n (1975)), that co l loca t ion procedures wi th ce r t a in p iecewise polynomials for the n u m e r i c a l solution of boundary value problems have a higher order of accuracy i f Gauss ian points are used as co l loca t ion points . One w i l l expect, that a s i m i l a r statement applies to the finite difference equations d i scussed i n the previous sect ion. However , i t may be su rp r i s ing that the points z^, for which an extra high order of consis tency is attained, do not genera l ly coincide wi th the Gauss ian points . Never the less , these spec ia l points can be cha rac te r i zed r e l a t ive ly eas i ly as w i l l now be shown. Cons ider again the bas is functions wX(x) of P . , defined by 6 v r+s+m-1 3 equations (2.13) and (2.14). This l i n e a r l y independent set may be extended to f o r m a basis of P , , by adding a po lynomia l w (x) e r+s+m ^ } v ' •r^r+s+m that vanishes at the mesh points . Thus this po lynomia l has the f o r m -20-r+s+m, . k = l w (x) = m-1 r+s T T^-t,). FT j . r + k ) k=0 m-1 r+s r r t w v TT (t m-^. r t k) with t m * t k , (.1 < k < m-1), and t^_ + x., (-r<j< s). E x p a n s i o n of y(x) i n t e r m s of the polynomials w (x) gives r+s+m y(x) = I c,w k(x) + 0 ( h r + s + m + 1 ) . k=0 Substitution of this e x p r e s s i o n into the truncation e r r o r m T j s E V j + i • E e i f ( z i } y i e l d s i= l i = - r r+s+m k=0 m E d i w ( xj+i) - E e i L w <zi> i = - r i = i + 0 ( h r + s + m " n + 1 ) The quantity between square brackets equals z e r o for 0 < k < _ r + s + m - 1, Moreover, w r + s + m ( X j _ ^ ) = 0, (-r < i < s). Hence m V T r+s+m, . ~ ,,r+s+m-n+l> T. = -c , , > e.Lw (z.) + 0( h ) . j r+s+m u i i' i = l E s t i m a t e s as i n the proof of theorem (2. 1) again r e v e a l that T\ = 0 ( h r ^ " S ^ " m n ) . However, i t is also c l e a r f r o m the above expression, that an additional o r d e r of con s i s t e n c y i s obtained, i f the c o l l o c a t i o n points z. , , ,n r+s+m, » • j -i-i. • i. v. T r+s+m, , d w (x) _ , coinc i d e with points where L Q W (x) = — s—-. = 0. Thus the dx following has been shown: T h e o r e m 3. 1. L e t the c o e f f i c i e n t s d. and e. and the n o r m a l i z i n g f a c t o r E be g i v e n by (2.10), (2.11) and (2. 12) r e s p e c t i v e l y . A s s u m e that E + 0 ~ r+s+m m-1 r+s and that r + s > n. L e t w (x) = ~| |~ (x-t, ) 1 T (x - x . , ) and a s s u m e k=i K k=o J " r + K s*** -f s T in. that LpW ( z ^ = 0,(1 < i < m ) . Then the o r d e r of c o n s i s t e n c y of the f i n i t e d i f f e r e n c e a p p r o x i m a t i o n (2.4) i s at l e a s t e q u a l to r + s.+ m - n + 1. If m - 1 and the d i f f e r e n c e a p p r o x i m a t i o n i s compact, i . e . , r + s = n, ~ r+s+m n then w (x) = ~| ["(x-x. , ). I n t h i s c a s e t h e r e i s only one k=0 J _ r + c o l l o c a t i o n p o i n t f o r w h i c h a h i g h e r o r d e r of c o n s i s t e n c y i s obtained. If m = 1 and r + s > n then t h e r e a r e r + s + 1 - n p o s s i b l e c h o i c e s of thi s point. ( A s s u m i n g E ± 0.) F o r the c a s e m > 1, r + s = n t h e r e i s a (m-1) - p a r a m e t e r f a m i l y of p o i n t s z^, (1 < i < m ) , f o r w h i c h the i m p r o v e d o r d e r of c o n s i s t e n c y i s obtained. The p a r a m e t e r s i n q u e s t i o n a r e the points t ^ i n the d e f i n i t i o n of w (x) i n t h e o r e m (3. 1). P a r t i c u l a r e x a m p l e s of these s p e c i a l c o l l o c a t i o n p o ints w i l l be i n c l u d e d i n the next s e c t i o n . 4. E X A M P L E S O F F I N I T E D I F F E R E N C E A P P R O X I M A T I O N S E x a m p l e 4. 1. L e t n = 1, r = 0 and s = 1; i . e . , we want to c o n s t r u c t a two p o i n t f i n i t e d i f f e r e n c e a p p r o x i m a t i o n to the equation (4.1) L y ( x ) = y'(x) + a°(x)y(x) = f(x) Take m = 1 and zi = x^+± = _-(x.. + X j + 1 ) . Then the b a s i s f u n c t i o n s w 1(x) as def i n e d by (2.13) and (2.14) a r e w°(x) = -(x-x. .) / h., and w X ( x ) = (x-x.)/h., J+1 J J ' J The d i f f e r e n c e a p p r o x i m a t i o n has the f o r m •22-where the coef f i c i e n t s d Q , and e^ can be obtained f r o m (2. 10) and (2. 11), v i z . A T 0 , V "I , 1 0 d Q = Lw ( Z l) = + 2 a. +± , and A T 1 / \ 1 , 1 0 d 1 = Lw ( Z j L) = — + 2 a, , i h. J 1+1 The differ e n c e approximation can then be written as (4.2) { uj +r uj ) / hj + i aj+i (V uj+i ) = • which i s r e c o g n i z e d as the centered E u l e r equation, or Box scheme of K e l l e r (1969). A c c o r d i n g to t h e o r e m (2.1) the o r d e r of consistency i s at l e a s t equal to one. That the o r d e r of c o n s i s t e n c y actually equals two follows f r o m theorem (3.1), because z^ is the c r i t i c a l point of w 2(x) = ( x - X j ) ( x - x . . + 1 ) . A n equally s i m p l e diff e r e n c e equation, the t r a p e z o i d a l method, i s d e r i v e d when taking m = 2, with z, = x. and z_ = x. , T h i s equation has 1 j 2 j+1 the f o r m d_.u. + d . u. , . = e . f. + e 0 f . ... ' , : 0 j 1 J+1 1 J 2 j+1 2 2 With w (x) = (x - x _ . ) ( x - X j ^ ) / h. the c o e f f i c i e n t s are 1 E 0 2 Lw (z^) L w (z^) L w ° ( z 2 ) L w 2 ( z 2 ) d l = E 1 2 L w (z^) Lw (z^) 1 2 L w ( z 2 ) L w ( z 2 ) -23-1 2 -1 2 2 2 e_ = __" L w (z_.) a n d e2 = ~ET " L i W ^ z i ) ' w i t n E = L Q w ^ z l ^ ~ L 0 W ^ Z 2 ^ A N C I i L^y(x) = y (x). The equation thus obtained is (u., .-u.) / h. + |(a°u.+a°, ,u. .) = _-(f'.+f. ,). J + 1 J J J J J+1 J+1 2 V J j+l' By theorem (2. 1) the order of consistency of this difference equation is equal to two. Next, consider the general two point difference approximation with m = 2. That such a difference approximation is always consistent follows from theorem (2.4). Theorem (2.1) implies that the order of consistency is equal to two. From theorem (3. 1) it follows that the order becomes ~ 3 three if the collocation points coincide with the cr i t i c a l points of w (x) = (x-t, )(x-x.)(x-x. , ,). Let t, = x. + rih. , z. = x. + £,h. and z_=x. + £_h.. 1 J J+1 1 J J 1 J 1 J 2 j 2 j Then the special values of ^ and £ 2 that give a third order formula are i the roots of p (£) = 0, where p(£) = £(£-l)(£-n). These roots are frz— t _ V+ 1 ± v f) -T7 + 1 5 - • • •• 3 If n = •§• then the roots are £ = \ ± ^  <y~3 . In that case the collocation points are placed symmetrically in the interval [ x ^ x ^ ^ ] . This results in another order of consistency being gained, so that the actual order of consistency is then equal to four. The coefficients of this approximation are -j - 1 , 1 , 0 , . , 0, 1 , 0, . 0, . d Q = — + (a ( z ^ + a ( z 2 ) ) - y-r h.a ( z ^ a (z 2), j J d«_ = jj- + 2f ( a (z_) + a ( z 2 ^ + T7 h i a < z l ) a (z2*' -24-e 1 = 2 " h j a ( z 2) a n d e 2 ~ 2 + h j i 2 ~ a ( z l ^ It should be p o i n t e d out here, that the c o l l o c a t i o n points i n t h i s d i f f e r e n c e a p p r o x i m a t i o n a r e the G a u s s i a n points w i t h r e s p e c t to the i n t e r v a l [XJ , X j ^ ] . The p r o c e d u r e of S e c t i o n 2 w i t h m G a u s s i a n p o i n t s w i l l y i e l d a f o r m u l a of o r d e r 2 m. T h i s i s h i g h e r than p r e d i c t e d by t h e o r e m (3. 1). However, t h i s r e s u l t does not g e n e r a l i z e to h i g h e r o r d e r equations. The c o r r e s p o n d i n g d i f f e r e n c e equations a r e e q u i v a l e n t to g e n e r a l i z e d m-stage i m p l i c i t R unge-Kutta p r o c e s s e s . (See B u t c h e r (1964), W r i g h t (1970).) It i s a l s o e q u i v a l e n t to a c o l l o c a t i o n p r o c e d u r e w i t h s p l i n e s of degree m (W e i s s (1974), Hulme( 1971). ) E x a m p l e 4.2. W i t h r = 0, s = 2, m = 1, z , = x. , _ and h. ,,. = h., _ = h the * 1 j + 2 j+1 j+2 p r o c e d u r e of S e c t i o n 2 g i v e s a three point (two step) equation of the f o r m d 0 U j + d l U j + l + d 2 U j + 2 = f ( x j + 2 } ' w h e r e d^ = L w 1 ( z ^ ) , 0 < i < 2 . The d i f f e r e n c e a p p r o x i m a t i o n to equation (4.1) i s then found to be r— {u.- 4u. . + 3u. } + a?. u. = f , 2h L j j+1 j+2 J j+2 j+2 j + 2 w h i c h i s one of Gear's f o r m u l a s , s p e c i a l l y s u i t e d f o r s t i f f equations. (See G e a r (1971), p. 217.) H i g h e r o r d e r Gear's methods a r e obtained w i t h r = 0, s > 2, m = 1 and z . = x. , . The o r d e r of c o n s i s t e n c y i s e q u a l to s. 1 j + s y -i E x a m p l e 4.3. L e t n = 2 and r = s = 1, i . e . , we want to c o n s t r u c t a c o m p a c t -25-finite difference a p p roximation to the equation (4.3) Ly(x) == y"(x) + a 1(x)y'(x) + a°(x)y(x) = f(x). With m = l and z^=x^ the coe f f i c i e n t s i n the dif f e r e n c e equation a r e given by d - l U j - l + d 0 U j + d l u j + l = £ j • d , = L w (x.) d n = Lw (x.) = 0 y 2-h.,, a. ,1+1 .1 h j ( h . + h. + 1) • • 2 ( h r h i+ i ) a i , J J + 1 J and 2 2+h.a. d, = Lw (x.) = z J j 1 y V V V ^ F r o m theorem (2. 1) it follows that the o r d e r of consistency of this d i f f e r e n c e approximation i s equal to one. The o r d e r of con s i s t e n c y ~ 3 becomes two i f one lets c o i n c i d e with the i n f l e c t i o n point of w (x) = (x-x. )^ (x-x.)(x-x. ^ ) . Th i s point i s given by z^ = x. + (h. ^-h.) / 3. The coefficients d^ of the c o r r e s p o n d i n g finite d i f f e r e n c e equation are given by 1 ~ h.(h.+ h. ± () J J J + 1 *0 ~ h.h., . J J+1 •2 + (Zh.+h.^Jfh.+Zh.^) a (z a (z, -26-1 h j + 1 ( h j + h ) 2+ 3 J a (z A)+ 1 — ^ a ( z ^ E x a m p l e 4.4. A t h i r d o r d e r three p o i n t f o r m u l a f o r equation (4. 3) i s obtained w i t h m=3, z, =x. ., z~=x. and z =x.,,. The c o e f f i c i e n t s of the 1 J-1 2 j 3 j+1 d i f f e r e n c e equation a r e somewhat c o m p l i c a t e d i n thi s c a s e . However, i f the d i f f e r e n t i a l equation i s g i v e n by (4.4) L y ( x ) = y"(x) + a°(x)y(x) = f(x) , then the r e s u l t i n g f i n i t e d i f f e r e n c e equation i s i d e n t i c a l to the " M e h r s t e l l e n v e r f a h r e n " m e n t i o n e d i n S e c t i o n 1. ( E q u a t i o n (1. 7). ) T h e r e the m e s h i s a s s u m e d to be u n i f o r m . If t h i s i s not the c a s e then the d i f f e r e n c e equation b e c o m e s u. , ,- u. u.- u. 4 _ J±J 1 _ •) . l - 1 (4:5) — J + w ( A ) = w ( f ) , where W. w. = a.w. , + (l-a.-b.)w. + b.w. ,, , h J J J-1 J J J J J+1 w i t h v . 2 i h i + l a j ~ 6 6h.(h.+ h.__.) J J J+1 and -27-2 , h. b = i _ J j 6 6h. x 1(h.+ h. x.) • J+ 1 J J + 1 F o r a u n i f o r m m e s h the o r d e r of c o n s i s t e n c y i s equal to f o u r . E x a m p l e 4. 5. In thi s example, a number of c h o i c e s of the c o l l o c a t i o n points z. w i l l be given, f o r w h i c h the o r d e r of c o n s i s t e n c y of the c o r r e s p o n d i n g f i n i t e d i f f e r e n c e equations i s g r e a t e r than p r e d i c t e d by t h e o r e m (2.1). A g a i n the d i f f e r e n t i a l equation under c o n s i d e r a t i o n i s g i v e n by (4. 3) and the f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s a r e a s s u m e d to be c o m p a c t w i t h r = s = l. The m e s h i s taken u n i f o r m . F i r s t l e t m = 2. T h i s g i v e s a d i f f e r e n c e a p p r o x i m a t i o n of o r d e r two. L e t z ^  = X j + £^h, z^ = X j + ^ h i and t^ = x.. + r j h . F r o m t h e o r e m (3. 1) i t f o l l o w s that the o r d e r b e c o m e s three i f the v a l u e s of e_ and e_ a r e taken to be the r o o t s of p (£) = 0, w h e r e p(£) = (£-n)(£-l) £(£+1). The c o l l o c a t i o n points w i l l be s y m m e t r i c a l l y p l a c e d i f one l e t s n = 0. T h i s y i e l d s z. = x. - and z_ = x- + . B e c a u s e of s y m m e t r y and the u n i f o r m m e s h the o r d e r of thi s equation a c t u a l l y equals f o u r . Next, c o n s i d e r the c a s e m = 4. T h i s w i l l i n g e n e r a l l e a d to a f o u r t h o r d e r f o r m u l a . The o r d e r b e c o m e s f i v e i f the c o l l o c a t i o n p o ints a r e the i n f l e c t i o n points of w^(x) = (x-t^Hx-t^Hx-t^Mx-x j ^ ) ( x - X j ) ( x - x ^ + ^ ) . Now r e q u i r e the c o l l o c a t i o n p o ints to have the f o r m Z l = X j - | " ^ z 3 = X j + i - eh z 2 = x . _ i + eh . z 4 = x j + 1 + eh -28-This choice has the advantage that the co l loca t ion points of consecutive difference equations par t ly coincide , thereby l i m i t i n g the. .necessary total number of function evaluat ions. Upon some computation i t i s found that these points are defined by the value So wi th either value of £ and wi th z. as above one obtains a finite ' 1 difference approximat ion for which the o rde r of consis tency i s at least equal to five according to theorem (3.1) . The actual order equals s ix because of symmet ry and the un i fo rm m e s h . Example 4. 6. In this example some five point finite difference a p p r o x i -mat ions to (4.3) w i l l be d i scussed . So let r = s = 2. F i r s t take m = l , w i th z^ = x j - The o rder of consis tency of the corresponding equation i s equal to three. Th is follows f rom theorem (2.1). If the mesh is un i fo rm then this equation i s iden t ica l to (2. 5). B y theorem (3. 1) the order i s then equal to four. If m = 5 and i f the co l loca t ion points z.., ( l < i < 5), coincide wi th the meshpoints, then the order of the resu l t ing equation i s equal to seven i f the mesh is not un i form and equal to eight i f the mesh i s un i fo rm. Of course the coefficients are ve ry compl ica ted i n this case . However , i f the different ia l equation does not involve the a (x)y (x) t e rm, i . e . , when the equation i s given by (4.4), then the finite difference approximat ion is the same as the five point "Mehrs te l l enver fahren" of C o l l a t z (I960), p . 502, and can be expressed as or s (31u + 128u where -29 -W , w - = (23w. _ + 688w. , + 2358w. + 688w.. , + 23w.^ o ) /3780 . h J 3-2 J - l J J+1 J + 2 " A s i n example (4.5) one may w i s h to determine the spec ia l points z^, for which the finite difference equation has an ex t ra high order of cons is tency. If m = 2 and i f the mesh is un i fo rm, then a computation s i m i l a r to those i n the previous example shows that the order is equal to s ix i f one lets z4 = x . - £h, z„ = x . + £h, and i f the value of £ i s either 1 J 2 J 5. F I N I T E D I F F E R E N C E A P P R O X I M A T I O N S T O B O U N D A R Y C O N D I T I O N S To define a complete difference scheme, approximations to the boundary conditions are r e q u i r e d . Le t a boundary condit ion be given by k-1 (5.1) By(0) = y ( k ) ( 0 ) + 2 b Y ( l ) ( 0 ) = b, i=0 wi th k < n . To construct a finite difference approximat ion to (5.1) one may proceed i n a fashion that resembles the const ruct ion of difference equations i n Section 2. Le t d = s + m and assume that a po lynomia l p(x) e P ^ is uniquely determined by the equations (5.2) p(x.) = u . , 0 < i < s (5.3) L p ( Z i ) = f(z.) , 1 < i < m -30-T h e n the equation (5.4) Bp(O) = b generates a r e l a t i o n of the f o r m m (5.5) B.u,, = ). d.u. = e„b + ), e.f(z.) h 0 . ^ 1 1 0 .u. l v r i=0 1=1 If the basi s functions w 1(x), ( 0 < i < d ) , of a r e chosen as i n Section 2, (Equations (2.13) and (2.14) with r = 0 . ) , then one can show that the coe f f i c i e n t s can be r e p r e s e n t e d by the following determinants: (5 . 6 ) d i = E; Bw1(0) Bw s + 1(0) . . Bw d(0) JLw 1(zp J_.wS + 1(z^) . . J_.wC1(zj) L w ^ z ) J L w S + 1 ( z ) . . Lw^(z ) m m m 0 < i < s (5.7) '0 E, Lw S + 1 ( Z j ) Lw S + '(z ) m Lw d(z 1) Lw^(z ) m (5.8) e = ^ -1 E 0 <3+ 1 B w S (0) L w (z .) L w S + 1 ( z . ,) l - 1' L w S + 1 ( z i + 1 ) L w S + 1 ( z ) m B w d ( 0 ) L w d ( z )^ L w (z^ ^ L w d ( z i + 1 ) L w d ( z ) m 1 < i < m It i s convenient to define the n o r m a l i z i n g factor E N by (5.9) E 0 = (-i) m+1 T s + 1 / \ L Q w (z{) L Q w (zj) L w S + 1 ( z ) ' ' L A w d ( z ) 0 m 0 m where L Q y ( x ) = y ^ ( x ) . If m = 0 define E Q = 1 . The proofs of the fol lowing theorems c lose ly fol low those given i n Sect ion 2 . T h e o r e m 5. 1 Assume that E Q ^ 0 . If m = 0 let s >k; otherwise i f m > 0 , let s + m > n . Then at least k + 1 of the coefficients d. are non-— l ze ro for a l l s m a l l enough h , and the t runcat ion e r r o r m T q S B h y ( 0 ) - e Q b - I e.f(z.) i = l of the difference approximat ion (5.5) sat isf ies T Q < C h s+m-k+1 -32-i . e . , the order of consis tency of (5.5) i s at least equal to s + m - k + 1 . Proof . Cons ider the equations (5.10) s m J d.vAx.) - YJ e i L w i ( z . ) = e 0 B w i ( x 0 ) , 0<i<d . i=0 i= l Here the polynomials w (x) are the same as those used i n the defini t ion of the coeff ic ients . These equations can a lso be wr i t t en as I - A 0 - C l0 m = e, B w " ( x 0 ) B w d ( x 0 ) where I is the (s+1) X (s+1) identi ty m a t r i x , L w V j . . • L w " ( z ) 1 m L w (z^) L w (z ) m and C = L w S + 1 ( z j ) L w d ( z j) L w S + ' ( z ) m L w d ( z ) rrr -33-Thus the equations (5.10) a d m i t a unique s o l u t i o n o n l y i f det(C) + 0. In p a r t i c u l a r , the a s s u m p t i o n that E Q + 0 i m p l i e s that det(C) + 0, p r o v i d e d that h i s s m a l l enough. Next l e t d^, ( 0 < i < s ) , be d e f i n e d as the c o e f f i c i e n t s d^ i n (5.6), (k) but w i t h L Q and B Q r e p l a c i n g J_/ and B. H e r e BQy(O) = y (0). A s s u m e that l e s s than k+1 of the c o e f f i c i e n t s d^ a r e n o n z e r o . M o r e p r e c i s e l y suppose that d^ #0, (1 < v < k^ < k+1), and d^ = 0 o t h e r w i s e . L e t v k l q(x) = "I ( x - x 0 ) c l _ ( x ) ' w h e r e q^(x) e ^ has l e a d i n g c o e f f i c i e n t l / k ! . v,= l v " 1 T h e n q j x j e P ^ , and by a s s u m p t i o n k < d . T h e r e f o r e s m Y ^ ( x . ) - Y ^LQqCz.) = e 0 B 0 q ( x 0 ) . i=0 i = l The c o n s t r u c t i o n of q(x) m akes the e n t i r e l e f t h a n d side v a n i s h . A l s o , BQq(xQ) = 1. Hence eQ = 0. B u t t h i s i s a c o n t r a d i c t i o n s i n c e eQ = ( - l ) m + 1 . Thus i f E _ j. 0 then at l e a s t k + 1 of the c o e f f i c i e n t s d. a r e n o n z e r o . If h 0 i i s s m a l l the same c a n be s a i d of the c o e f f i c i e n t s d.. T h i s p r o v e s the f i r s t statement of the t h e o r e m . The p r o o f of the second a s s e r t i o n c l o s e l y f o l l o w s that of p a r t of t h e o r e m (2. 1) and w i l l be o m i t t e d . T h e o r e m 5. 2. The n o r m a l i z i n g f a c t o r i s n o n z e r o i f m = 0 or i f s = k. (In the l a t t e r c ase the d i f f e r e n c e equation i s s a i d to be compact. ) P r o o f . A s a l r e a d y stated, one takes E Q = 1 i f m = 0. If s = k and E Q = 0, then i t f o l l o w s f r o m (5.9) that t h e r e a r e c o n s t a n t s c., ( s + l < i < d ) , s u c h s+m 1 ut that the p o l y n o m i a l q(x) = V c.w 1(x) s a t i s f i e s L f tq(z.) = 0, ( l < i < m ) . B 1 = 8+1 L Q ^ X ) e P s + m _ n a n d s + m - n< m b e c a u s e s = k < n . Hence L Q q ( x ) = q W ( x ) = 0. So q(x) = 0. T h i s however i s i m p o s s i b l e , s i n c e not a l l c a r e equal -34-to zero and since the polynomials w1(x) are linearly independent. As was the case for finite difference approximations to the differential equation, it is possible to identify collocation points z^, for which the approximation to the boundary condition attains a higher order of consistency, than predicted by theorem (5.1). A minor modification of the proof of theorem (3.1) shows the following. Theorem 5. 3. Let the coefficients d. and e. and the normalizing factor E Q be given by equations (5.6), (5.7), (5.8) and (5.9). Assume that E . + 0 and that s+m > n. Suppose w S * m + ^ ( x ) e P satisfies 0 _ r-r- s+m+1 ~ s+m+1. . _ _ ^  . ^  w (xi) =0, 0 < I < s , and jk ~ s+m+l/ri, d w (0) _ A k ' dx Also assume that the collocation points are chosen such that LQW S + m +*( z.) =0, (1 < i < m). Then the order of consistency of the finite difference approxi-mation (5.5) is at least equal to s + m - k + 2. 6. E X A M P L E S OF APPROXIMATIONS TO BOUNDARY CONDITIONS Example 6.1. Let n = 2 and k = 1. Thus the differential equation is given by (6.1) Ly(x ) E= y"(x) + a^ xjyV) + a°(x)y(x) = f(x), and the boundary condition by -35-(6. 2) By(0) = y (0) + b 0y(0) = b. L e t m = 0, so that the finite d i f f e r e n c e approximation to (6. 2) does not involve the d i f f e r e n t i a l equation (6.1). With s = l the c o n s t r u c t i o n p r o c e d u r e of the previous s e c t i o n y i e l d s e^ = E Q = 1 , d Q = Bw (0) = B = ~h + V and d 1 = Bw (0) = B (x-x )/h 1 1"" h The approximation i s therefore u r u o — h ~ ~ + b o u o = b ] which has o r d e r of c o n s i s t e n c y equal to one. If s = 2 and i f the m e s h i s uniform, then one obtains an equally well-known 3-point equation, v i z . u 2 - 4 u 1 + 3 u Q 2h + b Q u 0 = b, which is a second o r d e r formula. E x a m p l e 6.2. Now let m = l and s = l , i s compact and t h e r e f o r e .consistent for 1 Z -(x-x.)/h, w (x) = ( x - x n ) / h and w (x) = ( so that the approximation to (6. 2) any choice of z^. L e t w^(x) = x - x 0 ) ( x - X l ) / h . T h en the -36-2 2 n o r m a l i z i n g factor i s given by E ^ = L^w (z^) = — . The co e f f i c i e n t s h of the difference a p proximation ( 6 . 3 ) do uo + d l u l = e o b + e l f ( z l ^ ' a r e then equal to eo = ET L w 2( z I > = 1 + I a l ( z i ) -0 d 0 = E " 0 Bw°(0) Bw 2(0) L w " ( z j ) L w ^ ( z j ) = T - a ( Z l) + b Q + 2-b Qa ( z p , and B w V ) L w (z 1' Bw 2(0) Lw (z^) k + a (z ^  + ^ a ( z ^ The o r d e r of con s i s t e n c y i s equal to two. It i s now of i n t e r e s t to determine for which choice of the c o l l o c a t i o n point z^ the o r d e r i n c r e a s e s to three. F r o m theorem ( 5 . 3 ) i t follows that this happens i f z. i s the x h ~ 3 2 1 1 i n f l e c t i o n point of w (x) = x (x-x^). T h i s point i s z^ = — = — . Example 6 . 3 . L e t s = 2 and m = 1, so that the d i f f e r e n c e a p proximation to ( 6 . 2 ) has the f o r m -37-(6.5) d Q u 0 + dini + d 2 u 2 = e Q b + e ^(z ^  . A s s u m e that the m e s h i s u n i f o r m and l e t z^ = X Q + £h. Then E Q = 3 L Q W ( X Q + |h) and E Q = 0 i f and only i f £, = 1. Hence one should not l e t z ^  = x i » The o r d e r of c o n s i s t e n c y i n c r e a s e s to fou r i f z ^  i s an i n f l e c t i o n p o i n t of w (x) = x ( x - x ^ x - x ^ ) . T h i s happens i f i = (9 - \ / 3 3 ) / l 2 = 0. 27 or % =(9 + N/"33)/12 = 1. 23. 7. T H E S T A B I L I T Y O F F I N I T E D I F F E R E N C E S C H E M E S . In th i s s e c t i o n the s t a b i l i t y t h e o r y of K r e i s s (1972) w i l l be used to i n v e s t i g a t e the s t a b i l i t y of f i n i t e d i f f e r e n c e schemes based upon c o n s i s t e n t f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s i n t r o d u c e d i n p r e v i o u s s e c t i o n s . F o r this p u r p o s e the m e s h i s a s s u m e d to be u n i f o r m , i . e . , h^ = h = l / J , (1 < j < J). . .• . - . c. . . '.. , . .' A d i f f e r e n c e scheme c o n s i s t s of equations of the f o r m s m (7. 1) L, u. = V d. .u. ,. = V e. . f ( z . .) = f., r < j < J-s , h j U J , i J+i u j , i v j , i ' j - J -i=- r i = l together w i t h b o u n d a r y c o n d i t i o n s s k ( 0 ) m k ( 0 ) (7.2) B h ) k ( 0 ) u 0 , £ d k ) . ( 0 ) u . = e k ) 0 ( 0 ) b k ( 0 ) + 2 e f c > . (0 ) f ( z k > . (0)) i=0 i = l = b k ( 0 ) , l < k < n 0 , 0 m k { i ) <7'3) B n t ^ J s . l | n d k , i < 1 ) u J + i = e k , 0 < 1 ) V 1 ) + Z 1 ^ . i ^ k , ! * 1 " i = - r k ( l ) i = l b, (1), 1 < k < n - n r -38-In addition, i f (7. 1) i s not compact, i . e . , i f r + s > n then r + s - n extra boundary conditions a r e r e q u i r e d i n o r d e r to match the number of equations and the number of unknowns. Although this is not necessary, i t w i l l be ass u m e d that these extra equations a r e also consistent with the di f f e r e n t i a l equation and given by s. m. (7.4) L, u. = V d. .u. = X e. .f(z. .) = £., k„ < j < r - l , k n > 0 , v h j u j , i i u J , I X J , I ' y 0 - J - 0 — i=0 i = l 0 m j 1 d j , i U J + i = I e),ii{zj,i] 55 V J " B + I < j < J-n+kQ. i=-r. i = l J T L e t Uj. = (UQ, ' " , U J ) . Then the equations (7. 1) through (7.5) can be compactly written as (7. 5) L. u. = H J (7.6) L h u h = f h . H e r e f^ e R J + X i s the app r o p r i a t e right hand side vector and L ^ i s a ~ J + 1 e a s i l y seen to i m p l y that f.-f(x.) < ch. A l s o , let e^ e R be the (J+1) X (J+1) matrix. Consistency of the equations (7.1), (7.4) and (7. 5) is I . n. ISO c ,  iJ 1 y T error vector, i.e., e^ = (e Q, "',ej) with e^  = y. - u.., and let be the vector of truncation errors, T n = ( T_(0)> "" > T n T k ' " ' TJ-n+k ' T ^ D , • • • , T ( i : ) ) T , with l n-n Q Tk< 0 ) 5 5 B h , k ( 0 ) y 0 " V°>« 1 < k < V T j = L h Y j " V k Q < j < J - n + k Q , -39-and s B h , k < 1 ) y J - V 1 ^ 1 < k < n - n Q F o r w, e R J + 1 l e t I lw, j j = m a x |w. |. If A, i s a (J+1) X (J+1) m a t r i x 0 < j < J J I I I I i l l ! ^ ^ l l ^ l l ^ ^ then I j A, | is the induced matrix norm, i.e., | |A, | | = max . w h*0 ||wh|| The f i n i t e d i f f e r e n c e scheme (7.6) i s s a i d to be c o n s i s t e n t i f | |T^| | < k^h _ i as h - * 0 and stable i f e x i s t s f o r a l l s m a l l enough h and | | T ^ | j < k^. H e r e k^ and k^ a r e constants that do not depend on h. The s t a b i l i t y p r o p e r t y i s e s s e n t i a l l y d e t e r m i n e d by the d i f f e r e n c e a p p r o x i -m a t i o n to the h i g h e s t d e r i v a t i v e . To be m o r e s p e c i f i c , the n o t i o n of c h a r a c t e r i s t i c p o l y n o m i a l of (7. 1) i s needed. L e t E and d^ be as i n Se c t i o n 2, (Equations (2.12) and (2.16)), and define c(t) ^ h n I d . t r + i . i = - r It i s not d i f f i c u l t to show that the equation c(t) = 0 m u s t have a r o o t t = l of m u l t i p l i c i t y n, be c a u s e of c o n s i s t e n c y . Hence one c a n w r i t e c(t) = ( t - l ) n c ( t ) , w h e r e c(t) i s d e f i n e d to be the c h a r a c t e r i s t i c p o l y n o m i a l a s s o c i a t e d w i t h the d i f f e r e n c e equation (7. 1). A s s u m e that c(t) i s e x p l i c i t l y g i v e n by N (7.7) c(t) = V a.t 1, w i t h N = r + s - n. i=0 1 If the f i n i t e d i f f e r e n c e scheme (7. 6) i s not compact, then one a l s o has r + s - n c h a r a c t e r i s t i c p o l y n o m i a l s a s s o c i a t e d w i t h the e x t r a b o u n d a r y c o n d i t i o n s (7.4) and (7.5). It i s a s s u m e d that these have the f o r m -40-N. J (7.8) c.(t) = £ a.^.t 1 , k Q < j < r - 1 , i=0 and (7. 9) c . ( t ) = J a t 1 , J - s + 1 < j < J - n + k„. i=0 F i n a l l y , c o n s i d e r the homogeneous d i f f e r e n c e equation N (7. 10) YJ V j + i = ° ' 0 < j < oo. i=0 w i t h boundary c o n d i t i o n s N. J (7. 11) Y a- i v - = °» kQ ! J < r " x> S U P l v - I < const. , 0 < j < oo , . r . J ' 1 1 0< j < oo J 1=0 —J — and N (7.12) Y a N - i V J - j - i = °' 0 < J < o o > i=0 subject to N. J const. (7.13) Y a i N ^ . = 0, J - s + l < j < J-n+k Q, sup Ivj.Jl i=0 j 0 < j < c o Th e n one c a n state the f o l l o w i n g t h e o r e m due to K r e i s s (1972). T h e o r e m 7.1. L e t the homogeneous p r o b l e m c o r r e s p o n d i n g to (1.1) and (1. l a , b) o n l y a d m i t the t r i v i a l s o l u t i o n . A s s u m e that the d i f f e r e n c e scheme (7.6) i s c o n s i s t e n t and that a l l r o o t s t of the c h a r a c t e r i s t i c equation c(t) = 0 -41-s a t i s f y t. + 1. 1 l F u r t h e r , i f the d i f f e r e n c e scheme i s not c o m p a c t then a l s o a s s u m e that the d i f f e r e n c e equations (7.10) w i t h b oundary c o n d i t i o n s (7.11) and the d i f f e r e n c e equations (7. 12) w i t h b oundary c o n d i t i o n s (7. 13) o n l y have the t r i v i a l s o l u t i o n . T h e n (7.6) has a unique s o l u t i o n f o r a l l s m a l l enough h and the d i f f e r e n c e scheme i s s t a b l e . We now i n v e s t i g a t e the s t a b i l i t y p r o p e r t i e s of the a p p r o x i m a t i o n (7. 1). If equation (7. 1) i s c o m p a c t (and co n s i s t e n t ) then c o n s i s t e n c y of the boundary c o n d i t i o n s (7.2) and (7.3) i s s u f f i c i e n t to guarantee s t a b i l i t y . If (7.1) i s not co m p a c t then the c h a r a c t e r i s t i c p o l y n o m i a l c(t) s a i d to be I"-!" S — II 1 s y m m e t r i c i f c(t) = t c(—) and s t r i c t l y d i a g o n a l l y dominant i f the r+s -n d e g r e e of c(t) i s even and i f | a J g | > _ ^ l a j j - ( H e r e I = (r+s-n) /2. ) To m o t i v a t e the l a s t d e f i n i t i o n , c o n s i d e r the g e n e r a l f o r m of a fi v e p o i n t f o r m u l a that i s c o n s i s t e n t w i t h the second d e r i v a t i v e . T h i s 2 2 2 a p p r o x i m a t i o n c a n be w r i t t e n as K. D. u. w here D, u. = (u. , , - 2u.+ u. 4 ) / h v v h h j h j j+1 j and K,w. = a„w. , + a.w. + a 0w... The c h a r a c t e r i s t i c equation i s h j 0 j-1 1 j 2 j+1 > 1 2 c(t) = a^ + a^t + a^t and c(t) i s d i a g o n a l l y dominant i f and only i f the o p e r a t o r ( m a t r i x ) i s d i a g o n a l l y dominant. N o r m a l l y one w i l l c o n s t r u c t the d i f f e r e n c e a p p r o x i m a t i o n s (7. 1) su c h that c(t) i s s y m m e t r i c . (But the c h a r a c t e r i s t i c equations a s s o c i a t e d w i t h the e x t r a b o u n d a r y c o n d i t i o n (7.4) and (7.5) need not be.) L e m m a 7.2. A s s u m e that c(t) i s s y m m e t r i c . If the degree of c(t) i s odd then c ( - l ) = 0. -42-If the degree of c(t) i s even and i f c(t) i s s t r i c t l y d i a g o n a l l y dominant w i t h p o s i t i v e c o e f f i c i e n t s , then t h e r e a r e no r o o t s of c(t) = 0 on the unit c i r c l e . P r o o f . The p r o o f of the f i r s t statement i s i m m e d i a t e . If the d e g r e e of c(t) equals 2N^, then 2 N N , - l . N fi l c ( e ) | = L a k e = l ( a N + h a k ( e + e } ) e I k=0 v 1 k=0 N f l N -i N A-1 = I a, T + 2 V a. cos k0 I > a. T - 2 I V a, c o s k9 I > a M -2 V a. > 0. ' N l k^O k N l k^O k N l k^O k E x a m p l e 7. 3. In the s p e c i a l c a s e w h e r e c(t) = a^t + ( l - 2 a Q ) t + a^ the a s s u m p t i o n s of the l e m m a h o l d i f 0 < a^ < A s i m p l e c o m p u t a t i o n shows that i n fact t here a r e no r o o t s on the unit c i r c l e i f f - co < a^ < T h i s shows that the a s s u m p t i o n s i n the l e m m a a r e not s t r i c t l y n e c e s s a r y , although p e r h a p s d e s i r a b l e . The effect of the c h o i c e of the c o l l o c a t i o n p o i n t s z^, (1 < i < m),on the r o o t s of the c h a r a c t e r i s t i c equation c(t) = 0 w i l l be i l l u s t r a t e d by means of the next e x a m p l e . 2 E x a m p l e 7.4. L e t c(t) = a^t + (l-2aQ)t + aQ and l e t n = 2, i . e . , we c o n s i d e r s y m m e t r i c 5-point a p p r o x i m a t i o n s to the second d e r i v a t i v e . F o r c(t) to be s y m m e t r i c i t i s n e c e s s a r y to p l a c e the c o l l o c a t i o n points s y m m e t r i c a l l y i n the i n t e r v a l X j - 2 ' Xj+2. |. Hence i f m = 1 t h e r e i s one c o l l o c a t i o n p o i n t z^ = Xy The d i f f e r e n c e a p p r o x i m a t i o n to y (x) = f(x) i s i n t h i s c a s e g i v e n by f - u . 0+16u. .-30u.+ 16u., . - u. ,~^ / 1 2 h 2 = f(x.). The l e f t hand s i d e c a n y \ 3-2 J-1 J J+1 J+2/ ' v y a l s o we w r i t t e n as 1_ 12b 2 so that the c h a r a c t e r i s t i c p o l y n o m i a l of t h i s a p p r o x i m a t i o n i s g i v e n by 1 2 14 1 c(t) = - J2" t + t - Yz ' c h a r a c t e r i s t i c equation has r o o t s L = 7 ± 4 \/T, so that | t | ,£ 1. Next, i f m = 2 l e t z, = x. - qh and z_ = x. + qh, q > 0. T h i s generates 1 J ^ J a d i f f e r e n c e a p p r o x i m a t i o n to y " ( x ) = f(x) of the f o r m ( a 0 u . _ 2 + ( l ^ a ^ u . ^ + (-2 + 6 a 0 ) u . + ( l - 4 a 0 ) u . + 1 + a Q u . + 2 ) / h 2 = ±(f( z J + f ( z 2 ) ) . 2 The c o r r e s p o n d i n g c h a r a c t e r i s t i c p o l y n o m i a l i s c(t) = a^t + (l-2aQ)t + aQ. 2 It i s e a s i l y shown that a^ and q a r e r e l a t e d by a Q = (6q - 1)/ 12. B y e x a m p l e (7. 3) t h e r e a r e no r o o t s of c(t) = 0 on the unit c i r c l e i f and only i f -co < aQ < 4- I* 1 t e r m s of q t h i s c o n d i t i o n b e c o mes 0 < q < ~ 0. 816. r e If the c o l l o c a t i o n p o ints a r e g i v e n by z^ = x^-q^h and z 2 = x^ + qQh, whe q_ = J1- \j-rp ~ 0. 38 or q„ = Ji + N/T? « 1. 36, then the o r d e r of c o n s i s t e n c y 0 15 U 15 of the c o r r e s p o n d i n g f i n i t e d i f f e r e n c e a p p r o x i m a t i o n i s e q u a l to s i x . (See e x ample (4.6).) T h e r e f o r e the t h e o r y of K r e i s s guarantees s t a b i l i t y only f o r the f i r s t c a se, w h e r e qg = sj 1 - ^Y^ • Of c o u r s e , t h i s does not i m p l y that a f i n i t e d i f f e r e n c e a p p r o x i m a t i o n b a s e d upon the second v a l u e of q^ i s n e c e s s a r i l y u nstable. That s u c h an a p p r o x i m a t i o n may l e a d to a sta b l e s cheme i s s u p p o r t e d by the data o b t a i n e d i n the next s e c t i o n . (See e x ample (8. 1), e x p e r i m e n t number 15.) The f i r s t v a l u e of q^ a p p e a r s to give a b e t t e r e r r o r constant, however. F i n a l l y , take m = 5 and l e t the c o l l o c a t i o n p o i n t s z^, (1 < i < 5), c o i n c i d e w i t h the m e s h p o i n t s xj +j» ( - 2 < i < 2 ) . T h i s generates the f i v e p o i n t M e h r s t e l l e n v e r f a h r e n of C o l l a t z (I960), page 502, v i z . - 4 4 -1 + 128u 3 1 8 u . + 1 2 8 u . , t + 3 1 u J J+1 2 52h' J - 1 1 3f. , + 688 f . . + 2 3 5 8 f . + 6 8 8 f . , . + 23 f . _ V J " 2 J " 1 J J+1 J + 2 / 3780 T h e a s s o c i a t e d c h a r a c t e r i s t i c e q u a t i o n i s 1 ( 3 1 t 2 + 190t + 31) = 0 , 2 5 2 w i t h r o o t s t^ ~ - 0 . 1 6 8 a n d t^ ~ - 5 . 9 6 . H e n c e t h e s e c o l l o c a t i o n p o i n t s l e a d to a s t a b l e d i f f e r e n c e a p p r o x i m a t i o n . N o w c o n s i d e r t h e e x t r a b o u n d a r y c o n d i t i o n s ( 7 . 4 ) a n d ( 7 . 5 ) . A g a i n , a s s u m e tha t t h e c h a r a c t e r i s t i c p o l y n o m i a l c(t ) o f (7 . 1) i s s y m m e t r i c a n d t h a t the d e g r e e o f c( t ) i s e v e n a n d e q u a l to 2 N ^ . If i n a d d i t i o n c(t ) i s s t r i c t l y d i a g o n a l l y d o m i n a n t w i t h p o s i t i v e c o e f f i c i e n t s t h e n t h e c h a r a c t e r -i s t i c e q u a t i o n c(t ) = 0 h a s e x a c t l y N ^ r o o t s i n s i d e t h e u n i t c i r c l e a n d N ^ r o o t s o u t s i d e the u n i t c i r c l e . A n e c e s s a r y c o n d i t i o n f o r the d i f f e r e n c e e q u a t i o n s ( 7 . 1 0 ) , ( 7 . 1 1 ) a n d ( 7 . 1 2 ) , ( 7 . 1 3 ) to a d m i t the z e r o s o l u t i o n o n l y i s t h e n tha t the n u m b e r o f e x t r a b o u n d a r y c o n d i t i o n s a t x = 0 i s the s a m e a s the n u m b e r a t x = 1 a n d e q u a l to N ^ = (r + s - n) / 2 . ( H e n c e = ( r - s + n ) / 2 i n ( 7 . 4 ) a n d ( 7 . 5 ) ) . It i s a l s o r e a s o n a b l e to a s s u m e n o w t h a t t h e c h a r a c t e r i s t i c e q u a t i o n s o f t h e e x t r a e q u a t i o n s a t x = 0 a n d x = 1 a r e r e l a t e d b y (7 . 14) ° J + r - s - j (t) = t J c (-), ( r + s - n ) / 2 = k Q < j < r - l , i . e . , the c o n d i t i o n s a t x = 1 a r e the " r e f l e c t i o n " o f t h o s e a t x = 0 . F o r s t a b i l i t y i t i s t h e n s u f f i c i e n t to s h o w t h a t the d i f f e r e n c e e q u a t i o n s (7 . 10) -45-subject to (7.11) have the z e r o s o l u t i o n only. F o r thi s p u r p o s e define H e r e N = r + s - n and the c o e f f i c i e n t s a^ a r e the same as those of the c h a r a c t e r i s t i c p o l y n o m i a l c(t) i n (7.6). If (7.10) subj e c t to (7.11) a d m i t s a n o n t r i v i a l s o l u t i o n then i t i s e a s i l y seen that the p o l y n o m i a l s c^(t), ((r+s-n) / 2 = kQ < j < r - l ) , and p^(t), (o < j < max N.-N), a r e l i n e a r l y dependent. Hence we have shown the f o l l o w i n g T h e o r e m (7. 5). L e t the homogeneous p r o b l e m c o r r e s p o n d i n g to (1.1) and ( l . l a , b ) only have the t r i v i a l s o l u t i o n and l e t the d i f f e r e n c e s c heme (7.6) be c o n s i s t e n t . A s s u m e that c(t) i s s y m m e t r i c . A l s o suppose that the degree of c(t) i s even and that c(t) i s s t r i c t l y d i a g o n a l l y dominant w i t h p o s i t i v e c o e f f i c i e n t s . L e t the c h a r a c t e r i s t i c p o l y n o m i a l s of the e x t r a b o u n d a r y c o n d i t i o n s be r e l a t e d as i n E q u a t i o n (7.14). If the p o l y n o m i a l s c..(t) and Pj(t) d e f i n e d by (7. 8) and (7. 15) r e s p e c t i v e l y a r e l i n e a r l y i n d e -pendent then the d i f f e r e n c e scheme i s s t a b l e . Hence th e r e e x i s t s a con-stant K independent of h su c h that | | e^ | | < K| | T ^ | |. 2 1 E x a m p l e 7.6. L e t c(t) = a^t + (l-2aQ)t + a^ w i t h - co < a^ < If the degree of the c h a r a c t e r i s t i c p o l y n o m i a l c^(t) of the e x t r a boundary c o n d i t i o n at x = 0 i s a l s o e q u a l to two then the d i f f e r e n c e scheme i s s t a b l e i f c.^(t) ^  c ( t). If the d e g r e e of c^(t) i s three then s t a b i l i t y i s guaranteed 3 2 i f t h e r e i s no constant b s u c h that c , ( t ) = b ( a n t + ( l - 2 a n ) t +a„t)+(l-b) p o l y n o m i a l s p-(t) by N (7. 15) p (t) EE I a t ^ 1 , j > 0 i=0 8. N U M E R I C A L E X A M P L E S The m a i n p u r p o s e of the n u m e r i c a l e x a m p l e s g i v e n i n t h i s s e c t i o n i s to check the c o r r e c t n e s s of statements i n p r e v i o u s s e c t i o n s . They a l s o give some i n d i c a t i o n as to what the r e l a t i v e a c c u r a c y of v a r i o u s d i s c r e t i z a t i o n s i s . A l l c o m p u t a t i o n s w e r e c a r r i e d out on an I B M 370/168, u s i n g double p r e c i s i o n a r i t h m e t i c . No attempt was made to o p t i m i z e the e f f i c i e n c y of the com p u t a t i o n s , so that t h e r e w i l l be no c o n c l u s i o n s about the r e l a t i v e m e r i t of v a r i o u s f i n i t e d i f f e r e n c e s chemes. E x a m p l e 8.1. L e t the d i f f e r e n t i a l equation be g i v e n by Ly ( x ) = y " + y' - 2y = 2 ( l - 6 x ) e X w i t h boundary c o n d i t i o n s y(0) = y ( l ) = 0 , The s o l u t i o n of thi s p r o b l e m i s y(x) = 2 x ( l - x ) e X . The r e s u l t s of some n u m e r i c a l c o m p u t a t i o n s a r e g i v e n i n T a b l e (8.1). In th i s table r, s and m a r e as i n the f i n i t e d i f f e r e n c e a p p r o x i m a t i o n (2.4). So the w i d t h of the a p p r o x i m a t i o n equals to = r + s + 1 and m denotes the number of c o l l o c a t i o n p o i n t s . The l e t t e r s A, B, C, D and E i n the c o l u m n s headed b y " c " a r e a code i n d i c a t i n g the l o c a t i o n of the c o l l o c a t i o n p o i n t s , v i z . A: The points z^, ( l < i < m ) , a r e o p t i m a l i . e . , the o r d e r of c o n -s i s t e n c y i s as h i g h as p o s s i b l e f o r the p a r t i c u l a r v a l u e s of r, s and m c o n s i d e r e d . B: The points z^, (1 < i < m ) , a r e o p t i m a l under the r e s t r i c t i o n s that e a ch s u b i n t e r v a l c o n t a i n the same number of c o l l o c a t i o n p o i n t s , and that these points a r e p l a c e d s y m m e t r i c a l l y w i t h r e s p e c t to the sub-i n t e r v a l . (See the l a s t c a s e d i s c u s s e d i n E x a m p l e 5 of S e c t i o n 4. ) C: The c o l l o c a t i o n p o ints c o i n c i d e w i t h the m e s h p o i n t s and m = r + s + 1 . D: m = 1 and z, = x.. 1 J E: The p l a c e m e n t of the c o l l o c a t i o n p o i n t s i s not o p t i m a l , but they a r e p l a c e d s y m m e t r i c a l l y i n the i n t e r v a l x. , x., The c o l u m n headed by "o" g i v e s the o r d e r of c o n s i s t e n c y of the f i n i t e d i f f e r e n c e a p p r o x i m a t i o n as p r e d i c t e d by t h e o r e m s i n t h i s c h a p t e r . C o l u m n s 7 through 11 define the f i n i t e d i f f e r e n c e equations (2.4) f o r r <_ j < J - s . If the w i d t h of t h i s a p p r o x i m a t i o n i s equal to 5 then a s p e c i a l f i n i t e d i f f e r e n c e equation m u s t be d e f i n e d f o r j = 1. T h i s i s done i n c o l u m n s 2-6. (The s p e c i a l equation n e c e s s a r y f o r j = J - l i s a s s u m e d to be the " r e f l e c t i o n " of the one f o r j = 1.) The m e s h i s taken u n i f o r m i n t h i s example, so that h. = h = 4. F o r a number of v a l u e s of J the J J o b s e r v e d m a x i m u m e r r o r , i . e . , max |u.-y ( X j ) |, i s given. The n o t a t i o n - 1 -1 0.438 means 0.438 10 etc. In the f i n a l c o l u m n , headed by " a " , the o b s e r v e d a s y m p t o t i c o r d e r of a c c u r a c y of the g i v e n f i n i t e d i f f e r e n c e scheme i s l i s t e d . If i t i s not c l e a r f r o m the n u m e r i c a l r e s u l t s what thi s o r d e r i s then the expected o r d e r i s g i v e n between b r a c k e t s . M o s t of the r e s u l t s that appear i n the table a r e s e l f e x p l a n a t o r y . The f i r s t s e ven e x p e r i m e n t s i n v o l v e c o m p a c t d i f f e r e n c e a p p r o x i m a t i o n s . F o r e x p e r i m e n t s 5 and 6 the c o l l o c a t i o n p o i n t s a r e g i v e n b y z. = x. • i_ ± £h. T h e r e a r e two v a l u e s of £ f o r w h i c h the o p t i m a l o r d e r of c o n s i s t e n c y i s r e a c h e d . (See the l a s t c a s e d i s c u s s e d i n E x a m p l e (4. 5)). T h e s e v a l u e s a r e £ = N Y2 ~ 2 45~ ' u s e a - ^ n E x p e r i m e n t 5, and £ = V - j - ^ " + 2" 45~ ' u s e d i n 6 E x p e r i m e n t s 8-12 show the effect that v a r i o u s c h o i c e s of the e x t r a boundary conditions have on the o v e r a l l a c c u r a c y . Note that even i f the o r d e r of c o n s i s t e n c y of the e x t r a f i n i t e d i f f e r e n c e equations i s only equal to two, then t h e o r d e r of a c c u r a c y of the scheme r e m a i n s f o u r . T h i s phenomenon i s a l s o e x p l a i n e d i n the p a p e r of K r e i s s (1972). (See a l s o B r a m b l e and H u b b a r d (1964) and S h o o s m i t h (1975). ) The a c t u a l a c c u r a c y h o wever i s s e r i o u s l y e f f e c t e d . ( A m a z i n g l y th i s i s not the c a s e f o r e x p e r i m e n t 9. but th i s m u s t be c o n s i d e r e d e x c e p t i o n a l . ) In 13 and 14 the c o l l o c a t i o n p o i n t s of the m a i n f i n i t e d i f f e r e n c e equations a r e z ^  = x.. - £h and = X j + £h. A g a i n , as has b e e n m e n t i o n e d p r e v i o u s l y i n E x a m p l e (4. 6), t h e r e a r e two v a l u e s of £ f o r w h i c h the o r d e r of c o n s i s t e n c y becomes s i x . T h e s e a r e £ = v 1 - v T T used i n 13 and 1 e = 1 + N TT used i n 14. E x a m p l e 8. 2. C o n s i d e r the equation 11 1 x y + x y - (1 + x)y = -(2 + x)e w i t h boundary c o n d i t i o n s y'(0) = y ( l ) = 0 -49-T A B L E 8. 1 # r s m c o J = 4 J = 8 J = 16 J .= 32 J = 64 a 1 1 1 2 D 2 .438" 1 . 1 1 1 " 1 . 2 8 1 " 2 .704" 3 . 1 7 6 " 3 2 2 1 1 2 B 2 . 1 1 1 " 1 . 2 7 9 " 2 . 7 0 0 " 3 . 1 7 5 " 3 ' -4 .439 2 3 1 1 2 A 4 -4 . 830 -5 . 534 ° . 3 3 7 ~6 . 2 1 1 " 7 . 1 3 2 " 8 4 4 1 1 3 C 4 .516" 3 -4 . 325 -5 '. 203 . 128"6 .798" 8 4 ' 5 1 1 4 B 6 . 2 1 3 ~ 6 . 3 3 7 - 8 . 5 2 8 " 1 0 - 12 '.817 C * * 6 6 1 1 4 B 6 .704" 6 . 106" 7 . 1 6 5 " 9 : 257" 1 1 6 7 1 1 5 E 6 . 156" 6 . 2 4 9 " 8 . 3 9 0 " 1 0 -12 .628 L 6 8 1 1 1 D 2 2 2 1 D 4 . 2 9 1 " 1 . 2 6 5 ~ 2 . 2 0 4 " 3 . 1 4 7 " 4 . 9 8 7 " 6 4 9 1 3 1 D 3 2 2 1 D 4 . 4 4 0 " 2 -3 . 233 -5 . 897 . 3 1 0 "6 . 1 6 0 " 7 (4) 10 1 4 1 D 4 2 2 1 D 4 -4 . 369 . 4 5 8 "5 .326" 6 _7 . 207 4 11 1 1 2 A 4 2 2 1 D 4 -3 . 533 _4 . 719 _5 .518 . 3 3 6 " 6 _ 7 . 212 4 12 1 1 3 C 4 2 2 1 D 4 .884" 3 -4 . 781 . 5 2 8 ~5 . 3 3 8 " 6 _7 . 212 4 13 1 1 5 E 6 2 2 2 A 6 _7 .571 ' . 1 8 5 " 8 . 3 6 7 " 1 0 - 12 .601 L * * 6 14 1 1 5 E 6 2 2 2 A 6 -4 . 260 . 7 4 5 "6 -7 . 122 . 1 6 0 "9 . 2 1 5 " 1 1 (6) C o n t a m i n a t e d by roundoff -50-The solution to this problem is y(x) = (l-x)e . In this example only compact approximations to the differential equations are considered. Numerical test calculations are performed with i various difference approximations to the boundary condition y (0) = 0. Results appear in Table (8. 2). The notation used is the same as in the previous example. The approximation to the boundary condition is defined in columns 2-6, while the finite difference approximation to the differential equation is defined in the next five columns. In Experiment 9 the collocation point for the boundary condition is z^ = X g + y. For this value of z^ the order of consistency is three. (See Example (6.2).) In Experiments 13 and 14 this collocation point is z^ = X Q + £,h. In 13 the value of £ is (9 - N/ 33)/12 and in 14 this value is (9 + N/33)/12. For these collocation points the order of consistency is equal to four rather than three. (See Example (6. 3). ) Note that the order of accuracy is not greater than the order of consistency of the discrete boundary condition. This differs from observations made about the extra boundary conditions in Example (8.1). Example 8.3. The purpose of this example is to illustrate the usefulness of allowing nonuniform meshes. Let the differential equation be given by 8 2 y " - x 2y= 2 s 2 - x 4 - e e " X / 2 s , £ = 0.001, with boundary conditions -1 y'(0) = 0 and y ( l ) = 1 + e 2 S « 1. -51-T A B L E 8. 2 # r s mc o r s" mc o ' J = '4 J = 8 ' J =16 J = 32 J = 64 a 1, 0 1 0 - 1 1 1 1 D 2 . 117° .494-1 .226" 1 . 108 - 1 .526" 2 1 2 0 2 0 - 2 1 1 1 D 2 .240" 1 .492" 2 . 110" 2 .261" 3 _4 .633 2 3 0 1 1 * 2 1 1 1 D 2 -1 . 250 .631" 2 . 158" 2 .396" 3 -4 .989 . •2 4 0 1 1 **2 1 1 1 D 2 .240" 1 .492" 2 . n o - 2 .261" 3 _4 .633 2 5 0 1 0 - 1 1 1 2 A 4 .106° .466" 1 .219" 1 . 106" 1 .521" 2 1 6 0 2 0 - 2 1 1 2 A 4 .362" 1 .786" 2 . 183" 2 - 3 .443 . 109" 3 2 7 0 1 1 * 2 1 1 2 A 4 .152- 1 .359" 2 .876" 3 .216" 3 -4 .538 2 8 0 1 1 **2 1 1 2 A 4 .338" 1 .764" 2 . 181" 2 .440" 3 . 108" 3 2 9 0 1 1 A 3'; 1. 1 2; A 4. .531"3' -4 . 60.2 * - 5-. 71'! .862"6 . 106" 6 3 10 0 1 2 C 3 1 1 2 A 4 . 146" 2 - 3 . 171 -4 . 207 .254"5 .316" 6 3 11 0 4 0 - 4 1 1 2 A 4 .344" 2 . 164" 3 .898" 5 .526" 6 _7 .317 4 12 0 1 3 E 4 1 1 2 A 4 .482 -5 . 324 0 .208"6 -7 . 132 .828 - 9 4 13 0 2 1 A4 1 1 2 A 4 . 165" 4 . 108" 5 _7 . 730 .472" 8 .300" 9 4 14 (0 2. 1 A^4 I 1 1 2 A 4 .563" 3 _4 . 249 -5 . 129 -7 .734 .437" 8 4 z l z l -52-A s was the case i n the p r e v i o u s examples, this p r o b l e m has been constructed so i t would have a known solution. H e r e the solution i s . . - x 2 / 2 c 2 y(x) = e ' + x , which has a boundary l a y e r at x = 0. One could equally w e l l impose the boundary condition y(0) = 1, but f r o m n u m e r i c a l point of view the condition i y (0) = 0 makes it c o n s i d e r a b l y m o r e d i f f i c u l t to obtain a reasonably a c c u r a t e solution. Hence the given condition i s the m o re i n t e r e s t i n g one here. R e s u l t s a r e given i n Table (8. 3),. In the nonuniform case the m e s h i s taken such that 3 1 there are J meshpoints, equally spaced, to the left of x = 0. 1 and J mesh-points, also u n i f o r m l y spaced, to the r i g h t of x = 0. 1. Hence for each e x p e r i -ment there i s only one finite d i f f e r e n c e approximation that needs to take into account the nonuniformity of the mesh. The actual finite d i f f e r e n c e equations i n v o l v e d a r e the second case d i s c u s s e d i n E x a m p l e (4. 3), (with a^(x) = 0), and the equation given i n E x a m p l e (4.4). T A B L E 8.3 U n i f o r m M e s h # r s m c o r s m c o J=8 J=16 J=32 J=64 I 1 2 0 2 0 - 2 0 4 0 - 4 1 1 1 D 2 1 1 3 C 4 1.000° . 987° . 954° . 883° .486° •909" 1 .923" 1 _ i . 782 Nonuniform M e s h # r s m c o r s m c o J=8 J=16 J=32 J=64 3 4 0 2 0 - 2 0 4 0 - 4 1 1 1 A 2 1 1 3 C 3* .110° . 934" 1 . 127" 1 .718" 2 .508" 2 _ 2 . 223 .203" 2 -3 . 725 3 4 i f the m e s h i s l o c a l l y u n i f o r m E x a m p l e 8.4. L e t the p r o b l e m be the same as i n E x a m p l e (8.2). A p p r o x i m a t e the d i f f e r e n t i a l equation b y a c o m p a c t f i n i t e d i f f e r e n c e equation w i t h r = s = m = l . L e t = x j - T n e m e s h i s taken u n i f o r m . A i one p a r a m e t e r set of a p p r o x i m a t i o n s to the b o u n d a r y c o n d i t i o n y (0) = 0 i s g i v e n by u 0 " a U l + ( a " 1 ) u 2 (a - 2)h = 0 . T h i s a p p r o x i m a t i o n i s c o n s i s t e n t i f a + 2 and the o r d e r of c o n s i s t e n c y i s 4 one. If a - then the o r d e r becomes two and the f i n i t e d i f f e r e n c e equation i s i n that case i d e n t i c a l to the one obtained w i t h the p r o c e d u r e of S e c t i o n 5, l e t t i n g r = 0, s = 2 and m = 0. It w i l l be shown i n S e c t i o n 6 of C h a p t e r II that i f a > 2 then the f i n i t e d i f f e r e n c e o p e r a t o r L ^ has an eig e n v a l u e \ •— + oo as h —*• 0. T h i s m a k e s the a p p r o x i m a t i o n u s e l e s s f o r t i m e dependent p r o b l e m s . F r o m the g e n e r a l t h e o r y of K r e i s s (1972) i t f o l l o w s that the same a p p r o x i m a t i o n when used i n the c u r r e n t b o u n d a r y v a l u e p r o b l e m , does not l e a d to an i n s t a b i l i t y . T h i s i s supported b y the f o l l o w i n g data, obtained w i t h a = 4. T A B L E 8.4 J = 4 J = 8 J = 16 J = 32 J = 64 a .784° : .248° .100° . 4 5 3 " 1 .216" 1 1 -54-Chapter II INITIAL V A L U E PROBLEMS 1. INTRODUCTION This chapter is concerned with the construction of finite difference approximations to initial value problems. In Section 2 the construction procedure of the previous chapter will be extended to systems of linear first order ordinary differential equations subject to initial conditions. Examples of this very general class of difference approximations appear in Section 3. In particular, a constant coefficient problem is discussed in detail, (Example (3.4)), in view of its application to the stability analysis of parabolic equations, (Section 6). Well-known concepts in the stability analysis of finite difference approximations to initial value problems in ordinary differential equations are summarized in Section 4. Particular attention is paid to methods for "stiff" problems, because of their applicability to the numerical solution of parabolic equations. In Section 5 it is shown how the difference methods for boundary value problems and initial value problems may be combined to yield quite general difference schemes for the solution of linear parabolic equations in one space variable. The numerical stability of such schemes is investigated in Section 6. Particular attention is paid to the effect of finite difference approximations to the boundary conditions. Several examples are given. Finally, in Section 7 the results of some numerical experiments are given. -55-2 . I N I T I A L V A L U E P R O B L E M S I N 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 In t h i s s e c t i o n the c o n s t r u c t i o n of f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s to the f i r s t o r d e r l i n e a r s y s t e m ( 2 . 1 ) L w(t) = A ( t ) w ' ( t ) - B(t)w(t) = f(t), 0 < t < T , w i l l be d i s c u s s e d . H e r e w(t) = ^w^(t), w 2 ( t ) , ,Wj _ ^ ( t ) ^ T f(t) = (f^t), f 2 ( t ) , "", fj_l))T and A(t) and B(t) a r e continuous (J-1) X (J-1) m a t r i c e s . It i s a s s u m e d that the above s y s t e m , together w i t h the i n i t i a l c o n d i t i o n w (0 ) = w ° , a d m i t s a unique s o l u t i o n w(t), (0 < t < T), that i s as many t i m e s c o n t i n u o u s l y d i f f e r e n t i a b l e as n e c e s s a r y . The c o n s t r u c t i o n p r o c e d u r e , that i s an e x t e n s i o n of the p r o c e d u r e d i s c u s s e d i n the p r e v i o u s c h a p t e r , a p p l i e s e q u a l l y w e l l to h i g h e r o r d e r s y s t e m s , s u b j e c t to m o r e g e n e r a l boundary c o n d i t i o n s . However, these w i l l not be c o n s i d e r e d h e r e . D e f i n e a m e s h S^ by S h = { t . : 0 = t 0 <t' r< < t N = T } . The m e s h s p a c i n g s a r e ( A t ) = t - t ., ( l < n < N ) . L e t A t = max(At) , c °— n. n n-1 — — n n so that (At) = At = T/N f o r u n i f o r m m e shes. A s s u m e that m in(At) > ^ A t , n ' n — K fo r some constant K ;> 1 . F i n i t e d i f f e r e n c e a p p r o x i m a t i o n s to ( 2 . 1) a r e a s s u m e d to have the f o r m 0 u, ( 2 . 2 ) L . u n = V D .un+£ = 7, E , f ( r ), ,,l<n<N, v ' h Li n, I U n, SL y b n , I" ' — — l = -P 1=1 n •56 with initial condition = w ° . Here u' n , n n (u^, u^, u n j ,^ is the approximation to w ( t n ) - For each n the points t, ., (1 <i < u ), are distinct points, that will normally lie in the interval n-P n n The coefficients D . and E „ are (J-l) n, I n, I x ' X (J-l) matrices. These coefficients may be determined by requiring that 0 p. !=-P .P+P-" 1=1 ,J-1 for a l l p(t)e -^ Here is the space of all (J-l)-vector valued polynomials of degree at most d. (For notational convenience the sub-script n will frequently be omitted.) This requirement leads to the matrix equation (2.3) I 0 0 0 0 I 0 0 0 0 I 0 -L -L -L -L T 0,1 T 1,1 T P, 1 T P+l, 1 0 0 0 -L P+H-l, 1 •L •L. -L •L •L 0, \s. 1, [L T P, [i. T P+l, u. P+U--1, j x J T D-P 0 T D - P + l 0 T Do — 0 T E l 0 E T _ M- _ 0 where I and 0 are the (J-l) X (J-l) identity and zero matrix respectively. T L. . is the transpose of L. ., where L i , i = L 0 w i < ^ ) A < V - " i a i ) B ( ^ ) , -57-w i t h L,^(t) = j^f^  • T h e s c a l a r p o l y n o m i a l s c / ( t ) , (0 < i < P + u . -1 ) , a r e l i n e a r l y i n d e p e n d e n t a n d c h o s e n s u c h t h a t w l ( tn-P+i ) = 6U> °<i> l< P , a n d o J 1 ( t n _ p + ^ ) = 0 , P + l < i < P + L x - l , 0 < 4 < I F o r e x a m p l e , t h e s e p o l y n o m i a l s m a y b e d e f i n e d i n a s i m i l a r f a s h i o n a s t h e p o l y n o m i a l s w X ( x ) i n S e c t i o n ( I . 2). ( E q u a t i o n s (1.2.1,3) a n d ( 1 . 2 . 1 4 ) . ) T h e c o e f f i c i e n t m a t r i c e s a n d E ^ i n ( 2 . 3 ) a r e u n i q u e l y d e f i n e d up to a c e r t a i n n o r m a l i z a t i o n , p r o v i d e d t h a t t h e m a t r i x C EE P + l , 1) Li. 'P+H-1,1 L P + l , (Ji Li, ' P + f J L - 1 , (J. h a s r a n k e q u a l to ( u - 1 )(J - 1 ) . T h i s i s t h e c a s e f o r a l l s m a l l e n o u g h A t , p r o v i d e d t h a t t he m a t r i x L g o / ^ a ^ A c ^ ) L 0co P +^- 1(C 1)A(^) Lo«p+1<yA(V ' ' ' L o " p + f X _ 1 ( y A ( y h a s r a n k (JJL-1) ( J - 1 ) . T h i s i n t u r n i s t r u e i f S + 0 f o r a l l s m a l l e n o u g h A t . H e r e £ i s t h e a n a l o g u e o f the n o r m a l i z i n g f a c t o r E i n the p r e v i o u s c h a p t e r , ( E q u a t i o n ( 1 . 2 . 1 2 ) ) , a n d c a n b e e x p r e s s e d a s (2.4) e = (-D f -58-1 L Q C / - 1 " 1 ^ ) • • • ^ t f ^ i i j ) P+l, A l t h o u g h i t i s not i m p o r t a n t f o r t h e o r e t i c a l p u r p o s e s i n what m a n n e r the d i f f e r e n c e equation (2.2) i s n o r m a l i z e d , i t may be a s s u m e d , that t h i s i s done i n such a way that c < II E E J I < c > i=0 f o r a l l s m a l l A t . H e r e c and C a r e p o s i t i v e c o n s t a n t s that do not depend on A t . The p r o o f s of the f o l l o w i n g statements a r e then v e r y m u c h l i k e those of T h e o r e m (I. 2.1), T h e o r em( I. 2 . 4), and T h e o r e m (I. 3.1) and w i l l be o m i t t e d . T h e o r e m 2.1. A s s u m e that P > 1 and that £ + 0. Then the o r d e r of c o n s i s t e n c y of the f i n i t e d i f f e r e n c e a p p r o x i m a t i o n (2.2) i s at l e a s t equal to P + p-1. In p a r t i c u l a r , i f p. = 1 or i f P = 1 then 8^0. If the p o i n t s a r e ~ p + H i l l 0_ 1 the c r i t i c a l points of co ^(t) = | | (t-y.) | ( t - t .) then the o r d e r i s at 1=1 i = - p n + l l e a s t p+u.. (Thus t h e r e i s a (u.-l)-parameter f a m i l y of s u c h p o i n t s w i t h p a r a m e t e r s \^ , (1 < i < p.-l). ) 3. E X A M P L E S O F F I N I T E D I F F E R E N C E A P P R O X I M A T I O N S In t h i s s e c t i o n a number of s p e c i f i c e x a m p l e s of f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s to equation (2.1) w i l l be d i s c u s s e d . Some of these have been m e n t i o n e d b e f o r e i n S e c t i o n (1.4), f o r the one d i m e n s i o n a l c a s e . E x a m p l e 3.1. W i t h p = 1, u, = 1 and = t i = i ( t + t ,) one obtains the B o x scheme of K e l l e r (1969). w r i t t e n as -59-T h i s d i f f e r e n c e a p p r o x i m a t i o n m a y be w h e r e D - l = A l A ^ ~ ^ n - ^ -and D 0 = AT A ^ " • That the o r d e r of c o n s i s t e n c y i s two f o l l o w s f r o m T h e o r em (2. 1). Next i f p = 1 and i i = 2, w i t h I, = t . and L. = t then the m a t r i x f ' =1 n _ l '2 n equation (2. 3) becomes (3.1) I 0 0 I 0 0 •L •L T D - l T D 0 E T T 0,1 - LT  L 0 , 2 T 1,1 - LT h,2 T T T E l 2,1 " L 2 , 2 _ T LE2-0 0 0 w h e r e L. . = L nco 1(t ,)A(t ,) - oi X(t ,)B(t ,), 0 < i < 2 i , l 0 n - r v n - r v n - r x n - r — — and L. = L n o ) X ( t )A(t ) - (/(t )B(t ), 0 < i < 2 l, 2 0 n * n n x n — — W i t h o>0(t) = " ( t - t n ) / A t , (oV) = ( t - t n _ ! ) / A t a n d ^ 2 ( t ) = ( t - t n _ 1 ) ( t - t n ) / A t 2 , i t i s found that L 0 , l = A - t ^ W " B < W • L 0 , 2 = A l Mtn) L l , l = A T ^ n - l ) • L l , 2 = AT A( f cn> " B ^ • •60-L 2,1 AT A ( t n - i ) L i 2,2 AT A ( t ) v n F r o m (3.1) i t f o l l o w s that E,A(t ,) = E _ A ( t ).' Take E, =1. Then E. i s 1 n-1 2 v n 2 1 the s o l u t i o n of the m a t r i x equation E,A(t ,) = A ( t ). C l e a r l y , i t i s m o r e 1 n-1 n 1 d i f f i c u l t now to give e x p l i c i t r e p r e s e n t a t i o n s of the c o e f f i c i e n t s E^ than was the c a s e f o r s c a l a r equations i n Cha p t e r I. Of c o u r s e , one c a n alw a y s o b t a i n the c o e f f i c i e n t s n u m e r i c a l l y . In the s p e c i a l c a s e w h e r e A ( t n ^ ) and A ( t n ) commute, i t i s p o s s i b l e to state e x p l i c i t l y what the c o e f f i c i e n t s of the above f i n i t e d i f f e r e n c e equation a r e . In p a r t i c u l a r , i f A i s a constant m a t r i x , then w i t h E_> = 2 I i t i s found that E^ * h D - i = A T A f B ( t ,) and D n = - f - A - | B ( t ). * n-1 0 At ^ n Thus the d i f f e r e n c e equation may be w r i t t e n as A , n n-1. (u -u ) At ^ ( B ^ u ^ B ^ ^ u1 1 - 1 ) = | ( f ( t n ) + m^j) w h i c h i s r e c o g n i z e d as the u s u a l t r a p e z o i d a l method. The o r d e r of c o n s i s t e n c y of this e quation i s a l s o equal to two. E x a m p l e 3.2. W i t h p = 2, |_L = 1 and = t ^ + € A t the d i f f e r e n c e a p p r o x i -m a t i o n to (2.1) becomes A ^ At |(2M)un-2(e-Dun-1 + l(2e-3)u n- 2 + l(£-l)(£-2) B(^1) ie(e-i)un-e(e-2)u-x + n-1 (The m e s h i s u n i f o r m i n this example. ) The o r d e r of c o n s i s t e n c y of t h i s a p p r o x i m a t i o n i s equal to two. F r o m T h e o r e m (2.1) i t f o l l o w s that the o r d e r of c o n s i s t e n c y i n c r e a s e s to three i f one takes £ =^-±\^3 ~ 1-58. •61-If £ = 2 then the above f o r m u l a r e d u c e s to v n' A t 3 n _ n-1 , i n-2 — u - 2u + 2 U B ( t ) u n = f(t ) , n n w h i c h i s one of Gear's b a c k w a r d d i f f e r e n t i a t i o n f o r m u l a s . (See G e a r (1971), p. 217.) E x a m p l e 3. 3. C o n s i d e r the constant c o e f f i c i e n t p r o b l e m (3. 2) A w (t) - Bw(t) = f(t) , 0 < t < T , (3. 2a) w(0) = w , w h e r e A i s n o n s i n g u l a r and w h e r e A and B commute. T h i s o c c u r s f o r exa m p l e when A = I. L e t p. = 2 and p = 1. The c o n s t r u c t i o n p r o c e d u r e of S e c t i o n 2 a g a i n l e a d s to the m a t r i x equation (3.1) g i v e n i n E x a m p l e (3.1), now f o r g e n e r a l c h o i c e of the p o i n t s and t,^. L e t E1 = -L 0co 2(£ 2) A + ^{t,^ B. Then, u s i n g the f a c t that A and B commute, i t i s found that E 2 = L 0co Z(q)A - co 2(?, 1)B,. D _ j _ = E 1{L 0co°(c: i)A - u°(t>1)B}+ E 2{L 0co°(t: 2)A-co 0(C 2)B} , and D Q - E ^ L Q W ^ C ^ A - W ^ ^ B } + E.,{L.Q b i 1[t> 2 )A-( l i\t, 2)B}, where LQto(t) = co'(t). In p a r t i c u l a r i f one uses the G a u s s i a n p o i n t s -62-= ( t n - 1 + t n ) / 2 - AtN / 3 /6 and C 2 = ( ^ . i + t j / 2 + A t ^ 3 / 6 then the d i f f e r e n c e equation b e c o m e s (3. 3) { ^ A 2 - ± A B - § B 2 } u n " 1 + { ^ A 2 - l A B + | B 2 } u n = {|A + ^ I - B } f(Cl) 4- {±A - ^ - B } f(&2) The o r d e r of c o n s i s t e n c y of this equation i s equal to four. If f(t) = 0 then the f i n i t e d i f f e r e n c e a p p r o x i m a t i o n (3. 3) c a n a l s o be obtained f r o m a p a r t i c u l a r P a de a p p r o x i m a t i o n to the s o l u t i o n w(t) = B w ^ of (3.2), (3.2a). T hese a p p r o x i m a t i o n s , o r i g i n a l l y due to P a d e (1892) have been e x t e n s i v e l y s t u d i e d i n c o n n e c t i o n w i t h t h e i r s t a b i l i t y p r o p e r t i e s . (See e. g. , V a r g a (1961).) E x a m p l e 3.4. The p r e c e d i n g e xample suggests that the g e n e r a l p-step f i n i t e d i f f e r e n c e a p p r o x i m a t i o n to (3. 2) w i t h p. c o l l o c a t i o n p o i n t s has the f o r m i.=-p v = 0 1=1 v=0 where and (3^ a r e c o n s t a n t s that depend on the c h o i c e of the p o i n t s t,^, (1 < SL < p.). JL JL To v e r i f y t h i s and to show how the c o e f f i c i e n t s av and |3 m a y be found, c o n s i d e r the s c a l a r equation (3. 5) Lw(t) = w'(t) - X w(t) = f(t) , -63-The f i n i t e d i f f e r e n c e a p p r o x i m a t i o n to (3.5) i s g i v e n by (3.6) l = - p n+4 1=1 where (3.7) e. i-il |JL+4+l T P +1 / / X Lw (£ ) and (3.8) 6^  = J L w H + i ( C k ) e k , -P< k< 0 k=l The n o r m a l i z i n g f a c t o r S i s g i v e n by (2.4). C o n s i d e r i n g (3.5), (3.7) and (3. 8) i t i s o b s e r v e d that the c o e f f i c i e n t s have the f o r m [X u.-l v=0 v=0 4 4 f o r c e r t a i n constants a and 6 In o r d e r to show that these co n s t a n t s a r e the same as the cons t a n t s i n equation (3.4), note that by d e f i n i t i o n •64-Y e | L a ) k ( ^ ) = 0 , p + l < k < p + HL-1 4=1 i . e . , .-1 4=1 v=0 fo r a l l X and f o r p + l < k < p + [ i - l = o , So q^(X.)= 0 and t h e r e f o r e ° x k ( C ) = 0 f o r a l l m a t r i c e s C. In p a r t i c u l a r w i t h C = A X B one obtains H u-1 4=1 v=0 V L Q C O ^ ^ J I - C O ^ C ^ A ^ B 0 . Upon m u l t i p l i c a t i o n b y A ^ and u s i n g the a s s u m p t i o n that A and B commute one obtains p. u-1 4=1 v=0 L Q c o k ( ^ ) A - o o k ( ^ ) B 0 , p + l < k < p + [ i - l . Thus i f one define s .-1 v=0 then the l a s t u.-2 r e l a t i o n s i n the m a t r i x equation (2.3) a r e s a t i s f i e d . F u r t h e r , the f i r s t p+l r e l a t i o n s of (2.3) a r e a u t o m a t i c a l l y s a t i s f i e d by l e t t i n g P-D - S I E v L 0cop +^(c : V ) A(c: v)- w p + i(c V ) B(c: v) T h e r e f o r e i t has been shown, that i n c a s e of the constant c o e f f i c i e n t p r o b l e m (3.2), w i t h c o m m u t i n g A and B, one c a n e x p l i c i t l y g i v e the c o e f f i c i e n t m a t r i c e s and of the d i f f e r e n c e a p p r o x i m a t i o n (2.2). 4. T H E S T A B I L I T Y O F F I N I T E D I F F E R E N C E A P P R O X I M A T I O N S TO  I N I T I A L V A L U E P R O B L E M S In t h i s s e c t i o n some w e l l - k n o w n concepts i n the s t a b i l i t y a n a l y s i s of f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s to s y s t e m s of f i r s t o r d e r o r d i n a r y d i f f e r e n t i a l equations w i l l be s u m m a r i z e d . To define s t a b i l i t y of a f i n i t e d i f f e r e n c e a p p r o x i m a t i o n one only has to c o n s i d e r the f o r m this d i f f e r e n c e i equation takes when a p p l i e d to the equation w (t) = 0 . H e r e w(t) i s a s c a l a r f u n c t i o n . Thus f o r the type of a p p r o x i m a t i o n s d i s c u s s e d i n S e c t i o n 2 the r e l e v a n t d i f f e r e n c e e quation i s 0 (4.1) Y, &i u n + i = 0 , p < n < N , i = - p w h e r e v ^ v v p + 1v • • V ^ X w i t h LQCo(t) = to'(t). The m e s h i s taken to be u n i f o r m i n t h i s s e c t i o n , so that ( A t ) n = At = 1/N, (1 < n < N). The a p p r o x i m a t i o n i s s a i d to be s t a b l e p r o v i d e d that f o r a l l g i v e n i n i t i a l data u 1 1, (0 < n < p-1), the s o l u t i o n to (4.1) s a t i s f i e s 5 i = (-1Y -66-p-1 | u n | < const Y I ' p < n < N. 1=0 T h i s i s the c a s e i f a l l r o o t s r)^ of the c h a r a c t e r i s t i c e quation 0 (4.2) <r0(rj) = At £ \ r f + l = 0 i = - p s a t i s f y (4. 3) |r). | < 1 , 1 < i < p , and (4.3a) |n.J = 1 i m p l i e s r\^ i s a s i m p l e r o o t . T h e s e two c o n d i t i o n s a r e u s u a l l y r e f e r r e d to as the r o o t c o n d i t i o n . F o r p r a c t i c a l p u r p o s e s i t i s d e s i r a b l e that the a p p r o x i m a t i o n i s s t r o n g l y s t a b l e , i . e . , the r o o t s s a t i s f y 1j_ = 1 . r? i | < 1 , 2 < i < p . E x a m p l e 4.1. If p = 1 then the c o e f f i c i e n t s i n (4.1) a r e g i v e n by 6 Q = and 5 i ~ ^7/ • T n e c h a r a c t e r i s t i c equation has one r o o t r) = 1 and thus thi s a p p r o x i m a t i o n i s s t a b l e . If p = 2 and p = 1 w i t h t, = t + £ At, then the d i f f e r e n c e 1 n-2 --67-approximation to w'(t) = 0 has the f o r m T n-2 , -r n-1 -r- n n O Q u + u + 5., u = 0 , where the coefficients 6^  a r e given by S Q = (2£-3)/2 At, 6^  = (4-4£)/2At and S_, = (2£,-l)/2At. The roots of the c h a r a c t e r i s t i c equation (4.2) a r e ri1 = 1 and r\ 2 = (2£-3)/(2£-l). Now | r j 2 | < 1 i f and only if £ > 1. Thus the approximation i s strongly stable i f and only i f £ > 1. In p a r t i c u l a r , the second o r d e r f o r m u l a of Gear, for which % = 2, and the t h i r d o r d e r f o r m u l a obtained i f % - 1 + -|- N/3~ a r e strongly stable. (See example (3.2).) Now let p = 2 and p. = 2 with C,^  = t ^ + £]_ A t and £ 2 = t 2 + ^2 In this case the c o e f f i c i e n t s a r e 6(^+e 2) - 6 ^ | 2 - 8 -3(e 1+e 2) + 6e xe 2 + 2. and 6 Q = - bY - 6 2 . The c h a r a c t e r i s t i c equation °"Q( " ) = 0 has roots n ^  = 1 , and * i -9( e x + e 2) + 6 +14 § 2 -3(e x + e 2) + 6 ^ e 2 + 2 F o r the approximation to be strongly stable i t i s n e c e s s a r y that |i?2|< !• 2 If £ 2 = 2 then this i s the case p r o v i d e d that ^ < 0 or 1^ > "J* ^ ^2 = ^ 5 14 then i t i s n e c e s s a r y that > and i f = 0 then one needs — < ^ < 2. 62 = b2 n 1 A t -68-Another well-known concept i s that of stability r e g i o n of a finite d i f f e r e n c e approximation. If i n the s y s t e m (2.1) the m a t r i x A has eigenvalues with l a r g e negative r e a l part, then a l a r g e c l a s s of d i f f e r e n c e approximations r e q u i r e s At to be v e r y s m a l l for the approximation to be stable, even though this r e q u i r e m e n t may not be n e c e s s a r y for reasons of a c c u r a c y . Thus i t i s of i n t e r e s t to c h a r a c t e r i z e those approximations that do not suffer f r o m this shortcoming. F o r investigating the effect of the eigenvalues upon the st a b i l i t y i t i s c u s t o m a r y to apply the finite d i f f e r e n c e approximations under c o n s i d e r a t i o n to the p r o b l e m j_ ^w(t)= w (t) - \w(t) = 0, where w(t) i s a s c a l a r function. The a n a l y s i s r e m a i n s v a l i d for the s y s t e m Aw'(t) = Bw(t), where w i s a vector and A and B a r e constant m a t r i c e s , p r o v i d e d that A ^B has a complete set of eigenvectors. F i n i t e d i f f e r e n c e approximations to the s c a l a r equation w'(t) = X.w(t) have the f o r m 0 n+l Y 6L u n ' * = 0 , P < n < N , where the c o e f f i c i e n t s 6^  a r e given by (-l)1"-L ) A. U. F o r s t a b i l i t y i t i s again n e c e s s a r y that the roots rj. of the c h a r a c t e r i s t i c 0 equation q(rj) = At Y 8.r] P = 0 sa t i s f y the root condition (4.3), (4.3a). i=-P F r o m (4.4) i t follows that q(r?) has the f o r m -69-(4.5) ' q(n) = <r0(n) + £ ( At X ) V cr^n) , v=l w h e r e ^ ( r j ) i s g i v e n b y (4.2) and w h e r e cr (r?) e P . y p The s t a b i l i t y r e g i o n S a s s o c i a t e d w i t h a f i n i t e d i f f e r e n c e a p p r o x i -m a t i o n of the f o r m (2.2) i s now defined as that r e g i o n of the c o m p l e x A t X - p l a n e f o r w h i c h a l l r o o t s rj^(AtX) of q(rj) = 0 s a t i s f y l 1 ? ^ ^ 1- T hus i f a g i v e n d i f f e r e n c e a p p r o x i m a t i o n i s a p p l i e d to the equation w'(t) = Xw(t), and i f AtX eS, then the a p p r o x i m a t e s o l u t i o n u 1 1 w i l l tend to z e r o as n —••co fo r any c h o i c e of i n i t i a l data. F o r " s t i f f p r o b l e m s " , i . e . , f o r s y s t e m s that have an eigenvalue X w i t h Re(X) < < 0, i t i s t h e r e f o r e d e s i r a b l e that the s t a b i l i t y r e g i o n S i n c l u d e a c o n s i d e r a b l e p o r t i o n of the negative h a l f plane. In fact, f o r c e r t a i n a p p l i c a t i o n s i t i s n e c e s s a r y that S c o n t a i n the e n t i r e negative r e a l a x i s . (See S e c t i o n 6.) Note that an a p p r o x i m a t i o n s a t i s f i e s the r o o t c o n d i t i o n i f 0 e 5S.. H e r e 6S i s the b o u n d a r y of S. A n e c e s s a r y c o n d i t i o n f o r the negative r e a l a x i s to be c o n t a i n e d i n S i s of c o u r s e that the method i s s t a b l e at oo, i . e . , a l l r o o t s of cr (r\) = 0 l i e on o r i n s i d e the unit c i r c l e . (If the r o o t s a r e s t r i c t l y c o n t a i n e d i n the unit c i r c l e then the method w i l l be c a l l e d s t r i c t l y s t a b l e at co. ) T h e r e a r e d i f f e r e n c e a p p r o x i m a t i o n s however, that a r e s t r i c t l y s t a b l e at co, but f o r w h i c h the negative r e a l a x i s i s not c o n t a i n e d i n S. (See V a r a h (1975).) A r e l a t e d concept, i n t r o d u c e d by Widlund (1967), i s that of A{a)  s t a b i l i t y . A d i f f e r e n c e method i s s a i d to be A(a) - s t a b l e i f S c o n t a i n s a r e g i o n of the f o r m S = { z : TT-a < a r g z < t r + a > z # 0 } . -70-Thus i f an a p p r o x i m a t i o n i s A ( a ) - s t a b l e , then the negative r e a l a x i s i s c o n t a i n e d i n S and hence i n S. a The s t a b i l i t y p r o p e r t i e s of e x p r e s s i o n s of the f o r m (4. 5), ( " f i n i t e d i f f e r e n c e f o r m s " ) , have been e x t e n s i v e l y i n v e s t i g a t e d for the c a s e u. = 1 by, f o r example, D a h l q u i s t (1959 and 1963), G e a r (1971) and V a r a h (1975). The g e n e r a l c a s e has been s t u d i e d b y R e i m e r (1968). The m o t i v a t i o n f o r c o n s i d e r i n g g e n e r a l f i n i t e d i f f e r e n c e f o r m s i s the f a c t that the s t a b i l i t y a n a l y s i s of many d i f f e r e n c e methods le a d s to studying such f o r m s . E x a m p l e s of these i n c l u d e methods based upon Pade r a t i o n a l a p p r o x i m a t i o n to the e x ponential, (see V a r g a (1961)), Runge K u t t a methods and second d e r i v a t i v e methods, ( E n r i g h t (1974).) In th i s s e c t i o n i t has been shown how these h i g h e r o r d e r f i n i t e d i f f e r e n c e f o r m s a r i s e i n the s t a b i l i t y a n a l y s i s of the v e r y g e n e r a l type of f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s c o n s i d e r e d i n S e c t i o n s 2 and 3. It i s not the p u r p o s e of t h i s s e c t i o n to e x t e n s i v e l y c o n t r i b u t e to the i n v e s t i g a t i o n s r e f e r r e d to above, but by means of some e x a m p l e s the e f f e c t that the c h o i c e of the points t,^ has on the s t a b i l i t y r e g i o n S w i l l be i l l u s t r a t e d . E x a m p l e 4.2. L e t p = 1 and u. = 1 w i t h 2,^  = t ^ + £ At. So the d i f f e r e n c e a p p r o x i m a t i o n to (_ ^ w ( t ) = w'(t) - Xw(t) = 0 has the f o r m ~ n-1 ~ n _ - l u 0 U = ' w h e r e the c o e f f i c i e n t s 6 and 8^ a r e found to be equal to * - i = ii - x a n d ~6o = it - x ^ • T h e r e f o r e the c h a r a c t e r i s t i c p o l y n o m i a l i s g i v e n by -71-q(") = 6QT7 + 6_1 = tr Q(rj) + A t X ^ ( n ) , where <rQ(r)) = 17 - 1 and cr^rj) = -£"+1-1. The r o o t of 0^(77) = 0 i s 77^  = (£-l)/£, so that |rj^|< 1 i f and o n l y i f £ > \. Thus this d i f f e r e n c e a p p r o x i m a t i o n i s st a b l e at 00 i f £ It i s easy to ch e c k that the a p p r o x i m a t i o n i s i n fa c t A(a)-stable f o r suc h £. The method i s s t r i c t l y s t a b l e at 00 only i f £ > \. When a p p l i e d to the g e n e r a l s y s t e m (2.1) the c h o i c e s £ = 0, £ = \ and £ = 1 c o r r e s p o n d to the s t a n d a r d e x p l i c i t one step E u l e r a p p r o x i m a t i o n , the B o x scheme, ( K e l l e r (1969) ), and the t o t a l l y i m p l i c i t one step a p p r o x i m a t i o n . E x a m p l e 4.3. If one takes p = 2 and u. = 1 w i t h ^ = t + £ At , then q(r?) = °- 0(r7) + AtXo-^rj), w h e r e <rQ{n) = (2^-1)^ + 4(1-£)n + 2£-3 and 0^(77) = r)Z + Zi{i-2)r\ - (£ 2-3£+2). In E x a m p l e (4.1) i t was found that th i s d i f f e r e n c e a p p r o x i m a t i o n i s st a b l e at z e r o i f and only i f £ > 1, i . e . , the r o o t s of °"Q(") = 0 l i e o n o r i n s i d e the unit c i r c l e i f and only i f £ >.l. The a p p r o x i m a t i o n i s stable at 00 i f the r o o t s of 0^(77) = 0 l i e on o r i n s i d e the unit c i r c l e . T h ese r o o t s a r e g i v e n by Some co m p u t a t i o n r e v e a l s that 177.. | < 1 p r o v i d e d £ > 1 + \ \l2 « 1. 7. F o r example, the second o r d e r f o r m u l a of G e a r i s s t r i c t l y s t a b l e at 00, but the t h i r d o r d e r f o r m u l a c o r r e s p o n d i n g to £ = 1 + - j \l3 , (see E x a m p l e (3.2)), i s not stabl e at 00 . 5. I N I T I A L B O U N D A R Y V A L U E P R O B L E M S . T h i s s e c t i o n i s c o n c e r n e d w i t h the c o n s t r u c t i o n of f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s to l i n e a r s e c o n d o r d e r p a r a b o l i c p a r t i a l d i f f e r e n t i a l equations of the f o r m -72-(5.1) p(x, t) y f c(x, t) = L(t) y(x, t) + f(x, t ) , 0 < x < 1 , 0 < t < T , where L(t) y ( x , t ) = y (x, t) + q(x, t) y (x, t) + r(x, t) y(x, t), XX X with i n i t i a l condit ion (5.1a) y(x, 0) = g(x) , 0 < x < 1 , and boundary conditions (5.1b) a 0 y(0 , t ) + b 0 y x ( 0 , t ) = gQ(t) , 0 < t < T , and (5.1c) a l Y ( l , t) + b j y ^ l , t) = g :(t), 0 < t < T . It i s assumed that 0 < p , <p(x, t) < p < co, so that ( l . i i ) i s parabol ic i n the sense of P e t r o v s k i i (1938), and that the above p r o b l e m has a unique solution, that i s as many t imes continuously differentiable as necessa ry . Define a m e s h R^(T) by (5. 2) R h ( T ) = { ( X j , t n ) : 0 = x Q < X ; L < - < X j = 1 ; 0 = t Q < tj < ' • < t N = T }. The mesh spacings are h . = x . - x . ^, (1 < j < J ) , and ( A t ) n = t n ~ t n y (l<n<N). L e t h = max h. and At = max (At) . So for u n i f o r m meshes h. = h= -r and j J n n J J 1 n ( A t ) ^ = At = j^r . F o r functions w(x>t) defined on R^(T) let w. = w ( x j > tn)» n n , n n w, = (w,, w_, w n .T J - l ' ' n w. = max h J - l i n , 2 i 2 |w^| and ||w£||2= { J (w*) } . F u r t h e r , for m a t r i c e s || • || and || • || ^  denote the m a t r i x norms induced by the c o r r e s p o n d i n g v e c t o r norms. In p r i n c i p l e i t i s p o s s i b l e to extend the co n s t r u c t i o n method of Chapter I to p a r t i a l d i f f e r e n t i a l equations. However, this approach w i l l not be taken here, but instead d i s c r e t i z a t i o n i n space and time w i l l be done separately. D i s c r e t i z i n g i n space f i r s t r e s u l t s i n a s y s t e m of l i n e a r f i r s t o r d e r o r d i n a r y d i f f e r e n t i a l equations. T h i s s y s t e m can then be solved by a di f f e r e n c e method of the type d i s c u s s e d i n Sections 2.3 and 4. M o r e p r e c i s e l y this s y s t e m i s given by (5.3) K h ( t ) p ( X j , t)w (x.,t) = L h(t)w(x., t) + K h(t)f(x., t), 1 < j < J - l , where and with i n i t i a l data L h ( t ) w(x., t) = J d. j.(t)w(x. + i, t) , i = - r . J m J K (t)w(x t) = V e .(t)w(z t) , 1 1 J J' 1 J >1 (5. 3a) w ( X j , 0) = g(x.) , 0 < j < J . -74-F o r e ach t the c o e f f i c i e n t s cL(t) and e^(t) a r e d e f i n e d i n p r e c i s e l y the same manner as the c o e f f i c i e n t s d^ and e^ i n C h a p t e r I. Thus these c o e f f i c i e n t s a r e g i v e n by equations (1.2.10) and (1.2.11) w i t h L(t) r e p l a c i n g L. (To s i m p l i f y the notation, the s u b s c r i p t j i n d^ .., e^ ^ etc. w i l l u s u a l l y be omitted.) The i m p o r t a n t a s s u m p t i o n i s made that the c o l l o c a t i o n p oints c o i n c i d e w i t h the m e s h points x. ., i . e . i f z. i s a c o l l o c a t i o n p o i n t then x J + 1 x f z^ = f° r some i n t e g e r I . F o r example, i f m = 1 one may have K h ( t ) w ( x j t t ) = e 1(t) w ( X j , t ) = w ( X j , t ) , 1 < j < J - l , and i f m = 3 1 K h ( t ) w ( x . ) t ) = £ e.(t) w(x. +.,t) , l < j < J - l . i = - l T h i s r e s t r i c t i o n on the c h o i c e of the c o l l o c a t i o n p oints i s n e c e s s a r y f o r (5.3) to c o n s t i t u t e a p r o p e r s y s t e m of o r d i n a r y d i f f e r e n t i a l equations. The t r u n c a t i o n e r r o r T(t) a s s o c i a t e d w i t h (5. 3) i s def i n e d as m s m r(t)£ T e.(t) p(z., t) y t ( z . , t ) - Y d . ( t ) y ( x ., t) - V e.(t) f(z., t) , i=l i = - r J i = l where y(x, t) i s the exact s o l u t i o n of (5.1) su b j e c t to (5.1a,b,c). T h e r e -f o r e , u s i n g the d i f f e r e n t i a l equation, m s T(t) = I e.(t) L y ( z . , t ) - Y d i ( t ) y(x, i, t ) . i = l i = - r J Since the c o e f f i c i e n t s d^(t) and e^(t) a r e d e f i n e d as i n S e c t i o n d . 2) one may -75-a p p l y T h e o r e m (I. 2.1) and T h e o r em (I. 2 . 4) to a r r i v e at the f o l l o w i n g . T h e o r e m 5. 2. F o r each t, (0 < t < T), l e t d.(t), e^t) be defined b y equations (I. 2. 10) and(I. 2. 11) r e s p e c t i v e l y , w i t h the n o r m a l i z i n g f a c t o r E as i.n(I.2.12). A s s u m e that f o r each i , ( l < i < m ) , z^ = x j + ^ _ f° r some i n t e g e r S.^. A l s o a s sume that r+s > n and that E ± 0. Then t h e r e e x i s t s a constant C that does not depend on h, s u c h that | T ( t ) | < C h r + s + m " 2 In p a r t i c u l a r , i f r+s = n o r i f m = 1 then E + 0. F o r each t, (0 < t < T), there a r e a l s o a p p r o x i m a t i o n s to the bounda r y c o n d i t i o n s (5.1b, c) of the f o r m s 0 (5.3b) YJ d Q ) . w(x., t) = g Q ( t ) , 0 < t < T , i=0 and 0 (5.3c) Y d J , i w ( x J + i ' f c ) = Ei(t) , 0 < t < T . i = " r J T h e se a p p r o x i m a t i o n s may be defi n e d as i n S e c t i o n (I. 5)» E q u a t i o n (1.5. 5)> w i t h m = 0. Thus these d i s c r e t e boundary c o n d i t i o n s a r e a s s u m e d to be independent of the d i f f e r e n t i a l equation. It i s e a s y to show that i n t h i s c a s e the c o e f f i c i e n t s d^ Q and d j Q a r e non z e r o . T h i s a l l o w s e l i m i n a t i o n of the unknowns W(XQ, t) and w ( x j , t ) i n the s y s t e m (5.3), so that one i s l e f t w i t h a s y s t e m of J-1 o r d i n a r y d i f f e r e n t i a l equations i n the J-1 v a r i a b l e s w ( X j , t), (1 <_ j <_ J - l ) . T h i s r e d u c e d s y s t e m w i l l be c o m p a c t l y w r i t t e n as (5.4) p y t ) w h ( t ) = L h ( t ) w h ( t ) + f h ( t ) , 0 < t < T w h e r e w-^(t) = (w(x^, t), "", w ( x j y t)) and f^(t) i s the a p p r o p r i a t e inhomogeneous t e r m . The i n i t i a l data a r e (5.4a) w h(0) = g h = (g(x 1), g ( x J _ 1 ) ) T The d i s c r e t i z a t i o n c a n now be c o m p l e t e d by a p p l y i n g one of the methods f o r s y s t e m s of o r d i n a r y d i f f e r e n t i a l equations d i s c u s s e d i n S e c t i o n s 2, 3 and 4. If the s p a t i a l d i s c r e t i z a t i o n as w e l l as the d i s c r e t i z a t i o n i n time a r e c o n s i s t e n t , then i t f o l l o w s f r o m the eq u i v a l e n c e t h e o r e m of L a x , (see R i c h t m e y e r and M o r t o n (1967), p. 45), that the s o l u t i o n of (5.4), (5.4a) c o n v e r g e s t o the s o l u t i o n of the continuous p r o b l e m (5.1), (5.1a, b, c) i f and only i f the a p p r o x i m a t i o n i s s t a b l e . The s t a b i l i t y p r o b l e m i s d i s -c u s s e d i n the next s e c t i o n , w h e r e a s s p e c i f i c e x a m p l e s of f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s w i t h n u m e r i c a l tes t s a r e g i v e n i n S e c t i o n 7. 6. T H E S T A B I L I T Y O F F I N I T E D I F F E R E N C E A P P R O X I M A T I O N S TO  I N I T I A L B O U N D A R Y V A L U E P R O B L E M S In t h i s s e c t i o n the n u m e r i c a l s t a b i l i t y of f i n i t e d i f f e r e n c e a p p r o x i -m a t i o n s to l i n e a r p a r a b o l i c p a r t i a l d i f f e r e n t i a l equations i n one space v a r i a b l e w i l l be i n v e s t i g a t e d . F o r s i m p l i c i t y , we c o n s i d e r the d i f f u s i o n equation -77-(6.1) y, (x, t) = L y ( x , t) = y (x, t) , 0 < x < 1, 0 < t < T w i t h i n i t i a l c o n d i t i o n (6.1a) y(x, 0) = g(x), 0 < x < 1 and boundary c o n d i t i o n s (6. lb) y(0, t) = 0 , 0 < t < T , and (6.1c) y (1, t) = 0 , 0 < t < T D i s c r e t i z i n g i n space f i r s t as i n S e c t i o n 5 g i v e s the equations m. Si (6.2) V e~. w'(z , t) = Y d w(x , t), 1 < j < J-1, 0 < t < T , i=l i = - r . J w h e r e the c o e f f i c i e n t s d. . and e. . a r e g i v e n by equation (1.2.16) and j , i j , i (1-2.17) r e s p e c t i v e l y , and w h e r e the c o l l o c a t i o n points z. . a r e r e q u i r e d J > 1 to c o i n c i d e w i t h the m e s h p o i n t s . F u r t h e r , the m e s h i s a s s u m e d to be 1 1 u n i f o r m , so that h. = h = — and ( A t ) n = A t = |^ • The boundary c o n d i t i o n (6.1c) i s a p p r o x i m a t e d by 0 Y b. w ( x J + . , t) = 0 , 0 < t < T . i = " r J A s i n the p r e v i o u s s e c t i o n , one may use the d i s c r e t e b o u n d a r y c o n d i t i o n s - 7 8 -to eliminate W ( XQ , t) and w(xj,t) f r o m (6.2). The reduced s y s t e m of o r d i n a r y d i f f e r e n t i a l equations i s then w r i t t e n as (6. 3) K h w^t) = L h w h(t) , with i n i t i a l data (6.3a) w h(0) = g h . B e f o r e d i s c r e t i z i n g i n time sufficient conditions for the solution of (6. 3), (6.3a) to converge to the solution of (6.1), (6.1a,b, c) w i l l be stated. T h e o r e m 6.1. The solution w ^ ( t ) of (6. 3), (6. 3a) converges to the solution of (6.1), (6.1a,b,c), and this d i s c r e t i z a t i o n i s stable i n the i^-n o r m with r e s p e c t to perturbations i n the i n i t i a l data, p r o v i d e d that there exists, a constant c^.; that,.does , no.tv depend on h, such that for a l l s m a l l enough h the following conditions are s a t i s f i e d . 1. The spa t i a l d i s c r e t i z a t i o n (6.3) i s consistent with (6.1), (6.1b,c). 2. K, i s nonsingular. h 3. The eigenvalue p r o b l e m L ^ v ^ = ^ h ^ h v h a d m i * - s a complete set of eigenvectors. 4. The eigenvalues X., ., (1 < i < J-1), of the above eigenvalue p r o b l e m n, I s a t i s f y Re(\, .) < 0. ' h, x — 5. L e t be the m a t r i x of eigenvectors. Then ||||^ j |"V | | ^ — ° 2' Proof. L e t eh(t) = y h ( t ) - w h(t). Then K h e^(t) = L h e h ( t ) + T h ( t ) . T h e r e f o r e t K T 1 L, t (t-r)K~lU e h(t) = e n n e h ( 0 ) + / e & T h ( r ) dr t A , t ( t - r ) A , , = Vh 6 Veh<°» + / 0 V h e h V ; 1 T h ( r ) d r > w h e r e A ^ i s the d i a g o n a l m a t r i x of e i g e n v a l u e s . If n g i s the o r d e r of c o n s i s t e n c y of (6. 3), then the e s t i m a t e |K(t) || 2 < c,{ || e (0)|| + T -max.. || T , ( t ) n c c i i c o < t < T f o l l o w s i m m e d i a t e l y f r o m the a s s u m p t i o n s . A f u l l y d i s c r e t i z e d f i n i t e d i f f e r e n c e a p p r o x i m a t i o n to (6. 1 ) , (6. l a , b, c) i s obtained b y a p p l y i n g one of the c l a s s of methods d i s c u s s e d i n S e c t i o n s 2, 3 and 4 to the s y s t e m (6.3). F r o m the d i s c u s s i o n i n E x a m p l e (3.4) i t f o l l o w s that i f and commute then th i s a p p r o x i m a t i o n has the f o r m 0 _ u. J A I t -.11 T i l .A % = 0 , 1 < n < N V ^h l=-p v=0 The c o e f f i c i e n t s a £ a r e a s s u m e d to be independent of n, so that the same d i f f e r e n c e a p p r o x i m a t i o n i s used at each t i m e - s t e p . The i n i t i a l data a r e u^ = g^. F u r t h e r , i f p > 1 then s u f f i c i e n t l y a c c u r a t e a p p r o x i m a t i o n s to uj^ , ( 1 <_ n < P - l ) , m u s t be g i v e n a l s o . T h e s e may be obtained f o r example by i n i t i a l l y u s i n g a one step method, (p = l ), having the same o r d e r of c o n s i s t e n c y as (6.4). The d i f f e r e n c e a p p r o x i m a t i o n (6.4) i s s a i d to be stab l e i n the l ^-n o r m i f t h e r e e x i s t s a co n s t a n t K, that does not depend on h and A t , such that -80-Ki i 2 < KiZ"1 K i i 2 > 1 / 2 • p < n < N ' 1=0 f o r a l l s m a l l enough h and A t . S u f f i c i e n t c o n d i t i o n s f o r s t a b i l i t y , and hence f o r conv e r g e n c e , a r e g i v e n i n the f o l l o w i n g t h e o r e m . T h e o r e m 6.2. The s o l u t i o n u^ of (6.4), w i t h a p p r o p r i a t e i n i t i a l c o n d i t i o n s , c o n v e r g e s to the s o l u t i o n of (6. 1 ), (6.1a, b, c), and thi s d i f f e r e n c e s c heme i s s t a b l e i n the i ^ - n o r m w i t h r e s p e c t to p e r t u r b a t i o n i n the i n i t i a l data, p r o v i d e d that there e x i s t p o s i t i v e constants c^ and c^, that do not depend on h, s u c h that f o r a l l s m a l l enough h the f o l l o w i n g c o n d i t i o n s a r e s a t i s f i e d . 1. The d i f f e r e n c e scheme (6.4) i s c o n s i s t e n t . 2. i s n o n s i n g u l a r . 2a. K, and L, commute, h h 3. The eigenvalue p r o b l e m L*^v^ = X^K^v^ a d m i t s a c o m p l e t e set of e i g e n v e c t o r s . 4. The eig e n v a l u e s \, . of the above eigenvalue p r o b l e m l i e i n a r e g i o n n, I T of the f o r m a T = {z e C : T - a < a r e z < ir +a ; z ± 0 \ , w h i c h has the p r o p e r t y that i f z e T , then the r o o t s f)A,z) °f the c h a r a c t e r i s t i c equation L=-p v=0 s a t i s f y - 8 1 -\rii(z)\< 1 , and 1 ^ ( 2 ) 1 <1 - ci , 2 < l < p . 5. Le t be the m a t r i x of e igenvectors . Then l v h l 2 I I ^ ' I U < c 2 Proof . By assumption 3 one can wr i te u^ = c^. Substitution into (6.4), a lso using 2 and 2a, gives 0 p. K h L I % A t Kfa L h V h C h 4=-p v=0 0 p. AT K h L L a v A t v h A h c h 4=-p v=0 0 p. X t K h v h L h % A t A h c h = 0 -4 = -p v=0 where is the diagonal m a t r i x of eigenvalues. Hence each component cl^ j of satisfies the difference equation 0 ^ (6.5) 2 YJ \, i'ATK i ) V c h + i = °> P < n < N = T / A t . l=-p v=0 v n ' x ' J 1 +1 ^ Le t c!1 . denote the vector (c-P . , c^ f , " ' , c5* P + 1 ) and let - h , J h , j h , j h , j P"1 2 | II • | L = { E I c S " m I } 2 - Define the p X p m a t r i x A(z) by h>3 L m=o A' J -82 -A ( z ) = • p p - l ( z ) P p ( z ) 1 0 " P p - 2 ( Z )  P p ( z ) 0 1 P p (z ) 0 0 1 0 m v where p (z) = 7. v = one step difference equation ^ a z' , 0 < ^ m < ^ p . Then (6.5) can be wr i t t en as a 0 :JJ . = A(AtX, -) c ? " . 1 , p < n < N , 1 < j < J - l The eigenvalues I j ^ A t X ^ j) of A ( A t X ^ ..) are p r e c i s e l y the roots of the cha rac te r i s t i c equation 0 n The definit ion of T i m p l i e s that i f X, . e T then AtX, . e T . F r o m a r h , j a h , j a assumption 4 i t now follows that |n ; 1(Atx h J)| < i , and |n (AtX, .)| < 1-c, , 2 < m < p , 1 m h , j 1 — 1 ' — for a l l j , (1 < j < J - l ) . It is w e l l known, (See R ich tmeyer and M o r t o n (1967), p . 86), that -83-this i m p l i e s that the f a m i l y of m a t r i c e s JA(z) : z e T } s a t i s f i e s || ( A ( z j ) n | | z < K x , p < n < N . Thus under these conditions one has n-p+1 l c S . j l S b S . , 1 2 = l K , i ) ^ . j l l z ^ l K j I l z and t h e r e f o r e J=1 ' J j=l ' J j=l m=0 ' J m=0 so that K » 2 = l | V h c J J I I 2 < | V h l l 2 l | c S l l 2 < I I V J I . K , ! J J c ™ \ \ l f m=0 = "Viz Il^lll}*! I I V J U I V ^ K / Z w^wif m=0 m=0 m=0 T h i s proves the theorem. Bef o r e proceeding to analyze a number of examples, some r e m a r k s ar e i n o r d e r . F i r s t , it should be noted, that conditions 3, 4 and 5 of T h e o r e m (6.1) and Theorem(6. 2) are i n g e n e r a l d i f f i c u l t to v e r i f y a n a l y t i c a l l y . However, a good check on whether a proposed d i s c r e t i z a t i o n leads to a stable d i f f e r e n c e scheme, is obtained by n u m e r i c a l computation of the eigenvalues and the - 8 4 -c o n d i t i o n n u m b e r K (V^) = I  V j J | 2 I  ^ i / l l 2 a ^ e w v a ^ n e s °^ k . T h e a s s u m p t i o n that a l l e i g e n v a l u e s l i e i n the lef t h a l f p l ane i s not s t r i c t l y n e c e s s a r y , but r e a s o n a b l e i n v i e w of the d i f f e r e n t i a l e q u a t i o n u n d e r c o n s i d e r a t i o n . In p r o b l e m s w h e r e t h e r e a r e s o m e , but not m o r e t h a n a f i x e d f i n i t e n u m b e r of e i g e n v a l u e s w i t h p o s i t i v e r e a l p a r t , one n e e d o n l y a d d the r e q u i r e m e n t that the m o d u l u s of these e i g e n v a l u e s c a n be u n i f o r m l y b o u n d e d . T h e e s s e n t i a l d i f f e r e n c e b e t w e e n the s t a b i l i t y a n a l y s i s i n t h i s s e c t i o n and the s t a b i l i t y a n a l y s i s of d i f f e r e n c e m e t h o d s fo r s y s t e m s of o r d i n a r y d i f f e r e n t i a l equa t ions i n S e c t i o n 4 i s the fac t that the c u r r e n t s y s t e m , (equa t ion ( 6 . 3 ) ) , does not have a f i x e d d i m e n s i o n . T h i s r e s u l t s i n the r e q u i r e m e n t s that K ( V J ^ be bounded a n d that the s t a b i l i t y r e g i o n a d m i t a r e g i o n T as i n T h e o r e m ( 6 . 2 ) . A m u l t i s t e p m e t h o d that a d m i t s s u c h a r e g i o n w i l l be c a l l e d s t r o n g l y A(a)-stable. (See S e c t i o n 4 f o r a d e f i n i t i o n of A ( a ) s t a b i l i t y . ) T h u s i f a m e t h o d i s s t r o n g l y A ( a ) - s t a b l e , t hen i t i s A ( a ) - s t a b l e , but the c o n v e r s e n e e d not be t r u e . A l l A ( a ) - s t a b l e one s tep m e t h o d s a r e s t r o n g l y A ( a ) - s t a b l e , b e c a u s e the c h a r a c t e r i s t i c e q u a t i o n q(n) = 0 has o n l y one r o o t i n that c a s e . O t h e r e x a m p l e s of s t r o n g l y A ( a ) - s t a b l e m e t h o d s a r e the f o r m u l a s of G e a r , (see e x a m p l e (3 .2) a n d e x a m p l e ( 1 . 4 . 2 ) ) , a n d the s e c o n d d e r i v a t i v e m e t h o d s g i v e n i n E n r i g h t (1974). T h e s e m e t h o d s a r e u s e d i n the n u m e r i c a l e x a m p l e s of S e c t i o n 7 . C o n d i t i o n 4 of T h e o r e m (6.2) i s c l o s e l y r e l a t e d to the n o t i o n of p a r a b o l i c i t y of a d i f f e r e n c e s c h e m e , i n t r o d u c e d by W i d l u n d (1965 and 1966). In these p a p e r s the s t a b i l i t y of v e r y g e n e r a l m u l t i s t e p m e t h o d s f o r l i n e a r p a r a b o l i c s y s t e m s w i t h o u t b o u n d a r y c o n d i t i o n s i s c o n s i d e r e d , w i t h F o u r i e r t r a n s f o r m s r e p l a c i n g the e i g e n v a l u e - e i g e n v e c t o r a p p r o a c h i n t h i s s e c t i o n . (See a l s o H a k b e r g (1970) a n d N o r d m a r k (1974).) T h e m a j o r d i s advan t age -85-of the c u r r e n t e i g e n s y s t e m a n a l y s i s i s per h a p s the f a c t that the a s s u m p t i o n of a complete set of e i g e n v e c t o r s m a y not be s a t i s f i e d . A n advantage i s that the e f f e c t of v a r i o u s a p p r o x i m a t i o n s to the boundary c o n d i t i o n s c a n be i n v e s t i g a t e d n u m e r i c a l l y . T h e s e e f f e c t s have a l s o been i n v e s t i g a t e d by V a r a h f o r the s c a l a r d i f f u s i o n equation, (1970, 1971), and f o r p a r a b o l i c s y s t e m s , (1971a). F o r a p p l i c a t i o n s of the s t a b i l i t y t h e o r y of V a r a h see G o t t l i e b and G u s t a f s s o n (1975). E a r l i e r w o r k on the s t a b i l i t y of p a r a b o l i c s y s t e m s w i t h boundary c o n d i t i o n s was done by Osher (1968). The s u f f i c i e n t c o n d i t i o n s f o r s t a b i l i t y of one step methods, as g i v e n by V a r a h , a r e f i r s t that the scheme m u s t be stabl e when a p p l i e d to a pure i n i t i a l v a l u e p r o b l e m . (Thus the w o r k of W i d l u n d c a n be u s e d to c h e c k t h i s c o n d i t i o n . ) F u r t h e r the a m p l i f i c a t i o n m a t r i x Q = [ Q Q ] °f the d i f f e r e n c e scheme Q g u j ^ = Q_^u^ ^ sh o u l d have no e i g e n v a l u e s z w i t h |z| > 1, and i n a d d i t i o n , some c o n d i t i o n on the s p e c t r u m of Q near z = 1 m u s t be s a t i s f i e d . It i s a l s o shown by V a r a h how one c a n c h a r a c t e r i z e e i g e n v a l u e s z w i t h |z| > 1. E s s e n t i a l l y the same c h a r a c t e r i z a t i o n i s obt a i n e d when a n a l y z i n g the c a s e w h ere the eigenvalue p r o b l e m L ^ v ^ = ^ n ^ - n v n °^ t h i s s e c t i o n has e i g e n -v a l u e s z w i t h Re(z) > 0. T h i s w i l l be c o n s i d e r e d a g a i n i n E x a m p l e (6.3) and (6.4). It i s not c l e a r how the c o n d i t i o n on the s p e c t r u m of Q r e l a t e s to the c u r r e n t e i g e n s y s t e m a n a l y s i s . The p u r pose of the next e x a m p l e i s to i n d i c a t e why i t i s d e s i r a b l e to r e s t r i c t to d i s c r e t i z a t i o n s i n t i m e w h i c h have a s t a b i l i t y r e g i o n that i n c l u d e s at l e a s t the e n t i r e negative r e a l a x i s . E x a m p l e 6.3. It f o l l o w s f r o m the g e n e r a l t h e o r y of K r e i s s (1972), that e a c h eigenvalue of the eige n v a l u e p r o b l e m (6.6) L y = y " = Xy , y(0) = y'(l) = 0 , -86-is approximated by an eigenvalue of the discrete eigenvalue problem ( 6 - 7 ) Vh - W h . as h —- 0, provided that L^u-^ = K^f^ defines a consistent and stable finite difference approximation to L,y=f. (See Section (I. 7) J Since the eigen-values of (6.6) are given by \ f c = - L"5 1 ' k - *» ^ f o l l o w s t h a t there are eigenvalues X.^  of (6.7) such that \^ -*• -co as h — 0. Hence the restriction to discretizations in time with an unbounded region T as in 0 a Theorem (6.2) is desirable. If T were bounded then a restriction on the a size of the time step would have to be imposed. For example, consider the case where L, u. = -— (u. ,-2u.+u. ,,) and K. u. = u. , (1 < i < J-1). The J h 2 1 j-1 j j+r h j J - J -discrete boundary conditions are U Q = 0 and ( U j - U j _ ^ ) / h = 0. Then the eigenvalues of (6.7) are easily found to be given by Thus if the region T includes only a part (-(3,0) of the negative real axis a Q 2 then At must satisfy At < -*r h in order to guarantee that At L . e T . ' — 4 to h, x a The stability of particular finite difference approximations w i l l be investigated in the next two examples as well as in numerical examples in the next section. Since one can always apply a strongly A(a)~stable method to the system (6.3), it is necessary only to study the effect of various spatial difference approximations on the eigenvalues of the eigen-value problem (6.7). Specifically, it is important to check that this dis-crete eigenvalue problem does not allow eigenvalues X.^  that tend to 00 in the right half plane and that the condition number K-("v"^ ) is bounded. - 8 7 -Example 6 . 4 . -fAlthough stability and consistency of the spatial difference approximation imply that the eigenvalues of (6.6) are approximated by those of the discrete eigenvalue problem ( 6 . 7 ) , it need not be true that a l l eigenvalues of ( 6 . 7 ) converge to eigenvalues of ( 6 . 6 ) . Consider for example the case where L^ and are as in the previous example, but where the boundary condition at x = 1 is approxi-mated by (2-c) h = 0 This approximation is consistent with y (1) = 0 , provided that c + 2 . The 4 order of consistency is equal to one. If c =-j then the order of consistency is equal to two. (See also Example (1 . 8.4').) For this discretization the operator in ( 6 . 3 ) is the identity and the operator L ^ can be written as = P j j l ^ , where " - 2 1 P h = 2-c and 1 - 2 1 i ~ Let v^ = v^. Then the eigenvalue problem T^v^ = X^K^v^ becomes JL ~ JL ~ ~ P£ -J-hPj^h = ^h vh' w ^ c ^ admits a complete set of orthonormal eigenvectors with real eigenvalues. The matrices and containing the eigenvectors ~ -1 ~ v^ ^ and v^ ^ respectively by column are then related by = P ^ h " Hence K ( V h ) " V h 2 V h 2 - " 'h " 2 11 "h " 2 ^2 I U II 7 < m a x { | c - 2 | , T J - T } , JL )2 h " 2 c - 2 -88 -provided c + 2. Th i s shows that the last condi t ion of T h e o r e m (6.1) and T h e o r e m (6.2) i s sa t is f ied for this d i s c r e t i za t i on i n space when c + 2. Next consider the e igenvalues . App l i ca t i on of the G e r s g o r i n theorem, (Lancaster (1969), p . 226), shows that a l l eigenvalues are negative i f and only i f c < 2. If c = 2 then there is an eigenvalue at ze ro and i f c > 0 then there is a s t r i c t l y posi t ive eigenvalue. So assume that c > 2, X.^ > 0 and let v^ be the eigenvector assoc ia ted wi th X.^. Thus v^ sat isf ies v . , , - (2 + h 2X.)v. + v . , = 0, 1 < j < J - 1 . j+1 J J - 1 - -The solut ion to these difference equations is v . = c . z i + c_z,^ , 0 < j < J . j l h 2 h ' — — where z. and 1/ z, are the roots of h h (6.8) z 2 - (2 + h 2 X ) z + 1 = 0 But v^ must a lso sat isfy the d iscre te boundary condi t ions . In pa r t i cu l a r , the condit ion V Q = 0 i m p l i e s that v~ = C j ( z ^ + z ^ 1 ) , (0 < j < J), a n d the boundary condit ion at x = 1 leads to the equation (6.9) ( z h ~ z h J ) " c ( z h " 1 _ z h J + 1 ) + ( c - l ) ( z ^ " 2 - z ^ J + 2 ) = 0. Cons ider now the equation obtained f r o m (6.9) by omit t ing a l l negative powers , v i z . z. 2 - cz , + c-1 = 0. h h Th i s is the cha rac t e r i s t i c equation of the boundary condi t ion at x = 1. Its roots are z^ = 1 and z 2 = c - 1 . Since by assumption c-1 > 1 i t fol lows that there is a root z ^ of (6.9) such that z ^ — c - i as h -*• 0. The - 8 9 -c or responding eigenvalue i s obtained f r o m ( 6 . 8 ) and a s y m p t o t i c a l l y given by h h 2 z h h 2 ( c - l ) Thus X^ —> + oo as h —» 0. T h i s eigenvalue computation i s e s s e n t i a l l y equivalent to the technique of V a r a h when a p p l i e d to the c u r r e n t d i s -c r e t i z a t i o n (see V a r a h (1970), p. 33, and V a r a h (1971).) T h e r e the effect of the d i s c r e t e boundary conditions i s r e l a t e d d i r e c t l y to the eigenvalues of the a m p l i f i c a t i o n m a t r i x , whereas here these effects are r e l a t e d to the eigenvalues of K ^ L ^ . If the d i s c r e t i z a t i o n i s completed by applying a i m u l t i s t e p method to the s y s t e m K^u^ = J - ^ ^ > a n d i f this m u l t i s t e p method i s strongly A(a)-stable, but the positive r e a l axis i s not contained i n its s t a b i l i t y region, then the r e s u l t i n g scheme i s stable i f and only i f c < 2. The s i m p l e s t example of such a m u l t i s t e p method is the t r a p e z o i d a l r u l e . (See E x a m p l e (1.4. 1).) F o r this method the s t a b i l i t y r e g i o n i s the entire negative half - plane. That this scheme i s unstable for c > 2 can be checked by applying it to the i n i t i a l data g^ = v^, where v ^ i s the eigen-vector a s s o c i a t e d with the positive eigenvalue T h e r e are d i s c r e t i z a t i o n s i n time for which the complete difference scheme w i l l be stable f o r s m a l l enough h, even i f c > 2. In fact, a l l methods that are s t r i c t l y stable at co have this p r o p e r t y . E x ample 6.5. A g a i n c o n s i d e r the eigenvalue p r o b l e m ( 6 . 6 ) and its d i s c r e t e v e r s i o n (6.7) with L, v. = -4- (v. , - 2v. + v..,) n J h 2 3 3 3 -90-and K, v . = -r— v. - + -pr v . + -pr v . , , h j 1 2 j - 1 1 2 j 1 2 j + 1 Thus this eigenvalue p r o b l e m a r i s e s i n the s tabi l i ty analys is when using the fourth order spat ia l difference approximat ion d i scussed i n Example ( 1 . 4 . 4 ) . Let the approximat ion to the boundary condit ion y ' ( l ) = 0 have the f o r m ( 6 . 1 0 ) 2 b i u J + i = 0 1 = _ r J It i s now diff icult to check ^ (V^) , so only the eigenvalues w i l l be cons idered . Proceeding as i n the previous example , i t is found that the eigenvalues are given by ( 6 . 1 1 ) x. = 2 1 2 2 ( Z - 1 ) 2  1 1 l l ( z * + 1 0 Z + 1 ) for each z which is a root of 0 ( 6 . 1 2 ) 2 " z " ( J + 1 ) ) = 0 . i = - P j Note that i f z is a root of ( 6 . 1 2 ) then so i s \ . The cha rac te r i s t i c po ly-nomia l of the discre te boundary condit ion ( 6 . 1 0 ) i s given by cj(z) = . Z b. z J 1 = _ r J Now equation ( 6 . 1 2 ) has a root z sat isfying | z | > i + e , J l a rge , i f and only i f the cha rac te r i s t i c equation C j ( z ) = 0 has a root |z[ > 1 . F u r t h e r , i f z i s a root of ( 6 . 1 2 ) w i th | z | = 1 , then an easy computation shows that the c o r r e s p o n d i n g eigenvalue i s r e a l and ne g a t i v e . T hus i t i s s u f f i c i e n t to c h e ck that none of the r o o t s of the c h a r a c t e r i s t i c e q u a t i o n of the d i s -c r e t e boundary c o n d i t i o n exceeds one i n a b s o l u t e v a l u e . A s i n the p r e v i o u s e x a m p l e the s e c o n d o r d e r a p p r o x i m a t i o n h ( 2 U J " 2 u J - l + 2 u J - 2 ) = ° s a t i s f i e s the ei g e n v a l u e c o n d i t i o n because the r o o t s of i t s c h a r a c t e r i s t i c e q u ation a r e = 1 and z^ S i m i l a r l y , the c h a r a c t e r i s t i c e q u a t i o n of the t h i r d o r d e r a p p r o x i m a t i o n 6 S < l l u J - 1 8 u J - l . + 9 u J - 2 " 2 u J - 3 } = ° has r o o t s z ^  = 1, z^ = ^ z 3 = 2 2 — w r f c h \z^ \= | z | = 0. 426, The c h a r a c t e r i s t i c e q uation of the f o u r t h o r d e r a p p r o x i m a t i o n T 2 h ( 2 5 u J " 4 8 u J - l + 3 6 u J - 2 - l 6 u J - 3 + 3 u J - 4 ) = 0 has r o o t s zy = 1, z 2 = 0.381, z^ = 0.269 + 0.492 i and z^ = 0.269 - 0.492 w i t h \z^\ = I z^ I = 0.561. T h e r e f o r e , these d i f f e r e n c e equations a l s o s a t i s f y the eigenvalue c o n d i t i o n . 7. N U M E R I C A L E X A M P L E S E x a m p l e 7.1. The purpose of t h i s e x a m p l e i s to i n d i c a t e the r e l a t i v e a c c u r a c y of v a r i o u s f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s . F o r t h i s purpose c o n s i d e r the s i m p l e p r o b l e m = ^ x x + ( 1 + """2)y ' 0 < x < 1 > ° < t < 3 , w i t h i n i t i a l c o n d i t i o n -92-y(x,0) = s i n i r x , 0 < x < 1 , and boundary c o n d i t i o n s y(0,t) = y ( l , t ) = 0 , 0 < t < 3 . The s o l u t i o n i s y(x,t) = e ^ s i n f i r x ) . The m a x i m u m value t h i s s o l u t i o n 3 1 r e a c h e s i s y{j,3) = e = 20. U n i f o r m m e s hes a r e u s e d w i t h h = — and A t = j | . F i n i t e d i f f e r e n c e a p p r o x i m a t i o n s have the f o r m 0 u l=-p 1-^ =0 v n+i n u. = 0 h w h e r e V j . =.2 d i V i i=-r J and m K,u! = ). e. u n ( z . . ) . h J i = l 1 J ' 1 The c o l l o c a t i o n p o ints z. . c o i n c i d e w i t h the m e s h p o i n t s . If m = 1 then z. . = x. and i f m = 3 then z. . = x. _,. . (1 < i < 3). F o r the g i v e n J , i J 3,1 J-2+i - -p r o b l e m the c h o i c e m = 3 c o r r e s p o n d s to the " M e h r s t e l l e n v e r f a h r e n " of C o l l a t z . (See equation (1. 1. 7) and E x a m p l e (1.4.4).) R e s u l t s appear i n T a b l e (7.1). The s p a t i a l d i f f e r e n c e a p p r o x i m a t i o n i n e x p e r i m e n t 5 i s not compact. Thus near the b o u n d a r i e s other d i f f e r e n c e -93-T A B L E 7. 1 # r s m o P c o J = 5 J = 10 J = 20 J=40 N 1 1 i 1 2 1 C 2 + 2 .392 ^ .772+1 .279 + 1 .170+1 6 .320 + 2 .585+1 .149 + 1 .517° 15 . 311 .560+1 . 132+1 .360° 30 . 127 + 1 .321° 60 .311° 120 2 1 1 1 2 2 G 2 .599 + 3 .861+2 +2 .523 +^ .459+2 6 .755 + 2 .185+2 .106+2 .888+1 15 .475+2 .108+2 .516+1 .391+1 30 .382+2 .795+1 .306+1 + 1 . 197 + 1 60 .343+2 .669+1 .212+1 .111+1 120 + 2 .325 +^ .609 + 1 .168+1 .701° 240 3 1 1 3 4 1 P 4 .371° . 189" 1 .380" 2 .521" 2 6 .376° .241" 1 . 138" 2 -4 .399 15 .242"1 . 150" 2 -4 .859 30 4 1 1 3 4 4 G 4 .308° .161° .151° .150° 6 .299° .293"1 . 118" 1 . 108"1 15 .223"1 .224" 2 .986" 3 30 .229" 1 .149" 2 . 154" 3 60 , -2 -4 .147 .959 120 -4 .932 240 5 2 2 1 4 4 G 4 + 3 . 120 .158° .151° .150° 6 + 1 .212 .245"1 -1 .135 . 109" 1 15 .243+1 .167"1 .4H" 2 . 112"2 30 . 169" 1 .348" 2 .299" 3 60 .353 - 2 .245" 3 120 .245" 3 240 T A B L E 7.1 -94-(Continued) # r s m o P c o J = 5 J = 10 J = 20 J = 40 N 6 1 1 3 4 2 S 4 .298° .527 c .136" 1 . 148" 1 6 .351° .222"1 .978 - 3 _3 .343 15 .234"1 .143" 2 .631" 4 30 .148" 2 -4 .909 60 7 1 1 3 4 1 P 6 . 377° .243" 1 .152"2 . 108" 3 6 .377° .242"1 -2 . 151 -4 .947 4 15 _4 .944 30 8 1 1 3 4 6 G 6 - 1 .351 .536" 2 .342" 2 .330" 2 6 .238° . 156" 1 .H6" 2 .268"3 15 .307° . 198" 1 .124"2 -4 .824 * 30 -4 .858 60 9 1 1 3 4 4 S 6 .186° . 116" 1 .334"3 .370" 3 6 .300° . 194"1 .120" 2 -4 .719 15 .218"1 .136" 2 _4 .848 30 10 1 1 3 4 1 P 8 . 377° .242"1 -2 .151 .663" 4 6 . 377° .242"1 . 151" 2 . 109" 3 15 -3 . 101 30 -95-f o r m u l a s must be used. The a p p r o x i m a t i o n at x = x^ i s defined by the values r = 1, s = 4 and m = 1, while f o r the approximation at x = Xj. ^ these values a r e r = 4, s = 1 and m = 1. The o r d e r of c o n s i s t e n c y of these f o r m u l a s i s equal to four, as i s the o r d e r of the m a i n finite d i f f e r e n c e a pproximation. Stab i l i t y of this p a r t i c u l a r d i s c r e t i z a t i o n w i l l be c o n s i d e r e d i n Example (7.5). The d i s c r e t i z a t i o n i n time i s defined by the number of steps "p", the o r d e r of c o n s i s t e n c y "o", and a code " c " , indica t i n g the type of method used. "C" i s the t r a p e z o i d a l r u l e . (See E xample (1.4. 1) and Example (3.1).) T h i s approximation i s u s u a l l y r e f e r r e d to as a C r a n k - N i c o l s o n d i s c r e t i z a t i o n i n t i m e . I "G" denotes a method of G e a r . The coe f f i c i e n t s a of these v methods can be found i n G e a r (1971), p. 217. (See a l s o , Example (1.4.2) and E x a m p l e (3.2).) "P" stands f o r a Pade method, based upon a diagonal entry i n the Pade table of r a t i o n a l approximations to the exponential function. The coef f i c i e n t s of this d i s c r e t i z a t i o n i n time may be obtained f r o m V a r g a (1961). (See E xample (3.3) for a fo u r t h o r d e r Pade' formula.) F i n a l l y , "S" denotes a s o - c a l l e d second d e r i v a t i v e method. T h e s e methods have the f o r m (7.1) with u. = 2 and can be found i n E n r i g h t (1974). I The coe f f i c i e n t s a^ a r e d e t e r m i n e d f r o m c o n s i s t e n c y and A ( a ) - s t a b i l i t y r e q u i r e m e n t s . E x ample 7.2. In this example u n i f o r m as w e l l as non-uniform meshes are employed to solve the equation y t = f x u ^ - x ( l - x ) u x + x u , 0 < x < 1 , t > 0 , wi th i n i t i a l c o n d i t i o n -96-u(x ,0 ) = s i n TTX , 0 < x < 1 , and boundary condi t ions u(0,t ) = u ( l , t ) = 0 . T h e d i s c r e t i z a t i o n i n space c o r r e s p o n d s to the s e c o n d c a s e d i s c u s s e d i n E x a m p l e (1.4.3). If the m e s h is u n i f o r m th is f o r m u l a r e d u c e s to the s t a n d a r d s e c o n d o r d e r t h r e e - p o i n t f o r m u l a . If the m e s h i s not u n i f o r m , ( h j + r h ) then the c o l l o c a t i o n point z . , = x . + ~ J - does not c o i n c i d e w i th x . . T o c i r c u m v e n t th is d i f f i c u l t y , the o p e r a t o r is r e d e f i n e d i n t e r m s of the L a g r a n g i a n i n t e r p o l a t i o n f o r m u l a . M o r e p r e c i s e l y , is g i v e n by K h u . = w°(zLi)M._1 + w 1 ( z j i l ) u . + w 2 ( z j j l ) u . + 1 w h e r e the L a g r a n g i a n b a s i s funct ions w 1 (x) a r e de f ined by 0 (x -x )(x-x. + 1 ) 1 (x -x . t ) ( x - x . + 1 ) w ( x ) = h . i . + h . ' ) » w ( x ) = — h ! h . + i — a n d 2 (x-x. t ) (x-x.) w (x) = -r y~— r*- . T h i s y i e l d s the e x p r e s s i o n ( h i - h i + 1 ) ( 2 h i + 1 + h i ) ( 2 h . + h . + 1 ) (h . + 2 h . + 1 ) (h. + 1 - h . ) ( 2 h . + h . + 1 ) V j = 9 h . ( h j + h j + 1 ) U j - 1 9h. h j + 1 V 9h. + 1 ( h . + h . + 1) V l i T h e r e s u l t i n g s p a t i a l d i f f e r e n c e a p p r o x i m a t i o n K ^ u ^ = L ^ u ^ r e m a i n s s e c o n d o r d e r . T h e d i s c r e t i z a t i o n is c o m p l e t e d by a p p l y i n g the t r a p e z o i d a l r u l e i n t i m e . T h e c o m p l e t e d i f f e r e n c e s c h e m e is then -97-, n n- i " (u, - u, ) . , xr h h 1 T , n n-1. K h At = 2 + U h } ' 0 . , . u. = s i n ( i r x . ) . J J Thus this scheme i s a g e n e r a l i z a t i o n of the well-known C r a n k - N i c o l s o n scheme to non-uniform meshes. T h i s g e n e r a l i z a t i o n i s by no means unique. Other second o r d e r f o r m u l a s for non-uniform meshes that use m o r e than one evaluation of each c o e f f i c i e n t function may be d e r i v e d a l s o . F i g u r e s (7. l a , b, c, d) show the approximate solution i n the case 1 1 of a u n i f o r m mesh, with h. = 47J and At = — . P l o t t e d are the i n i t i a l function along with the approximate solution at v a r i o u s t i m e s . The approximate solution i s i n t e r p o l a t e d with a cubic s p l i n e . R e s u l t s of a computation on a non-uniform m e s h ar e given i n F i g u r e s (7. 2a, b, c, d, and e). Now there a r e only 31 meshpoints d i s t r i b u t e d as i n d i c a t e d . The t i m e s t e p i s as before. C l e a r l y , this computation i s able to follow the shock p r o f i l e f o r a longer time than the computation on the u n i f o r m mesh. Example 7.3. In this example the s t a b i l i t y of a number of finite d i f f e r e n c e schemes f o r solving the d i f f u s i o n equation y t(x,t) = Ly(x,t) = y n ( x , t ) , with y(x, 0) = g(x) and y(0,t) = y ( l , t ) = 0, i s t e s t e d n u m e r i c a l l y . F r o m T h e o r e m ( 6 . 2) it follows that for this purpose one has to investigate the behavior of the eigenvalues X., of L. v. = l K , v , , A l s o , the condition & h h h h h h ' number KCV^) must be bounded. In T a b l e (7. 2) a numbei of d i s c r e t i z a t i o n s F i g u r e 7. l a -99-0*91 U d 0'91 EES"EI Z.99'07 Q'8 s i x d - n E E E ' S Z.99'2 O ' D - 100-Q'QZ O'ST D'Ol O'S o L D ca u •d SuO •i—i j i i -i • i • i • * in Ln (M J I it Ji J I o " a LTS2 o'oe i 0'S7 9IXH-H O ' O l O'S O'O F i g u r e 7.Id CD. <n m CO r~ L D CO 01 IT) I m r— T x x x x .0 if it X X X X X X X X 0.125 0.25 0.375 X ¥ X X X ¥ ¥ X in_l i 0.5 X - f l X I 5 0.625 • _, in T -4.0 cn y - p y F i g u r e 7.2d 0(.0 — 1 5 . 0 30,0 _ J u - n x i s 45.0 60.0 _J 75.0 90.0 _J 4 O Q CD CD J10-JL 45.0 EO.Q 75.0 H II 00 90.0 -901--107-T A B L E 7. 2 j = 1 j = 2 3 < j < J - 3 # r s m o r s m o r s m o 1 1 2 1 2 2 2 1 4 2 2 1 4 2 1 4 1 4 •2 2 1 4 2 2 1 4 3 1 4 1 4 2 3 1 4 3 3 1 6 4 1 3 4 6 2 2 5 8 2 2 5 8 i n space of the f o r m (6. 2) i s given, for which the eigenvalues X^ were observed to be r e a l with < 0. In fact, a l l eigenvalues converged to those of the continuous p r o b l e m Lw(x) = X w(x), w(0) = w(l) = 0. T h e s e a l g e b r a i c eigenvalue computations were p e r f o r m e d with h = ^ , . A l s o computed was K(V^) and as h -*• 0 nop si g n i f i c a n t i n c r e a s e i n its value was observed. Thus the d i s c r e t i z a t i o n s defined i n .Table (7.2) can be completed by applying any s t r o n g l y A( a)-stable multistep method. L i s t e d i n this table a r e the main di f f e r e n c e equation at x^, (3 < j < J-3), and the equations at x^ and which may d i f f e r f r o m the m a i n approximation. The equations at X j and X j ^ a r e the " r e f l e c t i o n " of those at x^ and x^ r e s p e c t i v e l y . If m = 1 the c o l l o c a t i o n point i s Z j = X j . In experiment 4 the c o l l o c a t i o n points for the approximation at x^ a r e z^ = x^ ^, (1 < i < 4). The o r d e r of c o n s i s t e n c y i s given under "o". The lower o r d e r equations near the boundaries do not affect the o v e r a l l a c c u r a c y of the schemes. Example 7.4. C o n s i d e r again the a l g e b r a i c eigenvalue p r o b l e m of the previous example. L e t the d i s c r e t i z a t i o n i n space be defined as follows. -108-F o r each approximat ion let there be only one co l loca t ion point = x „ (So m = 1.) F o r the m a i n difference approximat ion take r = s = 3. Thus this equation i s a seven point formula of order s i x . The approximat ion at x^ i s defined by the values r = 2 and s = 5. A t x^ these values are r = 1 and s = 6. The difference equations at X j and X j ^ are once more the " ref lec t ion" of those at x^ and x^ . Thus the extra difference equations near the boundaries are also s ixth o rde r . A g a i n the eigenvalues \ ^ were computed for h = ^ , y^- , . F o r each h a l l eigenvalues were r e a l and negative, wi th the exception of two conjugate pa i r s of complex eigenvalues l i s ted below: h = | - : \ = -14. 9 ± 7. 18 i , \ = -20. 5 ± 6. 50 i , h = : \i = -66. 1 ± 24. 8 i , \ = -67. 8 ± 25. 1 i , h = : \ i = -150. 5 ± 55. 8 i , X = -150. 8 ± 56. 0 i . So apparently there i s a double complex root as h —- 0, that i s a p p r o x i -mately given by \ = ~Y ( - ° - 2 6 ± ° - ° 9 6 i ) h Therefore , i f the d i s c r e t i za t i on i s completed by applying a mult is tep method i n t ime, then this method should be s t rongly A{a)-stable wi th tan a > ^6* ' ° r a > 2 - l 2 ° - There are s t rongly A(a)-s table methods that do not satisfy this requirement . (See V a r a h (1975), p . 19). Example 7 . 5 . Cons ider the diffusion equation of Example (7. 3) , but now let the boundary conditions be given by y (0, t) = y(l> t) = 0. The approxi-mat ion to the boundary condition at x = 0 i s assumed to have the f o r m s Q -109-YJ b.w.(t) = 0 w i t h o r d e r of c o n s i s t e n c y e q u a l to s n . In T a b l e (7.3) i=0 1 1 U t h r e e s p a t i a l d i s c r e t i z a t i o n s a r e g i v e n f o r w h i c h the d i s c r e t e e i g e n v a l u e p r o b l e m L ^ v ^ = ^ h ^ n v h w a s o b s e r v e d to have only r e a l and negative 1 1 1 e i g e n v a l u e s . The c o m p u t a t i o n s w e r e p e r f o r m e d f o r h = -g , , -^g. The c o n d i t i o n number K ( V ^ ) was a l s o c o m puted and no i n c r e a s e of s i g n i f i c a n c e was seen f o r s m a l l e r h. W i t h the e x c e p t i o n of the boundary c o n d i t i o n at x = 0 the d i s c r e t i z a t i o n s i n T a b l e (7.3) and T a b l e (7.4) below a r e i d e n t i c a l to those of T a b l e (7.2). T A B L E 7.3 j : = 1 j = 2 3 < j < J-3 # s o r s m o r s m o r s m o i 4 1 2 1 2 2 2 1 4 2 2 1 4 2 4 1 ' 4 1 4 2 2 1 4 2 2 1 4 3 6 1 4 1 4 2 3 1 4 3 3 1 6 Two d i s c r e t i z a t i o n s f o r w h i c h the a l g e b r a i c e i g e n v a l u e p r o b l e m was o b s e r v e d to have c o m p l e x e i g e n v a l u e s a r e g i v e n i n T a b l e (7.4). In e x p e r i m e n t 4 t h e r e a r e two d i s t i n c t conjugate p a i r s of e i g e n v a l u e s a p p r o x i -m a t e l y g i v e n by \ j = - ^ ( - 0 . 178 ± 0.049 i) and \ = ( - 0 . 26l± 0.097 i ) . In h h j e x p e r i m e n t 5 t h e r e i s o n l y one conjugate p a i r g i v e n by \. =—r (-0.246 ± 0 .044 i ) , T he m u l t i s t e p d i f f e r e n c i n g i n t i m e s h o u l d t h e r e f o r e be s t r o n g l y A(a)-stable w i t h a > a . . The v a l u e a . as c o m puted f r o m the c o m p l e x e i g e n -— m i n m m r r ° v a l u e s above i s a l s o g i v e n i n the t a b l e , -110-T A B L E 7.4 j = 1 J = 2 3 < j < J -3 # s o r s m o r s m o r s m o a m i n 4 6 1 6 1 6 2 5 1 6 3 3 1 6 21° 5 8 1 3 4 6 2 2 5 8 2 2 5 8 11° Example 7 .6 , The purpose of this example i s to show the effect of usinjj an imprope r approximat ion to a boundary condi t ion. The equation under considera t ion i s y f c(x, t) = ( J^T ) y x x ( x , t ) , 0 < x < 1 , t > 0 , wi th i n i t i a l condit ion 3TT y(x, 0) = cos(—r- x) , 0 < x < 1 , and boundary conditions y x ( 0 , t ) = y ( l , t) = 0 , t > 0 . -1 3TF The solution i s y(x, t) = e cosC-^- x) . L e t the d i s c r e t i za t i on i n space be given by K ^ w (x^ t) = L ^ w f X j , t) , 1 < j < J-1 2 where K h w . = W j . j + Q W j , + w . + 1 and L h w . = --j(w. _ i-2w.+w j + . -1 11-Th e order of consis tency of this approximat ion i s four. (See Example (I. 4. 4). ) A one parameter f ami ly of approximations to the boundary condit ion at x = 0 i s given by 88w Q - (5a + 144)w 1 + (8a +72)w 2 - (3a +16)w 3 = 0 This approximat ion i s consistent i f a + 24. The order of consis tency i s equal to two, except i f a = 0, when the o rder i s three. The cha rac t e r -i s t i c po lynomia l of the d iscre te boundary condi t ion i s given by c Q ( z ) = 8 8 z 3 - (5a +144)z 2 + (8a + 72)z - (3a + 16). If a > .24 then the cha rac te r i s t i c equation c^fz) = 0 has a r e a l root > 1, F o r such a the eigenvalue p rob lem K ^ v ^ = X ^ L ^ v ^ has a pos i t ive eigen-value that i s a sympto t i ca l ly given by 12(z - l ) 2 h h 2 ( z 2 + l O z j + 1) This may lead to an in s t ab i l i t y . (See Example (6.5).) F o r example, let the d i sc re t i za t ion i n t ime be the fourth order method of Gear . (See Gear (1971), p . 217.) This leads to the difference scheme TTST Kh(25% - *K'1 + 3 K" 2 - 'K " 3 + K*) = V S • 2 Of the r e a l axis only the i n t e rva l [0, 10—] does not l i e i n the s tabi l i ty reg ion of this mul t is tep method i n t ime . (See Gear (1971), p . 216.) Two numer i ca l computations were per formed wi th a = 30 and At = 0.025. The f i r s t computat ion wi th h = 0.050 and the second wi th h = 0 .025. The values of A t X^ for these data are equal to 2.03 and 8.12 r e spec t ive ly . Hence both schemes are unstable. The approximate solution at t ime t = 0. 5, interpolated wi th a cubic spl ine, as w e l l as the exact solut ion are shown i n F igu re s (7. 3a, b). The same computation was repeated wi th a = 20. This i s stable. Resul ts at t = 0. 5 are shown i n F i g u r e s (7 .4a ,b ) . F i g u r e 7.3a •J .0 0.0 U-RXIS 2.0 o > H II r i -ll II o O o o O r o r o o o 5.0 cn ro Ln -1,0 *1 P CD cr L.O 2.0 3.0 4.0 5.0 U-RXIS -o.o CJ (TQ 0) _ST T --9TT--117-B I B L I O G R A P H Y G. B i r k h o f f and S. G u l a t i , O p t i m a l few-point d i s c r e t i z a t i o n s of l i n e a r s o u r c e p r o b l e m s , S I A M J . Numer. A n a l . 11, 1974, pp.700-728. C. de B o o r and B. S w a r t z , C o l l o c a t i o n at G a u s s i a n p o i n t s , S I A M J . Numer. A n a l . 10, 1973, pp. 582-606. J. H. B r a m b l e and B. E. Hubbard, On a f i n i t e d i f f e r e n c e analogue of an e l l i p t i c b o u n d a r y v a l u e p r o b l e m w h i c h i s n e i t h e r d i a g o n a l l y dominant nor of the nonnegative type, J . M a t h , and P h y s . 43, 1964, pp. 117-132. J. C. B u t c h e r , I m p l i c i t R u n g e - K u t t a p r o c e s s e s , Math, of Comp. 18, 1964, pp. 50-64. P. G. C i a r l e t , M. H. S c h u l t z and R. S. V a r g a , N u m e r i c a l methods of hi g h o r d e r a c c u r a c y f o r n o n l i n e a r b o u n d a r y v a l u e p r o b l e m s , N u m er. Math. 9, 1967, pp. 394-430. L. C o l l a t z , N u m e r i s c h e Behandlung von D i f f e r e n t i a l - g l e i c h u n g e n , S p r i n g e r -V e r l a g , B e r l i n , I960. G. D a h l q u i s t , S t a b i l i t y and e r r o r bounds i n the n u m e r i c a l i n t e g r a t i o n of o r d i n a r y d i f f e r e n t i a l equations, K u n g l . Tekn. Hogsk. St o c k h o l m , no. 130, 1959, 87 pp. G. D a h l q u i s t , A s p e c i a l s t a b i l i t y p r o b l e m f o r l i n e a r m u l t i s t e p methods, B I T , 3, 1963, pp. 27-43. W. H. E n r i g h t , Second d e r i v a t i v e m u l t i s t e p methods f o r s t i f f o r d i n a r y d i f f e r e n t i a l equations, S I A M J. Numer. A n a l . 11, 1974, pp. 321-331. C. W. Gear, N u m e r i c a l i n i t i a l v a l u e p r o b l e m s i n o r d i n a r y d i f f e r e n t i a l equations, P r e n t i c e - H a l l , 1971. R. D. G r i g o r i e f f , D i e K o n v e r g e n z des Rand-und E i g e n w e r t p r o b l e m s l i n e a r e r g e w o h n l i c h e r D i f f e r e n z e n g l e i c h u n g e n , N u m e r . Math., 1970, pp. 15-48. D a v i d G o t t l i e b and B e r t i l G u s t a f s s o n , G e n e r a l i z e d Du F o r t - F r a n k e l methods fo r p a r a b o l i c i n i t i a l b o u n d a r y v a l u e p r o b l e m s , I C A S E R e p o r t no. 75-5, NASA's L a n g l e y R e s e a r c h C e n t e r , Hampton, V A . B. Ha k b e r g , U n i f o r m l y m a x i m u m n o r m stab l e d i f f e r e n c e s chemes, B.I.T. 10, 1970, pp. 266-276. B. L. Hulme, D i s c r e t e G a l e r k i n and r e l a t e d one step methods f o r o r d i n a r y d i f f e r e n t i a l equations, R e s e a r c h r e p o r t , Sandia L a b o r a t o r i e s , A l b u q u e r q u e , N. M. , 1971. H. B. K e l l e r , N u m e r i c a l methods f o r i two-point b o u n d a r y v a l u e p r o b l e m s , B l a i s d e l l , London, 1968. H. B. K e l l e r , A c c u r a t e d i f f e r e n c e methods f o r 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 s y s t e m s s u b j e c t to l i n e a r c o n s t r a i n t s , S L A M J . Numer. A n a l . , 6, 1969, pp. 8-30. -118-H. B. K e l l e r , N u m e r i c a l s o l u t i o n of boundary v a l u e p r o b l e m s f o r o r d i n a r y d i f f e r e n t i a l equations, A. K. A z i z , ed. , A c a d e m i c P r e s s , 1975, pp. 27-88. H. O. K r e i s s , D i f f e r e n c e a p p r o x i m a t i o n s f o r b o u n d a r y and e i g e n v a l u e p r o b l e m s f o r o r d i n a r y d i f f e r e n t i a l equations, Math. o f C o m p . 26, 1972, pp. 605-624. P. L a n c a s t e r , T h e o r y of m a t r i c e s , A c a d e m i c P r e s s , 1969. S. N o r d m a r k , U n i f o r m s t a b i l i t y of a c l a s s of p a r a b o l i c d i f f e r e n c e o p e r a t o r s , B .I.T. 14, 1974, pp. 314-325. S. J. Osher, M a x i m u m n o r m s t a b i l i t y f or p a r a b o l i c d i f f e r e n c e s c h e m e s i n h a l f - s p a c e , H y p e r b o l i c equations and waves, B a t e l l e I n s t i t u t e , Seattle, Washington, 1968, pp. 61-75. H. Pade, Sur l a r e p r e s e n t a t i o n a p p r o c h e e d'une f o n c t i o n p a r des f r a c t i o n s r a t i o n e l l e s , t h e s i s , Ann. de L ' E c . Nor., (3), 9, 1892. I. G. P e t r o v s k i i , U b e r das C a u c h y s c h e P r o b l e m f u r e i n S y s t e m l i n e a r e r p a r t i e l l e r D i f f e r e n t i a l g l e i c h u n g e n i m Gebiete ^der N i c h t a n a l y t i s c h e n F u n c t i o n e n , B u l l . U n i v . d'Etat M o s c o u , 1A, No. 7, 1938, pp. 1-74. M. R e i m e r , F i n i t e d i f f e r e n c e f o r m s c o n t a i n i n g d e r i v a t i v e s of h i g h e r o r d e r , S I A M J . Numer. A n a l . 5, 1968, pp. 725-738. R. D. R i c h t m e y e r and K. W, M o r t o n , D i f f e r e n c e methods f o r i n i t i a l v a l u e p r o b l e m s , I n t e r s c i e n c e P u b l i s h e r s , 1967. R. D. R u s s e l l and L. F. Shampine, A c o l l o c a t i o n method f o r boundary v a l u e p r o b l e m s , N u mer. Math. 19, 1972, pp. 1-28. R. D. R u s s e l l and J . M. V a r a h , A c o m p a r i s o n of g l o b a l methods f o r l i n e a r two-point boundary v a l u e p r o b l e m s , Math, of Comp. 29, 1975, pp. 1-13. M. H. S c h u l t z , S p l i n e A n a l y s i s , P r e n t i c e - H a l l , 1973. J . N. Shoosmith, A h i g h e r o r d e r f i n i t e d i f f e r e n c e method f o r the s o l u t i o n of two-point boundary v a l u e p r o b l e m s on a u n i f o r m mesh, N u m e r i c a l s o l u t i o n of b o u n d a r y v a l u e p r o b l e m s f o r o r d i n a r y d i f f e r e n t i a l equations, A. K. A z i z , ed. , A c a d e m i c P r e s s , 1975, pp. 355-369. G. S t r a n g and G. J. F i x , A n a n a l y s i s of the f i n i t e element method, P r e n t i c e - H a l l , 1973. B. K. S w a r t z , The c o n s t r u c t i o n and c o m p a r i s o n of f i n i t e d i f f e r e n c e analogs of some f i n i t e e l ement schemes, r e p o r t , L o s A l a m o s S c i e n t i f i c L a b o r a t o r y , L o s A l a m o s , N . M. , 1974. J. M. V a r a h , M a x i m u m n o r m s t a b i l i t y of f i n i t e d i f f e r e n c e a p p r o x i m a t i o n s to the m i x e d i n i t i a l b o undary v a l u e p r o b l e m f o r the heat equation, Math, of Comp. 24, 1970, pp. 31-44. J . M. V a r a h , S t a b i l i t y of h i g h o r d e r a c c u r a t e d i f f e r e n c e methods f o r p a r a b o l i c equations w i t h boundary c o n d i t i o n s , S I A M J . Numer. A n a l . 8, 1971, pp. 598-615. - 1 1 9 -J . M. V a r a h , S t a b i l i t y of h i g h e r o r d e r a c c u r a t e d i f f e r e n c e methods f o r p a r a b o l i c equations w i t h boundary conditions, S I A M J . Numer. A n a l . 8, 1971a, pp. 569-574. J. M. V a r a h , S t i f f l y s t a b l e l i n e a r m u l t i s t e p methods of extended o r d e r , T e c h n i c a l R e p o r t , Dept.of Comp. Sc., U n i v e r s i t y of B r i t i s h C o l u m b i a , 1975. R. S. V a r g a , On h i g h e r o r d e r sta b l e i m p l i c i t methods f o r s o l v i n g p a r a b o l i c p a r t i a l d i f f e r e n t i a l equations, J . M a t h and P h y s . , 40, 1961, pp. 220-231. R. S. V a r g a , H e r m i t e i n t e r p o l a t i o n type R i t z methods f o r two-point b o u n d a r y v a l u e p r o b l e m s , N u m e r i c a l s o l u t i o n of p a r t i a l d i f f e r e n t i a l equations, J . H. B r a m b l e , ed. , A c a d e m i c P r e s s , 1966, pp. 365-373. R. W e i s s , The a p p l i c a t i o n of i m p l i c i t R u n g e - K u t t a and c o l l o c a t i o n methods to b o u n d a r y v a l u e p r o b l e m s , Math, of Comp. 28, 1974, pp. 449-464. O. B. Widlund, On the s t a b i l i t y of p a r a b o l i c d i f f e r e n c e schemes, M a t h , of Comp. 19, 1965, pp. 1-13. O. B. . Widlund, S t a b i l i t y of p a r a b o l i c d i f f e r e n c e s chemes i n the m a x i m u m norm, N u m e r . Math. 8, 1966, pp. 186-202. K. Wright, Some r e l a t i o n s h i p s between i m p l i c i t Runge-Kutta, c o l l o c a t i o n and L a n c z o s T methods, and t h e i r s t a b i l i t y p r o p e r t i e s , B . I.T. 10, 1970, pp. 217-227. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080139/manifest

Comment

Related Items