Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations

# Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations

Philippe Chartier INRIA-Rennes Bretagne Atlantique, IPSO Project    Nicolas Crouseilles INRIA-Rennes Bretagne Atlantique, IPSO Project    Mohammed Lemou CNRS and IRMAR, Université de Rennes 1 and INRIA-Rennes Bretagne Atlantique, IPSO Project    Florian Méhats IRMAR, Université de Rennes 1 and INRIA-Rennes Bretagne Atlantique, IPSO Project
###### Abstract

This work is devoted to the numerical simulation of nonlinear Schrödinger and Klein-Gordon equations. We present a general strategy to construct numerical schemes which are uniformly accurate with respect to the oscillation frequency. This is a stronger feature than the usual so called “Asymptotic preserving” property, the last being also satisfied by our scheme in the highly oscillatory limit. Our strategy enables to simulate the oscillatory problem without using any mesh or time step refinement, and the orders of our schemes are preserved uniformly in all regimes. In other words, since our numerical method is not based on the derivation and the simulation of asymptotic models, it works in the regime where the solution does not oscillate rapidly, in the highly oscillatory limit regime, and in the intermediate regime with the same order of accuracy. In the same spirit as in [5], the method is based on two main ingredients. First, we embed our problem in a suitable “two-scale” reformulation with the introduction of an additional variable. Then a link is made with classical strategies based on Chapman-Enskog expansions in kinetic theory despite the dispersive context of the targeted equations, allowing to separate the fast time scale from the slow one. Uniformly accurate (UA) schemes are eventually derived from this new formulation and their properties and performances are assessed both theoretically and numerically.

## 1 Introduction

This work is concerned with the numerical solution of highly-oscillatory differential equations in an infinite dimensional setting. Our main two applications here are the nonlinear Schrödinger equation and the nonlinear Klein-Gordon equation, although, prior to addressing them specifically, we envisage the more general situation of an abstract differential equation in a Hilbert space. To be a bit more specific, we shall consider equations of the form

 ddtuε(t)=F(t,t/ε,uε(t)),t∈[0,T],uε(0)=u0, (1.1)

where the vector field is supposed to be periodic of period with respect to the variable (we shall denote ). The parameter is supposed to have a positive real value in an interval of the form for some . However, is not necessarily vanishing and may be as well thought of as being close to : this means we can consider equation (1.1) simultaneously in different regimes, namely highly-oscillatory for small values of or smooth for larger values of , and our aim is to design a versatile numerical method, capable of handling these two extreme regimes as much as all intermediate ones.

Generally speaking, standard numerical methods for equation (1.1) exhibit errors of the form for some positive and . The user of such methods is thus forced to restrict the step-size to values less than in order to obtain some accuracy. This becomes an unacceptable constraint for vanishing values of . Whenever equation (1.1) admits a limit model, Asymptotic-Preserving (AP) schemes [10] have been designed to overcome this restriction: the methods we construct obey the corresponding requirement, i.e. they degenerate into a consistent numerical scheme for the limit model whenever tends to zero.

As favorable as this property may seem, the error behavior of an AP scheme may deteriorate for “intermediate” regimes where is neither very small nor large. The derivation of asymptotic models for (1.1) has been the subject of many works – see e.g. [3, 4, 15, 16] for time-averaging techniques and [1, 8, 14] for homogenization techniques – and a hierarchy of averaged vector fields and models at different orders of can be classically written from asymptotic expansions of the solution. However, these asymptotic models are valid only when is small enough and any numerical methods based on the direct approximation of such averaged vector fields introduce a truncation index and a corresponding incompressible error .

In sharp contrast, our strategy in this paper consists in developing numerical schemes that solve directly (1.1) for a wide range of -values with uniform accuracy. The main output of our work are numerical methods for highly oscillatory equations of type (1.1), which are uniformly accurate (UA) with respect to the parameter . These methods, as we shall demonstrate, are able to capture the various scales occurring in the system, while keeping numerical parameters (for instance ) independent of the degree of stiffness ().

The main idea underlying our strategy (see also [5]) consists in separating the two time scales naturally present in (1.1), namely the slow time and the fast time . To this aim, we embed the solution into a two-variable function while imposing that coincides with on the diagonal . Clearly, this implies that

 ddtuε(t)=∂tUε(t,t/ε)+1ε∂τUε(t,t/ε)=F(t,t/ε,Uε(t,t/ε)).

