High-order energy-conserving Line Integral Methods for charged particle dynamics

TL;DR

This paper develops high-order energy-conserving numerical methods based on Line Integral Methods for accurately simulating charged particle dynamics while preserving energy, extending the framework of Hamiltonian Boundary Value Methods.

Contribution

It introduces a new class of high-order energy-conserving methods for charged particles within the LIM framework, expanding the applicability of HBVMs.

Findings

Methods preserve energy accurately over long simulations

Numerical tests confirm theoretical properties

High-order methods outperform lower-order counterparts

Abstract

In this paper we study arbitrarily high-order energy-conserving methods for simulating the dynamics of a charged particle. They are derived and studied within the framework of Line Integral Methods (LIMs), previously used for defining Hamiltonian Boundary Value Methods (HBVMs), a class of energy-conserving Runge-Kutta methods for Hamiltonian problems. A complete analysis of the new methods is provided, which is confirmed by a few numerical tests.

Equations

$$\dot{q}=p,\dot{p}=p\times L(q)-\nabla U(q),q(0)=q_0,p(0)=p_0\in {\mathbb{R}}^3,$$

$$H(q,p)=\frac{1}{2}{p}^{\top }p+U(q),$$

$$\frac{\mathrm{d}}{\mathrm{d}t}H(q,p)=\nabla U(q{)}^{\top }\dot{q}+p^{\top }\dot{p}=p^{\top }\left(p\times L(q)\right)=0.$$

$$u(0)$$

$$u(h)$$

$$H(q_1,p_1)-H(q_0,p_0)$$

$$\dot{q}=p,\dot{p}=B(q)p-\nabla U(q),q(0)=q_0,p(0)=p_0\in {\mathbb{R}}^3,$$

$$B(q)=\left(\begin{array}[]{rrr}0&-\ell_{3}(q)&\ell_{2}(q)\\ \ell_{3}(q)&0&-\ell_{1}(q)\\ -\ell_{2}(q)&\ell_{1}(q)&0\end{array}\right)=-B(q)^{\top}.$$

$$\ddot{q}=B(q)\dot{q}-\nabla U(q),q(0)=q_0,\dot{q}(0)=p_0\in {\mathbb{R}}^m,B(q{)}^{\top }=-B(q),$$

$$\mathrm{deg}{P}_i=i,{\int }_0^1{P}_i(x)P_j(x)\mathrm{d}x={\delta }_{ij},\forall i,j=0,1,\dots ,$$

$$\dot{q}(ch)=\sum _{j\ge 0}{P}_j(c){\gamma }_j(p),\dot{p}(ch)=\sum _{j\ge 0}{P}_j(c)\left[B(q(ch)){\gamma }_j(p)-{\eta }_j(q)\right],c\in [0,1],$$

$${\gamma }_j(p)={\int }_0^1{P}_j(\tau )p(\tau h)\mathrm{d}\tau ,{\eta }_j(q)={\int }_0^1{P}_j(\tau )\nabla U(q(\tau h))\mathrm{d}\tau .$$

$$P_j(c)B(q(ch))=\sum _{i\ge 0}{P}_i(c){\rho }_{ij}(q),j=0,1,\dots ,$$

$${\rho }_{ij}(q)\equiv {\rho }_{ji}(q)={\int }_0^1{P}_i(\tau )P_j(\tau )B(q(\tau h))\mathrm{d}\tau ,i,j=0,1,\dots .$$

$$\dot{q}(ch)=\sum _{i\ge 0}{P}_i(c){\gamma }_i(p),\dot{p}(ch)=\sum _{i\ge 0}{P}_i(c)\left[-{\eta }_i(q)+\sum _{j\ge 0}{\rho }_{ij}(q){\gamma }_j(p)\right],c\in [0,1].$$

$${\int }_0^1{P}_j(c)c^{i}g(ch)\mathrm{d}c=O(h^{j-i}),i=0,\dots ,j.$$

$$c^{i}g(ch)=\sum _{r\ge 0}\frac{g^{(r)}(0)}{r!}{h}^r{c}^{r+i}.$$

