An exponential integrator for a highly oscillatory Vlasov equation
AAn exponential integrator for a highly oscillatory Vlasov equation
Emmanuel Fr´enod ∗ Sever A. Hirstoaga † Eric Sonnendr¨ucker ‡ Abstract
In the framework of a Particle-In-Cell scheme for some 1D Vlasov-Poisson system de-pending on a small parameter, we propose a time-stepping method which is numericallyuniformly accurate when the parameter goes to zero. Based on an exponential time differ-encing approach, the scheme is able to use large time steps with respect to the typical sizeof the fast oscillations of the solution.
The problem of studying a charged particle beam focused by external electric or magnetic fieldsis important in many applications. The motion of such particle beams is governed by theinteractions between the electric field generated by the particles themselves and by the externalfocusing electromagnetic field. The modelling framework consists in coupling a kinetic equationwith Maxwell equations. Disregarding the collisions between particles, the kinetic modelling isperformed by means of the Vlasov equation. In this paper we will consider only non-relativisticlong and thin beams. Therefore, instead of studying the phenomenon by means of the full Vlasov-Maxwell system, we can use its paraxial approximation (see [4] for mathematical modelling andnumerical simulation of focused particle beams dynamics). In this framework, the effects ofthe self-consistent magnetic and electric fields can be both taken into account by solving asingle Poisson equation. Finally, we are led to solve a two dimensional phase space Vlasov-Poisson system with a parameter ε . This small parameter acts on the time variable (in fact thelongitudinal variable of a thin beam, see [6] for details and physical meaning of the parameter),producing the highly oscillatory behaviour in phase space. In this framework, the aim of thispaper is to propose a numerical scheme able to study efficiently the evolution of a beam over alarge number of fast periods.Although more precise alternatives exist, we have chosen in this work to perform the nu-merical solution of Vlasov equation by particle methods, which consist in approximating thedistribution function by a finite number of macroparticles. The trajectories of these particlesare computed from the characteristic curves of the Vlasov equation, whereas the self-consistentelectric field is computed on a mesh in the physical space. This method allows to obtain satisfy-ing results with a small number of particles (see [1]). The contribution of this paper is to proposea new numerical method for solving the characteristic curves, or equivalently, for computing the ∗ Universit´e Europ´eenne de Bretagne, LMBA (UMR CNRS 6205), Universit´e de Bretagne-Sud & Inria Nancy-Grand Est, CALVI Project & IRMA (UMR CNRS 7501), Universit´e de Strasbourg, France † Inria Nancy-Grand Est, CALVI Project & IRMA (UMR CNRS 7501), Universit´e de Strasbourg, France ‡ Max Planck Institute for Plasma Physics, EURATOM Association, Boltzmannstr. 2, 85748 Garching, Ger-many a r X i v : . [ m a t h . NA ] J un acroparticles’ trajectories. Namely, we are faced with solving the following stiff differentialequation u (cid:48) ( t ) = u ( t ) /ε + F ( t, u ( t )) , u (0) = u , (1.1)for several small values of the parameter ε and where F represents a nonlinear term whichplays the role of the self-consistent electric field. The difficulty arising in the numerical solutionfor this equation relies on the ability of the scheme to be uniformly stable and accurate when ε goes to zero. Following the survey article [7], we are encouraged to use efficient numericalmethods like exponential integrators in order to describe the dynamics of equation (1.1). Thebasic idea behind these methods is to make use of the solution’s exact representation given bythe variation-of-constants formula u ( t ) = e t/ε u + (cid:90) t e ( t − τ ) /ε F ( τ, u ( τ )) dτ. (1.2)Applying this formula from one time step to another has the merit to solve exactly the linear(stiff) part. Classical numerical schemes fail to capture the stiff behaviour regardless of the sizeof the time step with respect to the small parameter. One may consult [7] for construction,mathematical analysis and implementation of exponential integrators in the two classical typesof stiff problems encapsulated in equation (1.1). As a specific implementation, one may cite[9], where, in the context of laser-plasma interactions, an exponential integrator is used for aParticle-In-Cell (PIC) method in order to model the high-frequency plasma response. Once thestiff part is exactly solved, one may use an exponential time differencing (ETD) method (see[3]) for the specific numerical treatement of the nonlinear term. The ETD schemes turn out tooutperform many other schemes when treating problems like (1.1); see [3, 8] for comparisons ofETD against the Implicit-Explicit method, the Integrating Factor method, etc.In the present paper we construct and implement an exponential integrator in order to solvethe characteristics of the highly oscillatory Vlasov-Poisson problem (2.2). The aim is to usea scheme with large time steps compared to the fast period that arises from the linear termwithout loosing the accuracy when the small parameter vanishes. The novelty of this method isin the numerical approximation of the integral term in (1.2). More precisely, when the time stepis much larger than the rapid period, the idea of the algorithm is the following: we first finelysolve the ODEs over one fast period by means of a high-order solver (we have used explicit 4thorder Runge-Kutta). Then, thanks to formula (1.2), we may compute an approximation of thesolution over a large integral number of periods. We have also found that using a more accurateperiod, instead of the period of the solution to the system (1.1) without the nonlinear term,leads to more accurate simulations. In addition, we have checked if the scheme gives accuratesolutions starting with an initial condition which lies on the slow manifold or not. We cite [2] fora “definition” of the slow manifold: “The slow manifold is that particular solution which variesonly on the slow time scale; the general solution to the ODE contains fast oscillations also.”The remainder of the paper is organized as follows. Following [6], we briefly recall in Section2 the paraxial approximation together with the axisymmetric beam assumption. In Section 3we describe the PIC method for the Vlasov-Poisson system in which we are interested. Then,Section 4 is devoted to the construction of the new numerical scheme as an exponential integratorfor solving the time-stepping part of the PIC algorithm. Eventually, in Section 5, we implementand test our method on several test cases related to the Vlasov-Poisson system.2 The paraxial model: an axisymmetric Vlasov equation
The paraxial approximation relies on a scaled Vlasov-Poisson system in a phase space of di-mension four, 2 D for space variable x and 2 D for velocity variable v . This simplified model ofthe full Vlasov-Maxwell system is particularly adapted to the study of long and thin beams ofcharged particles and it describes their evolution in the plane transverse to their direction ofpropagation.Subject of many research investigations, the paraxial model was derived in a number ofpapers, see e.g. [4, 6]. In this work we are interested in solving numerically a paraxial modelwith some additional hypotheses, see (2.2) below. The solution of system (2.2) is represented bya beam of particles in phase space. The beam evolves by rotating around the origin in the phasespace, and in long times a bunch forms around the center of the beam from which filaments ofparticles are going out. These filaments are difficult to capture with classical numerical methods.We now introduce the paraxial model that we aim to solve. In the additional axisymmetricbeam assumption (i.e. invariant beam under rotation in R for x ), we are led to change the x coordinate in the polar frame. We thus write the model in polar coordinates ( r, θ ), where r = | x | and θ ∈ [0 , π ) is such that x = r cos θ and x = r sin θ . Then we use new velocity variables v r = v · x /r and v θ = v · x ⊥ /r , where x ⊥ = ( − x , x ). Assuming in addition, as in [6], that f ε is concentrated in angular momentum, i.e. rv θ = 0, the paraxial model becomes ∂f ε ∂t + 1 ε v r ∂f ε ∂r + (cid:0) E εr + Ξ εr (cid:1) ∂f ε ∂v r = 0 , r ∂ ( r E εr ) ∂r = ρ ε , ρ ε ( t, r ) = (cid:90) R f ε ( t, r, v r ) dv r ,f ε ( t = 0 , r, v r ) = f ( r, v r ) , (2.1)where the external force Ξ εr writesΞ εr ( t, r ) = (cid:16) − ε H + H (cid:0) tε (cid:1)(cid:17) r, with H > H ( t ) some 2 π -periodic function with zero mean value. In this paper we assumethat there is no time oscillation in the external field but only the strong uniform focusing. Wethus take H = 1 and H = 0 and the Vlasov-Poisson system in which we are interested writes ∂f ε ∂t + 1 ε v r ∂f ε ∂r + (cid:16) E εr − rε (cid:17) ∂f ε ∂v r = 0 , r ∂ ( r E εr ) ∂r = (cid:90) R f ε ( t, r, v r ) dv r ,f ε ( t = 0 , r, v r ) = f ( r, v r ) . (2.2)In order to test the numerical method that we propose, we first consider two numerically simplertest cases where the electric field is not issued from Poisson equation, but it has analytical forms:in the first case E εr ( t, r ) = − r and in the second one E εr ( t, r ) = − r . Unlike the second case, forthe first one, we can analytically compute the solution to (2.2)(a).3 A Particle-In-Cell method for Vlasov-Poisson system
We solve the Vlasov-Poisson system (2.2) by using a Particle-In-Cell (PIC) method [1]. We thusintroduce the following Dirac mass sum approximation of f ε f N p ε ( r, v, t ) = N p (cid:88) k =1 w k δ ( r − R k ( t )) δ ( v − V k ( t )) , where N p is the number of macroparticles and (cid:0) R k ( t ) , V k ( t ) (cid:1) is the position in phase space ofmacroparticle k moving along a characteristic curve of equation (2.2)(a). Therefore, the problemis to find the positions and velocities ( R n +1 k , V n +1 k ) at time t n +1 from their values at time t n , bysolving R (cid:48) ( t ) = V ( t ) /ε,V (cid:48) ( t ) = − R ( t ) /ε + E ε ( t, R ( t )) ,R ( t n ) = R nk , V ( t n ) = V nk . (3.1)In this case, the standard PIC algorithm writes as follows: (1) deposit particles on a spatial grid,leading to the grid density; (2) solve Poisson equation on the grid, leading to the grid electricfield; (3) interpolate the grid electric field in each particle; (4) push particles with the previouslyobtained electric field.The first three steps are classically treated. The first one deals with the computing of thegrid density by convoluting ρ N p ε defined by ρ N p ε ( r, t ) = N p (cid:88) k =1 w k δ ( r − R k ( t )) , with a first order spline (this corresponds to the Cloud-in-Cell method in [1]). Then, we solvePoisson equation (2.2)(b) on a uniform one-dimensional grid by using finite differences. Inour case this amounts to only discretize some space integral. We have done this by usingthe trapezoidal rule. As for the third step, we used the same convolution function as for thedeposition step in order to get the particle electric field (this corresponds to a linear interpolationof the grid electric field on each cell).Eventually, the major issue is the fourth step of the PIC algorithm, consisting in the numer-ical integration of system (3.1). Here is the main focus of this paper, taking into account thatwe want to propose a stable and accurate scheme using large time steps with respect to the fastoscillation. To this end, we introduce in the next section a method based on exponential timedifferencing. We now describe the exponential numerical integrator that we have implemented to solve thefourth step of the PIC algorithm. We first write down the exponential time differencing (ETD)method in the case of the stiff ODE system we are interested in. Then, in Section 4.2, we developan algorithm based on the exponential time differencing in the framework of the PIC method.4 .1 Exponential time differencing for a highly oscillatory ODE
The so-called exponential time differencing scheme arose originally in the field of computationalelectrodynamics but has been reinvented many times over the years (see [3] and the referencestherein). We take details of the ideas behind the various ETD schemes from the comprehensivepaper by Cox and Matthews [3].Recall the stiff system of ODEs that we have to solve: R (cid:48) ( t ) = 1 ε V ( t ) ,V (cid:48) ( t ) = − ε R ( t ) + E ε ( t, R ( t )) , (4.1)with some initial condition ( R , V ). In this section we assume that the electric field E ε is given.As exposed in [3], the stiffness comes from the two scales on which the solution evolves: therapid oscillations due to the linear term and a slower evolution due to the nonlinear (electric)term. Thus, while any explicit scheme is limited to a small time step, of order ε , a fully implicitone requires nonlinear problems to be solved and is therefore slow. A suitable time-steppingscheme for (4.1) should be able to avoid the small time steps when treating the stiffness. Theessence of the ETD methods is to solve the stiff linear part exactly and to derive appropriateapproximations when integrating numerically the slower nonlinear term [3, 8].To derive the exponential time differencing (ETD) method for this system we first apply r ( − t/ε ) to (4.1) and then integrate the obtained equation from t n to t n +1 = t n + dt to deducethat (cid:18) R n +1 V n +1 (cid:19) = r (cid:16) dtε (cid:17) (cid:18) R n V n (cid:19) + r (cid:16) dtε (cid:17) (cid:90) t n +1 t n r (cid:16) t n − τε (cid:17) (cid:18) E ε ( τ, R ( τ )) (cid:19) dτ, (4.2)where r ( τ ) = (cid:18) cos ( τ ) sin ( τ ) − sin ( τ ) cos ( τ ) (cid:19) . (4.3)This formula is exact. In addition, it is also useful to write (4.2) by replacing ( t n , t n +1 ) by anycouple ( s, t ) such that s < t . More precisely, we have (cid:18) R ( t ) V ( t ) (cid:19) = r (cid:16) t − sε (cid:17) (cid:18) R ( s ) V ( s ) (cid:19) + r (cid:16) t − sε (cid:17) (cid:90) ts r (cid:16) s − τε (cid:17) (cid:18) E ε ( τ, R ( τ )) (cid:19) dτ. (4.4)Now the main question is how to derive approximations to the integral term in (4.2). All theETD schemes are results of this process. In this spirit, it was shown in [3] that ETD methodscan be extended to any order by using multistep or Runge-Kutta methods and explicit formulaefor such arbitrary order ETD methods were derived. In particular, explicit coefficients for ETDRunge-Kutta methods of order up to four have been computed. The authors also illustrated onseveral ODEs and PDEs that ETD is superior over the Integrating Factor and Implicit-Explicitmethods, two other classical schemes able to avoid the small time step. Nevertheless, in theform written in [3], a high-order ETD scheme (e.g. ETDRK4) suffers from numerical instabilityas explained in [8]. The problem has been solved in [8] by using a contour integral method forevaluating the coefficients. Then, this modified ETDRK4 scheme has been tested in [8] againstfive other 4th order schemes on several PDEs, and ETDRK4 has been found the best in termsof errors. These results encouraged us to use an ETD scheme in order to solve stiff ODEs.5n this paper we do not use a high-order ETD method as described in the references abovebut merely the formula (4.2) for justifying the derivation of Algorithm 4.1. More precisely, weadopt a different approach for approximating the integral term in (4.2), which is justified bythe remark in the next paragraph and by our aim to use very large time steps compared to thefast oscillations, say dt = 100 ε when ε = 10 − . This is the core of the method described in thefollowing section.Figure 1: A nonlinear case: the global Euclidean error at time 3.14 for two initial conditions:on (left) and off (right) the slow manifoldHowever, we have done tests with ETDRK2 (see formula (22) in [3]) against symplecticVerlet and RK4 schemes, in the case E ε ( t, R ) = − R with ε = 10 − (see the test case in Section5.2). Thus, we have noticed that this ETD scheme has smaller errors for some range for thetime step and that, unlike its competitors, it is stable when dt ≥ ε (see Fig. 1). But thesesimulations also show that when dt is too large with respect to ε (e.g. dt ≥ ε ) the global error issignificant for some initial condition (particles off the slow manifold, see Section 5). The reasonis that for dt large enough with respect to ε , the nonlinear force term cannot be accurately takeninto account by any high-order Runge-Kutta solver. In this section we describe our time-stepping scheme starting from the equation (4.2). The mainidea is to remark that the time for one fast grand tour is approximated by 2 πε , the period ofthe solution to (4.1) without nonlinear term. Therefore, since we want to build a scheme witha time step dt much larger than the fast oscillation, we first need to find the unique positiveinteger N and the unique real ∆ ∈ [0 , πε ) such that dt = N · (2 πε ) + ∆ . (4.5)Thus the integral term in (4.2) becomes (cid:90) t n +1 t n dτ = N − (cid:88) j =0 (cid:90) t n +2 πε ( j +1) t n +2 πε j dτ + (cid:90) t n +1 t n +1 − ∆ dτ, (4.6)6hat we approximate by N (cid:90) t n +2 πεt n dτ + (cid:90) t n +1 t n +1 − ∆ dτ. (4.7)In this approximation, we make some assumptions that we develop in Remark 4.2. Computingthe integrals I = (cid:90) t n +2 πεt n r (cid:16) t n − τε (cid:17) (cid:18) E ε ( τ, R ( τ )) (cid:19) dτ and I = (cid:90) t n +1 t n +1 − ∆ r (cid:16) t n − τε (cid:17) (cid:18) E ε ( τ, R ( τ )) (cid:19) dτ knowing ( R n , V n ) will thus lead to ( R n +1 , V n +1 ). Next we describe this method in 3 steps. Theelectric field is discretized explicitely, i.e. it can be computed at any time solving the Poissonequation (2.2)(b) with a right-hand side obtained from depositing the particles on the grid. : Using (4.4) with s = t n and t = t n + 2 πε and since r (2 π ) = Id, we have I = (cid:18) R ( t n + 2 πε ) − R ( t n ) V ( t n + 2 πε ) − V ( t n ) (cid:19) . (4.8)Now, we finely solve (4.1) with initial conditions ( R n , V n ) by a 4th order Runge-Kutta scheme, inorder to get an accurate approximation of (cid:0) R ( t n + 2 πε ) , V ( t n + 2 πε ) (cid:1) . This is done by followingthe four steps of the standard PIC algorithm as exposed in Section 3. The time step used in thefourth step for the RK4 solver is √ ε ε (we explain this choice in Section 5). : We have to calculate an approximation of (cid:0) R ( t n + N · πε ) , V ( t n + N · πε ) (cid:1) since it will be usefull in the computations to do in the 3rd step. Using (4.4) with s = t n and t = t n + N · πε we obtain, since r ( N · π ) = Id, (cid:18) R ( t n + N · πε ) V ( t n + N · πε ) (cid:19) ≈ (cid:18) R n V n (cid:19) + N · I , (4.9)where the right integral (the first one in the right-hand side of (4.6)) was approximated as donein (4.7) by N · I . : Now, we compute I by using (4.4) with s = t n + N · πε = t n +1 − ∆ and t = t n +1 I = r (cid:16) − ∆ ε (cid:17) (cid:32) (cid:101) R ( t n +1 ) (cid:101) V ( t n +1 ) (cid:33) − (cid:18) R ( t n + N · πε ) V ( t n + N · πε ) (cid:19) , (4.10)where (cid:0) (cid:101) R ( t n +1 ) , (cid:101) V ( t n +1 ) (cid:1) may be found as done above for the approximation of (cid:0) R ( t n + 2 πε ) ,V ( t n + 2 πε ) (cid:1) : we follow the steps of the standard PIC algorithm, using for the particles’ pusha 4th order Runge-Kutta solver with initial conditions ( R ( t n + N · πε ) , V ( t n + N · πε )) andwith total time ∆ (as for the 1st step above, we choose a small time step, √ ε ε ).Finally, we replace (4.9) in (4.10) which we put in (4.2): the term N · I will cancel and atthe end we have R n +1 = (cid:101) R ( t n +1 ) and V n +1 = (cid:101) V ( t n +1 ) , meaning that the vector (cid:0) (cid:101) R ( t n +1 ) , (cid:101) V ( t n +1 ) (cid:1) calculated within the 3rd step above is an approxi-mation of the solution at time t n +1 .In summary, writing dt as in (4.5), the implementation of our time-stepping algorithm followsthree steps: 7 lgorithm 4.1. Assume that ( R n , V n ) the solution of (4.1) at time t n is given. Then1. compute ( R, V ) at time t n + 2 πε by using a fine Runge-Kutta solver with initial condition ( R n , V n ) .2. compute ( R, V ) at time t n + N · πε by the following rule (cid:18) R ( t n + N · πε ) V ( t n + N · πε ) (cid:19) = (cid:18) R n V n (cid:19) + N (cid:18) R ( t n + 2 πε ) − R n V ( t n + 2 πε ) − V n (cid:19) . (4.11)
3. compute ( R, V ) at time t n +1 by using a fine Runge-Kutta solver with initial condition ( R, V ) obtained at the previous step. In the framework of the PIC method, Algorithm 4.1 is applied to each particle.We will call the modified ETD-PIC scheme , the Algorithm 4.1 where 2 πε is replaced in thefirst two steps with a more accurate fast time for particles. The reason for this choice is explainedin Sections 5.1 and 5.4. Remark 4.2.
In the approximation (4.7) , we have made the assumptions that the solution’speriod does not change significantly in time and the same for the electric field E ε . Nevertheless,although in this section E ε in supposed to be given, in the case of the Vlasov-Poisson coupling (2.2) , E ε is the self-consistent electric field. Therefore, we expect an almost periodic trajectory ofthe charged particles to involve a similar time behaviour for the electric field that they generate.Thus, in such a case, we make only the assumption that the particles’ period does not varysignificantly in time. We now validate our algorithms in the test cases presented in Section 2. We illustrate theirnumerical performance and desired properties exposed in the Introduction with global errorcurves. Thus, we have checked if the method is numerically stable and accurate when ε becomessmaller. Since we solve a stiff system, we compute the errors for both types of initial particles,on and off the slow manifold (see [2, 3]). In fact, for numerical reasons, we should rather use thedesignations “close or far from the slow manifold” than “on or off the slow manifold”. Never-theless, in order to be in line with the literature we cite, we keep in this paper the designations“on or off the slow manifold”.In order to solve the Vlasov-Poisson model (2.2), we choose as initial distribution functionthe following semi-Gaussian beam (see [6]) f ( r, v ) = n √ π v th exp (cid:16) − v v (cid:17) χ [ − . , . ( r ) , (5.1)where the thermal velocity is v th = 0 . χ [ − . , . ( r ) = 1 if r ∈ [ − . , . N p =10000 macroparticles with equal weights w k = 1 /N p .8 .1 A linear case As a first test of validation of our scheme we consider a case where the solution is analyticallyknown. It is the linear case where the self-consistent electric field in (4.1) is given by E ε ( t, R ) = − R . The solution satisfying the initial condition R ( s ) = R s , V ( s ) = V s is R ( t ) = V s √ ε sin (cid:16) √ εε ( t − s ) (cid:17) + R s cos (cid:16) √ εε ( t − s ) (cid:17) V ( t ) = V s cos (cid:16) √ εε ( t − s ) (cid:17) − R s √ ε sin (cid:16) √ εε ( t − s ) (cid:17) (5.2)for all t ≥ s . We note that the trajectory in the phase space is an ellipse since we have ∀ t ≥ s R ( t ) + (cid:16) V ( t ) √ ε (cid:17) = R s + (cid:16) V s √ ε (cid:17) . In this particular case, we can compute exactly the time for one rapid grand tour. It is t p = 2 πε √ ε (5.3)and this value (and not 2 πε ) was used in our simulation for the push of particles within Algorithm4.1. In fact, using 2 πε instead of the right rapid time makes our method unstable when the finaltime is large enough. The reason is that 2 πε − t p is big enough so that the accumulated in timeerrors issued from the rule (4.11) lead the particles to drift outward in the phase space. Wealso notice that, unlike general case, the rapid time does not depend on the initial condition.Thus, the beam of particles is rotating in the phase space without spiraling. This was checkednumerically for our scheme.In addition, in this case, the slow manifold can easily be determined. Knowing that it isthat particular solution which varies only on the slow time scale ([2]), we deduce that the fastoscillations in the solution given by (5.2) can be removed when R = 0 and V = 0. Thestationary solution ( R, V ) = (0 ,
0) is the slow manifold in this particular case. Numerically, onand off the slow manifold particles were arbitrarly chosen with R = 0 . V ∼ · − (close to (0 , R = 0 . V ∼ . dt = √ ε ε . More precisely, we have first noticed the need for the choice of a uniformly smallratio dt/ε with respect to ε in order to get an accurate solution. Indeed, we have computed areference solution with a time step dt = ε/
10, when ε = 10 − . The outcome of this simulation atfinal times t = 1 . t = 3 . dt = ε/
10 is notsufficiently small so that good accuracy for the reference solution be reached. Second, we haveexperimented different uniformly small time steps, dt = ε / , ε / , ε / , ε . Thus, computingsolutions with time steps dt ≤ ε / do not lead to considerably smaller error, for small ε , and inaddition it requires large CPU time. We have finally declared acceptable, at worst of order 10 − ,the errors obtained with dt = ε / , and therefore, this time step was chosen in the following testcases when computing the reference solution. 9igure 2: Computing a reference solution in phase space: use of an explicit Runge-Kutta solverwith dt = ε/
10 when ε = 10 − at final time 1.0 (left) and 3.5 (right) We now take the case of E ε ( t, R ) = − R . The solution to system (4.1) can be writen in termsof Jacobi elliptic functions and it cannot be put in an analytical form. We therefore compute areference solution using a fine explicit Runge-Kutta solver. Unlike the nonlinear case in Section5.3, we do not have numerical errors due to the computation of the electric field.Figure 3: First nonlinear case: the global Euclidean error at time 3.5 for an initial particle on(left) and off (right) the slow manifoldWe have used the initial conditions, on and off the slow manifold, mentioned in the previoussection, since ( R, V ) = (0 ,
0) is still a stationary solution to problem (4.1). We have appliedAlgorithm 4.1 for each particle. Denoting by ( R ( t ) , V ( t )) the result of the ETD-PIC methodand by ( R ref ( t ) , V ref ( t )) the reference solution of problem (4.1) for some initial condition, the10lobal errors (the max value in time of the Euclidean norm of ( R ( t ) − R ref ( t ) , V ( t ) − V ref ( t ))) atfinal time 3 . ε . The error curves in the left panel arederived with the above initial condition on the slow manifold.Figure 4: Vlasov-Poisson case: the global Euclidean error at time 3.5 for an initial particle on(left) and off (right) the slow manifold No analytical solution to system (4.1) is available in the case where E ε is the solution of Poissonequation (2.2)(b). We solve numerically this equation by using a trapezoidal formula for theintegral in r with 128 cells for the space interval. As above, the time step of the RK4 solver forcomputing the reference solution is dt = √ ε ε and we have used the same initial conditions, onand off the slow manifold. The global errors of the ETD-PIC method are shown in Fig. 4 fordifferent values of ε . Remark 5.1.
1. From Figs. 3 and 4 one can see the announced property of the scheme, thatis the uniform accuracy when ε goes to . We also observe that the order of accuracydecreases uniformly when ε goes to .2. Note the very large time steps with respect to ε that the ETD scheme allows to use. Forinstance, when the time step is . this means that it goes from ε when ε = 10 − to ε when ε = 10 − . First, we notice that for each nonlinear case above, the errors are larger when the initialcondition is off the slow manifold. Then, when ε = 10 − , we can see in Fig. 5 that the differencebetween errors for off and on the slow manifold particles is more significant in the first nonlineartest case.Second, still when ε = 10 − , we have compared the two nonlinear test cases. When simulationis done with an initial particle on the slow manifold (ONSM), the error is much more significant11igure 5: Global errors at time 3.5 for initial particles on and off the slow manifold in bothnonlinear cases (1st case is the one in Section 5.2, 2nd case is the Vlasov-Poisson case)for the second nonlinear case, as expected because of the additional numerical errors due toPoisson solver (see Fig. 5). Surprisingly, for an initial particle off the slow manifold (OFFSM)the error behaves conversely, although they are of the same order of magnitude.We think the reason for both comments above is in the use of 2 πε as the fast period for allthe particles, which is more or less accurate. Precisely, we have experimentally found that whilethe period of an ONSM particle in the first nonlinear case is very close to 2 πε (thus justifyingthe much smaller error at the right of Fig. 5), in the second nonlinear case this period is ratherclose to 2 πε + 3 ε . In addition, for an OFFSM particle our experiments lead to fast times closeto 2 πε − ε in the first nonlinear case and to 2 πε + ε in the second nonlinear case. Even ifthese numbers seem to be small to induce such significant errors, we recall that in the lineartest case above, 2 πε − t p is very close to 3 ε and the use of a not accurate fast time for particlesleads to an unstable simulation. That is why we test in the following section modified ETD-PICschemes , where we use a more accurate fast time than 2 πε . Finally, we present the numerical results obtained with Algorithm 4.1 where 2 πε was replacedwith a more accurate period. More precisely, we first compute numerically an approximation ofthe period for each particle, based on the fact that particles trajectories in time are sinusoid-likealmost periodic functions. The trajectories are computed with the same Runge-Kutta solverused for the reference solution’s calculation. Thus, the particles are pushed during some verysmall time until all the particles reach their trajectory’s third extremum. The criterion forfinding these extrema is the velocity’s change of sign. We then stated that each particle’s periodis the time interval between the first and the third extremum. Eventually, we compute the meanof these periods. It is by this mean that we replace 2 πε in the first two steps of Algorithm 4.1.Summarising, the implementation of the modified ETD-PIC scheme starts with the finding ofthe mean period and then uses this value in Algorithm 4.1 all along the simulation.Our experiments show that, at best, in the coupling with Poisson case, the error is smaller12igure 6: First nonlinear case, Algorithm 4.1 with the particles mean period: the global Eu-clidean error at time 3.5 for an initial particle on (left) and off (right) the slow manifoldby a factor of 10 than that computed with period 2 πε (see Figs. 4 and 7). For the nonlinear casein Section 5.2, as expected (see Section 5.4), the errors are slightly larger than those computedwith period 2 πε for an initial ONSM particle, but smaller for an initial OFFSM particle.At the end, we also illustrate the results of ETD-PIC scheme vs. modified ETD-PIC scheme by representing the solution to Vlasov equation with an electric field, first given as in Section5.2 (Fig. 8) and second as solution of the Poisson equation (Fig. 9). These beams of particleswere obtained with the larger time step used for the errors calculus, dt = 0 . t = 3 . In this paper we have introduced a new numerical scheme for solving some stiff (highly oscil-latory) differential equation. This scheme is based on exponential time differencing and canaccurately handle large time steps with respect to the fast oscillation of the solution. It is ap-plied in the framework of a Particle-In-Cell method for solving some Vlasov-Poisson equation.Since the numerical results are encouraging, several ways to explore in the future may be usefullto improve some points in the scheme development.We have seen that the use of 2 πε as fast time within the first step of Algorithm 4.1 may leadto an unstable simulation. In addition, even in a stable simulation, the use of the particles meanperiod gives smaller errors than those obtained with 2 πε . Therefore, we think it is importantto find theoretically a more accurate approximation of the fast time.It will be interesting to see if our numerical scheme preserves the two-scale asymptotic limit,meaning that ε = 0 in the numerical scheme leads to a consistent discretization of the two-scalelimit model. This remark is based on the fact that the ETD discretization we have used is veryclose to an explicit discretization of the limit model in Theorem 1.1 of [5].Last but not least, the ETD-PIC schemes proposed in this paper need to be tested on othersystems where different types of stiff differential equations are to be solved.13igure 7: Vlasov-Poisson case, Algorithm 4.1 with the particles mean period: the global Eu-clidean error at time 3.5 for an initial particle on (left) and off (right) the slow manifold References [1]
C. K. Birdsall, A. B. Langdon , Plasma physics via computer simulation, Institute ofPhysics, Bristol (1991).[2]
J. P. Boyd , Chebyshev and Fourier Spectral Methods, Dover, New York (2001).[3]
S. M. Cox, P. C. Matthews , Exponential Time Differencing for Stiff Systems , J. Comput.Phys. (2002), 430-455.[4]
F. Filbet, E. Sonnendr¨ucker , Modeling and numerical simulation of space charge dom-inated beams in the paraxial approximation , Math. Models Methods Appl. Sci. -5 (2006),763-791.[5] E. Fr´enod , Application of the averaging method to the gyrokinetic plasma , Asympt. Analysis (2006), 1-28.[6] E. Fr´enod, F. Salvarani, E. Sonnendr¨ucker , Long time simulation of a beam in aperiodic focusing channel via a two-scale PIC-method , Math. Models Methods Appl. Sci. -2 (2009), 175-197.[7] M. Hochbruck, A. Ostermann , Exponential integrators , Acta Numer. (2010), 209-286.[8] A.-K. Kassam, L. N. Trefethen , Fourth-order time-stepping for stiff PDEs , SIAM J.Sci. Comput. -4 (2005), 1214-1233.[9] T. T¨uckmantel, A. Pukhov, J. Liljo, M. Hochbruck , Three-Dimensional RelativisticParticle-in-Cell Hybrid Code Based on an Exponential Integrator , IEEE Trans. Plasma Sci. -9 (2010), 2383-2389. 14igure 8: First nonlinear case, ε = 10 − : phase space solution computed with a time step 8750 ε at final time 3.5, using for particles period 2 πε (at left) and the mean period (at right)Figure 9: Vlasov-Poisson case, ε = 10 − : phase space solution computed with a time step 8750 ε at final time 3.5, using for particles period 2 πεπε