By virtue of the “separation” principle, we then consider the equation over the whole -domain, i.e.

 ∂tUε(t,τ)+1ε∂τUε(t,τ)=F(t,τ,Uε(t,τ)). (1.2)

An observation of paramount importance is that no initial condition for (1.2) is evident, since only the value is prescribed: consequently, as such, the transport equation (1.2) is not a Cauchy problem and may have many solutions. This apparent obstacle is in fact the way out to our numerical difficulties: given that for any smooth initial condition satisfying , we can recover the solution from the values of on the diagonal , the missing Cauchy condition should be regarded as an additional degree of freedom.

Now, it turns out that for some specific choice of , it is possible to prove that and its time-derivatives are bounded on uniformly w.r.t. . The point is, in this two-scale formulation (1.2) of (1.1), that stiffness is confined in the sole term . Interpreting this singularly perturbed term as a “collision” operator, we can derive the asymptotic behavior of through a Chapman-Enskog expansion (see for instance [6]) from which averaged models (first and second order) can easily be obtained. The initial datum is then chosen so as to satisfy this expansion at , a requirement compatible with . Two numerical schemes are then proposed for this augmented problem, following the strategy in [5]. In the present work, these schemes are proved to be uniformly accurate with respect to : they have respectively orders and uniformly in . These properties are assessed by numerical experiments on the nonlinear Klein-Gordon and Schrödinger equations.

This paper is organized as follows. In Section 2.1, we present the two-scale formulation in a general framework and perform in Subsection 2.3 the Chapman-Enskog expansion of . The question of the choice of the initial datum for this augmented equation (1.2) is addressed. In Section 3, a first-order numerical scheme is introduced and analyzed while a second-order one is similarly studied in Section 4. Finally, Section 5 is devoted to a series of numerical tests which confirm the theoretical properties of our schemes when applied to the Schrödinger and Klein-Gordon equations and demonstrate the relevance of our strategy.

## 2 Two-scale formulation of the oscillatory equation

In this section, we formulate and analyze the equation obtained by decoupling the slow variable and the fast one .

### 2.1 Setting of the problem

Given , we consider the following highly-oscillatory evolution problem

 ∂tuε=F(t,t/ε,uε),t≥0,ε∈]0,ε0],uε(0)=u0, (2.1)

where the unknown is a smooth map onto a Sobolev space (either or with ) and the vector-field is a smooth map, -periodic w.r.t. (). Let us emphasize that may also depend on , although we shall not reflect specifically this dependence: whenever this is the case, all bounds on and its derivatives then implicitly hold uniformly in . In order to work in Banach algebras, we require that , a condition whose necessity will become apparent for the numerical schemes.

As described in the Introduction section, we envisage as the diagonal solution of the following transport equation which constitutes our starting point:

 ∂tUε+1ε∂τUε=F(t,τ,Uε),Uε(0,τ)=Uε0(τ), (2.2)

where the unknown is now the function . The choice of the Cauchy condition is discussed below, but it is already clear that and coincide provided that .

For our purpose, we shall need that the vector field obeys the following assumption, where each derivation w.r.t. or typically costs 2 derivatives in the space variable. Indeed, for applications to nonlinear Klein-Gordon or Schrödinger equation – see (5.7) and (5.12) –, one has in mind vector fields of the form

 F(t,τ,Uε)=f(eitΔUε,eiτΔUε),

where is a smooth function.

###### Assumption (A)

For all , and , for all such that , the functional is continuous and locally bounded from to

 L(Hσ×…×Hσγ%times,Hσ−2(α+β)).

### 2.2 Bounds in Hσ of the solution of the transport equation

Similar related transport equations will occur in our analysis, with possibly other functions than and initial conditions with various regularities. A somehow preliminary result thus concerns the existence and uniqueness of the solution of the general Cauchy problem

 ∂tΦε+1ε∂τΦε=G(t,τ,Φε),Φε(0,τ)=Φε0(τ)∈Hσ, (2.3)

where (possibly) depends on and is assumed to be not identically zero.

###### Proposition 2.1

Let , let and suppose that is a locally Lipschitz continuous map from into , that it admits derivatives and which are continuous and locally bounded from into, respectively, and . If is uniformly bounded in with respect to the norm then, for any , there exists such that for all , equation (2.3) has a unique solution and we have

 ∀t∈[0,Tκ],supε∈]0,ε0]∥Φε(t,⋅)∥L∞τ(Hσ)≤κsupε∈]0,ε0]∥Φε0(⋅)∥L∞τ(Hσ). (2.4)