$$\int_{0}^{1}P_{j}(c)c^{i}g(ch)\mathrm{d}c=\sum_{r\geq 0}\frac{g^{(r)}(0)}{r!}h^{r}\int_{0}^{1}P_{j}(c)c^{r+i}\mathrm{d}c=\sum_{r\geq j-i}\frac{g^{(r)}(0)}{r!}h^{r}\int_{0}^{1}P_{j}(c)c^{r+i}\mathrm{d}c=O(h^{j-i}).\mbox{~{}$\Box{~{}}$}$$

$${\gamma }_j(\sigma ),{\eta }_j(\sigma )=O(h^j),{\rho }_{ij}(\sigma )=O(h^{|i-j|}).\forall i,j=0,1,\dots .$$

$${\rho }_{ij}(\sigma )=-{\rho }_{ij}(\sigma {)}^{\top },\forall i,j=0,1,\dots .$$

$$\dot{u}(ch)=\sum _{i=0}^{s-1}{P}_i(c){\gamma }_i(v),\dot{v}(ch)=\sum _{i=0}^{s-1}{P}_i(c)\left[-{\eta }_i(u)+\sum _{j=0}^{s-1}{\rho }_{ij}(u){\gamma }_j(v)\right],c\in [0,1],$$

$$u(ch)$$

$$v(ch)$$

$$q_1:=u(h)\equiv q_0+h{\gamma }_0(v),p_1:=v(h)\equiv p_0-h\left[{\eta }_0(u)-\sum _{j=0}^{s-1}{\rho }_{0j}(u){\gamma }_j(v)\right],$$

$$H(q_1,p_1)-H(q_0,p_0)=H(u(h),v(h))-H(u(0),v(0))={\int }_0^h\frac{\mathrm{d}}{\mathrm{d}t}H(u(t),v(t))\mathrm{d}t$$

$$\dot{y}=f(y),\qquad y=\left(\begin{array}[]{c}q\\ p\end{array}\right),$$

$$\frac{\partial }{\partial y_0}y(t,t_0,y_0)=\Phi (t,t_0),\frac{\partial }{\partial t_0}y(t,t_0,y_0)=-\Phi (t,t_0)f(y_0),$$

$$\dot{\Phi }(t,t_0)=f^{\prime }(y(t,t_0,y_0))\Phi (t,t_0),\Phi (t_0,t_0)=I.$$

$$q_1=q(h)+O(h^{2s+1}),p_1=p(h)+O(h^{2s+1}),$$

$$w=\left(\begin{array}[]{c}u\\ v\end{array}\right),\quad w_{1}:=w(h)\equiv\left(\begin{array}[]{c}q_{1}\\ p_{1}\end{array}\right),\quad w_{0}:=w(0)\equiv\left(\begin{array}[]{c}q_{0}\\ p_{0}\end{array}\right),\quad\Phi(t,t_{0})\equiv\left[\Phi_{1}(t,t_{0}),\,\Phi_{2}(t,t_{0})\right],$$

$$v(ch)=\sum _{j=0}^s{P}_j(c){\gamma }_j(v),\nabla U(u(ch))=\sum _{i\ge 0}{P}_i(c){\eta }_i(u),$$

$$B(u(ch))v(ch)=\sum _{i\ge 0}{P}_i(c)\sum _{j=0}^s{\rho }_{ij}(u){\gamma }_j(v),c\in [0,1],$$

Full text

**Keywords: ** Charged particle dynamics, Lorentz force system, plasma physics, Energy-conserving methods, Line Integral Methods, LIMs, Hamiltonian Boundary Value Methods, HBVMs.

**MSC: ** 65L05, 65P10.

1 Introduction

We shall here be concerned with the dynamics of a charged particle, which is described by the following system of ODEs,

[TABLE]

where $-\nabla U(q)$ and $L(q)$ are the electric and magnetic fields, respectively.333The equations (1) are sometimes referred to as Lorentz force system (see, e.g., [21, 29]). Hereafter, we shall assume both of them to be time independent and suitably smooth functions of $q$. For this motion, the energy,

[TABLE]

turns out to be conserved, along the solution of (1) since, by the chain rule,

[TABLE]

Problem (1) is a relevant one in plasma physics, since many important phenomena in plasmas can be understood and analyzed in terms of the motion of a single-particle [2]. This is due to the fact that collisions occur infrequently in hot plasmas and, in fact, (1) provides a collisionless model of the plasma. However, despite the simplicity of the model, generic numerical methods are not appropriate for its long-time simulation, since an (unphysical) numerical drift in the energy (2) is experienced, giving rise to a complete wrong solution orbit (see, e.g., [30]).

