Uniformly accurate time-splitting methods for the semiclassical Schrödinger equation Part 1 : Construction of the schemes and simulations
UUNIFORMLY ACCURATE TIME-SPLITTING METHODS FORTHE SEMICLASSICAL SCHRÖDINGER EQUATIONPART 1 : CONSTRUCTION OF THE SCHEMES ANDSIMULATIONS
PHILIPPE CHARTIER, LOÏC LE TREUST, AND FLORIAN MÉHATS
Abstract.
This article is devoted to the construction of new numerical methodsfor the semiclassical Schrödinger equation. A phase-amplitude reformulation ofthe equation is described where the Planck constant ε is not a singular parameter.This allows to build splitting schemes whose accuracy is spectral in space, of upto fourth order in time, and independent of ε before the caustics. The second-order method additionally preserves the L -norm of the solution just as the exactflow does. In this first part of the paper, we introduce the basic splitting schemein the nonlinear case, reveal our strategy for constructing higher-order methods,and illustrate their properties with simulations. In the second part, we shallprove a uniform convergence result for the first-order splitting scheme applied tothe linear Schrödinger equation with a potential. Introduction
In this paper, we are concerned with the numerical approximation of the solution Ψ ε : R + × R d → C , d ≥ , of the nonlinear Schrödinger (NLS) equation iε∂ t Ψ ε = − ε ε + | Ψ ε | Ψ ε , (1.1)where ε > is the so-called semiclassical parameter. The initial datum is assumedto be of the form Ψ ε (0 , · ) = A ( · ) e iS ( · ) /ε such that (cid:107) A (cid:107) L ( R d ) = 1 . (1.2)Note that the L -norm, the energy and the momentum of Ψ ε (0 , · ) , namelyMass: (cid:107) Ψ ε ( t, · ) (cid:107) L ( R d ) , (1.3)Energy: (cid:90) R d (cid:0) ε |∇ Ψ ε ( t, x ) | + | Ψ ε ( t, x ) | (cid:1) dx, (1.4)Momentum: ε Im (cid:90) R d Ψ ε ( t, x ) ∇ Ψ ε ( t, x ) dx, (1.5)are all preserved by the flow of (1.1), whenever Ψ ε (0 , · ) ∈ H ( R d ) .Owing to its numerous occurrences in a vast number of domains of applicationsin physics, equation (1.1) has been widely studied: existence and uniqueness results Mathematics Subject Classification.
Key words and phrases. nonlinear Schrödinger equation, semiclassical limit, numerical simula-tion, uniformly accurate, Madelung transform, splitting schemes, eikonal equation.This work was supported by the ANR-FWF Project Lodiquas ANR-11-IS01-0003 and the ANR-10-BLAN-0101 Grant (L.L.T.) and by the ANR project Moonrise ANR-14-CE23-0007-01. a r X i v : . [ m a t h . A P ] M a y P. CHARTIER, L. LE TREUST, AND F. MÉHATS can be found for instance in [15]. In particular, it has been frequently used for thedescription of Bose-Einstein condensates, as well as for the study of the propagationof laser beams (see [30] for a detailed presentation of the physical context). In thesemiclassical regime where the rescaled Planck constant ε is small, its asymptoticstudy allows for an appropriate description of the observables of Ψ ε through thelaws of hydrodynamics. We refer to [10] for a detailed presentation of the semiclas-sical analysis of the NLS equation and to [23] for a review of both theoretical andnumerical issues.1.1. Motivation.
Generally speaking, numerical methods for equation (1.1) ex-hibit an error of size ∆ t p /ε r + ∆ x q /ε s , where ∆ t and ∆ x are the time and spacesteps and p, q, r, s strictly positive numbers. For time-splitting methods in the linearcase for instance, the error on the wave function behaves like ∆ x/ε + ∆ t p /ε [2, 18],while in the nonlinear framework, numerical experiments [3, 13, 19] suggest largervalues of r and s after the appearance of caustics . Even if we content ourselveswith observables in the linear case , the error of a splitting method of Bao, Jin, andMarkowich [2] grows like ∆ x/ε + ∆ t p . Now, achieving a fixed accuracy for varyingvalues of ε requires to keep both ratios ∆ t/ε r/p and ∆ x/ε s/q constant, and becomes prohibitively costly when ε → . Our aim, in this article, is thus to develop new numerical schemes that are Uniformly Accurate (UA) w.r.t. ε , i.e. whose accuracydoes not deteriorate for vanishing ε . In other words, schemes for which r, s = 0 . Thisseems highly desirable as all available methods with the exception of [6], namelyfinite difference methods [1, 17, 24, 32], splitting methods [5, 18, 19, 25, 28, 31],relaxation schemes [4] and symplectic methods [29] fail to be UA.It is the belief of the authors that, prior to the construction of UA-schemes, it isnecessary to reformulate (1.1) as in [6] and we now describe how this can be done.1.2. Reformulation of the problem.
In the spirit of the
Wentzel-Kramers-Brillouin(WKB) techniques, we decompose Ψ ε as the product of a slowly varying amplitudeand a fast oscillating factor Ψ ε ( t, · ) = A ε ( t, · ) e iS ε ( t, · ) /ε . (1.7)From this point onwards, various choices are possible, depending on whether A ε iscomplex or not : taking A ε ∈ C leads to the following system [21] ∂ t S ε + |∇ S ε | | A ε | = 0 , (1.8a) ∂ t A ε + ∇ S ε · ∇ A ε + A ε S ε = iε ∆ A ε (1.8b) We refer to [19] for a study of local error estimates for the semiclassical NLS equation. These authors performed extensive numerical tests in both linear and nonlinear cases [2, 3]. Considering the WKB-ansatz (1.7) transforms the invariants (1.3) into respectively (cid:107) A ε (cid:107) L ( R d ) , (cid:90) R d (cid:0) | ε ∇ A ε + iA ε ∇ S ε | + | A ε | (cid:1) dx and Im (cid:90) R d A ε ( ε ∇ A ε + iA ε ∇ S ε ) dx. (1.6) The Madelung transform [26] relates the semiclassical limit of (1.1) to hydrodynamic equations Ψ ε ( t, · ) = (cid:112) ρ ε ( t, · ) e iS ε ( t, · ) /ε and amounts to choosing A ε ∈ R + . However, this formulation leads to both analytical andnumerical difficulties in the presence of vacuum, i.e. whenever ρ ε vanishes [12, 16]. with S ε (0 , · ) = S ( · ) and A ε (0 , · ) = A ( · ) , whose analysis leans on symmetrizablequasilinear hyperbolic systems (the first existence and uniqueness result has beenobtained by Grenier [21]). Under appropriate smoothness assumptions, ( A ε , S ε ) ∈ C × R converges when ε → to the solution ( A , S ) of ∂ t S + |∇ S | | A | = 0 , (1.9a) ∂ t A + ∇ S · ∇ A + A S = 0 . (1.9b)Notice that ( ρ, v ) = ( | A | , ∇ S ) is then solution of the compressible Euler system ∂ t v + v · ∇ v + ∇ ρ = 0 , (1.10a) ∂ t ρ + div( ρv ) = 0 . (1.10b)Now, an important drawback of (1.8) stems from the formation of caustics in finitetime [10]: the solution of (1.8) may indeed cease to be smooth even though Ψ ε isglobally well-defined for ε > . In order to obtain global existence for ε > (at leastin the 1D-case), Besse, Carles and Méhats [6] suggested an alternative formulationby introducing an asymptotically-vanishing viscosity term in the eikonal equation(1.8a). Therein, system (1.8) is replaced by ∂ t S ε + |∇ S ε | | A ε | = ε ∆ S ε , (1.11a) ∂ t A ε + ∇ S ε · ∇ A ε + A ε S ε = iε ∆ A ε − iεA ε ∆ S ε (1.11b)where S ε (0 , x ) = S ( x ) , A ε (0 , x ) = A ( x ) and where x ∈ R d . Let us emphasize thatboth (1.8) and (1.11) are equivalent to (1.1) in the following sense: as long as thesolution ( S ε , A ε ) of (1.8) (resp. (1.11)) is smooth, the function Ψ ε defined by (1.7)solves (1.1). The following existence and uniqueness result given from [21, 6] is thusperfectly satisfactory for our purpose. Theorem 1.1 (Grenier, Besse-Carles-Méhats ) . Let ε max > . Assume that ( S , A ) belongs to H s +2 ( R d ) × H s ( R d ) where s > d/ . Then, there exist T > , inde-pendent of ε ∈ (0 , ε max ] and a unique solution ( S ε , A ε ) ∈ C ([0 , T ]; H s +2 ( R d ) × H s ( R d )) of system of equations (1.8) , resp. (1.11) . Moreover, ( S ε , A ε ) is bounded in C ([0 , T ]; H s +2 ( R d ) × H s ( R d )) uniformly in ε and Ψ ε ( t, · ) = A ε ( t, · ) e iS ε ( t, · ) /ε is the unique solution of (1.1) withthe initial datum (1.2) . The main advantage of the WKB reformulation (1.11) over (1.1) is apparent: thesemiclassical parameter ε does not give rise to singular perturbations . Hence, it Besse et al. [6] also proved that in the one-dimensional case, the solution of (1.11) is global intime under the assumptions of Theorem 1.1 and that ( S ε , A ε ) is globally uniformly bounded. The Cole-Hopf transformation [20, Section . . ] w ε = exp (cid:18) − S ε ε (cid:19) − (1.12)transforms (1.11a) into ∂ t w ε − | A ε | ε ( w ε + 1) = ε ∆ w ε for which the regularizing effect of theviscosity term becomes arguably more apparent. P. CHARTIER, L. LE TREUST, AND F. MÉHATS constitutes a good basis for the development of UA schemes (at least prior to theappearance of caustics), as witnessed by the methods introduced later in this paper.1.3.
Content of the paper.
First and only (up to our knowledge) UA schemes arebased on the formulation (1.11) introduced in [6]. Nevertheless, these schemes arestill subject to CFL stability conditions and are of low order in time and space. Inthis paper, we consider, in lieu of finite differences as in [6], time-splitting methods,for they enjoy the following favorable features:(i) they do not suffer from stability restrictions on the time step;(ii) they are easy to implement;(iii) they can be composed to attain high-order of convergence in time while re-maining spectrally convergent in space.The third point requires further explanation: standard compositions (for ordershigher or equal to ) are inappropriate here, as they necessarily involve negativecoefficients [7]. At the same time, parabolic terms prevent the corresponding equa-tions from being solved backward in time. Negative coefficients are thus forbidden.To bypass this apparent contradiction, it has become customary to resort to com-plex coefficients with positive real parts [8, 9, 14, 22]. Now, the situation is hererendered even more involved by the presence of Schrödinger terms, which are in-compatible with complex coefficients. The way out will consist in distributing thereal and complex coefficients to the different parts of the splitting in a clever way.The strategy will be discussed thoroughly in Section 3.The outline of the remaining of the paper is as follows. In Section 2, we willintroduce first and second order (in time) splitting schemes, preserving exactly the L -norm, and in Section 4, we will present extensive numerical experiments com-paring our methods to the Strang splitting method studied in [2, 3]. The analysis ofthe uniform order of convergence of the splitting scheme is postponed to Part 2 ofthis work and, given the technicality of the proofs, will be envisaged only for orderone and for the linear Schrödinger equation with a potential. Second-order numerical scheme
The UA scheme that we now introduce is built upon splitting techniques (see forinstance [27] for a general exposition). The first building block is the resolution ofthe eikonal equation. In the sequel, h denotes the time step.2.1. The eikonal equation.
In this section, we introduce two numerical schemesfor equation (5.1).2.1.1.
A semi-Lagrangian scheme for the eikonal equation.
Let us rewrite equation(5.1) ∂ t S + H ( ∇ S ) = 0 , where H ( p ) = | p | / . Thanks to the method of characteristics, we get for x ∈ R d and h ≥ small enough that S ( h, x ( h, x )) = S (0 , x ) + h | p ( h, x ) | where ( x, p ) solves the Hamilton equation, (see [20, Section 3.2]) ∂ h x ( h, x ) = ∂ p H ( p ( h, x )) = p ( h, x ) ,∂ h p ( h, x ) = − ∂ x H ( p ( h, x )) = 0 with x (0 , x ) = x and p (0 , x ) = ∇ S (0 , x ) , i.e x ( h, x ) = x + h ∇ S (0 , x ) and p ( h, x ) = ∇ S (0 , x ) . Hence, defining implicitly for any x ∈ R d , the function y ( h, x ) = Γ( x, h, y ( h, x )) where Γ( x, h, y ( h, x )) = x − h ∇ S (0 , y ( h, x )) , we get S ( h, x ) = S (0 , y ( h, x )) + h |∇ S (0 , y ( h, x )) | . Let us define y ( h, x ) = Γ( x, h, x ) ,y ( h, x ) = Γ( x, h, Γ( x, h, x )) . and for i = 1 , , S i ( h, x ) = S (0 , x − h ∇ S (0 , y i ( h, x ))) + h |∇ S (0 , y i ( h, x )) | . After long but straightforward calculations, we can prove that S and S definemethods of order and in time respectively that do not suffer from CFL restric-tions. Remark 2.1.
The authors think that y k ( h, x ) = Γ( x, h, y k − ( h, x )) gives rise to amethod of order k in time. In the case where the space is also discretized, the functions S (0 , · ) and ∇ S (0 , · ) are interpolated thanks to the discrete Fourier transform.Most of the numerical schemes introduced for equation (5.1) (ENO and WENOtype or semi-lagrangian methods) are designed to deal with low regularity solutions.In general, their precisions is of low order in time and space and the procedure toget high order schemes is very involved. The finite difference methods even sufferfrom CFL restrictions.The methods we introduced are of spectral precision in space and can be of anyorder in time. This can be achieved in our case since we deal with regular solutions.2.1.2. A splitting scheme for the eikonal equation.
Let us remind that the Cole-Hopftransformation w ε = exp (cid:18) − S ε ε (cid:19) − ensures that solving ∂ t S ε + |∇ S ε | ε ∆ S ε , S ε (0 , · ) = S ( · ) , (2.1)is equivalent to solving the heat equation ∂ t w ε = ε ∆ w ε , w ε (0 , · ) = exp (cid:18) − S ( · )2 ε (cid:19) − , (2.2)(see Remark 1.12). P. CHARTIER, L. LE TREUST, AND F. MÉHATS
Due to the limitations of the floating-point representation, it is completely inac-curate to solve (2.1) using (2.2) for small values of ε . To overcome this difficulty,we split the nonviscous eikonal equation (5.1) into two subequations, at each timestep; the parabolic term ε ∆ S will be dealt with separately. The key idea is to allow S to be a complex-valued function despite the fact that the solutions of equations(1.8a), (1.11a), (2.1) and (5.1) take their values in R . First flow: let us define φ h as the exact flow at time h ∈ R of equation ∂ t S + ∇ S · ∇ S − i ∆ S = 0 (2.3)where ∇ S · ∇ S = (cid:80) dk =1 ( ∂ k S ) . This equation can be solved thanks to the followingmodified Cole-Hopf transformation w = exp (cid:18) iS (cid:19) − leading to i∂ t w = − ∆ w, w (0 , · ) = exp (cid:18) iS (0 , · )2 (cid:19) − , (2.4)solved in the Fourier space and S ( h, · ) − S (0 , · ) = − i log (cid:18) w ( h, · ) − w (0 , · ) w (0 , · ) + 1 (cid:19) . Remark 2.2.
This formula is well-defined whenever (cid:13)(cid:13)(cid:13)(cid:13) w ( h, · ) − w (0 , · ) w (0 , · ) + 1 (cid:13)(cid:13)(cid:13)(cid:13) L ∞ < , a condition ensured as soon as h is small enough.Second flow: φ h is the exact flow at time h ∈ R of the free Schrödinger equation ∂ t S + i ∆ S = 0 . (2.5)We are now able to define time-splitting schemes for equation (5.1). In particular,the Lie-Trotter splitting formula , Re φ h ◦ φ h (2.6)gives us an approximation of first-order in time of the solution of (5.1) (which isreal-valued) while the Strang splitting formula Re φ h/ ◦ φ h ◦ φ h/ (2.7)provides a second-order method.2.2. Numerical schemes for system of equations (1.11) . We are now in posi-tion to introduce our numerical schemes. To this aim, we split system (1.11) ∂ t S ε + |∇ S ε | | A ε | = ε ∆ S ε ,∂ t A ε + ∇ S ε · ∇ A ε + A ε S ε = iε ∆ A ε − iεA ε ∆ S ε , into four subsystems. First flow:
Let us define ϕ h as the approximate flow at time h ∈ R of the systemof equations: ∂ t S + |∇ S | , (2.8a) ∂ t A + ∇ S · ∇ A + A S = i ∆ A . (2.8b)The eikonal equation (2.8a) is solved according to Sec. 2.1. Equation (2.8b) is dealtwith by noticing that w = A exp ( iS ) satisfies the free Schrödinger equation i∂ t w = −
12 ∆ w. Second flow: ϕ h is the exact flow at time h ∈ R of ∂ t S = 0 , (2.9a) ∂ t A = i ( ε −
1) ∆ A (2.9b)solved in the Fourier space. Third flow: ϕ h is the exact flow at time h ∈ R of ∂ t S = −| A | , (2.10a) ∂ t A = 0 . (2.10b) Fourth flow: ϕ h is the exact flow at time h ∈ R + of ∂ t S = ε ∆ S, (2.11a) ∂ t A = − iεA ∆ S. (2.11b)Equation (2.11a) is solved in the Fourier space and the solution of (2.11b) is A ( h, · ) = exp (cid:0) − iε − ( S ( h, · ) − S (0 , · )) (cid:1) A (0 , · ) . Remark that ϕ h gathers the terms of (1.11) which are not in (1.8) and can thus beviewed as a regularizing flow. Remark 2.3.
Let us stress that ϕ h is not defined for h such that Re h < . As amatter of fact, the propagator e z ∆ is well-defined, in the distributional sense, if andonly if Re( z ) ≥ . We consider now the following methods for (1.11) ϕ h ◦ ϕ h ◦ ϕ h ◦ ϕ h (2.12)and ϕ h/ ◦ ϕ h/ ◦ ϕ h/ ◦ ϕ h ◦ ϕ h/ ◦ ϕ h/ ◦ ϕ h/ . (2.13) Remark 2.4.
It is worth mentioning that both schemes preserve exactly the L norm of A since all ϕ h , ϕ h , ϕ h and ϕ h do so. P. CHARTIER, L. LE TREUST, AND F. MÉHATS Fourth-order numerical scheme
The splitting of (1.11) into the four flows (2.8), (2.9), (2.10), (2.11) proposed inthe previous section is incompatible with splitting methods of order higher than with real-valued coefficients. Indeed, it is known that such methods involve at leastone negative time step for each part of the splitting (see for instance [7]). Therefore,we cannot built such a scheme for (1.11) because of its time irreversibility.To circumvent this difficulty, it is possible to use splitting methods with complexcoefficients [14, 8, 9, 22]. Let us point out the main restrictions on the coefficientsin order for the methods to be well-defined. For obvious consistency reasons, if aflow ϕ h is used with complex coefficients, then both ϕ αh and ϕ βh with Re α > and Im β < will appear. Hence, the flows ϕ h of (2.8) and ϕ h of (2.9) containingSchrödinger type terms have to be integrated with only real-valued coefficients,otherwise some parts of the splitting would be ill-posed. Moreover, ϕ h of (2.11)contains parabolic terms and it should be used with coefficients with nonnegativereal part. In this section, we introduce a four-flow complex splitting method takinginto account all these constraints.The main remaining problem originates from the non-analytic character of thenonlinearity appearing in the flow ϕ h of (2.10). To overcome it, we split the realand imaginary parts of A ε as in [8].Although we content ourselves in the sequel with a fourth-order scheme, let usemphasize that the strategy adopted here is amenable to higher orders.3.1. The new splitting scheme.
We commence from the original system (1.11) ∂ t S ε + |∇ S ε | | A ε | = ε ∆ S ε ,∂ t A ε + ∇ S ε · ∇ A ε + A ε S ε = iε ∆ A ε − iεA ε ∆ S ε and rewrite it (following the steps exposed in [8]) in term of the unknowns A ε =Re A ε , A ε = Im A ε and S ε : ∂ t S ε + |∇ S ε | A ε ) + ( A ε ) = ε ∆ S ε , (3.1a) (cid:18) ∂ t + ∇ S ε · ∇ + ∆ S ε (cid:19) (cid:18) A ε A ε (cid:19) = − iεσ (cid:18) ∆2 − ∆ S ε (cid:19) (cid:18) A ε A ε (cid:19) (3.1b)in order for the nonlinearity ( A ε ) + ( A ε ) to be an analytic function of A ε and A ε .The matrix σ is here the second Pauli matrix σ = (cid:18) − ii (cid:19) , so that P σ P = σ with σ = (cid:18) − (cid:19) and P = 1 √ (cid:18) − ii − (cid:19) = 1 √ σ + σ ) . (3.2)Let V ε = (cid:18) v ε v ε (cid:19) = P (cid:18) A ε A ε (cid:19) . System (3.1) becomes ∂ t S ε + |∇ S ε | − iv ε v ε = ε ∆ S ε , (3.3a) (cid:18) ∂ t + ∇ S ε · ∇ + ∆ S ε (cid:19) V ε = − iεσ (cid:18) ∆2 − ∆ S ε (cid:19) V ε . (3.3b)We are now in position to define a four-flow splitting which is compatible withcomplex coefficients. First flow:
Let us define (cid:101) ϕ h as the approximate flow at time h ∈ R of: ∂ t S + ∇ S · ∇ S , (3.4a) ∂ t V + ∇ S · ∇ V + ∆ S V = i ∆ V . (3.4b)To solve (3.4a), we use either the semi-lagrangian method of order of Section 2.1.1or the following fourth-order time-splitting from [33] φ α h ◦ φ α h ◦ φ α h ◦ φ α h ◦ φ α h ◦ φ α h ◦ φ α h , (3.5)with coefficients α , . . . α defined by α = α = 12(2 − / ) , α = α = 0 . − α ,α = α = 12 − / , α = 1 − α (3.6)and where the numerical flows φ h and φ h are those introduced in Sec. 2.1.2. Tosolve (3.4b), we proceed as for (2.8b). Second flow: (cid:101) ϕ h is the exact flow at time h ∈ R of ∂ t S = 0 , (3.7a) ∂ t V = − i ( εσ + σ ) ∆ V (3.7b)solved in the Fourier space. Here σ denotes the identity matrix. Third flow: (cid:101) ϕ h is the exact flow at time h ∈ R of ∂ t S = 2 iv v , (3.8a) ∂ t V = 0 . (3.8b) Fourth flow: (cid:101) ϕ h is the exact flow at time h ∈ { z ∈ C , Re z ≥ } of ∂ t S = ε ∆ S, (3.9a) ∂ t V = iεσ V ∆ S. (3.9b)We will see below that in the complex splitting method that we use, the coefficientsrelated to (cid:101) ϕ h are complex so that S , v and v are complex-valued functions. Splitting scheme of fourth-order for system (1.11) . We define below acomplex splitting method whose coefficients related to (cid:101) ϕ h have positive real partwhereas those associated with (cid:101) ϕ h and (cid:101) ϕ h are real-valued.The simplest way to get a fourth-order time-splitting scheme for four flows is tocompose several times a fourth-order time-splitting for two flows: using the sametime-splitting scheme as in (3.5), we define the following fourth-order schemes (for h ∈ R ): (cid:101) ϕ h = (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h for the system of equations ∂ t S + ∇ S · ∇ S ,∂ t V + ∇ S · ∇ V + ∆ S V = − iεσ ∆ V . and (cid:101) ϕ h = (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h ◦ (cid:101) ϕ α h , for the system of equations ∂ t S + ∇ S · ∇ S − iv v = 0 ,∂ t V + ∇ S · ∇ V + ∆ S V = − iεσ ∆ V . The coefficients α , . . . , α are defined by (3.6). Since (cid:101) ϕ h is not reversible, we cannotuse the scheme (3.5) anymore to define our four-flow method, given that α , α and α are negative. To avoid this problem, we use a complex splitting method of Blanes et al. [8]: (cid:101) ϕ h = (cid:101) P (cid:0) (cid:101) ϕ β h ◦ (cid:101) ϕ β h ◦ (cid:101) ϕ β h ◦ (cid:101) ϕ β h ◦ (cid:101) ϕ β h ◦ (cid:101) ϕ β h ◦ (cid:101) ϕ β h ◦ (cid:101) ϕ β h ◦ (cid:101) ϕ β h (cid:1) (cid:101) P (3.10)where (cid:101) P = (cid:18) P (cid:19) ,P is defined in (3.2) and β = β = 0 . − . i,β = β = 0 . ,β = β = 0 . . i,β = β = 0 . − β = 0 . ,β = 1 − β − β = 0 . − . i. Observe that all the coefficients β , β , β , β and β for the irreversible flow (cid:101) ϕ h have a positive real part and all the coefficients β , β , β and β for the flow (cid:101) ϕ h containing all the Schrödinger terms are real-valued.Following [8], we state without proof the following proposition. Proposition 3.1.
Assume that ( A , S ) ∈ H s ( R d ) × H s +2 ( R d ) for a large enough s . Then, the following error bound holds true (cid:107) (cid:101) ϕ h ( S , Re A , Im A ) − ( S ε ( h, · ) , Re A ε ( h, · ) , Im A ε ( h, · )) (cid:107) L ≤ Ch where C does not depend on ε ∈ [0 , ε max ] and ( S ε , A ε ) is the solution of system ofequations (1.11) . Let us remark that since S ε takes complex values, the L norm of A ε is not exactlyconserved by the flows (cid:101) ϕ h and (cid:101) ϕ h . We stress that the function S ε is not projectedon the set of real-valued functions after each flow, as it was done in Section 2.1.2,since it would reduce the order of convergence. Numerical experiments
In this part, we illustrate the behavior of the schemes (2.13) and (3.10) introducedin Sections 2 and 3 and compare their properties to those of the Strang splittingmethod [3] for (1.1). As mentioned in the introduction, quadratic observables havesome peculiarities for this problem. For this reason, the convergence properties ofthe different schemes will be illustrated separately, on the one hand for the functions S ε , A ε (resp. Ψ ε for the Strang splitting scheme) and, on the other hand, for thedensity ρ ε = | A ε | (resp. ρ ε = | Ψ ε | ). We restrict ourselves to the one-dimensionalperiodic setting in which the equations studied remain unchanged.We consider the following initial data: A ( x ) = sin( x ) , S ( x ) = sin( x ) / , Ψ ε (0 , · ) = A ( · ) e iS ( · ) /ε , (4.1)where x ∈ T = R / π Z , for which caustics appear numerically at time T c = 0 . . Inour simulations, the semiclassical parameter ε varies from to − .The numerical solutions ( S ε , A ε ) , resp. Ψ ε , are compared to corresponding refer-ence solutions ( S εref , A εref ) , resp. Ψ εref , which, in the absence of analytical solutions,are respectively obtained thanks to our fourth order splitting method (3.10) andthanks to a splitting scheme of order for (1.1) (see [33, 3]), with very small timeand space steps. More precisely, to compute ( S εref , A εref ) , we have taken N x = 2 and h = 2 − T f , and to compute Ψ εref , in order to fit with the constraints on thetime step and on the space step h (cid:28) ε and ∆ x (cid:28) ε, the space interval [0 , π ] is discretized with N x = 2 points and the time step is h = 2 − T f .The various errors that are represented on the figures below are defined as follows: err ρ ε ( T ) = (cid:107) ρ εref ( T ) − ρ ε ( T ) (cid:107) L (cid:107) ρ εref ( T ) (cid:107) L , err Ψ ε ( T ) = (cid:107) Ψ εref ( T ) − Ψ ε ( T ) (cid:107) L (cid:107) ψ εref ( T ) (cid:107) L , and err ( S ε ,A ε ) ( T ) = (cid:32) (cid:107) S εref ( T ) − S ε ( T ) (cid:107) L + (cid:107) A εref ( T ) − A ε ( T ) (cid:107) L (cid:107) S εref ( T ) (cid:107) L + (cid:107) A εref ( T ) (cid:107) L (cid:33) / , where (cid:107) u (cid:107) L = ∆ x N x − (cid:88) k =0 | u k | , (cid:107) u (cid:107) L = (cid:118)(cid:117)(cid:117)(cid:116) ∆ x N x − (cid:88) k =0 | u k | , and ρ εref ( T ) = | Ψ εref ( T ) | . As far as the Strang splitting scheme is concerned, ρ ε ( T ) = | Ψ ε ( T ) | whereas ρ ε ( T ) = | A ε ( T ) | for our methods.We first study qualitatively the dynamics, in order to guess what is the time ofappearance of the caustics. Figures 1a and 1b represent the density | A ε | and the phase S ε at times T f = 0 . , . , . , . for ε = 2 − . The caustics appear around t = 0 . . At time t = 0 . , oscillations at other scales than those of the phase can beobserved in | A ε | whereas S ε ceases to be smooth. These figures are obtained byusing our scheme (2.13) with N x = 2 and N t = T f /h = 2 .Let us now illustrate the behavior of the Strang splitting scheme for (1.1) at time T f = 0 . i.e. before the caustics. On Figures 2 and 3, errors on ρ ε and Ψ ε withrespect to the time step h , for fixed N x = 2 , are represented and on Figures 4 to 5,errors with respect to ∆ x , for fixed N t = h/T f = 2 , are represented. Regardingthe observable ρ ε = | Ψ ε | , 2a, 2b, 4a and 4b corroborate the fact that the errorbehaves as O (cid:0) h + C ε,N ∆ x N (cid:1) where N > and C ε,N → + ∞ as ε → [2, 3, 11]. This is in agreement withthe results obtained by Carles [11] in the weakly nonlinear case before the caustics;however, our simulations suggest that this behavior persists in the supercritical case.If we observe the wave function, the situation is completely different: the Strangsplitting scheme is not UA any more when h → . Figures 3a, 3b, 5a and 5b indeedsuggest that the error of Ψ ε behaves like O (cid:18) h ε + C ε,N ∆ x N (cid:19) where N > and C ε,N → + ∞ as ε → .Let us now focus on the experiments performed with our second and fourth-ordermethods, in the same situation. We start with the second-order scheme (2.13).Figures 6 and 7 represent the errors on ρ ε and ( S ε , A ε ) w.r.t. the time step h for a fixed N x = 2 . Figures 8 and 9 represent the errors w.r.t. ∆ x for fixed N t = h/T f = 2 . All these figures illustrate the fact that our scheme is UA withrespect to ε , for the quadratic observables as well as for the whole unknown ( S ε , A ε ) itself. Figures 6 and 7 show that (2.13) is uniformly of order 2 in time, whereasFigures 8 and 9 show that the convergence is uniformly spectral in space.Figures 10 to 13 illustrate the behavior of our fourth-order scheme (3.10) at T f = 0 . : here again, it appears that, before the caustics, our method is UA withan order 4 in time and with spectral in space accuracy.Finally, let us explore the behavior of the splitting methods after caustics, byobserving the error on the density ρ ε . Figures 14, 15, 16 and 17 present the samesimulations as Figures 2, 4, 6 and 8, except that the final time is now T f = 0 . , i.e. we illustrate the behaviors of Strang splitting method and of scheme (2.13) after thecaustics. In that case, it appears that none of these methods is UA, neither in h ,nor in ∆ x , with respect to ε . Concerning the Strang splitting scheme, this behaviorwas already reported in [13, 23]. Notice that, although it is not UA any longer, ourscheme (2.13) still has second-order accuracy in time and spectral accuracy in space(with ε -dependent constants). Recall that the same scheme written on (1.8) wouldnot be usable in the same situation, since S ε ceases to be regular for ε > , afterthe formation of caustics.All the experiments have been performed with the methods of both Sections 2.1.1and 2.1.2 for the eikonal equation (5.1) but we just present here the graphs obtainedwith the semi-lagrangian methods of Section 2.1.1. Final remarks
Let us emphasize that, as a by-product, we have also derived a new numericalscheme based on splitting techniques to approximate the solution of the Hamilton-Jacobi (eikonal) equation ∂ t S + |∇ S | (5.1)based on the Cole-Hopf transform. Finally, it is interesting to mention that, al-though we have chosen to focus here on the supercritical regime of (1.1) (see [10]),our approach is also relevant in other semiclassical regimes, whether the linearSchrödinger equation with a given potential, or the weakly nonlinear geometricoptics, where | Ψ ε | Ψ ε is replaced with ε | Ψ ε | Ψ ε . References [1]
G. D. Akrivis , Finite difference discretization of the cubic Schrödinger equation , IMA J.Numer. Anal., 13 (1993), pp. 115–124.[2]
W. Bao, S. Jin, and P. A. Markowich , On time-splitting spectral approximations for theSchrödinger equation in the semiclassical regime , J. Comput. Phys., 175 (2002), pp. 487–524.[3]
W. Bao, S. Jin, and P. A. Markowich , Numerical study of time-splitting spectral dis-cretizations of nonlinear Schrödinger equations in the semiclassical regimes , SIAM J. Sci.Comput., 25 (2003), pp. 27–64.[4]
C. Besse , Relaxation scheme for time dependent nonlinear Schrödinger equations , in Mathe-matical and numerical aspects of wave propagation (Santiago de Compostela, 2000), SIAM,Philadelphia, PA, 2000, pp. 605–609.[5]
C. Besse, B. Bidégaray, and S. Descombes , Order estimates in time of splitting meth-ods for the nonlinear Schrödinger equation , SIAM J. Numer. Anal., 40 (2002), pp. 26–40(electronic).[6]
C. Besse, R. Carles, and F. Méhats , An asymptotic preserving scheme based on a newformulation for nls in the semiclassical limit , Multiscale Modeling and Simulation, 11 (2013),pp. 1228–1260.[7]
S. Blanes and F. Casas , On the necessity of negative coefficients for operator splittingschemes of order higher than two , Appl. Numer. Math., 54 (2005), pp. 23–37.[8]
S. Blanes, F. Casas, P. Chartier, and A. Murua , Splitting meth-ods with complex coefficients for some classes of evolution equations
Optimized high-order splitting methods for some classes of parabolic equations , Math.Comp., 82 (2013), pp. 1559–1576.[10]
R. Carles , Semi-classical analysis for nonlinear Schrödinger equations , World Scientific,2008.[11] ,
On Fourier time-splitting methods for nonlinear Schrödinger equations in the semi-classical limit , SIAM J. Numer. Anal., 51 (2013), pp. 3232–3258.[12]
R. Carles, R. Danchin, and J.-C. Saut , Madelung, Gross-Pitaevskii and Korteweg , Non-linearity, 25 (2012), pp. 2843–2873.[13]
R. Carles and L. Gosse , Numerical aspects of nonlinear Schrödinger equations in thepresence of caustics , Math. Models Methods Appl. Sci., 17 (2007), pp. 1531–1553.[14]
F. Castella, P. Chartier, S. Descombes, and G. Vilmart , Splitting methods withcomplex times for parabolic equations , BIT, 49 (2009), pp. 487–508.[15]
T. Cazenave , Semilinear Schrödinger equations , vol. 10, AMS Bookstore, 2003.[16]
P. Degond, S. Gallego, and F. Méhats , An asymptotic preserving scheme for theSchrödinger equation in the semiclassical limit , C. R. Math. Acad. Sci. Paris, 345 (2007),pp. 531–536.[17]
M. Delfour, M. Fortin, and G. Payre , Finite-difference solutions of a nonlinearSchrödinger equation , J. Comput. Phys., 44 (1981), pp. 277–288. x | A † | Tf = 0.1Tf = 0.3Tf = 0.5Tf = 0.6 (a) Evolution of the density | A ε | for ε = 2 − . x S † Tf = 0.1Tf = 0.3Tf = 0.5Tf = 0.6 (b) Evolution of the phase S ε for ε = 2 − . Figure 1.
Evolution of the density and of the phase −4 −3 −2 −9 −8 −7 −6 −5 −4 h E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − Slope 2 (a) err ρ ε ( T f = 0 . w.r.t h , N x = 2 −2 −1 −9 −8 −7 −6 −5 −4 † E rr o r N t =1024 N t =512 N t =256 N t =128 N t =64 N t =32 N t =16 N t =8 N t =4 N t =2 (b) err ρ ε ( T f = 0 . w.r.t ε , N x = 2 Figure 2.
Error on the density ρ ε for the Strang splitting schemefor (1.1) before the caustics: dependence on ε and on the time step h . −4 −3 −2 −9 −8 −7 −6 −5 −4 −3 −2 h E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − Slope 2 (a) err Ψ ε ( T f = 0 . w.r.t h , N x = 2 −2 −1 −9 −8 −7 −6 −5 −4 −3 −2 † E rr o r N t =1024 N t =512 N t =256 N t =128 N t =64 N t =32 N t =16 N t =8 N t =4 N t =2 (b) err Ψ ε ( T f = 0 . w.r.t ε , N x = 2 Figure 3.
Error on the wave function Ψ ε for the Strang splittingscheme for (1.1) before the caustics: dependence on ε and on thetime step h . −1 −12 −10 −8 −6 −4 −2 ∆ x E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − (a) err ρ ε ( T f = 0 . w.r.t ∆ x , N t = 2 −2 −1 −12 −10 −8 −6 −4 −2 † E rr o r N x =128 N x =64 N x =32 N x =16 N x =8 (b) err ρ ε ( T f = 0 . w.r.t ε , N t = 2 Figure 4.
Error on the density ρ ε for the Strang splitting schemefor (1.1) before the caustics: dependence on ε and on ∆ x . −1 −12 −10 −8 −6 −4 −2 ∆ x E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − (a) err Ψ ε ( T f = 0 . w.r.t ∆ x , N t = 2 −2 −1 −12 −10 −8 −6 −4 −2 † E rr o r N x =128 N x =64 N x =32 N x =16 N x =8 (b) err Ψ ε ( T f = 0 . w.r.t ε , N t = 2 Figure 5.
Error on the wave function Ψ ε for the Strang splittingscheme for (1.1) before the caustics: dependence on ε and on ∆ x . −4 −3 −2 −9 −8 −7 −6 −5 −4 h E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − Slope 2 (a) err ρ ε ( T f = 0 . w.r.t h , N x = 2 −3 −2 −1 −9 −8 −7 −6 −5 −4 −3 † E rr o r N t =1024 N t =512 N t =256 N t =128 N t =64 N t =32 N t =16 N t =8 N t =4 N t =2 (b) err ρ ε ( T f = 0 . w.r.t ε , N x = 2 Figure 6.
Error on the density ρ ε for the splitting scheme (2.13)of order before the caustics: dependence on ε and on h . −4 −3 −2 −9 −8 −7 −6 −5 −4 h E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − Slope 2 (a) err ( S ε ,A ε ) ( T f = 0 . w.r.t h , N x = 2 −3 −2 −1 −9 −8 −7 −6 −5 −4 −3 † E rr o r N t =1024 N t =512 N t =256 N t =128 N t =64 N t =32 N t =16 N t =8 N t =4 N t =2 (b) err ( S ε ,A ε ) ( T f = 0 . w.r.t ε , N x = 2 Figure 7.
Error on ( S ε , A ε ) for the splitting scheme (2.13) of order before the caustics: dependence on ε and on h . −1 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 ∆ x E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − (a) err ρ ε ( T f = 0 . w.r.t ∆ x , N t = 2 −3 −2 −1 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 † E rr o r N x =128 N x =64 N x =32 N x =16 N x =8 (b) err ρ ε ( T f = 0 . w.r.t ε , N t = 2 Figure 8.
Error on the density ρ ε for the splitting scheme (2.13)of order before the caustics: dependence on ε and on ∆ x . −1 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 ∆ x E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − (a) err ( S ε ,A ε ) ( T f = 0 . w.r.t ∆ x , N t = 2 −3 −2 −1 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 † E rr o r N x =128 N x =64 N x =32 N x =16 N x =8 (b) err ( S ε ,A ε ) ( T f = 0 . w.r.t ε , N t = 2 Figure 9.
Error on ( S ε , A ε ) for the splitting scheme (2.13) of order before the caustics: dependence on ε and on ∆ x . −4 −3 −2 −13 −12 −11 −10 −9 −8 −7 h E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − Slope 4 (a) err ρ ε ( T f = 0 . w.r.t h , N x = 2 −3 −2 −1 −13 −12 −11 −10 −9 −8 −7 † E rr o r N t =1024 N t =512 N t =256 N t =128 N t =64 N t =32 N t =16 N t =8 N t =4 N t =2 (b) err ρ ε ( T f = 0 . w.r.t ε , N x = 2 Figure 10.
Error on the density ρ ε for the splitting scheme (3.10)of order before the caustics: dependence on ε and on h . −4 −3 −2 −14 −13 −12 −11 −10 −9 −8 −7 h E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − Slope 4 (a) err ( S ε ,A ε ) ( T f = 0 . w.r.t h , N x = 2 −3 −2 −1 −14 −13 −12 −11 −10 −9 −8 −7 † E rr o r N t =1024 N t =512 N t =256 N t =128 N t =64 N t =32 N t =16 N t =8 N t =4 N t =2 (b) err ( S ε ,A ε ) ( T f = 0 . w.r.t ε , N x = 2 Figure 11.
Error on ( S ε , A ε ) for the splitting scheme (3.10) of order before the caustics: dependence on ε and on h . −1 −12 −10 −8 −6 −4 ∆ x E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − (a) err ρ ε ( T f = 0 . w.r.t ∆ x , N t = 2 −3 −2 −1 −12 −10 −8 −6 −4 † E rr o r N x =128 N x =64 N x =32 N x =16 N x =8 (b) err ρ ε ( T f = 0 . w.r.t ε , N t = 2 Figure 12.
Error on ρ ε for the splitting scheme (3.10) of order before the caustics: dependence on ε and on ∆ x . −1 −14 −12 −10 −8 −6 −4 ∆ x E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − (a) err ( S ε ,A ε ) ( T f = 0 . w.r.t ∆ x , N t = 2 −3 −2 −1 −14 −12 −10 −8 −6 −4 † E rr o r N x =128 N x =64 N x =32 N x =16 N x =8 (b) err ( S ε ,A ε ) ( T f = 0 . w.r.t ε , N t = 2 Figure 13.
Error on ( S ε , A ε ) for the splitting scheme (3.10) of order before the caustics: dependence on ε and on ∆ x . −3 −2 −1 −7 −6 −5 −4 −3 −2 −1 h E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − Slope 2 (a) err ρ ε ( T f = 0 . w.r.t h , N x = 2 −2 −1 −7 −6 −5 −4 −3 −2 −1 † E rr o r N t =1024 N t =512 N t =256 N t =128 N t =64 N t =32 N t =16 N t =8 N t =4 N t =2 (b) err ρ ε ( T f = 0 . w.r.t ε , N x = 2 Figure 14.
Error on ρ ε for the Strang splitting scheme for (1.1)after the caustics, dependence on ε and on h . −1 −10 −8 −6 −4 −2 ∆ x E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − (a) err ρ ε ( T f = 0 . w.r.t ∆ x , N t = 2 −2 −1 −10 −8 −6 −4 −2 † E rr o r N x =128 N x =64 N x =32 N x =16 N x =8 (b) err ρ ε ( T f = 0 . w.r.t ∆ x , N t = 2 Figure 15.
Error on ρ ε for the Strang splitting scheme for (1.1)after the caustics, dependence on ε and on ∆ x . [18] S. Descombes and M. Thalhammer , An exact local error representation of exponentialoperator splitting methods for evolutionary problems and applications to linear Schrödingerequations in the semi-classical regime , BIT, 50 (2010), pp. 729–749.[19] ,
The Lie-Trotter splitting for nonlinear evolutionary problems with critical parameters:a compact local error representation and application to nonlinear Schrödinger equations in thesemiclassical regime , IMA J. Numer. Anal., 33 (2013), pp. 722–745.[20]
L. C. Evans , Partial differential equations , Providence, Rhode Land: American MathematicalSociety, 1998.[21]
E. Grenier , Semiclassical limit of the nonlinear Schrödinger equation in small time , Proc.Amer. Math. Soc., 126 (1998), pp. 523–530.[22]
E. Hansen and A. Ostermann , High order splitting methods for analytic semigroups exist ,BIT, 49 (2009), pp. 527–542.[23]
S. Jin, P. Markowich, and C. Sparber , Mathematical and computational methods forsemiclassical Schrödinger equations , Acta Numer., 20 (2011), pp. 121–209.[24]
O. Karakashian, G. D. Akrivis, and V. A. Dougalis , On optimal order error estimatesfor the nonlinear Schrödinger equation , SIAM J. Numer. Anal., 30 (1993), pp. 377–400.[25]
C. Lubich , On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödingerequations , Math. Comp., 77 (2008), pp. 2141–2153.[26]
E. Madelung , Quanten theorie in hydrodynamischer form , Zeit. F. Physik, 40 (1927),pp. 322–326.[27]
R. I. McLachlan and G. R. W. Quispel , Splitting methods , Acta Numer., 11 (2002),pp. 341–434.[28]
D. Pathria and J. L. Morris , Pseudo-spectral solution of nonlinear Schrödinger equations ,J. Comput. Phys., 87 (1990), pp. 108–125.[29]
J. M. Sanz-Serna and J. G. Verwer , Conservative and nonconservative schemes for thesolution of the nonlinear Schrödinger equation , IMA J. Numer. Anal., 6 (1986), pp. 25–42.[30]
C. Sulem and P.-L. Sulem , The nonlinear Schrödinger equation , vol. 139 of Applied Math-ematical Sciences, Springer-Verlag, New York, 1999. Self-focusing and wave collapse.[31]
J. A. C. Weideman and B. M. Herbst , Split-step methods for the solution of the nonlinearSchrödinger equation , SIAM J. Numer. Anal., 23 (1986), pp. 485–507.[32]
L. Wu , Dufort-Frankel-type methods for linear and nonlinear Schrödinger equations , SIAMJ. Numer. Anal., 33 (1996), pp. 1526–1533.[33]
H. Yoshida , Construction of higher order symplectic integrators , Phys. Lett. A, 150 (1990),pp. 262–268.
E-mail address : [email protected]
INRIA Rennes, IRMAR and ENS Rennes, IPSO Project Team, Campus de Beaulieu,F-35042 Rennes
E-mail address : [email protected] IRMAR, Université de Rennes 1 and INRIA, IPSO Project
E-mail address : [email protected] IRMAR, Université de Rennes 1 and INRIA, IPSO Project −3 −2 −1 −7 −6 −5 −4 −3 −2 −1 h E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − Slope 2 (a) err ρ ε ( T f = 0 . w.r.t h , N x = 2 −2 −1 −7 −6 −5 −4 −3 −2 −1 † E rr o r N t =1024 N t =512 N t =256 N t =128 N t =64 N t =32 N t =16 N t =8 N t =4 N t =2 (b) err ρ ε ( T f = 0 . w.r.t ε , N x = 2 Figure 16.
Error on ρ ε for the splitting scheme (2.13) of order after the caustics, dependence on ε and on h . −1 −10 −8 −6 −4 −2 ∆ x E rr o r † = 2 − † = 2 − † = 2 − † = 2 − † = 2 − (a) err ρ ε ( T f = 0 . w.r.t ∆ x , N t = 2 −2 −1 −10 −8 −6 −4 −2 † E rr o r N x =128 N x =64 N x =32 N x =16 N x =8 (b) err ρ ε ( T f = 0 . w.r.t ε , N t = 2 Figure 17.
Error on ρ ε for the splitting scheme (2.13) of order after the caustics, dependence on ε and on ∆ xx