Moreover, has first derivatives w.r.t. both and which are functions of . If in addition, satisfies the estimate

 ∀(t,τ)∈[0,T]×T,∀v∈Hσ,∥G(t,τ,v)∥Hσ≤CG∥v∥Hσ+DG

for some positive constants and , then equation (2.3) has a unique solution in satisfying

 ∀t∈[0,T],supε∈]0,ε0]∥Φε(t,⋅)∥L∞τ(Hσ)≤(supε∈]0,ε0]∥Φε0(⋅)∥L∞τ(Hσ)+DGt)etCG.

Proof. Considering a smooth solution of (2.3) and denoting , it is easy to check that

 ∂tφε(t,τ)=∂tΦε(t,τ+t/ε)+1ε∂τΦε(t,τ+t/ε)=G(t,τ+t/ε,φε(t,τ)),

so that the smooth function , parametrized by , is then solution of the ordinary differential equation

 ∂tφε(t,τ)=G(t,τ+t/ε,φε(t,τ)),φε(0,τ)=Φε0(τ). (2.5)

According to Cauchy-Lipschitz theorem in (a Banach space), equation (2.5) has a unique maximal solution on an interval of the form (when ) or a solution on , which furthermore satisfies the following inequality

 ∥φε(t,τ)∥Hσ≤∥Φε0(τ)∥Hσ+∫t0∥G(θ,τ+θ/ε,φε(θ,τ)∥Hσdθ.

Denote

 R=supε∈]0,ε0]∥Φε0(τ)∥Hσ

and

 Mκ=sup{∥G(t,τ,u)∥Hσ,0≤t≤T,τ∈T,∥u∥Hσ≤κR}.

Now, as long as , we have

 ∥φε(t,τ)∥Hσ≤∥Φε0(τ)∥Hσ+tMκ,

so that

 Tεmax≥Tκ:=min(T,(κ−1)RMκ)>0

and estimate (2.4) holds. Now, since (resp. ) is a continuous and locally bounded function from to (resp. to ), then is the unique solution on of the linear differential equation in

 ∂tψε(t,τ) =(∂τG)(t,τ+t/ε,φε(t,τ))+∂uG(t,τ+t/ε,φε(t,τ))ψε(t,τ), ψε(0,τ) =∂τΦε0(τ)∈Hσ−2.

Hence has first derivatives in and in . Finally, since , we have

 ∂τΦε(t,τ)=∂τφε(t,τ−t/ε) and ∂tΦε(t,τ)=∂tφε(t,τ−t/ε)−1ε∂τφε(t,τ−t/ε)

so that has also first derivatives in . Finally, the proof of the subsequent assertions in Proposition 2.1 can be done in the same way, the last estimate being a consequence of the Gronwall lemma.

###### Remark 2.2

From previous formulae, it appears that exists in but is not necessarily uniformly bounded in . In order to get a solution with uniformly bounded first derivatives, we have to consider an appropriate -dependent initial condition . In the next two subsections, we shall consider a formal expansion of in so as to determine how this initial condition should be prescribed.

### 2.3 A formal Chapman-Enskog expansion

In this subsection, we analyze formally the behavior of (2.2) in the limit under the assumption that its solution has uniformly bounded (in ) derivatives up to order . Following [5], we thus consider the linear operator , defined for all periodic (regular) function by

 Lh=∂τh.

This operator is skew-adjoint with respect to the scalar product and its kernel is the set of constant functions. The -projector on this kernel is the averaging operator

 Πh:=1P∫P0h(τ)dτ

which obviously satisfies . On the set of functions with vanishing average, is invertible with inverse defined by

 (L−1h)(τ)=(\rm I−Π)∫τ0h(θ)dθ.

In order to alleviate notations, we further introduce which operates on the set of periodic functions onto the set of zero-average periodic functions.

The Chapman-Enskog expansion (see for instance [6]) consists in writing the solution in the form

 Uε(t,τ)=U––ε(t)+hε(t,τ), (2.6)

where

 U––ε(t)=Π(Uε(t,τ)),Πhε=0,

and then, under some regularity assumptions on with respect to and , one seeks the correction as an expansion in powers of :

 hε(t,τ)=εh1(t,τ,U––ε(t))+ε2h2(t,τ,U––ε(t))+…. (2.7)

Inserting the decomposition (2.6) into (2.2) leads to

 ∂tU––ε+∂thε+1εLhε=F(t,τ,U––ε+hε). (2.8)