Though several approaches have been recently proposed for numerically solving (1) (see, e.g., [23, 24, 25, 31, 32, 33]), a widely used method is the Boris method [3], which is symmetric and second-order accurate. In fact, because of its ease of implementation, such a method has become a de facto standard for simulating (1). The analysis of its good behaviour, in turn, has been the subject of many investigations (see, e.g., [20, 22, 30]). The Boris method, however, is not energy-conserving and, in fact, a numerical drift may still be observed [22].444In particular, the talk given by Ernst Hairer at the 2019 RSME Congress, presenting the results in [22], has inspired the investigations reported in this paper. On the other hand, a second-order energy conserving method based on a line-integral approach has been recently proposed in [29]. Nevertheless, to the best of our knowledge no arbitrarily high-order energy-conserving methods are known for the simulation of problem (1), and this motivates the present paper. The new energy-conserving methods will be derived within the framework of Line Integral Methods [26, 27, 28], already used to devise a number of energy-conserving methods for various conservative problems [1, 4, 8, 15, 16, 6, 7], with the main instance provided by Hamiltonian Boundary Value Methods (HBVMs), for the numerical solution of Hamiltonian problems [11, 12, 14, 17]. We also refer to the monograph [9] and to the review paper [10], for an overview.

The basic idea Line Integral Methods rely on is that of rewriting the conservation property (3) at $t=h$ in integral form, and defining approximate paths $u(t)\approx q(t)$, and $v(t)\approx p(t)$ such that:555This procedure defines, indeed, the very first step of application of a one-step method.

[TABLE]

For our purposes, it will be convenient to rewrite (1) in the equivalent form

[TABLE]

where, setting L(q)=\left(\begin{array}[]{ccc}\ell_{1}(q),&\ell_{2}(q),&\ell_{3}(q)\end{array}\right)^{\top},

[TABLE]

In so doing, the arguments that we shall use for solving (7)-(8) can be naturally extended for solving more general problems, in the form

[TABLE]

which possess the same invariant (2), with $\dot{q}=p$.

With this premise, the structure of the paper is as follows: in Section 2 we describe the framework in which the methods will be derived; in Section 3 we provide a fully discrete method; in Section 4 its actual implementation is studied; in Section 5 we present a few numerical tests confirming the theoretical findings; at last, in Section 6 we give some concluding remarks.

2 Derivation of the method

Following the approach in [14], let us now rewrite the problem (7)-(8) by expanding the right-hand sides along a suitable orthonormal basis, which we choose as the orthonormal Legendre polynomial basis $\{P_{j}\}$ on the interval $[0,1]$,

[TABLE]

with $\delta_{ij}$ the Kronecker symbol. We then obtain, at first:

[TABLE]

with

[TABLE]

We observe that, in (11), $B(q(ch))$ has not yet been expanded. Next, consider the expansions

[TABLE]

with

[TABLE]

As a result, from (11) we eventually arrive at:

[TABLE]

The following properties hold true.666Lemma 1 is a generalization of [14, Lemma 1].

Lemma 1

Assume $g:[0,h]\rightarrow V$, with $V$ a vector space, admit a Taylor expansion at 0. Then, for all $j=0,1,\ldots:$

[TABLE]

Proof By the hypotheses on $g$, one has:

[TABLE]

Consequently, for all $i=0,\dots,j$, by virtue of (10) it follows that:

[TABLE]

Corollary 1

With reference to (12) and (13), for any suitably regular path $\sigma:[0,h]\rightarrow\mathbb{R}^{3}$ one has:777 The path would become $\sigma:[0,h]\rightarrow\mathbb{R}^{m}$, in the case of the more general problem (9).

[TABLE]

Lemma 2

With reference to (13), for any path $\sigma:[0,h]\rightarrow\mathbb{R}^{3}$ one has:888 Footnote 7 applies also here.

[TABLE]

Proof The statement follows from the definition (13) and the skew-symmetry of matrix $B$ in (8).  $\Box{~{}}$

