A posteriori error estimates for leap-frog and cosine methods for second order evolution problems
Emmanuil H. Georgoulis, Omar Lakkis, Charalambos Makridakis, Juha M. Virtanen
aa r X i v : . [ m a t h . NA ] N ov A POSTERIORI ERROR ESTIMATESFOR LEAP-FROG AND COSINE METHODSFOR SECOND ORDER EVOLUTION PROBLEMS
EMMANUIL H. GEORGOULIS ∗ , OMAR LAKKIS † , CHARALAMBOS G. MAKRIDAKIS ‡ , AND
JUHA M. VIRTANEN § Abstract.
We consider second order explicit and implicit two-step time-discrete schemes forwave-type equations. We derive optimal order a posteriori estimates controlling the time discretiza-tion error. Our analysis, has been motivated by the need to provide a posteriori estimates for thepopular leap-frog method (also known as
Verlet ’s method in molecular dynamics literature); it isextended, however, to general cosine-type second order methods. The estimators are based on anovel reconstruction of the time-dependent component of the approximation. Numerical exper-iments confirm similarity of convergence rates of the proposed estimators and of the theoreticalconvergence rate of the true error.
1. Introduction.
This work is concerned with second order explicit and im-plicit two-step time-discrete schemes for wave-type equations. Our objective is toderive optimal order a posteriori estimates controlling the time-discretization error.To the best of our knowledge, error control for wave equations, discretized by pop-ular methods is limited so far to first order schemes [9, 13]. Despite the importanceof such wave-type problems, the lack of error control for time-discretizations usedextensively in applications is probably due to the two-step character of these meth-ods, and the associated technical issues. Our analysis, has been motivated by theneed to provide a posteriori estimates for the leap-frog method or, as often termed,
Verlet ’s method in the molecular dynamics literature. It extends, however to gen-eral cosine-type second order methods [5, 6].Adaptivity and a posteriori error control for parabolic problems have been devel-oped in, e.g., [12, 25, 21, 15, 19, 8, 10, 17]. In particular, as far as time discretizationis concerned, all implicit one-step methods can be treated within the framework de-veloped in [2, 20, 3, 4, 18]. Although some of these results apply (directly or afterappropriate modifications) to the wave equation also, when written as a first ordersystem and discretised by implicit Runge-Kutta or Galerkin schemes, this frame-work does not cover popular two-step implicit or explicit time-discretisation meth-ods. The recent results in [9, 13] cover only first order time discrete schemes; see also[1] for certain estimators to standard implicit time-stepping finite element approx-imations of the wave equation. For earlier works on adaptivity for wave equationsfrom various perspectives we refer, e.g., to [16, 7, 23, 24]. ∗ Department of Mathematics, University of Leicester, Leicester LE1 7RH, England UK, andSchool of Applied Mathematical and Physical Sciences, National Technical University of Athens,Athens 15780, Greece. Email:
[email protected] † Department of Mathematics, University of Sussex, Brighton BN1 9QH, England UK. ‡ Department of Mathematics, University of Sussex, Brighton BN1 9QH, England UK. Email:
[email protected] § Department of Mathematics, University of Leicester, Leicester LE1 7RH, England UK. Email: [email protected] odel problem and notation.. Let ( H, h· , ·i ) be a Hilbert space and A :[0 , T ] → D ( A ), positive definite, self-adjoint, linear operator on D ( A ), the domainof A , which is assumed to be dense in H . For time t ∈ (0 , T ], we consider the linearsecond order hyperbolic problem: find u : [0 , T ] → D ( A ), such that u ′′ ( t ) + A u ( t ) = f ( t ) for 0 < t ≤ T,u (0) = u ,u ′ (0) = v , (1.1)where f : [0 , T ] → H , u , v ∈ H . Leap-frog time-discrete schemes..
We shall be concerned with the popularleap-frog time-discrete scheme for (1.1). We consider a subdivision of the timeinterval (0 , T ] into disjoint subintervals ( t n , t n +1 ], n = 0 , . . . , N −
1, with t = 0and t N = T , and we define k n := t n +1 − t n , the time-step. For simplicity ofthe presentation, we shall assume that k n = k is constant, although this is not arestriction of the analysis below. Despite being two-step, the schemes consideredherein can be formulated for variable time steps also, with their consistency andstability properties then being influenced accordingly, cf. [22, 11]; the study of suchextensions is out of the scope of this work. We shall use the notation t n +1 / :=( t n +1 + t n ) / U n +1 ∈ D ( A ) of the exact values u n +1 := u ( t n +1 ), suchthat: ∂ U n +1 + A U n = f n , n = 1 , . . . , N − , (1.2)where f n := f ( t n ) ∈ H , ∂ U n +1 := ∂U n +1 − ∂U n k = U n +1 − U n + U n − k , (1.3)with ∂U n +1 := U n +1 − U n k , assuming knowledge of U and U . We set U := u and we define U by ∂U − v k + 12 A U = 12 f , (1.4)where f := f (0) and ∂U := ( U − U ) /k . This is a widely used and remarkablemethod in many ways: it is the only two-step explicit scheme for second order prob-lems which is second order accurate, it has important conservation and geometricproperties, as it is symplectic, and it is very natural and simple to formulate andimplement. We refer to the review article [14] for a thorough discussion. Explicitschemes, such as (1.2) are suited for the discretization of wave-type partial differ-ential equations, since their implementation requires a mild CFL-type conditionof the form k/h ≤ C (in contrast to parabolic problems,) where h stands for thespace-discretization parameter. osine methods.. The leap-frog scheme is a member of a general class oftwo-step methods for second order evolution problems, which are based on the ap-proximation of cosine and are used extensively in practical computations. In atwo-step cosine method, for n = 1 , . . . , N − , we seek approximations U n +1 suchthat ∂ U n +1 + (cid:2) q A U n +1 − p A U n + q A U n − (cid:3) = (cid:2) q f n +1 − p f n + q f n − (cid:3) , (1.5)where we assume that p = q − for second order accuracy; we refer to [5, 6] for adetailed discussion and analysis of general multi-step cosine schemes. In this case,the rational function r ( x ) = (1 + p x ) / (1 + q x ) is a second order approximationto the cosine, in the sense that for | x | sufficiently small, | r ( x ) − cos( x ) | ≤ C x . (1.6)When q = 0 the above methods are explicit, and the condition p = q − impliesthat the only explicit second order member of this family is the leap-frog method(1.2).In this work, we derive a posteriori error bounds in the L ∞ ()-in-time-energy-in-space norm of the error. The derived bounds are of optimal order, i.e., of the sameorder as the error (which is known to be second order) for the class of the schemesconsidered [5, 6]. This is verified by the numerical experiments presented herein.Our approach is based on the following ingredients: first, we rewrite the schemes asone-step system on staggered time grids. In turn, this can be seen as a second orderperturbation of the staggered midpoint method. Further, by introducing appropriateinterpolants, we arrive to a form which can be viewed as perturbation of (1.1) writtenas a first-order system. Finally, we employ an adaptation of the time reconstructionfrom [2], yielding the desired a posteriori error estimates. An interesting observationis that our estimates hold without any additional time-step assumption, which atthe fully discrete level would correspond to a CFL-type restriction. Thus, in aposteriori analysis, standard stability considerations of time discretisation schemesmight influence the behaviour of the estimator, but are not explicitly required; thepossible instability is sufficiently reflected by the behaviour of the estimator, seeSection 4. Although not done here, by employing space reconstruction techniquesit would be possible to derive error estimates for fully discrete schemes in variousnorms, using ideas from [19, 17, 13].The remaining of this work is organized as follows. In § § § § . Reformulation of the methods. It will be useful for the analysis to re-formulate the methods as a system in two staggered grids.
Starting with the leap-frog method, we introduce the auxiliaryvariable V n +1 / := ∂U n +1 , (2.1)for n = 0 , , . . . , N −
1, and we set V − / := 2 v − V / . (Note that, then, v =( V − / + V / ) / U − := U − kV − / and we observe that wehave v = U − U − k . Further, we introduce the notation ∂V n +1 / := V n +1 / − V n − / k , n = 0 , , . . . , N − , (2.2)noting that the identity ∂V / = 2( ∂U − v ) /k also holds.We can now write the method (1.2) as a system in the staggered form consideredin [14]: ∂U n +1 − V n +1 / = 0 ,∂V n +1 / + A U n = f n , (2.3)for n = 0 , , . . . , N − U : [ − k, T ] → D ( A ) to be the piecewise linear interpolant ofthe sequence { U n } Nn = − , at the points { t n } Nn = − , with t − := − k . In addition, let V : [ − k/ , t N − / ] → D ( A ) be the piecewise linear interpolant of the sequence { V n +1 / } N − n = − , at the points { t n +1 / } N − n = − . Using the notation U n +1 / := U ( t n +1 / ) ,V n := V ( t n ) , n = 0 , . . . , N − . (2.4)we, then, have U n +1 / = 12 ( U n +1 + U n ) , V n = 12 ( V n +1 / + V n − / ) , (2.5)for n = 0 , , . . . , N − ∂U n +1 −
12 ( V n +1 + V n ) = −
14 ( V n +3 / − V n +1 / + V n − / ) ,∂V n +1 / + 12 A ( U n +1 / + U n − / ) = f n + 14 A ( U n +1 − U n + U n − ) , (2.6)for n = 0 , . . . , N −
1. Upon defining the piecewise constant residuals R U ( t ) | ( t n − / ,t n +1 / ] ≡ R nU := 14 A ( U n +1 − U n + U n − ) , V ( t ) | ( t n ,t n +1 ] ≡ R n +1 / V := −
14 ( V n +3 / − V n +1 / + V n − / ) , it is easy to check that, given that the leap-frog method is second order (in both U n and V n +1 / ), we have R nU = O ( k ) and R n +1 / V = O ( k ) . Hence, (2.6) can be viewedas a second order perturbation of the staggered mid-point method for (1.1) writtenas first order system u ′ − v = 0 ,v ′ + A u = f . (2.7)In what follows, it will be useful to rewrite (2.6) as a perturbation of the continuoussystem (2.7). To this end, we introduce two time interpolants onto piecewise linearfunctions defined on the staggered grids: we define U : [0 , T ] → D ( A ) to be thepiecewise linear interpolant of the sequence { U n +1 / } N − n = − and V : [0 , t N − ] → D ( A )to be the piecewise linear interpolant of the sequence { V n } N − n =0 . Then, (2.6) can bewritten as U ′ − I V = R V ,V ′ + A ˜ I U = ˜ I f + R U , (2.8)where we define the interpolators˜ I : piecewise constant midpoint interpolator on { ( t n − / , t n +1 / ] } N − n =0 ,I piecewise constant midpoint interpolator on { ( t n − , t n ] } N − n =1 . (2.9)This formulation will be the starting point of our analysis in the next section. We shall see that cosine methods(1.5) can be reformulated in a similar way as a staggered system. As in the leap-frog case we introduce the auxiliary variable V n +1 / := ∂U n +1 , (2.10)and we let ∂V n +1 / := V n +1 / − V n − / k , n = 0 , , . . . , N − , (2.11)Then the methods (1.5) can be rewritten in system form: ∂U n +1 − V n +1 / = 0 ,∂V n +1 / + (cid:2) q A U n +1 − p A U n + q A U n − (cid:3) = (cid:2) q f n +1 − p f n + q f n − (cid:3) , (2.12) or n = 0 , , . . . , N −
1. Using the same notation and conventions as in the leap-frogcase, we observe, respectively, (cid:2) q A U n +1 − p A U n + q A U n − (cid:3) = 12 A ( U n +1 / + U n − / ) − A ( U n +1 / + U n − / ) + (cid:2) q A U n +1 − p A U n + q A U n − (cid:3) = 12 A ( U n +1 / + U n − / ) − A ( U n +1 / + U n − / ) + A U n + (cid:2) q A U n +1 − q A U n + q A U n − (cid:3) = 12 A ( U n +1 / + U n − / ) − (cid:2) A U n +1 − A U n + A U n − (cid:3) + (cid:2) q A U n +1 − q A U n + q A U n − (cid:3) = 12 A ( U n +1 / + U n − / ) − (1 − q )4 (cid:2) A U n +1 − A U n + A U n − (cid:3) , (2.13)where we used the fact that p = q − . Therefore, as before we conclude that ∂U n +1 −
12 ( V n +1 + V n ) = −
14 ( V n +3 / − V n +1 / + V n − / ) ,∂V n +1 / + 12 A ( U n +1 / + U n − / ) = ˜ f n + (1 − q )4 A ( U n +1 − U n + U n − ) , (2.14)where ˜ f n = (cid:2) q f n +1 − p f n + q f n − (cid:3) and n = 0 , . . . , N −
1. Let us now define R cos U ( t ) | ( t n − / ,t n +1 / ] ≡ R cos ,nU := (1 − q )4 A ( U n +1 − U n + U n − )+ q (cid:2) f n +1 − f n + f n − (cid:3) ,R cos V ( t ) | ( t n ,t n +1 ] ≡ R cos ,n +1 / V := −
14 ( V n +3 / − V n +1 / + V n − / ) . As in the leap-frog case, it is easy to check that given that the method is secondorder, we have R cos ,nU = O ( k ) and R cos ,n +1 / V = O ( k ) . Hence, (2.6) can be seen asa second order perturbation of the staggered mid-point method for (2.7). Further,still using the same notation for time interpolants as the leap-frog case, we obtain U ′ − I V = R cos V ,V ′ + A ˜ I U = ˜ I ˜ f + R cos U . (2.15)where ˜ I ˜ f | ( t n − / ,t n +1 / ] = ˜ f n . It is interesting to compare (2.15) to (2.8).
We briefly discuss an alternativeformulation of cosine methods. This time we let V n +1 / := (cid:0) I + k q A (cid:1) ∂U n +1 , (2.16)and, as before, ∂V n +1 / := V n +1 / − V n − / k , n = 0 , , . . . , N − . (2.17) hen, using again p = q − , we rewrite the methods (1.5) as ∂U n +1 − V n +1 / = − k q A ∂U n +1 ,∂V n +1 / + A U n = (cid:2) q f n +1 − p f n + q f n − (cid:3) , (2.18)for n = 0 , , . . . , N −
1. Using the same notation and conventions as in the leap-frogcase, we finally conclude ∂U n +1 −
12 ( V n +1 + V n ) = − k q A ∂U n +1 −
14 ( V n +3 / − V n +1 / + V n − / ) ,∂V n +1 / + 12 A ( U n +1 / + U n − / ) = ˜ f n + 14 A ( U n +1 − U n + U n − ) , (2.19)where ˜ f n = (cid:2) q f n +1 − p f n + q f n − (cid:3) and n = 0 , . . . , N −
1. Upon defining theperturbations as R cos , U ( t ) | ( t n − / ,t n +1 / ] ≡ R cos , ,nU := 14 A ( U n +1 − U n + U n − )+ q (cid:2) f n +1 − f n + f n − (cid:3) ,R cos , V ( t ) | ( t n ,t n +1 ] ≡ R cos , ,n +1 / V := − k q A ∂U n +1 −
14 ( V n +3 / − V n +1 / + V n − / ) , it is easy to check, again, that the method is second order ( R cos , ,nU = O ( k ) and R cos , ,n +1 / V = O ( k ) , ) and, thus, (2.6) can be interpreted as a second order pertur-bation of the staggered mid-point method for (2.7). Still, using the same notationas before, we have U ′ − I V = R cos , V ,V ′ + A ˜ I U = ˜ I ˜ f + R cos , U . (2.20)
3. A posteriori error bounds.
We have seen that all above schemes can bewritten in the form, V ′ + A ˜ I U = ˜ I ( f + ρ U ) ,U ′ − I V = I ρ V , (3.1)where ˜ I , I are defined in (2.9) and ˜ I ρ U equal to R U , R cos U , or R cos , U , for the leap-frog, the first and second cosine method formulations, respectively; similarly, I ρ V is equal to R V , R cos V , or R cos , V for each of the respective 3 formulations; cf., (2.8),(2.15) and (2.20). It is possible, in principle, to consider non-constant ρ U , ρ V oneach time-step; nevertheless the, easiest to implement, constant ones consideredhere suffice to deliver optimal estimator convergence rates, as will be highlighted inthe numerical experiments below. We continue by defining appropriate time reconstruc-tions, cf., [2]. To this end, on each interval ( t n − / , t n +1 / ], for n = 0 , . . . , N −
1, wedefine the reconstruction ˆ V of V byˆ V ( t ) := V n − / + Z tt n − / ( −A U + ˜ I f + ρ U ) d t, where ˜ I is a piecewise linear interpolant on the mesh { ( t n − / , t n +1 / ] } N − n =1 , suchthat ˜ I f ( t n ) = f n . We observe that ˆ V ( t n − / ) = V n − / , andˆ V ( t n +1 / ) = V n − / + k ( −A U ( t n ) + ˜ I f ( t n ) + ρ U ( t n )) = V n +1 / , sing the mid-point rule to evaluate the integral and the first equation in (2.6).Also, on each interval ( t n − , t n ], for n = 1 , . . . , N −
1, we define the reconstructionˆ U of U by ˆ U ( t ) := U n − + Z tt n − ( V + ρ V ) d t. Again, we observe that ˆ U ( t n − ) = U n − and thatˆ U ( t n ) = U n − + k ( V ( t n − / ) + ρ V ( t n − / )) = U n , using the mid-point rule. Notice that each of the above reconstructions is similar inspirit to the Crank-Nicolson reconstruction of [2]; however, we note that, althoughˆ U , ˆ V are both globally continuous functions, their derivatives jump alternatingly atthe nodes of the staggered grid. Setting ˆ e U := u − ˆ U and ˆ e V := u ′ − ˆ V ,we deduce ˆ e ′ V + A ˆ e U = R + R f ˆ e ′ U − ˆ e V = R , (3.2)with R := −A ( ˆ U − U ) − ρ U , R := ˆ V − V − ρ V , R f := f − ˜ I f. (3.3)For Φ = ( φ , φ ), Ψ = ( ψ , ψ ) ∈ D ( A ) × H , we define the bilinear form hh Φ , Ψ ii := hA / φ , A / ψ i + h φ , ψ i . It is evident that hh· , ·ii is an inner product on [ D ( A / ) × H ] . This is the standardenergy inner product and the induced norm, denoted by |k·|k , i.e., |k Φ |k = ( ||A / φ || + || φ || ) / , is the natural energy norm for (1.1).Then, the a posteriori error estimates will follow by applying standard energyarguments to the error equation (3.2). More specifically, in view of (3.2), we have12 dd t |k (ˆ e U , ˆ e V ) |k = hh (ˆ e ′ U , ˆ e ′ V ) , (ˆ e U , ˆ e V ) ii = hA ˆ e ′ U , ˆ e U i + h ˆ e ′ V , ˆ e V i = hA ˆ e V , ˆ e U i + hAR , ˆ e U i − hA ˆ e U , ˆ e V i + hR , ˆ e V i + hR f , ˆ e V i = hAR , ˆ e U i + hR , ˆ e V i + hR f , ˆ e V i , using the self-adjointness of A . Hence, using the Cauchy-Schwarz inequality, wearrive to 12 dd t |k (ˆ e U , ˆ e V ) |k ≤ |k ( R , R + R f ) |k|k (ˆ e U , ˆ e V ) |k . Integrating between 0 and τ , with 0 ≤ τ ≤ t N such that |k (ˆ e U , ˆ e V )( τ ) |k = sup t ∈ [0 ,t N ] |k (ˆ e U , ˆ e V )( t ) |k , e arrive to12 |k (ˆ e U , ˆ e V )( τ ) |k ≤ |k (ˆ e U , ˆ e V )(0) |k + |k (ˆ e U , ˆ e V )( τ ) |k Z τ |k ( R , R + R f ) |k d t, which implies |k (ˆ e U , ˆ e V )( τ ) |k ≤ |k (ˆ e U , ˆ e V )(0) |k + 4 (cid:16) Z τ |k ( R , R + R f ) |k d t (cid:17) . This already gives the following a posteriori bound.
Theorem 3.1.
Let u be the solution of (1.1) , ˆ e U := u − ˆ U and ˆ e V := u ′ − ˆ V .
Then, the following a posteriori error estimate holds sup t ∈ [0 ,t N ] |k (ˆ e U , ˆ e V )( t ) |k ≤ |k (ˆ e U , ˆ e V )(0) |k + 4 (cid:16) Z t N |k ( R , R + R f ) |k d t (cid:17) , where R , R , and R f are defined in (3.3) . An immediate Corollary from Theorem3.1 is an posteriori bound for the error sup [0 ,t N ] |k ( u − U, u ′ − V ) |k which can betrivially deduced through a triangle inequality. Remark 3.2.
Notice that due to the two-mesh stagerring, the computation ofthe ‘last’ V used in the above estimate, V N − / , requires the computation of U N +1 . This can be obtained by advancing one more time step in the computation beforeestimating.
4. Numerical Experiments.4.1. Fully discrete formulation..
Although the focus of the present work isin time-discretization, we shall introduce a fully discrete version of the time-steppingschemes for the numerical experiments bellow. To this end, we consider Let Ω ⊂ R d , d = 2 , ∂ Ω. We consider the initial-boundary valueproblem for the wave equation: find u ∈ L ∞ (0 , T ; H (Ω)) such that, u tt + A u = f in Ω × (0 , T ] , (4.1) u = u in Ω × { } (4.2) u t = v in Ω × { } (4.3) u = g on ∂ Ω × (0 , T ] , (4.4)where, for simplicity, we take A = − c ∆ , c = 0 , and g ∈ H / ( ∂ Ω). Further, for each n , we consider the standard, conforming finite element space S ph ⊂ H (Ω), based ona quasiuniform triangulation of Ω consisting of finite elements of polynomial degree p , with h denoting the largest element diameter. Focusing on time-discretizationissues, we shall use the same spatial discretization for all time-steps. The respectivediscrete spatial operator is denoted by A ph . The fully-discrete leap-frog method isthen defined as follows: for each n = 2 , . . . , N , find U n +1 ∈ S ph , such that U n +1 = 2 U n − U n − + k ( ¯ f n − A ph U n ) , (4.5)and U ∈ S ph such that U = U + kV + k A ph U − ¯ f ); (4.6) ere ¯ f n ( · ) := Π f ( · , t n ) for each t n , where Π : L (Ω) → S ph denotes a suitableinterpolation/projection operator onto the finite element space S ph of the sourcefunction f . We also set U := Π u and V := Π v . Note that V can be calculatedas above through U , or can be computed as follows: find V n +1 / ∈ S ph such that V n +1 / = V n − / + k ( ¯ f n − A ph U n ); (4.7)(4.7) can be used to overcome the difficulty of evaluating estimators defined onstaggered time mesh and depending on the term V n +3 / .To assess the time-error estimator, we replace A by its approximation A ph inthe a posteriori estimators discussed above, and in the (( · , · ))-inner product and ||| · ||| -norm. For brevity, we introduce the notation: e R := (ˆ e U , ˆ e V ) and e L := e R + ( U − ˆ U , V − ˆ V ) = ( u − U, u t − V ) . The objective is to study the performance of the a posteriori estimator η := ||| e R (0) ||| + 4 (cid:18)Z T ||| ( R , R + R f ) ||| (cid:19) ! / , (4.8)from Theorem 3.1. For t ∈ [0 ,
1] and Ω := (0 , , we consider the modelproblem (4.1) - (4.3) in the following set up: A := − c ∆ and f = 0, for c > u of problem (4.1) - (4.3) is given by: u ( t, x ) := ∞ X j,k =1 sin( kπx ) sin( jπx ) ( α k,j cos( ξ k,j πt ) + β k,j sin( ξ k,j πt )) (4.9)where α k,j > β k,j > ξ k,j := c p k + j . To illustrate the estimator’sbehaviour, we have chosen the following sets of parameters as numerical examples: ( c := 1 . ,α , = β , = 15 . , β k,j = α k,j = 0 for k, j = 1 , (4.10) ( c := 1 . ,α , = β , = 1 . , β k,j = α k,j = 0 for k, j = 3 , (4.11) ( c := 5 . ,α , = β , = 15 . , β k,j = α k,j = 0 for k, j = 1 , (4.12)The solutions of (4.10) - (4.12) are all smooth, but (4.12) oscillates much fastertemporally, while (4.11) has greater space-dependence of the error. In the numericalexperiments, the C++ library FEniCS/dolfin 1.2.0 and PETSc/SuperLU were usedfor the finite element formulation and the linear algebra implementation. For eachof the examples, we compute the solution of (4.5) using finite element spaces ofpolynomial degree p = 2, and time step size k = Ch r / ( p + 1) , r = 1 ,
2, for someconstant
C >
0, with h > S ph . The sequences of meshsizes considered (with respective olouring in the figures below) are: h = 1 / / / (4 √
2) (yellow),1 / /
10 (purple), for Examples (4.10) and (4.12), and h = 1 / (4 √
2) (cyan),1 / /
10 (yellow), 1 /
12 (red), 1 /
14 (purple) for Example (4.11). Note thatthe CFL-condition required by the leap-frog method, is satisfied by the restrictionon the time step k ≤ C h ( p + 1) , (4.13)for sufficiently small C > η , and the errors e R and e L , as well as the effectivity indexover time on a sequence of uniformly refined meshes with mesh sizes given as pereach example. We also monitor the energy of the reconstructed solution: E reconstruction := ( ˆ U , ˆ V ) . (4.14)We define experimental order of convergence (EOC) of a given sequence of pos-itive quantities a ( i ) defined on a sequence of meshes of size h ( i ) byEOC( a, i ) = log( a ( i + 1) /a ( i ))log( h ( i + 1) /h ( i )) , (4.15)the inverse effectivity index ,IEI( k e k L ∞ (0 ,t m ; |||·||| ) , η ) = k e k L ∞ (0 ,t m ; |||·||| ) η . (4.16)The IEIhas the same information as the (standard) effectivity index and has theadvantage of relating directly to the inequality appearing in Theorem 3.1. Theresults of numerical experiments on uniform meshes, depicted in Figures 5.1 and 5.2,indicate that the error estimators are reliable and also efficient provided the timesteps are kept sufficiently small. In the last experiment, Figure 5.3, the behaviourof the estimator is displayed, for the case when the CFL condition (4.13) is violated;the estimator remains reliable in this case also: the instability is reflected in thebehaviour of estimator, cf., Figure 5.3. Indeed, the inverse effectivity idex IEI ofthe estimator oscillates around a constant value for all times in the course of thenumerically unstable behaviour.
5. Concluding remarks.
An a posteriori error bound for error measured inL ∞ - norm in time and energy norm in space for leap-frog and cosine-type time semi-discretizations for linear second order evolution problems was presented and studiednumerically. The estimator was found to be reliable, with the same convergencerate as the theoretical convergence rate of the error. In a fully discrete setting, thisestimator corresponds to the control of time discretization error. The estimatorswere also found to be sharp on uniform meshes provided that the time steps and,thus, the time-dependent part of the error, is kept sufficiently small. Investigationinto the suitability of the proposed estimators within an adaptive algorithm remainsa future challenge. Slimane Adjerid , A posteriori finite element error estimation for second-order hyperbolicproblems , Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 4699–4719.[2]
Georgios Akrivis, Charalambos Makridakis, and Ricardo H. Nochetto , A posteriorierror estimates for the Crank-Nicolson method for parabolic equations , Math. Comp., 75(2006), pp. 511–531 (electronic).[3] ,
Optimal order a posteriori error estimates for a class of Runge-Kutta and Galerkinmethods , Numer. Math., 114 (2009), pp. 133–160.[4] ,
Galerkin and Runge-Kutta methods: unified formulation, a posteriori error estimatesand nodal superconvergence , Numer. Math., 118 (2011), pp. 429–456.[5]
Garth A. Baker, Vassilios A. Dougalis, and Steven M. Serbin , High order accu-rate two-step approximations for hyperbolic equations , RAIRO Anal. Num´er., 13 (1979),pp. 201–226.[6] ,
An approximation theorem for second-order evolution equations , Numer. Math., 35(1980), pp. 127–142.[7]
Wolfgang Bangerth and Rolf Rannacher , Adaptive finite element techniques for theacoustic wave equation , J. Comput. Acoust., 9 (2001), pp. 575–591.[8]
A. Bergam, C. Bernardi, and Z. Mghazli , A posteriori analysis of the finite elementdiscretization of some parabolic equations , Math. Comp., 74 (2005), pp. 1117–1138 (elec-tronic).[9]
Christine Bernardi and Endre S¨uli , Time and space adaptivity for the second-order waveequation , Math. Models Methods Appl. Sci., 15 (2005), pp. 199–225.[10]
Christine Bernardi and R¨udiger Verf¨urth , A posteriori error analysis of the fully dis-cretized time-dependent Stokes equations , M2AN Math. Model. Numer. Anal., 38 (2004),pp. 437–455.[11]
M. P. Calvo and J. M. Sanz-Serna , The development of variable-step symplectic integra-tors, with application to the two-body problem , SIAM J. Sci. Comput., 14 (1993), pp. 936–952.[12]
Kenneth Eriksson and Claes Johnson , Adaptive finite element methods for parabolicproblems. II. Optimal error estimates in L ∞ L and L ∞ L ∞ , SIAM J. Numer. Anal., 32(1995), pp. 706–740.[13] Emmanuil H. Georgoulis, Omar Lakkis, and Charalambos Makridakis , A posteriori L ∞ ( L ) -error bounds for finite element approximations to the wave equation , IMA J.Numer. Anal., 33 (2013), pp. 1245–1264.[14] Ernst Hairer, Christian Lubich, and Gerhard Wanner , Geometric numerical integra-tion illustrated by the St¨ormer-Verlet method , Acta Numer., 12 (2003), pp. 399–450.[15]
Paul Houston and Endre S¨uli , Adaptive Lagrange-Galerkin methods for unsteadyconvection-diffusion problems , Math. Comp., 70 (2001), pp. 77–106.[16]
Claes Johnson , Discontinuous Galerkin finite element methods for second order hyperbolicproblems , Comput. Methods Appl. Mech. Engrg., 107 (1993), pp. 117–129.[17]
Omar Lakkis and Charalambos Makridakis , Elliptic reconstruction and a posteriori errorestimates for fully discrete linear parabolic problems , Math. Comp., 75 (2006), pp. 1627–1658 (electronic).[18]
Christian Lubich and Charalambos Makridakis , Interior a posteriori error estimates fortime discrete approximations of parabolic problems , Numer. Math., 124 (2013), pp. 541–557.[19]
Charalambos Makridakis and Ricardo H. Nochetto , Elliptic reconstruction and a pos-teriori error estimates for parabolic problems , SIAM J. Numer. Anal., 41 (2003), pp. 1585–1594 (electronic).[20] ,
A posteriori error analysis for higher order dissipative methods for evolution problems ,Numer. Math., 104 (2006), pp. 489–514.[21]
Marco Picasso , Adaptive finite elements for a linear parabolic problem , Comput. MethodsAppl. Mech. Engrg., 167 (1998), pp. 223–237.[22]
Robert D. Skeel , Variable step size destabilizes the St¨ormer/leapfrog/Verlet method , BIT,33 (1993), pp. 172–175.[23]
E. S¨uli , A posteriori error analysis and global error control for adaptive finite volume approx- mations of hyperbolic problems , in Numerical analysis 1995 (Dundee, 1995), vol. 344 ofPitman Res. Notes Math. Ser., Longman, Harlow, 1996, pp. 169–190.[24] Endre S¨uli , A posteriori error analysis and adaptivity for finite element approximations ofhyperbolic problems , in An introduction to recent developments in theory and numericsfor conservation laws (Freiburg/Littenweiler, 1997), vol. 5 of Lect. Notes Comput. Sci.Eng., Springer, Berlin, 1999, pp. 123–194.[25]
R. Verf¨urth , A posteriori error estimates for finite element discretizations of the heat equa-tion , Calcolo, 40 (2003), pp. 195–212. 13 ig. 5.1 . Examples (4.10) and (4.11). η )01234560 0.25 0.5 0.75 1EOC(||e R || L ∞ (0,T;||| ⋅ |||) )01234560 0.25 0.5 0.75 1EOC(||e L || L ∞ (0,T;||| ⋅ |||) ) 01000020000300004000050000600000 0.25 0.5 0.75 1E reconstruction R || L ∞ (0,T;||| ⋅ |||) , η )10 -7 -6 -5 -4 -3 -2 -1 η -7 -6 -5 -4 -3 -2 -1 R || L ∞ (0,T;||| ⋅ |||) -7 -6 -5 -4 -3 -2 -1 L || L ∞ (0,T;||| ⋅ |||) (a) Example (4.10). Errors, estimator and IEI are plotted on thetop row and EOCs and energy of the reconstructed solution onthe bottom row over time ( x -axis). Results are computed on thesequence of uniform meshes with mesh size h , fixed time step k =0 . h/ ( p +1) and p = 2. The IEI behaviour indicates that the erroris well estimated by the estimator and the convergence rate of theestimator remains near to EOC ≈
2, i.e., to that of the errors e L and e R . η )01234560 0.25 0.5 0.75 1EOC(||e R || L ∞ (0,T;||| ⋅ |||) )01234560 0.25 0.5 0.75 1EOC(||e L || L ∞ (0,T;||| ⋅ |||) ) 01000020000300004000050000600000 0.25 0.5 0.75 1E reconstruction R || L ∞ (0,T;||| ⋅ |||) , η )10 -7 -6 -5 -4 -3 -2 -1 η -7 -6 -5 -4 -3 -2 -1 R || L ∞ (0,T;||| ⋅ |||) -7 -6 -5 -4 -3 -2 -1 L || L ∞ (0,T;||| ⋅ |||) (b) Example (4.11). Errors, estimator and IEI are plotted on thetop row and EOCs and energy of the reconstructed solution onthe bottom row over time ( x -axis). Results are computed on thesequence of uniform meshes with mesh size h , fixed time step k =0 . h/ ( p +1) and p = 2. The IEI behaviour indicates that the erroris well estimated by the estimator and the convergence rate of theestimator remains near to EOC ≈
2, i.e., to that of the errors e L and e R . 14 ig. 5.2 . Example (4.12). η )01234560 0.25 0.5 0.75 1EOC(||e R || L ∞ (0,T;||| ⋅ |||) )01234560 0.25 0.5 0.75 1EOC(||e L || L ∞ (0,T;||| ⋅ |||) ) 01000020000300004000050000600000 0.25 0.5 0.75 1E reconstruction R || L ∞ (0,T;||| ⋅ |||) , η )10 -7 -6 -5 -4 -3 -2 -1 η -7 -6 -5 -4 -3 -2 -1 R || L ∞ (0,T;||| ⋅ |||) -7 -6 -5 -4 -3 -2 -1 L || L ∞ (0,T;||| ⋅ |||) (a) Example (4.12). Errors, estimator and IEI are depicted onthe top row and EOCs and energy of the reconstructed solutionon the bottom row over time ( x -axis). Results are computed onthe sequence of uniform meshes with mesh size h , time-step k =0 . h/ ( p +1) and p = 2. The IEI behaviour indicates that the erroris overestimated by the estimator. The convergence rate of theestimator is slightly below 2, probably due to the somewhat coarsetime-step for asymptotic convergence, cf., subfigure (b) below. η )01234560 0.25 0.5 0.75 1EOC(||e R || L ∞ (0,T;||| ⋅ |||) )01234560 0.25 0.5 0.75 1EOC(||e L || L ∞ (0,T;||| ⋅ |||) ) 01000020000300004000050000600000 0.25 0.5 0.75 1E reconstruction R || L ∞ (0,T;||| ⋅ |||) , η )10 -7 -6 -5 -4 -3 -2 -1 η -7 -6 -5 -4 -3 -2 -1 R || L ∞ (0,T;||| ⋅ |||) -7 -6 -5 -4 -3 -2 -1 L || L ∞ (0,T;||| ⋅ |||) (b) Example (4.12). Errors, estimator and IEI are depicted onthe top row and EOCs and energy of the reconstructed solutionon the bottom row over time ( x -axis). Results are computed onthe sequence of uniform meshes with mesh size h , time-step k =0 . h / ( p +1) , and p = 2. The effect of using slightly smaller time-step for fine h (resulting from the h -term) is evident compared tosubfigure (a) above, in that the EOC ≥ ig. 5.3 . Example (4.10), violation of the CFL condition. η )01234560 0.25 0.5 0.75 1EOC(||e R || L ∞ (0,T;||| ⋅ |||) )01234560 0.25 0.5 0.75 1EOC(||e L || L ∞ (0,T;||| ⋅ |||) ) 01000020000300004000050000600000 0.25 0.5 0.75 1E reconstruction R || L ∞ (0,T;||| ⋅ |||) , η )10 -7 -6 -5 -4 -3 -2 -1 η -7 -6 -5 -4 -3 -2 -1 R || L ∞ (0,T;||| ⋅ |||) -7 -6 -5 -4 -3 -2 -1 L || L ∞ (0,T;||| ⋅ |||) (a) Example (4.10). Errors, estimator and IEI are depicted onthe top row and EOCs and energy of the reconstructed solutionon the bottom row over time ( x -axis). Results are computed onthe sequence of uniform meshes with mesh size h and time step k = 2 . h/ ( p + 1) and pp