Projecting on the kernel of and taking into account that , we obtain

 ∂tU––ε=Π(F(t,τ,U––ε+hε)), (2.9)

and then subtracting from (2.8)

 ∂thε+1εLhε=(\rm I% −Π)(F(t,τ,U––ε+hε)). (2.10)

Since belongs to the range of , we get

 hε=εA(F(t,τ,U––ε+hε))−εL−1(∂thε). (2.11)

Therefore, provided and its first time derivative are uniformly bounded w.r.t. , we first deduce from this last equation that . Now if we additionally assume that the second and third time derivatives are uniformly bounded w.r.t. , then, by a simple induction on (2.11), we get

 hε(t,τ)=εh1(t,τ,U––ε)+ε2h2(t,τ,U––ε)+O(ε3),

with and defined by

 h1(t,τ,U)= AF(t,τ,U), (2.12) h2(t,τ,U)= A∂uF(t,τ,U)AF(t,τ,U)−A2(∂uF(t,τ,U)ΠF(t,τ,U)+∂tF(t,τ,U)). (2.13)

Inserting these corrections into equation (2.9) yields the first and second order averaged models

 ∂tU––ε =ΠF(t,τ,U––ε)+O(ε), ∂tU––ε =ΠF(t,τ,U––ε)+εΠ(∂uF(t,τ,U––ε)AF(t,τ,U––ε))+O(ε2).

Anticipating on next sections, let us now briefly address the crucial issue of the initial condition for (2.2). According to the above calculations, one expects to get a smooth solution of (2.2) if the initial condition follows the same expansion as above, i.e.

 Uε0(τ) =U––ε0+εAF0(τ,U––ε0) (2.14) +ε2(A∂uF0(τ,U––ε0)AF0(τ,U––ε0)−A2∂uF0(τ,U––ε0)ΠF0(τ,U––ε0)−A2(∂tF)0(τ,U––ε0)))

where we have denoted by a subindex the evaluation of functions at and where is chosen so as to be compatible with which is the initial condition for the original problem (2.1). Starting from and inserting successively higher-order terms in the previous equation, we can obtain the expression of and then of . For instance, we have

 U––ε0=u0+εΠ∫τ0(\rm I% −Π)F0(θ,u0)dθ+O(ε2),

so that

 Uε0(τ)=u0+ε∫τ0(\rm I−Π)F0(θ,u0)dθ+O(ε2)

which provides an initial condition for our first order numerical scheme (see Subsection 3.3). The explicit computation of second order terms is postponed to Subsection 4.3.

### 2.4 Estimates of time derivatives

In this subsection, we indeed prove that the initial condition (2.14) ensures that time derivatives of up to order are uniformly bounded in . In the sequel, the following functional space will be useful:

 Xσ=⋂0≤2ℓ<σ−d/2Cℓ(T;Hσ−2ℓ). (2.15)
###### Proposition 2.3

Suppose that satisfies Assumption (A) and let and . Consider the following initial condition

 ∀τ∈T,Uε(0,τ)=Uε0(τ)=U––ε0+(εh1+ε2h2)(0,τ,U––ε0)+ε3rε(τ), (2.16)

where is assumed to be uniformly bounded in , where and are given by (2.12) and (2.13), and where the remainder term is assumed to be bounded in uniformly in . Then the following holds:

1. is uniformly bounded in and there exists such that, for all , equation (2.2), subject to the initial condition (2.16), has a unique solution , which satisfies the uniform bound

 ∀t∈[0,Tκ],supε∈]0,ε0]∥Uε(t)∥L∞τ(Hs)≤κsupε∈]0,ε0]∥Uε0∥L∞τ(Hs). (2.17)
2. Moreover, for any for which (2.17) holds, the solution satisfies the following estimates

 ∀t∈[0,Tκ],supε∈]0,ε0]∥∂αt∂βτUε(t)∥L∞τ(Hs−2(α+β))≤C,α=0,1,2,3,β=0,1, (2.18)

for some constant .

Proof. We prove this proposition in several steps.

Existence of and uniform bound. Let us first estimate the initial condition defined by (2.16). From (2.12), (2.13), one gets

 Uε0=U––ε0+εAF0+ε2A∂uF0AF0−ε2A2∂uF0ΠF0−ε2A2(∂tF)0+ε3rε, (2.19)

where, for conciseness, we have further omitted the dependence111In the sequel, we explicitly mention the dependence while stands for and stands for . of in and . We notice that, by Assumption (A), we have

 F0,∂uF0AF0,∂uF0ΠF0∈Xs+2, and (∂tF)0∈Xs,