Next, in order to obtain polynomial approximations $u\approx q$ and $v\approx p$ of degree $s$, we truncate the infinite series in (14) after $s$ terms, thus getting:

[TABLE]

with $\gamma_{j}(v),\eta_{j}(u),\rho_{ij}(u)$ defined according to (12) and (13), by formally replacing $q,p$ with $u,v$, respectively. After that, we need to satisfy the Line Integral conditions (4)–(6), i.e.,

- requirement (4):

we impose the initial conditions $u(0)=q_{0}$, $v(0)=p_{0}$, thus obtaining, integrating both equations in (17),

[TABLE]

- requirement (5):

accordingly, we set

[TABLE]

where, by virtue of (10), we took into account that $\int_{0}^{1}P_{i}(c)\mathrm{d}c=\delta_{i0}$ ;

- requirement (6):

at last, next theorem states the property of energy-conservation.

Theorem 1

With reference to (2) and (17)–(19), one has $H(q_{1},p_{1})=H(q_{0},p_{0})$.

Proof In fact, one has:

[TABLE]

where the last equality follows from (16).  $\Box{~{}}$

We now study to what extent $q_{1},p_{1}$ approximate $q(h),p(h)$, respectively. For this purpose, in order to simplify the notation, we rewrite (7)-(8) as

[TABLE]

and denote by $y(t,t_{0},y_{0})$ the solution of such an equation, satisfying the initial condition $y(t_{0})=y_{0}$. We also recall the following known perturbation results:

[TABLE]

where $\Phi(t,t_{0})$ is the fundamental matrix solution of the associated variational problem,

[TABLE]

The following result then holds true, whose proof is based on the arguments used in [14, Theorem 1].

Theorem 2

With reference to (1) and (17)–(19), one has

[TABLE]

i.e., the approximation procedure has order $2s$.

Proof In fact, by using the notation (20), the perturbation results (21), setting

[TABLE]

with $\Phi_{i}(t,t_{0})\in\mathbb{R}^{6\times 3}$, $i=1,2$,999 $\Phi_{i}(t,t_{0})\in\mathbb{R}^{2m\times m}$ in the case of the more general problem (9). and considering that from (11)–(14) (see also (17) and (18)) it follows that

[TABLE]

one has:

[TABLE]

where we have set

[TABLE]

We observe that, by virtue of Lemma 1 and Corollary 1,

[TABLE]

Consequently, one has:

[TABLE]

and, as a result, one concludes that $(*)=O(h^{2s+1})$.  $\Box{~{}}$

3 Discretization

The approximation procedure (17)–(19) described in the previous section is not yet a ready to use numerical method. In fact, in order for this to happen, the integrals $\gamma_{i}(v),\eta_{i}(u),\rho_{ij}(u)$, $i,j=0,\dots,s-1$, defined in (12)-(13) need to be conveniently computed or approximated. We observe that, since $v\in\Pi_{s}$, then $\gamma_{i}(v)$ can be exactly computed by using the interpolatory quadrature formula of order $2s$ based at the zeros of $P_{s}$. If we denote $(\hat{c}_{\ell},\hat{b}_{\ell})$ the nodes and weights of such a quadrature,101010I.e., $P_{s}(\hat{c}_{\ell})=0$, $\ell=1,\dots,s$. one has then:

[TABLE]

We shall use the same quadrature for approximating $\rho_{ij}(u)$:

[TABLE]

with $\hat{\Delta}_{ij}(h)$ the quadrature error. Finally, as it has been done in the case of HBVMs [12], for approximating $\eta_{i}(v)$ we shall use a Gauss-Legendre quadrature of order $2k$, i.e., the interpolatory quadrature rule based at the zeros of $P_{k}$, for a convenient value $k\geq s$, with abscissae and weights $(c_{\ell},b_{\ell})$:111111I.e., $P_{k}(c_{\ell})=0$, $\ell=1,\dots,k$.

[TABLE]

with $\Delta_{i}(h)$ the quadrature error. Concerning the quadrature errors, one verifies that, for all $i,j=0,\dots,s-1:$

[TABLE]

