A filtered Boris algorithm for charged-particle dynamics in a strong magnetic field
aa r X i v : . [ m a t h . NA ] J u l Version of 16 July 2019
A filtered Boris algorithm for charged-particle dynamicsin a strong magnetic field
Ernst Hairer , Christian Lubich ,Bin Wang Abstract
A modification of the standard Boris algorithm, called filtered
Borisalgorithm, is proposed for the numerical integration of the equations of motionof charged particles in a strong non-uniform magnetic field in the asymptoticscaling known as maximal ordering. With an appropriate choice of filters,second-order error bounds in the position and in the parallel velocity, andfirst-order error bounds in the normal velocity are obtained with respect tothe scaling parameter. The proof compares the modulated Fourier expansionsof the exact and the numerical solutions. Numerical experiments illustrate theerror behaviour of the filtered Boris algorithm.
Keywords.
Charged particle, magnetic field, guiding center, Boris algorithm,filter functions, exponential integrator, modulated Fourier expansion.
Mathematics Subject Classification (2010):
In this paper we propose and analyse a numerical integrator for the equationsof motion of a charged particle in a strong inhomogeneous magnetic field,¨ x ( t ) = ˙ x ( t ) × B ( x ( t ) , t ) + E ( x ( t ) , t )with B ( x, t ) = 1 ε B ( εx ) + B ( x, t ) for 0 < ε ≪ . (1.1) Dept. de Math´ematiques, Univ. de Gen`eve, CH-1211 Gen`eve 24, Switzerland,E-mail: [email protected] Mathematisches Institut, Univ. T¨ubingen, D-72076 T¨ubingen, Germany,E-mail: [email protected] School of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R.China,E-mail: [email protected] E. Hairer, Ch. Lubich, B. Wang
This scaling is of interest in particle methods in plasma physics and is called maximal ordering in [2]; see also [11] for a careful discussion of scalings anda rigorous analysis of this model. It is assumed that | B (0) | ≥
1, that B , B and E are smooth functions that are bounded independently of ε on boundeddomains together with all their derivatives, and that the initial data x (0) = x ,˙ x (0) = v are bounded independently of ε .In (1.1), x ( t ) ∈ R represents the position at time t of a charged particle(of unit mass and charge) that moves in the magnetic field B and the electricfield E . The motion is composed of fast rotation around a guiding center (withthe Larmor radius proportional to ε ) and slow motion of the guiding center.The standard integrator for charged particles in a magnetic field is the Boris algorithm [1] (see also, e.g., [5]), which in the two-step formulation withstep size h is given by x n +1 − x n + x n − h = x n +1 − x n − h × B ( x n , t n ) + E ( x n , t n ) (1.2)with the velocity approximation v n = h (cid:0) x n +1 − x n − (cid:1) at time t n = nh .This algorithm does, however, not behave well for (1.1) with small ε . Here wepropose a modification, which we name filtered Boris algorithm . This modifiedintegrator allows us to obtain better accuracy with considerably larger timesteps, at minor additional computational cost. It is still a symmetric algorithm.We formulate and discuss this new algorithm in Section 2. It comes in differentvariants that depend on the choice of a suitable filter function and of thepositions where the magnetic field is evaluated, and we identify favourablechoices.In Section 3 we state the main theoretical result, Theorem 3.1, which givesan error bound that behaves favourably with respect to ε . While most filtersonly yield a first-order error bound in the positions, for a particular, non-trivialchoice of the filter a second-order error bound is obtained. A second-order errorbound is also obtained for the component of the velocity that is parallel tothe magnetic field. For the normal velocity approximation, there is only afirst-order error bound for any filter. The proof of Theorem 3.1 is based oncomparing the modulated Fourier expansions of the exact and the numericalsolutions, which are derived in Sections 4 and 5, respectively. Combining thoseresults, the proof of Theorem 3.1 is finally completed in Section 6.We remark that the differential equations for the coefficient functions ofthe modulated Fourier expansions derived explicitly up to O ( ε ) in Section 4also yield the motion of the guiding center up to O ( ε ). They coincide up to O ( ε ) with the guiding center equations of the numerical approximation givenby the filtered Boris integrator for an appropriate filter and for non-resonantstep sizes h ≤ Cε with a possibly large constant C . This does not hold truefor the standard Boris integrator.In Section 7 we describe a related, but different integrator, called two-pointfiltered Boris algorithm, which evaluates the magnetic field both in the currentposition and in the current guiding center approximation in each step, and filtered Boris algorithm for charged-particle dynamics 3 which has similar convergence properties to the previously considered filteredBoris method.In Section 8 we present the results of numerical experiments in which wecompare the standard and filtered Boris algorithms.In the Appendix we show how the filters are evaluated efficiently using aRodriguez-type formula. Using the velocity approximation at the mid-point, v n − / = 1 h (cid:0) x n − x n − (cid:1) = v n − h v n × B ( x n , t n ) − h E ( x n , t n ) , the Boris algorithm (1.2) is usually written and implemented as a one-stepmethod ( x n , v n − / ) ( x n +1 , v n +1 / ), v n − / = v n − / + h E ( x n , t n ) v n +1 / − − v n − / = h (cid:0) v n +1 / − + v n − / (cid:1) × B ( x n , t n ) v n +1 / = v n +1 / − + h E ( x n , t n ) x n +1 = x n + h v n +1 / . (2.1)To capture the high oscillations in the velocity more accurately, the secondline of (2.1) needs to be modified, and one should rather work with averagedvelocities v n +1 / ≈ h R t n +1 t n v ( t ) d t and possibly averaged positions. This canbe achieved with the help of filter functions like in [3,8] and [7, Section XIII.2].For a vector B = ( b , b , b ) ⊤ ∈ R we denote by | B | the Euclidean normof B and we use the common notation v × B = − b B v, b B = − b b b − b − b b . (2.2)For real-analytic functions Ψ ( ζ ) (such as exp( ζ )) we will form matrix func-tions Ψ ( − h b B ), which are efficiently computed by a Rodriguez-type formula asdescribed in the Appendix.We denote by x n ⊙ = x n + v n × B n | B n | (2.3)with B n = B ( x n , t n ) the guiding center approximation at time t n (cf. [10]).For the argument of B in the algorithm we choose a point on the straight lineconnecting x n and x n ⊙ : ¯ x n = θ n x n + (1 − θ n ) x n ⊙ (2.4) E. Hairer, Ch. Lubich, B. Wang with θ n = θ ( h | B n | ) for a real function θ . It turns out that there is a uniquechoice of θ such that a second-order error bound will be obtained: θ ( ξ ) = 1sinc( ξ/ , (2.5)where sinc( ξ ) = sin( ξ ) /ξ . We note that with the scaling (1.1), we have ¯ x n = x n + O ( ε ), provided that h | B n | is bounded away from non-zero integral mul-tiples of 2 π .We consider the following modification of the Boris algorithm. Algorithm 2.1 (Filtered Boris algorithm)
Given ( x n , v n − / ) , the algo-rithm computes ( x n +1 , v n +1 / ) as follows, with B n = B ( x n , t n ) , ¯ B n = B (¯ x n , t n ) with ¯ x n defined by (2.4) , and E n = E ( x n , t n ) : v n − / = v n − / + h Ψ ( h c B n ) E n v n +1 / − = exp (cid:0) − h c ¯ B n (cid:1) v n − / .v n +1 / = v n +1 / − + h Ψ ( h c B n ) E n x n +1 = x n + h v n +1 / , (2.6) where Ψ ( ζ ) = tanch( ζ/ with tanch( ζ ) = tanh( ζ ) /ζ .The velocity approximation v n is computed as v n = Φ ( h c ¯ B n ) x n +1 − x n − h − hΥ ( h c B n ) E n , (2.7) where Φ ( ζ ) = 1sinch( ζ ) with sinch( ζ ) = sinh( ζ ) ζ , and Υ ( ζ ) = Φ ( ζ ) − ζ .The starting approximation v / is computed from (2.10) below with n = 0 . For the choice θ n = 1, the algorithm is explicit and requires only matrix-vector multiplications that can be done efficiently with a Rodriguez-type for-mula (see the Appendix).For θ n = θ ( h | B n | ) with θ ( ζ ) from (2.5), the algorithm is implicit, because¯ x n then depends on v n and appears in the argument of ¯ B n in the second line.This can be solved by a rapidly convergent fixed-point iteration for ¯ x n , withan error reduction by a factor O ( ε ) in each iteration. We start the iterationwith ¯ x n = x n , then compute v n +1 / − from (2.6) and v n from (2.7) using h ( x n +1 − x n − ) = (cid:0) v n +1 / + v n − / (cid:1) = (cid:0) v n +1 / − + v n − / (cid:1) . (2.8)This then yields x n ⊙ from (2.3) and the new ¯ x n from (2.4). In practice, oneiteration is sufficient to get the improved accuracy. We note that all matrix-vector multiplications can be done with a Rodriguez-type formula.We mention that Algorithm 2.1 preserves volume in phase space exactlyin the case of constant B (and time-dependent B ( t )), but only approximatelyup to O ( hε ) in the general case of an inhomogeneous magnetic field (1.1). filtered Boris algorithm for charged-particle dynamics 5 Two-step formulation.
The filtered Boris algorithm has the symmetric two-step formulation x n +1 − x n + x n − h = 2 h tanh (cid:0) − h c ¯ B n (cid:1) x n +1 − x n − h + Ψ ( h c B n ) E n , (2.9)as is readily obtained by taking two consecutive steps and using (2.8). Thisformulation is the basis of our theoretical analysis. Starting value.
The starting value v / is chosen such that formulas (2.6)-(2.7) also hold for n = 0. We find, for arbitrary n , that v n ± / = ϕ (cid:0) ∓ h c ¯ B n (cid:1)(cid:16) v n + hΥ ( h c B n ) E n (cid:17) ± h Ψ ( h c B n ) E n , (2.10)where ϕ ( ζ ) = ( e ζ − /ζ . Note that, for given x and v , the vectors x n ⊙ and¯ x n are known, so that (2.7) provides an explicit formula for v / . One-step map ( x n , v n ) ( x n +1 , v n +1 ) . Using the last formula of (2.6) to-gether with (2.10) for relating x n +1 and x n , and (2.10) with n and + andwith n + 1 and − for relating v n +1 and v n , the filtered Boris algorithm can bewritten as x n +1 = x n + hΦ n + v n + h Ψ n + E n Φ n +1 − v n +1 = Φ n + v n + h Ψ n +1 − E n +1 + h Ψ n + E n , (2.11)where Φ n ± = ϕ ( ∓ h c ¯ B n ) and Ψ n ± = Ψ ( h c B n ) ± Φ n ± Υ ( h c B n ).The method is symmetric in the sense that exchanging n ↔ n + 1 and h ↔ − h gives the same formulas. The integrator in the case of a constant magnetic field.
For constant B , we note that ( Φ n +1 − ) − Φ n + = exp( − h b B ), and so (2.11) reduces to the expo-nential integrator (with the notation Ψ ± ( ζ ) = Ψ ( ζ ) ∓ ϕ ( ± ζ ) Υ ( ζ )) x n +1 = x n + hϕ ( − h b B ) v n + h Ψ + ( − h b B ) E n v n +1 = exp( − h b B ) v n + h (cid:0) Ψ ( − h b B ) E n + Ψ ( − h b B ) E n +1 (cid:1) (2.12)with Ψ ( ζ ) = Ψ + ( ζ ) /ϕ ( − ζ ) and Ψ ( ζ ) = Ψ − ( ζ ) /ϕ ( − ζ ). The method is exactfor a constant magnetic field B and vanishing electric field E , becauseexp (cid:18) tI − t b B (cid:19) = I t ϕ ( − t b B )0 exp( − t b B ) ! . (2.13)Since we have chosen Ψ ( ζ ) = tanch( ζ/ exact for constant B and E . This is seen as follows: For constant B , the variation-of-constantsformula for the system ˙ x = v, ˙ v = x × B + E ( x ) reads, in view of (2.13), x ( t n + h ) = x ( t n ) + hϕ ( − h b B ) v ( t n )+ h R (1 − s ) ϕ ( − (1 − s ) h b B ) E ( x ( t n + hs ))d s,v ( t n + h ) = exp( − h b B ) v ( t n ) + h R exp( − (1 − s ) h b B ) E ( x ( t n + hs ))d s. For constant E , this becomes (2.12), which yields Ψ ± ( ζ ) = ϕ ( ± ζ ), where ϕ ( ζ ) = ( e ζ − − ζ ) / ( ζ /
2) = R (1 − s ) ϕ ((1 − s ) ζ )d s . E. Hairer, Ch. Lubich, B. Wang
Our main theoretical result in this paper is the following error bound for thefiltered Boris algorithm. Here we denote, for the exact velocity v ( t ) = ˙ x ( t ), v k ( t ) = B ( x ( t ) , t ) | B ( x ( t ) , t ) | (cid:18) B ( x ( t ) , t ) | B ( x ( t ) , t ) | · v ( t ) (cid:19) , v ⊥ ( t ) = v ( t ) − v k ( t ) , and similarly for the numerical velocity v n , v n k = B ( x n , t ) | B ( x n , t n ) | (cid:18) B ( x n , t n ) | B ( x n , t n ) | · v n (cid:19) , v n ⊥ = v n − v n k . We then have the following result.
Theorem 3.1
We assume the following, with arbitrarily chosen positive con-stants c , C , M and T :1. The initial velocity satisfies an ε -independent bound | v | ≤ M. (3.1)
2. The exact solution x ( t ) of (1.1) stays in a bounded set K (independentof ε ) for ≤ t ≤ T .3. The step size satisfies h ≤ Cε and is such that the following non-resonancecondition is satisfied: (cid:12)(cid:12) sinc (cid:0) kh | B ( x ( t ) , t ) | (cid:1)(cid:12)(cid:12) ≥ c > for k = 1 , , . (3.2) If in the filtered Boris algorithm, – ¯ x n is given by (2.4) with the function θ of (2.5) , and – the filter functions Ψ and Υ are defined as in Algorithm 2.1,then the errors in the positions and the velocities are bounded by x n − x ( t n ) = O ( ε ) ,v n k − v k ( t n ) = O ( ε ) , v n ⊥ − v ⊥ ( t n ) = O ( ε ) . (3.3) For a different choice of the functions θ , Ψ and Υ , the error bounds are notbetter than O ( ε ) for general problems (1.1) . The constants in the O -notationare independent of ε and h and n with ≤ t n = nh ≤ T , but depend on T , onthe velocity bound M and the constants c and C , and on bounds of derivativesof B , B and E in a neighbourhood of the set K . We remark that in view of the error bounds, the non-resonance conditionmight be required along the numerical solution x n instead of the exact solution x ( t ) as in (3.2).The proof of this theorem will compare the modulated Fourier expansionof the exact solution (as given in Section 4) with that of the numerical ap-proximation (as given in Section 5). It will be given in Section 6. Remark 3.1
The proof also shows that the choice ¯ x n = x n is sufficient fororder 2 if the magnetic field satisfies, for all z ∈ C and x ∈ K and all times t ,Im ( z × ∂ x B ( x, t )¯ z ) · B ( x, t ) = O ( ε ) . filtered Boris algorithm for charged-particle dynamics 7 We write the solution of (1.1) as a modulated Fourier expansion x ( t ) ≈ X k ∈ Z z k ( t ) e i kφ ( t ) /ε (4.1)with coefficient functions z k ( t ) for which all time derivatives are bounded in-dependently of ε , where ˙ φ ( t ) /ε = (cid:12)(cid:12) B (cid:0) z ( t ) , t (cid:1)(cid:12)(cid:12) , and z ( t ) describes the motionof the guiding center. Such a formal expansion has first been considered in[9] for proving the existence of an adiabatic invariant (essentially the mag-netic moment | ˙ x × B ( x ) | / | B ( x ) | ). It has been used for a rigorous proofof the long-time near-conservation of the magnetic moment in [6], where thisapproach was extended to the numerical solution of a variational integrator,for which near-conservation of the magnetic moment and of the energy is rig-orously proved over long times that cover arbitrary negative powers of ε .Following [6], we diagonalize the linear map v v × B ( x, t ), which haseigenvalues λ = i | B ( x, t ) | , λ = 0, and λ − = − i | B ( x, t ) | . We denote thenormalized eigenvectors by v ( x, t ) , v ( x, t ) , v − ( x, t ), and remark that v ( x, t )is collinear to B ( x, t ). We let P j ( x, t ) = v j ( x, t ) v j ( x, t ) ∗ be the orthogonal pro-jections onto the eigenspaces. Furthermore, we write the coefficient functionsof (4.1) in the time-dependent basis v j (cid:0) z ( t ) , t (cid:1) , z k = z k + z k + z k − , z kj ( t ) = P j (cid:0) z ( t ) , t (cid:1) z k ( t ) . (4.2)Since x ( t ) is real, we assume z − k = z k for all k . Together with the fact that v − ( x, t ) = v ( x, t ) and v ( x, t ) is real, it follows z − k − = z k , z − k = z k , z − k = z k − . (4.3)The following result is a variant of Theorem 4.1 in [6], adapted to the presentcase of a strong magnetic field of the form (1.1). Note that B in this papercorresponds to B/ε in [6].
Theorem 4.1
Let x ( t ) be a solution of (1.1) with bounded initial velocity (3.1) that stays in a compact set K for ≤ t ≤ T . For an arbitrary truncationindex N ≥ we then have x ( t ) = X | k |≤ N z k ( t ) e i kφ ( t ) /ε + R N ( t ) , (4.4) where the phase function satisfies ˙ φ ( t ) = ε | B ( z ( t ) , t ) | = O (1) .(a) The coefficient functions z k ( t ) together with their derivatives (up toorder N ) are bounded as z j = O (1) for j ∈ {− , , } , z = O ( ε ) , z − − = O ( ε ) , z kj = O ( ε ) for | k | = 1 , j = k , and for the remaining ( j, k ) with | k | ≤ N , z kj = O ( ε | k | +1 ) . (4.5) They are unique up to O ( ε N +2 ) . Moreover, we have ˙ z × B ( z , t ) = O (1) . E. Hairer, Ch. Lubich, B. Wang (b) The remainder term and its derivative are bounded by R N ( t ) = O ( t ε N ) , ˙ R N ( t ) = O ( tε N ) for ≤ t ≤ T. (4.6) (c) The functions z , z ± , z , z − − satisfy the differential equations ¨ z = P ( z , t ) E ( z , t ) + 2 P ( z , t ) Re (cid:16) i ˙ φε z × B ′ ( z , t ) z − − (cid:17) + 2 ˙ P ( z , t ) ˙ z + ¨ P ( z , t ) z + O ( ε ) , (4.7)˙ z ± = ˙ P ± ( z , t ) z ± i ε ˙ φ P ± ( z , t ) E ( z , t ) + O ( ε ) , (4.8)˙ z ± ± = − ¨ φ ˙ φ z ± ± + O ( ε ) = O ( ε ) , (4.9) where we use the notation ˙ P j ( z , t ) = dd t P j (cid:0) z ( t ) , t (cid:1) and similar for ¨ P j ( z , t ) .All other coefficient functions z kj are given by algebraic expressions dependingon z , ˙ z , z , z − − .(d) Assuming φ (0) = 0 , initial values for the differential equations ofitem (c) are given by z (0) = x (0) + ˙ x (0) × B ( x (0) , | B ( x (0) , | + O ( ε ) , ˙ z (0) = P ( x (0) ,
0) ˙ x (0) + ˙ P ( x (0) , x (0) + O ( ε ) ,z ± ± (0) = ∓ i ε ˙ φ (0) P ± ( x (0) ,
0) ˙ x (0) + O ( ε ) . The constants symbolised by the O -notation are independent of ε and t with ≤ t ≤ T , but they depend on N , on the velocity bound M in (3.1) , on boundsof derivatives of B and E , and on T .Remark 4.1 We note that the guiding center motion of the system (1.1) isgiven by the non-oscillating term z ( t ) in the modulated Fourier expansion.By the uniqueness of the modulated Fourier expansion up to high powers of ε ,the equations in item (d) hold not only at time 0, but for all t ≤ T . Proof (a) and (b): Compared to Theorem 4.1 in [6], where a more generalstrong magnetic field is considered, the time interval of validity of the modu-lated Fourier expansion is here O (1) instead of just O ( ε ), and the bound (4.5)is improved by a factor ε . The improvement of the time scale comes about byobserving that a function x ∗ ( t ) that solves (1.1) up to a defect d ( t ), i.e.,¨ x ∗ ( t ) = ˙ x ∗ ( t ) × B ( x ∗ ( t ) , t ) + E ( x ∗ ( t ) , t ) + d ( t ) , satisfies an error bound, for 0 ≤ t ≤ T , | x ∗ ( t ) − x ( t ) | ≤ C (cid:16) | x ∗ (0) − x (0) | + | ˙ x ∗ (0) − ˙ x (0) | + Z t | d ( t ) | d t (cid:17) , filtered Boris algorithm for charged-particle dynamics 9 where C is independent of ε but grows exponentially with T . This is proved bydecomposing B ( x, t ) = ε − B (0)+ ε − ( B ( εx ) − B (0))+ B ( x, t ) and using thevariation-of-constants formula and the Gronwall inequality. The improvementof the bound (4.5) is a consequence of the fact that the derivatives of B ( x, t )are bounded independently of ε .(c): For the error bound of Section 6 we need precise formulas for thedominant terms of (4.4). Inserting the expansion (4.1) into the differentialequation (1.1) and comparing the coefficients of e i kφ ( t ) /ε yields¨ z k + 2i k ˙ φε ˙ z k + (cid:16) i k ¨ φε − k ˙ φ ε (cid:17) z k = F k , (4.10)where, using Taylor series expansion for the nonlinearities, F k = X k + k = k (cid:16) ˙ z k + i k ˙ φε z k (cid:17) × X m ≥ s ( α )= k m ! B ( m ) ( z , t ) z α + X m ≥ s ( α )= k m ! E ( m ) ( z , t ) z α . Here, B ( m ) ( x, t ) and E ( m ) ( x, t ) denote the m th derivative with respect to x , α = ( α , . . . , α m ) is a multi-index with α j ∈ Z \ { } , s ( α ) = α + . . . + α m , | α | = | α | + . . . + | α m | , and z α = ( z α , . . . , z α m ).From (4.10) it follows that the motion of the guiding center z ( t ) is givenby¨ z = ˙ z × B ( z , t ) + E ( z , t ) + 2 Re (cid:16) i ˙ φε z × B ′ ( z , t ) z − (cid:17) + O ( ε ) . (4.11)The solution z ( t ) is influenced by the functions z ± which, by (4.10), satisfy ±
2i ˙ φε ˙ z ± + (cid:16) ± i ¨ φε − ˙ φ ε (cid:17) z ± = (cid:16) ˙ z ± ± i ˙ φε z ± (cid:17) × B ( z , t ) + O ( ε ) . (4.12)Note that, whereas B ( z , t ) is of size O ( ε − ), its derivatives are bounded in-dependently of ε due to the special form (1.1).To get solutions with derivatives bounded uniformly in ε , one has to extractthe dominant terms. Multiplying (4.11) with P ( z , t ) eliminates the ε − -termthat is present in B ( z , t ), and the second derivative ¨ z becomes dominant.Differentiating the relation z = P ( z , t ) z with respect to time yields ¨ z = P ( z , t )¨ z + 2 ˙ P ( z , t ) ˙ z + ¨ P ( z , t ) z . This then gives (4.7). Note that, dueto the special form of B ( x, t ), the time derivatives of P j ( z , t ) are of size O ( ε ).A multiplication of (4.11) with P ± ( z , t ) gives P ± ( z , t )¨ z = ± i ˙ φε P ± ( z , t ) ˙ z + P ± ( z , t ) E ( z , t ) + O ( ε ) . Substituting P ± ( z , t ) ˙ z = ˙ z ± − ˙ P ± ( z , t ) z , and extracting ˙ z ± yields (4.8).Note that ˙ z ± = O ( ε ), so that also ¨ z ± = O ( ε ), and P ± ( z , t )¨ z = O ( ε ).Since ˙ φ/ε = | B ( z , t ) | , the ε − -terms cancel in (4.12) after projection with P ± ( z , t ). Therefore, the ε − -terms are dominant and we obtain (4.9). (d): Assuming φ (0) = 0, initial values are determined from (4.4) by x (0) = z (0) + (cid:0) z (0) + z − (0) (cid:1) + O ( ε )˙ x (0) = ˙ z (0) + (cid:0) ˙ z (0) + ˙ z − (0) (cid:1) + i ˙ φ (0) ε (cid:0) z (0) − z − (0) (cid:1) + O ( ε ) . (4.13)This is a nonlinear system for z (0) , ˙ z (0) , z (0) , z − − (0). We write the vectorsin the basis { v j ( z (0) , } , and we select the dominant terms in each equation.They are z (0) in the upper relation of (4.13), and ˙ z (0) , z (0) , z − − (0) in thelower relation. Fixed-point iteration then yields the stated equations for theinitial values. Note that the relation P j ( z (0) ,
0) = P j ( x (0) ,
0) + O ( ε ) hasbeen applied. ⊓⊔ We consider the two-step formulation (2.9) of the filtered Boris algorithm, andwe write the numerical approximation x n as x n ≈ X k ∈ Z z k ( t ) e i kφ ( t ) /ε , t = nh. (5.1)We use the same notation for the coefficient functions as in Section 4. Note,however, that for the numerical solution these functions are not the same anddepend on the additional parameter h . We again consider the basis { v j ( x, t ) } and the corresponding orthogonal projections P j ( x, t ), and we write the coeffi-cient functions z k as in (4.2), with the only difference that here the argument z ( t ) is the non-oscillating part of (5.1) and not that of (4.1). Theorem 5.1
Let { x n } be a numerical solution of the filtered Boris algorithmapplied to (1.1) with bounded initial velocity (3.1) , and suppose that it staysin a compact set K for ≤ nh ≤ T . We assume the non-resonance condition (cid:12)(cid:12) sinc (cid:0) kh | B ( x n , t n ) | (cid:1)(cid:12)(cid:12) ≥ c > for k = 1 , . . . , N + 1 , (5.2) for a fixed, but arbitrary truncation index N ≥ , and (for convenience ofpresentation) also the bound η = h/ε ≤ C . Moreover, we assume that the filterfunction Ψ in (2.9) is bounded by | Ψ (i ξ ) | ≤ C | tanc( ξ ) | for all real ξ , where tanc( ξ ) = tan( ξ ) /ξ . Then, we have that x n = X | k |≤ N z k ( t ) e i kφ ( t ) /ε + R N ( t ) , t = nh, (5.3) where the phase function is given by ˙ φ ( t ) = ε | B ( z ( t ) , t ) | .(a) and (b) The coefficient functions z k ( t ) together with their derivatives(up to order N ) as well as the remainder term and its derivative satisfy thebounds of items (a) and (b) of Theorem 4.1. filtered Boris algorithm for charged-particle dynamics 11 (c) The functions z , z ± , z , z − − satisfy the differential equations (with θ ( ξ ) used in the definition of ¯ x n in (2.4) ) ¨ z = P ( z , t ) E ( z , t ) + 2 P ( z , t ) Re (cid:16) i ˙ φε z × B ′ ( z , t ) z − − (cid:17) θ ( η ˙ φ ) sinc( η ˙ φ/ + 2 ˙ P ( z , t ) ˙ z + ¨ P ( z , t ) z + O ( ε ) , (5.4)˙ z ± = ˙ P ± ( z , t ) z ± Ψ (cid:0) i η ˙ φ (cid:1) tanc (cid:16) η ˙ φ (cid:17) i ε ˙ φ P ± ( z , t ) E ( z , t ) + O ( ε ) , (5.5)˙ z ± ± = − (cid:16) η ˙ φ (cid:17) ¨ φ ˙ φ z ± ± + O ( ε ) = O ( ε ) . (5.6) All other coefficient functions z kj are given by algebraic expressions depend-ing on z , ˙ z , z , z − − .(d) Assuming φ (0) = 0 , initial values for the differential equations ofitem (c) are given by the same equations as for the exact solution, up to O ( ε ) , z (0) = x (0) + ˙ x (0) × B ( x (0) , | B ( x (0) , | + O ( ε ) , ˙ z (0) = P ( x (0) ,
0) ˙ x (0) + ˙ P ( x (0) , x (0) + O ( ε ) , (5.7) z ± ± (0) = ∓ i ε ˙ φ (0) P ± ( x (0) ,
0) ˙ x (0) + O ( ε ) . The constants symbolised by the O -notation are independent of ε and n with ≤ nh ≤ T , but they depend on N , on the velocity bound M in (3.1) , onbounds of derivatives of B and E , and on T .Proof (a) and (b) We do not present the details of the proof of the existence ofthe modulated Fourier expansion and the bounds for the coefficient functionsand the remainder term, since this uses the same kind of arguments as inprevious such proofs, e.g. in [7,4,6]. In particular, for | k | = 1 , j = k and for | k | ≥ z kj is multiplied by 4 η sin (cid:16) kη ˙ φ (cid:17) sin (cid:16) ( k − j ) η ˙ φ (cid:17) . Under the non-resonance assumption (5.2) this expression is bounded frombelow by a positive constant, so that an algebraic relation for z kj can be ex-tracted.By construction of the coefficient functions the truncated series of (5.3)satisfies the two-step relation (2.9) up to a defect of size O ( ε N ). A standarddiscrete Gronwall argument then gives the bounds on the remainder.(d) The initial values are obtained from x (0) = z (0) + (cid:0) z (0) + z − (0) (cid:1) + O ( ε ) , (5.8) which is a consequence of (5.1), and from˙ x (0) = Φ (cid:0) h b B ( x (0) , (cid:1) ˙ z (0) + i ˙ φ (0) ε z (0) − i ˙ φ (0) ε z − − (0) − hΥ (cid:0) h b B ( x (0) , (cid:1) E ( x (0) ,
0) + O ( ε ) , (5.9)which follows from (2.7) and Lemma 5.1. As in the proof of Theorem 4.1this constitutes a nonlinear system for the values z (0) , ˙ z (0) , z (0) , z − − (0).The relation (5.8) yields z (0). Multiplication of (5.9) with P j (cid:0) z (0) , (cid:1) = P j (cid:0) x (0) , (cid:1) + O ( ε ) gives ˙ z (0) for j = 0 and z ± ± (0) for j = ±
1, where weuse in addition that Φ ( h b B ( x (0) , Φ ( h b B ( z (0) , O ( ε ) and P ± ˙ z =˙ z ± − ˙ P ± z = O ( ε ). Remarkably we get, up to terms of size O ( ε ), the sameformulas for the initial values as for the exact solution.By the uniqueness of the modulated Fourier expansion (up to O ( ε N )), theserelations hold not only at time 0, but for arbitrary times t ≤ T , except for aphase factor e ∓ i φ ( t ) in the equation for z ± ± . This phase factor did not appearin (5.7) because we chose φ (0) = 0.(c) To derive the differential equations for the coefficient functions we firstexpand the perturbed argument of B ( x, t ) in the filtered Boris algorithm as¯ x n ≈ X k ∈ Z ζ k ( t ) e i kφ ( t ) /ε , t = nh. (5.10)The coefficient functions ζ k ( t ) are obtained as follows: inserting the modu-lated Fourier expansion (5.1) into (2.7), using Lemma 5.1 below, and replacing Φ (cid:0) h c ¯ B n (cid:1) by Φ (cid:0) h b B ( z ( t n ) , t n ) (cid:1) yields with t = nhv n = ˙ z ( t ) + i ˙ φ ( t ) ε z ( t ) e i φ ( t ) /ε − i ˙ φ ( t ) ε z − − ( t ) e − i φ ( t ) /ε + O ( ε ) , see also the more detailed computation in Section 6. Since we have ˙ z ( t ) = P (cid:0) z ( t ) , t (cid:1) z ( t ) + O ( ε ) and B n = b B (cid:0) z ( t n ) , t n (cid:1) + O ( ε ), this implies v n × B n | B n | = − z ( t ) e i φ ( t ) /ε − z − − ( t ) e − i φ ( t ) /ε + O ( ε ) , and consequently x n ⊙ = z ( t n )+ O ( ε ), which shows that x n ⊙ is an excellent ap-proximation of the non-oscillating part of the numerical solution x n . Togetherwith the definition (2.4) of ¯ x n we find the dominating terms of the expansion(5.10) as ζ ( t ) = z ( t ) + O ( ε ) , ζ ± ( t ) = θ ( t ) z ± ( t ) + O ( ε ) , (5.11)where θ ( t ) = θ ( h ˙ φ ( t ) /ε ) = θ (cid:0) h | B ( z ( t ) , t ) | (cid:1) . filtered Boris algorithm for charged-particle dynamics 13 After this preparation, we insert the expansions (5.1) for x n and (5.10)for ¯ x n into the two-step formulation (2.9) of the filtered Boris algorithm. Us-ing Lemma 5.1 below, expanding the nonlinearities around ζ and z , andcomparing the coefficients of e i kφ ( t ) /ε yields X l ≥ ε l − d kl d l d t l z k = X k + k = k (cid:16) X m ≥ s ( α )= k m ! T ( m ) b B ( ζ , t ) ζ α (cid:17)(cid:16) X l ≥ ε l − c k l d l d t l z k (cid:17) + X k + k = k (cid:16) X m ≥ s ( α )= k m ! Ψ ( m ) b B ( z , t ) z α (cid:17)(cid:16) X m ≥ s ( α )= k m ! E ( m ) ( z , t ) z α (cid:17) , where T ( m ) b B ( x, t ) denotes the m th derivative of T b B ( x, t ) = h tanh (cid:0) − h b B ( x, t ) (cid:1) with respect to x and, similarly, Ψ ( m ) b B ( x, t ) is the m th derivative of Ψ b B ( x, t ) = Ψ (cid:0) − h b B ( x, t ) (cid:1) with respect to x . These derivatives are bounded under theassumption that η = h/ε ≤ c and the non-resonance condition (3.2).For k = 0 we obtain¨ z = T b B ( ζ , t ) ˙ z + Ψ b B ( z , t ) E ( z , t )+ 2 Re (cid:16)(cid:0) T ′ b B ( ζ , t ) ζ − (cid:1) i εη sin( η ˙ φ ) z (cid:17) + O ( ε ) , (5.12)and for k = ± ε − d ± z ± + ε − d ± ˙ z ± = T b B ( ζ , t ) (cid:0) ε − c ± z ± + c ± ˙ z ± (cid:1) + O ( ε ) . (5.13)Because of (5.11), the argument ζ can be replaced by z in these equations.In the limit h →
0, i.e., η → P ( z , t ) T b B ( z , t ) = 0 , P ± ( z , t ) T b B ( z , t ) = ± i 2 h tan (cid:16) h ˙ φ ε (cid:17) P ± ( z , t ) ,P ( z , t ) Ψ b B ( z , t ) = P ( z , t ) , P ± ( z , t ) Ψ b B ( z , t ) = Ψ (cid:0) ± i h ˙ φε (cid:17) P ± ( z , t ) . Multiplying the equation (5.12) with P ( z , t ) and applying the differentiationformula of Lemma 5.2 yields the differential equation P ( z , t ) ¨ z = P ( z , t ) E ( z , t )+ 2 P ( z , t ) Re (cid:18) − η ˙ φ tan (cid:16) η ˙ φ (cid:17)(cid:16) b B ′ ( z , t ) ζ − − (cid:17) i ˙ φε sinc( η ˙ φ ) z (cid:19) ) + O ( ε ) . Using (cid:0) b B ′ ( x, t ) ∆x (cid:1) v = − v × B ′ ( x, t ) ∆x , which follows from differentiation of b B ( x, t ) v = − v × B ( x, t ), the trigonometric identity sin(2 α ) = 2 sin( α ) cos( α ),and the second relation of (5.11), this equation becomes (5.4). A multiplication of (5.12) with P ± ( z , t ) permits to extract the dominantfirst derivative ˙ z ± and gives (5.5).We next consider the equation (5.13). The ε − -terms in the left and rightsides are contained in − ε η sin (cid:16) η ˙ φ (cid:17) z ± and ± εh tanh (cid:16) − h b B ( z , t ) (cid:17) i η sin( η ˙ φ ) z ± . After multiplication with P ± ( z , t ) these terms cancel because of the aboveformula for P ± ( z , t ) T b B ( z , t ). The remaining terms lead to˙ z ± ± = − (cid:16) cos( η ˙ φ ) + tan (cid:16) η ˙ φ (cid:17) sin( η ˙ φ ) (cid:17) ¨ φ η (cid:16) sin( η ˙ φ ) − tan (cid:16) η ˙ φ (cid:17) cos( η ˙ φ ) (cid:17) z ± ± + O ( ε ) , which simplifies to (5.6). ⊓⊔ In the above proof we referred to the following lemmas.
Lemma 5.1 ([6])
For smooth functions φ ( t ) and z k ( t ) let y k ( t ) = e i kφ ( t ) /ε z k ( t ) ,and denote η = h/ε . The finite differences of y k ( t ) then satisfy δ h y k ( t ) = y k ( t + h ) − y k ( t − h )2 h = e i kφ ( t ) /ε X l ≥ ε l − c kl ( t ) d l d t l z k ( t ) δ h y k ( t ) = y k ( t + h ) − y k ( t ) + y k ( t − h ) h = e i kφ ( t ) /ε X l ≥ ε l − d kl ( t ) d l d t l z k ( t ) , where c j = 0 , c j +1 = η j / (2 j + 1)! , and d = 0 , d j = 2 η j − / (2 j )! , d j +1 =0 . The leading coefficients are c k ( t ) = i η sin (cid:0) kη ˙ φ ( t ) (cid:1) − ε kη (cid:0) kη ˙ φ ( t ) (cid:1) ¨ φ ( t ) + O ( ε ) c k ( t ) = cos (cid:0) kη ˙ φ ( t ) (cid:1) + O ( ε ) d k ( t ) = − η sin (cid:16) kη ˙ φ ( t )2 (cid:17) + i ε k cos (cid:0) kη ˙ φ ( t ) (cid:1) ¨ φ ( t ) + O ( ε ) d k ( t ) = 2 i η sin (cid:0) kη ˙ φ ( t ) (cid:1) + O ( ε ) . (5.14)Note that these coefficients depend on η , ε , and t via derivatives of φ ( t ). Proof
Expanding φ ( t ± h ) and z k ( t ± h ) into Taylor series around t yields thestated formulas. ⊓⊔ filtered Boris algorithm for charged-particle dynamics 15 Lemma 5.2
Let T b B ( x, t ) = h tanh (cid:0) − h b B ( x, t ) (cid:1) , and let P j ( x, t ) be the orthog-onal projections onto the eigenspace of b B ( x, t ) corresponding to the eigenvalues λ = 0 and λ = − i | B ( x, t ) | = − i ˙ φ ( x, t ) /ε , and λ − = i | B ( x, t ) | = i ˙ φ ( x, t ) /ε ,respectively. Omitting the argument ( x, t ) , we then have with η = h/ε , P (cid:16) T ′ b B ∆x (cid:17) P ± = ∓ η ˙ φ tan (cid:16) η ˙ φ (cid:17) P (cid:16) b B ′ ∆x (cid:17) P ± , where prime indicates the derivative with respect to x .Proof Writing tanh as a Taylor series with coefficients γ l and differentiatingterm by term, we obtain P (cid:16) T ′ b B ∆x (cid:17) P ± = 2 h X l ≥ γ l (cid:16) − h (cid:17) l P (cid:16) b B ′ ∆x (cid:17) b B l − P ± = 2 h X l ≥ γ l (cid:16) − h (cid:17) l P (cid:16) b B ′ ∆x (cid:17)(cid:16) ∓ i ˙ φε (cid:17) l − P ± = 2i η ˙ φ tanh (cid:16) ± i η ˙ φ (cid:17) P (cid:16) b B ′ ∆x (cid:17) P ± . This proves the statement of the lemma. ⊓⊔ Theorems 4.1 and 5.1 show that the coefficient functions z k ( t ) (and also ˙ z ( t ))of the modulated Fourier expansions of the exact and numerical solutionscoincide up to O ( ε ) for the choice (2.5) and Ψ ( ζ ) = tanch( ζ/ φ (with ˙ φ ( t ) = ε | B ( z ( t ) , t ) | ) differ only by O ( ε ), respectively. Since all coefficient functions z k of the modulated Fourierexpansion with the exception of z are of size O ( ε ) or smaller, this yields thatall summands z k ( t )e i kφ ( t ) /ε still differ only by O ( ε ). So we obtain the O ( ε )error bound for the positions as stated in Theorem 3.1.We now turn to the error bound for the velocities. By Theorem 4.1, usingthat ˙ z ± ± = O ( ε ) and z kj = O ( ε ) for | k | = 1 and k = j and for | k | ≥ j = − , ,
1, together with their derivatives, the velocity of the exact solutionsatisfies v ( t ) = ˙ x ( t ) = ˙ z ( t ) + i ˙ φ ( t ) ε z ( t ) e i φ ( t ) /ε − i ˙ φ ( t ) ε z − − ( t ) e − i φ ( t ) /ε + O ( ε ) . (6.1)We shall show below that the numerical solution admits the same expan-sion with functions φ ( t ) , ˙ z ( t ), z ( t ), z − − ( t ) that correspond to the modulatedFourier expansion (5.1) of the numerical solution and not to (4.1) of the exactsolution. By Theorems 4.1 and 5.1, these functions differ only by O ( ε ). Be-cause of the denominator ε in the second and third terms on the right-handside of (6.1), this yields v n − v ( t n ) = O ( ε ) , but v n k − v k ( t n ) = O ( ε ) , (6.2) and proves the statement of Theorem 3.1.Using Lemma 5.1 and ¨ φ ( t ) = O ( ε ), we have, with t = nh , that x n +1 − x n − h = ˙ z ( t )+sinc (cid:0) η ˙ φ ( t ) (cid:1) i ˙ φ ( t ) ε (cid:16) z ( t )e i φ ( t ) /ε − z − − ( t )e − i φ ( t ) /ε (cid:17) + O ( ε ) . A consequence of the maximal ordering in (1.1) is that Φ (cid:0) h b B (¯ x n , t n ) (cid:1) = Φ (cid:0) h b B ( z ( t n ) , t n ) (cid:1) + O ( ε ), and Υ (cid:0) h b B ( x n , t n ) (cid:1) = Υ (cid:0) h b B ( z ( t n ) , t n ) (cid:1) + O ( ε ).Splitting Φ ( · ) ˙ z into ˙ z + (cid:0) Φ ( · ) − I (cid:1) ˙ z and using Υ ( ζ ) = (cid:0) Φ ( ζ ) − (cid:1) /ζ , wetherefore have v n = Φ (cid:0) h b B (¯ x n , t n ) (cid:1) x n +1 − x n − h − hΥ (cid:0) h b B ( x n , t n ) (cid:1) E ( x n , t n )= ˙ z ( t ) + i ˙ φ ( t ) ε (cid:16) z ( t )e i φ ( t ) /ε − z − − ( t )e − i φ ( t ) /ε (cid:17) + hΥ (cid:0) h b B ( z ( t ) , t ) (cid:1)(cid:16) b B (cid:0) z ( t ) , t (cid:1) ˙ z ( t ) − E (cid:0) z ( t ) , t (cid:1)(cid:17) + O ( ε ) . Since Υ (0) = 0 we have Υ (cid:0) h b B ( z ( t ) , t ) (cid:1) P ( z ( t ) , t ) = 0. On the other hand P ± ( z ( t ) , t ) (cid:16) b B (cid:0) z ( t ) , t (cid:1) ˙ z ( t ) − E (cid:0) z ( t ) , t (cid:1)(cid:17) = O ( ε )which follows from (5.5) for Ψ (i y ) = tanc( y/ Algorithm 2.1 evaluates the magnetic field B at ¯ x n given by (2.4)–(2.5), whichcan be far from both x n and the guiding center approximation x n ⊙ of (2.3)when h | B ( x n ) | is close to a nonzero integral multiple of 2 π . In the followingwe propose an alternative filtered Boris algorithm with the same second-orderconvergence properties as in Theorem 3.1, which evaluates the magnetic fieldat the two points x n and x n ⊙ . Algorithm 7.1 (Two-point filtered Boris algorithm)
Given ( x n , v n − / ) ,the algorithm computes ( x n +1 , v n +1 / ) as follows, with B n = B ( x n , t n ) , B n ⊙ = B ( x n ⊙ , t n ) and E n = E ( x n , t n ) : v n − / = v n − / + h Ψ ( h b B n ) E n Φ ( h b B n ⊙ )( v n +1 / − − v n − / ) = h Φ ( h b B n ) (cid:0) v n +1 / − + v n − / (cid:1) × B n v n +1 / = v n +1 / − + h Ψ ( h b B n ) E n x n +1 = x n + h v n +1 / , (7.1) filtered Boris algorithm for charged-particle dynamics 17 where Ψ ( ζ ) = tanch( ζ/ and Φ ( ζ ) = 1sinch( ζ ) are as in Algorithm 2.1, and Φ ( ζ ) = 1sinch( ζ/ .The velocity approximation v n is again computed by (2.7) , with B n insteadof ¯ B n . For constant B , Algorithms 2.1 and 7.1 are identical and explicit. In thegeneral case, both methods are implicit, but this time the fixed-point itera-tion for x n ⊙ requires not only the evaluation of matrix functions by the Ro-driguez formula, but in addition the solution of a linear system with the3 × Φ ( h b B n ⊙ ) + h b B n Φ ( h b B n ). We further note that in the caseof a vanishing electric field, E n = 0, Algorithm 2.1 preserves the velocitynorm | v n +1 / | = | v n − / | , which is satisfied only approximately up to O ( hε )by Algorithm 7.1. While these properties are unfavourable for Algorithm 7.1,our numerical experiments indicate that it yields higher accuracy than Al-gorithm 7.1 for stepsizes such that h | B | is large, and in particular it is lesssensitive to near-resonances where h | B | is close to an integral multiple of 2 π .The two-step formulation of Algorithm 7.1 is x n +1 − x n + x n − h (7.2)= Φ ( h b B (¯ x n , t n )) − (cid:16) Φ ( h b B n ) x n +1 − x n − h × B n (cid:17) + Ψ ( h b B n ) E n . The starting value v / is chosen such that formulas (7.1) and (2.7) alsohold for n = 0. With the abbreviations Λ n = Φ ( h b B n ⊙ ) − Φ ( h b B n ) , Ψ n = Ψ ( h b B n ) , Υ n = Υ ( h b B n ) ,Φ n ± = ( I ∓ Λ n h b B n ) sinch( h b B n ) ,Ψ n ± = Ψ n ± Φ n ± Υ n , we find, for arbitrary n , that like in (2.10), v n ± / = Φ n ± v n ± h Ψ n ± E n , and the one-step map ( x n , v n ) ( x n +1 , v n +1 ) is then again given by (2.11)with these modified matrices Φ n ± and Ψ n ± .For the two-point filtered Boris algorithm, the second-order convergenceresult of Theorem 3.1 in x and v k and the first-order convergence in v ⊥ remainvalid, as can be shown by an adaptation of the proof of Theorem 5.1, for whichwe omit the details. As an illustrative numerical experiment, we consider the charged-particle mo-tion in the magnetic field B ( x, t ) = ∇ × ε x + ∇ × x x = 1 ε + − x x , and the electric field E ( x, t ) = −∇ x U ( x ) with the potential U ( x ) = 1 p x + x . The initial values are chosen as x (0) = ( , , ) ⊺ and v (0) = ( , , ⊺ . Wesolve this problem for 0 ≤ t ≤ h = ǫ, ǫ, ǫ and compare the numericalerrors of the following methods: – the standard Boris algorithm, – Exp-A: the filtered Boris method of Algorithm 2.1 with θ = 1 in (2.4)(where ¯ x n = x n and the method is explicit), – Imp-A: the filtered Boris method of Algorithm 2.1 with θ of (2.5), – Two P-A: the two-point filtered Boris method of Algorithm 7.1.The errors in x and v k , v ⊥ against different ǫ = 1 / j are displayed in Fig. 8.1,where j = 4 , . . . ,
13. Then we fix ǫ = 1 / and show the errors at t = 1against h/ǫ in Fig. 8.2. It is observed that all three filtered Boris methods im-prove considerably over the standard Boris method, and the optimally filteredmethods Imp-A and Two P-A show second order, whereas method Exp-A onlyshows first order. Methods Exp-A and Two P-A behave very similar away fromstepsize resonances, but method Two P-A appears more robust near stepsizeresonances. For the implicit methods Imp-A and Two P-A, the error behaviourremains essentially unchanged after just one fixed-point iteration. Appendix: Implementation
The filtered Boris algorithm requires the computation of matrix functionsapplied to a vector. This can be done very efficiently with a Rodriguez-likeformula. Consider a vector B = ( b , b , b ) ⊤ ∈ R and the skew-symmetricmatrix b B of (2.2), and let b = | B | . Assume that the function ϕ ( ζ ) can beexpanded into a Taylor series at the origin with real coefficients c n , and write ϕ (i y ) = ϕ (0) + i yϕ ( y ) − y ϕ ( y )with ϕ ( y ) = P j ≥ c j +1 ( − y ) j and ϕ ( y ) = P j ≥ c j +2 ( − y ) j . The fact that b B = − b b B filtered Boris algorithm for charged-particle dynamics 19 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -10-8-6-4-20 l og ( G E ) The global errors of x with h= ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -10-8-6-4-20 l og ( G E ) The global errors of v || with h= ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -10-8-6-4-20 l og ( G E ) The global errors of v ⊥ with h= ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -8-6-4-202 l og ( G E ) The global errors of x with h=4 ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -8-6-4-202 l og ( G E ) The global errors of v || with h=4 ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -8-6-4-202 l og ( G E ) The global errors of v ⊥ with h=4 ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -8-6-4-202 l og ( G E ) The global errors of x with h=16 ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -8-6-4-202 l og ( G E ) The global errors of v || with h=16 ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1 -4 -3.5 -3 -2.5 -2 -1.5 -1 log ( ǫ ) -8-6-4-202 l og ( G E ) The global errors of v ⊥ with h=16 ǫ Exp-AImp-ATwo P-ABorisslope 2slope 1
Fig. 8.1
The logarithm of the global error against the logarithm of ǫ . implies that ϕ ( b B ) = ϕ (0) I + ϕ ( b ) b B + ϕ ( b ) b B . (8.1)This permits us to compute ϕ ( b B ) v by evaluating the scalars ϕ (0) , ϕ ( b ) , ϕ ( b ),and by forming twice a product of b B with a vector. Note that b Bv = B × v .For the case that ϕ ( ζ ) has only even powers of ζ , we have ϕ ( y ) = 0,and the formula simplifies. Similarly, for the case where only odd powers of ζ are present, we have ϕ (0) = 0 and ϕ ( y ) = 0. For the matrix functions of h/ ǫ -6-4-20 l og ( G E ) The global errors of x for Boris π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of v || for Boris π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of v ⊥ for Boris π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of x for Exp-A π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of v || for Exp-A π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of v ⊥ for Exp-A π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of x for Imp-A π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of v || for Imp-A π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of v ⊥ for Imp-A π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of x for Two P-A π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of v || for Two P-A π π π π h/ ǫ -6-4-20 l og ( G E ) The global errors of v ⊥ for Two P-A π π π π Fig. 8.2
The logarithm of the global error at t = 1 against h/ǫ for ǫ = 1 / and h = 1 /k ,where k = 60 , , . . . , Algorithm 2.1 we thus haveexp( − h b B ) = I − sin( hb ) b b B + 1 − cos( hb ) b b B ,Ψ ( h b B ) = I + 1 − tanc( hb/ b b B ,Φ ( h b B ) = I + 1 − sinc( hb ) − b b B ,Υ ( h b B ) = 1 − sinc( hb ) − hb b B. filtered Boris algorithm for charged-particle dynamics 21 Acknowledgement
C.L. thanks Eric Sonnendr¨ucker for stimulating discussions during the Ober-wolfach workshop 2019-04. This work was supported by the Fonds NationalSuisse, Project No. 200020-159856, by Deutsche Forschungsgemeinschaft, SFB1173, and by the Humboldt Foundation.
References Boris, J. P.
Relativistic plasma simulation-optimization of a hybrid code.
Proceedingof Fourth Conference on Numerical Simulations of Plasmas (November 1970), 3–67.2.
Brizard, A. J., and Hahm, T. S.
Foundations of nonlinear gyrokinetic theory.
Rev.Modern Phys. 79 , 2 (2007), 421–468.3.
Garc´ıa-Archilla, B., Sanz-Serna, J. M., and Skeel, R. D.
Long-time-step methodsfor oscillatory differential equations.
SIAM J. Sci. Comput. 20 (1999), 930–963.4.
Hairer, E., and Lubich, C.
Long-term analysis of the St¨ormer-Verlet method forHamiltonian systems with a solution-dependent high frequency.
Numer. Math. 34 (2016), 119–138.5.
Hairer, E., and Lubich, C.
Energy behaviour of the Boris method for charged-particledynamics.
BIT 58 (2018), 969–979.6.
Hairer, E., and Lubich, C.
Long-term analysis of a variational integrator for charged-particle dynamics in a strong magnetic field.
Numerische Mathematik, to appear (2019).7.
Hairer, E., Lubich, C., and Wanner, G.
Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations , 2nd ed. Springer Series inComputational Mathematics 31. Springer-Verlag, Berlin, 2006.8.
Hochbruck, M., and Lubich, C.
A Gautschi-type method for oscillatory second-orderdifferential equations.
Numer. Math. 83 (1999), 403–426.9.
Kruskal, M.
The gyration of a charged particle.
Rept. PM-S-33 (NYO-7903), Prince-ton University, Project Matterhorn (1958).10.
Northrop, T. G.
The adiabatic motion of charged particles . Interscience Tracts onPhysics and Astronomy, Vol. 21. Interscience Publishers John Wiley & Sons New York-London-Sydney, 1963.11.
Possanner, S.
Gyrokinetics from variational averaging: existence and error bounds.