with norms uniformly bounded w.r.t. . Hence, observing that and are bounded operators on , for all , one deduces that belongs to and is uniformly bounded w.r.t. . Hence, according to Proposition 2.1, exists on an interval independent of , and satisfies

 ∀0≤t≤Tκ,supε∈]0,ε0]∥Uε(t,⋅)∥L∞τ(Hs)≤κsupε∈]0,ε0]∥Uε0(⋅)∥L∞τ(Hs).

Furthermore, derivatives of w.r.t. and exist and are functions with values in .

Estimate of the first derivative in . The first derivative satisfies the equation

 ∂tVε+1ε∂τVε=∂uF(t,τ,Uε)Vε+∂tF(t,τ,Uε) (2.20)

with initial condition

 Vε0=F0(Uε0)−1εLUε0.

From (2.19) and , , we obtain

 Vε0=F0(Uε0)−F0+ΠF0−ε(\rm I−Π)∂uF0AF0+εA∂uF0ΠF0+εA(∂tF)0−ε2Lrε.

Taylor-Lagrange expansions with integral remainder at order one and two give222The notation is used here for terms uniformly bounded in with the appropriate -norm.

 F0(Uε0)−F0=∫10∂uF0(U––ε0+μ(Uε0−U––ε0))dμ(Uε0−U––ε0)=OXs(ε) =∂uF0(U––ε0)(Uε0−U––ε0)+∫10(1−μ)∂2uF0(U––ε0+μ(Uε0−U––ε0))dμ(Uε0−U––ε0,Uε0−U––ε0) =ε∂uF0AF0+OXs(ε2),

where we used that is uniformly bounded in . Therefore333Notice that is continuously embedded in .

 Vε0 =ΠF0+OXs−2(ε) (2.21) =ΠF0+ε∂uF0AF0−ε(\rm I−Π)∂uF0AF0+εA∂uF0ΠF0+εA(∂tF)0+OXs−2(ε2) =ΠF0+εΠ∂uF0AF0+εA∂uF0ΠF0+εA(∂tF)0+OXs−2(ε2). (2.22)

In particular, and is uniformly bounded in w.r.t. . According to the second part of Proposition 2.1 (for ) with

 G(t,τ,V)=∂uF(t,τ,Uε(t,τ))V+∂tF(t,τ,Uε(t,τ))

which is a map from into , we thus have an estimate of the form

 ∀t∈[0,Tκ],supε∈]0,ε0]∥Vε(t,⋅)∥L∞τ(Hs−2)≤(supε∈]0,ε0]∥Vε0(τ)∥L∞τ(Hs−2)+DGTκ)eTκCG.

Estimate of the second derivative in . We proceed in an analogous way for by considering

 ∂tWε+1ε∂τWε= ∂2uF(t,τ,Uε)(Vε,Vε)+2∂t∂uF(t,τ,Uε)Vε +∂2tF(t,τ,Uε)+∂uF(t,τ,Uε)Wε. (2.23)

The initial condition for can be obtained from (2.20) at and (2.22)

 Wε0 =−1εLVε0+∂uF0(Uε0)Vε0+(∂tF)0(Uε0) =−1εL(εA∂uF0ΠF0+εA(∂tF)0)+∂uF0(Uε0)Vε0+(∂tF)0(Uε0)+OXs−4(ε) =−(\rm I−Π)∂uF0ΠF0−(\rm I−Π)(∂tF)0+∂uF0ΠF0+(∂tF)0+OXs−4(ε) =Π∂uF0ΠF0+Π(∂tF)0+OXs−4(ε), (2.24)

which is uniformly bounded in the -norm w.r.t. both and . By Proposition 2.1 applied to with , one gets that is uniformly bounded in .

Estimate of the third derivative in . Finally, we derive the equation for , which reads

 ∂tYε+1ε∂τYε =∂3uF(t,τ,Uε)(Vε,Vε,Vε)+3∂t∂2uF(t,τ,Uε)(Vε,Vε) +3∂2uF(t,τ,Uε)(Vε,Wε)+3∂2t∂uF(t,τ,Uε)Vε +3∂t∂uF(t,τ,Uε)Wε+∂3tF(t,τ,Uε)+∂uF(t,τ,Uε)Yε. (2.25)

We then extract from (2.4):

 Yε0=−1εLWε0+∂2uF0(Uε0)(Vε0,Vε0)+2(∂t∂uF)0(