Consequently, it is straightforward to prove that the results of Corollary 1 and Lemma 2 continue formally to hold for $\hat{\gamma}_{i}(v),\hat{\eta}_{i}(u),\hat{\rho}_{ij}(u)$ defined in (22)–(24) (and, of course, $\hat{\rho}_{ij}(u)=\hat{\rho}_{ji}(u)$). In so doing, the polynomials (17)–(18) respectively become 121212For sake of brevity, we shall continue to denote such polynomials by $u$ and $v$, respectively.

[TABLE]

and

[TABLE]

with the new approximations given by

[TABLE]

in place of (19). Clearly, the Line Integral requirements (4) and (5) are satisfied by the new polynomial paths $u$ and $v$ defined by (22)–(28). Concerning the requirement (6), the following result holds true.

Theorem 3

With reference to (2) and (22)–(28), one has either $H(q_{1},p_{1})=H(q_{0},p_{0})$, if  $U\in\Pi_{\nu}$ with $\nu\leq 2k/s$, or $H(q_{1},p_{1})=H(q_{0},p_{0})+O(h^{2k+1})$, otherwise.

Proof In fact, one has, by considering that  $\gamma_{i}(v)=\hat{\gamma}_{i}(v)$  and  $\hat{\rho}_{ij}(u)=\hat{\rho}_{ji}(u)=-\hat{\rho}_{ij}(u)^{\top}$:

[TABLE]

[TABLE]

The statement then follows from (25) and considering that $\gamma_{i}(u)=O(h^{i})$.  $\Box{~{}}$

Remark 1

As is clear from Theorem 3, an exact energy conservation is obtained in the polynomial case, by choosing $k$ large enough. However, also in the non-polynomial case, one can always gain a practical energy conservation, by choosing $k$ large enough so that the $O(h^{2k+1})$ energy error falls within the round-off error level.

Next results states that the order $2s$ of the approximation procedure (17)–(19) is retained by the new one.

Theorem 4

With reference to (1) and (22)–(28), for all $k\geq s$ one has

[TABLE]

Proof In fact, by using the same notations and preliminary results used in the proof of Theorem 2, one has:

[TABLE]

where the last equality follows from the proof of Theorem 2, in which formally the same terms were involved, except those including the quadrature errors (25). Concerning these latter terms, by recalling that

[TABLE]

one has:

[TABLE]

Consequently, the statement follows.  $\Box{~{}}$

Definition 1

Hereafter, we shall refer to the method defined by (22)–(28) as Line Integral Method with parameters $(k,s)$, in short LIM$(k,s)$, for solving problem (7)–(8).

We then conclude that, according to Theorem 4, the LIM$(k,s)$ has order $2s$, for all allowed values of $s$ and $k\geq s$. Moreover, according to Theorem 3 and Remark 1, it is energy-conserving, either exactly or practically, by choosing $k$ suitably large. This, in turn, will be not a drawback, since the discrete problem that one has to solve has dimension $s$, independently of $k$, as we are going to see in the next section.

4 The discrete problem

We now study the efficient implementation of the method (22)–(28). For this purpose, we need the following matrices, defined by the Legendre polynomial basis (10) and the nodes and weights of the Gauss-Legendre quadratures of order $2s$, $(\hat{c}_{\ell},\hat{b}_{\ell})$, and $2k$, $(c_{\ell},b_{\ell})$:

[TABLE]

Moreover, we need to introduce the vector  {\bm{1}}_{r}=\left(\begin{array}[]{ccc}1,&\dots,1\end{array}\right)^{\top}\in\mathbb{R}^{r}, and the block vectors and matrices, with reference to (22)–(24):

[TABLE]

We observe that, for computing what in (51), we need to evaluate $v(\hat{c}_{\ell}h),u(\hat{c}_{\ell}h)$, $\ell=1,\dots,s$, and $u(c_{\ell}h)$, $\ell=1,\dots,k$. For this purpose, to have a more compact notation, we define the vectors:

[TABLE]

A similar notation will be used for the functions in (1)–(8), when a block vector argument is specified. With reference to (43) and (51), one then obtains

[TABLE]

and

[TABLE]

Hereafter, $I$ will denote the identity matrix having the size of the vectors $q_{0}$ and $p_{0}$ (i.e., 3 for problem (1), or $m$ for problem (9)). The last equation can be rewritten, by defining the block vector

[TABLE]

as:

[TABLE]

From (28), one then obtains that the new approximations are given by

[TABLE]

Consequently, we need to compute the three block vectors $\hat{\bm{\gamma}}(v)$, $\hat{\bm{\eta}}(u)$, and ${\bm{\varphi}}(u,v)$, respectively defined in (51) and (53). By using the matrices defined in (43), and the formulae (22)–(24), we then obtain

[TABLE]

and

[TABLE]

having set

[TABLE]

Nevertheless, by taking into account that $\hat{\cal P}_{s}^{\top}\hat{\Omega}=\hat{\cal P}_{s}^{-1}$ and the first equation in (56), one obtains that (57) can be rewritten as

[TABLE]

Finally, by considering the expressions (52) and (54), we eventually obtain the following set of equations:

[TABLE]

where we have removed the arguments of the (block) vectors, since now this is only an algebraic system of equations which, moreover, can be further simplified. In fact, we observe, at first, by taking into account

[TABLE]

and

[TABLE]

that (59) can be rewritten as

[TABLE]

Next, by plugging the right-hand side of (64) in those of (60)-(61), and taking into account that

[TABLE]

one obtains:

[TABLE]

Further, by defining the block vector

[TABLE]

subtracting (66) from (4) provides us with the (block) vector equation

[TABLE]

Once (69) has been solved, one easily computes the first block entry of (64) by taking into account (63), so that the new approximations (55) become

[TABLE]

Remark 2

As is clear from (70), in order to obtain the new approximation $q_{1}$, one must have $s\geq 2$, so that all LIM$(k,s)$ methods, with $k\geq s$, have order $2s\geq 4$. Consequently, the methods here derived are completely different from the second-order method studied in [29].

Moreover, in the case where the magnetic field is zero (i.e., ${\cal B}=O$ in (69)), then the discrete problem (69)-(70) reduces to that generated by a HBVM$(k,s)$ method applied to the special second order problem  $\ddot{q}+\nabla U(q)=0$,  $q(0)=q_{0}$,  $\dot{q}(0)=p_{0}$  (see, e.g., [9, Chapter 4.1.2]).

Also, we observe that, introducing the notation

[TABLE]

then (69) can be rewritten in a form closer to that of the original problem (1), i.e.,

[TABLE]

Last, but not least, we stress that, in order to obtain an energy-conserving method of order $2s$, we need to solve, at each integration step, the set of $3s$ algebraic equations (69) (i.e., (71)), whose dimension is remarkably independent of $k$, as previously anticipated.

Next result guarantees the existence and uniqueness of the solution of (69).

Theorem 5

Let $B(q)p$ and $\nabla U(q)$ in (7) be continuous and satisfy a Lipschitz condition with constant $\mu$ (with respect to $p$ and $q$) given by

[TABLE]

Then, if the stepsize $h$ is small enough to satisfy

[TABLE]

there exists a unique solution of (69), and the fixed-point iteration

[TABLE]

converges to it.

Proof First, note that the Lipschitz conditions ensure that the same functions with block vector arguments are also Lipchitz with the same constant. Then, it is straightforward to prove that the function defined by the right-hand side of (69) is a contraction with respect to ${\bm{\psi}}$ if $h$ satisfies (72). Consequently, the Fixed-Point Theorem ensures the existence and uniqueness of the solution, as well as the convergence of the iteration (73).  $\Box{~{}}$

4.1 Symmetry

The compact structure of the discrete problem (69)-(70) allows us to prove the important property of symmetry of the method, under the following symmetry condition of the abscissae,

[TABLE]

which is clearly satisfied by the choice of the Gauss-Legendre abscissae $P_{s}(\hat{c}_{\ell})=0$, $\ell=1,\dots,s$, $P_{k}(c_{\ell})=0$, $\ell=1,\dots,k$. For interpolatory quadrature formulae, in turn, (74) implies the symmetry of the weights:

[TABLE]

We also need to define the following matrices,

[TABLE]

and related preliminary results.

Lemma 3

With reference to the matrices defined in (43), (63), and (76), and the vector ${\bm{e}}_{1}$ defined in (62), one has:

[TABLE]

Proof The properties on the first line follows by the fact that $\hat{D},\hat{P},P$ are symmetric and orthogonal, and from (63). The properties on the subsequent two lines derive from the following symmetries of the Legendre polynomials:

[TABLE]

At last, the properties on the last line follow from (75).  $\Box{~{}}$

We recall that the method (69)-(70) is symmetric if, when starting from $(q_{1},p_{1})$ with timestep $-h$, it brings back to $(q_{0},p_{0})$.131313The proof of symmetry will be akin to that done for HBVMs in [9, Theorem 3.11].

Theorem 6

Under the hypotheses (74), the method (69)-(70) is symmetric.

Proof By using the method (69)-(70) for solving (7)-(8) starting from $(q_{1},p_{1})$ with stepsize $-h$, we obtain the equations

[TABLE]

We have then to prove that

[TABLE]

Preliminarily, we observe that the symmetry conditions (74) can be respectively rewritten as

[TABLE]

Moreover, from the last two equations in (81) and (70), one has

[TABLE]

Then, (82) follows if we show that (see (68))

[TABLE]

We shall prove this statement by verifying that the (block) vector ${\bm{\psi}}^{*}$ satisfies the very same equation (81) which defines $\bar{\bm{\psi}}$. From (69) and Lemma 3, one has:

[TABLE]

Next, upon observing that (see (63), (70), and (84))

[TABLE]

we obtain, taking into account again that $p_{0}=p_{1}-h\psi_{0}^{*}$,

[TABLE]

The statement then follows by considering that (see (43), (63), (65), and Lemma 3)

[TABLE]

and, similarly,  \left[{\bm{c}}{\bm{e}}_{1}^{\top}-{\cal I}_{s}X_{s}^{\top}\right]={\cal I}_{s}X_{s}.\,\mbox{~{}\Box{~{}}}

4.2 Solving the discrete problem

When the stepsize $h$ satisfies the conditions of Theorem 5, we proved that the fixed-point iteration (73) converges to the solution of (69). However, sometimes a simplified Newton iteration may be more appropriate for solving the equation

[TABLE]

Neglecting the $O(h^{2})$ terms in the simplified Jacobian, and taking into account (63), one then obtains the iteration

[TABLE]

This iteration, though straightforward, requires the factorization of a matrix having dimension $3s$. This complexity can be reduced, e.g., by using variants akin to those used for HBVMs [13, 5]. In particular, we consider the blended iteration associated to (86), which, by setting 141414We refer to [18] for a more comprehensive analysis of blended methods.

[TABLE]

reads

[TABLE]

In so doing, only the $3\times 3$ matrix $\Theta$ in (87) needs to be computed.

Sometimes, the electric field may be much stronger than the magnetic one. In such a case, the simplified Newton iteration (86) becomes, by neglecting the contribution of the magnetic field, and considering that ${\cal P}_{s}^{\top}\Omega{\cal I}_{s}=X_{s}$, as defined in (63),

[TABLE]

We can again use a corresponding blended iteration [9, 19],

[TABLE]

where $\rho_{s}$ is the same parameter defined in (87), and

[TABLE]

Consequently, also in this case, only a $3\times 3$ matrix need to be factored.

Remark 3

We observe that, in the case of the more general problem (9), the previous arguments remains formally the same, with the only difference that, now, the matrices $\Theta$ and $\Theta_{1}$ defined in (87) and (91), respectively, have dimension $m\times m$.

Remark 4** (The case of a constant $B$)**

When the magnetic field $L(q)$ is uniform, then with reference to matrix $B(q)$ defined in (8), one has:

[TABLE]

In this case, taking into account (62)-(63), one has that the discrete problem (69) simplifies to

[TABLE]

Moreover, also the blended iteration (87)-(88) greatly simplifies, since the matrix $\Theta$ turns out to be given by  $\Theta=(I-h\rho_{s}B)^{-1}$  and, therefore, it is constant for all time-steps.

5 Numerical tests

We here report a few numerical tests, aimed at assessing the theoretical findings of the previous sections. We also compare the LIM$(k,s)$ method (69)-(70) with the Boris method.151515An efficient implementation of this latter method is explained in [22]. All numerical tests have been done on a 2.8GHz Intel core i7 computer with 16GB of memory, running Matlab 2017b.

Example 1.

We start considering the problem (1)–(8) with:161616Hereafter, $q_{i}$ will denote the $i$th entry of $q$.

[TABLE]

