Twizell, E H

English

Brunel University Research Archive

                    TR/10/86                                                   July 1986  A family of numerical methods for the numerical solution of Blasius' equation.  E.H. Twizell.                                                                          w9259194   1 1. INTRODUCTION I t  i s  w e l l  k n o w n  t h a t  t h e  g e n e r a l  i n i t i a l  v a l u e  p r o b l e m  o f  o r d e r  N  given by ,0xx),)1N(y,...,'y,y,x(f)x((N)y >−=     (1)  wi th  ini t ia l  condi t ions   ),1N,...,2,1r(rz)0x()r(y,0y)0x(y −===     (2)  may be expressed  as the system of N first order equations  ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡+−=−)Nu,...,2u,1u,x(f0..00Nu1Nu...2u1u0Nu' 1Nu..'2u'1u10......1010 (3) where u1 (x) = y(x) and us (x)=y(s - 1 )   (x) for s  = 2,  .  .  .  ,N.  The init ial  conditions associated with (3) are given by     Τ−= ]1Nz,...,2z,1z,0y[)0(u ,     (4) T  d e n o t i n g  t r a n s p o s e .   C l e a r l y ,  t h e  s o l u t i o n  v e c t o r  i s  u =  u ( x )  .]Nu,...,2u,1u[T=  A  p a r t i c u l a r  e x a m p l e  o f  ( 1 )  i s  B l a s i u s '  e q u a t i o n         0xx,)x("y)x(y21)x('"y >−=        (5) w i t h  I n i t i a l  c o n d i t i o n s   .2z)0x("y,1z)0x('y,0y)0x(y ===         (6) Blasius’  equat ion was or iginal ly  a  boundary value problem: numerical  results were reported by Howarth [1] and the problem was formulated as a n  i n i t i a l  v a l u e  p r o b l e m i n  a  r e c e n t  p a p e r  b y  R a d o k  a n d  C h a n  [ 3 ] .  Al l ied to  (3) ,  the  Blasius  ini t ia l  value problem can be expressed as  the  f i rs t  order  system  2  ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡2u1u21003u2u1u000100010'3u'2u'1u        (7) with initial conditions Τ=≡ ]2z,1z,0y[0u)0(u                     (8) o r  as    0u)x(u;)u(F)x(uM)x('u)x(uD =+=≡                     (9) where D = diag{d/dx} is of order three,  ,000100010M⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=    a n d   .]3u1u21,0,0[)u(F Τ−=  ]u,u,u[)x(uu 321== 2. THE FAMILY OF NUMERICAL METHODS Numerical  methods for  the  solut ion of  the  Blasius  problem (5) ,  (6)  wil l  be  obtained by approximat ing the exponent ia l  term in  ).x(u)hDexp()hx(u =+       (10) In  equat ion (10) ,  h  is  an increment  in  x  so that  the  independent  var iable  x  i s  d i sc re t i sed  in  such  a  way tha t  x n =x 0 +nh (n=0,1 ,2 , . . . ) .  Using the approximation ],hD)1(I[)hDI(~)hDexp( 1 θ−+θ−− −        (11) whe re   i s  a  p a r a me t e r  a n d  I  i s  t he  i d e n t i t y  ma t r i x  o f  o r d e r  10 ≤θ≤t h r e e ,  l e a d s  t o  )x(U]hD)1(I[)hx(U)hDI( θ−+=+θ−       (12) which, using (9),  gives the numerical method  .nFh)1(nU]hM)1(I[1nFh1nU)hMI( θ−+θ−+=+θ−+θ−     (13)  3  In (13) )nx(UnU = denotes the solution of the numerical method and )nU(FnF = with  n  =  0 ,1 ,2 ,… Equation (13) describes a family of numerical  methods for the s o l u t i o n  o f  ( 9 )  a n d  t h u s  o f  t h e  B l a s i u s  i n i t i a l  v a l u e  p r o b l e m  ( 5 ) ,  ( 6 ) .  Choos ing   =  0  i s  equiva len t  to  the  use  of  the  (0 ,1)  Padé  approximant  θto  rep lace  exp(hD)  in  (10)  and  leads  to  the  wel l  known expl ic i t  Euler  me t h o d  f o r  s o l v i n g  ( 9 ) ;   c h o o s i n g  θ = 1  i s  e q u i v a l e n t  t o  t h e  u s e  o f  t h e  (1 ,0)  Padé  approximant  to  exp(hD)  and  leads  to  the  fu l ly  impl ic i t  o r  backward Euler method for solving (9);  and choosing θ  = 21   is equivalent to  the  use  of  the  (1 ,1)  Padé  approximant  to  exp(hD)  and  g ives  the  modi f ied  Euler  or  t rapezoida l  ru le  for  so lv ing  (9) .  Implementing the family of methods is  easy; with =0 the solution θa t  x n + 1  i s  g ive n  by   ,n2hUn1U1n1U +=+   ,n2hUn1U1n2U +=+         (14)  nUn1hUnU1nU 32133 −=+   f o r  n  =  0 ,1 ,2 , . . .  .  Wi th  0  <  θ  ≤  1  t h e  s o l u t i o n  i s  o b t a i n e d  b y  s o l v i n g  a  nonl inear  a lgebra ic  sys tem of  the  form 0)U( 1n1n =φ≡φ ++ ,         (15) where 0  is the zero vector of order three and the elements of 1n +φ a r e   ,hU)1(UhUU)U( n2n11n21n11n11n1 θ−−−θ−=φ=φ ++++  ,hU)1(UhUU)U( n3n21n31n21n21n2 θ−−−θ−=φ=φ ++++                               (16) .UhU)1(UUhUU)U( n3n1n31n31n11n31n31n3 2121 θ−+−θ+=φ=φ +++++   The elements   ( r ,s  = 1,2,3)  of  the Jacobian of  1nrsJ+ 1n+φ are   01n13J,h1n12J,11n11J =+θ−=+=+   h1n23J,12n22J,01n21J θ−=+=+=+         (17)  1n1hU11n33J,01n32J,1n3hU1n31J 2121 +θ+=+=++θ=+  4  and the solut ion vector  1nU +  i s  found from (15)  by applying the Newton-Raphson method     ,J1n)k(1n)k(z)k(++ φ−=    (18)  ,)k(zUU1n)k(1n)1k( += +++    (19) i n  w h i c h  t h e  s u b s c r i p t  k =  0 , 1 , 2 , . . .  i s  t h e  i t e r a t e  n u m b e r .  I n  ( 1 8 )  a n d  ( 1 9 )  )k(z  i s  t h e  c o r r e c t i o n  v e c t o r  w h i c h  i s  d e t e r mi n e d  b y  s o l v i n g  t h e  l i n e a r  a l g e b r a i c  s y s t e m  ( 1 8 ) .   I t  i s  a n  e a s y  m a t t e r  t o  w r i t e  a  computer  program to implement  (13)  f rom (14)  or  (16) ,  (17) ,  (18) ,  (19) .    The local  t runcation error  vector  ]h;)x(u[L  of (12)/(13) is  givenby  )x(u...}4D4h)(3D3h)(2D2h){(]h;)x(u[L61241216121 +θ+θ−+θ−= −  (20)  f rom which  i t  fo l lows  tha t  the  method  i s  f i r s t  o rder  accura te  wi th  e r ro r  constant C2 = θ−21  (see Lambert [2]) for .21|=θ  but is second order accura te  wi th  e r ror  cons tan t  C 3  =- 121  for  21=θ   The  g loba l  t runca t ion  vector  ]h;)x(u[G is given by  ).x(u...}Dh)(Dh)(hD){(]h;)x(u[G 4361241322161221 +θ−+θ−+θ−=   (21)  S tab i l i ty  o f  the  fami ly  o f  methods  (13)  i s  inves t iga ted  in  re la t ion  to  the  usua l  s ing le  t e s t  equa t ion  00 u)x(u;)x(w)x('w =λ=      (22)  in  which  is  real .  Equat ion (10)  thus becomes 0<λ )x(w)hexp()hx(w λ=+     (23) and it follows from (11), using W to represent the solution of a numer ica l  method ,  tha t  nWh1h)1(11nW ⎟⎠⎞⎜⎝⎛λθ−λθ−+−+    (24) Def in ing  hh =  ,   so  tha t  0h <  ,   the  te rm       h1h)1(1S θ−θ−+=θ     (25) 5  is  var iously cal led the ampli f icat ion factor  or  ampli f icat ion symbol  or  symbol  of  the numerical  method def ined by (10) ,  (11) .   In  the case of    (9 ) ,   i s  usua l ly  t aken  to  be  the  rea l  pa r t  ( a ssumed  nega t ive )  o f  an  λe igenva lue  o f  the  Jacob ian  o f )u(F)x(uM + .   For  λ< 0  rea l ,  the  s t ab i l i ty  interval  of  the numerical  method is  the  range of  values  of  h  for  which | |  <  1 .  I t  i s  e a s y  t o  v e r i fy  t h a t  ( 1 2 ) / ( 1 3 )  h a s  a  f i n i t e  s t a b i l i t y  θSin te rva l  fo r  210 <θ≤  and  tha t  (12) / (13)  i s  abso lu te ly  s t ab le  o r  A-s tab le  f o r  121 ≤θ≤ .  F u r t h e r m o r e ,  i t  i s  e a s y  t o  v e r i f y  t h a t  θ  =  1  i s  t h e  o n l y  L-stable  member of  the family descr ibed by (13) .   Clear ly  the methods a re  app l i cab le  a l so  to  genera l  in i t i a l  va lue  p rob lems  o f  h igher  o rder .   3.   GLOBAL EXTRAPOLATION OF THE NUMERICAL SOLUTION Introducing a  s l ight  change of  notat ion,  suppose that  the  discret izat ion of  x  used so far  is  cal led Grid 1 and consis ts  of  the points  )1(nx  (n=0 ,1 ,2 , . . . )  and  suppose  fu r the r  tha t  the  numer ica l  so lu t ion  i s  sough t  a t  some  f ixed  po in t   Here ,  the  supersc r ip t  r e fe r s  to  .Qh0x)1(QxX +==Grid 1 and the chosen numerical method is used Q times to integrate from x =X 0  to  x  = x0+ Qh.  The global  t runcat ion error  a t  X is  given by (21) .   Suppose now that  the interval  of  integrat ion is  divided into 2Q subintervals each of width 21 h giving a second discretization to be called Grid 2 consisting of the points ).Q2,...,1,0(ih210x)2(ix +=  Clearly the points  (r  =0,2,4, . . . ,2Q).  of Grid 2 are coincident with the points ( )2rxof Grid 1 .   The global error for Grid 2 has the form  ).x(u...}Dh)(Dh)(hD)({]h);x(u[G 436124181322161412212121)2( +θ−+θ−+θ−=  (26)  S u p p o s e ,  w i t h  a n o t h e r  s l i g h t  c h a n g e  o f  n o t a t i o n ,  t h a t  Q)1(U  a n d  Q2)2(U  a r e  t h e  s o l u t i o n s  o b t a i n e d  a t  X  o n  G r i d s  1  a n d  2  r e s p e c t i v e l y .  Then i t  may be shown that,  when ,21|=θ  the globally extrapolated solution Q)1(U)1(Q2)2(U)E(U α−+= ,   (27) 6 where  is  some parameter,  is  second order accurate provided  = 2,  for α αthen  the  t e rm in  h  in  the  assoc ia ted  g loba l  ex t rapo la ted  e r ro r  func t ion  )1(G)1()2(G)E(G α−+α=          (28) vanishes  (here  )1(G  i s  g iven  by  (21) ) .  I t  i s  easy  to  show tha t  the  g loba l  ex t rapo la t ion  p rocedure  jus t  descr ibed ,  when appl ied  to  the  method wi th  21=θ ,  g ives  a  th i rd  o r d e r  m e t h o d  p r o v i d e d  34=α  Global extrapolation thus improves the numerical methods developed  in §2 for the numerical solution of Blasius’ equation by one order of accuracy .   4 . NUMERICAL RESULTS Numerical  resul ts  were obtained for  the  Blasius  ini t ia l  value problem 33206.0)0("y,0)0('y)0(y;"yy"'y 21 ===−=         (29) which was used by Radok and Chan [3], Four numerical methods were used; these were  Method A: the first  order L-stable method obtained by writing θ  = 1      in  (11) ,   Method B:  the second order extrapolation of Method A,  Method C: the second order A-stable method obtained by writing 21=θ      in (11),  Method D:  the third order extrapolation of Method C. The solution vector Τ= ]"y,'y,y[U was computed for x = 0(h)9 using step sizes h = 1, 0.5, 0.25, 0.2, 0.1, 0.05, 0.025, 0.02, 0.01. Following Radok and Chan [3], the results of Haworth [1] were taken to be the theoretical  solution of the problem. Haworth showed that the solution y(x) increases l inearly (approximately) for x > 4,  that  y '(x) →1 (from below) as x increases, and that y"(x) →  0 (from above) as x increases  7 indicating that Methods A and B should model the behaviour of y"(x) fairly accurately. It was found that Methods A,B,C,D each gave accurate representations of all  three components y,  y ' ,  y" of the solution vector U , the accuracy increasing as h was decreased. Comparing the results obtained with those of Haworth [1] and the best results of Radok and Chan [3] which are given in the Appendix, i t  was found that for x  and 5x ≥ ( i)  for h=0.01, Method A gave results for y(x) which were almost as  accurate  as  those of  Radok and Chan,  while  the resul ts  for  y1(x)   and  y" (x)  were  more  accura te .   Resu l t s  fo r  h -0 .01  a re  g iven  in   T a b l e  1 ,   (i i)  for h ≤  0.25, Method B gave results for y(x)> y '(x),  y"(x) which  were closer to those of Haworth than those of Radok and Chan.  Resul ts  for  h=0.25 are  given in  Table  2;   (iii) for h ≤  0.25, Method C, which is of the same order as Method B,  also gave more accurate results than the method reported by Radok  and Chan in [3];  ( iv )  for all  values of h tested, Method D gave results for all  of  y, y', y" which were closer to those of Haworth than those of Radok  and Chan.  Results for h =1 are given in Table 3. Overall, the methods proposed in the present paper were found to give more accurate results for larger values of x than the method of Radok and Chan [3], when tested on the model Blasius problem (29). Following Radok and   Chan, all results are given to five decimal places .  5. SUMMARY A family of numerical methods has been developed for the numerical solution of the Blasius initial value problem y'' '(x) = 21− y(x) y"(x), x >x0 with y(x0), y'(x0), y"(x0) given.  The methods were tested on a problem from the  8 l i t e ra tu re  and  were  seen  to  g ive  good  accuracy  fo r  h igher  va lues  o f  x .   The family of methods may be used to solve general initial value p r o b l e m s  o f  h i g h e r  o r d e r .   REFERENCES 1.  L. Howarth, On the solution of the laminar boundary layer equations,  Proc. Roy. Soc. London A 164 (1938) 547-579. 2.  J. D. Lambert, Computational Methods in Ordinary Differential Equations,  John Wiley and Sons, Chichester, 1973. 3.  R. Radok and P.S. Chan, A linearizing algorithm for nonlinear differential  equations, Comp. Meth. Appl. Mech. Engng. 54 (1986) 245-253.                        9 Table 1:  Values of y, y'  and y" for x = 0(1)9 using h = 0.01    X A  B C D y           y'           y" 0 1 2 3 4 5 6 7 8 9  0 1 2 3 4 5 6 7 8 9  0 1 2 3 4 5 6 7 8 9 0.00000 0.16719 0.65266 1.39890 2.30551 3.27975 4.27258 5.25868 6.26514 7.26162  0.00000 0.32965 0.62882 0.84373 0.95223 0.98799 0.99543 0.99640 0.99648 0.99649  0.33206 0.32274 0.26597 0.16060 0.06403 0.01599 0.00246 0.00023 0.00001 0.00000  0.00000 0.16557 0.65004 1.39683 2.30578 3.28331 4.27966 5.27928 6.27925 7.27925  0.00000 0.32978 0.62977 0.84605 0.95552 0.99154 0.99898 0.99992 1.00000 1.00000  0.33206 0.32301 0.26675 0.16136 0.06423 0.01591 0.00240 0.00022 0.00001 0.00000 0.00000 0.16557 0.65003 1.39681 2.30576 3.28329 4.27964 5.27927 6.27925 7.27925  0.00000 0.32978 0.62977 0.84605 0.95552 0.99155 0.99898 0.99993 1.00000 1.00000  0.33206 0.32301 0.26675 0.16136 0.06423 0.01591 0.00240 0.00022 0.00001 0.00000 0.00000 0.16557 0.65003 1.39682 2.30576 3.28330 4.27965 5.27927 6.27925 7.27926  0.00000 0.32978 0.62977 0.84605 0.95552 0.99155 0.99898 0.99993 1.00000 1.00000  0.33206 0.32301 0.26675 0.16136 0.06423 0.01591 0.00240 0.00022 0.00001 0.00000  10 Table  2:   Values  of   y ,  y '  and y"  for  x=0(1)9 using h=0.25  X A B C D y           y'           y"  0 1 2 3 4 5 6 7 8 9  0 1 2 3 4 5 6 7 8 9  0 1 2 3 4 5 6 7 8 9 0.00000 0.20461 0.70876 1.43459 2.28331 3.18142 4.09391 5.00949 5.92556 6.84170  0.00000 0.32478 0.60243 0.78677 0.87590 0.90676 0.91449 0.91592 0.91612 0.91614  0.33206 0.31559 0.24794 0.14510 0.06049 0.01795 0.00386 0.00061 0.00007 0.00001 0.00000 0.16633 0.65385 1.40410 2.31380 3.28978 4.28396 5.28140 6.27919 7.27701  0.00000 0.33062 0.63149 0.84652 0.95404 0.98939 0.99681 0.99776 0.99782 0.99782  0.33206 0.32339 0.26632 0.15996 0.06357 0.01589 0.00237 0.00017 0.00001 0.00000 0.00000 0.16543 0.64903 1.39434 2.30214 3.27937 4.27585 5.27569 6.27590 7.27613  0.00000 0.32951 0.62893 0.84515 0.95521 0.99168 0.99921 1.00016 1.00023 1.00023  0.33206 0.32276 0.26659 0.16167 0.06459 0.01596 0.00236 0.00021 0.00001 0.00000 0.00000 0.16557 0.65003 1.39682 2.30576 3.28329 4.27965 5.27927 6.27925 7.27925  0.00000 0.32978 0.62977 0.84605 0.95552 0.99155 0.99898 0.99993 1.00000 1.00000  0.33206 0.32301 0.26675 0.16136 0.06423 0.01591 0.00240 0.00022 0.00001 0.00000  11 Table  3:  Values  of  y ,  y '  and y" for  x  = 0(1)9 using h=1.   X A B C D y           y '            y"   0 1 2 3 4 5 6 7 8 9  0 1 2 3 4 5 6 7 8 9  0 1 2 3 4 5 6 7 8 9 0.00000 0.29001 0.78805 1.40818 2.08803 2.79280 3.50663 4.22336 4.94093 5.65873  0.00000 0.29001 0.49804 0.62013 0.67985 0.70477 0.71383 0.71674 0.71757 0.71779  0.33206 0.29001 0.20804 0.12208 0.05973 0.02492 0.00905 0.00291 0.00084 0.00022 0.00000 0.18896 0.71417 1.48225 2.37882 3.32546 4.28658 5.25073 6.21528 7.17982  0.00000 0.34249 0.64097 0.83502 0.92604 0.95669 0.96384 0.96475 0.96466 0.96456  0.33206 0.32411 0.25501 0.14668 0.05812 0.01528 0200184 0.00051 0.00040 0.00000 0.00000 0.16278 0.63417 1.35891 2.25021 3.22229 4.22080 5.22435 6.22829 7.23222  0.00000 0.32557 0.61721 0.83228 0.95031 0.99386 1.00316 1.00394 1.00393 1.00392  0.33206 0.31908 0.26420 0.16594 0.07012 0.01699 0.00161 0.00004 0.00000 0.00000 0.00000 0.16570 0.65000 1.39634 2.30514 3.28281 4.27915 5.27866 6.27855 7.27846  0.00000 0.32973 0.62954 0.84588 0.95557 0.99151 0.99886 0.99983 0.99991 0.99991  0.33206 0.32298 0.26675 0.16146 0.06417 0.01585 0.00246 0.00022 0.00001 0.00000  12 APPENDIX Values of y, y’ and y" for x=0(1)9;  H :Howarth's solution [1],  RC : Radok and Chan [3].  X H RC y           y'           y" 0 1 2 3 4 5 6 7 8 9  0 1 2 3 4 5 6 7 8 9  0 1 2 3 4 5 6 7 8 9 0.00000 0.16557 0.65003 1.39682 2.30576 3.28329 4.27964 5.27926 6.27923 7.27923  0.00000 0.32979 0.62977 0.84605 0.95552 0.99155 0.99898 0.99992 1.00000 1.00000  0.33206 0.32301 0.26675 0.16136 0.06424 0.01591 0.00240 0.00022 0.00001 0.00000 0.00000 0.16557 0.65003 1.39682 2.30576 3.28290 4.27487 5.26736 6.26023 7.25318  0.00000 0.32978 0.62977 0.84605 0.95552 0.98955 0.99193 0.99276 0.99295 0.99298  0.33206 0.32301 0.26675 0.16136 0.06419 0.00646 0.00152 0.00023 0.00003 0.00000  

A family of numerical methods for the numerical solution of Blasius' equation

Abstract

Similar works

Full text

Available Versions

Brunel University Research Archive