Global extrapolation procedures, in space and time are considered for the numerical Solution of linear partial differential equations. Global extrapolation procedures in time only are reviewed.
The procedures are tested on three problems from the literature, one of which has a nonlinear source term

Twizell, E H

English

Brunel University Research Archive

               TR/11/87          November 1987 G l o b a l  e x t r a p o l a t i o n  p r o c e d u r e s  f o r  l i n e a r  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 .   E.H. Twizell.  (0)   ABSTRACT  Global extrapolation procedures, in space and time are considered for the numerical Solution of linear partial differential equations. Global extrapolation procedures in time only are reviewed.  The procedures are tested on three problems from the literature, one of which  has a nonlinear source  term. (1)  INTRODUCTION  Consider the  linear partial differential equation  (PDE)  ∞<≤<Ω∈+−∂∂ΩΩ Tt0,x,),t,x(fuAtuvv    (1) in which v=l ,2, Ω= ),t,x(uu ,  is a spatial domain in R, R2 or R3 with boundary , and Ω∂ΩA  denotes a linear differential operator of order q which differentiates the function u with respect to the space variables. Obviously, equation (1) describes a first order hyperbolic equation when v=l and q=l, a second order hyperbolic equation when v=2 and q=2, a second order parabolic equation when v=l and q=2, and a fourth order parabolic equation when v=2 and q=4. Associated with the PDE (1) are the initial conditions Ω∂Ωε−==∂∂ Ux,)1v()1(0r,u)0,x(tu )r(0rr     (2) and the boundary conditions ,Tt0,x,)t(fuA <<Ω∂ε= Ω∂Ω∂        (3) where  is a linear differential operator of order less than q , which also Ω∂Adifferentiates u with respect to the space variables and acts on the boundary Ω∂ . It must be noted that  and  may all depend on Ω∂ΩΩ A,f,A Ω∂A x  . A popular method ( cf  [l ,3,4,5,7,8,9,10,13,14,15,16,17]) of solving the initial-boundary value problem {(1),(2),(3)} is the so-called method of lines (MOL) in which the space domain Ω is discretized in some way and, via some finite difference or finite element approximation to AΩ, the PDE is transformed into a time-continuous system of ordinary differential equations (ODEs) of order v , which has the form (2)  ,Tt0,)t(hf)t(hUhA)t(hUVD <<+=      (4) with the associated initial conditions .)1v()1(0r,)r(0u)0(hUrD −==       (5) In (4) and (5), upper case U has  been introduced to distinguish the theoretical solution of an approximating method from the theoretical solution (lower case u) of the initial-boundary value problem, h is the parameter of a grid in )t(hUΩ∂ΩU  is an N-dimensional vector the elements of which are the approximations to the unknown dependent variable u at the grid points, the N-dimensional vector )t(hf  arises from Ωf  and , the N-vector(s) Ω∂f ))1v()1(0r()r(0u −=  arise from (2), Ah is a time-independent matrix of order N which arises from AΩ and , and D = diag{d/dt} is of order N. Ω∂AThe MOL approach then solves the initial value problem {(4), (5)} by dividing the time-interval 0 < t < T into Q subintervals (time-steps) each of length l  , say, so that , and then employing a k→step (k>v) ODE solver to  integrate  from t=0 to t=T. TQ =lIn the case of PDEs with v=2, the initial conditions are used to determine an approximation to ,)(U h l thus providing enough starting values for a two-step ODE solver. The MOL thus determines the solution in a recursive manner on a grid G1 in [ ] [ ]Tt0 <<×Ω∂ΩU  which has a total of Q+l time levels. Under conditions such as those detailed in [11,Th.4.1] the full global error at each of the N grid points of G1 at time t=T=Q  is given by the quantity El 1, which has the form E1 = hs C  +  Pl B + K          (6) In (6), p is the order of the ODE solver and s is the order of the approximations to AΩ  and  (the reader is referred to the work in [11] on order reduction and to Ω∂A[2,12,18] for the reduction in the accuracy of the time integration in the presence of time-dependent boundary conditions). The quantities C and B are independent of h, l  and T and the quantity K is 0(hs* + l P*) where p* > p and s* > s are integers. (3)  It will be assumed that the space and time increments h and l  satisfy any restriction of the form                (7) ,v/qa,ah =λ<lwhere λ is a  fixed positive constant,  which must be  imposed for stability.   2.  GLOBAL EXTRAPOLATION IN TIME  A number of authors ( cf [3,4,5,7,8,9]) have used local extrapolation methods to integrate (4) (with v=l) and (5) from time t to time t+2 l , t+3 l  or t+4 l . The merits of local extrapolation in time were shown in [17] to be overshadowed somewhat by those of global extrapolation in time. Using a half-time step of length l21 , the time interval 0 < t < T is divided into 2Q subintervals and, if the same space step h is retained, a new grid G2 is constructed which has N interior grid points of G2 at time )21(Q2Tt = l=  is given by the quantity E2 which has the form .KBPP2Csh2E +−+= l         (8) The grid parameters h and l21  of G2 clearly satisfy (7) if the parameters h and  lof G1 do so.  Introducing some new notation, suppose that )T(U 1,h denotes the computed solution vector at time T on G1, and that )T(U 2,h is the associated vector on G2. Consider now the approximation )T(U)1()T(U)T(V 1,h2,h α−+α=        (9) and the associated error Ev defined by .E)1(EE 12V α−+α=         (10) It is easy to show that the term in l P in (10) vanishes when the parameter α takes the value .)12(/11with)12(/2 PP −−=α−−=α       (11) This  g loba l  ex t rapola t ion  in  t ime only  us ing  gr ids  G 1  and  G 2  has  thus  produced an approximation )T(V defined by (9) which is  provided α takes the )Ph(0 *s l+  (4) value given in (11). 3.  GLOBAL EXTRAPOLATION IN SPACE AND TIME Suppose now, that, in addition to halving the time step  , the space parameter h is lalso halved. In the special case where Ω is an interval XεR with N interior points on Grids 1 and 2, X is now divided into 2N+2 subintervals. A third grid G3 is thus constructed which has 2N+1 grid points at each of its 2Q+1 time levels. The full global error at the 2N+1 grid points of G3 at time t=T=2Q( l21 ) is given by the quantity E3 which has the form .KB2ch2E ppss3 ++= −− l         (12) Replacing  and h in (7) by  l21  and 21 h, respectively, gives v/qa,2/h 1aa =λ< −l         (13)  so that in a space-time global extrapolation procedure the parameters h and l  used on G1, must satisfy (13) instead of (7). Suppose, now, that 3,h21U (T) denotes the computed solution vector at time T on G3, then 3,h21U (T) has 2N+1 elements. Let hh21I be an operator which isolates the N elements of 3,h21U (T) corresponding to the N elements of 1,hU (T) and 2,hU (T):  that is, hh21I 3,h21U (T) has N elements. Consider, next, the approximation )T(U)1()T(U)T(,UI)T(W 1,h2,h3h21hh21 γ−β−+γ+β=     (14) and the associated error EW defined by .E)1(EEE 123W γ−β−+γ+β=         (15) It may be shown that the terms in hs and P in (15) vanish when the parameters β land γ take the values ])12()12([/)22(,)12(/2 pspsss −−−=γ−=β      (16) (5) so that .)12(/11 p −−=γ−β−          (17) This global extrapolation in both space and time using grids G1 ,G 2 and G3 has thus produced an approximation )T(W which is provided the parameters β and )Ph(0 **s l+γ take the values given in (16). In special cases when s=p the parameter γ vanishes and thus only two grids, G1 and G 3, are needed to obtain the global space-time extrapolation. A notable example of this is the Crank-Nicolson method, for which s=p=2, for solving second order parabolic equations (v=l,q=2)(see Problem 2, in §4). As noted in [17] it is of course theoretically possible to extrapolate to arbitrarily high orders: this remark is applicable to space-time extrapolation as well as time only [17]. However, in the belief that the extra orders achieved in §§2 and 3 of the present paper are high enough for PDE's, no further extrapolations will be considered. The economics of the extrapolations are easy to compare. When compared with the result 2,hU (T) computed on grid G2, the computation of the time-only extrapolation vector V defined by (9) requires an additional computation effort of 50%. When compared with the result 3,h21U (T) computed on grid G3, the computation of the space-time extrapolation vector W  defined by (14) requires 1.75 times as many operations. This factor is reduced to 1.25 when s=p, for then γ=0 and grid G 2 is not required. 4.  NUMERICAL RESULTS The global extrapolation algorithms outlined in Sections 2 and 3 were tested on three problems from the literature, two from the literature on linear second order parabolic equations and one from the literature on nonlinear second order hyperbolic equations. (The last problem is not, of course, of the type described by (1), but the authors feel that the results obtained were sufficiently interesting to report in the present paper.) Problem 1 (Lawson and Morris [9], Gourlay and Morris [3], Twizeli and Khaliq [15]). ,0t,2x0,xutu22><<∂∂=∂∂  (6) wi th  boundary  condi t ions  u(0 , t )  =  u(2 , t )  =  0 ,  t  >  0  and  in i t ia l  condi t ions  u(x ,0)  =  1, 0 ≤  x ≥  2. This problem is of continuing interest because it has a discontinuity between initial and boundary conditions. The theoretical solution is given by [ ] .)ts41(exp)xs21(sins2)1(1)t,x(u 221ss π−ππ−−= ∑∞= Lawson and Morris [9] transformed this problem into the system  )t(UAdt/)t(Ud h=  by the MOL, where ,2112112112hhA2⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−−−= −             (18) and  then used the fully implicit,  L0-stable method  )t(U)t(U)hAI( =+− ll         (19) to obtain a numerical solution; for (19),  p=l   and s=2. The behaviour of (19) and comparisons of it with other methods are well documented (see, for instance, [3,9,15,17]). The purpose of the present paper is to show that the global extrapolation procedures described in Sections 2 and 3 improve accuracy. Four numerical experiments were carried out in which the space and time increments on grid G1, were given by (h, l )  = (0.1,0.1), (0.1,0.05), (0.05,0.002), (0.05,0.05), respectively, the increments on G2 and G3 following in an obvious manner. It was observed in [9] that (19) produces a smooth solution to Problem 1, with maximum errors occurring where x=l. The errors at this point at time t=1.0 using (19) on grids G1, G2 and G3 individually, the time-only global extrapolation of (18) on G1, and G2 (defined by (9) with α=2), and the space time global extrapolation of (19) on G1, G2 and G3 (defined by (14) with β=4/3 and γ=2/3), are given in Table 1. It is easy to see from Table 1 that the two extrapolation procedures produce the improved accuracy predicted in Sections 2 and 3. Problem 20t,1x0;)2x(exutu 2t22><<+−∂∂=∂∂ −  (7) with initial conditions u(x,0) = l + x2, 0 ≤ x ≤ l, and boundary conditions u(0,t) = 1, u(l ,t) = 1+e-t, t > 0. The theoretical solution is given by u(x,t) = 1 +e-tx2. This problem is the one-space dimensional form of Problem (5.2) in Verwer and de Vries [17].    Using the MOL, Problem 2 may be transformed into the equivalent first order initial value problem  )t(hw)t(UhAdt/)t(Ud +=  in which Ah is as defined  by (18)  and  T)]t(h,Nw,...,)t(hw,)t(hw[)t(hw ,2,1=  (T, here, denoting transpose) with   ,)2h(eh)t(hw2t2,1+−= −−  1N,...,3,2m,)2hm(e)t(h,mw22t −=+−= −  .)2hN(e)e1(h)t(h,Nw22tt2 +−+= −−−The A0 -stable implicit method  )t(U)hA21I()t(U)hA21I( +=+− ll  ,)t(hw)hA21I()t(hW)hA21I21⎥⎦⎤⎢⎣⎡ +++−+ llll   (20) for which p=s=2, was used to compute the solution at time t=l .0. Six numerical experiments were carried out for which the space and time steps on grid G1 were taken to be h, ) = (0.1,0.1), (0.1,0.05), (0.1,0.002), (0.05,0.1), (0.05,0,05) and (0.05,0.02) l(respectively. The errors at the grid points where x = 0.3,0.6,0,9 are given in Table 2 for grids G1 and G3 individually and for the space-time global extrapolation of (20) on G1 and G3 (defined by (14) with β=4/3 and γ=0). The errors on G2 were found to be similar to those on G1, and, consequently, the errors for the time-only extrapolation to be similar to those for the space-time extrapolation experiments (because α=β=4/3 and γ=0). It may be seen from Table 2 that errors were larger in modulus near the time-dependent boundary where some of the errors on G1, had a dominant effect on the errors after extrapolation. (8) Problem 3 0t,bx0,u)u,x(xu)x(gdtu 2222241 ≥≤≤λ+∂∂=∂∂  with boundary conditions  0t,0)0,b(tu)0,x(xu >=∂∂=∂∂  and initial conditions .bx0,)b/x(cos)gd()0,x(tu,)b/x(sin)0,x(u 21≤≤ππ−=∂∂π=   (21) In this problem (which does not conform to the type described by equation (1)), d(x) is a depth function given by d(x) = d*[2+cos(2πx/b)], g is the acceleration due to gravity, and λ(x,u) is the coefficient of bottom friction defined by λ = g | u |/(C2d), C being the Chezy coefficient. Using the MOL, this second order hyperbolic problem may be converted into the second order initial value problem  0t,))t(U,t(f))t(U(M])C4(/g[)t(UhMBgdt)t(Ud 24222>≡ϕ+= −   (22) with initial conditions U(0) and dU(0)/dt derived from (21). In (22), M = diag{di} is a diagonal matrix in which di = d(xi) for i=0,1,...,N, and xi=ih with h=b/N. The matrix Bh is of order N+l and is given by   ;2212112122hB 2h⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−−−= −the vectors U(t) = [U0(t),U,(t),....UN(t)]T, T10 ]N,...,,[ ϕϕϕ=ϕ with  ϕ=ϕ ))t(U(ii = i=0,l,...,N, and f = [f 0,f, ,f 2,...,fN]T are of order N+l.  ),t(U|)t(U| i2iThe simple explicit method U(t+ ) = 2l U(t) – U(t- l ) + l 2 f (t,U(t)) ; t= ,2 l ,...             (23) lfor which p=s=2, was used to compute the solution on grid G1. In using (23) U(0), derived from (21), and U( ) = l U(0) + 2l f(0,U(0)) were used as starting values. Following van der Houwen and Sommeijer [6], the various parameters associated with the problem were given the values (9) h = 10,   31=l  (grid G1) ; b= 100,   g=9.81. d*=10, C=50 and the solution was computed to  time t=3600.  Taking as reference solution that used by van der Houwen and Sommeijer [6], the second order Runge-Kutta-Nystrdm method given by their equation (3.14), the solution in the interval 3567 ≤ t ≤ 3600 at the point x=80 is depicted in Figure 1. The graphs in Figure 1 refer to the reference solution, the solution obtained on G1 using (23), the time-only extrapolation defined by (9) with α=4/3, and the space-time extrapolation defined by (14) with β=4/3 and γ=0. It is seen from Figure 1 that the extrapolation procedures clearly improve accuracy, as predicted in §§2 and 3, a marked reduction in phase-lag being evident (cf van der Houwen and Sommeijer [6]). 5.  SUMMARY This paper has considered global extrapolation procedures, in time only and in space and time, for the numerical solution of linear partial differential equations. The procedures were tested on three problems from the literature, one of which was nonlinear. It was seen that the extrapolation procedure involving both space and time produced notable reductions in error. REFERENCES 1. J. R. Cash, "Two new finite difference schemes for parabolic equations", SINUM, 21(3), 443-446, 1984.  2. G. Fairweather and A. R. Mitchell, "A new computational procedure for ADI methods",  SINUM, 4, 163-170, 1967.  3. A. R. Gourlay and J. Ll. Morris, "The extrapolation of  first order methods for parabolic partial differential equations. II", SINUM, 17(5), 641-655, 1980.  4. A.  R. Gourlay and J.  Ll. Morris, "Linear combinations of generalized Crank Nicolson schemes", IMA J.  Numer, Anal., 1, 347-357,  1981, (10) 5.  A. R. Gourlay and J. Ll. Morris, "The extrapolation of first order methods for variable coefficient inhomogeneous parabolic partial differential equations", Arab Gulf J. scient.  Res., 1(2), 599-622, 1983. 6. P. J. van der   Houwen and B. P. Sommeijer,  "Explicit Runge-Kutta(-Nystrom) methods with reduced phase errors for computing oscillating solutions", SINUM, 24(3), 595-617, 1987. 7. A. Q. M  Khaliq and E. H, Twizell, "Backward difference replacements of the space derivative in first order hyperbolic equations",  Comp. Meth. in Appl. Mechs. and Engng., 43(1),  45-56, 1984. 8. A. Q. M. Khaliq and E. H. Twizell,  "L0 -stable splitting methods for the simple heat equation in two space dimensions with homogeneous boundary conditions", SINUM, 23(3),  473-484,  1986. 9. J. D. Lawson and J. Ll. Morris, "The extrapolation of first order methods for parabolic partial differential equations.I.", SINUM, 15(6), 1212-1224, 1978. 10. J. Oliger, "Fourth order difference methods for the initial-boundary value problem for hyperbolic equations", Math. Comp. 28(125), 15-25, 1974. 11. J. M. Sanz-Serna, J. G. Verwer and W. H. Hundsdorfer, "Convergence and order reduction of Runge-Kutta schemes applied to evolutionary problems  in partial differential equations", Numer. Math., 50, 405-418, 1987. 12. B. P. Sommeijer, P. J. van der Houwen and J. G. Verwer, "On the treatment of time-dependent   boundary conditions in splitting methods for parabolic differential equations", Intern. J. Numer.  Meth, Engng., 17., 335-346, 1981. 13. E. H. Twizell, "An explicit difference method for the wave equation with extended stability range", BIT, 19(3), 378-383, 1979. 14. E. H. Twizell and A. Q. M. Khaliq, "A difference scheme with high accuracy in time for fourth order parabolic equations", Comp. Meth. in Appl. Mech.  Eng., 41(1), 91-104, 1983. 15. E. H. Twizell and A. Q. M. Khaliq, "A family of numerical methods for diffusion and reaction-diffusion equations", Num. Meth. Partial Diff. Equations, 2, 31-45, 1986. 16. E.  H. Twizell and S. I. A. Tirmizi, "A family of methods for the wave equation in one- and two-space dimensions", Num. Meth. Partial Diff. Equations, 1(2), 105-125, 1985. 17. J. G. Verwer and H. B. de Vries, "Global extrapolation of a first order splitting method", SISSC, 6(3), 771-780, 1985. 18. Ye. G. D'Yakonov, "Some difference schemes for solving boundary problems", USSR Comput. Math. Math. Phys., 1, 55-77,  1963. (11) T a b l e 1 .  Errors  U-u where x=1.0 at  t ime t=1.0 for  the four  experiments  o f  P r o b l e m  1 .  Increments on G1  Errors on grids  h l  G1 G2 G3 G1,G2 G1,G2,G30.1       0.1 0.33E-1 0.17E-1 0.16E-1 0.37E-3 0.31E-3 0.1 0.05 0.17E-1 0.85E-2 0.83E-2 0.14E-3 0.73E-4   0.05   0.002 0.66E-2 0.34E-2 0.33E-2 0.29E-4 0.87E-5    0.05  0.05 0.16E-1 0.83E-2 0.82E-2 0.82E-2 0.62E-4  Table 2. Errors U-u at the points where x=0.3,0.6,0.9 at t ime t=1.0  for the six experiments of Problem 2.  Increments on G1  Errors on grids  h l  x G1 G3 G1,G3  0.1 0.1 0.3 0.16E-2 0.45E-3 0.77E-4    0.6 0.47E-2 0.36E-3 -0.11E-2    0.9 -0.54E-1 -0.27E-1 -0.18E-1  0.1 0.05 0.3 0.46E-3 0.11E-3 0.30E-6    0.6 0.41E-3 0.11E-3 0.15E-4    0.9 -0.22E-1 -0.83E-4 0.73E-2  0.1 0.002 0.3 0.73E-4 0.18E-4 0.30E-7    0.6 0.73E-4 0.18E-4 0.60E-7    0.9 -0.36E-2 0.18E-4 0.12E-2  0.05 0.1 0.3 0.15E-2 0.45E-3 0.10E-3    0.6 0.54E-2 0.35E-3 -0.13E-2    0.9 -0.19E+0 -0.16E-1 0.43E-1  0.05 0.05 0.3 0.45E-3 0.11E-3 0.17E-5    0.6 0.36E-3 0.11E-3 0.30E-4    0.9 -0.27E-1 0.30E-2 0.13E-1  0.05 0.02 0.3 0.73E-4 0.18E-4 -0.64E-6    0.6 0.73E-4 0.18E-4 -0.77E-6    0.9 0.64E-3 0.14E-3 -0.33E-4    (12)  Figure 1: Numerical solutions of Problem 3 with h =10, 31=l  on grid G1.     .......... Reference solution on G1.   Method (23) on G1.                Time-only extrapolation of (23).               Space-time extrapolation of (23). 

A difference scheme with high accuracy in time for fourth order parabolic equations&quot;,

A family of methods for the wave equation in one- and two-space dimensions&quot;,

A family of numerical methods for diffusion and reaction-diffusion equations&quot;,

A new computational procedure for ADI methods&quot;,

An explicit difference method for the wave equation with extended stability range&quot;,

Backward difference replacements of the space derivative in first order hyperbolic equations&quot;,

Convergence and order reduction of Runge-Kutta schemes applied to evolutionary problems in partial differential equations&quot;,

de Vries, &quot;Global extrapolation of a first order splitting method&quot;,

Explicit Runge-Kutta(-Nystrom) methods with reduced phase errors for computing oscillating solutions&quot;,

Fourth order difference methods for the initial-boundary value problem for hyperbolic equations&quot;,

L0 -stable splitting methods for the simple heat equation in two space dimensions with homogeneous boundary conditions&quot;,

Linear combinations of generalized Crank Nicolson schemes&quot;,

On the treatment of time-dependent boundary conditions in splitting methods for parabolic differential equations&quot;,

Some difference schemes for solving boundary problems&quot;,

The extrapolation of first order methods for parabolic partial differential equations.

The extrapolation of first order methods for parabolic partial differential equations.I.&quot;,

The extrapolation of first order methods for variable coefficient inhomogeneous parabolic partial differential equations&quot;,

Two new finite difference schemes for parabolic equations&quot;,

Global extrapolation procedures for linear partial differential equations

Global extrapolation procedures for linear partial differential equations

Abstract

Similar works

Full text

Available Versions

Brunel University Research Archive