for which in [22] it has been shown that the Boris algorithm exhibits a $O(th^{2})$ drift in the energy. For this problem, the LIM$(2s,s)$ method turns out to be energy-conserving, according to the result of Theorem 3, since $U(q)$ is a polynomial of degree 4. The corresponding order is $2s$, as stated in Theorem 4. The drift for the Boris method is confirmed by the left-plot in Figure 2, where the Hamiltonian error for such method, using a stepsize $h=10^{-2}$ over the interval $[0,3\cdot 10^{4}]$, is shown: the obtained plot, upon scaling by $h^{2}$ and reversing the sign, perfectly fits those in [22, Fig. 1], thus confirming the $O(th^{2})$ energy drift. In the right-plot of the same figure, there is the Hamiltonian error for the LIM(4,2) method using the same stepsize: as one may see, in such a case energy-conservation is gained (clearly, up to round-off errors). It must be observed, however, that the execution time for the Boris method is about 11 times smaller than that for the LIM(4,2) method (101 sec vs. 1149 sec). Nevertheless, the availability of arbitrarily high-order energy-conserving line-integral methods (i.e., LIM$(2s,s)$, $s=2,3,\dots$), coupled with the efficient nonlinear iteration (87)-(88) studied in Section 4.2, allows to recover on the efficiency of the methods, when we compare the solution error vs. the execution time. For this, we fix a shorter interval, say $[0,100]$, and we compare various energy-conserving line integral methods, along with the Boris method, in terms of solution error vs. execution time. The obtained results are plotted in Figure 2, from which one realizes that the higher order LIMs are more effective than the lower order ones and, in any case, all of them perform better than the Boris method, which is the lowest order method, among those considered.

Example 2.

Next, we consider the problem defined by (93)-(96), and

[TABLE]

for which the energy error of the Boris algorithm defines a random walk [22]. For this problem, since $U(q)\in\Pi_{4}$ and $L(q)\in\Pi_{1}$, one has that the quadratures (22)–(24) are exact for all LIM$(2s,s)$ methods, so that the discrete method (27)-(28) coincides with the approximation procedure (18)-(19). In particular, all of such methods are energy-conserving and have order $2s$. In Table 1 we list the obtained results when solving the problem on the interval $[0,25]$ with stepsize $h=0.05/n$: as one may see, the expected convergence order is confirmed.

Example 3.

Finally, we consider the 2D dynamics of a charged particle in a static, non-uniform electromagnetic field. The model is an important application in the study of the single particle motion and the guiding center dynamics [29, 30]. The problem that we consider is defined by:

[TABLE]

and $L(q)$ defined as in (98). In this case, the motion, depicted in the upper-left plot of Figure 3 for $t\in[0,10^{3}\pi]$, occurs in the $(q_{1},q_{2})$-plane, and another invariant is given by the angular momentum [23],

[TABLE]

We solve numerically this problem by using the Boris and the LIM$(2s,s)$ methods, $s=2,3,4,5$, with stepsize $h=\pi/10$. The numerical solution of the Boris method is depicted in the upper-right plot of Figure 3, whereas in the lower plots one finds the numerical solutions computed by the LIM(4,2) method (left) and LIM(6,3) method (right): as is clear, the higher the order of the method, the better the numerical solution. For completeness, in Table 2 we also list the maximum solution, Hamiltonian, and angular momentum errors for all methods, along with the measured execution times, again confirming that the higher order methods are more effective. In particular, it is worth mentioning that all considered LIMs are practically energy-conserving and, remarkably, their execution times are almost the same.

6 Conclusions

In this paper we have derived symmetric and arbitrarily high-order energy-conserving methods for the numerical simulation of the dynamics of a charged particle, within the framework of Line Integral Methods. A complete analysis of the methods, including their efficient implementation, has been given. Possible extensions to more general problems have also been sketched. Numerical results on significant test problems, including the guiding center dynamics, duly confirm the theoretical findings.

Acknowledgements

This work has been finalized during the visit of the first author at the I.U.M.A., Universidad de Zaragoza, whose financial support is acknowledged (MINECO project MTM2016-77735-C3-1-P). The first author also acknowledges the interesting discussions with Ernst Hairer, during the RSME meeting in Santander, concerning the Boris method.