# Closing the gap between trigonometric integrators and splitting methods for highly oscillatory differential equations

Closing the gap between trigonometric integrators and splitting methods for highly oscillatory... Abstract An error analysis of symmetric trigonometric integrators applied to highly oscillatory linear second-order differential equations is given. Second-order convergence is shown uniformly in the high frequencies under a finite-energy condition on the exact solution. The main novelty is the concept of the proof itself, which is based on the interpretation of trigonometric integrators as splitting methods for averaged differential equations. This allows one to combine techniques for splitting methods with those for trigonometric integrators. For the bound of the global error, cancellations in the error accumulation have to be studied carefully. 1. Introduction Many ordinary differential equations in science can be written as   u′=F1(u)+F2(u) (1.1) with linear or nonlinear vector fields $$F_1$$ and $$F_2$$, and with the property that solving the two sub-problems   v′=F1(v)andw′=F2(w) exactly or approximately is numerically much more efficient than solving (1.1) directly. Typical examples include many-body problems and other Hamiltonian systems with a certain structure, but also semidiscretizations of time-dependent partial differential equations such as, e.g., linear or nonlinear Schrödinger equations, linear or nonlinear wave equations, or Maxwell’s equations. In such a situation, splitting methods are easy to implement and very efficient, because these integrators compute a numerical approximation of the full problem by a composition of the flows of the two sub-problems. In many cases, the numerical solution inherits certain geometric properties of the exact solution, such as norm conservation or symplecticity of the flow. The order conditions and the accuracy of splitting methods have been analysed in many papers such as, e.g., Strang (1968), Jahnke & Lubich (2000), Lubich (2008), Hansen & Ostermann (2009), Thalhammer et al. (2009), Koch & Lubich (2011), Holden et al. (2013), Einkemmer & Ostermann (2014), Einkemmer & Ostermann (2015), Faou et al. (2015), Hochbruck et al. (2015) and references therein. Overviews have been given, e.g., in McLachlan & Quispel (2002), Hairer et al. (2006), Hundsdorfer & Verwer (2007), Holden et al. (2010) and Blanes & Casas (2016). The classical procedure is to estimate the local error by proving bounds for certain iterated commutators between the two vector fields. The error bound for the global error then follows by standard arguments like Gronwall inequalities or Lady Windermere’s fan. Unfortunately, the accuracy typically suffers when splitting methods are applied to problems with highly oscillatory solutions. In order to obtain a reasonable approximation, the step size must be considerably smaller than the inverse of the highest frequency, which reduces the efficiency considerably. We remark, however, that certain invariants of the exact flow are often conserved over long times even in the presence of oscillations (see e.g., Hairer et al., 2006; Faou, 2012; Cohen et al., 2015). In the present article, we study the discretization error on finite time intervals and do not consider the additional numerical resonances that show up on these longer time scales. A particular class of problems with highly oscillatory solutions takes the semilinear form   q″(t)=−Ω2q(t)+g(q(t)),t>0,q(0)=q0,q′(0)=q0′, (1.2) where $${\it{\Omega}}$$ is a symmetric positive semidefinite matrix of arbitrary large norm, and where $$g$$ is smooth and bounded. For such problems, trigonometric integrators have been constructed and analysed in GarcÍa-Archilla et al. (1999), Hochbruck & Lubich (1999), Hairer et al. (2006) and Grimm & Hochbruck (2006) under a finite-energy condition. These methods involve filter functions which are chosen in such a way that oscillatory parts of the local error do not sum up in the global error. To prove this in the error analysis is a delicate matter, and usually excludes to apply the technique of Lady Windermere’s fan in a standard way. An exception is semilinear wave equations with polynomial or analytic nonlinearities, see Gauckler (2015). Hence, it seems that splitting methods and trigonometric integrators are two different approaches, based on different ideas, having different properties, and requiring different techniques for their analysis. Nevertheless, there is a bridge between these two worlds: it is long known that symmetric trigonometric integrators applied to (1.2) can be interpreted as a Strang splitting method applied to an averaged version of the first-order formulation of (1.2). This interesting link raises a number of questions: Is it possible to gain a better understanding by considering trigonometric integrators as splitting methods for averaged equations? Does this relation allow to apply the techniques developed for one class of methods to the other one? And, most importantly, is it then possible to construct and analyze efficient numerical integrators for fully nonlinear problems? These questions are our motivation for the analysis below. In this article, we will restrict ourselves to the linear variant of (1.2), i.e., $$g(q) = Gq$$ with a matrix $$G$$ with a moderate norm $$\| G \| \ll \| {\it{\Omega}}\|$$. Here and in the following, $$\| \cdot \|$$ denotes the Euclidean vector norm or its induced matrix norm, respectively. The analysis in Lubich (2008) shows that for splitting methods the calculus of Lie derivatives allows, at least to some extent, to carry over techniques developed for linear differential equations to nonlinear ones. After reformulation as a first-order problem, we will derive the corresponding averaged equation and prove a bound for the difference between the solution of both equations; cf. Theorem 4.1. Then, we present an error analysis for the classical Strang splitting applied to the averaged equation, which yields a result very similar, but not equivalent to the one obtained in Grimm & Hochbruck (2006). We stress that the novelty of our analysis is not the error bound itself, but the fact that it is proven by techniques which, to the best of our knowledge, have so far not been considered in the context of trigonometric integrators. 2. Problem setting, assumptions and notation The situation we have in mind is that the ordinary differential equation   q″(t)=−Ω2q(t)+Gq(t),0≤t≤tend,q(0)=q0,q′(0)=q0′ (2.1) stems from a spatial discretization of a linear wave equation with finite elements, finite differences or spectral methods on a family of finer and finer meshes. In this situation, the matrix $$-{\it{\Omega}}^2$$ stems from a spatial discretization of the Laplacian or a more general differential operator, and $${\it{\Omega}}$$ is a symmetric, positive definite matrix (possibly after shifting $${\it{\Omega}} \rightarrow \sqrt{{\it{\Omega}}^2+I}$$ and $$G\rightarrow G+I$$). Then $$\left\lVert {{\it{\Omega}}} \right\rVert$$ becomes arbitrarily large if the spatial mesh width tends to zero, but $$\lVert {\it{\Omega}}^{-1} \rVert$$ remains uniformly bounded independently of the discretization. This motivates the following assumption. Assumption 2.1 Let $${\it{\Omega}} \in \mathcal{F}$$, where $$\mathcal{F}$$ is a family of symmetric, positive definite matrices such that there is a constant $$C_\text{inv}$$ with   ‖Ω−1‖≤Cinvfor all Ω∈F. (2.2) For the matrix $$G$$ we assume that its norm is bounded independently of $${\it{\Omega}}$$. Our aim is to prove error estimates which are uniform for all matrices in the family $$\mathcal{F}$$, which means that they are independent of $$\left\lVert {{\it{\Omega}}} \right\rVert$$ (i.e., independent of the spatial discretization). On the other hand, the constant $$C_\text{inv}$$ only depends on the coercivity constant of the corresponding differential operator. For the solution we rely on the following assumption. Assumption 2.2 Let the solution $$q:[0,t_{\text{end}}] \rightarrow \mathbb R^d$$ of (2.1) fulfill the finite-energy condition  ‖Ωq(t)‖2+‖q′(t)‖2≤K2,0≤t≤tend (2.3) with a constant $$K>0$$ on a finite time interval of length $$t_{\text{end}}$$. In fact, it can be shown that (2.3) is true on bounded time intervals if the initial data satisfy the bound $$\left\lVert {{\it{\Omega}} q(0)} \right\rVert^2 + \left\lVert {q'(0)} \right\rVert^2 \leq K_0^2$$ with a sufficiently small $$K_0 \leq K$$; cf. Lemma 4.2 below. In order to apply a splitting scheme, we formulate (2.1) as a first-order problem. We define the new variable   u=[qΩ−1q′],u0=[q0Ω−1q0′] which solves the differential equation   u′=Au+Bu,u(0)=u0 (2.4) with matrices   A=[0Ω−Ω0],B=[00Ω−1G0]. Since $$A$$ is skew symmetric, the exponential   exp⁡(tA)=[cos⁡(tΩ)sin⁡(tΩ)−sin⁡(tΩ)cos⁡(tΩ)] is unitary and thus   ‖etA‖=1,t∈R. (2.5) The finite-energy condition (2.3) is equivalent to   ‖Au(t)‖≤K,0≤t≤tend. (2.6) The solution of (2.4) could, in principle, be approximated with the classical Strang splitting, but in our setting this method only yields an acceptable accuracy if $$\tau \ll \|{\it{\Omega}}\|^{-1}$$. Such a severe step-size restriction is not acceptable in practice. We will show that much better approximations are obtained if the Strang splitting is applied to an averaged version of (2.4). The methods considered below involve analytic (mostly trigonometric) matrix functions of $$\tau{\it{\Omega}}$$. We assume that products of such matrix functions with a given vector can be computed efficiently, e.g., via fast Fourier transforms or (rational) Krylov subspace methods (see e.g., Grimm & Hochbruck, 2008). 3. Trigonometric integrators as splitting methods Symmetric one-step trigonometric integrators (Hairer et al., 2006, Section XIII.2.2) and Grimm & Hochbruck (2006) written in terms of the variables $$u_n=[q_n, {\it{\Omega}}^{-1}q'_n] \approx u(t_n)=[q(t_n), {\it{\Omega}}^{-1}q'(t_n)]$$ at $$t_n = n\tau$$ are given by   un+1 =eτAun+τ2[sin⁡(τΩ)Ω−1G~qncos⁡(τΩ)Ω−1G~qn+Ω−1G~qn+1]. (3.1) Here, $$\tau>0$$ denotes the step-size, and the matrix $$\widetilde{G}$$ is defined as   G~=ΨSGΦ,Φ=ϕ(τΩ),ΨS=ψS(τΩ) (3.2) with appropriate scalar filter functions $$\phi$$ and $$\psi_S$$. Note that $$\psi_S$$ was denoted by $$\psi_1$$ in Hairer et al. (2006) and Grimm & Hochbruck (2006). In these references, it was shown that under the assumptions of Section 2 and certain (sufficient) conditions on the filter functions $$\phi$$ and $$\psi_S$$, these methods yield a global accuracy of $$O(\tau^2)$$ in the positions $$q_n$$ and $$O(\tau)$$ in the velocities $$q_n'$$. The constants in the error bounds are independent of $$\| {\it{\Omega}} \|$$, in spite of the oscillatory nature of the solution. It is well known that symmetric trigonometric integrators can be interpreted as splitting schemes. With   B~=[00Ω−1G~0], (3.3) the method (3.1) can be written as   (I−τ2B~)un+1 =eτA(I+τ2B~)un. (3.4) Now we use that   eθB~=I+θB~for all θ∈R (3.5) due to the particular structure of $$\widetilde{B}$$. Substituting (3.5) into (3.4) gives   un+1=Sun,S=eτ2B~eτAeτ2B~. (3.6) Hence, the trigonometric integrator (3.1) is equivalent to the exponential Strang splitting method applied to the averaged equation   u~′=Au~+B~u~,u~(0)=u~0=u0 (3.7) with $$\widetilde{B}$$ defined in (3.3) and (3.2). This means that applying the Strang splitting to the averaged equation yields very good approximations for step sizes, where the Strang splitting applied to the original problem (2.4) fails if the norm of $${\it{\Omega}}$$ is large. We point out that not only the numerical solution $$u_n$$ depends on the step-size $$\tau$$, but also the exact solution $$\widetilde{u}$$ of the averaged equation via the filter functions $${\it{\Phi}}=\phi(\tau {\it{\Omega}})$$ and $${\it{\Psi}}_S=\psi_S(\tau{\it{\Omega}})$$. Note that for semilinear problems (1.2) the trigonometric integrators can still be interpreted as a Strang splitting applied to a (semilinear) averaged equation. The interpretation of trigonometric integrators as second-order Strang splitting methods is well known (see GarcÍa-Archilla et al., 1999 or Hairer et al., 2006). Recently, it has been shown in Blanes (2015) that also certain methods of higher order designed for oscillatory differential equations of the form (1.2) can be interpreted as higher order splitting methods. The methods considered in Blanes (2015) are exponentially fitted Runge–Kutta–Nyström methods and extended Runge–Kutta–Nyström methods (see also the monographs Ixaru & Vanden Berghe, 2004; Wu et al., 2013; Wu et al., 2015 and the review article Paternoster, 2012). We emphasize, however, that higher order methods cannot be expected to converge with order higher than two under the finite-energy condition (2.6). Instead, third or higher order error bounds are expected to require a stronger condition on the exact solution of the form $$\left\lVert A^q u(t) \right\rVert \leq K$$ with $$q\ge 2$$, which amounts in particular to a small-energy assumption if the matrix $${\it{\Omega}}$$ contains only high frequencies. Our aim is to prove second-order error bounds for the trigonometric integrator (3.1) on the basis of the interpretation as a splitting method (3.6) and under the finite-energy condition (2.6), using a carefully adapted Lady Windermere’s fan argument familiar from splitting integrators. In order to analyse the error of the splitting scheme, we first study the properties of the solution of the averaged problem (3.7). This allows us to bound the error which results from solving the averaged equation instead of (2.4) (Section 4). To analyse the error of the splitting scheme for the averaged problem, we first give a new representation of the local error. Unfortunately, this local error still contains a term which is not uniformly of third order in $$\tau$$ and requires a more careful investigation (Section 5). The error analysis given below relies on the following assumption on the filter functions. Assumption 3.1 The filter functions $$\chi = \phi$$ or $$\chi = \psi_S$$ are even analytic functions with the properties   χ(0) =1, (3.8a)  |χ(x)| ≤M1, (3.8b)  |xχ(x)| ≤M2, (3.8c)  |cot⁡(x2)xχ(x)| ≤M3 (3.8d) for certain constants $$M_j$$ uniformly for all $$x \in \mathbb R$$. Even functions $$\chi$$ guarantee that the scheme is symmetric. A popular example used in trigonometric integrators is $$\chi(x) = \mathop{\mathrm{sinc}}(x)$$. Note that (3.8a) and the condition that $$\chi$$ is even analytic imply   |x−2(1−χ(x))|≤M4. (3.8e) In the following $$C$$ denotes a generic constant which may have different values at different occurrences. Our main result is stated in the following theorem. Theorem 3.2 (Main result) Let Assumptions 2.1, 2.2 and 3.1 be fulfilled. Then the global error of (3.6) applied to (2.4) is bounded by   ‖un−u(tn)‖≤Cτ2,0≤tn=nτ≤tend with a constant $$C$$ that only depends on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_j$$, $$j=1,\ldots,4$$ and $$t_{\text{end}}$$ but not on $$\left\lVert {{\it{\Omega}}} \right\rVert$$. Proof. We combine the results from Theorem 4.1 and Theorem 5.3 below to obtain   ‖un−u(tn)‖≤‖un−u~(tn)‖+‖u~(tn)−u(tn)‖≤Cτ2,0≤tn=nτ≤tend, where the constant $$C$$ only depends on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_j$$, $$j=1,\ldots,4$$ and $$t_{\text{end}}$$. □ 4. Properties of the averaged equation In this section, we prove that the solutions of the original problem (2.4) and the averaged equation (3.7) differ only by $$O(\tau^2)$$. Theorem 4.1 Let the assumptions (3.8a) and (3.8b) and the finite-energy condition (2.3) be fulfilled. Let $$u$$ be the solution of (2.4) and $$\widetilde{u}$$ be the solution of the averaged system (3.7). Then it holds   ‖u(t)−u~(t)‖≤Cavτ2,0≤t≤tend, where $$C_\text{av}$$ only depends on $$C_\text{inv}$$, $$\|G\|$$, $$K$$, $$M_1$$, $$M_4$$ and $$t_{\text{end}}$$. Proof. The variation-of-constants formula yields   u(t) =etAu0+∫0te(t−s)ABu(s)ds, (4.1a)  u~(t) =etAu0+∫0te(t−s)AB~u~(s)ds. (4.1b) We define the block diagonal matrices   Ψ^=[ΨS00ΨS],Φ^=[Φ00Φ]. (4.2) Obviously, they satisfy $$A\widehat{{\it{\Psi}}}=\widehat{{\it{\Psi}}} A$$ and $$A\widehat{{\it{\Phi}}}=\widehat{{\it{\Phi}}} A$$. We also have $$\widetilde{B} = \widehat{{\it{\Psi}}} B \widehat{{\it{\Phi}}}$$ so that   u(t)−u~(t) =∫0te(t−s)A((I−Ψ^)Bu(s)+Ψ^B(I−Φ^)u(s)+Ψ^BΦ^(u(s)−u~(s)))ds =I1(t)+I2(t)+∫0te(t−s)AB~(u(s)−u~(s))ds with   I1(t) =∫0te(t−s)A(I−Ψ^)Bu(s)ds,I2(t)=∫0te(t−s)AΨ^B(I−Φ^)u(s)ds. For the first term integration by parts yields   I1(t) =−τ2e(t−s)A(τA)−2(I−Ψ^)ABu(s)|0t+τ2∫0te(t−s)A(τA)−2(I−Ψ^)ABu′(s)ds. With (2.2) and (2.6) it follows that   ‖u‖ =‖A−1Au‖≤CinvK, (4.3a)  ‖u′‖ ≤‖Au‖+‖Bu‖≤K+Cinv2‖G‖K. (4.3b) Using in addition (2.5), (3.8e) and   ‖AB‖=‖G‖, (4.4) this yields the estimate $$\lVert I_1(t) \rVert \leq C \tau^2$$ with a constant $$C$$ depending only on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_4$$ and $$t_{\text{end}}$$. For the second term, we insert the variation of constants formula once more and obtain   I2(t) =∫0te(t−s)AΨ^B(I−Φ^)u(s)ds =∫0te(t−s)AΨ^B(I−Φ^)esAu0ds+∫0t∫0se(t−s)AΨ^B(I−Φ^)e(s−σ)ABu(σ)dσds. Integration by parts yields   I2(t) =e(t−s)AΨ^B(I−Φ^)A−1esAu0|0t+∫0te(t−s)AAΨ^B(I−Φ^)A−1esAu0ds −∫0te(t−s)AΨ^B(I−Φ^)A−1e(s−σ)ABu(σ)|0sds +∫0t∫0se(t−s)AΨ^B(I−Φ^)A−1e(s−σ)ABu′(σ)dσds =τ2e(t−s)AΨ^B(I−Φ^)(τA)−2esAAu0|0t+τ2∫0te(t−s)AΨ^AB(I−Φ^)(τA)−2esAAu0ds −τ2∫0te(t−s)AΨ^B(I−Φ^)(τA)−2e(s−σ)AABu(σ)|0sds +τ2∫0t∫0se(t−s)AΨ^B(I−Φ^)(τA)−2e(s−σ)AABu′(σ)dσds. By (2.2), (2.5), (3.8), (4.4) and (4.3) we obtain $$\lVert I_2(t) \rVert \leq C \tau^2$$, where $$C$$ depends on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_1$$, $$M_4$$ and $$t_{\text{end}}$$. Hence, we have   ‖u(t)−u~(t)‖ =‖I1(t)+I2(t)+∫0te(t−s)AB~(u−u~)ds‖≤Cτ2+∫0tβ‖u(s)−u~(s)‖ds with   β=Cinv‖G‖M12≥‖B~‖. (4.5) An application of Gronwall’s Lemma proves the desired result. □ We finally show that the solution of the averaged problem (3.7) inherits the finite-energy condition (2.3) of the original problem. Lemma 4.2 Let the assumption (3.8b) and the finite-energy condition (2.3) be fulfilled. Then the solution $$\widetilde{u}$$ of the averaged system (3.7) satisfies the finite-energy condition   ‖Au~(t)‖≤K~,0≤t≤tend, (4.6) where $$\widetilde{K}$$ depends only on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_1$$ and $$t_{\text{end}}$$. Proof. Using (2.5) and the variation of constants formula (4.1b) yields   ‖Au~(t)‖ ≤‖Au0‖+∫0t‖AB~‖‖A−1‖‖Au~(s)‖ds. By definition (3.3) and (3.8b), it holds   ‖AB~‖=‖G~‖≤M12‖G‖ (4.7) and the statement thus follows from (2.2) and Gronwall’s Lemma. □ 5. Finite-time error analysis of the splitting scheme In this section we finally study the error of the splitting scheme. We use the notation   S=eτ2B~eτAeτ2B~ defined in (3.6) for the numerical flow over a time step $$\tau$$ and   T=eτ(A+B~) for the exact flow of the averaged equation (3.7). By (3.3) and (3.5) we have   max{‖eθB~‖,‖AeθB~A−1‖}≤1+θβ≤eθβ (5.1) for $$\theta\ge 0$$ with $$\beta$$ from (4.5). Together with (2.5), this provides the stability estimate   ‖Sn‖≤(1+τ2β)2n≤eβtend,0≤tn≤tend. (5.2) We next consider the local error. Lemma 5.1 (Local error) Assume that the filter functions are bounded as stated in (3.8b) and (3.8c). If the finite-energy condition (4.6) holds true, then the local error at time $$t_n=n\tau$$, $$0 \leq t_n \le t_{\text{end}}-\tau$$, of the splitting method (3.6) as an approximation to the averaged system (3.7) is given by   δn=(S−T)u~(tn)=δ^n+Dn, where   δ^n=12∫0τ∫0ξ∫0θe(ξ−σ)ALe(σ+τ−ξ)Au~(tn)dσdθdξ,L=A2B~+B~A2 (5.3) and $$\lVert D_n \rVert \leq C \tau^3$$. The constant $$C$$ depends on $$C_\text{inv}$$, $$\|G\|$$, $$\widetilde{K}$$, $$M_1$$, $$M_2$$ and $$t_{\text{end}}$$. Proof. We start with the following representation of the local error   δn =∫0τddξ(eξ2B~eξAeξ2B~e(τ−ξ)(A+B~))u~(tn)dξ =12∫0τeξ2B~([B~,eξA]eξ2B~+2eξA[A,eξ2B~])e(τ−ξ)(A+B~)u~(tn)dξ, (5.4) where $$\left[ A_1,A_2 \right]=A_1A_2 - A_2A_1$$ denotes the commutator of two matrices $$A_1$$ and $$A_2$$. The commutators in $$\delta_n$$ are written as integrals   [B~,eξA] =∫0ξddθ(e(ξ−θ)AB~eθA)dθ=−∫0ξe(ξ−θ)A[A,B~]eθAdθ and   [A,eξ2B~] =∫0ξddθ(eξ−θ2B~Aeθ2B~)dθ=12∫0ξeξ−θ2B~[A,B~]eθ2B~dθ. This yields   δn =12∫0τeξ2B~∫0ξ(−e(ξ−θ)A[A,B~]eθAeξ2B~+eξAeξ−θ2B~[A,B~]eθ2B~)e(τ−ξ)(A+B~)u~(tn)dθdξ. We add and subtract the term $${\it \,{e}}^{\xi A} [A,\widetilde{B}] {\it \,{e}}^{\frac{\xi}2 \widetilde{B}}$$ inside the brackets, and the local error can be split into $$\delta_n = \delta_n^{(1)} + D_n^{(1)}$$, where   δn(1) =12∫0τeξ2B~∫0ξ(eξA[A,B~]−e(ξ−θ)A[A,B~]eθA)eξ2B~e(τ−ξ)(A+B~)u~(tn)dθdξ, (5.5a)  Dn(1) =12∫0τeξ2B~∫0ξeξA(eξ−θ2B~[A,B~]eθ2B~−[A,B~]eξ2B~)e(τ−ξ)(A+B~)u~(tn)dθdξ. (5.5b) The term in brackets in (5.5a) can be written as   eξA[A,B~]−e(ξ−θ)A[A,B~]eθA =−∫0θddσ(e(ξ−σ)A[A,B~]eσA)dσ =∫0θe(ξ−σ)A[A,[A,B~]]eσAdσ. For the double commutator, we have   [A,[A,B~]] =L−R,L=A2B~+B~A2,R=2AB~A, and we split $$\delta_n^{(1)}$$ accordingly into $$\delta_n^{(1)} = \delta_n^{(2)}+D_n^{(2)}$$ with   δn(2) =12∫0τeξ2B~∫0ξ∫0θe(ξ−σ)ALeσAeξ2B~e(τ−ξ)(A+B~)u~(tn)dσdθdξ,Dn(2) =−12∫0τeξ2B~∫0ξ∫0θe(ξ−σ)AReσAeξ2B~e(τ−ξ)(A+B~)u~(tn)dσdθdξ. To extract the term $$\widehat \delta_n$$ defined in (5.3) from $$\delta_n^{(2)}$$, we use (3.5) and once more the variation-of-constants formula to obtain   δn(2) =δ^n+Dn(3)+Dn(4)+Dn(5), where   Dn(3) =12∫0τ∫0ξ∫0θe(ξ−σ)ALeσA∫0τ−ξe(τ−ξ−ν)AB~eν(A+B~)u~(tn)dνdσdθdξ,Dn(4) =14∫0τ∫0ξ∫0θe(ξ−σ)ALeσAξB~e(τ−ξ)(A+B~)u~(tn)dσdθdξ,Dn(5) =14∫0τξB~∫0ξ∫0θe(ξ−σ)ALeσAeξ2B~e(τ−ξ)(A+B~)u~(tn)dσdθdξ. In this way, we end up with the decomposition   δn=δ^n+Dn,Dn=Dn(1)+Dn(2)+Dn(3)+Dn(4)+Dn(5) of the local error. To prove the statement of the lemma, we have to estimate the terms $$D_n^{(i)}$$, $$i=1,\dots,5$$. The properties (2.2), (3.8b), (3.8c), (4.4) and (4.7) imply that   ‖RA−1‖ =2‖AB~‖≤2M12‖G‖, (5.6a)  ‖τLA−1‖ =‖τAAB~A−1+τB~A‖=‖(τAΨ^)(AB)Φ^A−1+Ψ^B(Φ^τA)‖≤2CinvM1M2‖G‖, (5.6b) where we use the notation (4.2). Writing   Dn(2) =−12∫0τeξ2B~∫0ξ∫0θe(ξ−σ)A(RA−1)eσA(Aeξ2B~A−1)Au~(tn+τ−ξ)dσdθdξ,Dn(3) =12τ∫0τ∫0ξ∫0θe(ξ−σ)A(τLA−1)eσA∫0τ−ξe(τ−ξ−ν)A(AB~)A−1Au~(tn+ν)dνdσdθdξ,Dn(4) =14τ∫0τξ∫0ξ∫0θe(ξ−σ)A(τLA−1)eσA(AB~)A−1Au~(tn+τ−ξ)dσdθdξ,Dn(5) =14τ∫0τξB~∫0ξ∫0θe(ξ−σ)A(τLA−1)eσA(Aeξ2B~A−1)Au~(tn+τ−ξ)dσdθdξ, the above estimates (5.6) and the bounds (2.2), (2.5), (4.5), (4.6), (4.7) and (5.1) yield   ‖Dn(2)‖ ≤16eτβ‖G‖M12K~τ3, ‖Dn(3)‖ ≤124Cinv2M13M2‖G‖2K~τ3,‖Dn(4)‖ ≤116Cinv2M13M2‖G‖2K~τ3, ‖Dn(5)‖ ≤116βeτβCinvM1M2‖G‖K~τ3. It remains to bound the integral $$D_n^{(1)}$$ defined in (5.5b). Using the representation   eξ−θ2B~[A,B~]eθ2B~−[A,B~]eξ2B~ =∫θξddν(eν−θ2B~[A,B~]eθ−ν2B~)eξ2B~dν =12∫θξeν−θ2B~[B~,[A,B~]]eξ+θ−ν2B~dν, we can write   Dn(1)=14∫0τeξ2B~∫0ξeξA∫θξeν−θ2B~[B~,[A,B~]]A−1(Aeξ+θ−ν2B~A−1)Au~(tn+τ−ξ)dνdθdξ. For the double commutator, we get from (2.2) and (4.7)   ‖[B~,[A,B~]]A−1‖ =‖2B~AB~A−1−B~2−AB~2A−1‖≤4Cinv2M14‖G‖2, which yields with (2.5) and (5.1)   ‖Dn(1)‖≤124eτβK~‖[B~,[A,B~]]A−1‖τ3≤16eτβCinv2M14‖G‖2K~τ3. The assertion now follows from $$D_n = D_n^{(1)} + D_n^{(2)} + D_n^{(3)} + D_n^{(4)} + D_n^{(5)}$$. □ We now investigate the local error further. Lemma 5.2 The dominating part $$\widehat \delta_n$$ of the local error $$\delta_n$$ defined in (5.3) satisfies   δ^n=δ^n(1)+δ^n(2),δ^n(1)=AΨ^Z1u~(tn),δ^n(2)=Z2Φ^A2u~(tn), where   ‖Z1‖≤112M1‖G‖τ3,‖Z2‖≤112CinvM1‖G‖τ3,‖AZ2‖≤112M1‖G‖τ3. Proof. The matrix $$L$$ can be written as   L =(AΨ^)ABΦ^+Ψ^BΦ^A2, since $$\widetilde{B}=\widehat{{\it{\Psi}}} B \widehat{{\it{\Phi}}}$$ with the notation (4.2). Since the filter functions $${\it{\Psi}}_S$$ and $${\it{\Phi}}$$ are independent of the integration variables, we have by (5.3)   δ^n =AΨ^Z1u~(tn)+Z2Φ^A2u~(tn), where   Z1 =12∫0τ∫0ξ∫0θe(ξ−σ)AABΦ^e(σ+τ−ξ)Adσdθdξ,Z2 =12∫0τ∫0ξ∫0θe(ξ−σ)AΨ^Be(σ+τ−ξ)Adσdθdξ. The desired bounds now follow immediately from (2.2), (2.5), (3.8a) and (4.4). □ Although the matrices $$Z_{1}$$, $$Z_2$$ and $$AZ_2$$ contain a factor of $$\tau^3$$, this is not sufficient to use a standard Lady Windermere’s fan argument to prove that the global error is of second order. The reason is the additional factors of $$A$$ in $$\widehat \delta_n^{(1)}$$ and $$\widehat \delta_n^{(2)}$$ which would yield a constant depending on $$\lVert {\it{\Omega}} \rVert$$. The proof of the following result treats the two terms defined in Lemma 5.2 separately in an appropriate way by using summation by parts. Theorem 5.3 (Global error of the averaged problem) Let the assumptions (3.8) and the finite-energy condition (4.6) be fulfilled. Then the global error of the splitting scheme (3.6) as an approximation to the solution of the averaged system (2.4) is bounded by   ‖un−u~(tn)‖≤Cτ2,0≤tn=nτ≤tend, (5.7) where $$C$$ only depends on $$C_\text{inv}$$, $$\| G \|$$, $$\widetilde{K}$$, $$M_j$$, $$j = 1,2,3$$ and $$t_{\text{end}}$$. Proof. By a telescopic identity, the global error can be written as   e~n=un−u~(tn)=(Sn−Tn)u0=∑j=0n−1Sn−j−1(S−T)Tju0=∑j=0n−1Sn−j−1δj with the local errors $$\delta_j$$ of Lemma 5.1. Lemmas 5.1 and 5.2 motivate to split the error into   e~n=e~n(1)+e~n(2)+e~n(3), (5.8) where   e~n(1)=∑j=0n−1Sn−j−1δ^j(1),e~n(2)=∑j=0n−1Sn−j−1δ^j(2),e~n(3)=∑j=0n−1Sn−j−1Dj. From the stability bound (5.2) and $$\|D_j\| \leq C \tau^3$$ by Lemma 5.1, we thus obtain   ‖e~n(3)‖=‖∑j=0n−1Sn−j−1Dj‖≤Cτ2(nτ)eβtend≤Cτ2. To bound $$\widetilde{e}_n^{(1)}$$ and $$\widetilde{e}_n^{(2)}$$, we write $$\tau A \widehat{\chi}$$ for $$\widehat{\chi} = \widehat{{\it{\Phi}}}$$ and $$\widehat{\chi} = \widehat{{\it{\Psi}}}$$ (see (4.2)) as   τAχ^=τχ^A=(eτA−I)Θ^χ=Θ^χ(eτA−I) (5.9) with   Θ^χ=θ^χ(τΩ),θ^χ(x)=12[cot⁡(x2)−11cot⁡(x2)][xχ(x)00xχ(x)] This representation follows from the identities   sin⁡(x)=cot⁡(x2)(1−cos⁡(x))andcos⁡(x)+1=cot⁡(x2)sin⁡(x). Note that by assumptions (3.8c) and (3.8d) on the filter functions $$\chi$$, we have   ‖Θ^χ‖≤C. (5.10) Next we use summation by parts. With   Ej=∑k=0j−1Sk,Fj=∑k=0j−1Tk, it holds   e~n(1)=Enδ^0(1)+∑j=0n−2En−j−1(δ^j+1(1)−δ^j(1)). Using Lemma 5.2 and (5.9) to replace $$\widehat \delta_{j}^{(1)}$$ by $$(e^{\tau A} - I) \widehat{{\it{\Theta}}}_\psi \tfrac1{\tau} Z_1 T^j u_0$$, this implies   e~n(1) =En(eτA−I)Θ^ψ1τZ1u0+∑j=0n−2En−j−1(eτA−I)Θ^ψZ11τ(T−I)Tju0. (5.11) Now we estimate the matrices in this expression. The matrices $$\widehat{{\it{\Theta}}}_\psi$$ and $$Z_1$$ can be estimated with (5.10) and Lemma 5.2, respectively. To bound $$E_{n-j-1} (e^{\tau A} - I)$$ in (5.11), we start from   Ej(eτA−I)=Ej(eτA−S)+Ej(S−I)=Ej(eτA−S)+Sj−I, and we use the stability bound (5.2) and $$\| e^{\tau A} - S \| \leq C \tau$$ by (2.5), (3.5) and (5.1) to show that   ‖Ej(eτA−I)‖≤C. To bound $$\tfrac1{\tau} (T-I) T^ju_0$$ in (5.11), we start from   1τ(T−I)Tju0 =1τ(eτA−I)Tju0+1τ(T−eτA)Tju0 =(eτA−I)(τA)−1Au~(tj)+1τ∫0τe(τ−ξ)AB~A−1Au~(tj+ξ)dξ by the variation-of-constants formula (4.1b), and we use (2.2), the finite-energy condition (4.6) and $$\|(e^{\tau A}-I) (\tau A)^{-1}\| = \| \int_0^1 {\it \,{e}}^{\sigma \tau A} {\rm d}\sigma\| \leq 1$$ by (2.5) to show that   ‖1τ(T−I)Tju0‖≤C. This yields the estimate $$\|\widetilde{e}_n^{(1)}\| \leq C \tau^2$$ for (5.11) since $$\lVert Z_1 \rVert \leq C \tau^3$$ by Lemma 5.2. Analogously, we have by Lemma 5.2 and (5.9)   e~n(2) =1τZ2Θ^ϕA(eτA−I)Fnu0+∑j=0n−2Sj1τ(S−I)Z2Θ^ϕA(eτA−I)Fn−j−1u0. Using the variation-of-constants formula (4.1b) yields   A(eτA−I)Fju0 =A(eτA−T)Fju0+A(Tj−I)u0 =−∫0τe(τ−ξ)AAB~eξ(A+B~)Fju0dξ+A(u~(tj)−u0). The finite-energy condition (4.6) shows   ‖τeξ(A+B~)Fju0‖=‖τ∑k=0j−1u~(tk+ξ)‖≤tendCinvK~. By (4.7) we thus obtain   ‖A(eτA−I)Fju0‖≤C. Moreover, by Lemma 5.2 we have   ‖1τ(S−I)Z2‖ ≤‖1τ(S−eτA)Z2‖+‖(eτA−I)(τA)−1AZ2‖≤Cτ3 since $$\|S-e^{\tau A}\|\le C\tau$$ and $$\| (e^{\tau A} - I) (\tau A)^{-1}\|\le 1$$. Together with the stability estimate (5.2) and the bound (5.10), this proves $$\| \widetilde{e}_n^{(2)} \| \leq C \tau^2$$, and thus (5.7) by (5.8). □ Remark It is also possible to prove Theorem 5.3 if (5.3) in Lemma 5.1 is replaced by the representation of the local error which has been derived in Jahnke & Lubich (2000, Theorem 2.1). While they used quadrature errors to bound the local error, we used the fundamental theorem of calculus as in (5.4). An advantage of the new representation is that $$\left\lVert G \right\rVert$$ appears only quadratically in the bounds, while it appears cubically in Jahnke & Lubich (2000). We also believe that the representation (5.3) is more suitable for a future extension of our approach to nonlinear problems. 6. Discussion and comparison Error bounds for trigonometric integrators applied to second-order oscillatory differential equations have been previously derived in Grimm & Hochbruck (2006) and, for the case of a single high frequency, in Theorems XIII.4.1 and XIII.4.2 of Hairer et al. (2006). They differ from the error bounds of the present paper (Theorem 3.2) in their statement, in the required assumptions on the filter functions and, most notably, in the technique of proof. The difference in the statement concerns the error bound for the velocities. While the error bounds of Grimm & Hochbruck (2006) and Hairer et al. (2006) for the velocities $$q'$$ are of first order in $$\tau$$, the error bound of Theorem 3.2 (in the linear case) is of second order in the rescaled velocities $${\it{\Omega}}^{-1} q'$$. For small time step sizes $$\tau \left\lVert {{\it{\Omega}}} \right\rVert=O(1)$$, the latter bound implies the former. Besides this difference in the statement, there are differences in the assumptions on the filter functions. The splitting method (3.6) with filter functions $$\phi$$ and $$\psi_S$$ coincides with the trigonometric integrator in Grimm & Hochbruck (2006) and Hairer et al. (2006) if the filter functions $$\phi$$, $$\psi$$, $$\psi_1$$ and $$\psi_0$$ in Grimm & Hochbruck (2006) and Hairer et al. (2006) are chosen in the following way:   ϕ =ϕ,ψ1=ψS,ψ=sinc⁡(⋅)ψS,ψ0=cos⁡(⋅)ψS. Now, we can compare the conditions (11)–(16) in Grimm & Hochbruck (2006), and the conditions (XIII.4.1) and (XIII.4.8) in Hairer et al. (2006) to our conditions (3.8). Our condition (3.8b) on the boundedness of the filter functions coincides identically with $$(11)$$ of Grimm & Hochbruck (2006). The conditions (13)–(16) of Grimm & Hochbruck (2006) imply that $$\chi=\phi$$ and $$\chi=\psi_S$$ are zero whenever $$\sin(\frac{\cdot}2)$$ is zero, meaning for $$x = 2k\pi$$, $$k \in \mathbb Z$$. This behaviour can also be found in our condition (3.8d). Still, an analogon to our condition (3.8c) is missing in Grimm & Hochbruck (2006) since this condition requires that $$\chi$$ decreases at least like $$x^{-1}$$ for $$x \to \infty$$. In comparison to Hairer et al. (2006), our conditions (3.8b)–(3.8d) are implied by the conditions (XIII.4.1) and (XIII.4.8) used there, which require in particular that $$|\chi(x)|$$ and $$|x\sin(\frac12 x)^{-1}\chi(x)|$$ are bounded. Note that all sets of conditions are proven to be sufficient but not necessary. But the main difference in comparison to Grimm & Hochbruck (2006) and Hairer et al. (2006) is the technique of proof. It will be interesting to see whether the technique developed in the present paper, namely a Lady Windermere’s fan that takes cancellations in the error accumulation into account, can help to gain further insight into the error behaviour of trigonometric integrators and splitting methods, in particular for nonlinear problems. Acknowledgements The authors thank Roland Schnaubelt for helpful discussions on the representation of the local error. Funding We gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through RTG 1294, CRC 1114, CRC 1173 and project GA 2073/2-1. References Blanes, S. ( 2015) Explicit symplectic RKN methods for perturbed non-autonomous oscillators: splitting, extended and exponentially fitting methods. Comput. Phys. Commun.,  193, 10– 18. Google Scholar CrossRef Search ADS   Blanes, S. & Casas, F. ( 2016) A concise introduction to geometric numerical integration . A Chapman & Hall book: Monographs and research notes in mathematics. Boca Raton: CRC Press. Google Scholar CrossRef Search ADS   Cohen, D., Gauckler, L., Hairer, E. & Lubich, C. ( 2015) Long-term analysis of numerical integrators for oscillatory Hamiltonian systems under minimal non-resonance conditions. BIT,  55, 705– 732. Google Scholar CrossRef Search ADS   Einkemmer, L. & Ostermann, A. ( 2014) Convergence analysis of Strang splitting for Vlasov-type equations. SIAM J. Numer. Anal.,  52, 140– 155. Google Scholar CrossRef Search ADS   Einkemmer, L. & Ostermann, A. ( 2015) Overcoming order reduction in diffusion-reaction splitting. Part 1: Dirichlet boundary conditions. SIAM J. Sci. Comput.,  37, A1577–A1592. Faou, E. ( 2012) Geometric numerical integration and Schrödinger equations.  Zurich Lectures in Advanced Mathematics. Zürich: European Mathematical Society (EMS), pp. xiii+138. Faou, E., Ostermann, A. & Schratz, K. ( 2015) Analysis of exponential splitting methods for inhomogeneous parabolic equations. IMA J. Numer. Anal.,  35, 161– 178. Google Scholar CrossRef Search ADS   García-Archilla, B., Sanz-Serna, J. M. & Skeel, R. D. ( 1999) Long-time-step methods for oscillatory differential equations. SIAM J. Sci. Comput.,  20, 930– 963. Google Scholar CrossRef Search ADS   Gauckler, L. ( 2015) Error analysis of trigonometric integrators for semilinear wave equations. SIAM J. Numer. Anal.,  53, 1082– 1106. Google Scholar CrossRef Search ADS   Grimm, V. & Hochbruck, M. ( 2006) Error analysis of exponential integrators for oscillatory second-order differential equations. J. Phys. A,  39, 5495– 5507. Google Scholar CrossRef Search ADS   Grimm, V. & Hochbruck, M. ( 2008) Rational approximation to trigonometric operators. BIT,  48, 215– 229. Google Scholar CrossRef Search ADS   Hairer, E., Lubich, C. & Wanner, G. ( 2006) Geometric numerical integration : structure-preserving algorithms for ordinary differential equations . Springer series in computational mathematics; 31, 2. edn. Berlin: Springer. Hansen, E. & Ostermann, A. ( 2009) Exponential splitting for unbounded operators. Math. Comp.,  78, 1485– 1496. Google Scholar CrossRef Search ADS   Hochbruck, M., Jahnke, T. & Schnaubelt, R. ( 2015) Convergence of an ADI splitting for Maxwell’s equations. Numer. Math.,  129, 535– 561. Google Scholar CrossRef Search ADS   Hochbruck, M. & Lubich, C. ( 1999) A Gautschi-type method for oscillatory second-order differential equations. Numer. Math.,  83, 403– 426. Google Scholar CrossRef Search ADS   Holden, H., Karlsen, K. H., Lie, K.-A. & Risebro, N. H. ( 2010) Splitting methods for partial differential equations with rough solutions . EMS Series of Lectures in Mathematics. Zürich: European Mathematical Society (EMS), pp. viii+226. Google Scholar CrossRef Search ADS   Holden, H., Lubich, C. & Risebro, N. H. ( 2013) Operator splitting for partial differential equations with Burgers nonlinearity. Math. Comp.,  82, 173– 185. Google Scholar CrossRef Search ADS   Hundsdorfer, W. H. & Verwer, J. G. ( 2007) Numerical solution of time dependent advection diffusion reaction equations . Springer series in computational mathematics; 33, corr. 2. print. edn. Berlin, Heidelberg: Springer. Ixaru, L. G. & Vanden Berghe, G. ( 2004) Exponential fitting . Mathematics and its Applications, vol. 568. Kluwer Academic Publishers, Dordrecht, pp. xiv+308. Google Scholar CrossRef Search ADS   Jahnke, T. & Lubich, C. ( 2000) Error bounds for exponential operator splittings. BIT,  40, 735– 744. Google Scholar CrossRef Search ADS   Koch, O. & Lubich, C. ( 2011) Variational-splitting time integration of the multi-configuration time-dependent Hartree-Fock equations in electron dynamics. IMA J. Numer. Anal.,  31, 379– 395. Google Scholar CrossRef Search ADS   Lubich, C. ( 2008) On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations. Math. Comp.,  77, 2141– 2153. Google Scholar CrossRef Search ADS   McLachlan, R. I. & Quispel, G. R. W. ( 2002) Splitting methods. Acta Numerica,  11, 341– 434. Google Scholar CrossRef Search ADS   Paternoster, B. ( 2012) Present state-of-the-art in exponential fitting. A contribution dedicated to Liviu Ixaru on his 70th birthday. Comput. Phys. Commun.,  183, 2499– 2512. Google Scholar CrossRef Search ADS   Strang, G. ( 1968) On the construction and comparison of difference schemes. SIAM J. Numer. Anal.,  5, 506– 517. Google Scholar CrossRef Search ADS   Thalhammer, M., Caliari, M. & Neuhauser, C. ( 2009) High-order time-splitting Hermite and Fourier spectral methods. J. Comput. Phys.,  228, 822– 832. Google Scholar CrossRef Search ADS   Wu, X., You, X. & Wang, B. ( 2013) Structure-preserving algorithms for oscillatory differential equations . Springer, Heidelberg; Science Press Beijing, Beijing, pp. xii+236. Google Scholar CrossRef Search ADS   Wu, X., Liu, K. & Shi, W. ( 2015) Structure-preserving algorithms for oscillatory differential equations. II . Springer, Heidelberg; Science Press Beijing, Beijing, pp. xv+298. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Closing the gap between trigonometric integrators and splitting methods for highly oscillatory differential equations

, Volume 38 (1) – Jan 1, 2018
18 pages

/lp/ou_press/closing-the-gap-between-trigonometric-integrators-and-splitting-EjgcoxJkyE
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx007
Publisher site
See Article on Publisher Site

### Abstract

Abstract An error analysis of symmetric trigonometric integrators applied to highly oscillatory linear second-order differential equations is given. Second-order convergence is shown uniformly in the high frequencies under a finite-energy condition on the exact solution. The main novelty is the concept of the proof itself, which is based on the interpretation of trigonometric integrators as splitting methods for averaged differential equations. This allows one to combine techniques for splitting methods with those for trigonometric integrators. For the bound of the global error, cancellations in the error accumulation have to be studied carefully. 1. Introduction Many ordinary differential equations in science can be written as   u′=F1(u)+F2(u) (1.1) with linear or nonlinear vector fields $$F_1$$ and $$F_2$$, and with the property that solving the two sub-problems   v′=F1(v)andw′=F2(w) exactly or approximately is numerically much more efficient than solving (1.1) directly. Typical examples include many-body problems and other Hamiltonian systems with a certain structure, but also semidiscretizations of time-dependent partial differential equations such as, e.g., linear or nonlinear Schrödinger equations, linear or nonlinear wave equations, or Maxwell’s equations. In such a situation, splitting methods are easy to implement and very efficient, because these integrators compute a numerical approximation of the full problem by a composition of the flows of the two sub-problems. In many cases, the numerical solution inherits certain geometric properties of the exact solution, such as norm conservation or symplecticity of the flow. The order conditions and the accuracy of splitting methods have been analysed in many papers such as, e.g., Strang (1968), Jahnke & Lubich (2000), Lubich (2008), Hansen & Ostermann (2009), Thalhammer et al. (2009), Koch & Lubich (2011), Holden et al. (2013), Einkemmer & Ostermann (2014), Einkemmer & Ostermann (2015), Faou et al. (2015), Hochbruck et al. (2015) and references therein. Overviews have been given, e.g., in McLachlan & Quispel (2002), Hairer et al. (2006), Hundsdorfer & Verwer (2007), Holden et al. (2010) and Blanes & Casas (2016). The classical procedure is to estimate the local error by proving bounds for certain iterated commutators between the two vector fields. The error bound for the global error then follows by standard arguments like Gronwall inequalities or Lady Windermere’s fan. Unfortunately, the accuracy typically suffers when splitting methods are applied to problems with highly oscillatory solutions. In order to obtain a reasonable approximation, the step size must be considerably smaller than the inverse of the highest frequency, which reduces the efficiency considerably. We remark, however, that certain invariants of the exact flow are often conserved over long times even in the presence of oscillations (see e.g., Hairer et al., 2006; Faou, 2012; Cohen et al., 2015). In the present article, we study the discretization error on finite time intervals and do not consider the additional numerical resonances that show up on these longer time scales. A particular class of problems with highly oscillatory solutions takes the semilinear form   q″(t)=−Ω2q(t)+g(q(t)),t>0,q(0)=q0,q′(0)=q0′, (1.2) where $${\it{\Omega}}$$ is a symmetric positive semidefinite matrix of arbitrary large norm, and where $$g$$ is smooth and bounded. For such problems, trigonometric integrators have been constructed and analysed in GarcÍa-Archilla et al. (1999), Hochbruck & Lubich (1999), Hairer et al. (2006) and Grimm & Hochbruck (2006) under a finite-energy condition. These methods involve filter functions which are chosen in such a way that oscillatory parts of the local error do not sum up in the global error. To prove this in the error analysis is a delicate matter, and usually excludes to apply the technique of Lady Windermere’s fan in a standard way. An exception is semilinear wave equations with polynomial or analytic nonlinearities, see Gauckler (2015). Hence, it seems that splitting methods and trigonometric integrators are two different approaches, based on different ideas, having different properties, and requiring different techniques for their analysis. Nevertheless, there is a bridge between these two worlds: it is long known that symmetric trigonometric integrators applied to (1.2) can be interpreted as a Strang splitting method applied to an averaged version of the first-order formulation of (1.2). This interesting link raises a number of questions: Is it possible to gain a better understanding by considering trigonometric integrators as splitting methods for averaged equations? Does this relation allow to apply the techniques developed for one class of methods to the other one? And, most importantly, is it then possible to construct and analyze efficient numerical integrators for fully nonlinear problems? These questions are our motivation for the analysis below. In this article, we will restrict ourselves to the linear variant of (1.2), i.e., $$g(q) = Gq$$ with a matrix $$G$$ with a moderate norm $$\| G \| \ll \| {\it{\Omega}}\|$$. Here and in the following, $$\| \cdot \|$$ denotes the Euclidean vector norm or its induced matrix norm, respectively. The analysis in Lubich (2008) shows that for splitting methods the calculus of Lie derivatives allows, at least to some extent, to carry over techniques developed for linear differential equations to nonlinear ones. After reformulation as a first-order problem, we will derive the corresponding averaged equation and prove a bound for the difference between the solution of both equations; cf. Theorem 4.1. Then, we present an error analysis for the classical Strang splitting applied to the averaged equation, which yields a result very similar, but not equivalent to the one obtained in Grimm & Hochbruck (2006). We stress that the novelty of our analysis is not the error bound itself, but the fact that it is proven by techniques which, to the best of our knowledge, have so far not been considered in the context of trigonometric integrators. 2. Problem setting, assumptions and notation The situation we have in mind is that the ordinary differential equation   q″(t)=−Ω2q(t)+Gq(t),0≤t≤tend,q(0)=q0,q′(0)=q0′ (2.1) stems from a spatial discretization of a linear wave equation with finite elements, finite differences or spectral methods on a family of finer and finer meshes. In this situation, the matrix $$-{\it{\Omega}}^2$$ stems from a spatial discretization of the Laplacian or a more general differential operator, and $${\it{\Omega}}$$ is a symmetric, positive definite matrix (possibly after shifting $${\it{\Omega}} \rightarrow \sqrt{{\it{\Omega}}^2+I}$$ and $$G\rightarrow G+I$$). Then $$\left\lVert {{\it{\Omega}}} \right\rVert$$ becomes arbitrarily large if the spatial mesh width tends to zero, but $$\lVert {\it{\Omega}}^{-1} \rVert$$ remains uniformly bounded independently of the discretization. This motivates the following assumption. Assumption 2.1 Let $${\it{\Omega}} \in \mathcal{F}$$, where $$\mathcal{F}$$ is a family of symmetric, positive definite matrices such that there is a constant $$C_\text{inv}$$ with   ‖Ω−1‖≤Cinvfor all Ω∈F. (2.2) For the matrix $$G$$ we assume that its norm is bounded independently of $${\it{\Omega}}$$. Our aim is to prove error estimates which are uniform for all matrices in the family $$\mathcal{F}$$, which means that they are independent of $$\left\lVert {{\it{\Omega}}} \right\rVert$$ (i.e., independent of the spatial discretization). On the other hand, the constant $$C_\text{inv}$$ only depends on the coercivity constant of the corresponding differential operator. For the solution we rely on the following assumption. Assumption 2.2 Let the solution $$q:[0,t_{\text{end}}] \rightarrow \mathbb R^d$$ of (2.1) fulfill the finite-energy condition  ‖Ωq(t)‖2+‖q′(t)‖2≤K2,0≤t≤tend (2.3) with a constant $$K>0$$ on a finite time interval of length $$t_{\text{end}}$$. In fact, it can be shown that (2.3) is true on bounded time intervals if the initial data satisfy the bound $$\left\lVert {{\it{\Omega}} q(0)} \right\rVert^2 + \left\lVert {q'(0)} \right\rVert^2 \leq K_0^2$$ with a sufficiently small $$K_0 \leq K$$; cf. Lemma 4.2 below. In order to apply a splitting scheme, we formulate (2.1) as a first-order problem. We define the new variable   u=[qΩ−1q′],u0=[q0Ω−1q0′] which solves the differential equation   u′=Au+Bu,u(0)=u0 (2.4) with matrices   A=[0Ω−Ω0],B=[00Ω−1G0]. Since $$A$$ is skew symmetric, the exponential   exp⁡(tA)=[cos⁡(tΩ)sin⁡(tΩ)−sin⁡(tΩ)cos⁡(tΩ)] is unitary and thus   ‖etA‖=1,t∈R. (2.5) The finite-energy condition (2.3) is equivalent to   ‖Au(t)‖≤K,0≤t≤tend. (2.6) The solution of (2.4) could, in principle, be approximated with the classical Strang splitting, but in our setting this method only yields an acceptable accuracy if $$\tau \ll \|{\it{\Omega}}\|^{-1}$$. Such a severe step-size restriction is not acceptable in practice. We will show that much better approximations are obtained if the Strang splitting is applied to an averaged version of (2.4). The methods considered below involve analytic (mostly trigonometric) matrix functions of $$\tau{\it{\Omega}}$$. We assume that products of such matrix functions with a given vector can be computed efficiently, e.g., via fast Fourier transforms or (rational) Krylov subspace methods (see e.g., Grimm & Hochbruck, 2008). 3. Trigonometric integrators as splitting methods Symmetric one-step trigonometric integrators (Hairer et al., 2006, Section XIII.2.2) and Grimm & Hochbruck (2006) written in terms of the variables $$u_n=[q_n, {\it{\Omega}}^{-1}q'_n] \approx u(t_n)=[q(t_n), {\it{\Omega}}^{-1}q'(t_n)]$$ at $$t_n = n\tau$$ are given by   un+1 =eτAun+τ2[sin⁡(τΩ)Ω−1G~qncos⁡(τΩ)Ω−1G~qn+Ω−1G~qn+1]. (3.1) Here, $$\tau>0$$ denotes the step-size, and the matrix $$\widetilde{G}$$ is defined as   G~=ΨSGΦ,Φ=ϕ(τΩ),ΨS=ψS(τΩ) (3.2) with appropriate scalar filter functions $$\phi$$ and $$\psi_S$$. Note that $$\psi_S$$ was denoted by $$\psi_1$$ in Hairer et al. (2006) and Grimm & Hochbruck (2006). In these references, it was shown that under the assumptions of Section 2 and certain (sufficient) conditions on the filter functions $$\phi$$ and $$\psi_S$$, these methods yield a global accuracy of $$O(\tau^2)$$ in the positions $$q_n$$ and $$O(\tau)$$ in the velocities $$q_n'$$. The constants in the error bounds are independent of $$\| {\it{\Omega}} \|$$, in spite of the oscillatory nature of the solution. It is well known that symmetric trigonometric integrators can be interpreted as splitting schemes. With   B~=[00Ω−1G~0], (3.3) the method (3.1) can be written as   (I−τ2B~)un+1 =eτA(I+τ2B~)un. (3.4) Now we use that   eθB~=I+θB~for all θ∈R (3.5) due to the particular structure of $$\widetilde{B}$$. Substituting (3.5) into (3.4) gives   un+1=Sun,S=eτ2B~eτAeτ2B~. (3.6) Hence, the trigonometric integrator (3.1) is equivalent to the exponential Strang splitting method applied to the averaged equation   u~′=Au~+B~u~,u~(0)=u~0=u0 (3.7) with $$\widetilde{B}$$ defined in (3.3) and (3.2). This means that applying the Strang splitting to the averaged equation yields very good approximations for step sizes, where the Strang splitting applied to the original problem (2.4) fails if the norm of $${\it{\Omega}}$$ is large. We point out that not only the numerical solution $$u_n$$ depends on the step-size $$\tau$$, but also the exact solution $$\widetilde{u}$$ of the averaged equation via the filter functions $${\it{\Phi}}=\phi(\tau {\it{\Omega}})$$ and $${\it{\Psi}}_S=\psi_S(\tau{\it{\Omega}})$$. Note that for semilinear problems (1.2) the trigonometric integrators can still be interpreted as a Strang splitting applied to a (semilinear) averaged equation. The interpretation of trigonometric integrators as second-order Strang splitting methods is well known (see GarcÍa-Archilla et al., 1999 or Hairer et al., 2006). Recently, it has been shown in Blanes (2015) that also certain methods of higher order designed for oscillatory differential equations of the form (1.2) can be interpreted as higher order splitting methods. The methods considered in Blanes (2015) are exponentially fitted Runge–Kutta–Nyström methods and extended Runge–Kutta–Nyström methods (see also the monographs Ixaru & Vanden Berghe, 2004; Wu et al., 2013; Wu et al., 2015 and the review article Paternoster, 2012). We emphasize, however, that higher order methods cannot be expected to converge with order higher than two under the finite-energy condition (2.6). Instead, third or higher order error bounds are expected to require a stronger condition on the exact solution of the form $$\left\lVert A^q u(t) \right\rVert \leq K$$ with $$q\ge 2$$, which amounts in particular to a small-energy assumption if the matrix $${\it{\Omega}}$$ contains only high frequencies. Our aim is to prove second-order error bounds for the trigonometric integrator (3.1) on the basis of the interpretation as a splitting method (3.6) and under the finite-energy condition (2.6), using a carefully adapted Lady Windermere’s fan argument familiar from splitting integrators. In order to analyse the error of the splitting scheme, we first study the properties of the solution of the averaged problem (3.7). This allows us to bound the error which results from solving the averaged equation instead of (2.4) (Section 4). To analyse the error of the splitting scheme for the averaged problem, we first give a new representation of the local error. Unfortunately, this local error still contains a term which is not uniformly of third order in $$\tau$$ and requires a more careful investigation (Section 5). The error analysis given below relies on the following assumption on the filter functions. Assumption 3.1 The filter functions $$\chi = \phi$$ or $$\chi = \psi_S$$ are even analytic functions with the properties   χ(0) =1, (3.8a)  |χ(x)| ≤M1, (3.8b)  |xχ(x)| ≤M2, (3.8c)  |cot⁡(x2)xχ(x)| ≤M3 (3.8d) for certain constants $$M_j$$ uniformly for all $$x \in \mathbb R$$. Even functions $$\chi$$ guarantee that the scheme is symmetric. A popular example used in trigonometric integrators is $$\chi(x) = \mathop{\mathrm{sinc}}(x)$$. Note that (3.8a) and the condition that $$\chi$$ is even analytic imply   |x−2(1−χ(x))|≤M4. (3.8e) In the following $$C$$ denotes a generic constant which may have different values at different occurrences. Our main result is stated in the following theorem. Theorem 3.2 (Main result) Let Assumptions 2.1, 2.2 and 3.1 be fulfilled. Then the global error of (3.6) applied to (2.4) is bounded by   ‖un−u(tn)‖≤Cτ2,0≤tn=nτ≤tend with a constant $$C$$ that only depends on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_j$$, $$j=1,\ldots,4$$ and $$t_{\text{end}}$$ but not on $$\left\lVert {{\it{\Omega}}} \right\rVert$$. Proof. We combine the results from Theorem 4.1 and Theorem 5.3 below to obtain   ‖un−u(tn)‖≤‖un−u~(tn)‖+‖u~(tn)−u(tn)‖≤Cτ2,0≤tn=nτ≤tend, where the constant $$C$$ only depends on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_j$$, $$j=1,\ldots,4$$ and $$t_{\text{end}}$$. □ 4. Properties of the averaged equation In this section, we prove that the solutions of the original problem (2.4) and the averaged equation (3.7) differ only by $$O(\tau^2)$$. Theorem 4.1 Let the assumptions (3.8a) and (3.8b) and the finite-energy condition (2.3) be fulfilled. Let $$u$$ be the solution of (2.4) and $$\widetilde{u}$$ be the solution of the averaged system (3.7). Then it holds   ‖u(t)−u~(t)‖≤Cavτ2,0≤t≤tend, where $$C_\text{av}$$ only depends on $$C_\text{inv}$$, $$\|G\|$$, $$K$$, $$M_1$$, $$M_4$$ and $$t_{\text{end}}$$. Proof. The variation-of-constants formula yields   u(t) =etAu0+∫0te(t−s)ABu(s)ds, (4.1a)  u~(t) =etAu0+∫0te(t−s)AB~u~(s)ds. (4.1b) We define the block diagonal matrices   Ψ^=[ΨS00ΨS],Φ^=[Φ00Φ]. (4.2) Obviously, they satisfy $$A\widehat{{\it{\Psi}}}=\widehat{{\it{\Psi}}} A$$ and $$A\widehat{{\it{\Phi}}}=\widehat{{\it{\Phi}}} A$$. We also have $$\widetilde{B} = \widehat{{\it{\Psi}}} B \widehat{{\it{\Phi}}}$$ so that   u(t)−u~(t) =∫0te(t−s)A((I−Ψ^)Bu(s)+Ψ^B(I−Φ^)u(s)+Ψ^BΦ^(u(s)−u~(s)))ds =I1(t)+I2(t)+∫0te(t−s)AB~(u(s)−u~(s))ds with   I1(t) =∫0te(t−s)A(I−Ψ^)Bu(s)ds,I2(t)=∫0te(t−s)AΨ^B(I−Φ^)u(s)ds. For the first term integration by parts yields   I1(t) =−τ2e(t−s)A(τA)−2(I−Ψ^)ABu(s)|0t+τ2∫0te(t−s)A(τA)−2(I−Ψ^)ABu′(s)ds. With (2.2) and (2.6) it follows that   ‖u‖ =‖A−1Au‖≤CinvK, (4.3a)  ‖u′‖ ≤‖Au‖+‖Bu‖≤K+Cinv2‖G‖K. (4.3b) Using in addition (2.5), (3.8e) and   ‖AB‖=‖G‖, (4.4) this yields the estimate $$\lVert I_1(t) \rVert \leq C \tau^2$$ with a constant $$C$$ depending only on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_4$$ and $$t_{\text{end}}$$. For the second term, we insert the variation of constants formula once more and obtain   I2(t) =∫0te(t−s)AΨ^B(I−Φ^)u(s)ds =∫0te(t−s)AΨ^B(I−Φ^)esAu0ds+∫0t∫0se(t−s)AΨ^B(I−Φ^)e(s−σ)ABu(σ)dσds. Integration by parts yields   I2(t) =e(t−s)AΨ^B(I−Φ^)A−1esAu0|0t+∫0te(t−s)AAΨ^B(I−Φ^)A−1esAu0ds −∫0te(t−s)AΨ^B(I−Φ^)A−1e(s−σ)ABu(σ)|0sds +∫0t∫0se(t−s)AΨ^B(I−Φ^)A−1e(s−σ)ABu′(σ)dσds =τ2e(t−s)AΨ^B(I−Φ^)(τA)−2esAAu0|0t+τ2∫0te(t−s)AΨ^AB(I−Φ^)(τA)−2esAAu0ds −τ2∫0te(t−s)AΨ^B(I−Φ^)(τA)−2e(s−σ)AABu(σ)|0sds +τ2∫0t∫0se(t−s)AΨ^B(I−Φ^)(τA)−2e(s−σ)AABu′(σ)dσds. By (2.2), (2.5), (3.8), (4.4) and (4.3) we obtain $$\lVert I_2(t) \rVert \leq C \tau^2$$, where $$C$$ depends on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_1$$, $$M_4$$ and $$t_{\text{end}}$$. Hence, we have   ‖u(t)−u~(t)‖ =‖I1(t)+I2(t)+∫0te(t−s)AB~(u−u~)ds‖≤Cτ2+∫0tβ‖u(s)−u~(s)‖ds with   β=Cinv‖G‖M12≥‖B~‖. (4.5) An application of Gronwall’s Lemma proves the desired result. □ We finally show that the solution of the averaged problem (3.7) inherits the finite-energy condition (2.3) of the original problem. Lemma 4.2 Let the assumption (3.8b) and the finite-energy condition (2.3) be fulfilled. Then the solution $$\widetilde{u}$$ of the averaged system (3.7) satisfies the finite-energy condition   ‖Au~(t)‖≤K~,0≤t≤tend, (4.6) where $$\widetilde{K}$$ depends only on $$C_\text{inv}$$, $$\lVert G \rVert$$, $$K$$, $$M_1$$ and $$t_{\text{end}}$$. Proof. Using (2.5) and the variation of constants formula (4.1b) yields   ‖Au~(t)‖ ≤‖Au0‖+∫0t‖AB~‖‖A−1‖‖Au~(s)‖ds. By definition (3.3) and (3.8b), it holds   ‖AB~‖=‖G~‖≤M12‖G‖ (4.7) and the statement thus follows from (2.2) and Gronwall’s Lemma. □ 5. Finite-time error analysis of the splitting scheme In this section we finally study the error of the splitting scheme. We use the notation   S=eτ2B~eτAeτ2B~ defined in (3.6) for the numerical flow over a time step $$\tau$$ and   T=eτ(A+B~) for the exact flow of the averaged equation (3.7). By (3.3) and (3.5) we have   max{‖eθB~‖,‖AeθB~A−1‖}≤1+θβ≤eθβ (5.1) for $$\theta\ge 0$$ with $$\beta$$ from (4.5). Together with (2.5), this provides the stability estimate   ‖Sn‖≤(1+τ2β)2n≤eβtend,0≤tn≤tend. (5.2) We next consider the local error. Lemma 5.1 (Local error) Assume that the filter functions are bounded as stated in (3.8b) and (3.8c). If the finite-energy condition (4.6) holds true, then the local error at time $$t_n=n\tau$$, $$0 \leq t_n \le t_{\text{end}}-\tau$$, of the splitting method (3.6) as an approximation to the averaged system (3.7) is given by   δn=(S−T)u~(tn)=δ^n+Dn, where   δ^n=12∫0τ∫0ξ∫0θe(ξ−σ)ALe(σ+τ−ξ)Au~(tn)dσdθdξ,L=A2B~+B~A2 (5.3) and $$\lVert D_n \rVert \leq C \tau^3$$. The constant $$C$$ depends on $$C_\text{inv}$$, $$\|G\|$$, $$\widetilde{K}$$, $$M_1$$, $$M_2$$ and $$t_{\text{end}}$$. Proof. We start with the following representation of the local error   δn =∫0τddξ(eξ2B~eξAeξ2B~e(τ−ξ)(A+B~))u~(tn)dξ =12∫0τeξ2B~([B~,eξA]eξ2B~+2eξA[A,eξ2B~])e(τ−ξ)(A+B~)u~(tn)dξ, (5.4) where $$\left[ A_1,A_2 \right]=A_1A_2 - A_2A_1$$ denotes the commutator of two matrices $$A_1$$ and $$A_2$$. The commutators in $$\delta_n$$ are written as integrals   [B~,eξA] =∫0ξddθ(e(ξ−θ)AB~eθA)dθ=−∫0ξe(ξ−θ)A[A,B~]eθAdθ and   [A,eξ2B~] =∫0ξddθ(eξ−θ2B~Aeθ2B~)dθ=12∫0ξeξ−θ2B~[A,B~]eθ2B~dθ. This yields   δn =12∫0τeξ2B~∫0ξ(−e(ξ−θ)A[A,B~]eθAeξ2B~+eξAeξ−θ2B~[A,B~]eθ2B~)e(τ−ξ)(A+B~)u~(tn)dθdξ. We add and subtract the term $${\it \,{e}}^{\xi A} [A,\widetilde{B}] {\it \,{e}}^{\frac{\xi}2 \widetilde{B}}$$ inside the brackets, and the local error can be split into $$\delta_n = \delta_n^{(1)} + D_n^{(1)}$$, where   δn(1) =12∫0τeξ2B~∫0ξ(eξA[A,B~]−e(ξ−θ)A[A,B~]eθA)eξ2B~e(τ−ξ)(A+B~)u~(tn)dθdξ, (5.5a)  Dn(1) =12∫0τeξ2B~∫0ξeξA(eξ−θ2B~[A,B~]eθ2B~−[A,B~]eξ2B~)e(τ−ξ)(A+B~)u~(tn)dθdξ. (5.5b) The term in brackets in (5.5a) can be written as   eξA[A,B~]−e(ξ−θ)A[A,B~]eθA =−∫0θddσ(e(ξ−σ)A[A,B~]eσA)dσ =∫0θe(ξ−σ)A[A,[A,B~]]eσAdσ. For the double commutator, we have   [A,[A,B~]] =L−R,L=A2B~+B~A2,R=2AB~A, and we split $$\delta_n^{(1)}$$ accordingly into $$\delta_n^{(1)} = \delta_n^{(2)}+D_n^{(2)}$$ with   δn(2) =12∫0τeξ2B~∫0ξ∫0θe(ξ−σ)ALeσAeξ2B~e(τ−ξ)(A+B~)u~(tn)dσdθdξ,Dn(2) =−12∫0τeξ2B~∫0ξ∫0θe(ξ−σ)AReσAeξ2B~e(τ−ξ)(A+B~)u~(tn)dσdθdξ. To extract the term $$\widehat \delta_n$$ defined in (5.3) from $$\delta_n^{(2)}$$, we use (3.5) and once more the variation-of-constants formula to obtain   δn(2) =δ^n+Dn(3)+Dn(4)+Dn(5), where   Dn(3) =12∫0τ∫0ξ∫0θe(ξ−σ)ALeσA∫0τ−ξe(τ−ξ−ν)AB~eν(A+B~)u~(tn)dνdσdθdξ,Dn(4) =14∫0τ∫0ξ∫0θe(ξ−σ)ALeσAξB~e(τ−ξ)(A+B~)u~(tn)dσdθdξ,Dn(5) =14∫0τξB~∫0ξ∫0θe(ξ−σ)ALeσAeξ2B~e(τ−ξ)(A+B~)u~(tn)dσdθdξ. In this way, we end up with the decomposition   δn=δ^n+Dn,Dn=Dn(1)+Dn(2)+Dn(3)+Dn(4)+Dn(5) of the local error. To prove the statement of the lemma, we have to estimate the terms $$D_n^{(i)}$$, $$i=1,\dots,5$$. The properties (2.2), (3.8b), (3.8c), (4.4) and (4.7) imply that   ‖RA−1‖ =2‖AB~‖≤2M12‖G‖, (5.6a)  ‖τLA−1‖ =‖τAAB~A−1+τB~A‖=‖(τAΨ^)(AB)Φ^A−1+Ψ^B(Φ^τA)‖≤2CinvM1M2‖G‖, (5.6b) where we use the notation (4.2). Writing   Dn(2) =−12∫0τeξ2B~∫0ξ∫0θe(ξ−σ)A(RA−1)eσA(Aeξ2B~A−1)Au~(tn+τ−ξ)dσdθdξ,Dn(3) =12τ∫0τ∫0ξ∫0θe(ξ−σ)A(τLA−1)eσA∫0τ−ξe(τ−ξ−ν)A(AB~)A−1Au~(tn+ν)dνdσdθdξ,Dn(4) =14τ∫0τξ∫0ξ∫0θe(ξ−σ)A(τLA−1)eσA(AB~)A−1Au~(tn+τ−ξ)dσdθdξ,Dn(5) =14τ∫0τξB~∫0ξ∫0θe(ξ−σ)A(τLA−1)eσA(Aeξ2B~A−1)Au~(tn+τ−ξ)dσdθdξ, the above estimates (5.6) and the bounds (2.2), (2.5), (4.5), (4.6), (4.7) and (5.1) yield   ‖Dn(2)‖ ≤16eτβ‖G‖M12K~τ3, ‖Dn(3)‖ ≤124Cinv2M13M2‖G‖2K~τ3,‖Dn(4)‖ ≤116Cinv2M13M2‖G‖2K~τ3, ‖Dn(5)‖ ≤116βeτβCinvM1M2‖G‖K~τ3. It remains to bound the integral $$D_n^{(1)}$$ defined in (5.5b). Using the representation   eξ−θ2B~[A,B~]eθ2B~−[A,B~]eξ2B~ =∫θξddν(eν−θ2B~[A,B~]eθ−ν2B~)eξ2B~dν =12∫θξeν−θ2B~[B~,[A,B~]]eξ+θ−ν2B~dν, we can write   Dn(1)=14∫0τeξ2B~∫0ξeξA∫θξeν−θ2B~[B~,[A,B~]]A−1(Aeξ+θ−ν2B~A−1)Au~(tn+τ−ξ)dνdθdξ. For the double commutator, we get from (2.2) and (4.7)   ‖[B~,[A,B~]]A−1‖ =‖2B~AB~A−1−B~2−AB~2A−1‖≤4Cinv2M14‖G‖2, which yields with (2.5) and (5.1)   ‖Dn(1)‖≤124eτβK~‖[B~,[A,B~]]A−1‖τ3≤16eτβCinv2M14‖G‖2K~τ3. The assertion now follows from $$D_n = D_n^{(1)} + D_n^{(2)} + D_n^{(3)} + D_n^{(4)} + D_n^{(5)}$$. □ We now investigate the local error further. Lemma 5.2 The dominating part $$\widehat \delta_n$$ of the local error $$\delta_n$$ defined in (5.3) satisfies   δ^n=δ^n(1)+δ^n(2),δ^n(1)=AΨ^Z1u~(tn),δ^n(2)=Z2Φ^A2u~(tn), where   ‖Z1‖≤112M1‖G‖τ3,‖Z2‖≤112CinvM1‖G‖τ3,‖AZ2‖≤112M1‖G‖τ3. Proof. The matrix $$L$$ can be written as   L =(AΨ^)ABΦ^+Ψ^BΦ^A2, since $$\widetilde{B}=\widehat{{\it{\Psi}}} B \widehat{{\it{\Phi}}}$$ with the notation (4.2). Since the filter functions $${\it{\Psi}}_S$$ and $${\it{\Phi}}$$ are independent of the integration variables, we have by (5.3)   δ^n =AΨ^Z1u~(tn)+Z2Φ^A2u~(tn), where   Z1 =12∫0τ∫0ξ∫0θe(ξ−σ)AABΦ^e(σ+τ−ξ)Adσdθdξ,Z2 =12∫0τ∫0ξ∫0θe(ξ−σ)AΨ^Be(σ+τ−ξ)Adσdθdξ. The desired bounds now follow immediately from (2.2), (2.5), (3.8a) and (4.4). □ Although the matrices $$Z_{1}$$, $$Z_2$$ and $$AZ_2$$ contain a factor of $$\tau^3$$, this is not sufficient to use a standard Lady Windermere’s fan argument to prove that the global error is of second order. The reason is the additional factors of $$A$$ in $$\widehat \delta_n^{(1)}$$ and $$\widehat \delta_n^{(2)}$$ which would yield a constant depending on $$\lVert {\it{\Omega}} \rVert$$. The proof of the following result treats the two terms defined in Lemma 5.2 separately in an appropriate way by using summation by parts. Theorem 5.3 (Global error of the averaged problem) Let the assumptions (3.8) and the finite-energy condition (4.6) be fulfilled. Then the global error of the splitting scheme (3.6) as an approximation to the solution of the averaged system (2.4) is bounded by   ‖un−u~(tn)‖≤Cτ2,0≤tn=nτ≤tend, (5.7) where $$C$$ only depends on $$C_\text{inv}$$, $$\| G \|$$, $$\widetilde{K}$$, $$M_j$$, $$j = 1,2,3$$ and $$t_{\text{end}}$$. Proof. By a telescopic identity, the global error can be written as   e~n=un−u~(tn)=(Sn−Tn)u0=∑j=0n−1Sn−j−1(S−T)Tju0=∑j=0n−1Sn−j−1δj with the local errors $$\delta_j$$ of Lemma 5.1. Lemmas 5.1 and 5.2 motivate to split the error into   e~n=e~n(1)+e~n(2)+e~n(3), (5.8) where   e~n(1)=∑j=0n−1Sn−j−1δ^j(1),e~n(2)=∑j=0n−1Sn−j−1δ^j(2),e~n(3)=∑j=0n−1Sn−j−1Dj. From the stability bound (5.2) and $$\|D_j\| \leq C \tau^3$$ by Lemma 5.1, we thus obtain   ‖e~n(3)‖=‖∑j=0n−1Sn−j−1Dj‖≤Cτ2(nτ)eβtend≤Cτ2. To bound $$\widetilde{e}_n^{(1)}$$ and $$\widetilde{e}_n^{(2)}$$, we write $$\tau A \widehat{\chi}$$ for $$\widehat{\chi} = \widehat{{\it{\Phi}}}$$ and $$\widehat{\chi} = \widehat{{\it{\Psi}}}$$ (see (4.2)) as   τAχ^=τχ^A=(eτA−I)Θ^χ=Θ^χ(eτA−I) (5.9) with   Θ^χ=θ^χ(τΩ),θ^χ(x)=12[cot⁡(x2)−11cot⁡(x2)][xχ(x)00xχ(x)] This representation follows from the identities   sin⁡(x)=cot⁡(x2)(1−cos⁡(x))andcos⁡(x)+1=cot⁡(x2)sin⁡(x). Note that by assumptions (3.8c) and (3.8d) on the filter functions $$\chi$$, we have   ‖Θ^χ‖≤C. (5.10) Next we use summation by parts. With   Ej=∑k=0j−1Sk,Fj=∑k=0j−1Tk, it holds   e~n(1)=Enδ^0(1)+∑j=0n−2En−j−1(δ^j+1(1)−δ^j(1)). Using Lemma 5.2 and (5.9) to replace $$\widehat \delta_{j}^{(1)}$$ by $$(e^{\tau A} - I) \widehat{{\it{\Theta}}}_\psi \tfrac1{\tau} Z_1 T^j u_0$$, this implies   e~n(1) =En(eτA−I)Θ^ψ1τZ1u0+∑j=0n−2En−j−1(eτA−I)Θ^ψZ11τ(T−I)Tju0. (5.11) Now we estimate the matrices in this expression. The matrices $$\widehat{{\it{\Theta}}}_\psi$$ and $$Z_1$$ can be estimated with (5.10) and Lemma 5.2, respectively. To bound $$E_{n-j-1} (e^{\tau A} - I)$$ in (5.11), we start from   Ej(eτA−I)=Ej(eτA−S)+Ej(S−I)=Ej(eτA−S)+Sj−I, and we use the stability bound (5.2) and $$\| e^{\tau A} - S \| \leq C \tau$$ by (2.5), (3.5) and (5.1) to show that   ‖Ej(eτA−I)‖≤C. To bound $$\tfrac1{\tau} (T-I) T^ju_0$$ in (5.11), we start from   1τ(T−I)Tju0 =1τ(eτA−I)Tju0+1τ(T−eτA)Tju0 =(eτA−I)(τA)−1Au~(tj)+1τ∫0τe(τ−ξ)AB~A−1Au~(tj+ξ)dξ by the variation-of-constants formula (4.1b), and we use (2.2), the finite-energy condition (4.6) and $$\|(e^{\tau A}-I) (\tau A)^{-1}\| = \| \int_0^1 {\it \,{e}}^{\sigma \tau A} {\rm d}\sigma\| \leq 1$$ by (2.5) to show that   ‖1τ(T−I)Tju0‖≤C. This yields the estimate $$\|\widetilde{e}_n^{(1)}\| \leq C \tau^2$$ for (5.11) since $$\lVert Z_1 \rVert \leq C \tau^3$$ by Lemma 5.2. Analogously, we have by Lemma 5.2 and (5.9)   e~n(2) =1τZ2Θ^ϕA(eτA−I)Fnu0+∑j=0n−2Sj1τ(S−I)Z2Θ^ϕA(eτA−I)Fn−j−1u0. Using the variation-of-constants formula (4.1b) yields   A(eτA−I)Fju0 =A(eτA−T)Fju0+A(Tj−I)u0 =−∫0τe(τ−ξ)AAB~eξ(A+B~)Fju0dξ+A(u~(tj)−u0). The finite-energy condition (4.6) shows   ‖τeξ(A+B~)Fju0‖=‖τ∑k=0j−1u~(tk+ξ)‖≤tendCinvK~. By (4.7) we thus obtain   ‖A(eτA−I)Fju0‖≤C. Moreover, by Lemma 5.2 we have   ‖1τ(S−I)Z2‖ ≤‖1τ(S−eτA)Z2‖+‖(eτA−I)(τA)−1AZ2‖≤Cτ3 since $$\|S-e^{\tau A}\|\le C\tau$$ and $$\| (e^{\tau A} - I) (\tau A)^{-1}\|\le 1$$. Together with the stability estimate (5.2) and the bound (5.10), this proves $$\| \widetilde{e}_n^{(2)} \| \leq C \tau^2$$, and thus (5.7) by (5.8). □ Remark It is also possible to prove Theorem 5.3 if (5.3) in Lemma 5.1 is replaced by the representation of the local error which has been derived in Jahnke & Lubich (2000, Theorem 2.1). While they used quadrature errors to bound the local error, we used the fundamental theorem of calculus as in (5.4). An advantage of the new representation is that $$\left\lVert G \right\rVert$$ appears only quadratically in the bounds, while it appears cubically in Jahnke & Lubich (2000). We also believe that the representation (5.3) is more suitable for a future extension of our approach to nonlinear problems. 6. Discussion and comparison Error bounds for trigonometric integrators applied to second-order oscillatory differential equations have been previously derived in Grimm & Hochbruck (2006) and, for the case of a single high frequency, in Theorems XIII.4.1 and XIII.4.2 of Hairer et al. (2006). They differ from the error bounds of the present paper (Theorem 3.2) in their statement, in the required assumptions on the filter functions and, most notably, in the technique of proof. The difference in the statement concerns the error bound for the velocities. While the error bounds of Grimm & Hochbruck (2006) and Hairer et al. (2006) for the velocities $$q'$$ are of first order in $$\tau$$, the error bound of Theorem 3.2 (in the linear case) is of second order in the rescaled velocities $${\it{\Omega}}^{-1} q'$$. For small time step sizes $$\tau \left\lVert {{\it{\Omega}}} \right\rVert=O(1)$$, the latter bound implies the former. Besides this difference in the statement, there are differences in the assumptions on the filter functions. The splitting method (3.6) with filter functions $$\phi$$ and $$\psi_S$$ coincides with the trigonometric integrator in Grimm & Hochbruck (2006) and Hairer et al. (2006) if the filter functions $$\phi$$, $$\psi$$, $$\psi_1$$ and $$\psi_0$$ in Grimm & Hochbruck (2006) and Hairer et al. (2006) are chosen in the following way:   ϕ =ϕ,ψ1=ψS,ψ=sinc⁡(⋅)ψS,ψ0=cos⁡(⋅)ψS. Now, we can compare the conditions (11)–(16) in Grimm & Hochbruck (2006), and the conditions (XIII.4.1) and (XIII.4.8) in Hairer et al. (2006) to our conditions (3.8). Our condition (3.8b) on the boundedness of the filter functions coincides identically with $$(11)$$ of Grimm & Hochbruck (2006). The conditions (13)–(16) of Grimm & Hochbruck (2006) imply that $$\chi=\phi$$ and $$\chi=\psi_S$$ are zero whenever $$\sin(\frac{\cdot}2)$$ is zero, meaning for $$x = 2k\pi$$, $$k \in \mathbb Z$$. This behaviour can also be found in our condition (3.8d). Still, an analogon to our condition (3.8c) is missing in Grimm & Hochbruck (2006) since this condition requires that $$\chi$$ decreases at least like $$x^{-1}$$ for $$x \to \infty$$. In comparison to Hairer et al. (2006), our conditions (3.8b)–(3.8d) are implied by the conditions (XIII.4.1) and (XIII.4.8) used there, which require in particular that $$|\chi(x)|$$ and $$|x\sin(\frac12 x)^{-1}\chi(x)|$$ are bounded. Note that all sets of conditions are proven to be sufficient but not necessary. But the main difference in comparison to Grimm & Hochbruck (2006) and Hairer et al. (2006) is the technique of proof. It will be interesting to see whether the technique developed in the present paper, namely a Lady Windermere’s fan that takes cancellations in the error accumulation into account, can help to gain further insight into the error behaviour of trigonometric integrators and splitting methods, in particular for nonlinear problems. Acknowledgements The authors thank Roland Schnaubelt for helpful discussions on the representation of the local error. Funding We gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through RTG 1294, CRC 1114, CRC 1173 and project GA 2073/2-1. References Blanes, S. ( 2015) Explicit symplectic RKN methods for perturbed non-autonomous oscillators: splitting, extended and exponentially fitting methods. Comput. Phys. Commun.,  193, 10– 18. Google Scholar CrossRef Search ADS   Blanes, S. & Casas, F. ( 2016) A concise introduction to geometric numerical integration . A Chapman & Hall book: Monographs and research notes in mathematics. Boca Raton: CRC Press. Google Scholar CrossRef Search ADS   Cohen, D., Gauckler, L., Hairer, E. & Lubich, C. ( 2015) Long-term analysis of numerical integrators for oscillatory Hamiltonian systems under minimal non-resonance conditions. BIT,  55, 705– 732. Google Scholar CrossRef Search ADS   Einkemmer, L. & Ostermann, A. ( 2014) Convergence analysis of Strang splitting for Vlasov-type equations. SIAM J. Numer. Anal.,  52, 140– 155. Google Scholar CrossRef Search ADS   Einkemmer, L. & Ostermann, A. ( 2015) Overcoming order reduction in diffusion-reaction splitting. Part 1: Dirichlet boundary conditions. SIAM J. Sci. Comput.,  37, A1577–A1592. Faou, E. ( 2012) Geometric numerical integration and Schrödinger equations.  Zurich Lectures in Advanced Mathematics. Zürich: European Mathematical Society (EMS), pp. xiii+138. Faou, E., Ostermann, A. & Schratz, K. ( 2015) Analysis of exponential splitting methods for inhomogeneous parabolic equations. IMA J. Numer. Anal.,  35, 161– 178. Google Scholar CrossRef Search ADS   García-Archilla, B., Sanz-Serna, J. M. & Skeel, R. D. ( 1999) Long-time-step methods for oscillatory differential equations. SIAM J. Sci. Comput.,  20, 930– 963. Google Scholar CrossRef Search ADS   Gauckler, L. ( 2015) Error analysis of trigonometric integrators for semilinear wave equations. SIAM J. Numer. Anal.,  53, 1082– 1106. Google Scholar CrossRef Search ADS   Grimm, V. & Hochbruck, M. ( 2006) Error analysis of exponential integrators for oscillatory second-order differential equations. J. Phys. A,  39, 5495– 5507. Google Scholar CrossRef Search ADS   Grimm, V. & Hochbruck, M. ( 2008) Rational approximation to trigonometric operators. BIT,  48, 215– 229. Google Scholar CrossRef Search ADS   Hairer, E., Lubich, C. & Wanner, G. ( 2006) Geometric numerical integration : structure-preserving algorithms for ordinary differential equations . Springer series in computational mathematics; 31, 2. edn. Berlin: Springer. Hansen, E. & Ostermann, A. ( 2009) Exponential splitting for unbounded operators. Math. Comp.,  78, 1485– 1496. Google Scholar CrossRef Search ADS   Hochbruck, M., Jahnke, T. & Schnaubelt, R. ( 2015) Convergence of an ADI splitting for Maxwell’s equations. Numer. Math.,  129, 535– 561. Google Scholar CrossRef Search ADS   Hochbruck, M. & Lubich, C. ( 1999) A Gautschi-type method for oscillatory second-order differential equations. Numer. Math.,  83, 403– 426. Google Scholar CrossRef Search ADS   Holden, H., Karlsen, K. H., Lie, K.-A. & Risebro, N. H. ( 2010) Splitting methods for partial differential equations with rough solutions . EMS Series of Lectures in Mathematics. Zürich: European Mathematical Society (EMS), pp. viii+226. Google Scholar CrossRef Search ADS   Holden, H., Lubich, C. & Risebro, N. H. ( 2013) Operator splitting for partial differential equations with Burgers nonlinearity. Math. Comp.,  82, 173– 185. Google Scholar CrossRef Search ADS   Hundsdorfer, W. H. & Verwer, J. G. ( 2007) Numerical solution of time dependent advection diffusion reaction equations . Springer series in computational mathematics; 33, corr. 2. print. edn. Berlin, Heidelberg: Springer. Ixaru, L. G. & Vanden Berghe, G. ( 2004) Exponential fitting . Mathematics and its Applications, vol. 568. Kluwer Academic Publishers, Dordrecht, pp. xiv+308. Google Scholar CrossRef Search ADS   Jahnke, T. & Lubich, C. ( 2000) Error bounds for exponential operator splittings. BIT,  40, 735– 744. Google Scholar CrossRef Search ADS   Koch, O. & Lubich, C. ( 2011) Variational-splitting time integration of the multi-configuration time-dependent Hartree-Fock equations in electron dynamics. IMA J. Numer. Anal.,  31, 379– 395. Google Scholar CrossRef Search ADS   Lubich, C. ( 2008) On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations. Math. Comp.,  77, 2141– 2153. Google Scholar CrossRef Search ADS   McLachlan, R. I. & Quispel, G. R. W. ( 2002) Splitting methods. Acta Numerica,  11, 341– 434. Google Scholar CrossRef Search ADS   Paternoster, B. ( 2012) Present state-of-the-art in exponential fitting. A contribution dedicated to Liviu Ixaru on his 70th birthday. Comput. Phys. Commun.,  183, 2499– 2512. Google Scholar CrossRef Search ADS   Strang, G. ( 1968) On the construction and comparison of difference schemes. SIAM J. Numer. Anal.,  5, 506– 517. Google Scholar CrossRef Search ADS   Thalhammer, M., Caliari, M. & Neuhauser, C. ( 2009) High-order time-splitting Hermite and Fourier spectral methods. J. Comput. Phys.,  228, 822– 832. Google Scholar CrossRef Search ADS   Wu, X., You, X. & Wang, B. ( 2013) Structure-preserving algorithms for oscillatory differential equations . Springer, Heidelberg; Science Press Beijing, Beijing, pp. xii+236. Google Scholar CrossRef Search ADS   Wu, X., Liu, K. & Shi, W. ( 2015) Structure-preserving algorithms for oscillatory differential equations. II . Springer, Heidelberg; Science Press Beijing, Beijing, pp. xv+298. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations