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

DIFFERENCE METHODS FOR ORDINARY D I F F E R E N T I A L  EQUATIONS  WITH APPLICATIONS TO P A R A B O L I C  EQUATIONS  by E u s e b i u s 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 , 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 ,  1971 1973  A T h e s i s S u b m i t t e d i n P a r t i a l F u l f i l m e n t of the R e q u i r e m e n t s for the Degree of D o c t o r of P h i l o s o p h y i n the  Institute of  Applied Mathematics  We accept this t h e s i s as c o n f o r m i n g to the required  standard  T H E UNIVERSITY OF BRITISH C O L U M B I A January,  1976  In  presenting  an  advanced degree  the I  Library  further  for  this  shall  agree  thesis  in  at  University  the  make  that  it  freely  permission  s c h o l a r l y p u r p o s e s may  by  his  of  this  representatives. thesis  for  partial  financial  is  for  The  of  University  t  of  gain  ^l/^/jC  British  2075 W e s b r o o k P l a c e V a n c o u v e r , Canada V6T 1W5  of  Columbia,  British  by  shall  Columbia  f  for  extensive the  understood  written permission.  Department  of  available  be g r a n t e d  It  fulfilment  requirements  reference copying of  Head o f  that  not  the  I  agree  and  be a l l o w e d  that  study.  this  thesis  my D e p a r t m e n t  copying or  for  or  publication  without  my  ABSTRACT  The of  first  c h a p t e r o f the  finite difference  nth  order  thesis  approximations  ordinary differential  collocation procedure  to  a method  of u n d e t e r m i n e d  of  these finite  minants states  certain normalization factor compact  difference  The  new  difference  conditions Kreiss based also  are  i s applied  expressed  given.  treated i n detail.  to i n v e s t i g a t e  The  basic  the d e t e r theorem  equations with only be  improved  by  of known,  as  Approximations stability  number  one  suitable well  as  to b o u n d a r y  t h e o r y of H.  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 A  a  T h i s i s the c a s e f o r  Several examples  are  A  as  p r o v i d e d o n l y that  for difference  upon these a p p r o x i m a t i o n s .  O. schemes  of n u m e r i c a l e x a m p l e s  are  given.  In the f i r s t linear  the s e c o n d chapter  first  'included and is  be  o r d e r of c o n s i s t e n c y m a y  approximations also  can  a r e consistent,  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 .  based  which i s equivalent  s m a l l dimension.  and  in linear  It i s s h o w n that the c o e f f i c i e n t s  d o e s not v a n i s h .  equations  c o l l o c a t i o n point.  with polynomials,  approximations  that t h e s e a p p r o x i m a t i o n s  construction  This construction is  coefficients.  of m a t r i c e s of r e l a t i v e l y  w i t h the  to b o u n d a r y v a l u e p r o b l e m s  equations.  upon a l o c a l  difference  is concerned  can  order  c h a p t e r i t i s s h o w n how be  extended  to i n i t i a l v a l u e p r o b l e m s  ordinary differential  the w e l l - k n o w n  the c o n s t r u c t i o n m e t h o d  stability  equations.  Specific  of  for systems examples  theory f o r these d i f f e r e n c e  of  are  equations  summarized.  It i s t h e n s h o w n how linear first  parabolic discretizing  partial  these d i f f e r e n c e methods  differential  i n space  by  a  e q u a t i o n s i n one  suitable method  from  may space  be  applied  variable  the f i r s t  to  after  chapter.  -111-  The  stability  investigated  of  using  e f f e c t of v a r i o u s The  such difference an  approximations  Numerical  for parabolic  eigenvalue-eigenvector  r e l a t i o n of t h i s a n a l y s i s  indicated.  schemes  to the  to the  examples  are  analysis.  boundary  stability also  equations i s In p a r t i c u l a r ,  conditions  theory  included.  of  J.  is M.  the  considered. V a r a h is  -ivTABLE  OF  CONTENTS  ABSTRACT TABLE  i  OF CONTENTS  iv  ACKNOWLEDGMENTS  v  Chapter I.  II.  BOUNDARY VALUE  PROBLEMS  1  1.  Introduction  1  2.  The C o n s t r u c t i o n of F i n i t e Difference A p p r o x i m a t i o n s  6  3.  Effect of the C h o i c e of C o l l o c a t i o n Points  19  4.  E x a m p l e s of F i n i t e Difference A p p r o x i m a t i o n s  21  5.  F i n i t e Difference A p p r o x i m a t i o n s to B o u n d a r y Conditions  29  6.  E x a m p l e s of A p p r o x i m a t i o n s to B o u n d a r y Conditions  34  7.  The S t a b i l i t y of F i n i t e Difference Schemes  37  8.  Numerical Examples  46  INITIAL  VALUE  PROBLEMS  1.  Introduction  2.  Initial Value Problems  54 54 in O r d i n a r y D i f f e r e n t i a l  Equations  55  3.  E x a m p l e s of F i n i t e Difference A p p r o x i m a t i o n s  58  4.  The S t a b i l i t y of F i n i t e Difference A p p r o x i m a t i o n s to Initial Value Problems  65  5.  Initial Boundary Value Problems  71  6.  The S t a b i l i t y of F i n i t e Difference 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 Value P r o b l e m s  76  Numerical Examples  91  7.  BIBLIOGRAPHY  117  -v-  ACKNOWLEDGMENTS  The  author  and  guiding this  her  patience.  w i s h e s to thank  research. Finally,  the  J . M.  Varah for stimulating  H e i s a l s o t h a n k f u l to h i s w i f e  he i s i n d e b t e d  e x c e l l e n t t y p i n g of the m a n u s c r i p t proofreading  Professor  manuscript.  Adrienne for  to M r s . V i v i a n Davies  a n d to D r . F a u s t o  f o r her  Milinazzofor  -1-  Chapter  BOUNDARY  1.  VALUE  PROBLEMS  INTRODUCTION  In to  I  this  chapter  the l i n e a r  (1.1)  the c o n s t r u c t i o n  differential  Ly(x) = y  ( n )  (x)  of f i n i t e d i f f e r e n c e  approximations  equation  +  n  j  1  a (x)y k  ( k )  ( x ) = f(x), 0 < x <  1 ,  k=0 is  investigated.  D e f i n e a m e s h S. = {x. : 0 = x „ < x, < ••• < x h l 0 1 J  T  L  spacings  h. = x. - x. ^ , (1 < j < J ) .  assumed  that m i n  on  S, h  J  l e t w. = J Finite  > J  > -  1  L e t h = max  = 1 } with J  h..  mesh-  It w i l l a l w a y s  be  j  h f o r some  c o n s t a n t K.  F o r a n y f u n c t i o n w(x)  w(x.). J  difference approximations  the p o i n t x = x^ a r e a s s u m e d  to the d i f f e r e n t i a l  e q u a t i o n (1.1) a t  to h a v e the f o r m  s.  J  (1. 2)  L u.^ h  I i=-r.  <L . u  j +  . = f., k  < j< J - n + k  Q  Q  , k  > 0 ,  Q  J  where  f. i s s o m e  approximation  to f..  J precisely  Thus  J - n + 1 such  equations.  J approximation i s said  The  m e s h S,  J T h e w i d t h co.. of e a c h  to be c o m p a c t .  f o r s o m e i n t e g e r co max b  local  J  truncation  error  there a r e  Im-  a p p r o x i m a t i o n i s d e f i n e d a s co. s r . + s. + 1.  co. < co j — max  for each  J  finite  difference  If co. = n + 1 t h e n the  J The assumption  that d o e s  i s made  that  n o t d e p e n d o n S, . h r  i s g i v e n b y T . = L ^ y ^ - f..  If t h e r e i s a  constant C  n of h such that |T.|<Ch s with n > 0, J s  independent  high as possible,  then the difference approximation  n  as s  i s said to be consistent  the order of c o n s i s t e n c y i s equal to n s.  with the differential equation and To  and  specify y(x) completely it i s n e c e s s a r y  conditions, which a r e a s s u m e d to take the  to i m p o s e boundary  form  n (0) k  (1.1a)  B (0)y(0)=  £  k  b ^ .(0)y (0) = b ( 0 ) ,  l<k<n ,  b .(l)y  l<k<n-n ,  (l)  k  0  i=0 and n (l) k  (1.1b)  B (l)y(l)= k  ( l )  k >  (l) = b (l), k  0  i=0 where  n (0)  <  k  It w i l l always be  n  and  n  _r.(l)  •  a s s u m e d that the d i f f e r e n t i a l equation (1.1) subject  to the boundary conditions ( l . l a , b ) admits many times  n  a unique  solution that i s as  continuously differentiable as needed.  Similarly, approximations  to define a finite difference scheme that determines  u. to y. completely one has  approximations  conditions (1.1a, b) that look like  s  (1.2a)  B  h j k  (0)u EE 0  k< > 0  YJ  i=0 and  d k  , i  (  0  )  u  i  =  b  k  ( 0 )  '  i  <  k  <  n 0  '  the  to the boundary  -30 (1.2b)  B  h  (  k  ( l ) u  j  S  - 2  d  k >  .(l)u  J + i  =b (l).  l < k < n - n  k  0  .  i=-r (l) k  H e r e b (0) k  The  and b ( l ) k  a r e approximations to b (0) k  and b ( l )  respectively.  k  definitions of truncation e r r o r and c o n s i s t e n c y are  similar  the previous definitions and the equations are compact if s (0) = ( 0 ) n  k  r (l)  = n (l).  k  k  order (1.2)  A finite difference  s_ if the approximate and (1.2a,b) and  k  to and  scheme is said to be (accurate) of  solution Uy  |y^-u. | < K ^ h  S  for  (0 < j  < J),  is uniquely defined by  some constant  that does not  depend on h. N u m e r o u s texts and papers a r e partly or entirely c o n c e r n e d with finite difference approximations to (1.1) books by Collatz (I960) and K e l l e r  and ( l . l a , b ) .  In p a r t i c u l a r ,  the  (1968) need to be mentioned h e r e . The  general theory of difference approximations to systems of boundary value p r o b l e m s is contained i n papers by G r i g o r i e f f (1970) and K r e i s s (1972). An  extensive survey of finite difference methods for boundary value  p r o b l e m s can be found i n a paper by K e l l e r  (1975).  V a r i o u s finite difference approximations of the f o r m (1.2) the l i t e r a t u r e .  (1.3)  To the second o r d e r differential  y"(x)  + a^xjy'fx) + a°(x)y(x)  = f(x)  equation  ,  the s i m p l e s t difference approximation is perhaps given by  (1.4)  appear  in  -4H e r e D. u. = (u h J J+1 uniform, i . e . , h. = h.  -u. ,) / 2h and the m e s h i s The  order of consistency i s equal to two,  as i s  the order of a c c u r a c y of a difference scheme based upon (1.4), p r o v i d e d that the approximations  to the boundary conditions are consistent.  F r e q u e n t l y i t i s d e s i r a b l e to use equations for which the o r d e r of consistency i s m o r e than two.  One  way  of obtaining such equations i s by  allowing the width of the approximation to be greater than n+1. example,  (1.5)  The  a fourth o r d e r finite difference approximation to (1.3) i s  L u. h  For  EE  width of this approximation i s equal to five.  If a difference a p p r o x i -  mation i s not compact then an entire difference scheme cannot be upon that equation alone. used.  here.  N e a r the boundaries other equations need to be  In (1.5) for example,  p e r i o d i c i t y i s assumed.  based  L^u. i s not defined i f j = 1 or j = J-1 unless  Hence other approximations need to be  E x a m p l e s of these are given by Shoosmith (1975).  used  In this chapter  such extra boundary conditions w i l l be d i s c u s s e d i n Sections 7 and  8.  It is not n e c e s s a r y to i n c r e a s e the width of the approximation to obtain higher o r d e r equations, for one can also r e q u i r e m o r e " l o c a l i n f o r m a t i o n " to be available.  In this case each difference equation  involves the values of the coefficient functions and the inhomogeneous t e r m at m o r e than one point.  F o r example,  to  (1.6)  y"(x) + a°(x)y.(x) = f(x)  a fourth order approximation  -5i s given by  (1.7)  D*u.  + K ( f ) (a°u.) = K ( f ) f . , h  h  with  K (c)w. h  The  c o n s t r u c t i o n and  s  i|£. w.^  + cw.  + f w. 1  £  + 1  a n a l y s i s of such higher o r d e r  mations i s also included i n this  .  compact a p p r o x i -  chapter.  U s u a l l y the statement that a c e r t a i n difference approximation given order Systematic  of consistency can be justified by  unless one  methods such as G a l e r k i n , Ritz and  e.g.,  V a r g a (1966),  Shampine (1972), Strang and  F i x (1973) and  Swartz (1973),  Schultz (1973)).  compare v a r i o u s compact equations,  with  piecewise  ( F o r a d i s c u s s i o n of these,  C i a r l e t Schultz and V a r g a (1967),  deBoor and  appear i n f r e -  c o n s i d e r s finite element  collocation procedures  as finite difference methods.  but no  a  T a y l o r expanding.  methods of d e r i v i n g finite difference equations  quently i n the l i t e r a t u r e however,  polynomials  simply  has  see  R u s s e l l and  as w e l l as the books of B i r k h o f f and  Gulati: (1974)  attempt i s made to m e t h o d i c a l l y  construct such compact difference approximations. a p r o c e d u r e for d e r i v i n g difference approximations  Swartz (1974) d i s c u s s e s that r e s e m b l e s  the  p r o c e d u r e given i n this thesis, but no complete a n a l y s i s i s given.  The  m a i n purpose of the present  c l a s s of finite difference approximations their o r d e r  of consistency.  chapter  i s then to d e r i v e a l a r g e  to the p r o b l e m (1.1) and  T h i s i s done i n the next two  to give  sections.  Section 4 a number of examples i s given, including known as w e l l as  In new  -6difference equations. approximations  I n 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  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 .  t h e s e a p p e a r i n S e c t i o n 6.  g i v e n i n this chapter.  THE To  C O N S T R U C T I O N OF  u p o n the  i n S e c t i o n 8 the  are given.  FINITE D I F F E R E N C E  APPROXIMATIONS  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 o f the f o r m  the d i f f e r e n t i a l e q u a t i o n (1.1) one i n d e x j l e t z , < z_ < • • • < z 1 2 m J  satisfying  Finally,  | z . . - z , | >—-  h,  be  may r  i # k\  proceed  as f o l l o w s .  For  (1.2) to given  p o i n t s i n the i n t e r v a l Tx. ,x.. 1, j - r j+s L  (To  J  s i m p l i f y the n o t a t i o n the  subscript  K  1  j d e n o t i n g d e p e n d e n c e o n 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 omitted here,  so that,  e.g.,  exceeding  P_j d e n o t e the s p a c e d=r  + s + m - l .  be  z. . b e c o m e s z.). J> 1  Let  of  (1972) i s  of d i f f e r e n c e s c h e m e s b a s e d  r e s u l t s of s o m e n u m e r i c a l e x p e r i m e n t s  2.  Examples  I n S e c t i o n 7 the t h e o r y of K r e i s s  a p p l i e d to i n v e s t i g a t e t h e s t a b i l i t y difference approximations  difference  1  of a l l p o l y n o m i a l s w i t h d e g r e e  Assume  t h a t f o r s o m e i n t e g e r i,  not (-r < l< s ) ,  SL  a p o l y n o m i a l p (x) e P ^ (2. 1) P (x )_ =  is determined u j+i  uniquely by  i +  r < i < s , (2.2)  Lp (z.)  e q u a t i o n s (2.2) w i l l b e  the e q u a t i o n (2.3)  equations  t  f(z.) 1 < i < m  (The  the r + s + m  referred  to as  . "collocation"  equations.)  Then  relates  the q u a n t i t i e s j ^ »  (-r<i<s),  u  by an e x p r e s s i o n  (difference  e q u a t i o n ) of the f o r m  m (2.4)  u.  =  V  d.u.,. =  V i=l  l =-r  In o r d e r  to b e a b l e  and  e^ i n (2.4),  so  that one c a n  e.f(z.),  to e x h i b i t e x p l i c i t  expressions  l e t the p o l y n o m i a l s w ( x ) , 1  f o r the c o e f f i c i e n t s d^  (0 < i < d), f o r m  a basis  of  P^,  write  P*(x)  =  Yi  c w (x), 1  i  i=0 Also,  define  a matrix  r  0, (x. , ,) j-r+1  w (2.5)  D. Lw  x  (x. , ) j+s  Lw"(z  If the r e l a t i o n s (2.1),  ing  w  w  (x. , ,) j-r+1'  w  r  v  (x. , )  d. . (x. ) J-r'  u  v  (x. , ,) j-r+1'  u. , , J-r+1  x  w  (x. , ) j + s'  U  m  j  + B  f(z ) A  Lw^z^)  )  Lw*(z  m  f(z ) 2  )  Lw (z u  (2.2) a n d (2.3) a r e s a t i s f i e d  the d e t e r m i n a n t b y the l a s t  factor E  (2.6)  w  (x. j _ ')  (z^  Lw^(z^)  Evaluating  w  (Xj_ )  w  w  by  m'  )  f(z ) m v  t h e n d e t ( D ^ ) = 0.  c o l u m n and introducing  a  normaliz-  o n e f i n d s that the c o e f f i c i e n t s d^ a n d e^ i n (2.4) a r e g i v e n  d  i  =  c  o  f  L [ V i  ]  /  E  '  by  -8-  and  (2.7)  e. = - c o f [ f ( z . ) ] / E L  where cof-^f • ] i s the cofactor of the given element i n the m a t r i x D ^ , w i t h 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 m E i =l  = "  E  where  c  o  f  L  as  '  L y ( x ) EE y ^ x ) . Q  If the b a s i s functions w (x) are chosen such that 1  (2.8)  ( j.  w  x  r + k  )  =  6  f i  k  o  r  0 < i < r + s ,  and  (2.9)  w ( j x  + k  ) =  0  f  o  r r + s + l < i < d ,  r  -r < k < s ,  then j r+i. * , r+s + l \ Lw . ( z ^ . , lLw.u , (z^,;, , T  (2. 10)  d  m i "  , -r < i < s ,  E j r+i. Lw (z  and  Lw^(z ,  m  > )  r+s + l , . Lw (z ) m  T  x  Lw (z d  m  )  -9T  Lw  Lzll  e,  Lw (z d  T r+s+1, . Lw (z. ,) l-l j r+s+1, . (z )  m+i+1  (2.11)  r+s+1, . ( z ^  E  Lw  ,  r  W i t h the above c h o i c e normalizing  factor  of b a s i s  + 1  (-1)  related  )  e x p r e s s the  T  Q  Z l  m ,  = 1 then E  m  one c a n a l s o  t  r+s + 1, .. (z ) m/  1 L„w 0  If m  Lw^(z  r+s + 1, > . . . r+s+m-1, . (z ) L w ( )  j  Q  =  '  as  ,  E  "  functions,  1 L w  (2.12)  m  m  d  )  S  1 < i <  Lw (z. )  i + 1  Lw ^ "'"^(z  (z^ ^)  = -1  T  T  a n d e, = 1.  Further,  1  r+s+m-1, . L w (z 0 rrr n  v  v  the c o e f f i c i e n t s  d. a n d e. a r e i  i  by  m d  Remark. show  to  L k=l  =  L , w r + 1  T h e above f o r m u l a s  the f o r m  products  i  numerical problems,  Thereafter, s u c h a way,  z  e  collecting  functions one m a y  that t h e s e  "  r  -  -  1  In p a r t i c u l a r ,  that w i l l first  the v a r i o u s  the n u m e r i c a l  '  S  '  f o r d^ a n d e.. a r e e s s e n t i a l l y  of t h e s e c o e f f i c i e n t s .  of c o e f f i c i e n t  symbolically,  ^ k^ k  evaluate  the type o f  F o r actual  application  the d e t e r m i n a n t s  multiplicative  evaluation  coefficients  appear.  note  i n t e n d e d to  and additive  of the d^ a n d e^ m a y  constants.  be o r g a n i z e d i n  a r e o b t a i n e d i n the s m a l l e s t  possible  10number of a r i t h m e t i c  operations.  B a s i s functions w (x) s a t i s f y i n g (2.8) and (2.9) a r e 1  constructed.  In fact,  easily  i t i s convenient to a s s u m e that the w (x) a r e given 1  by r +s  rr  (x-x. ,. ) j-r+k  k=0 w (x) = k ^ i r+s (x. ,.-x. ,.) ' J-r+i J-r+k'  (2. 13)  0 < i < r + s ,  1  k = 0  k#i  and  m-1  r+s  TT(*-t (2. 14)  w  r + S + 1  (x) =  TT (x-x.j - r + k '  t  k=0 r+s  k=l k^a m- 1  TT<t_-t ) k  k=l M i  H e r e the t., I  (1 < i < m - 1 ) , — —  Lo  1 < i < m - 1. 1  J  "  r + k  1  are d i s t i n c t points i n the i n t e r v a l Tx. , x . , 1. j - r j+s L  r  J  It i s not difficult to check that the points t\ c a n be chosen i n such a way that the p o l y n o m i a l s w (x), 1  The p r o c e d u r e  (1 < i < d),  f o r m a l i n e a r l y independent  set.  leading up to the definition of the coefficients d^ and  e^ i s equivalent to a method of undetermined this c o n s i d e r the following e x a m p l e .  coefficients.  Let n = 2 ,  To i l l u s t r a t e  r = s = 1 and m = 2; i . e.,  we a r e 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 a difference a p p r o x i m a t i o n of the  form  -11-  (2.15)  d  _! j_l u  to the differential  +  d  u 0  + l j+i d  u  =  e  i ( i) + 2 ( 2 £  z  e  f  z  '  )  equation  Ly(x) = y  (x) + a (x)y (x) + a (x)y(x) = f(x)  T h i s can be accomplished polynomials  j  p(x) e P^.  .  by r e q u i r i n g that (2.15) be satisfied for a l l (Replacing f(z^) by Lp(z.)).  as i n (2. 13) and (2. 14), this requirement  With w (x), 1  ( 0 < i < 3),  leads to the following l i n e a r  system of equations:  1  0  0  •Lw^(z  0  1  0  -Lw (z^)  0  0  1  -Lw  0  0  0  -Lw  0"  L  •Lw^(z.  0  (z p  -Lw  CO  (z ,)  -Lw^(z )  2 3  (z ) 2  0'  2  A s s u m e that the rank of the s y s t e m i s equal to four, so that one c a n r e w r i t e the equations as  1  0  0  .  0  1  0  .  -Lw'(z ^)  0  0  1  .  2 -Lw (z^)  0  0  o  •  3 -Lw (z^)  where the m a t r i x  - l  Lw^(z^) e^  d  o  Lw (z )e  2  d  l  Lw (z )e  2  e  l _  d  1  2  2  2  3 _  on the left i s nonsingular,  Lw  (z )e 2  i . e . , Lw  2  (z.)±0.  Partition  12this  matrix  into  I  A'  0  B  as  indicated. -1  by  elementary  cofactors  of  linear  B  T  .  algebra  In  this  1 = ^g^g ^ C ,  B  simple  I 0  Then  example  B  A B  -.  -1  where  -1  C  Q  =  /  Lw (z )]e  2  ,  =  [Lw (z )  -  Lw (z )Lw (z )  /  Lw (z )]e  2  ,  [Lw (z )  -  Lw (z  1  2  and  3  [ -1/Lw  •Lw 2  =  3  1  3  1  2  (z )]Lw  3  4  Then  E  the  (  =  [Lw^(zj)Lw^(z_)  2  3  1  More the  matrix  2  z 2  )  1  3  1  Lw (z  e  3  2 >  -  e 2  coefficients  [Lw (z )Lw (z )  in  )Lw (z ) /  3  (z^)  dp  is  2  3  1  =  Lw  =  -Lw  agreement  generally  (z_)  /  (z ^) /  with  the  2  are  equal  to  Lw^(z )Lw (zp] / E 3  [ Lw"(Zj)Lw^(z_)  which  2  of  Henc e  Lw°(z )Lw (z ) 1  and  (z ^)  -  2  ,  matrix  [Lw°(z )  2  e  the  =  2  Let  is  -1 Lw  d  I - A B 0 B"  2  -  Lw^(z )Lw (z^)]  /  Lw (z )Lw (z^)]  /  3  2  2  3  2  E  ,  E  E , E  the  above  equations  (2.10)  procedure  is  and  (2.11).  well-defined  provided  that  -13-  T r + s + 1. > Lw (z^)  • • • Lw  j r + s+1, . Lw (z ) m  . . .  v  has  rank  zero  e q u a l to m-1.  h a s r a n k m-1  W i t h the d e f i n i t i o n s g i v e n  Theorem be and  given  of t h i s  2. 1  this  for a l l small  factor E  i s non-  e n o u g h h.  far, i ti s possible  to s t a t e the b a s i c  L e t the c o e f f i c i e n t s d. a n d e. a n d the n o r m a l i z i n g  that r + s > n.  (2. 11) a n d (2. 12) r e s p e c t i v e l y . If h i s s m a l l  d. a r e n o n z e r o  Assume  factor  that E ± 0  e n o u g h t h e n a t l e a s t n + 1 of the  a n d the o r d e r  l  difference approximation  Proof.  r+s+m-1, > Lw (z ) m  chapter.  b y (2. 10),  coefficients  T  r+s+m-1/ (z  I n p a r t i c u l a r i f the n o r m a l i z i n g  t h e n the m a t r i x  theorem  T  of consistency  o f the f i n i t e  1  (2.4) i s a t l e a s t e q u a l to r + s +  m-n.  Define  (2. 16)  d. = c o f  (2. 17)  e. =  i  T  L  -cof  L  [u.,.l / E, j+i  J  '  [f(z.)] / E ,  T J  -r < i< s, -  -  1 < i < m  ,  0  m so  that  Y  e  i  - 1-  i =l First coefficients  itwill  be s h o w n  d. a r e n o n z e r o . 1  l < k < k ^ < n  +  _  that i f E # 0 t h e n a t l e a s t n + 1 o f the F o r s u p p o s e o n the c o n t r a r y  1 a n d d^ = 0 o t h e r w i s e .  that d . + 0, k  Let ,  E  141 p(x)  EE  has  degree  L  certainly  k=i n and  n < d.  q(x), w h e r e q(x) e P , i s c h o s e n s u c h that p(x) k-l 1 leading coefficient — . B y a s s u m p t i o n n < r + s, s o Hence by  construction  m I i = -r But  from  the d e f i n i t i o n  Therefore  _  ^  e^ = 0,  provided  The  -E i oP< i> i=l e  j + i  L  z  of p(x) i t f o l l o w s  °-  =  that  ^  which  truncation  LQP(X)=1,  +  r  is a contradiction.  So  i f E + 0 then n + 1 or  _ d. / d. = 1 + 0(h), the s a m e l ' l  Since  h is small  Pj j_ = 0 a n d  d^  i=-  _ i=l d. a r e n o n z e r o . l  more d^,  m  ip  d  is true  f o r the  enough.  error  T. i s d e f i n e d  as  J m r  where  j -  V j  s  ~ I i ( i> i=l e  conditions  = I Vj+i i = -r  z  y(x) i s the s o l u t i o n  boundary  - I i i=l e  £  (  z  i  •  )  to the d i f f e r e n t i a l e q u a t i o n  (1. l a , b).  d T.  f  m  y  Taylor  <w  expand  to get  r- d  m  (1.1) s u b j e c t  y  =  (k) k!  k=0  i=-r  +  V i = -r  d  k=0  d.(x. ,.-x.  k=0  £=1  ,d+l y  (  d  +  1  )  (-)  (d+l)! -  ™ -  I i=l  ^  (z.-x.)  .  1  J  d+l  y  ( d + 1 )  (t  i :  (d+l)!  to the  -15-  ^ V  A  +  > t-i i=-r  for  some  the  quantity  i  xd+1  y  d.(x.,.-x.) l j+i j  between  square brackets  follows such  w i t h the u n d e r l y i n g that t h e r e  equals  zero  By  f o r e a c h k.  d. a n d e. a n d the b a s i s 1  functions  construction F r o m the w (x), 1  1  v  assumption  a r e constants  ,  v  s. b e t w e e n x. a n d x., . a n d t- b e t w e e n x. a n d z. . i J J+i i J i  d e f i n i t i o n s o f the c o e f f i c i e n t s together  1  v _ i _ r / ^+1 - > e. — , , , ... L(z.-x.) Li I (d+1)! I y i=l  i  [S }  — . ,, ,, , — (d+1)!  C^,  that ^  < tu < h,  and  (1 <_ j <_ J), i t  that do not d e p e n d  o n h,  that  d^[  e  < C^h  i I —  ,  n  ^2  -r < i < s ,  '  1  < i  <  m  ,  and  T / ^d+1I ^ „ , d + l - n L(z-x.) < C.h , J — 3  z. < z < z 1 — — m  1  W i t h these  estimates  the s e c o n d  conclusion  o f the t h e o r e m  follows  i m m e dia tely.  By  means  of two e x a m p l e s  the n e c e s s i t y  of the c o n d i t i o n  that E  be  n o n z e r o w i l l be i l l u s t r a t e d .  Example  2. 2  Consider  the c a s e w h e r e n = 2, m  T h u s we  a r e 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 II  a p p r o x i m a t i o n to the p r o b l e m the y-' J  p o i n t x., that i n v o l v e s J y-iJ-1 J' i '  y-x?J+^  This  Ly = y  of a finite d i f f e r e n c e  1  F  '  0  (x) + a (x)y (x) + a (x)y(x)  the a p p r o x i m a t i o n s F  = 2, r = 1, a n d s = 2.  = f(x) a t  u. <, u., u.,., u . , J-l J J+1 J+2  i s a c c o m p l i s h e d b y the p r o c e d u r e  o f this  0 y  to y . J - l '  section,  -16with  two  collocation points.  m e s h is uniform. (2. 14).  The  In t h i s c a s e  For  basis  they  simplicity  i t i s also assumed  f u n c t i o n s w (x)  a r e defined  X  + 1  = ( x - x _ ) ( x - x . ) ( x - x . ) / 2h  w (x)  = -(x-x._ )(x-x.)(x-x  j  2  4  ,  /  assumed = 3,  to be  i  e  h  a  i s i n this  so £  2  ± S_^.  z^ a n d  x  points.  z  x  One factor  c a n not a l l e v i a t e E.  For  d  The  2  distinct,  n  z  2  =  X  /  / !  6h  3  )(x-x.  j-1  example  3) / h .  2  j  ) /  ,  + 2  ,  4 )/h .  by  (x^ + j_|__) / 2, , w h i c h w o u l d  normalizing the  E  +  ) / 2h  V  + 1  (<_, -<r;^)(£ + 2  W  (x) = ( x - x ^ . _ ) ( x - x ) ( x - x .  i . e . , w h e n the p o i n t s  midpoint these  =  j+ 2  1  defined  ,  + 2  +  X  LQW^(Z^)  W  1  factor  1  1  1 = j-1  normalizing  +  (x-Xj_ )(x-x^)(x-Xj  the c o l l o c a t i o n p o i n t s be  The  1  1  (x) =  w  Let  ,  + 2  w^x)  3,  as i n (2. 13) a n d  c a n be w r i t t e n a s  w ° ( x ) = -(x-x.)(x-x. ) ( x - x . ) / 6h  w  that the  g i v e n by E  E  = L^w  seem  the m o s t  E  2  z^ a r e a l w a y s  symmetrically  the p r o b l e m b y with  (z ) -  - 0 i f and only i f ^  are placed  example,  h  c o l l o c a t i o n points  Hence 2  ^2  +  =  n a t u r a l p l a c e m e n t of r e d e f i n i n g the  l/h , ^  Q  ( ) Z l  L w Q  = 1 and  ( ) Z l  -1 < i <  d. = h' l  T  L w Q  \  (z ) 2  4. . L w (z ) T  Q  2  2  a b o u t the  formula  L w  + i_,  2 ,  = 2  -17gives  d_  —  2  e^ y  =  1  2  -2/h  , d  4 -2  and  (-2u.j . l  +  2h  ("Vl  -  one  +  -  3 U  j  take m basis  As  = 3.  So  functions  1  + 2u  +  U  in  y  w (x),  a basis  J  o f P^.  h  =  3  - 2f(  x  ( h  )  +l  f  i (x) = f (x),  4  2  .  Further,  so that f o r the  )  j+ 1  problem  - 2f(x.).  "  E  f  (  than y  /  '  h  J +  (x) = f ( x ) .  l e t n = 2,  i n example  Then  )  }  II  rather  example  J  j  X  collocation points.  }  -  .). J +  Now  L w (z ) E  ( Z ^ ) = 2,  r =  1 and  s = 2,  but  I n a d d i t i o n to the  (2.2) take  (x) = ( x - x . •. ) ( x - x . ) ( x - x . A ) ( x - x .  ^  Q  W  ) /h*  2 ) /  +  (0<i<4),  1  x. _ i = \ ( x . + x . J+2  j  a r e three  J  where  j + 2  2  = 2/h  2  becomes  i n the p r e v i o u s  there  w  l  +  which i s consistent with  2. 3  equation  +  d  4  -h L Q  =  and  obtains  j  3 U  6u.  2  = -6/h 2  e^  (x) = f ( x ) t h e d i f f e r e n c e  Example  ,  —  = h L Q W (Z^) =  Dividing by  2  - 6/h  Q  J  2  . ) ( x - x . _) / h  i  +1  the p o l y n o m i a l s  ,  +C  w (x), 1  (0<i<5),  form  4  5  1  i s given  by  4  L w (z ) 5  Q  2  L Q W  *  -  5 •  L Q W  '(ZJ)  L Q W  (Z^)  L ,  (Z )  L w (z )  Q  W  ( Z ^  ( Z J )  = L w Q  If the m e s h distributed  4  (z ) 3  L w Q  i s taken  5 (z )  to be  symmetrically,  then after s o m e  L w  3  4  Q  uniform i . e . , z^  computation  (z ) 3  L Q W  5  L w (z )  4  Q  3  5  2  Q  and i f the c o l l o c a t i o n points =  x  j _ ^ i . - £h,  i t i s found  that E  z  = j x  2  = c£  3  (4£  X  i '  ?  z  3  =  x  2  are j+- *  5 -3)/h , where  c is  -18some constant. zero.  Since the points z^ a r e d i s t i n c t ,  § cannot be equal to  Hence i f the points z.. a r e p l a c e d 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 . 8 6 6 . Remark.  A n e c e s s a r y c o n d i t i o n for the n o r m a l i z i n g factor E to be n o n -  z e r o i s that r + s > n .  F o r example,  the width of a difference a p p r o x i -  m a t i o n to a second o r d e r d i f f e r e n t i a l equation i s at l e a s t t h r e e . should not be s u r p r i s i n g . d w (x) / dx n  1  11  =0,  To p r o v e the a s s e r t i o n suppose r + s < n . T h e n  (0 < i <.r + s), because the f i r s t r + s + l  have degree r + s.  This  basis  functions  Hence the f i r s t c o l u m n i n the determinant defining d..  c o n s i s t s of z e r o s for each i .  Thus d.. = 0,  (-r < i < s), and f r o m  the  proof of t h e o r e m (2. 1) i t follows that E = 0. Some sufficient conditions for E to be nonzero a r e given i n the following  theorem.  T h e o r e m 2.4 or i f r + s = n.  The n o r m a l i z i n g factor E i s n o n z e r o i f m = 1 and r + s > n, ( i . e . , when the difference equation i s c o m p a c t . )  If m = 1 then E = - 1 .  Proof.  If- r + s = n, m > l  and E = 0,  then i t  follows f r o m (2.12) that there a r e constants c , (r + s < i < d), not a l l d z e r o , such that the p o l y n o m i a l q(x) = ( __, c^w (x)^ e P ^ s a t i s f i e s 1  1  i=r+s+l q^ \z/) n  -  c r + g  >  (1 < i < m ) .  satisfies q ^ ( z ^ ) =0.  Thus the p o l y n o m i a l q(x) = q(x) - c _  x /n! n  + g  B u t this i s i m p o s s i b l e since q^ ^(x) € P _ _ ^ because n  n  r + s =n. Remark.  If the a p p r o x i m a t i o n i s compact then for each n the coefficients  dj, a r e uniquely d e t e r m i n e d , points.  independent of the number of c o l l o c a t i o n  F o r e x a m p l e , i f n = 2 and r + s = 2 then d ^  VVVi  1  -19-  — d.g =  -2  — and d^ = ^  2 (h J . h f ' j j+1 j+l j j+1 second d e r i v a t i v e i s always g i v e n by u., . - u.  __±i  h..  1  s  o  t  n  approximation  e  u. - u... ,  _ J Lli h. J  J+1  3.  ^h.at  E F F E C T OF THE CHOICE OF COLLOCATION It i s w e l l known, (see,  collocation procedures  e.g.,  POINTS  R u s s e l l and V a r a n (1975)),  w i t h c e r t a i n p i e c e w i s e p o l y n o m i a l s for  that  the  n u m e r i c a l solution of boundary value p r o b l e m s have a higher o r d e r of a c c u r a c y i f G a u s s i a n points a r e used as c o l l o c a t i o n p o i n t s . expect,  One w i l l  that a s i m i l a r statement applies to the finite difference  d i s c u s s e d i n the p r e v i o u s s e c t i o n .  equations  H o w e v e r , i t m a y be s u r p r i s i n g that  the points z^, for w h i c h an e x t r a h i g h o r d e r of c o n s i s t e n c y i s do not g e n e r a l l y c o i n c i d e w i t h the G a u s s i a n p o i n t s .  attained,  Nevertheless,  these  s p e c i a l points c a n be c h a r a c t e r i z e d r e l a t i v e l y e a s i l y as w i l l now be shown. C o n s i d e r again the b a s i s functions w (x) of P X  6  equations  (2.13) and (2.14).  extended to f o r m a b a s i s of P  v  3  T h i s l i n e a r l y independent , , by adding a r+s+m  •r^r+s+m that vanishes at the m e s h p o i n t s . form  . , defined by r+s+m-1 set m a y be  polynomial w ^ }  v  (x) e '  Thus this p o l y n o m i a l has  the  -20-  m-1 w  r+s  TT^-t,).  r+s+m, . (x) =  rr t with  t  m  FT  k=l m-1  j.  k=0 r+s  w  r + k  )  TT ( t - ^ . )  v  m  rtk  * t , (.1 < k < m-1), a n d t^_ + x., ( - r < j < s ) . E x p a n s i o n  o f y(x) i n  k  t e r m s o f the p o l y n o m i a l s w (x) g i v e s r+s+m = I  y(x)  c,w (x) + 0 ( h k  r  +  s  +  m  +  1  )  .  k=0 S u b s t i t u t i o n of t h i s T  j  E  s  e x p r e s s i o n i n t o the t r u n c a t i o n e r r o r m  Vj+i  •  E  e  i  f  (  z  i  y  }  i e l d s  i=l  i = -r  r+s+m  m E  k=0  The  w  r  +  s  +  m  j  ( j+i) - E x  =  m V > u  i  e  L  +  < i>  w  z  0 ( h  r  +  s  +  m  "  n  +  1  )  i=i  square  ( X j _ ^ )  = -c , , r+s+m  T.  i  w  i =-r  quantity between  Moreover,  d  0,  brackets  equals  (-r < i < s ) .  zero  f o r 0 < k < _ r + s + m - 1,  Hence  T r+s+m, . ~ ,,r+s+m-n+l> e.Lw (z.) + 0 ( h ) . i i'  i=l Estimates 0(h ^" ^" r  S  m  a s i n the p r o o f n  ) .  However,  of t h e o r e m  (2. 1) a g a i n  it is also clear  from  r e v e a l that T\ =  the a b o v e  an additional order  expression,  that  of c o n s i s t e n c y i s o b t a i n e d , i f the c o l l o c a t i o n p o i n t s z. , , ,n r+s+m, » •j -i-i. • i. v. T r+s+m, , d w (x) _ , c o i n c i d e with points where L Q W (x) = — —-. = 0. T h u s the dx f o l l o w i n g h a s b e e n shown: s  Theorem E  3. 1.  L e t the c o e f f i c i e n t s d. a n d e. a n d t h e n o r m a l i z i n g f a c t o r  b e g i v e n b y (2.10),  and  (2.11) a n d (2. 12) r e s p e c t i v e l y .  ~ r+s+m Let w (x) =  t h a t r + s > n.  Assume  m-1 r+s ~| |~ ( x - t , ) 1 T ( x - x .  k=i  k=o  K  J  "  that E + 0  , ) and a s s u m e r + K  s*** -f s T in. t h a t LpW  ( z ^ = 0,(1 < i < m ) . T h e n the o r d e r  difference approximation If m  (2.4) i s a t l e a s t e q u a l t o r + s.+ m  - n + 1.  - 1 a n d 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 c o m p a c t , i . e . , r + s = n,  ~ r+s+m n then w (x) = ~| ["(x-x. k=0 J  , ). _  r  I n t h i s c a s e t h e r e i s only one  +  collocation point f o r which a higher m  of c o n s i s t e n c y of the f i n i t e  order  of c o n s i s t e n c y i s o b t a i n e d .  = 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  point.  ( A s s u m i n g E ± 0.)  -parameter  F o r the c a s e m  f a m i l y o f p o i n t s z^, (1 < i < m ) ,  of c o n s i s t e n c y i s o b t a i n e d . i n t h e d e f i n i t i o n of w these 4.  of t h i s  > 1, r + s = n t h e r e i s a f o r w h i c h the i m p r o v e d  (m-1) order  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 ^ (x) i n t h e o r e m (3. 1).  Particular  e x a m p l e s of  s p e c i a l c o l l o c a t i o n p o i n t s w i l l be i n c l u d e d i n the n e x t s e c t i o n .  EXAMPLES  E x a m p l e 4. 1. two  If  OF  FINITE DIFFERENCE  L e t n = 1, r = 0 a n d s = 1;  point finite difference approximation  (4.1)  Take m  APPROXIMATIONS  i . e . , we  w a n t to c o n s t r u c t a  to the e q u a t i o n  L y ( x ) = y'(x) + a°(x)y(x) = f(x)  = 1 and z  i  = x^ ± +  = _-(x.. + X j  +  1  ).  T h e n the b a s i s f u n c t i o n s w ( x ) a s 1  d e f i n e d b y (2.13) a n d (2.14) a r e w°(x) = - ( x - x . .) / h., a n d w ( x ) = ( x - x . ) / h . , J+1 J J' J T h e d i f f e r e n c e a p p r o x i m a t i o n h a s the f o r m X  •22where  the c o e f f i c i e n t s  d ,  a n d e^ c a n b e  Q  obtained  from  (2. 10) a n d  (2. 11), v i z .  A  0  T  d  Q  = Lw  1  = Lw  ,  "I  V  ( )  =  Z l  ,  1  +  2  0 a. ± , +  and  1/ \ 1 , 1 0 ( ) = h. — + 2 a, , i 1+1  T  A  d  ZjL  J  The  difference  approximation  (4.2)  {u  c a n then be w r i t t e n as  j r j j j+i V j+i u  )/h  +ia  (  u  •  ) =  +  which is recognized Keller least  (1969). equal  as the c e n t e r e d  According  to one.  That  follows  from  w (x) =  (x-Xj)(x-x.. ).  2  to t h e o r e m the o r d e r  (3.1),  because  equation,  or Box  (2.1) the o r d e r  scheme  of  of c o n s i s t e n c y i s a t  of c o n s i s t e n c y a c t u a l l y z^ i s the c r i t i c a l  equals  two  p o i n t of  + 1  An  equally simple  d e r i v e d when the  theorem  Euler  form  taking m  difference  equation,  = 2, w i t h z, = x. a n d z_ = x. , 1 j 2 j+1  d_.u. + d . u. , . = e . f. + e f . ... 0 j 1 J+1 1 J 2 j+1  1 E  2  This  equation  ' , :  0  With w  the t r a p e z o i d a l m e t h o d , i s  2 (x) = ( x - x _ . ) ( x - X j ^ ) / h. the c o e f f i c i e n t s a r e  Lw  0  (z^)  Lw  2 (z^) d  Lw°(z ) 2  Lw (z ) 2  2  l  =  E  Lw  1 (z^)  Lw  2 (z^)  Lw  1 (z )  Lw  2 (z )  2  2  has  -23-  e  _  1 2 = __" L  w  L^y(x) = y  (_.) z  i  a  (x).  n  d  e  2  -1 ~ET "  =  L i W  The equation  2 ^ i)' z  w  i  t  n  = Q  E  L  ^ l^  w  2  z  ~ 0 L  2 ^ 2^  W  Z  A N C I  thus obtained i s  (u., .-u.) / h. + |(a°u.+a°, ,u. .) = _-(f'.+f. ,). J+ 1 J J J J J+1 J+1 J j+l' 2 V  B y theorem (2. 1) the o r d e r  of consistency of this difference equation i s  equal to two.  Next, consider the general two point difference approximation m  = 2.  That such a difference approximation  f r o m theorem (2.4). i s equal to two.  i s always consistent follows  T h e o r e m ( 2 . 1 ) i m p l i e s that the o r d e r  From  with  of consistency  theorem (3. 1) i t follows that the o r d e r  becomes ~ 3  three i f the collocation points coincide with the c r i t i c a l points of w (x-t, )(x-x.)(x-x. , ,). 1 J J+1  (x) =  L e t t, = x. + r i h . , z. = x. + £,h. and z_=x. + £_h.. 1 J J 1 J 1 J 2 j 2 j  T h e n the special values  of ^  i  and £  2  that give a third o r d e r  the roots of p (£) = 0, where p(£) = £(£-l)(£-n).  formula are  These roots a r e  rz—  f  t  _  5  -  V+ 1 ± v ••  ••  f) -T7 + 1  3  If n = •§• then the roots a r e £ = \ ± ^ <y~3 . points a r e p l a c e d s y m m e t r i c a l l y  In that case the collocation  i n the i n t e r v a l [ x ^ x ^ ^ ] .  i n another o r d e r of consistency being gained, consistency i s then equal to four.  This  so that the actual o r d e r of  The coefficients of this  approximation  are -j d  Q  =  - 1 —  , 1 , 0 , . , 0, + (a ( z ^ + a ( z ) )  1 , 0, . 0, . - y-r h.a ( z ^ a ( z ) ,  2  2  j  d«_ = j j - + 2f (  J  a  ( _) + z  a  ( 2^ z  +  results  T7 i h  a  < l z  ) a  ( 2*' z  -24-  e  = 2 "  1  j  h  X j ^ ] .  [XJ,  The  However,  a n d  procedure 2 m.  ~ 2 +  e 2  of S e c t i o n 2 w i t h m  (1974),  E x a m p l e 4.2. *  With  (See B u t c h e r  r = 0,  s = 2, m  0 j U  +  d  l j +l U  +  d^ = L w ( z ^ ) ,  (1964), with  difference approximation  m  d  2 j+2 U  =  f ( x  j+2  }  E x a m p l e 4.3.  points w i l l  yield  t h e o r e m (3. 1).  equations.  The  Wright  (1970).)  p. 217.)  Higher  z . = x. , . 1 j+ s  L e t n=2  and  h. ,,. = h., _ = h j+1 j+2  The  form  .  to  be  = f , j+ 2  specially suited for stiff equations.  order  the  '  J  of G e a r ' s f o r m u l a s ,  = 1 and  interval  to e q u a t i o n (4.1) i s t h e n f o u n d  L  s > 2,  difference  step) e q u a t i o n of the  r— { u . - 4u. . + 3u. } + a?. u. 2h j j+1 j+2 j+2 j+2  G e a r (1971),  l ^  s p l i n e s of d e g r e e  = 1, z , = x. , _ a n d 1 j+ 2  0 < i < 2  1  i s one  z  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 - s t a g e  of S e c t i o n 2 g i v e s a t h r e e p o i n t (two  where  which  (  a  H u l m e ( 1971). )  d  The  Gaussian  t h i s r e s u l t d o e s n o t g e n e r a l i z e to h i g h e r o r d e r  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  procedure  i 2 ~  T h i s i s h i g h e r t h a n p r e d i c t e d by  i m p l i c i t Runge-Kutta processes.  (Weiss  j  p o i n t s w i t h r e s p e c t to the  c o r r e s p o n d i n g difference equations  m  h  t h a t the c o l l o c a t i o n p o i n t s i n t h i s  the G a u s s i a n  a f o r m u l a of o r d e r  z  2  It s h o u l d b e p o i n t e d out h e r e , approximation are  ( )  a  (See  G e a r ' s m e t h o d s a r e o b t a i n e d w i t h r = 0,  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. y  r = s = 1, i . e . , we  -i  w a n t to c o n s t r u c t a  compact  -25finite difference a p p r o x i m a t i o n  to the e q u a t i o n  (4.3)  1  L y ( x ) == y " ( x ) + a ( x ) y ' ( x ) + a°(x)y(x) = f ( x ) .  With m = l  a n d z^=x^  the c o e f f i c i e n t s  - l  d  are  U  j - l  +  d  0 j U  +  i n the d i f f e r e n c e  d  l  equation  j +l = j•  u  £  given by  d  2-h.,, a. ,1+1 .1 h h . h. ) •  , = L w (x.)  j (  d  = Lw  0 n  (x.) = y  •  +  + 1  r i+i i ,  2(h  h  )a  J J+1  J  and d,  = Lw  1  From  2 (x.) = z  2+h.a. J j  y V V V ^  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  difference approximation  i s e q u a l to o n e .  of c o n s i s t e n c y o f t h i s  The order  of c o n s i s t e n c y ~ 3  becomes (x-x.  two i f o n e l e t s  ^) (x-x.)(x-x. ^ ) .  coefficients  c o i n c i d e w i t h the i n f l e c t i o n p o i n t o f w  T h i s p o i n t i s g i v e n b y z ^ = x. + (h. ^-h.) / 3.  d^ of the c o r r e s p o n d i n g  finite  difference  J  a (z  ±(  J+1  (Zh.+h.^Jfh.+Zh.^) *0 ~ h.h., . J J+1  •2 +  The  equation a r e given by  1 ~ h.(h.+ h . ) J  (x) =  a (z,  -26-  1  h j  E x a m p l e 4.4.  +  1  ( h  A  j  +  h  )  2+  a (z )+  J  1  A  3  —  third order three point formula  obtained with m=3,  f o r e q u a t i o n (4. 3) i s  z, =x. ., z ~ = x . a n d z =x.,,. 1 J-1 2 j 3 j+1  difference equation a r e somewhat  ^ a ( z ^  The coefficients  c o m p l i c a t e d i n this case.  o f the  However, i f  the d i f f e r e n t i a l e q u a t i o n i s g i v e n b y  (4.4)  L y ( x ) = y " ( x ) + a°(x)y(x) = f(x) ,  t h e n t h e 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 e q u a t i o n i s i d e n t i c a l to t h e " 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. the m e s h i s a s s u m e d to b e u n i f o r m . difference equation  u. , ,- u. 1 _  u.- u. •)  —  4  . l -  1  + w ( A ) =w  J  (f),  where  W. w. = a.w. , + ( l - a . - b . ) w . + b.w. ,, , h J J J-1 J J J J J+1  with v. a  and  j  i ~ 6  There  I f t h i s i s n o t the c a s e t h e n t h e  becomes _J±J  (4:5)  ( E q u a t i o n (1. 7). )  2  i +l 6h.(h.+ h.__.) J J J+1 h  -272 h.  ,  b  For  a uniform  j  =  i _ 6  J  6h. (h.+ h. .) J+ J J+1 x1 1  m e s h the o r d e r  E x a m p l e 4. 5.  •  x  of c o n s i s t e n c y i s e q u a l to f o u r .  In this example, a number  p o i n t s z. w i l l b e g i v e n , corresponding  f o r w h i c h the o r d e r  finite d i f f e r e n c e equations  t h e o r e m (2.1).  of c h o i c e s  of the c o l l o c a t i o n  o f c o n s i s t e n c y o f the  i s g r e a t e r than p r e d i c t e d by  A g a i n 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 i s  g i v e n by (4. 3) a n d 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 compact with r = s = l.  The m e s h i s taken  F i r s t l e t m = 2.  This  L e t z ^ = X j + £^h,  two.  uniform.  gives a difference approximation  z ^ = X j + ^ h i a n d t ^ = x.. + r j h .  (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 are  a r e a s s u m e d to b e  t a k e n to be t h e r o o t s of p  of o r d e r  F r o m theorem  t h r e e i f the v a l u e s  o f e_ a n d e_  (£) = 0, w h e r e p(£) = (£-n)(£-l) £(£+1).  c o l l o c a t i o n p o i n t s 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. yields  z. = x. -  uniform  m e s h the o r d e r  Next, fourth order are Now  a n d z_ = x- +  equals  four.  T h i s w i l l i n g e n e r a l l e a d to a  The order becomes  f i v e i f the c o l l o c a t i o n  t h e i n f l e c t i o n p o i n t s o f w^(x) = ( x - t ^ H x - t ^ H x - t ^ M x - x j  z  l  3  =  X  =  X  "  j - |  j  +  i  ^  - eh  z  2  = x . _ i + eh .  z  4  = x  j + 1  +  eh  points  ^)(x-Xj )(x-x^ ^).  r e q u i r e t h e c o l l o c a t i o n p o i n t s to h a v e the f o r m  Z  This  B e c a u s e of s y m m e t r y a n d the  of this equation a c t u a l l y  c o n s i d e r t h e c a s e m = 4. formula.  .  The  +  -28T h i s choice has the advantage that the c o l l o c a t i o n points of consecutive difference  equations  p a r t l y c o i n c i d e , thereby l i m i t i n g the. .necessary  number of function e v a l u a t i o n s .  total  Upon some computation i t i s found that  these points a r e defined by the value  s  or  So w i t h either value of £ and w i t h z. as above one obtains a finite '  1  difference a p p r o x i m a t i o n for w h i c h the o r d e r of c o n s i s t e n c y i s at equal to five a c c o r d i n g to t h e o r e m ( 3 . 1 ) .  least  The a c t u a l o r d e r equals  six  because of s y m m e t r y and the u n i f o r m m e s h . E x a m p l e 4. 6.  In this example some five point finite difference  m a t i o n s to (4.3) w i l l be d i s c u s s e d . z^  =  x  j -  So let r = s = 2.  approxi-  F i r s t take m = l ,  with  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 equation i s equal  to t h r e e .  T h i s follows f r o m t h e o r e m ( 2 . 1 ) .  this equation i s i d e n t i c a l to (2. 5).  If the m e s h i s u n i f o r m then  B y t h e o r e m (3. 1) the o r d e r i s then  equal to four. If m = 5 and i f the c o l l o c a t i o n points z.., ( l < i < 5), c o i n c i d e w i t h the meshpoints,  then the o r d e r of the r e s u l t i n g equation i s equal to seven  i f the m e s h i s not u n i f o r m and equal to eight i f the m e s h i s u n i f o r m . c o u r s e the coefficients a r e v e r y c o m p l i c a t e d i n this c a s e . the d i f f e r e n t i a l equation does not i n v o l v e the a (x)y (x) t e r m , the equation i s g i v e n by (4.4),  However, i f i.e.,  (31u  + 128u  as  when  then the finite difference a p p r o x i m a t i o n i s  the same as the five point " 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), and c a n be e x p r e s s e d  Of  p . 502,  -29where  W , w - = (23w. _ + 688w. , + 2358w. + 688w.. , + 2 3 w . ^ ) / 3 7 8 0 . h J 3-2 J-l J J+1 J+ 2" o  A s i n example (4.5) one m a y w i s h to determine the s p e c i a l points z^, for w h i c h the finite difference consistency.  equation has an e x t r a h i g h o r d e r of  If m = 2 and i f the m e s h i s u n i f o r m , then a computation  s i m i l a r to those i n the p r e v i o u s example shows that the o r d e r i s equal to s i x i f one lets z  5.  1  4  = x . - £h, z „ = x . + £h, and i f the value of £ i s J 2 J  either  FINITE 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 CONDITIONS T o define a complete difference  boundary conditions a r e r e q u i r e d .  scheme,  a p p r o x i m a t i o n s to the  L e t a boundary c o n d i t i o n be g i v e n by  k-1 (5.1)  By(0) = y  with k < n.  ( k )  (0) +  2 i=0  b Y ( 0 ) = b,  T o construct a finite difference a p p r o x i m a t i o n to (5.1)  m a y p r o c e e d i n a fashion that r e s e m b l e s equations  ( l )  the c o n s t r u c t i o n of  one  difference  i n S e c t i o n 2.  L e t d = s + m and a s s u m e that a p o l y n o m i a l p(x) e P ^ i s u n i q u e l y d e t e r m i n e d by the  equations  (5.2)  p(x.) = u . ,  0 < i < s  (5.3)  L p ( ) = f(z.) ,  1 < i < m  Z i  -30T h e n the e q u a t i o n  (5.4)  generates  Bp(O)  a r e l a t i o n of the f o r m  (5.5)  B.u,, = ). d.u. = e „ b + h 0 . ^ 1 1 0 i=0  If the b a s i s f u n c t i o n s w ( x ) , 1  (Equations coefficients  (0<i<d),  u  v  a r e chosen  Bw (0)  Bw  JLw (zp  J_.w (z^)  d  i  =  (0)  determinants:  . . Bw (0) d  S + 1  . . J_.w (zj) C1  0 < i< s  E; Lw^z  '0  s + 1  a s i n S e c t i o n 2,  t h e n one c a n s h o w t h a t the  c a n be r e p r e s e n t e d b y t h e f o l l o w i n g  1  (5.7)  m ), e.f(z.) .. l r 1=1  of  (2.13) a n d (2.14) w i t h r = 0 . ) ,  1  (5.6)  = b  Lw  S  Lw  S +  +  m  1  )  JLw  ( Z j )  S + 1  (z  m  ) . . Lw^(z  Lw (z ) d  1  E,  '(z  m  )  Lw^(z  m  )  m  )  1  <3+  Bw  e  = 1  ^ E  0  d  (0)  S  Lw  (5.8)  Bw (0)  L w ( z ^) d  (z .)  Lw  S + 1  Lw  S +  Lw  S + 1  (z. (z  1  (z  l -  ,)  L w (z^  ^  )  Lw (z  )  1'  i + 1  1 d  ) m  Lw (z  i + 1  )  d  m  It i s convenient to define the n o r m a l i z i n g f a c t o r E  T w L  s  +  1  Q  (5.9)  E  0  = (-i)  (z\)  L w  {  N  by  (zj)  Q  m+1  L w 0  where  /  < i < m  S + 1  (z  m  ) '  '  L w (z d  0 A  L y ( x ) = y ^ ( x ) . If m = 0 define E Q  ) m  = 1.  Q  T h e proofs of the f o l l o w i n g  t h e o r e m s c l o s e l y follow those g i v e n i n S e c t i o n 2 . T h e o r e m 5. 1  A s s u m e that E Q ^ 0 .  m > 0 , let s + m > n . —  If m = 0 let s > k ;  otherwise if  T h e n at least k + 1 of the coefficients d. a r e l  z e r o for a l l s m a l l enough h , and the t r u n c a t i o n  error  m T  q  S  B y(0) - e b h  Q  -  I i  e  .f(z.)  =l  of the difference a p p r o x i m a t i o n (5.5) s a t i s f i e s  T  Q  < C h s+m-k+1  non-  -32-  i.e., Proof.  (5.10)  the o r d e r  of c o n s i s t e n c y of (5.5) i s at least equal to s + m - k + 1 .  C o n s i d e r the equations s J  d.vAx.) -  i=0  m YJ e L w ( z . ) = i=l i  i  e Bw (x ), i  0  0  0<i<d .  H e r e the p o l y n o m i a l s w (x) are the same as those used i n the of the c o e f f i c i e n t s .  T h e s e equations c a n a l s o be w r i t t e n  l  I  -A  0  -C  Bw"(x ) 0  0  = e, Bw (x ) d  0  m  w h e r e I is the (s+1) X (s+1) identity  L w V j  1  L w (z^)  matrix,  . . •  Lw"(z m  L w (z  m  )  )  and Lw  S +  1  (zj)  Lw  S +  '(z  m  C = L w ( z j) d  Lw (z d  )  rrr  )  as  definition  -33-  T h u s the  e q u a t i o n s (5.10) a d m i t a u n i q u e s o l u t i o n o n l y i f d e t ( C ) + 0. a s s u m p t i o n t h a t E Q + 0 i m p l i e s t h a t d e t ( C ) + 0,  p a r t i c u l a r , the that h i s s m a l l  In  provided  enough.  N e x t 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  LQ  and  BQ  r e p l a c i n g J_/ a n d  that l e s s than k+1 s u p p o s e t h a t d^  of the  #0,  B.  H e r e BQy(O) = y  c o e f f i c i e n t s d^ a r e  (1 < v < k^ < k+1),  and  nonzero.  (0).  Assume  More precisely  d^ = 0 o t h e r w i s e .  Let  v k  l  q(x) = "I ( 0 ) l _ ( ) ' w h e r e q^(x) e v,= l T h e n q j x j e P ^ , and by a s s u m p t i o n k < d . x - x  c  ^ has l e a d i n g c o e f f i c i e n t l / k ! . "1 Therefore  x  v  s  m  Y  - Y  ^(x.)  i=0 The  ^LQqCz.) =  1.  H e n c e eQ = 0.  small  the  same can  be  the  But  T h u s i f E _ j. 0 t h e n at l e a s t k + 1 0 is  5. 2.  The  the l a t t e r c a s e the  As  0  entire lefthand side vanish.  already  then i t follows f r o m  normalizing  of the c o e f f i c i e n t s d. a r e i  The  p r o o f of the  one  w i l l be  nonzero.  This  proves  takes  .  the  second a s s e r t i o n c l o s e l y omitted.  i s s a i d to b e  E Q = 1 if m  m + 1  If h  f a c t o r i s n o n z e r o i f m = 0 o r i f s = k.  difference equation  stated,  Also,  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 )  f o l l o w s that of p a r t of t h e o r e m (2. 1) and  Proof.  0  s a i d of the c o e f f i c i e n t s d..  f i r s t s t a t e m e n t of the t h e o r e m .  Theorem  0  i=l  c o n s t r u c t i o n of q(x) m a k e s  BQq(xQ) =  e B q(x ).  (In  compact. )  = 0.  E Q = 0,  If s = k and  (5.9)  t h a t 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 = V c.w (x) s a t i s f i e s L q ( z . ) = 0, ( l < i < m ) . B ut 1  that the p o l y n o m i a l  q(x)  1  ft  1=  LQ^X) =  0.  e  P s  So  +  m  _  a n d n  q(x) = 0.  s+ m This  8+1  - n< m  because s = k<n.  however is impossible,  Hence L q(x) = Q  since not a l l c  are  q (x) W  equal  -34to zero and since the polynomials  As  w (x) a r e l i n e a r l y independent. 1  was the case for finite difference approximations to the  differential equation,  i t i s p o s s i b l e to identify collocation points z^,  which the approximation of consistency,  to the boundary condition attains a higher  than p r e d i c t e d by theorem (5.1).  for  order  A m i n o r m o d i f i c a t i o n of  the proof of theorem (3.1) shows the following.  Theorem E  Q  5. 3.  L e t the coefficients d. and e. and the n o r m a l i z i n g factor  be given by equations  (5.6),  E .0 + 0 and that s + m > _ n.  (5.7),  Suppose r-r- w * S  (5.8) and (5.9). m +  ~ s+m+1. . _ w (x ) = 0 ,  A s s u m e that  ^ ( x ) e P s+m+1 satisfies  _ ^.^ 0< I < s,  i  and j k ~ s+m+l , d w (0) _ /ri  A  '  k  dx A l s o a s s u m e that the collocation points a r e chosen such that L Q W (1 < i < m).  T h e n the order  S + m +  * ( . ) =0, z  of consistency of the finite difference approxi-  mation (5.5) i s at least equal to s + m - k + 2.  6.  EXAMPLES  Example 6.1.  OF APPROXIMATIONS  L e t n = 2 and k = 1.  TO BOUNDARY  CONDITIONS  Thus the differential equation i s  given by  (6.1)  and  Ly(x)E= "(x) + y  the boundary condition by  a^xjyV)  + a°(x)y(x) = f(x),  -35-  (6. 2)  Let  m  By(0)  = 0,  involve  so  0  that 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  the d i f f e r e n t i a l  procedure  = y (0) + b y ( 0 ) = b.  equation  of the p r e v i o u s  (6.1).  Q  = Bw  (0) =  d  1  = Bw  (0) = B  not  W i t h s = l the c o n s t r u c t i o n  section yields  d  to (6. 2) d o e s  e^ = E Q = 1 ,  B  = ~h  V  +  and  The  approximation  w h i c h has  order  w e l l - k n o w n 3-point  u  r o  —  h ~ ~  +  b  o o u  to  is  second o r d e r  compact  -(x-x.)/h,  Now  and w  1  h  b ]  one.  t h e n one  obtains an e q u a l l y  equation, v i z .  1 +  3u  2h  6.2.  =  i f the m e s h i s u n i f o r m ,  2  Example  1""  u  u -4u  which is a  1  )/h  is therefore  of c o n s i s t e n c y e q u a l  If s = 2 and  (x-x  Q  + b u Q  = b,  formula.  let m = l  and  s =l,  so  that the a p p r o x i m a t i o n  t h e r e f o r e . c o n s i s t e n t f o r any  (x) = ( x - x ) / h a n d n  0  w  Z  c h o i c e of z^.  (x) = (x - x ) ( x 0  X l  )/h  .  to (6. 2)  L e t w^(x)  T h en  the  =  -36-  2  normalizing  factor  i s given  by E ^ = L ^ w  2  (z^) = — h  .  The coefficients  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  (6.3)  are  d  o o u  +  d  l l u  =  e  o  b +  e  l ( l^ ' f  z  1  +I  t h e n e q u a l to  e  o = ET  L w 2  0  Bw°(0) d  0=E"  ( I> = z  a l  ( i) z  Bw (0) 2  =  T  -  a  (  Z l  ) +b  Q  + 2-b a Q  (zp,  0  Lw"(zj)  Lw^(zj)  B w V )  Bw (0)  and  2  k + a (z ^ Lw  The  order  of c o n s i s t e n c y  (z  for which choice  to  From  theorem ~3  inflection point of w  Example to  6. 3.  to two.  It i s now  o f i n t e r e s t to  o f the c o l l o c a t i o n p o i n t z^ the o r d e r  ( 5 . 3 ) i t follows 2  (x) = x (x-x^).  L e t s = 2 and m  ( 6 . 2 ) h a s the f o r m  L w (z^)  i s equal  determine three.  1'  + ^ a ( z ^  = 1,  This  increases  that t h i s  h a p p e n s i f z. i s the x h 1 1 p o i n t i s z^ = — = — .  so that the d i f f e r e n c e a p p r o x i m a t i o n  -37-  (6.5)  d u Q  + dn  0  i  + d u  i  2  A s s u m e that t h e m e s h i s u n i f o r m  = e b + e ^(z ^ .  2  Q  a n d l e t z ^ = X Q + £h.  Then E Q =  3  LQW z^  ( X Q + |h) a n d E Q = 0 i f a n d o n l y i f £, = 1.  = x  i»  The order  point of w or  of c o n s i s t e n c y i n c r e a s e s  (x) = x ( x - x ^ x - x ^ ) .  H e n c e one should  not l e t  to f o u r i f z ^ i s a n i n f l e c t i o n  T h i s h a p p e n s i f i = (9 - \ / 3 3 ) / l 2 = 0. 27  % =(9 + N/"33)/12 = 1. 23.  7.  THE STABILITY OF FINITE DIFFERENCE  I n t h 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  SCHEMES.  o f K r e i s s (1972) w i l l b e u s e d t o  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 s c h e m e s b a s e d upon c o n s i s t e n t finite difference approximations  introduced i n previous  p u r p o s e t h e m e s h i s a s s u m e d t o be u n i f o r m , . .•  .  A  - . c.  . '..  , . .'  d. .u. ,. = V e. . f ( z . .) = f., r < j < J - s , J , i J+i u j , i j,i' j i=l v  J  with boundary conditions s (0)  m (0)  k  (7.2)  B  h ) k  k  (0)u ,£  d  0  k )  .(0)u. e =  k ) 0  (0)b (0) k  +  i=0  2  e  k  3)  B  n  t  ^ J  s  . l  . (0)f(z  k >  . (0))  l<k<n , 0  0 7  fc>  i=l = b (0),  <'  of the f o r m  m  L, u. = V h j U i=-r  together  F o r this  i . e . , h^ = h = l / J , (1 < j < J ) .  d i f f e r e n c e s c h e m e c o n s i s t s of equations s  (7. 1)  .  sections.  m  k,i< ) J i i=-r (l) d  | n  1  u  +  =  e  k,0< V Z 1 )  1 ) +  k  { i )  1  ^ . i ^ k , ! *  i=l  k  b, ( 1 ) ,  1< k < n -n  r  1  "  -38In addition,  i f (7. 1) i s n o t c o m p a c t ,  boundary conditions  are required  e q u a t i o n s a n d the n u m b e r w i l l be a s s u m e d differential  v  A l t h o u g h this i s not n e c e s s a r y , i t  equations  a r e a l s o c o n s i s t e n t w i t h the  m.  V  L, u. = h j  d. .u. =  u  j , i i  L . u. = J  X u  e. . f ( z . .) = £., J,I J,I'  k >0, 0 —  k„ < j < r - l ,  y  X  0 -  J  n  -  i=l  0 H  of unknowns.  to m a t c h the n u m b e r o f  equation and given b y  i=0  (7. 5)  i n order  that t h e s e e x t r a  s.  (7.4)  i . e . , i f r + s > n then r + s - n e x t r a  j  m  1 j,i J i i=-r. d  U  +  J  i{z  e  V  55  = I ),i j,i i=l  ]  J"B+I<j<  J-n+k . Q  T Let  Uj. = (UQ,  .  '",UJ)  Then  (7. 1) t h r o u g h (7.5) c a n b e  the e q u a t i o n s  compactly written as (7.6)  Here  L  f^ e R  J  +  = f  h  i s the a p p r o p r i a t e  X  (J+1) X (J+1) matrix. easily  u  h  h  .  right hand  side vector  and L ^ is a  Consistency of the equations (7.1), (7.4) and (7. 5) i s  s e e n to i m p l y  that  ~ cn. f.-f(x.) I <. ch. J y T  AA lI sS oO,,  J+ e,^ e R i le ect e  1  b e the  1  e r r o r vector, i . e . , e^ = (e , " ' , e j ) Q  vector of truncation e r r o r s , T ^ D , • • • , T l n - n ( i : ) ) , with  T  = n  with e^ = y. - u.., and let  ( _(0)> "" > T  T  T n  k ' " ' J-n+k ' T  T  Q  T  T  k<  0 ) 55 B  h,k( y 0 )  j = h j " V L  Y  0  "  V°>«  1 < k<  V  k <j < J -n + k , Q  Q  be the  -39and  s  For  w, e R  J  +  1  B  h,k<  1 ) y  - V ^  J  1 < k <n - n  1  l e t I lw, j j = m a x |w. |. 0<j<J  Q  If A, i s a (J+1) X (J+1)  matrix  J  I  I  i  I I  l  l  !  ^  ^ l l ^ l l  ^ ^  then I j A, | i s the induced matrix norm, i . e . , | |A, | | = max w *0 h  The  . ||w || h  f i n i t e d i f f e r e n c e s c h e m e (7.6) i s s a i d to b e c o n s i s t e n t i f | |T^| | < k ^ h _i  as  h - * 0 and stable i f  | | T ^ | j < k^. The  exists for a l l s m a l l  H e r e k^ a n d k^ a r e c o n s t a n t s  stability property  t h a t do n o t d e p e n d o n h.  i s e s s e n t i a l l y d e t e r m i n e d b y the d i f f e r e n c e  m a t i o n t o the h i g h e s t d e r i v a t i v e . characteristic polynomial S e c t i o n 2, ( E q u a t i o n s  enough h a n d  T o be m o r e  o f (7. 1) i s n e e d e d .  (2.12) a n d (2.16)), c(t) ^ h  n  I i= -r  d.t  r + i  specific,  approxi-  the n o t i o n o f  L e t E a n d d^ be a s i n  and define .  It i s n o t d i f f i c u l t to s h o w that the e q u a t i o n  c(t) = 0 m u s t h a v e a r o o t t = l  of m u l t i p l i c i t y n, b e c a u s e o f c o n s i s t e n c y .  Hence  (t-l) c(t), n  where  c(t) i s d e f i n e d  one c a n w r i t e c(t) =  to b e t h e c h a r a c t e r i s t i c  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 e q u a t i o n  (7. 1).  polynomial  A s s u m e t h a t c(t) i s  explicitly given b y  (7.7)  c(t) =  N V i=0  a.t , 1  with  N = r + s - n.  1  If the f i n i t e d i f f e r e n c e s c h e m e (7. 6) i s not c o m p a c t , r + s - n characteristic polynomials c o n d i t i o n s (7.4) a n d ( 7 . 5 ) .  then one a l s o h a 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  It i s a s s u m e d  that t h e s e h a v e t h e f o r m  -40N. J  (7.8)  c.(t) = £  a.^.t  1  ,  k  ,  J - s +  Q  < j < r - 1 ,  i=0 and  (7. 9)  c.(t) =  J  a  t  1  1  <  j <  J - n + k„.  i=0 Finally,  consider  the h o m o g e n e o u s d i f f e r e n c e  equation  N  YJ V j + i  (7. 10)  =  ° '  0 < j < oo.  i=0 with boundary  (7. 11)  conditions  N. J Y - i . . ' 1=0 a  v  J  1  = 1  °» Q ! J < k  r  "> x  S U  P l - I < c o n s t . , 0 < j < oo , 0< j < oo —J — v  J  r  and N (7.12)  Y N-i J-j-i i=0 a  V  =  °'  0  < J  <  o  o  >  s u b j e c t to  (7.13)  N. J  Y  a  i  ^  N  i=0  . = 0, J - s + l < j < J - n + k ,  j  Q  sup 0<j<co  T h e n one c a n s t a t e t h e f o l l o w i n g t h e o r e m due to K r e i s s T h e o r e m 7.1.  Ivj.Jl (1972).  L e t the h o m o g e n e o u s p r o b l e m c o r r e s p o n d i n g  (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 .  const.  to (1.1) a n d  A s s u m e that the d i f f e r e n c e s c h e m e  (7.6) i s c o n s i s t e n t a n d 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 e q u a t i o n c(t) = 0  -41satisfy  t. + 1. l  1  F u r t h e r , i f the d i f f e r e n c e s c h e m e i s n o t c o m p a c t t h e n a l s o that the d i f f e r e n c e equations d i f f e r e n c e equations  (7.10) w i t h b o u n d a r y c o n d i t i o n s  assume  (7.11) a n d the  (7. 12) w i t h b o u n d a r y c o n d i t i o n s (7. 13) o n l y h a v e t h e  trivial solution.  T h e n (7.6)  has a unique s o l u t i o n for a l l s m a l l enough h a n d the  difference scheme i s stable.  We (7. 1).  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  If equation  approximation  (7. 1) i s c o m p a c t ( a n d c o n s i s t e n t ) t h e n c o n s i s t e n c y o f  the b o u n d a r y c o n d i t i o n s  (7.2)  a n d (7.3)  i s s u f f i c i e n t to g u a r a n t e e  If (7.1) i s n o t c o m p a c t t h e n 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 I"!" S — II  s y m m e t r i c i f c(t) = t  motivate  five point formula approximation v v  c(t) s a i d to b e  1  c(—) a n d 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 i f the  d e g r e e o f c(t) i s e v e n a n d i f | To  stability.  a  r+s -n g| _^ l jj>  J  a  the l a s t d e f i n i t i o n ,  ( H e r e I = (r+s-n) /2. )  consider  the g e n e r a l  f o r m of a  that i s c o n s i s t e n t w i t h the s e c o n d d e r i v a t i v e .  This  2 2 2 c a n b e w r i t t e n a s K. D. u. w h e r e D, u. = (u. , , - 2u.+ u. ) / h h h j h j j+1 j 4  and  K,w. = a „ w . , + a.w. + a w . . . h j 0 j-1 1 j 2 j+1 2  c(t)  = a ^ + a ^ t + a ^ t a n d c(t) i s d i a g o n a l l y d o m i n a n t i f a n d o n l y i f t h e  0  operator  (matrix)  i s diagonally  The c h a r a c t e r i s t i c  equation i s >1  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 s u 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  w i t h 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 (7.4) Lemma  7.2.  Assume  t h e n c ( - l ) = 0.  a n d (7.5)  that c(t) i s s y m m e t r i c .  equations  (7. 1)  associated  n e e d n o t be.) I f t h e d e g r e e o f c(t) i s o d d  -42If the d e g r e e of c(t) i s e v e n a n d 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 with positive coefficients,  t h e n t h e r e a r e no r o o t s o f c(t) = 0 o n  the u n i t c i r c l e .  Proof. c(t)  T h e p r o o f of the f i r s t s t a t e m e n t i s i m m e d i a t e .  equals  2N^,  then  2N l ( ) | c  =  e  L k=0  N,-l a  = l( N 1  e  a  k  '  Example  N  7. 3.  h k k=0  +  a  v  N  = I a,  I f the d e g r e e o f  f  (  e  +  l  e  }  N  + 2 V k^O  T  l  .  k  a. c o s k 0 I > a . l  e  I  V k  N -1 A  a, c o s k 9 I > a -2 l k^O M  N  I n the s p e c i a l c a s e w h e r e c(t) = a ^ t  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^ <  fi  -i  - 2I k^O  T  N  )  N  A  V a. > 0.  k  + ( l - 2 a Q ) t + a ^ the  simple  computation shows  t h a t i n f a c t t h e r e a r e no r o o t s o n the u n i t c i r c l e i f f - co < a ^ < s h o w s that the a s s u m p t i o n s i n the l e m m a a r e n o t s t r i c t l y  This  necessary,  although perhaps d e s i r a b l e .  The  e f f e c t o f the c h o i c e  the r o o t s o f t h e c h a r a c t e r i s t i c m e a n s o f the n e x t  of the c o l l o c a t i o n p o i n t s z^, (1 < i < m ) , o n e q u a t i o n c(t) = 0 w i l l b e i l l u s t r a t e d b y  example. 2  Example  7.4.  symmetric  L e t c(t) = a ^ t + ( l - 2 a Q ) t + aQ a n d l e t n = 2, i . e . , w e c o n s i d e r  5-point a p p r o x i m a t i o n s  be s y m m e t r i c i t i s n e c e s s a r y i n the i n t e r v a l |. j - 2 ' j+2. X  z^ = Xy by y  X  to the s e c o n d d e r i v a t i v e .  to p l a c e the c o l l o c a t i o n p o i n t s s y m m e t r i c a l l y H e n c e 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  The difference approximation  to y  (x) = f(x) i s i n t h i s c a s e  f - u . + 1 6 u . .-30u.+ 16u., . - u. ,~^ / 1 2 h = f ( x . ) . \ 3-2 J-1 J J+1 J+2/ ' y 2  0  a l s o we w r i t t e n as  F o r c(t) t o  v  given  The left hand side c a n  1_ 12b  2  s o t h a 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 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 b y 1 c(t)  2  = - J2" t  7 ± 4 \/T,  14 +  1 t - Yz  '  c  h  a  r  a  c  t  e  r  i  s  tic  equation has  so t h a t | t | ,£ 1.  Next,  i f m = 2 l e t z, = x. - q h a n d z_ = x. + qh, q > 0. 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 u._ 0  2  roots L =  + ( l ^ a ^ u . ^ + (-2 + 6 a ) u . + ( l - 4 a ) u . 0  0  + 1  + a u. Q  + 2  )/h  2  This  = ±(f(  z  generates  J + f(z )) . 2  2 T h e 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 - 2 a Q ) t + 2 It i s e a s i l y  s h o w n that a^ and  e x a m p l e (7. 3) t h e r e a r e no  q a r e r e l a t e d by a  Q  = (6q - 1 ) / 12.  aQ.  By  r o o t s of c(t) = 0 o n the u n i t c i r c l e i f a n d  only i f  -co < Q < 4I* t e r m s of q t h i s c o n d i t i o n b e c o m e s 0 < q < ~ 0. 816. If 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^-q^h a n d z = x^ + qQh, whe r e a  1  2  q_ = J10 of  \j-rp ~ 0. 38 o r q„ = Ji 15 U  1. 36,  t h e n the o r d e r of c o n s i s t e n c y  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 .  e x a m p l e (4.6).) for  + N/T? « 15  T h e r e f o r e the t h e o r y of K r e i s s  the f i r s t c a s e ,  w h e r e qg = sj 1 - ^Y^  •  guarantees  Of c o u r s e ,  scheme i s supported by  That such an a p p r o x i m a t i o n may  number  give a better e r r o r constant, take m  v a l u e of q ^ i s  l e a d to a  the d a t a o b t a i n e d i n the n e x t s e c t i o n .  e x a m p l e (8. 1), e x p e r i m e n t  Finally,  stability only  t h i s d o e s not i m p l y  t h a 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 b a s e d u p o n the s e c o n d n e c e s s a r i l y unstable.  15.)  The  (See  stable  (See  f i r s t v a l u e of q^ a p p e a r s to  however.  = 5 and  l e t the c o l l o c a t i o n p o i n t s z^,  c o i n c i d e w i t h the m e s h p o i n t s j j » x  +  (-2<i<2).  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  (1 < i < 5),  This generates  ( I 9 6 0 ) , p a g e 502, v i z .  the f i v e  -44-  1  +  128u J-1  2 52h'  318u. + 128u., + J J+1  1  3f. , J"  3780  The  associated  to  stable  a  t^ ~ - 0 . 1 6 8  (31t  that  that  degree  strictly istic  outside  the  unit  as  the  (7.  (7.10),  that  the  number  s + n)/2 the  stability  conditions it  and of  and  and  exactly  Hence  0  ,  these  conditions  A  equal  extra  c(t) to  positive N^  collocation points  roots  (7.5)).  N^  It  equations  of  is the  conditions  addition  then  the  unit  circle  the  admit  (7.5).  lead  for  the at  and  c(t)  is  character-  the  and  N^  difference  zero  solution  x = 0  is  the  = (r  + s - n) / 2 .  also  reasonable  to  assume  equations  at  x = 0  extra  Again,  symmetric  in  condition to  and is  If  inside  (7.13)  to  ( 7 . 1)  coefficients  boundary  equal  of  (7.4)  2N^.  necessary  (7.12),  x = 1 and  °J+r-s-j  the  =  is  (Hence  only same  =  and  now x = 1  by  14)  i.e.,  31)  polynomial  with  circle.  number  (7.4)  +  boundary  even  has  characteristic  related  is  (7.11)  at  in  extra  dominant = 0  then  are  c(t)  c(t)  is  that  of  190t  ~ -5.96.  characteristic  diagonally  equations  -  the  equation  roots  (r  the  t^  +  2  equation  approximation.  consider  assume the  and  difference  Now  1  characteristic  1  roots  + 688f. . + 2358f. + 6 8 8 f . , . + 23f. _ V J " J J+1 J+2/  2  252  with  31u  t  is  then  (t)  at  =  t  J  c (-),  x = 1 are  sufficient  to  (r+s-n)/2  the show  = k  "reflection" that  the  Q  < j < r - l  of  those  difference  ,  at  x = 0.  equations  For (7.  10)  -45s u b j e c t to (7.11) h a v e t h e z e r o  solution only.  F o r this p u r p o s e  define  p o l y n o m i a l s p-(t) b y N p (t)  (7. 15)  a t^ ,  I  EE  1  j > 0  i=0 H e r e N = r + s - n a n d t h e c o e f f i c i e n t s a^ a r e t h e s a m e characteristic polynomial  c(t) i n ( 7 . 6 ) .  a s those of the  I f (7.10) s u b j 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 s e e n that the p o l y n o m i a l s ((r+s-n) / 2  =k  dependent.  H e n c e we  Theorem  (7. 5).  Q  < j< r-l),  and  < j < m a x N.-N), a r e l i n e a r l y  h a v e s h o w n the f o l l o w i n g  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  ( l . l a , b ) only have the t r i v i a l be consistent.  p^(t), (o  c^(t),  Assume  to (1.1) a n d  s o l u t i o n a n d l e t the d i f f e r e n c e s c h e m e (7.6)  t h a t c(t) i s s y m m e t r i c .  Also  suppose that the  d e g r e e o f c(t) i s e v e n a n d that 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 positive coefficients. boundary conditions  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 be r e l a t e d a s i n E q u a t i o n ( 7 . 1 4 ) .  c..(t) a n d P j ( t ) d e f i n e d  s t a n t K i n d e p e n d e n t of h s u c h that  Example  7.6.  of the e x t r a  I f the p o l y n o m i a l s  b y (7. 8) a n d (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  pendent then the d i f f e r e n c e s c h e m e i s s t a b l e .  with  Hence there  inde-  exists a con-  | | e ^ | | < K | | T ^ | |.  2 L e t c(t) = a ^ t + ( l - 2 a Q ) t + a ^ w i t h  d e g r e e of 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  - co < a ^ <  c^(t) of the e x t r a  1  If the  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 t h e n t h e d i f f e r e n c e s c h e m e i s s t a b l e i f c.^(t) ^ c ( t ) .  I f t h e d e g r e e of c ^ ( t ) i s t h r e e  then stability i s guaranteed 3  i f t h e r e i s no c o n s t a n t  b s u c h that c , ( t ) = b ( a t n  2 + ( l - 2 a ) t +a„t)+(l-b) n  8.  NUMERICAL  The  EXAMPLES  m a i n p u r p o s e of the n u m e r i c a l  i s to c h e c k t h e c o r r e c t n e s s o f s t a t e m e n t s  examples given i n this section i n previous  sections.  a l s o g i v e s o m e i n d i c a t i o n a s to w h a t t h e r e l a t i v e a c c u r a c y discretizations u s i n g double  is.  A l l computations  precision arithmetic.  e f f i c i e n c y o f the c o m p u t a t i o n s ,  They  of v a r i o u s  w e r e c a r r i e d out o n a n I B M 370/168,  No a t t e m p t  w a s m a d e to o p t i m i z e the  so t h a t t h e r e w i l l be no c o n c l u s i o n s a b o u t  the r e l a t i v e m e r i t o f v a r i o u s f i n i t e d i f f e r e n c e s c h e m e s .  Example  8.1.  L e t the d i f f e r e n t i a l e q u a t i o n b e g i v e n b y  L y ( x ) = y " + y' - 2 y = 2 ( l - 6 x ) e  X  with boundary conditions  y(0) = y ( l ) = 0 ,  The  s o l u t i o n of this 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 s o m e n u m e r i c a l  this table r, s and m So  computations  a r e g i v e n i n T a b l e (8.1). I n  a r e a s 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  the w i d t h o f the a p p r o x i m a t i o n  n u m b e r of c o l l o c a t i o n p o i n t s .  e q u a l s to = r + s + 1 a n d m  T h e l e t t e r s A, B ,  C,  (2.4).  denotes the  D a n d E i n the  c o l u m n s h e a d e d b y " c " a r e a c o d e i n d i c a t i n g the l o c a t i o n o f t h e c o l l o c a t i o n points, v i z .  A: T h e  p o i n t s 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 and  m  B: T h e  of r , s  considered.  p o i n t s z^,  (1 < i < m ) ,  a r e o p t i m a l u n d e r the r e s t r i c t i o n s  that  e a c h s u b i n t e r v a l c o n t a i n the s a m e n u m b e r of c o l l o c a t i o n p o i n t s , t h a t t h e s e p o i n t s 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 interval. C: T h e  (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  and sub-  5 of S e c t i o n 4. )  c o l l o c a t i o n p o i n t s 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..  J  1  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 , b u t t h e y  placed  The  s y m m e t r i c a l l y i n the  c o l u m n headed by  "o"  difference approximation Columns  7 through  r <_ j < J - s .  g i v e s the o r d e r  as p r e d i c t e d by  be  (The  0.438  -1  observed  of the one  maximum  asymptotic  chapter. (2.4) f o r special  T h i s i s done i n  for j = J - l i s assumed m e s h i s taken u n i f o r m  a n u m b e r of v a l u e s  to  in  of J  the  The  notation  J  error,  m e a n s 0.438 10  For  The  finite  i s e q u a l to 5 t h e n a  d e f i n e d f o r j = 1.  f o r j = 1.)  so t h a t h. = h = 4.  scheme is listed. order  t h e o r e m s i n this  s p e c i a l equation n e c e s s a r y  J observed  , x.,  of c o n s i s t e n c y of the  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  the " r e f l e c t i o n "  this example,  x.  11 d e f i n e the f i n i t e d i f f e r e n c e e q u a t i o n s  f i n i t e d i f f e r e n c e e q u a t i o n m u s t be c o l u m n s 2-6.  interval  are  -1  order  i . e . , max etc.  | u . - y ( X j ) |, i s g i v e n .  I n the f i n a l c o l u m n , h e a d e d by  of a c c u r a c y  order  the  of the g i v e n f i n i t e d i f f e r e n c e  If i t i s n o t c l e a r f r o m the n u m e r i c a l  i s t h e n the e x p e c t e d  "a",  i s given between  r e s u l t s what this  brackets.  M o s t of the r e s u l t s t h a t a p p e a r i n the t a b l e a r e The  first seven experiments  For  experiments  5 and  6 the c o l l o c a t i o n p o i n t s a r e  given by  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  reached.  (See  the l a s t c a s e  ~ 2  45~ '  u s e a  Experiments  - ^  8-12  n  E x p e r i m e n t 5, s h o w the  to  and  These values are  £ = V -j-^" + 2"  45~ '  then t h e o r d e r  the o v e r a l l a c c u r a c y .  of a c c u r a c y  H u b b a r d (1964) and  of the s c h e m e r e m a i n s f o u r .  S h o o s m i t h (1975). )  The  experiment  considered  exceptional.)  equations  13 a n d are  = X j + £h.  z ^ = x.. - £h and  of c o n s i s t e n c y b e c o m e s s i x . 1 + N TT  Example  8. 2.  6  equal  also  for  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  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  e =  i n  actual accuracy  t h i s i s not the c a s e  In  d  This  (See  (Amazingly  b u t t h i s m u s t be  e  i s only  however i s seriously effected. 9.  s  N o t e t h a t e v e n i f the  p h e n o m e n o n 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). B r a m b l e and  u  e f f e c t t h a t v a r i o u s c h o i c e s of the e x t r a  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 e q u a t i o n s  two,  z. = x. • i_ ± £h. 1  of c o n s i s t e n c y i s  d i s c u s s e d i n E x a m p l e (4. 5)).  b o u n d a r y conditions have on order  explanatory.  involve compact difference approximations.  There are  £ = N Y2  self  used i n  These are  Again,  y  11  the  + xy  £ = v 1 - v T T u s e d i n 13  equation  1  - (1 + x ) y = -(2 + x ) ex  with boundary conditions  y'(0) =  y(l)  =  0  been mentioned  v a l u e s of £ f o r w h i c h the o r d e r  14.  Consider  as h a s  and  -49T A B L E 8. 1  r sm  #  c o  J = 8  J = 4  1  1 1 2 D 2  .438"  1  .111"  1  2  1 1 2 B  2  .111"  1  .279"  2  3  1 1 2 A  4  -4 . 830  4  1 1 3 C 4  .516"  3  ' 5  1 1 4 B 6  .213~  6  6  1 1 4 B 6  .704"  6  7  1 1 5 E 6  . 156"  1 1 1 D 2  2 2 1 D 4  .291"  9  1 3 1 D 3  2 2 1 D 4  .440"  10  1 4 1 D 4  2 2 1 D 4  11  1 1 2 A 4  2 2 1 D 4  -3 . 533  12  1 1 3 C 4  2 2 1 D 4  .884"  13  1 1 5 E 6  2 2 2 A 6  .571  14  1 1 5 E 6  2 2 2 A 6  -4 . 260  1  2  3  _7  J = 64  a  .281"  2  .704"  3  .176"  2  .700"  3  .175"  3  .337~  6  .211"  7  -5 '. 203  .337  .528"  - 8  . 106"  8  J .= 32  -5 ° -4 . 325  . 534  6  J = 16  .249"  8  .265~  2  7  -3 . 233 -4 . 369 _4 . 719 -4 . 781 ' . 185" .745"  6  9  .390"  1 0  .204"  3  .798"  6  - 12 C  : 257" .628  2 8  4 4  8  **  6 6  1 1  -12  6  L  4  .147"  4  .987"  -5 . 897  .310"  6  . 160"  .458"  .326"  6  . 207  .336"  6  . 212  6  . 212  4  * *  6  5  _5 .518  8  '.817  -4  . 132"  . 128"  1 0  . 165"  ' .439  3  .528~  5  .338"  .367"  1 0  .601  -7 . 122  6  7  _7  4  _ 7 _7  - 12  . 160"  L  9  .215"  C o n t a m i n a t e d by roundoff  (4)  1 1  4  (6)  -50The  solution to this p r o b l e m i s y(x) = (l-x)e .  In this example only compact approximations to the d i f f e r e n t i a l equations  a r e considered.  N u m e r i c a l test calculations a r e p e r f o r m e d with i  v a r i o u s difference approximations to the boundary condition y (0) = 0. Results appear i n Table (8. 2). previous  example.  i n columns 2-6,  The  The  notation used i s the same as i n the  approximation  to the boundary condition i s defined  while the finite difference approximation  to the d i f f e r e n t i a l  equation i s defined i n the next five columns. In E x p e r i m e n t 9 the collocation point for the boundary condition i s z^ = X g + y. Example  For  (6.2).)  XQ  + £,h.  (9  + N/33)/12.  this value of z^ the o r d e r In E x p e r i m e n t s  In 13 the value For  of £ i s (9  of a c c u r a c y  14 this collocation point i s z^  - N/ 33)/12 and  d i s c r e t e boundary condition.  Example  8.3.  The  Note that the  of consistency of the  T h i s d i f f e r s f r o m observations  the extra boundary conditions i n E x a m p l e  made about  (8.1).  purpose of this example i s to i l l u s t r a t e the  of allowing nonuniform meshes.  8 y" 2  L e t the d i f f e r e n t i a l equation be  - x y= 2  2s  2  - x  4  - ee"  X  /  2 s  ,  -1 and  y(l)  = 1 + e  2 S  «  usefulness given  £ = 0.001,  with boundary conditions  y'(0) = 0  =  of consistency i s  (See E x a m p l e (6. 3). )  i s not greater than the o r d e r  (See  i n 14 this value i s  these collocation points the order  equal to four rather than three. order  13 and  of consistency i s three.  1.  by  -51-  TABLE  r s mc  #  o  r s"m c  o  8. 2  J = 8'  ' J = '4  J  =16  J = 64  J = 32 . 108  1,  0 10 - 1  11 1 D 2  . 117°  .494-1  .226"  2  0 2 0 - 2  11 1 D 2  .240"  .492"  2  . 110"  2  .261"  3  .633  3  0 1 1 * 2  11 1 D 2  -1 . 250  .631"  2  . 158"  2  .396"  3  .989  4  0 1 1 **2  11 1 D 2  .240"  .492"  2  .  .261"  3  .633  5  0 10 - 1  11 2 A 4  .106°  6  0 2 0 - 2  11 2 A 4  .362"  1  .786"  2  7  0 1 1 * 2  11 2A 4  .152-  1  .359"  2  8  0 1 1 **2  11 2 A 4  .338"  1  .764"  2  9  0 1 1 A 3';  1. 12; A 4. .531" '  -4 . 60.2 *  10  0 1 2 C 3  11 2 A 4  . 146"  . 171  11  0 4 0 - 4  11 2A 4  .344"  12  0 1 3 E 4  11 2A 4  .482  13  0 2 1 A4  11 2 A 4  . 165"  14  (0 2. 1 A^4 I  11 2A 4  .563"  1  1  .466"  3  2  3  0  3  _4 . 249  5  2  3  . 181"  - 1  . 106"  1  .876"  -3  . 108"  2  . 183"  -5 . 324 4  -  .219"  1  . 164"  2  n o  1  2  - 5. 71'! -4 . 207  1  .521"  2 1  2  . 109"  3  -4  3  . 108"  3  .862"  6  . 106"  6  .254"  5  .316"  6  . 132  -7  .472"  8  -7  .734  .538  6  _7 .317 .828  -9  2 2  .440"  6  l  . •2  3  .208"  z  -4  .216"  .526"  l  2  _4  -3  5  z  1  2  _4  .443  .898"  _7 . 730 -5 . 129  .526"  a  2 3 3 4 4  .300"  9  4  .437"  8  4  -52As  was  the c a s e i n the p r e v i o u s e x a m p l e s , t h i s p r o b l e m h a s b e e n c o n s t r u c t e d  so i t w o u l d h a v e a k n o w n s o l u t i o n .  H e r e the s o l u t i o n i s  . . -x /2c 2 y(x) = e ' + x , 2  w h i c h has a b o u n d a r y l a y e r at x = 0.  One  could equally well impose  the  b o u n d a r y c o n d i t i o n y(0) = 1, b u t f r o m n u m e r i c a l p o i n t of v i e w the c o n d i t i o n i y (0) = 0 m a k e s i t 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 o b t a i n a r e a s o n a b l y solution.  accurate  H e n c e the g i v e n c o n d i t i o n i s the m o r e i n t e r e s t i n g one h e r e .  Results  a r e g i v e n i n T a b l e (8. 3),. In the n o n u n i f o r m c a s e the m e s h i s t a k e n s u c h that  3 there are  1 J m e s h p o i n t s , e q u a l l y s p a c e d , to the l e f t of x = 0. 1 and  p o i n t s , a l s o u n i f o r m l y s p a c e d , to the r i g h t of x = 0. 1. m e n t t h e r e i s o n l y one 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 a c c o u n t the n o n u n i f o r m i t y  of the m e s h .  The  J mesh-  Hence for each e x p e r i that n e e d s to take i n t o  a c t u a l finite d i f f e r e n c e equations  i n v o l v e d a r e 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 (4. 3), (with a^(x) = 0), a n d the e q u a t i o n g i v e n i n E x a m p l e (4.4). TABLE Uniform  I  8.3 Mesh  #  r  1  0 2 0  -2  1 1 1 D  2  1.000°  2  0 4 0  -4  1 1 3 C  4  . 987°  s m  c  o  r  s m  c  o  J=8  Nonuniform  #  r  3  0 2 0  -2  1 1 1 A  2  .110°  4  0 4 0  -4  1 1 3 C  3*  . 934"  s m  c  o  r  s m  c  J=64  . 954°  .486°  .923"  . 883°  •909"  1  1  _i . 782  Mesh  J=8  o  J=32  J=16  J=32  J=16 . 1  127"  .718"  2  1  .508"  J=64 2  _2 . 223  .203"  2  -3 . 725  4 i f the m e s h i s l o c a l l y  3  uniform  Example  8.4.  L e t the p r o b l e m be the s a m e a s 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 e q u a t i o n 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 with r = s = m = l .  Let  =  j -  x  T  n  m e s h i s taken u n i f o r m .  e  A  i  one  parameter  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  given by u  This approximation  0  "  a  U  l  ( " ) 2 (a - 2 ) h +  a  1  u  = 0 .  i s consistent 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 -  t h e n the o r d e r  b e c o m e s two a n d t h e f i n i t e d i f f e r e n c e  e q u a t i o n i s i n that c a s e  identical  to the one o b t a i n e d w i t h t h e p r o c e d u r e  S e c t i o n 5, l e t t i n g r = 0,  s = 2 a n d m = 0.  It w i l l be s h o w n i n S e c t i o n 6 o f C h a p t e r finite difference operator L ^ has an eigenvalue m a k e s the a p p r o x i m a t i o n  II that i f a > 2 t h e n the \ •— + oo a s h —*• 0.  u s e l e s s f o r time dependent p r o b l e m s .  g e n e r a l t h e o r y o f K r e i s s (1972) i t f o l l o w s t h a t t h e s a m e w h e n u s e d 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 , instability.  This i s supported  J =4 .784°  J =8 :  .248°  J = 16 .100°  .453"  1  F r o m the  approximation  o b t a i n e d w i t h a = 4.  8.4  J = 32  This  d o e s n o t l e a d to a n  b y the f o l l o w i n g d a t a ,  TABLE  of  J = 64 .216"  1  a 1  -54Chapter II  INITIAL V A L U E  1.  PROBLEMS  INTRODUCTION  This  chapter is c o n c e r n e d with the construction of finite difference  approximations to i n i t i a l value p r o b l e m s . procedure of the previous chapter w i l l be first order ordinary  i n detail,  extended to systems of l i n e a r conditions.  c l a s s of difference approximations appear  (Example (3.4)), i n view of its application to the equations,  " s t i f f " problems,  differential equations  P a r t i c u l a r attention i s paid to methods for  because of their applicability to the n u m e r i c a l solution  equations.  In Section  5 it i s shown how  methods for boundary value p r o b l e m s and combined to y i e l d quite general linear parabolic  stability  stability analysis of finite difference  approximations to i n i t i a l value p r o b l e m s i n o r d i n a r y s u m m a r i z e d i n Section 4.  the difference  i n i t i a l value p r o b l e m s may  equations i n one  space v a r i a b l e . 6.  The  numerical  given.  n u m e r i c a l experiments are  Finally, given.  stability  P a r t i c u l a r attention is paid  to the effect of finite difference approximations to the boundary examples are  be  difference schemes for the solution of  of such schemes i s investigated i n Section  Several  discussed  (Section 6 ) .  Well-known concepts i n the  of parabolic  construction  In p a r t i c u l a r , a constant coefficient p r o b l e m is  analysis of parabolic  are  2 the  differential equations subject to i n i t i a l  E x a m p l e s of this v e r y general i n Section 3.  In Section  conditions.  i n Section 7 the results of some  -552.  INITIAL V A L U E  PROBLEMS INORDINARY DIFFERENTIAL  I n t h i s s e c t i o n the 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 to the f i r s t o r d e r  (2.1)  linear  f(t) =  (f^t),  matrices.  H e r e w(t) = ^ w ^ ( t ) , j_ ))  f ( t ) , "",  f  0 < t < T ,  ,Wj _ ^ ( t ) ^ T  2  A(t) and  B(t) a r e  procedure,  c o n t i n u o u s (J-1) together  differentiable as necessary.  that i s a n e x t e n s i o n  boundary conditions.  order  However,  initial  that i s a s  The construction  of the p r o c e d u r e d i s c u s s e d  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  X (J-1)  w i t h the  (0 < t < T),  a d m i t s a u n i q u e s o l u t i o n w(t),  many times continuously  general  w (t),  I t i s a s s u m e d that the a b o v e s y s t e m ,  condition w(0) = w°,  chapter,  and  T  l  2  approximations  system  L w(t) = A ( t ) w ' ( t ) - B ( t ) w ( t ) = f(t),  w i l l be d i s c u s s e d .  EQUATIONS  systems,  i n the  s u b j e c t to  previous more  t h e s e w i l l not b e c o n s i d e r e d  here.  D e f i n e a m e s h S^ b y  {t. : 0 = t  S = h  The so  meshspacings are °—  (At)  c  t h a t (At)  for  n  0  <t' < r  <t  N  = T }.  = t -t ., ( l < n < N ) . n. n n-1 — —  = At = T / N for u n i f o r m meshes. '  Let  A t = max(At) , n n  A s s u m e that m i n ( A t ) > ^ A t , n—K  K ;> 1 .  some constant  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 h a v e  form 0  (2.2) v  u,  L . u = V h Li  D  n  '  l = -P  n  .u n, I  n+£  = 7, E ,f(r ), , , l < n < N , U n, SL n , I" ' — — 1=1 yb  the  •56  with i n i t i a l condition  =w°.  , n n (u^, u^,  H e r e u'n  each n the points t,  n uj  ^,  i s the approximation to ( t ) w  n  ., (1 < i < u ), a r e distinct points,  lie i n the i n t e r v a l  n-P  The  n  coefficients D  normally  . and E „ are (J-l) n, I n, I ' x  n X (J-l) matrices.  that w i l l  For  T h e s e coefficients may  be determined  by r e q u i r i n g  that  for a l l p(t)e polynomials  0  p.  !=-P  1=1  .P+P-"  ^-  Here  ,J-1 i s the space  of degree at most d.  ( F o r notational convenience  s c r i p t n w i l l frequently be omitted.) matrix  (2.3)  of a l l ( J - l ) - v e c t o r valued  T h i s requirement  the sub-  leads to the  equation  I  0  0  0  I  0  0  0  I  0  0  0  0  0  0  T 0,1 T -L 1,1  •L  -L  0, \s.  •L. 1,  T P, 1 T -L P+l, 1  [L  T P, [i. T •L P+l, u. -L  -L  -L P+H-l, 1  •L  D  D  T -P T -P l  T D  E  —  oT  i,i  =  L  0  w  i  <^  )  A  <V  - " a )B(^) , i  i  0 0  l  0  T  _  where I and 0 a r e the (J-l) X ( J - l ) identity and z e r o m a t r i x  L  0  +  E P+U--1, j x J _ M-  T L. . i s the transpose of L. ., where  0  respectively.  -57with  = ^jf^  L,^(t)  linearly  •  independent  w l ( t  and  The and  n-P+i  oJ (t _ 1  n  For  example,  these  the  polynomials  The  coefficient  certain  p +  )  scalar  chosen  U>  =  ^)  such  = 0 ,  in  X  normalization,  l<  may  and  provided  E^  (0 < i < P + u . - 1 ) ,  are  P ,  0 < 4 < I  P+l < i < P + L x - l ,  be  S e c t i o n ( I . 2).  matrices  c/(t),  that  °<i>  6  polynomials  w (x)  polynomials  in a  (Equations  in  that  defined  (2.3)  the  are  similar  fashion  (1.2.1,3) a n d uniquely  as  (1.2.14).)  defined  up  to  a  matrix  Li.  P + l , 1)  'P+H-1,1  C EE  L  has  rank  provided  equal that  to  the  P+l,  Li,  ' P + f J L - 1 , (J.  (Ji  (u-1)(J-1).  This  is  the  case  o« <y V  Here  rank £  is  (Equation  (JJL-1) ( J - 1 ) . the  L  p+1  This  in  analogue  (1.2.12)),  all  small  enough  At,  matrix  L co ^- (C )A(^)  Lgo/^a^Ac^)  has  for  and  of  A(  turn  the  can  be  is  P+  1  0  ''' true  normalizing expressed  if  L  o"  p + f X _ 1  S + 0 for  factor as  1  E  in  (y (y A  all small the  enough  previous  At.  chapter,  -58-  1  e  (2.4)  =  (-D  L Q C / -  1  "  1  ^ )  •  •  ^ t f ^ i i j )  •  f  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 w h a t m a n n e r the d i f f e r e n c e e q u a t i o n (2.2) i s n o r m a l i z e d , i s done i n s u c h a w a y  i t may  be a s s u m e d ,  that this  that  c  <  I  E  E  JI  <  >  c  i=0 for  a l l small At.  on A t .  The  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 t h a t do not  p r o o f s of the f o l l o w i n g s t a t e m e n t s  t h o s e of T h e o r e m (I. 2.1), T h e o r e m ( I. 2 . 4), a n d  depend  a r e then v e r y m u c h T h e o r e m ( I . 3.1) a n d  like  will  be  omitted.  T h e o r e m 2.1.  A s s u m e that P > 1 a n d  t h a t £ + 0.  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 P + p-1.  In p a r t i c u l a r ,  ~p the c r i t i c a l p o i n t s o f co  T h e n the o r d e r  (2.2) i s a t l e a s t e q u a l to  i f p. = 1 o r i f P = 1 t h e n 8 ^ 0 . +  ^(t) =  Hill | | (t-y.)  1=1  0_ | (t-t  i =-  of  If the p o i n t s  are 1  .) t h e n the o r d e r i s at n + l  p  l e a s t p+u..  ( T h u s t h e r e i s a ( u . - l ) - p a r a m e t e r f a m i l y of s u c h p o i n t s w i t h  parameters  \^ , (1 < i < p.-l). )  3.  E X A M P L E S OF In  FINITE D I F F E R E N C E  this s e c t i o n a n u m b e r of s p e c i f i c  approximations  to e q u a t i o n (2.1) w i l l be  APPROXIMATIONS 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  discussed.  b e e n 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 Example  3.1.  W i t h p = 1, u, = 1 a n d  = t  S o m e of t h e s e dimensional  i = i(t + t  ,) one  have  case. obtains  the  -59Box  s c h e m e of K e l l e r  (1969).  This difference approximation  may  be  w r i t t e n as  where  D  and  D  T h a t the o r d e r  - l  0  A  =  = AT  l  A  A  ^  ^  ~ ^ n - ^ -  "  •  o f 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 e m (2. 1).  N e x t i f p = 1 a n d i i = 2, w i t h I, = t . a n d L. = t t h e n the m a t r i x f ' =1 _ l '2 n e q u a t i o n (2. 3) b e c o m e s n  T  (3.1)  I  0  0  I  0  T 0,1 T •L 1,1 T •L 2,1  0  - l  D  -L 0, 2 T  L  D  -L h,2 T  E  0  T  l  E  T " 2, 2 _  0  T 0  T  0  T  L  L 2E  where  L. . = L c o ( t , ) A ( t ,) - oi (t , ) B ( t , ) , i,l 0 n-r n-r n-r n-r  and  L. = L o ) ( t ) A ( t ) - (/(t ) B ( t ), l, 2 0 n *n n n  1  v  v  x  0 < i< 2 — —  X  x  n  W i t h o> (t) = " ( t - t ) / A t , (oV) = ( t - t _ ! ) / 0  0 < i < 2 — —  X  n  n  A  t  a  n  d  n  ^ (t) = (t-t _ )(t-t )/A 2  n  1  n  i t i s f o u n d that  A-t^W  " <W  L  0,l  =  L  l,l  = AT^n-l) •  B  •  L  L  Mt )  0, 2 = A l  l,2  =  n  AT ( n> " A  fc  B  ^  •  2 t  ,  •60-  L  From  AT  2,1  A  Li  (t -i) n  A(t ) n v  (3.1) i t f o l l o w s t h a t E , A ( t  ,) = E _ A ( t ).' T a k e E , =1. T h e n E. i s 2 n 2 1 e q u a t i o n E , A ( t ,) = A ( t ). C l e a r l y , i t i s m o r e 1 n-1 n n-1  1  the  AT  2,2  s o l u t i o n of t h e m a t r i x  v  1  d i f f i c u l t n o w to g i v e e x p l i c i t r e p r e s e n t a t i o n s was  the c a s e f o r s c a l a r equations  a l w a y s obtain the c o e f f i c i e n t s  of the c o e f f i c i e n t s E^ than  i n C h a p t e r I.  Of course,  one c a n  numerically.  I n t h e s p e c i a l c a s e w h e r e A ( t ^) a n d A ( t ) c o m m u t e , n  n  i ti s possible  to s t a t e e x p l i c i t l y w h a t the c o e f f i c i e n t s o f t h e a b o v e f i n i t e d i f f e r e n c e equation  are.  In particular,  2 I i t i s f o u n d that E ^  *  h  D  T h u s the d i f f e r e n c e equation , n n-1. ) A (u -u At  which i s recognized consistency  Example  i f A i s a constant -i  =  AT  t h e n w i t h E_> =  fB(t ,) a n d D = - f - A - | B ( t ). * n-1 0 At ^ n n  A  m a y be w r i t t e n as  ^ ( B ^ u ^ B ^ ^ u  1  1  -  1  )  = |(f(t ) n  m^j)  +  as the u s u a l t r a p e z o i d a l m e t h o d .  of t h i s e q u a t i o n i s a l s o e q u a l  3.2.  matrix,  W i t h p = 2, |_L = 1 a n d  T h e o r d e r of  to t w o .  = t  ^ + € A t the d i f f e r e n c e approxi-  m a t i o n to (2.1) b e c o m e s  |(2M)u -2(e-Du - + l ( 2 e - 3 ) u n  A  ^  At  n  1  n  2  B(^ ) 1  n-1  ie(e-i)u-e(e-2)u- + n  x  + l(£-l)(£-2)  ( T h e m e s h i s u n i f o r m i n t h i s e x a m p l e . ) T h e o r d e r of c o n s i s t e n c y o f t h i s approximation order  i s equal  to t w o .  of c o n s i s t e n c y i n c r e a s e s  From  T h e o r e m (2.1) i t f o l l o w s t h a t t h e  to t h r e e i f o n e t a k e s  £ ^- \^3 =  ±  ~ 1-58.  •61If £ = 2 t h e n the a b o v e f o r m u l a  v  reduces  to  3 n _ n-1 , i n-2 — u - 2u + 2  B(t  U  n' A t  )u  = f(t ) , n  n  n  w h i c h i s one o f G e a r ' 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 . p.  (See G e a r  (1971),  217.)  Example  3. 3.  Consider  the c o n s t a n t c o e f f i c i e n t  Aw  (3. 2)  (t) - B w ( t ) = f(t) ,  w(0)  (3. 2a)  where A  = w  i s nonsingular and where A  example when A  = I.  problem  0 < t < T  ,  ,  and B  commute.  L e t p. = 2 a n d p = 1.  This occurs for  The construction procedure  of  S e c t i o n 2 a g a i n l e a d s to t h e m a t r i x e q u a t i o n (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 t h e p o i n t s  Let E B  commute,  E  = -L co (£ ) A 2  1  0  i t i s found  = L co (q)A Z  2  2  0  + ^{t,^  a n d t,^.  B.  Then,  u s i n g the f a c t t h a t A  that  - co (?, )B,. 2  1  D _ j _ = E {L co°(c: )A - u°(t )B}+ E {L co°(t: )A-co (C )B} , 0  1  and  D  Q  -  0  i  E ^ L Q W ^ C ^ A  >1  - W ^ ^ B }  2  +  0  E.,{ L.  2  Qbi  1  [t )A-( \t , )B}, >2  w h e r e LQto(t) = co'(t).  In p a r t i c u l a r  i f one u s e s the G a u s s i a n  2  points  li  2  and  -62-  = (t  n - 1  + t ) / 2 - AtN/3/6 n  and  C = 2  (^.i+tj/  2  + At^3/6  t h e n 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  {^A  (3. 3)  2  -±AB  - §B }u "  = {|A  +  2  n  ^I-B}  1  + {^A -lAB 2  f( ) 4- { ± A Cl  |  +  - ^ - B  2 B  }  } u  n  f( ) &2  T h e o r d e r of c o n s i s t e n c y o f t h i s e q u a t i o n i s e q u a l to f o u r . f(t) = 0 t h e 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  (3. 3) c a n a l s o be  f r o m a particular Pade approximation  to t h e s o l u t i o n w(t) =  (3.2),  o r i g i n a l l y due to P a d e  (3.2a).  These approximations,  been extensively studied i n connection with their (See  B  If obtained w^  of  (1892) h a v e  stability properties.  e. g. , V a r g a (1961).)  Example  3.4.  The preceding  finite difference approximation  example suggests  that the g e n e r a l  p-step  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 h a s the  form  i.=-p  where  v  =  1=1  0  v=0  a n d (3^ a r e c o n s t a n t s that d e p e n d o n the c h o i c e of the p o i n t s t,^,  (1 < SL < p.). JL  T o v e r i f y t h i s a n d to s h o w h o w found, (3. 5)  c o n s i d e r the s c a l a r  the c o e f f i c i e n t s a  equation  L w ( t ) = w'(t) - X w(t) = f(t) ,  v  JL  a n d |3  may  be  -63The  finite difference approximation  to (3.5) i s g i v e n b y  n+4  (3.6)  l=-p  1=1  where  (3.7)  e.  i-il  P+ // 1  T  |JL+4+l  X  (£ )  Lw  and  6^ = J  (3.8)  L  H + i w  (C )e k  ,  k  -P< k< 0  k=l The n o r m a l i z i n g factor  S i s g i v e n b y (2.4).  Considering  (3.5),  (3.7) a n d  (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 h a v e the f o r m u.-l  [X  v=0  for c e r t a i n constants In order i n equation  a  v=0  4  and 6  4  to s h o w that t h e s e c o n s t a n t s  (3.4),  note that b y d e f i n i t i o n  a r e t h e s a m e a s the c o n s t a n t s  •64-  Y  e La |  k )  ) = 0 ,  (^  p+l < k < p  + HL-1  4=1 i.e., .-1 = o ,  4=1 for  v=0  a l l X and f o r p + l < k < p  + [ i - l  So q^(X.)= 0 a n d t h e r e f o r e ° x ( C ) = 0 f o r a l l m a t r i c e s  C.  k  In  p a r t i c u l a r w i t h C = A B one obtains X  H  u-1 0 .  L Q C O ^ ^ J I - C O ^ C ^ A ^ B  4=1  v=0  V  U p o n m u l t i p l i c a t i o n b y A ^ a n d u s i n g the a s s u m p t i o n that A a n d B one  commute  obtains p.  u-1 L co (^)A-oo (^)B k  k  0  Q  4=1  ,  p + l < k < p +[i-l .  v=0  Thus i f one defines .-1 v=0  t h e n t h e l a s t u.-2 r e l a t i o n s i n t h e m a t r i x Further,  the f i r s t  equation  (2.3) a r e s a t i s f i e d .  p + l r e l a t i o n s o f (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 b y  letting P-  D -  S  I  E  L cop ^(c: )A(c: )+  v  0  V  v  p+i w  (c )B(c: ) V  v  T h e r e f o r e i t has  b e e n shown,  that i n c a s e  p r o b l e m (3.2), w i t h c o m m u t i n g A coefficient matrices  4.  THE  and  S T A B I L I T Y OF  INITIAL V A L U E  and  B,  of the c o n s t a n t c o e f f i c i e n t one  c a n e x p l i c i t l y g i v e the  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 ) .  FINITE D I F F E R E N C E APPROXIMATIONS  TO  PROBLEMS  In this section some w e l l - k n o w n concepts of finite 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 the s t a b i l i t y  analysis  to s y s t e m s o f 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  summarized.  d i f f e r e n c e a p p r o x i m a t i o n one  o n l y has  To  d e f i n e s t a b i l i t y of a f i n i t e  to c o n s i d e r the f o r m  this  difference  i  e q u a t i o n t a k e s w h e n a p p l i e d to the e q u a t i o n w (t) = 0 . scalar function.  T h u s f o r the t y p e of a p p r o x i m a t i o n s  H e r e w(t) i s a discussed i n Section  2 the r e l e v a n t d i f f e r e n c e e q u a t i o n i s 0 Y,  (4.1)  i =  &i u  n  +  i  = 0 ,  p < n < N ,  -p  where  5  i  (-1Y =  v^v v  p+1  v •• V^X  w i t h L Q C o ( t ) = to'(t).  T h e m e s h i s t a k e n to be u n i f o r m i n t h i s ( At)  n  = A t = 1/N,  (1 < n < N).  The  11  satisfies  so t h a t  a p p r o x i m a t i o n i s s a i d to b 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 , (4.1)  section,  (0 < n <  stable  p-1), the s o l u t i o n to  -66p-1 |u | n  < const  I '  Y  p < n <  N.  1=0  This  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 q u a t i o n 0  (4.2)  <r (rj) =  At  0  £  i  \  r f  +  =  l  0  =-p  satisfy  (4. 3)  |r). | < 1 ,  1 < i <  p ,  and  (4.3a)  |n.J  T h e s e two  conditions  are  =1  implies  the  roots  i ~ ^7/  If p = 1 t h e n the •  For  = 1 .  i  5  root condition.  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 ,  r? | < 1 ,  and  the  root.  satisfy  1j_  E x a m p l e 4.1.  simple  u s u a l l y r e f e r r e d to as  practical purposes i t i s desirable i.e.,  r\^ i s a  T  n  e  2 < i < p .  c o e f f i c i e n t s i n (4.1)  c h a r a c t e r i s t i c equation has  one  are  given by  r o o t r) = 1 and  this a p p r o x i m a t i o n i s stable.  If p = 2 and  p = 1 with  t, = t 1 n-2  + £ At, -  t h e n the  difference  6  Q  =  thus  -67approximation  to w'(t)  = 0 has  T OQ  where and  the c o e f f i c i e n t s  S_, = (2£,-l)/2At.  ri  = 1 and  1  r\  second  formula  6^  are  The  order  this c a s e  form  , -r n-1 + u  -r-  +  5., u  is strongly formula  Now  % = 2,  strongly  stable.  p. = 2 w i t h C,^ = t  2  2  =  n At  ^ +  Thus  In p a r t i c u l a r ,  and  the t h i r d  (See  £]_ A t and  order  e x a m p l e (3.2).)  £  = t  2  2  ^2  +  are  6(^+e ) - 6 ^ |  2  o n l y i f £ > 1.  for which  are  o n l y i f £ > 1.  2  are  6^ = ( 4 - 4 £ ) / 2 A t  e q u a t i o n (4.2)  | r j | < 1 i f and  s t a b l e i f and  of G e a r ,  the c o e f f i c i e n t s  b  n  = 0 ,  r o o t s of the c h a r a c t e r i s t i c  let p = 2 and  6  n  S Q = ( 2 £ - 3 ) / 2 At,  g i v e n by  o b t a i n e d i f % - 1 + -|- N/3~  Now In  n-2  = (2£-3)/(2£-l).  2  the a p p r o x i m a t i o n the  u  the  - 8  2  1  -3(e +e ) + 6 e e 1  2  x  2  + 2.  and  6  The  characteristic  equation  °"Q( " )  *i § For If £  the a p p r o x i m a t i o n  2  -9(  6  = 0 has  2  .  roots n ^ = 1 ,  e + e) + 6 x  strongly  >  2  and  2  + 2  stable i t i s n e c e s s a r y  p r o v i d e d that ^ 5  and  +14  2  x  to be  that  -  Y  -3(e + e ) + 6 ^ e  2  = 2 t h e n t h i s i s the c a s e  then i t i s n e c e s s a r y  b  = -  Q  if  <  0 or  = 0 t h e n one  ^1  >  that  2 "J*  ^  |i? |< !• 2  ^2  =  14 needs — < ^  ^ < 2.  -68Another  w e l l - k n o w n c o n c e p t i s that of s t a b i l i t y  difference approximation. eigenvalues  with large  approximations stable,  accuracy. do  A t to be  this  eigenvalues  from  requirement  then a l a r g e  this  may  shortcoming.  under  w h e r e w(t) i s a  consideration  scalar  function.  For  not be  necessary  Aw'(t) = Bw(t),  where w  p r o v i d e d that A  difference approximations  ^B  investigating  analysis  a complete  to the s c a l a r  of d i f f e r e n c e  approximations the  A  and  that  difference  f o r the  are  constant  set of e i g e n v e c t o r s .  e q u a t i o n w'(t)  of  (t) - \w(t) = 0,  valid B  be  e f f e c t o f the  j_ ^w(t)= w  remains  to  for reasons  to a p p l y the f i n i t e  i s a v e c t o r and  has  those  to the p r o b l e m  The  class  has  s m a l l f o r the a p p r o x i m a t i o n  u p o n the s t a b i l i t y i t i s c u s t o m a r y  approximations  matrices,  very  finite  (2.1) the m a t r i x A  T h u s 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  not s u f f e r  system  negative r e a l part,  requires  even though  If i n the s y s t e m  r e g i o n of a  = X.w(t) h a v e  Finite the  form  0  Y  where  the c o e f f i c i e n t s  6  6^  n+l  u ' * n  L  are  given  = 0 ,  P < n < N  ,  by  (-l) "1  L  ) A.  F o r stability i t i s again n e c e s s a r y  that the r o o t s r j . o f the  U.  characteristic  0  e q u a t i o n q(rj) =  At Y  8.r]  P  = 0 satisfy  the r o o t c o n d i t i o n  i=-P From  (4.4) i t f o l l o w s  that q(r?) h a s  the  form  (4.3),  (4.3a).  -69-  '  (4.5)  q(n)  =  <r (n) + £  ( At X)  0  cr^n) ,  V  v=l where ^ ( r j ) i s given by  (4.2)  and  w h e r e cr (r?) e P  The  stability  m a t i o n of the f o r m  . p  y  region S associated with a finite difference (2.2) i s now  d e f i n e d as  t h a t r e g i o n of the  A t X - p l a n e f o r w h i c h a l l r o o t s r j ^ ( A t X ) of q(rj) = 0 s a t i s f y if a given difference approximation i f A t X eS,  t h e n the a p p r o x i m a t e s o l u t i o n u  for  any  of i n i t i a l d a t a .  choice  For  eigenvalue  the  region S include a considerable  plane.  In fact,  entire negative satisfies  X w i t h Re(X)  < < 0,  real axis.  (See  S e c t i o n 6.)  the r o o t c o n d i t i o n i f 0 e 5S..  Here  o r i n s i d e the u n i t c i r c l e .  (If the  t h e n the m e t h o d w i l l be  w h i c h the n e g a t i v e  called  r e g i o n of the  strictly  that a r e  r e a l a x i s i s not c o n t a i n e d  r e l a t e d concept, A  w'(t) = Xw(t),  to z e r o  ",i.e.,  for  as  n —••co  systems  i t i s t h e r e f o r e d e s i r a b l e that p o r t i o n of the n e g a t i v e  introduced  N o t e that a n  half  by  approximation  contained of  cr (r\) = 0 l i e on  s t a b l e at co. ) strictly i n S.  Widlund  i n the There  s t a b l e at co, (See  (1967),  Varah  z <  tr +  unit are  but  A{a)  -stable i f S contains  a > z # 0  } .  for  (1975).)  i s that of  form  = { z : TT-a < a r g  A  i n S i s of  s t r i c t l y contained  d i f f e r e n c e m e t h o d i s s a i d to be A(a)  S  the  6S i s the b o u n d a r y of S.  i . e . , a l l roots  roots are  difference approximations however,  stability.  Thus  1  w i l l tend  r e a l a x i s to be  c o u r s e that the m e t h o d i s s t a b l e at oo,  A  l ? ^ ^ 1-  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  n e c e s s a r y c o n d i t i o n f o r the n e g a t i v e  circle  11  "stiff problems  that have an stability  complex  i s a p p l i e d to the e q u a t i o n  and  approxi-  a  -70Thus i f an a p p r o x i m a t i o n contained i n S  The  a  and  i s A ( a ) - s t a b l e , t h e n the n e g a t i v e r e a l a x i s i s  h e n c e i n S.  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  difference forms"),  of the f o r m (4. 5),  h a v e b e e n e x t e n s i v e l y i n v e s t i g a t e d f o r the c a s e  by,  for example,  D a h l q u i s t (1959  The  general case  has  and  1963),  been studied by  Gear  Reimer  (1971) a n d  (1968).  The  E x a m p l e s of t h e s e i n c l u d e m e t h o d s b a s e d u p o n P a d e r a t i o n a l  derivative methods, how  these higher o r d e r  a n a l y s i s of the v e r y  stability  forms. approximation  (1961)), R u n g e K u t t a m e t h o d s a n d  ( E n r i g h t (1974).)  In this  s e c t i o n i t has  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  (1975).  motivation for  a n a l y s i s of m a n y d i f f e r e n c e m e t h o d s l e a d s to s t u d y i n g s u c h  the e x p o n e n t i a l , (see V a r g a  u. = 1  Varah  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  to  ("finite  been  second shown  stability  g e n e r a l t y p e 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.  I t 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 a b o v e , but  by  m e a n s of s o m e e x a m p l e s the e f f e c t t h a t the c h o i c e of the p o i n t s t,^ h a s the s t a b i l i t y r e g i o n S w i l l be E x a m p l e 4.2. approximation  u. = 1 w i t h 2,^ = t  ^ + £ At.  to (_ ^ ( t ) = w'(t) - Xw(t) = 0 h a s  the  w  n-1 -l  w h e r e the c o e f f i c i e n t s  * - i Therefore  illustrated.  L e t p = 1 and  ~  =  ~  and  ii  n 0  u  6  on  -  U  a  form  _ '  =  8^ a r e f o u n d  x  So  n  d  to b e  e q u a l to  ~ o = it  - ^ •  6  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  x  by  the d i f f e r e n c e  -71q(") = 6QT7 + 6_ The  = tr (rj) + A t X ^ ( n ) , w h e r e <r (r)) = 17 - 1 a n d cr^rj) = -£"+1-1.  1  Q  Q  r o o t of 0^(77) = 0 i s 77^ = (£-l)/£, so t h a t |rj^|< 1 i f a n d o n l y i f £ > \.  T h u s 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 i s s t a b l e a t 00 i f £  I t i s e a s y to  check  £.  t h a t the a p p r o x i m a t i o n i s i n f a c t A(a)-stable  m e t h o d i s s t r i c t l y s t a b l e a t 00 o n l y i f £ > \. system  the B o x s c h e m e ,  to t h e s t a n d a r d (Keller  (1969) ),  t h e t o t a l l y i m p l i c i t one s t e p a p p r o x i m a t i o n .  E x a m p l e 4.3. q(r?)  The  W h e n a p p l i e d to the g e n e r a l  (2.1) the c h o i c e s £ = 0, £ = \ a n d £ = 1 c o r r e s p o n d  e x p l i c i t one s t e p E u l e r a p p r o x i m a t i o n , and  for such  I f one t a k e s p = 2 a n d u. = 1 w i t h ^ = t  = °- (r7) + A t X o - ^ r j ) , w h e r e <r {n) = (2^-1)^ 0  r)  Z  Q  + Zi{i-2)r\  - (£ -3£+2). 2  + £ At ,  then  + 4(1-£)n + 2£-3 a n d 0^(77) =  I n E x a m p l e (4.1) i t w a s f o u n d  that 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 s t a b l e a t z e r o i f a n d o n l y i f £ > 1, i . e . , t h e r o o t s o f °"Q(")  =  0 lie  o  n  o  r  i n s i d e the u n i t c i r c l e i f a n d o n l y i f £ >.l. T h e  a p p r o x i m a t i o n i s s t a b l e a t 00 i f t h e r o o t s o f 0^(77) = 0 l i e o n o r i n s i d e the unit c i r c l e .  These roots a r e given by  Some computation example,  r e v e a l s t h a t 177.. | < 1 p r o v i d e d £ > 1 + \ \l2  the s e c o n d  « 1. 7.  For  o r d e r f o r m u l a o f G e a r i s s t r i c t l y s t a b l e a t 00, b u t  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 t o £ = 1 + - j \l3  ,  (see E x a m p l e  (3.2)),  i s not stable at 00. 5.  INITIAL  BOUNDARY VALUE  This  section i s concerned  approximations  to l i n e a r  e q u a t i o n s o f the f o r m  second  PROBLEMS. w i t h the 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 order parabolic partial  differential  -72-  (5.1)  p(x, t) y (x, t) = L(t) y(x, t) + f(x, t ) ,  0 < x <1 ,  fc  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 initial condition  (5.1a)  y(x, 0) = g(x) ,  0 < x < 1 ,  and boundary conditions  (5.1b)  a y(0,t) + b y (0,t) 0  0  x  = g (t) , Q  0<t< T ,  and  (5.1c)  a ( l , t) + b j y ^ l , t) = g (t), l Y  :  It i s a s s u m e d that 0 < p , <p(x, t) < p the sense of P e t r o v s k i i (1938), solution,  0<t< T .  < co, so that ( l . i i ) i s p a r a b o l i c i n  and that the above p r o b l e m has a unique  that i s as many t i m e s c o n t i n u o u s l y differentiable as  necessary.  Define a m e s h R ^ ( T ) by  (5. 2)  R (T) = h  {(Xj,  t ) : 0 =x < n  Q  X ; L  < - <  X  j  = 1 ; 0 =t  The m e s h spacings a r e h . = x . - x . ^, (1 < j < J ) ,  Q  < tj < ' • < t  and ( A t ) = t ~ t n  n  n  N  = T }.  y (l<n<N).  Let (At)^  h = max  j  h. a n d A t = m a x  1  J  = A t = j^r .  n  (At) .  n  So f o r u n i f o r m  meshes  h. = h = -r a n d  n  J  J  n F o r f u n c t i o n s w(x>t) d e f i n e d o n R ^ ( T ) l e t w. = ( j > t )» w  x  n  J-l n , n w, = (w,,  Further,  n  n .T w J-l' '  w_,  for matrices  the c o r r e s p o n d i n g  In Chapter  n w. h  || • || a n d  vector  separately.  first  order  (5.3)  2  Discretizing  equations.  However,  i n space  first  o r d i n a r y differential equations.  precisely  this  norms  of the type  results This  this a p p r o a c h and time  system  c a n then  discussed i n Sections  h  h  h  1 < j<  where  w(x., t) =  J  d . . ( t ) w ( x . , t) , j  +i  i = -r. J  and  m  K  with initial  (5. 3a)  11  will  of l i n e a r be  2.3  and  s y s t e m i s g i v e n by  L (t)  by  w i l l be  in a system  K ( t ) p ( X j , t)w (x.,t) = L ( t ) w ( x . , t) + K ( t ) f ( x . , t),  h  induced  the c o n s t r u c t i o n m e t h o d of  but i n s t e a d d i s c r e t i z a t i o n i n space  solved by a difference method More  2  ||w£|| = { J (w*) } .  norms.  I to p a r t i a l d i f f e r e n t i a l  done  and  || • || ^ d e n o t e the m a t r i x  p r i n c i p l e i t i s p o s s i b l e to e x t e n d  not be t a k e n h e r e ,  4.  ma axx |w^|  =  i  n,2i  (t)w(x  t) =  J  J  V  e  J'  1  .(t)w(z J >  1  t) ,  data  w ( j , 0) = g(x.) , X  0 < j < J.  J-l ,  -74For  e a c h t t h e c o e f f i c i e n t s cL(t) a n d e^(t) a r e d e f i n e d i n p r e c i s e l y t h e  s a m e m a n n e r a s the c o e f f i c i e n t s d^ a n d e^ i n C h a p t e r I. coefficients a r e given by equations L.  (1.2.10) a n d  Thus  these  (1.2.11) w i t h L ( t ) r e p l a c i n g  ( T o 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 j i n d^ .., e^ ^ e t c .  will  usually be omitted.)  The  important  a s s u m p t i o n i s m a d e that the c o l l o c a t i o n p o i n t s  c o i n c i d e w i t h the m e s h xp o i n t s x. ., i . e . J  z^ =  f°  r  some integer I .  +  i f z. i s a c o l l o c a t i o nf p o i n t t h e n  1  x  F o r example, i f m  K ( t ) w ( x t ) = e (t) w ( , t ) = w ( X j , t ) , h  and  if m  j t  1  = 1 one m a y  have  1< j< J-l,  X j  = 3 1  K (t) w(x. h  ) t )  £  =  e.(t) w ( x . . , t )  ,  +  l<j<  J-l.  i=-l This  r e s t r i c t i o n o n the c h o i c e  (5.3)  to c o n s t i t u t e a p r o p e r  The  o f the c o l l o c a t i o n p o i n t s i s n e c e s s a r y f o 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  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  m  s t  J  fore,  using  ., t) - V e.(t) f ( z . , t) , i=l  s o l u t i o n of (5.1) s u b j e c t to ( 5 . 1 a , b , c ) .  the d i f f e r e n t i a l m  (5. 3) i s d e f i n e d a s m  r(t)£ T e.(t) p(z., t) y ( z . , t ) - Y d.(t)y(x i=l i=-r  w h e r e y ( x , t) i s t h e e x a c t  equations.  There-  equation, s  T(t) = I e.(t) L y ( z . , t ) - Y i ( t ) y(x, , t ) . i=l i=-r d  i  J  Since  the c o e f f i c i e n t s d^(t) a n d e^(t) a r e d e f i n e d a s i n S e c t i o n d . 2) o n e  may  -75apply  T h e o r e m (I. 2.1) and  Theorem  5. 2.  T h e o r e m (I. 2 . 4) to a r r i v e a t the  F o r e a c h t, (0 < t < T ) ,  l e t d.(t),  following.  e^t) b e defined b y  e q u a t i o n s ( I . 2. 10) a n d ( 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 i.n(I.2.12). i n t e g e r S.^.  Assume Also  that f o r e a c h i ,  (l<i<m),  z^ = j ^ _  f°  x  +  factor E as r  some  a s s u m e that r + s > n a n d t h a t E ± 0.  Then there  exists a constant  C t h a t d o e s n o t d e p e n d o n h, s u c h  that  < C  | T ( t ) |  r h  + + " s  m  2  I n p a r t i c u l a r , i f r + s = n o r i f m = 1 t h e n E + 0.  For  e a c h t, (0 < t < T ) ,  boundary conditions  (5.1b, c) o f the  a r e a l s o a p p r o x i m a t i o n s to t h e  form  0  s  (5.3b)  there  YJ d i=0  Q )  . w(x., t) = g ( t ) Q  ,  0 < t < T ,  and 0 Y  (5.3c) i  =  d  J , i( J w  x  + i  '  f c  )  = Ei(t) ,  0 < t < T  .  " J r  T h e s e a p p r o x i m a t i o n s m a y b e d e f i n e d a s 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 conditions are  i n d e p e n d e n t o f the d i f f e r e n t i a l e q u a t i o n .  It i s e a s y to s h o w t h a t 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 a n d d j Q a r e n o n z e r o .  This  of the  (5.3),  u n k n o w n s W ( X Q , t) a n d w ( x j , t ) i n the  l e f t w i t h a s y s t e m o f J-1 o r d i n a r y  a s s u m e d t o be  system  allows  elimination  so that o n e i s  d i f f e r e n t i a l e q u a t i o n s i n the  J-1 v a r i a b l e s  w ( X j , t),  (1 <_ j <_ J - l ) .  (5.4)  This reduced  p y t ) w (t) h  w h e r e - ^ ( t ) = (w(x^, t), "", w  inhomogeneous The  w i l l be c o m p a c t l y  = L ( t ) w (t) + f (t) , h  w  (  h  x  j y t))  w r i t t e n as  0 < t < T  h  a n d f ^ ( t ) i s the a p p r o p r i a t e  term.  i n i t i a l data a r e  (5.4a)  The  system  w (0) = g h  h  = (g(x ),  g(  1  d i s c r e t i z a t i o n c a n n o w be c o m p l e t e d  x J  _ ))  T  1  b y a p p l y i n g o n e o f the m e t h o d s  for  s y s t e m s o f 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 d i s c u s s e d i n S e c t i o n s 2, 3  and  4.  I f the s p a t i a l d i s c r e t i z a t i o n a s w e l l a s the d i s c r e t i z a t i o n i n t i m e  are consistent, Richtmeyer converges and  then i t follows f r o m  the e q u i v a l e n c e t h e o r e m o f L a x , ( s e e  a n d M o r t o n (1967), p. 45), t h a t the s o l u t i o n o f ( 5 . 4 ) , (5.4a) t o the s o l u t i o n o f the c o n t i n u o u s p r o b l e m  o n l y i f the a p p r o x i m a t i o n i s s t a b l e .  c u s s e d i n the n e x t s e c t i o n ,  whereas  (5.1), (5.1a, b, c) i f  The stability problem  specific  e x a m p l e s of finite  i s disdifference  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 t e s t s a r e g i v e n i n S e c t i o n 7.  6.  THE  INITIAL  STABILITY OF  FINITE D I F F E R E N C E APPROXIMATIONS  BOUNDARY VALUE  TO  PROBLEMS  I n 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 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 mations  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 e q u a t i o n s i n o n e  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 . equation  F o r simplicity,  we  space  c o n s i d e r the d i f f u s i o n  -77-  (6.1)  y, (x, t) = L y ( x , t) = y  (x, t)  , 0 < x < 1,  0 < t < T  with initial condition  y ( x , 0) = g(x),  (6.1a)  and b o u n d a r y  0 < x < 1  conditions  (6. l b )  y ( 0 , t) = 0 ,  0 < t< T ,  y (1, t) = 0 ,  0 < t< T  and  (6.1c)  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 gives m. V  (6.2)  the e q u a t i o n s  Si  e~.  w'(z  i=l  , t) = Y  d  w(x  , t),  1 < j < J-1,  0<t<T,  i=-r. J  where  the c o e f f i c i e n t s d. . a n d e. . a r e g i v e n b y e q u a t i o n j,i  (1.2.16) a n d  j , i  (1-2.17) r e s p e c t i v e l y , a n d w h e r e  the c o l l o c a t i o n p o i n t s z. . a r e r e q u i r e d J >  to c o i n c i d e w i t h t h e m e s h p o i n t s . uniform,  1 so that h. = h = —  Further,  and ( A t )  n  = At  1  the m e s h i s a s s u m e d =  1 |^ •  to b e  The boundary condition  (6.1c) i s a p p r o x i m a t e d b y 0 Y i  As  i n the p r e v i o u s  =  b. w ( x . , t) = 0 , J +  0 < t< T .  " J r  s e c t i o n , one m a y  u s e the d i s c r e t e b o u n d a r y  conditions  -78to  e l i m i n a t e W ( X Q , t) a n d w ( x j , t )  ordinary  d i f f e r e n t i a l equations  (6. 3)  K  with initial  h  w (t) , h  h  discretizing  6.1.  exists, a  h  .  s u f f i c i e n t c o n d i t i o n s f o r the s o l u t i o n of (6. 3),  i n time  to the s o l u t i o n of (6.1),  T h e s o l u t i o n ^ ( t ) of (6. 3), w  s o l u t i o n of (6.1), with  s y s t e m of  data  (6.3a) to c o n v e r g e  norm  h  w (0) = g  Theorem  The reduced  i s then w r i t t e n a s  w^t) = L  (6.3a)  Before  (6.2).  from  (6.1a,b,c),  and this  (6.1a,b, c) w i l l be s t a t e d .  (6. 3a)  converges  discretization  r e s p e c t to p e r t u r b a t i o n s i n the i n i t i a l  constant  c^.; that,.does , no.t  v  depend  to the  i s s t a b l e i n the i^-  data,  provided  o n h, s u c h  that t h e r e  that f o r a l l s m a l l  e n o u g h 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 spatial discretization  2.  K, i s n o n s i n g u l a r . h  3.  T h e eigenvalue  problem  (6.3) i s c o n s i s t e n t w i t h  L^v^  =  ^h^h h v  a  d  m  i * -  s  (6.1),  (6.1b,c).  complete  a  s e t of  eigenvectors. 4.  The eigenvalues  X., ., (1 < i < J-1), n, I  of the a b o v e  eigenvalue  problem  s a t i s f y Re(\, .) < 0. ' h, x — 5.  Let  Proof.  be the m a t r i x  of e i g e n v e c t o r s .  L e t e (t) = y ( t ) - w ( t ) .  Therefore  h  h  h  Then K  h  Then  ||||^ j | " V | | ^  e^(t) = L e ( t ) + T ( t ) . h  h  h  —  ° 2'  t K T L, e (t) = e e (0) + /  t e  1  n  n  h  h  =V  h  tA,  6  (t-r)K~ U l  &  h  t  Vh<°» / h e  T ( r ) dr  +  V  e  (t-r)A, h  0  ,  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 o f e i g e n v a l u e s . If n  g  i s the o r d e r o f c o n s i s t e n c y of (6. 3), t h e n the e s t i m a t e  |K(t) || < c,{ || ei i(0)|| + To-max.. || ,(t) n <t< T 2  c  c  T  c  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 o b t a i n e d b y a p p l y i n g one o f the c l a s s o f m e t h o d s d i s c u s s e d i n S e c t i o n s 2, 3 a n d 4 to the s y s t e m ( 6 . 3 ) . that i f  and  From  c o m m u t e then this a p p r o x i m a t i o n 0  _ u.  l=-p £  V T  ^h  v=0  i  l  %  .A  = 0 ,  i s used at each time-step.  Further, i fp > 1  then s u f f i c i e n t l y a c c u r a t e  uj^ , ( 1 <_ n < P - l ) , m u s t be g i v e n a l s o .  1  < n < N  T h e i n i t i a l data a r e approximations  (p = l ), h a v i n g  the s a m e  of c o n s i s t e n c y as (6.4).  The difference approximation  (6.4) i s s a i d to b e s t a b l e i n the l ^-  n o r m i f t h e r e e x i s t s a c o n s t a n t K, t h a t d o e s n o t d e p e n d o n h a n d A t , such  to  These m a y be obtained f o r  e x a m p l e by i n i t i a l l y u s i n g a one s t e p m e t h o d , order  h a s the f o r m  a r e a s s u m e d to b e i n d e p e n d e n t o f n, s o t h a t t h e s a m e  difference approximation u ^ = g^.  I t -.11  A  J  The coefficients a  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  -80-  K i i < KiZ"  1  2  Kii2>  1  /  2  •  p <  n  <  N  '  1=0 for  a l l s m a l l enough h and A t .  hence f o r convergence, Theorem  6.2.  The  Sufficient conditions for stability,  a r e g i v e n i n the f o l l o w i n g  s o l u t i o n u ^ of (6.4),  theorem.  with appropriate  initial conditions,  c o n v e r g e s to the s o l u t i o n of (6. 1 ), ( 6 . 1 a , b, c ) , a n d t h i s d i f f e r e n c e is  s t a b l e i n the i ^ - n o r m  provided o n h,  that t h e r e  with respect  to p e r t u r b a t i o n i n the i n i t i a l  exist p o s i t i v e constants  c^  and  scheme data,  a n d c^, that do n o t d e p e n d  s u c h that f o r a l l s m a l l e n o u g h h the f o l l o w i n g c o n d i t i o n s a r e  satisfied. 1.  The difference scheme  2.  i s nonsingular.  (6.4) i s c o n s i s t e n t .  2a. K, a n d L, c o m m u t e , h h 3.  The  eigenvalue  problem  L*^v^  = X^K^v^  admits a complete  s e t of  eigenvectors. 4.  The eigenvalues T  a  \, . of the a b o v e e i g e n v a l u e n, I o f the f o r m  T  =  {z e C : T-a  < a r e z < ir +a ; z ± 0 \  w h i c h h a s the p r o p e r t y that i f z e T  satisfy  liein a  region  ,  , t h e n 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  problem  v  =0  -81-  1 ,  \ri (z)\< i  and 1 ^ ( 2 ) 1 <1  5.  Let  - c  , 2 < l < p .  i  be the m a t r i x of e i g e n v e c t o r s .  l hl  II^'IU < c  v  Proof.  By assumption  a l s o using 2 and 2a, 0 K  2  L  A  Kfa  t  K  h  h  L  0  p.  L  L v  X t  K  h  V  v  a  h  C  A  t  h % 4 = -p v=0  the difference  t  0  ^  2  YJ \, i' K AT  l=-p =0  v  n  '  x  i)  V  c  h i ' +  J  = °>  1 L e t c! . denote the v e c t o r (c-P . , c^ f , -h,J h,j h,j 1  P"  • |L = { E >3 =o L  m  2  1  I  c  A  '  S " J  m  I  h  A  c  h  A  h h  =  c  0  -  Hence each component cl^ j  equation  v  h  A  is the diagonal m a t r i x of eigenvalues.  (6.5)  II  h  p.  h h L v  v  =0  0  satisfies  Substitution into (6.4),  p.  4=-p  of  c^.  gives  AT  where  2  3 one c a n w r i t e u ^ =  I % 4=-p v=0  h  Then  | } 2  "',  P < n< N = T / A t . +1 ^ c5* P ) and let h,j + 1  Define the p X p m a t r i x A ( z ) by  -82-  • p-l P (z) p  "Pp-2 (z)  ( z )  ( Z )  P p  p  P (z) p  1  0  0  0  1  0  A(z) =  1  v , ^ a m z' (z) = 7. v=0 v= one step difference equation  0<^m<^p.  where p  :JJ . = A ( A t X ,  -) c ? " .  1  ,  0  T h e n (6.5) c a n be w r i t t e n as a  p < n < N ,  1 < j < J-l  The eigenvalues I j ^ A t X ^ j) of A ( A t X ^ ..) a r e p r e c i s e l y the roots of the characteristic  equation 0  n  The definition of T  i m p l i e s that i f X, . e T then AtX, . e T . a h,j a h,j a a s s u m p t i o n 4 i t now follows that  From  r  |n (Atx )| < i , ;1  hJ  and |n 1  for  all j ,  (1  m  (AtX, .)| < 1-c,  h, j  1  —  , 1 '  2 < m < p , —  < j < J-l).  It is w e l l known, (See R i c h t m e y e r and M o r t o n (1967), p . 86), that  -83this  i m p l i e s that the f a m i l y of m a t r i c e s  || ( A ( z j ) | | n  Thus under  these  z  < K  x  JA(z) : z e T  ,  p < n < N  }  satisfies  .  c o n d i t i o n s one h a s n-p+1  lcS.jlSbS.,12 = l K , i ) and  ^.jllz^lKjIlz  therefore  J=1  '  j=l  J  '  j = l m=0  J  '  m=0  J  so that  K » 2  = l|V cJJII < | V l l l | c S l l < h  h  2  = "Viz  2  2  Il^lll}*!  I I V J I . K , !  J c ™ \ \ l f  J  m=0  I I V J U I V ^ K / Z  m=0  w^wif  m=0  m=0  T h i s p r o v e s the t h e o r e m .  B e f o r e p r o c e e d i n g to a n a l y z e a n u m b e r of e x a m p l e s , s o m e are i norder.  remarks  F i r s t , i t s h o u l d be n o t e d , t h a t c o n d i t i o n s 3, 4 a n d 5 o f T h e o r e m  (6.1) a n d T h e o r e m ( 6 . 2) a r e i n g e n e r a l d i f f i c u l t t o v e r i f y a n a l y t i c a l l y . a good check on whether a p r o p o s e d  However,  d i s c r e t i z a t i o n l e a d s to a s t a b l e d i f f e r e n c e  scheme, i s obtained by n u m e r i c a l computation  of the e i g e n v a l u e s a n d t h e  -84-  condition number The  K ( V ^ ) = II V j J |  II ^ i / l l 2  ^  a  e  w  v a  ^  °^ k .  n e s  a s s u m p t i o n t h a t a l l e i g e n v a l u e s l i e i n the l e f t h a l f p l a n e i s not  strictly necessary, under  2  but r e a s o n a b l e  consideration.  In p r o b l e m s w h e r e  than a fixed finite number n e e d o n l y a d d the  i n v i e w of the  differential equation  there are  some,  but not  more  of eigenvalues w i t h p o s i t i v e r e a l p a r t ,  requirement  t h a t the  modulus  one  of t h e s e e i g e n v a l u e s  can  be u n i f o r m l y b o u n d e d . The  essential difference  s e c t i o n a n d the  b e t w e e n the  stability analysis i n this  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  methods  for  systems  o r d i n a r y differential equations  i n S e c t i o n 4 i s the f a c t t h a t t h e  system,  not h a v e  i n the  (equation (6.3)),  requirements  admit a region T  that as  does K(VJ^  a fixed  dimension.  be b o u n d e d a n d t h a t the  i n T h e o r e m (6.2).  d e f i n i t i o n of A ( a ) s t a b i l i t y . ) then it is A ( a ) - s t a b l e , one  step methods  are  e q u a t i o n q(n) = 0 h a s  converse  the  Other  strongly A(a)-stable methods  are  a n d e x a m p l e ( 1 . 4 . 2 ) ) , a n d the  second derivative methods  (1974).  These  methods  are  the f o r m u l a s  u s e d i n the  In t h e s e p a p e r s the parabolic transforms  systems  scheme,  of G e a r ,  characteristic examples (see  of  example  (3.2)  given in Enright  numerical examples  of S e c t i o n 7 .  n o t i o n of  i n t r o d u c e d b y W i d l u n d (1965 a n d 1966).  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 without boundary conditions is considered,  r e p l a c i n g the  a  A l l A(a)-stable  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 t o the p a r a b o l i c i t y of a d i f f e r e n c e  admits  strongly A(a)-stable,  because  root i n that c a s e .  results  stability region  n e e d not be t r u e .  strongly A(a)-stable, o n l y one  This  (See S e c t i o n 4 f o r  Thus if a method is  but t h e  current  A m u l t i s t e p m e t h o d that  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.  of  eigenvalue-eigenvector  (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).)  linear  with Fourier  approach i n this The major  for  section.  disadvantage  -85of t h e 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 p e r h a p s the f a c t that t h e a s s u m p t i o n  of a c o m p l e t e s e t o f e i g e n v e c t o r s t h a t t h e e f f e c t of v a r i o u s investigated numerically.  may  n o t be s a t i s f i e d .  a p p r o x i m a t i o n s to the boundary conditions  (1971a).  (1970, 1971), a n d f o r p a r a b o l i c  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  G o t t l i e b a n d G u s t a f s s o n (1975).  of V a r a h s e e  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 b o u n d a r y c o n d i t i o n s w a s done b y O s h e r c o n d i t i o n s f o r s t a b i l i t y of one s t e p m e t h o d s , a s g i v e n the  (1968).  The  s c h e m e m u s t be s t a b l e w h e n a p p l i e d t o a p u r e i n i t i a l v a l u e  amplification matrix  Q_^u^ ^ s h o u l d  Q = [QQ]  problem.  u  z with  | z | > 1, a n d i n a d d i t i o n ,  n e a r z = 1 m u s t be s a t i s f i e d .  s h o w n b y V a r a h h o w 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 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 obtained where the eigenvalue values and  problem L ^ v ^ = ^ ^ n  z w i t h R e ( z ) > 0.  (6.4).  v n  n  z with  It i s a l s o | z | > 1.  when analyzing the case  °^ t h i s s e c t i o n h a s  T h i s w i l l be c o n s i d e r e d  again  eigensystem analysis.  p u r p o s e of the next e x a m p l e i s t o i n d i c a t e w h y i t i s d e s i r a b l e  at l e a s t the e n t i r e negative  E x a m p l e 6.3. each eigenvalue  (6.6)  eigen-  i n E x a m p l e (6.3)  to r e s t r i c t t o 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 includes  some  It i s n o t c l e a r h o w t h e c o n d i t i o n o n t h e s p e c t r u m o f Q r e l a t e s  to the c u r r e n t  The  Further  °f t h e d i f f e r e n c e s c h e m e Q g j ^ =  h a v e no e i g e n v a l u e s  c o n d i t i o n on t h e s p e c t r u m of Q  sufficient  by V a r a h , a r e f i r s t that  ( T h u s t h e w o r k o f W i d l u n d c a n be u s e d t o c h e c k t h i s c o n d i t i o n . ) the  c a n be  T h e s e effects have a l s o been investigated 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, systems,  A n advantage i s  real axis.  It f o l l o w s f r o m t h e g e n e r a l of the eigenvalue  Ly  = y"  theory  of K r e i s s  problem  = Xy  that  ,  y(0) = y'(l) = 0  ,  (1972), t h a t  -86is approximated by an eigenvalue of the discrete eigenvalue problem ( 6  -  - W h  Vh  7 )  .  as h —- 0, provided that L^u-^ = K^f^ defines a consistent and stable finite difference approximation to L,y=f. values of (6.6) are given by \  (See Section (I. 7) J  = - L"  1 '  5  fc  k  - *» ^  Since the eigenf o l l o w s  are eigenvalues X.^ of (6.7) such that \^ -*• -co as h — 0.  t h  a t there  Hence the  restriction to discretizations i n time with an unbounded region T 0  Theorem (6.2) is desirable.  If T  a  a  as i n  were bounded then a restriction on the  size of the time step would have to be imposed.  F o r example, consider  the case where L, u. = -— (u. ,-2u.+u. ,,) and K. u. = u. , (1 < i < J-1). J 2 j-1 j j+r h j J - 1  The  J  h  discrete boundary conditions are  UQ  =  0 and  (Uj-Uj_^)/h  = 0.  Then the  eigenvalues of (6.7) are easily found to be given by  Thus i f the region T  a  includes only a part (-(3,0) of the negative r e a l axis Q  then At must satisfy At < -*r h ' — 4 The  2 in order to guarantee that At L . e T . h, x a to  stability of particular finite difference approximations w i l l be  investigated i n the next two examples as w e l l as i n numerical examples i n 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 eigenvalue problem (6.7).  Specifically, it i s 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 -("v"^) is bounded. K  -87-  Example 6 . 4 .  -fAlthough stability and  approximation  imply that the eigenvalues of ( 6 . 6 ) are approximated  consistency of the spatial difference  those of the discrete eigenvalue p r o b l e m ( 6 . 7 ) ,  by  it need not be true that  a l l eigenvalues of ( 6 . 7 ) converge to eigenvalues of  (6.6).  C o n s i d e r for example the case where L ^ and  are as i n the  previous example, but where the boundary condition at x = 1 i s a p p r o x i mated by  =  (2-c) h  i s consistent with y (1) = 0 ,  T h i s approximation  0  p r o v i d e d that c + 2 .  The  4  order of consistency i s equal to one.  If c =-j  (See also E x a m p l e ( 1 . 8.4').)  is equal to two.  F o r this d i s c r e t i z a t i o n the  i n ( 6 . 3 ) i s the identity and the operator L ^ can be w r i t t e n as  operator = Pjjl^,  where " -2  P  then the order of consistency  1  and  h = 2-c  1  -2 1  i  ~  Let v ^ = JL ~  v^.  T h e n the eigenvalue p r o b l e m T ^ v ^ = X ^ K ^ v ^ becomes  JL ~  ~  P £ -J-hPj^h  =  ^h h' v  w  ^  c  ^  with r e a l eigenvalues.  admits The  a complete  matrices  set of o r t h o n o r m a l  and  containing the eigenvectors  ~ v ^ ^ and v ^ ^ r e s p e c t i v e l y by column are then r e l a t e d by JL  ^ )2I U h"2 2  K ( V  h  )  "  V  h  2  V  h  2 -  " 'h"2  11  "h  "2  eigenvectors  II 7 <  -1 ~ = P^h"  Hence  max{|c-2|, - }, c-2 J  T  T  -88p r o v i d e d c + 2.  T h i s shows that the l a s t c o n d i t i o n of T h e o r e m (6.1) and  T h e o r e m (6.2) i s s a t i s f i e d for this d i s c r e t i z a t i o n i n space when c + 2. Next c o n s i d e r the e i g e n v a l u e s .  A p p l i c a t i o n of the G e r s g o r i n t h e o r e m ,  ( L a n c a s t e r (1969), p . 226), shows that a l l eigenvalues a r e negative i f and only i f c < 2.  If c = 2 then there i s an eigenvalue at z e r o and i f c > 0  then there i s a s t r i c t l y p o s i t i v e e i g e n v a l u e .  So a s s u m e that c > 2, X.^ > 0  and let v ^ be the e i g e n v e c t o r a s s o c i a t e d w i t h X.^. v.,,  - (2 + h X.)v. + v . , = 0, 2  j+1  J  J-1  The s o l u t i o n to these difference equations v. = c . z i + c_z,^ , j l h 2 h '  Thus v ^ s a t i s f i e s  1 < j < J-1. -  -  is 0 < j < J. — —  w h e r e z. and 1/ z, are the roots of h h (6.8)  z  - (2 + h X ) z  2  +1=0  2  But v ^ must a l s o s a t i s f y the d i s c r e t e boundary c o n d i t i o n s . the c o n d i t i o n V Q  In p a r t i c u l a r ,  = 0 i m p l i e s that v~ = C j ( z ^ + z ^ ) , (0 < j < J ) , n d the 1  a  boundary c o n d i t i o n at x = 1 leads to the equation  (6.9)  ( h~ h ) " z  z  J  c  (  z  h"  1  _  z  h  J  +  1  )  +  (c-l)(z^" -z^ 2  J + 2  )  = 0.  C o n s i d e r now the equation obtained f r o m (6.9) by o m i t t i n g a l l negative powers, v i z .  z. - c z , + c-1 = 0. h h 2  T h i s i s the c h a r a c t e r i s t i c equation of the boundary c o n d i t i o n at x = 1. roots are z^ = 1 and z  2  = c-1.  Its  Since by a s s u m p t i o n c-1 > 1 it follows that  there is a root z ^ of (6.9) s u c h that z ^ — c - i as h -*• 0.  The  -89c or responding  eigenvalue  is obtained  from  (6.8)  and  asymptotically  given  by  h z  h  Thus  X^  —>  equivalent  +  h —» 0.  to the t e c h n i q u e  cretization of the  oo as  (see  h (c-l)  2  Varah  2  h  This  eigenvalue  of V a r a h w h e n a p p l i e d to the  (1970), p.  33, a n d  discrete boundary conditions  of the a m p l i f i c a t i o n m a t r i x , eigenvalues  of K ^ L ^ .  computation is essentially  (1971).)  There  i s r e l a t e d d i r e c t l y to the  whereas  If the  Varah  current  here  these  effects are  dis-  the  effect  eigenvalues  r e l a t e d to  d i s c r e t i z a t i o n i s c o m p l e t e d by  applying  the a  i  multistep  m e t h o d to the  system  K^u^  = J-^^ >  n  m e t h o d i s s t r o n g l y A(a)-stable,  but  in its stability  resulting scheme  c < 2.  The  rule.  (See  simplest  c h e c k e d by  example  half - plane.  applying  f o r w h i c h the  e n o u g h h,  co h a v e t h i s Example  6.5. (6.7)  2.  this  That this  i f this  multistep  p o s i t i v e r e a l a x i s i s not  scheme  m e t h o d i s the stability  is unstable  g ^ = v^,  trapezoidal  region  v^  There  d i f f e r e n c e s c h e m e w i l l be a l l m e t h o d s that a r e  only i f  i s the  for c > 2 can  where  positive eigenvalue  In f a c t ,  contained  i s stable i f and  m e t h o d the  i n i t i a l data  complete  even if c >  d  of s u c h a m u l t i s t e p For  i t to the  a s s o c i a t e d w i t h the  in time  version  t h e n the  E x a m p l e (1.4. 1).)  entire negative  vector  region,  the  a  i s the  be  eigen-  are discretizations stable f o r  small  s t r i c t l y s t a b l e at  property. Again consider  the  eigenvalue  problem  with  -4-  L, v. = n J h  2  (v. , - 2v. 3  3  +  v..,) 3  (6.6)  and  its discrete  -90-  and K, v . = -r— v . - + -pr v . + -pr v . , , h j 12 j-1 12 j 12 j+1  Thus this eigenvalue p r o b l e m a r i s e s i n the s t a b i l i t y a n a l y s i s when u s i n g the fourth o r d e r s p a t i a l difference a p p r o x i m a t i o n d i s c u s s e d i n E x a m p l e (1.4.4).  L e t the a p p r o x i m a t i o n to the boundary c o n d i t i o n y ' ( l ) = 0 have the  form  2  (6.10) 1  =  _  i  b  =  u J  +  i  0  J  r  It i s now difficult to check ^ ( V ^ ) ,  so only the eigenvalues w i l l be c o n s i d e r e d .  P r o c e e d i n g as i n the p r e v i o u s e x a m p l e , i t i s found that the  eigenvalues  a r e g i v e n by  x. =  (6.11)  for  2  ll  1 1  1 2  2  ( Z  -  1 ) 2  (z* +  10Z  +  1)  each z w h i c h i s a root of 0  2  (6.12)  i=Note that i f z i s a root of  " z"  then so i s \ .  (6.12)  j( ) = . Z z  1  Now  equation  (6.12)  =  ) =  0.  P j  n o m i a l of the d i s c r e t e boundary c o n d i t i o n  c  ( J + 1 )  _  r  b. z  (6.10)  The c h a r a c t e r i s t i c p o l y i s g i v e n by  J  J  has a root z s a t i s f y i n g | z | > i + e , J l a r g e , i f and  only i f the c h a r a c t e r i s t i c equation C j ( z ) = 0 has a root i f z i s a root of  (6.12)  with | z |  = 1 ,  |z[ > 1 .  Further,  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 a n d n e g a t i v e .  Thus i t i s sufficient  t o c h e c k t h a t n o n e of t h e r o o t s o f 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 of 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 e x c e e d s one i n a b s o l u t e v a l u e .  As  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  h  (  2  J  U  "  2  u  J-l  2 J-2  +  u  )  =  approximation  °  s a t i s f i e s t h e e i g e n v a l u e c o n d i t i o n b e c a u s e t h e r o o t s of i t s c h a r a c t e r i s t i c equation a r e  = 1 and z^  the t h i r d o r d e r  approximation  6S<  l  l  u  J-  1  8  J-l.  u  has  r o o t s z ^ = 1, z ^ =  The  characteristic  T2h  has with  ( 2 5 u  roots z  y  J  9 u  "  4  8  = 1, z  2  E x a m p l e 7.1.  "  2 u  ^  z  J - l  u  +  3 6 u  J-3  3  }  °  =  22—  =  J-2  -  l 6 u  J-3  +  3 u  w r f c  h  \z^ \= | z | = 0. 426,  approximation  J-4  )  =  also  condition.  EXAMPLES  T h e p u r p o s e of t h i s e x a m p l e i s t o i n d i c a t e t h e r e l a t i v e  c o n s i d e r the s i m p l e  ^xx  0  T h e r e f o r e , these difference equations  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 .  =  equation of  = 0.381, z^ = 0.269 + 0.492 i a n d z^ = 0.269 - 0.492  \z^\ = I z^ I = 0.561.  NUMERICAL  J-2  e q u a t i o n of t h e f o u r t h o r d e r  s a t i s f y the e i g e n v a l u e 7.  +  S i m i l a r l y , the c h a r a c t e r i s t i c  +  ( 1  F o r this  problem  +  """  2)y  '  0  <  x  <  1  >  ° < t < 3 ,  purpose  -92with initial  condition  y(x,0) = s i n ir x  and boundary  ,  0 < x <  conditions  y(0,t) = y ( l , t ) = 0 ,  The  solution i s y(x,t) = e ^ s i n f i r x ) .  reaches  i s y{j,3)  At = j | .  1 ,  = e  3  The  0 < t < 3 .  m a x i m u m value this  1 U n i f o r m m e s h e s a r e u s e d w i t h h = — and  = 20.  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 h a v e the 0  l=-p  solution  form  u 1-^=0  n+i u. h  v  =0 n  where  V j .  =.2  d  i=-r and  i V i J  m K , u ! = ). e. u ( z . . ) . h J i=l J' n  1  The  1  c o l l o c a t i o n p o i n t s z. . c o i n c i d e w i t h t h e m e s h p o i n t s .  z. . = x. a n d i f m J,i J problem Collatz.  = 3 t h e n z. . = x. _,. . (1 < i < 3). 3,1 J-2+i -  the c h o i c e m  If m  = 1 then  F o r the g i v e n  = 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 " o f  (See e q u a t i o n (1. 1. 7) a n d E x a m p l e (1.4.4).)  R e s u l t s a p p e a r i n T a b l e (7.1). in experiment  5 i s not c o m p a c t .  The  spatial difference  T h u s n e a r the b o u n d a r i e s  approximation other  difference  -93TABLE  #  r s m o  P c o  1  1 i  1 C 2  1 2  J = 5  7. 1  J = 20  J = 10  +2 .392 ^ .320 + 2  . 311  .772+  1  .279  +1  .585+  1  .149  +1  .560+  1  . 132+  J=40 .170+  1  . 127 +  1  2  1 1 1 2  2 G 2  .599  +3  .755  + 2  .475+ .382+  2  2  4  5  1 1 3 4  2 2 1 4  1 P 4  4 G 4  4 G 4  .516+  1  .795+  1  .306+  1  .669+  1  .212+  1  .609  .371°  . 189"  .376°  .241"  1  .242"  1  +  1 1 3 4  .108+2  .343+ +2 .325 ^ 2  3  .185+2  +2 .523 ^ .106+2  .861+2  +1  .161°  .299°  .293"  1  .223"  1  .229"  1  1  1  1  .167"  1  . 169"  .360°  30  .321°  60  .311°  120  .459+2  6  .888+  15  1  30  120  .168+  .701°  240  .380"  2  .521" -4 .399 -4 .859  1  1  2  . 138"  2  . 150"  2  .224"  .149" , -2 .147 2  .986"  1  . 154" -4 .959 -4 .932  15 30  3  3  60 120 240 6  . 109"  1  15  . 112"2  30  2  .299"  3  60  - 2  .245"  3  120  .245"  3  240  2  .348" .353  30  .150°  -1  .135 .4H"  15  6  . 108" 2  60  6  .150° 1  1  15  1  .151°  .245"  .517°  + 1  . 118"  .158°  6  1  .391+ +1 . 197 .111+  .151°  .308°  +3 . 120 +1 .212 .243+  +  N  -94-  TABLE  7.1 (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  .136"  6  .351°  .222"  1  .978  .234"  1  .143"  2  . 148" _3 .343 .631"  .148"  2  .152"  2  7  8  9  1 1 3 4  1 1 3 4  1 1 3 4  1 P 6  6 G 6  4 S 6  c  . 377°  .243"  1  .377°  .242"  1  -1  .351  .536"  .238°  . 156"  .307°  . 198"  1 1 3 4  1 P 8  1  1  .300°  . 194"  1  1  . 377°  .242"  1  -4  . 108" -4 .947 _4 .944  .H6"  .330"  2  3  .334"  3  .120"  2  .136"  2  6  2  .370" -4 .719 _4 .848  2  15 30 60 6  3  .663"  6 30  .268" -4 .824 * -4 .858  -2 .151 . 151"  30  15  3  2  15 60  .909  2  .124"  . 116"  .242"  4  4  .342"  .186°  . 377°  -3  . 151  2  1  1  -2  1  .218" 10  1  15 30 6  4  . 109" -3 . 101  3  15 30  -95formulas values  m u s t be  r = 1,  values are formulas  r = 4,  "p",  "C"  s = 1 and  (7.5).  The  the o r d e r  "G"  coefficients (1961).  as  = 1.  The  order  i s the o r d e r  discretization  of c o n s i s t e n c y "o",  at x = Xj. ^ t h e s e  of c o n s i s t e n c y of t h e s e  of the m a i n f i n i t e  difference  d i s c r e t i z a t i o n w i l l be  in time and  the  considered  i s d e f i n e d by the n u m b e r a code  "c", indicating  of  the  (See  is usually  E x a m p l e (1.4. 1) a n d  r e f e r r e d to as  a  Example  Crank-Nicolson  in time.  be  Example  Pade table  i s defined by  used.  denotes  can  "P"  m  i s the t r a p e z o i d a l r u l e .  discretization  at x = x^  = 1, w h i l e f o r the a p p r o x i m a t i o n  This approximation  methods  approximation  S t a b i l i t y of t h i s p a r t i c u l a r  t y p e of m e t h o d  and  m  i s e q u a l to f o u r ,  in Example  (3.1).)  The  s = 4 and  approximation.  steps  used.  a m e t h o d of G e a r .  found i n G e a r  (1971),  The p.  coefficients  217.  (See  a  also,  I v  of t h e s e E x a m p l e (1.4.2)  (3.2).)  stands f o r a of r a t i o n a l of t h i s  approximations  discretization  (See E x a m p l e  Finally,  Pade method,  "S"  based  u p o n a d i a g o n a l e n t r y i n the  t o the e x p o n e n t i a l f u n c t i o n .  in time  may  (3.3) f o r a f o u r t h o r d e r  denotes  m e t h o d s h a v e the f o r m  a  so-called  second  (7.1) w i t h u. = 2 a n d  be  obtained f r o m  The Varga  Pade' f o r m u l a . )  derivative  c a n be  method.  These  f o u n d i n E n r i g h t (1974).  I The  coefficients  a^  are  determined  from  c o n s i s t e n c y and  A(a)-stability  requirements. Example  7.2.  In t h i s  e m p l o y e d to s o l v e the y  t  = f  e x a m p l e u n i f o r m as w e l l as  non-uniform  meshes  equation x u ^  - x(l-x) u  x  + xu,  0 <x<  1 , t >  0  ,  are  with initial  -96-  condition  u(x,0)  and boundary  = s i n TTX  = u(l,t)  E x a m p l e (1.4.3).  If  the m e s h i s u n i f o r m t h i s  standard second order  three-point  t h e n the c o l l o c a t i o n p o i n t c i r c u m v e n t this  the  , = x.  difficulty,  K u.  =  h  where  the  0 x  =  )  the  +  J  1  w°(z )M._ +  )  h.i.+h.')  1  (x-x. ) (x-x.) = -r y~— r*- .  w  (  x  )  -  w (z 1  = —  h  t  !  in  r e d u c e s to the  the m e s h i s not u n i f o r m ,  d o e s not c o i n c i d e w i t h  j i l  precisely,  + w (z  are  defined by  +  2  1  j  j  l  )u.  +  1  )  —  + i  of  is g i v e n by  )u.  )(x-x.  h.  x..  is redefined in t e r m s  More  (x-x.  1  »  formula  operator  1  +  second case discussed  r  ~  formula.  Li  +  L a g r a n g i a n b a s i s functions w (x)  (x-x )(x-x.  2 w  z.  Lagrangian interpolation  (  .  formula. If ( h h ) j  w  = 0  d i s c r e t i z a t i o n i n s p a c e c o r r e s p o n d s to t h e  To  1 ,  conditions  u(0,t)  The  0 < x <  ,  a  n  d  t  (x)  (h -h i  =  V j  i +  1  )(2h  9h.(h  j +  h  i +  1 +  j + 1  )  T h i s y i e l d s the  h )  (2h. h.  i  +  U  j-1  + 1  expression  )(h. 2 h . )  9h. h  +  j  +  (h.  + 1  +  V  1  1  -h.)(2h. h.  9h.  +  +  1  (h. h. +  +  1  + 1  )  )  V l  i The  resulting  order. time.  spatial difference  approximation  K^u^ = L^u^ remains  T h e d i s c r e t i z a t i o n i s c o m p l e t e d b y a p p l y i n g the t r a p e z o i d a l The complete  difference  scheme is  then  second rule  in  -97, n n- i " (u, - u, ) h h At  xr K  h  0 u.  scheme  Other  m o r e than  one  +  U  , n-1. h ' }  J  i s a g e n e r a l i z a t i o n of the w e l l - k n o w n  s c h e m e to n o n - u n i f o r m unique.  , n  T  . , . sin(irx.).  =  J  T h u s this  . 1 2  =  meshes.  second  order  T h i s g e n e r a l i z a t i o n i s by  formulas  for non-uniform  b,  c,  d) s h o w the  approximate  1  f u n c t i o n a l o n g w i t h the approximate  and  At  approximate  = —  .  of a  computation  F i g u r e s (7. 2a,  b,  c,  d,  and  d i s t r i b u t e d as  indicated.  on  e).  The  i s a b l e to f o l l o w the  computation  on the u n i f o r m  In t h i s  use  derived also.  s o l u t i o n i n the  s o l u t i o n at v a r i o u s t i m e s .  a non-uniform  Now  there are  timestep  computation  be  P l o t t e d a r e the  solution i s interpolated with a cubic  Results  means  case  1  w i t h h. = 47J  mesh,  no  m e s h e s that  e v a l u a t i o n of e a c h c o e f f i c i e n t f u n c t i o n m a y  F i g u r e s (7. l a ,  of a u n i f o r m  Crank-Nicolson  i s as  The  spline. mesh are  o n l y 31  before.  shock profile  initial  given in  meshpoints Clearly,  this  for a longer time  than  the  mesh.  Example  7.3.  e x a m p l e the  schemes  f o r s o l v i n g the d i f f u s i o n  stability  of a n u m b e r  of f i n i t e  difference  equation  y (x,t) = Ly(x,t) = y ( x , t ) , n  t  w i t h y ( x , 0) = g(x) a n d  y(0,t) = y ( l , t )  T h e o r e m ( 6 . 2) i t f o l l o w s t h a t f o r t h i s behavior  of the e i g e n v a l u e s &  number  KCV^)  must  be  = 0,  is tested numerically.  purpose  one  has  to i n v e s t i g a t e  From the  X., of L. v. = l K , v , , A l s o , the c o n d i t i o n h h h h h h ' bounded. In T a b l e (7. 2) a n u m b e i of d i s c r e t i z a t i o n s  Figure  7. l a  -99-  0*91  U  d  0'91  EES"EI  Z.99'07  Q'8  sixd-n  EEE'S  Z.99'2  O'D  - 100-  o  Q'QZ  O'ST  D'Ol  O'S  LD  ca  j  u  i  •d  SuO  i -  •i—i  i • i • i •  *  in  Ln (M  JI  it Ji JI  LTS2  o'oe  i  0'S7  9IXH-H  O'Ol  o "a  O'S  O'O  Figure  7.Id  • in _,  CD.  <n m CO  r~ LD  CO  01  T -4.0 IT) I  m  r—  T x x x x  in_l i  .0  if it 0.125  X X X X X X X X  0.25  X ¥ X X  0.375  X ¥ ¥ X  0.5 X-flXI5  0.625  cn  y-py  Figure  7.2d  0.0 (  —  15.0  30,0  _J  u-nxis  45.0  60.0 _J  90.0 _J  75.0  4 OQ  CD  CD  H II  J10-JL  45.0  -901-  E O . Q  75.0  00  90.0  -107TABLE  in  space  observed  r  s  1  1  2  2  1  3 4  to be  m  Also was  completed  by  2  2  1  4  2  1  4  2  2  1  4  2  3  1  4  3  3  1  6  2  2  5  8  2  2  5  8  2  2  4  1  4  •2  1  4  1  4  1  3  4  6  was  (6. 2) i s g i v e n , < 0.  problem  K(V^)  Thus  i n this  (3 < j < J-3), and  4  2  and  the  respectively.  z^ = x^ lower  ^,  The  equations near  of the  w(0) = w ( l ) = 0.  to  These  with h = ^ ,  .  difference  and  If m  order  be  multistep method.  e q u a t i o n at x^,  which may and  differ  from  X j ^ a r e the  the  "reflection"  = 1 the c o l l o c a t i o n p o i n t i s Z j =  the c o l l o c a t i o n p o i n t s f o r the a p p r o x i m a t i o n  (1 < i < 4).  order  accuracy  4  were  d e f i n e d i n . T a b l e (7.2) c a n  equations at X j  of t h o s e at x ^  In e x p e r i m e n t  = X w(x),  the m a i n  The  Xj.  X^  h -*• 0 nop s i g n i f i c a n t i n c r e a s e i n i t s v a l u e  e q u a t i o n s at x^  x^  o  a l l eigenvalues converged  s t r o n g l y A ( a)-stable  table a r e  m  the e i g e n v a l u e s  were performed  main approximation. and  for which  In fact, Lw(x)  as  m  the d i s c r e t i z a t i o n s  a p p l y i n g any  Listed  1  1  eigenvalue computations  observed.  s  s  r e a l with  computed  r  r  t h o s e of the c o n t i n u o u s algebraic  o  o  of the f o r m  3 < j < J- 3  j = 2  j = 1  #  7. 2  of c o n s i s t e n c y i s g i v e n u n d e r  the b o u n d a r i e s  do  not a f f e c t the  at x^ "o".  are The  overall  schemes.  Example  7.4.  Consider  previous  example.  a g a i n the a l g e b r a i c  eigenvalue problem  L e t the d i s c r e t i z a t i o n i n s p a c e  be  d e f i n e d as  of the follows.  -108F o r each a p p r o x i m a t i o n l e t there be only one c o l l o c a t i o n point (So m = 1.)  = x„  F o r the m a i n difference a p p r o x i m a t i o n take r = s = 3.  this equation i s a seven point f o r m u l a of o r d e r s i x . at x ^ i s defined by the values r = 2 and s = 5.  the " r e f l e c t i o n " of those at x ^ and x ^ .  The a p p r o x i m a t i o n  A t x^ these values  The difference equations at X j  r = 1 and s = 6.  \ ^ w e r e computed for h = ^ , y^- ,  .  are  and X j ^ are once m o r e  Thus the e x t r a  equations near the b o u n d a r i e s are a l s o s i x t h o r d e r .  Thus  difference  A g a i n the  eigenvalues  F o r each h a l l eigenvalues w e r e  r e a l and negative, w i t h the exception of two conjugate p a i r s of c o m p l e x eigenvalues l i s t e d below: h = |-  :  \  = - 1 4 . 9 ± 7. 18 i ,  h =  :  \  h =  :  \  i  =  - 6 6 . 1 ± 24. 8 i ,  \  = - 1 5 0 . 5 ± 55. 8 i ,  i  \  = - 2 0 . 5 ± 6. 50 i , = - 6 7 . 8 ± 25. 1 i ,  X  = - 1 5 0 . 8 ± 56. 0 i .  So apparently there i s a double c o m p l e x root as h —- 0,  that i s a p p r o x i -  m a t e l y g i v e n by  = ~Y ( - ° h  \  Therefore,  6  ±  °-°96 i )  i f the d i s c r e t i z a t i o n i s c o m p l e t e d by applying a m u l t i s t e p  method i n t i m e , tan a >  2  then this method should be s t r o n g l y A{a)-stable  ^6* ' °  r  a  >  2  -l °2  There are strongly A(a)-stable  that do not satisfy this r e q u i r e m e n t . Example 7.5.  with methods  (See V a r a h (1975), p . 19).  C o n s i d e r the diffusion equation of E x a m p l e (7. 3), but now  let the boundary conditions be g i v e n by y (0, t) = y(l> t) = 0.  The a p p r o x i -  m a t i o n to the boundary condition at x = 0 i s a s s u m e d to have the  form  s  -109-  Q  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 t o s . In T a b l e (7.3) i=0 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 n  1  1  U  problem  L^v^ = ^ h ^  v n  h  w  a  observed  s  to have only r e a l and 1  eigenvalues.  The  K ( V ^ ) was  s e e n f o r s m a l l e r h.  a l s o c o m p u t e d and  i d e n t i c a l to those  of T a b l e  (7.3) a n d  j  observed  The  no i n c r e a s e of s i g n i f i c a n c e  Table  (7.4) b e l o w  are  s  i  4  1  2  2  4  1 ' 4  3  6  1  7.3  = 1  j = 2  m  r  4  3 < j <  J-3  m  o  o  r  s  1  4  2  2  1  4  2  1  4  2  2  1  4  3  1  4  3  3  1  6  o  r  s  1  2  2  2  1  4  2  1  4  2  m  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 to have c o m p l e x eigenvalues  e x p e r i m e n t 4 t h e r e a r e two mately  :  o  s  , -^g.  (7.2). TABLE  Two  1  W i t h the e x c e p t i o n of the b o u n d a r y 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  #  1  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 ,  condition number was  negative  are  given i n Table  distinct conjugate  problem  (7.4).  was  In  p a i r s of e i g e n v a l u e s  approxi-  \ j = - ^ ( - 0 . 178 ± 0.049 i ) a n d \ = ( - 0 . 26l± 0.097 i ) . I n 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 c o n j u g a t e p a i r g i v e n by \. = — r (-0.246 ± 0 .044 The  g i v e n by  multistep differencing i n time  with a > a . . — min  The  value a  mm  .  v a l u e s a b o v e i s a l s o g i v e n i n the  s h o u l d t h e r e f o r e be  strongly  as c o m p u t e d f r o m the c o m p l e x r  table,  r  A(a)-stable eigen°  i),  -110TABLE  7.4  j = 1  #  o  s  J = 2  3 < j <  J-3  r  s  m  o  r  s  m  o  r  s  m  o  a  min  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 i m p r o p e r a p p r o x i m a t i o n to a boundary c o n d i t i o n .  The equation  under  consideration is  y (x, t) = ( J^T ) fc  y x x  (x, t),  0 < x < 1 ,  t > 0 ,  with initial condition  3TT  y(x, 0) = cos(—r- x) ,  and boundary  0 < x < 1 ,  conditions  y ( 0 , t ) = y ( l , t) = 0 , x  -1 The s o l u t i o n i s y(x, t) = e  t > 0 .  3TF  cosC-^- x).  L e t the d i s c r e t i z a t i o n i n space be  given by  K^w  (x^  t)  = L ^ w f X j , t)  ,  1 < j  <  J-1  2 where K w . = h  Wj.j  + Q  W j  , +  w.  + 1  and L w . = h  --j(w. _-2w.+w i  j +  .  -1 11T h e o r d e r of c o n s i s t e n c y of this a p p r o x i m a t i o n i s four. A one p a r a m e t e r f a m i l y of a p p r o x i m a t i o n s  (See E x a m p l e (I. 4. 4). )  to the boundary c o n d i t i o n at  x = 0 i s g i v e n by  88w  - (5a + 144)w + (8a +72)w  Q  1  2  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 + 24. equal to two,  except i f a = 0,  Q  3  =0  3  The o r d e r of c o n s i s t e n c y i s  when the o r d e r i s three.  i s t i c p o l y n o m i a l of the d i s c r e t e  c (z) = 88z  - (3a +16)w  The  character-  boundary c o n d i t i o n i s given by  - (5a +144)z  2  If a > .24 then the c h a r a c t e r i s t i c  + (8a + 72)z - (3a + 16).  equation c^fz) = 0 has a r e a l root  > 1,  F o r such a the eigenvalue p r o b l e m K ^ v ^ = X ^ L ^ v ^ has a p o s i t i v e e i g e n value that i s a s y m p t o t i c a l l y g i v e n by 12(z - l ) h  h (z 2  2  2  + lOzj +  T h i s m a y lead to an i n s t a b i l i t y .  1)  (See E x a m p l e (6.5).)  F o r example,  the d i s c r e t i z a t i o n i n time be the fourth o r d e r method of G e a r . (1971),  p . 217.)  TTST  K  T h i s leads to the difference  h( % - *K' 25  1  +3  K" - 'K" 2  let  (See Gear  scheme  3 +  K*)  =VS•  2 Of the r e a l a x i s only the i n t e r v a l [0, 10—] does not l i e i n the s t a b i l i t y r e g i o n of this m u l t i s t e p method i n t i m e .  (See Gear (1971),  p . 216.)  Two  n u m e r i c a l computations w e r e p e r f o r m e d w i t h a = 30 and A t = 0 . 0 2 5 . The f i r s t c o m p u t a t i o n w i t h h = 0.050 and the second w i t h h = 0 . 0 2 5 . The  values of A t X ^ for these data a r e equal to 2 . 0 3 and 8.12 Hence both schemes are unstable.  The a p p r o x i m a t e  respectively.  s o l u t i o n at  time  t = 0. 5, interpolated w i t h a cubic s p l i n e , as w e l l as the exact s o l u t i o n a r e shown i n F i g u r e s (7. 3a, b). a = 20.  T h i s i s stable.  The same computation was repeated w i t h  R e s u l t s at t = 0. 5 a r e  shown i n F i g u r e s ( 7 . 4 a , b ) .  Figure  7.3a  •J .0  U-RXIS  0.0  5.0  2.0  >H ri-  o  II  ll  II  o  O  o  o ro  O ro  o o  *1  P CD  cr  cn ro Ln  -1,0  L.O  2.0  3.0  4.0  5.0  U-RXIS -o.o CJ  (TQ  0)  _  ST T  -  -9TT-  -117BIBLIOGRAPHY  G.  B i r k h o f f a n d S. G u l a t i , O p t i m a l f e w - p o i n t d i s c r e t i z a t i o n s o f 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 . N u m e r . A n a l . 11, 1974, pp.700-728.  C.  de B o o r a n d B. S w a r t z , C o l l o c a t i o n a t G a u s s i a n p o i n t s , S I A M J . N u m e r . A n a l . 10, 1973, pp. 5 8 2 - 6 0 6 .  J.  H.  B r a m b l e a n d B . E . H u b b a r d , O n a f i n i t e d i f f e r e n c e a n a l o g u e of a n elliptic boundary value p r o b l e m which i s neither diagonally d o m i n a n t n o r of the n o n n e g a t i v e t y p e , J . M a t h , a n d P h y s . 43, 1964, pp. 117-132.  J.  C.  Butcher, Implicit Runge-Kutta processes, pp. 50-64.  P.  G.  C i a r l e t , M . H. S c h u l t z a n d R. S. V a r g a , N u m e r i c a l m e t h o d s o f high order accuracy for nonlinear boundary value problems, N u m e r . M a t h . 9, 1967, pp. 3 9 4 - 4 3 0 .  L.  Collatz,  G.  Dahlquist,  S t a b i l i t y a n d 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, Kungl. Tekn. Hogsk. S t o c k h o l m , no. 130, 1959, 87 p p .  G.  Dahlquist,  A special stability problem for linear multistep B I T , 3, 1963, pp. 2 7 - 4 3 .  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 stiff 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 , S I A M J . N u m e r . A n a l . 11, 1974, pp. 3 2 1 - 3 3 1 .  C.  W.  Gear,  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 gewohnlicher Differenzengleichungen, Numer. Math., 1970, pp. 15-48.  David  M a t h , of C o m p .  18, 1964,  Numerische Behandlung von Differential-gleichungen, SpringerV e r l a g , B e r l i n , I960.  N u m e r i c a l initial value problems i n o r d i n a r y e q u a t i o n s , P r e n t i c e - H a l l , 1971.  methods,  differential  Gottlieb and B e r t i l Gustafsson, G e n e r a l i z e d Du F o r t - F r a n k e l methods for 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. 7 5 - 5 , N A S A ' s L a n g l e y R e s e a r c h C e n t e r , H a m p t o n , V A .  B. Hakberg,  Uniformly maximumnorm 10, 1970, p p . 2 6 6 - 2 7 6 .  stable difference schemes,  B.I.T.  B.  L. Hulme,  D i s c r e t e G a l e r k i n a n d r e l a t e d one s t e p 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 equations, R e s e a r c h report, 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 fori two-point B l a i s d e l l , L o n d o n , 1968.  H.  B. K e l l e r ,  Accurate difference methods f o r linear ordinary differential 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 . N u m e r . A n a l . , 6, 1969, p p . 8-30.  boundary value  problems,  -118H.  B.  Keller,  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 e q u a t i o n s , A . K. A z i z , ed. , A c a d e m i c P r e s s , 1975, pp. 2 7 - 8 8 .  H.  O.  Kreiss,  D i f f e r e n c e approximations f o r boundary and eigenvalue p r o b l e m s f o r o r d i n a r y differential equations, Math. o f C o m p . 26, 1972, pp. 605-624.  P.  Lancaster,  S.  Nordmark,  U n i f o r m s t a b i l i t y o f 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. 3 1 4 - 3 2 5 .  S.  J. Osher,  M a x i m u m n o r m stability for parabolic difference schemes i n half-space, H y p e r b o l i c equations and waves, B a t e l l e Institute, S e a t t l e , W a s h i n g t o n , 1968, pp. 6 1 - 7 5 .  H.  Pade,  I.  G.  T h e o r y of m a t r i c e s , A c a d e m i c  Press,  1969.  S u r 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 d e s 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 , A n n . de L ' E c . N o r . , (3), 9, 1892.  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 G e b i e t e ^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.  Reimer,  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 S I A M J . N u m e r . A n a l . 5, 1968, pp. 725-738.  R.  D.  R i c h t m e y e r a n d K. W, M o r t o n , D i f f e r e n c e m e t h o d s 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 m e t h o d f o 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 e r . M a t h . 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 m e t h o d s f o r l i n e a r t w o - p o i n t b o u n d a r y v a l u e p r o b l e m s , M a t h , of C o m p . 29, 1975, pp. 1-13.  M.  H.  Schultz, Spline Analysis,  Prentice-Hall,  order,  1973.  J.  N.  S h o o s m i t h , 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 m e t h o d f o r the s o l u t i o n of t w o - p o i n t b o u n d a r y v a l u e p r o b l e m s o n a u n i f o r m m e s h , 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 e q u a t i o n s , A . K. A z i z , ed. , A c a d e m i c P r e s s , 1975, pp. 3 5 5 - 3 6 9 .  G.  Strang and  G. J . F i x , A n a n a l y s i s o f the f i n i t e e l e m e n t P r e n t i c e - H a l l , 1973.  B.  K.  T h e c o n s t r u c t i o n a n d 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 a n a l o g s of s o m e f i n i t e e l e m e n t s c h e m e s , r e p o r t , L 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.  Swartz,  method,  Alamos  J.  M.  Varah,  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 t h e m i x e d 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 f o r the h e a t e q u a t i o n , M a t h , of C o m p . 24, 1970, pp. 31-44.  J.  M.  Varah,  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 m e t h o d s f o r p a r a b o l i c equations with boundary conditions, S I A M J. Numer. A n a l . 8, 1971, pp. 598-615.  -119-  J.  M.  Varah,  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 m e t h o d s f o r p a r a b o l i c equations w i t h boundary conditions, S I A M J . N u m e r . A n a l . 8, 1971a, p p . 5 6 9 - 5 7 4 . S t i f f l y stable l i n e a r m u l t i s t e p m e t h o d s of extended o r d e r , T e c h n i c a l R e p o r t , Dept.of C o m p . 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 higher o r d e r stable i m p l i c i t methods f o r solving parabolic 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, 1 9 6 1 , p p . 2 2 0 - 2 3 1 .  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 e q u a t i o n s , J . H. B r a m b l e , ed. , A c a d e m i c P r e s s , 1966, p p . 3 6 5 - 3 7 3 .  R.  Weiss,  T h e a p p l i c a t i o n of i m p l i c i t Runge-Kutta 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 , M a t h , o f C o m p . 28, 1974, pp. 4 4 9 - 4 6 4 .  O. B. W i d l u n d , O n the s t a b i l i t y o f p a r a b o l i c C o m p . 1 9 , 1965, pp. 1-13.  difference  schemes,  Math, of  O. B . . W i d l u n d , S t a b i l i t y o f 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 t h e m a x i m u m n o r m , N u m e r . M a t h . 8, 1966, p p . 186-202. K. W r i g h t ,  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 m e t h o d s , 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, p p . 2 1 7 - 2 2 7 .  

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