A fourth order convergent finite difference method is developed for the numerical solution of the nonlinear fourth order boundary value
problem y(iv)(x) = f(x,y), a<x <b, y(a) = A0 , y"(a) = B0 , y(b) = A1,
y" (b) = B1 .
The method is based on a second order convergent method which is used on two grids, fourth order convergence being obtained by considering a linear combination of the individual results relating to the two grids.
Special formulas are developed for application to grid points adjacent to the boundaries x = a and x = b , the principal parts of the local truncation errors of these formulas being the same as that of the second order method used at other points of each grid.
Modifications to these special formulas are noted for problems with boundary conditions of the form y (a) = Ao , y'(a) = Co , y(b) = A1,
y'(b) =c1

Twizell, E H

English

Brunel University Research Archive

                       TR/12/85           May   1985 A   two - grid,   fourth   order  method   for nonlinear   fourth   order   boundary   value problems. E. H.   Twizell 0 ABSTRACT A  fourth  order   convergent   finite   difference  method   is   developed   for the  numerical   solution  of   the  nonlinear   fourth   order   boundary    value problem  y(iv)(x)  = f(x,y),  a < x  < b ,   y(a)  =  A0 , y"(a)  = B0 ,  y(b)   =  A1, y" (b) = B1 . The  method is based  on a  second order  convergent   method  which  is used   on   two   grids,   fourth  order convergence  being  obtained  by considering a   linear  combination  of   the   individual   results   relating   to   the   two   grids. Special   formulas   are  developed   for  application   to   grid  points adjacent  to  the boundaries  x  = a  and  x   =  b ,   the  principal   parts   of   the local   truncation   errors   of   these   formulas  being  the  same  as  that  of   the second   order  method  used  at  other points  of  each  grid. Modifications   to  these  special  formulas  are  noted  for  problems with  boundary  conditions  of  the  form y (a)  = Ao , y ' ( a )  = Co , y(b) = A1, y'(b) =c1. Z1488908 1.       THE  SECOND  ORDER  METHOD Consider   the   general   non-linear   fourth  order  boundary  value  problem  given by y(iv) (x) = f(x,y)  ,  a < x < b , a,b,x ∈ RI  (1) with  the   functional  and   second  order  derivative  boundary  conditions y(a)   =  Ao , y"(a) = B o  ,  y(b) = A1 , y"(b) =B1 . (2) It   is  assumed  that   f(x,y)   is   real   and  continuous   on   [a,b]   with   ∂f/∂y < 0, and   that  A0,   A1,  B 0  , B1   are  real   finite  constants.     For  a  detailed discussion  of   the   existence  and   uniqueness   of   the  real  valued  function y(x)  which  satisfies   (1)   and   (2) ,    the  reader  is   referred   to  Agarval and   Akrivis   [1]. Suppose  the   interval  x є [a ,  b]    is  discretized   into  N+  1   subintervals each  of  width  h = (b- a) /( N+1),  where N ≥ 3 is a positive  integer.   The solution   y(x)    will   be   computed  at   the  points  x n  =  a + nh  (n = 1,2,...,N) and   the  notation  y    will  be  used  to  denote   the  solution  of  an  approximat- ing  difference   scheme  at   the  grid  point  xn  ;   clearly  y0 = A0   and  yN + 1 = A1. It  was  noted  by  Twizell   and  Tirmizi   [2]   that  the  solution  y(x) of (1)   and   (2)    satisfies    the    recurrence    relation y(x-2h) - Ry(x-h) + Sy(x) -  Ry(x+h)  + y(x+2h) = 0 , (3) where  the  operators  R  and   S  are  given  by R  =  exp(hD)   +  exp(-hD)   +  exp(ihD)   +  exp(-ihD) (4) and S   =  2 + {exp(hD)  + exp(-hD)}{exp(ihD)  +exp(-ihD)} (5) with  i  =  +√-1   and  D  =  d/dx.     One  of  the  numerical  methods   discussed  in [2]   was  developed  by  replacing  the  exponential   terms   in   (4)   and   (5)    by 2 their   (1,2)   Pade  approximates.     This   leads  to  the  numerical  method (6),0)2nf1n14fn51f1n14f2n(f814h2ny1n4yn6y1n4y2ny=+++++−+−−+++−+−−−where  f s   ≡  f(xs  ,ys   ), s = 1,2,...,N,   which  has  local   truncation   error      ...)n(xyh6480)n(xyh18nt(viii)8119(vi)61(1) −−−=                 (7)at   the  grid  point   X  =  x   . Equation   (6)    is   applicable   only   to  the  N-2  mesh  points xn (n = 2,...,N-1) and   clearly  special   formulas   are  needed  for  the  mesh  points  x1  and  xN. In  order   to  make   possible   the  application  of  the   procedure   to   be  dis- cussed   in   §2,   these   special  formulas  must  be  second order accurate  and their  local   truncation   errors   must   be   of  the  form )i(x(vi)y6h181(1)t i −=  +  0 (h8) for  i =  1, N.  These   requirements  are  met  by   the  pair  of   formulas )8(0f4h1520B2h02A)3f276f1(205f4h36013y24y15y+−=++−+−and (9),1Nf4h1522B2h12A)N205f1-N76f2-N(f4h3601N5y1-N4y2-Ny++−=++−+− for  which   the   local  truncation  errors  are,   respectively,                  ...)1(x(viii)y8h60480521)1(x(vi)y6h181(1)1t −−−=                     (10)and                          ...)N(x(viii)y8h60480521)N(x(vi)y6h181(1)Nt −−−=            (11) The second order  algorithm  on  which   the  procedure of   §2  is  based, is   therefore  described  by   {(8),(6),(9)}   and   the   solution   vector 3 ,T]Ny,...,2y,1[yT](1)Ny,...,(1)2y,(1)1[y(1)Y ≡=   T  denoting   transpose,   is  obtained  by   solving   a   nonlinear   algebraic   system   of   order  N  which  has the  form (1)r)(1)Y((1)fhM4h811(1)YhA =−           (12)In   (12)   the  matrices  Ah,  and Mh ,    are  both  quindiagonal  and   of   order  N and  are  given  by                       (13)⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡−−−−−−−−−−=5414641014641...............1464114640145hAand  ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=125.461.17225.01451141011451141...............1145114111451140225.01.17125.46hM          (14)The   vectors   f(1)  and r (1)   each  have  N   components   and  are   given  by    T](1)Nf,,(1)2f,(1)1[f(1)f K=       (15) where    andN,,...1,2,n),(1)ny,nf(x)(1)Y((1)nf ==  T]1Nf4h1521B2h12A,1N812h10,,...0,,0f814h0A,0f4h1520B2h0[2A(1)r ++−++−+−+−= fA      (16) 4 The  vector  y(1)   =  [y(x1 ,y(x2) ,. . . ,y  (xN)  ] T   clearly   satisfies   the equation    ,(1)t(1)r)(1)y((1)fM4h811(1)yAh +=−         (17) where    T](1)Nt,,(1)2t,(1)1[t(1)t K=  is   the  vector   of   order  N  of   local  truncation   errors.  Defining  ,T](1)Ne,,(1)2e,(1)1[e(1)Y(1)y(1)E K=−=  it   is   seen   that   E(1) satisfies f ( l ) (y ( l ) )   -  f ( 1 )  ( Y ( l ) )  = Fh ( y( l )) E(l )     (18)  where  Fh =  Fh  (y(1))   =  diag is   a   diagonal   matrix   of }(1)iy/(1)if{ ∂∂order  N.      It   follows   that  E(1) satisfies E (1) = Ph t (1) (19) where  Ph.    is   the  matrix  of  order  N  given  by 1hA1)hFhM1hA4h811h(IhP−−−−=       (20) and   Ih,    is   the   identity  matrix  of   order  N. It   is   shown   in  Usmani   and   Marsden   [3]    that          ,hZ2a)(b4384h5||1hA|| −≤∞−        (21) where  Zh, = (b-a)2 + 0.8h2  (the norm referred  to  is   the  L∞     norm  and   from this  point  onwards   the   subscript  will  be   omitted) ,   and  it  is   easy   to   see from   (14)    that    || Mh || ≤ 81.    Using   these  norms,  it was  verified  in  [2] that   {(8),(6),(9)}   is  a  second  order   convergent   method   provided U< 76.8/ {( b- a)2 Zh}            (22) where  U= max|∂f/∂y|.  5 2.      THE  FOURTH   ORDER  METHOD Suppos e  now  that   a   second,   finer,   grid  of  step  size  21 h  is  used.  The interval  a < x < b    is   thus  divided   into  2N+2   subintervals   each  of  width ½h  and   the  points  x n  (n = 1,2,...,N)   of   the  coarse  grid  of   §t  are  named x m (m=  2,4,..,2N)   with   respect   to  the  fine  grid   for  which   the   finite difference  formulas   (8),   (6)   and   (9)  are modifed to   (25),22Nf1204h1B2h4112A)12N2052N76f12N(f57604h12N5y2N4y12Nyand(24)2N,...,1m,0)2mf1m14fm51f1m14f2m(f12964h2my1m4ym6y1m4y2my(23),0f1204h0B2h4102A)3f276f1f(20557604h3y24y15y++−=+++−−++−−−=+++−+−−−−+++−+−−−+−=++−+−The  solution  vector  of   these   2N+2  nonlinear  algebraic  equations will  be  denoted  by  Y(2)  so   that Y(2)  satisfies   the  equation                            .(2)r)(2)Y((2)fh21M4h12961(2)yh21A =−     (26)The   forms  of   the  matrices  A21 H ,  M21 h .  and of  the vector  f (2)  are obvious from (13),  (14)  and (15)   while the  vector  r (2)   is   given  by T]22Nf1204h1B2h4112A,22Nf12962h1A0,,...0,,of12964h0A,0f1204h0B2h410[2A(2)r ++−++−+−+−= (27) The  matrix  analogous   to  Ph,    is  P21 h ,    and  this   is  defined  by   1h21A1)h21Fh21M1h21A4h12961h21(Ih21P−−−−=    (28)in  which  I21 h , is   the   identity  matrix  of  order  2N+1and }(2)iy/(2)if{diagh21F ∂∂=  is  of  order   2N+1. 6 The  local  truncation  errors  of  (23),  (24)  and  (25)  are,  respectively,   ,...)1(x(viii)y8h15482880521)1(x(vi)y6h11521(2)1t −−−=              (29) (30),2N)...,1,(m,...)m(x(viii)y8h1658880521)m(x(vi)y6h11521(2)mt=−−−= )31(,...)12N(x(viii)y8h15482880521)12N(x(vi)y6h11521(2)12Nt −+−+−=+ and   the  theoretical   solution  vector  y (2) ,   relating   to   the  fine   gird, satisfies   the  equation   ,(2)t(2)r)(2)y((2)f21M4h8116(1)yh21A +=− h         (32) where   .T](2) 12Nt,,(2)2t,(2)1[t(2)t += KSuppose  now  that  hh21I is  a  fine-to-coarse  grid   restriction   operator defined  by        .](2)2Ny,..,.(2)4y,(2)2[y(2)Yhh21I =                                     (33) Then  the  components of (2)Yhh21I give second approximations  to y(x) at the  N  points  x n  (n = 1,2,...,N) of the original  (coarse)  grid  used  in  § 1. To  develop  a numerical method which is fourth order convergent, a  para- meter  α  must  be  determined such that the vector E* of  order N defined  by   ,](1)Yα)(1(2)Yhh21I[y(1)Eα)(1(2)Ehh21αI*E ++−=++= α            (34)where  ,T](2) 12Ne,,(2)2e,(2)1[e(2)Y(2)y(2)E +=−= K  has  norm   ||  E*  ||  =0(h4). In  (34)   it  is  clear  that      (2)th21P(2)E =             (35) where .T](2) 12Nt,,(2)2t,(2)1[t(2)t += K  is   the  vector  of  local  truncation errors  given  by   (29),(30)   and   (31),   an  that .T](2)2Ne,,(2)2e,(2)4[e(2)Ehh21I K=  7 It  may  be   shown   from   (21)   that ,h21Z2a)(b424h5||1h21A|| −≤−             (36) where   h21Z   =   (b-a)2   +   0.2h2,   and   it   is   easy   to   see   that    ||  h21M  ||    =81. It   then   follows   from   (20)   and   (28)   that || h21P  || < || Ph ||                                                            (37) Defining  τh   =   ||t(l)   || /h4   and  τh21   =   || t( 2 )   ||/h4,  it   is  clear from  (34) that || E* ||   <   || Ph IK(α  τh21  (1 - α)  Th)                     (38) provided   (22)    is   satisfied.      Noting   from   (7),    (10),    (11)    that || t(1)   ||    =  181 h6V6,   and   from  (29),  (30),  (31)  that   ||  t (l) ||   =  11521   h6V6 , where  V6     =   bxamax≤≤      | y (vi)  (x)|,    it   follows   that ||  E*  ||    =   O(h4)                                          (39) when  a   = 34  .       The   numerical   formulation        (1)Y31(2)Yhh21I34(E)Y −=             (40)  where  Y(E).     is   avector   of  order  N,   is  therefore  a  fourth order convergent method, Undoubtedly   the  main   advantage   of   the   fourth   order  method   just dveloped   is    its   ease   of    implementation,   especially   in   comparison  with other   fourth   order   methods    (in   particular,   multiderivative   methods    [2]). The  novel   method   requires   only   two  applications   of   a   second   order  method, the  solution   being    obtained  by   taking  a   linear  combination  of   the   results relating   to    the    two    applications. 8 3.       FIRST  ORDER  DERIVATIVE   BOUNDARY   CONDITIONS In   this   section,   the  boundary  value  problem  consisting   of   the  differ- ential   equation   (1)   with   the  boundary   conditions y  (a)  =  A0 ,    y'(a)  =  C0 , y(b)  =  A1 , y1(b) = C1 (41) will   be   considered. The  special  formulas  (8)  and  (9)  must be  replaced by (42)0f4h6305303hC0A211)39f2608f1(1079f4h252013y2y2919y++=−+−+−  and (43)1Nf4h6305313hC1A211)N9f1071-N608f2-N(-9f4h25201Ny91-Ny292-Ny+++=++−+−for  which  the   local   truncation  errors  are,   respectively,    )44(,...)1(x(viii)y8h60480271)1(x(vi)y6h181(2)1t −−−=and    )45(,...)N(x(viii)y8h60480271)N(x(vi)y6h181(2)Nt −−−=The   second   order   algorithm   is   then   given  by   {(42),(6),(43)}. The  matrix  Ah,    becomes  ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡−−−−−−−−−−=92914641014641...............14641146401299hA  (46)9 for  which   it   is   known   [3]    that ,hZ4a)(b4384h1||1hA|| −≤−      (47)where,  now,  Zh = 1+8h3 (b- a) - 3 ,  and  the matrix Mh.   becomes ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=280971135684280811451141011451141...............114511411145114028081356842809711hM  (48)for  which  || Mh || = 81   as   before.      The   vector  r (1)   also   requires   modi- fication;      it  becomes T]1Nf4h630531hC31A211,1Nf814h1A0,,...0,,0f814h0A,0f4h630530hC30A211[(1)r ++−++−+−++=         (49) The   restriction  on  U   given   in   (22)    is   no   longer  applicable;      it   must be   replaced  by U< 384/ {(b- a)4  Zh } , (50) where   Zh  = 1 + 8h3(b- a) - 3 . Considering,   next,   the  application  of   the   second   order   algorithm to   the   fine  grid,   it   is   evident   that   equations   (23)   and   (25)   need  modi- fications.      They   become        (51)0f4h100805303hC0A211)39f2608f1(1079f403204h3y2y2919y++=−+−+−and 10 (52),1Nf4h100805313hC1A211)N9f1071-N608f2-N(-9f403204hNy91-Ny292-Ny+++=++−+− respectively,   for  which   the   local   truncation     errors   are )53(,...)1(x(viii)y8h15482880271)1(x(vi)y6h11521(2)1t −−−=and )54(,...)N(x(viii)y8h15482880271)N(x(vi)y6h11521(2)Nt −−−=With   respect   to   the   fine   grid,   the   second   order  method   is   defined  by {(51),(24),(52)}   and   the   elements   of   the  vector   t (2)  are  given  by (53),    (30)   and    (54). The   forms   of   the   revised   matrices   A½h  and  M½h   are   obvious   from (46)    and   (48) ,   and   the   vector   r (2)   is   now   seen   to   take   the   form (55).T]22Nf4h10080531hC231A211,22Nf12964h1A0,,..,.,00f12964h0A,0f4h10080530hC230A211[(2)r++++−+−++=             Clearly,   ,h21Z2a)(b424h5||1h21A|| −≤− where Z21 h = 1+h3(b- a)- 3,  and ||  M21 h  ||  ≤  81,   and  convergence  of  the  fourth  order method  given by  (40) is   established   for   the  boundary  value  problem { ( 1 ) , ( 4 1 ) }   as  in  §2. 4, NUMERICAL   EXPERIMENTS The   fourth  order  method  discussed  in  §2  was  tested  on  the  following problem  from  the   literature. Problem  1 y(vi)(x)   =  6exp[-4y(x)] -  12 (l +x)- 4 , 0 <x <1 with     boundary   conditions 11 y(0)   =  0   ,   y(1)   =  ln2 , y"(0)  = 1 , y"(1)  = -0.25 for  which   the  theoretical   solution   is y(x) = ln(1+x) The   interval    [0,1]   was   divided   into   N+1   equal   parts   each   of   width h = 2- m  with  m=  3,4,5   so   that  N = 2m -  1  =7 ,15,31,   respectively. The  value   of   || y-Y(E)  ||  ,   where  Y(E)   is   defined   in  equation   (40), was   calculated   for   each  value   of  N.      The  numerical   results   are   tabulated in  Table   1   which  also   includes   results   obtained   using   the   fourth   order methods   of  Agarwall   and  Akrivis  [1]  as  well  as  the  second   order  method {(8) , (6),(9)},   with   the  coarse  grid,   on  which   the  novel   method   is   based. Observing   the  contents  of  Table   1,   it   is   evident   that   the  fourth order  method  of   §2   gives   better  numerical   results   than   either  of   the other   fourth  order  methods   tested.      Bearing   in  mind   its   ease   of   imple- mentation,   the   method   is   clearly  an   economic   alternative   to   the   other fourth  order  methods. The  adaptation   in  §3  of   the  novel   fourth  order  method   to  problems with   functional  and   first   order   derivative  boundary   conditions  was also   tested.     The  problem  on  which   this  adaptation  was   tested   is  given  by Problem 2 y(iv)(x) = 6exp[-4y(x)] - 12(l +x)- 4 ,  0 < x < 1 , with  boundary   conditions y(0) = 0 , y(1) =ln2 ,   y'(0)  = 1 , y'(1) = 0.5 , for  which  the  theoretical   solution  is  y =ln(1+x)   as   in  Problem  1.     The three  discretizations   used   for  Problem  1  were  also  used  for  Problem 2 and  the  value  of   || y - Y(E)  ||  was  calculated   in  each  case. The  numerical   results   are  given    in Table  2  from which  it may be 12 observed   that   the  method  of   §3   retains   the  accuracy  attained  for   the problem  with  functional   and   second  order  derivative  boundary  conditions.          Only   three   applications  of    the   Newton-Raphson   method   for   nonlinear algebraic   systems  were  need  to  give  convergence   to  three   significant figures   for  both  problems, 5.        SUMMARY A  fourth  order  convergent  finite  difference  method  has  been   developed and  analyzed  for   the  numerical   solution  of   the  general   fourth   order nonlinear  boundary  value  problem  y(iv)  (x)   =  f(x,y)   with  boundary  con- ditions  given   in   the  form  of   two  functional  values   together  with (i)   two   second   order  derivatives,   or    (ii)   two   first   order   derivatives. The  method  was  based on a  second order convergent method which was  used  on  two  grids,  fourth  order  convergence  being  obtained by considering  a  linear  combination  of   the   results   determined   for  each grid   individually.      Special   formulas  were  developed  for   application to  points   adjacent   to   the  boundary,   the  principal   parts  of   the   local truncation  errors   of   these  formulas  being   the   same  as  that  of  the second  order  method  which  was  used  at  other  points  of   the   discretization. REFERENCES [1]      R. P.  Agarwall and G. Akrivis,  Boundary value problem occurring  in plate deflection   theory,   J.   Comput.  Appl.  Math.  8  (1982)  145-154, [2[     E.H.   Twizell  and  S.I.A.   Tiraizi,   Multiderivative  methods  for  non-linear fourth  order  boundary  value  problems,   Brunel University Department  of Mathematics  and Statistics  Technical  Report TR/07/85  (1985), [3]     R.A. Usmani and M.J. Marsden,  Convergence of a numerical procedure for  the  solution  of  a  fourth  order  boundary  value  problem,   Proc. Indian Acad.  Sci.88 (1979), 20-21. 13 Table   1:     Error  norms   for  Problem   1.    second  order    Agarwal  and  Akrivis fourth  order m N method Method  A Method  B method 3 7 0.19E-3 0.14E-4 0.14E-4 0.37E-5 4 15 0.46E-4 0.83E-6 0.83E-6 0.29E-6 5 31 0.11E-4 0.54E-7 0.54E-7 0.19E-7 Table   2:     Error  norms   for  Problem  2.     econd  order fourth  order m N method method 3 7 0.15E-3 0.22E-4 4 15 0.23E-4 0.42E-5 5 31 0.27E-5 0.67E-6  

A two - grid, fourth order method for nonlinear fourth order boundary value problems

Abstract

Similar works

Full text

Available Versions

Brunel University Research Archive