On a finite-difference scheme defining a birational non-quadratic map between time layers

Full Text

1.     Introduction

Dynamical systems with quadratic right–hand sides can be approximated by finite-difference schemes. Timestep of these schemes is a birational map [1–5]. Appelroth’s quadratization allows any dynamical system with a polynomial right-hand side to be transformed into a system with a quadratic right-hand side to which the Kahan scheme can be applied [6]. Consequently, a very broad class of continuous dynamical systems admits approximation by 𝑡-symmetric finite-difference schemes that define a birational map between time layers. Following [5], we refer to such schemes as reversible. This property is of fundamental physical relevance. The time evolution of a dynamical system should, in general, be reversible: every initial state corresponds to exactly one final state, and every final state comes from exactly one initial state. Surprisingly, although the discrete-time methods often preserve one-to-one correspondence, it is not always true for the original continuous systems over the different time intervals [5].

It is worth mentioning that in the papers above the time step was implemented using a quadratic Cremona transformation. However, in the 1980s, V. N. Abrashin proposed a multicomponent alternating-direction method, which was successfully applied by S. N. Sytova. This method was used for modeling the nonlinear dynamics of electromagnetic radiation generated by charged particle beams in multidimensional spatially periodic structures [7, n. 4]. Like Appelroth’s quadratization, this approach also introduces extra variables. However, instead of using a complicated substitution, it does so in a very simple and clear way: by duplicating every original variable. Because of this straightforward construction, we call it the system doubling method. A closer look shows that this method can be used almost for any dynamical system to build a finite-difference scheme where each time step is described by a Cremona transformation. In this paper, we give a general description of the doubling method and explain its key algebraic properties.

2.     Methods

2.1.      System doubling method

Let 𝔵 = (𝑥₁,…,𝑥_𝑛) and there is a dynamical system

𝑑𝔵

                                                                            = 𝔣(𝔵).                                                                                                (1)

𝑑𝑡

For simplicity, assume that the right-hand sides in (1) are polynomials.

Introduce a new set of variables 𝔵^′, which duplicate the original set of variables 𝔵, and form the doubled system:

                                                                            𝑑𝔵             _′

                                                                          ⎧ ^𝑑𝑡_′ = 𝔣(𝔵 ),                                                                                              (2)

⎨𝑑𝔵

⎩ 𝑑𝑡 = 𝔣(𝔵).

We formulate the equivalence of these systems as a theorem.

Theorem 9. Every solution 𝔵 = 𝔥(𝑡) of system (1) extends to a solution 𝔵 = 𝔵^′ = 𝔥(𝑡) of the doubled system (2). Conversely, any solution of system (2) that satisfies the initial condition 𝔵 = 𝔵^′ at some time 𝑡₀ satisfies the equation 𝔵(𝑡) = 𝔵^′(𝑡) for all 𝑡, and 𝔵 = 𝔵(𝑡) is a solution of the original system (1).

Proof. (i) Suppose 𝔵 = 𝔥(𝑡) is a solution of system (1). Then

𝑑𝔥

= 𝔣(𝔥).

𝑑𝑡

Substituting 𝔵(𝑡) = 𝔵^′(𝑡) = 𝔥(𝑡) into system (2), both equations are satisfied identically, since

                                             𝑑𝔵 = 𝑑𝔥 = 𝔣(𝔥) = 𝔣(𝔵′),                     𝑑𝔵′ = 𝑑𝔥 = 𝔣(𝔥) = 𝔣(𝔵).

                                          𝑑𝑡        𝑑𝑡                                            𝑑𝑡         𝑑𝑡

Thus, (𝔥(𝑡),𝔥(𝑡)) is a solution of the doubled system (2).

(ii) Conversely, let (𝔵(𝑡),𝔵^′(𝑡)) be a solution of (2) satisfying the initial condition 𝔵(0) = 𝔵^′(0). Define the differences

                                                        𝑢_𝑖(𝑡) = 𝑥_𝑖^′(𝑡) − 𝑥_𝑖(𝑡),                𝑖 = 1,…,𝑛,

and collect them into the vector 𝔲(𝑡) = 𝔵^′(𝑡) − 𝔵(𝑡). Differentiating gives

                                                        𝑑𝑢𝑖        𝑑𝑥𝑖′       𝑑𝑥𝑖                                 ′).

                                                               =            −              = 𝑓𝑖(𝔵) − 𝑓𝑖(𝔵

                                                         𝑑𝑡          𝑑𝑡         𝑑𝑡

Since each component 𝑓_𝑖 is a polynomial (by assumption), the difference 𝑓_𝑖(𝔵) − 𝑓_𝑖(𝔵^′) can be expressed as a linear combination of the differences 𝑥_𝑗^′ − 𝑥_𝑗 = 𝑢_𝑗 with polynomial coefficients. More precisely, there exist polynomials 𝑔_𝑖𝑗(𝔵,𝔵^′) such that

𝑛

𝑓_𝑖(𝔵) − 𝑓_𝑖.

Hence, 𝔲(𝑡) satisfies the linear homogeneous system

                                                 ^{𝑑𝑢𝑖               𝑛                                    ′}(𝑡))𝑢_𝑗,               𝑖 = 1,…,𝑛.

= ∑ 𝑔_𝑖𝑗(𝔵(𝑡),𝔵

𝑑𝑡

𝑗=1

This system has continuous (in fact, smooth) coefficients along the trajectory (𝔵(𝑡),𝔵^′(𝑡)), and the initial condition 𝔲(0) = 0 holds by hypothesis. By uniqueness of solutions to linear ODEs, it follows that 𝔲(𝑡) ≡ 0 for all 𝑡 in the interval of existence. Therefore,

                                                                   𝔵(𝑡) = 𝔵^′(𝑡)       for all 𝑡.

Substituting this identity into the first equation of (2) results in

𝑑𝔵

= 𝔣(𝔵),

𝑑𝑡

so 𝔵(𝑡) is indeed a solution of the original system (1).                                                                                                                  

Remark 1. The polynomials 𝑢_𝑖 arising in the proof are a natural generalization of Darboux polynomials [8–10].

System (2) is more complex than system (1), yet it admits a straightforward finite-difference discretization:

                                                                        𝔵 − 𝔵 = 𝔣(𝔵̂   ^′)𝛥𝑡,

                                                                        ^{𝔵^′̂ − 𝔵^′ = 𝔣(𝔵)𝛥𝑡.̂                                                                                           (3)

This scheme defines a birational map between the points (𝔵,𝔵^′) and (𝔵,̂ 𝔵^′̂ ).

In contrast to Kahan’s method, the scheme (3) is not 𝑡-symmetric [5]. However, it possesses a generalized analogue of 𝑡-symmetry: the equations (3) are invariant under the transformation

                                                        𝛥𝑡 ↦ −𝛥𝑡,            𝔵 ↦ 𝔵^′̂ ,             𝔵^′ ↦ 𝔵.̂

We shall refer to this property as generalized 𝑡-symmetry.

We have implemented the doubling-based scheme (3) in our fdm package for SAGE. All computations presented below were made using this package.

2.2.      Quadratic integrals

Theorem 10. Suppose that system (1) admits a homogeneous quadratic first integral

𝑣 = ∑𝑣𝑖𝑗 𝑥𝑖𝑥𝑗.

𝑖,𝑗

Then the doubled system (2) possesses the first integral

𝑤 = ∑𝑣𝑖𝑗 𝑥𝑖𝑥𝑗′,

𝑖,𝑗

which is exactly preserved by the discrete scheme (3), i.e.,

𝑤(𝔵,̂ 𝔵^′̂ ) = 𝑤(𝔵,𝔵^′).

Proof. We first observe that

𝑣 = 𝔵^T𝑉𝔵,

where 𝑉 is the symmetric matrix representing the quadratic form 𝑣. By assumption,

𝑑𝑣

= 0,

𝑑𝑡

along solutions of (1), i.e.,

𝔣(𝔵)^T𝑉𝔵 + 𝔵^T𝑉𝔣(𝔵) = 0.

Since 𝑉 is symmetric, this simplifies to

                                                                           𝔵^T𝑉𝔣(𝔵) = 0.                                                                                               (4)

This identity holds for all 𝔵 ∈ ℝ_𝑛, as it is satisfied along every solution of (1).

The quantity 𝑤 is the bilinear form associated with the quadratic form 𝑣, so we can write

𝑤 = 𝔵^T𝑉𝔵^′.

Differentiating 𝑤 along solutions of the doubled system (2) gives

𝑑𝑤 = 𝔣(𝔵′)T𝑉𝔵′ + 𝔵T𝑉𝔣(𝔵).

𝑑𝑡

Each term on the right-hand side vanishes by identity (4). Consequently, 𝑤 is a first integral of system (2).

We now turn to the finite-difference setting. Consider the increment

𝑤̂ − 𝑤 = 𝑤(𝔵,̂ 𝔵^′̂ ) − 𝑤(𝔵,𝔵^′) = 𝔵^T̂ 𝑉𝔵^′̂ − 𝔵^T𝑉𝔵^′. Adding and subtracting the term 𝔵^T̂ 𝑉𝔵^′ gives

𝑤̂ − 𝑤 = (𝔵^T̂ 𝑉𝔵^′̂ − 𝔵^T̂ 𝑉𝔵^′) + (𝔵^T̂ 𝑉𝔵^′ − 𝔵^T𝑉𝔵^′).

Using the discrete scheme (3), we have

                                                        𝔵^′̂ − 𝔵^′ = 𝔣(𝔵)𝛥𝑡,̂                𝔵 − 𝔵 = 𝔣(𝔵̂  ^′)𝛥𝑡,

and therefore

𝔵^T̂ 𝑉𝔵^′̂ − 𝔵^T̂ 𝑉𝔵^′ = 𝔵^T̂ 𝑉(𝔵^′̂ − 𝔵^′) = 𝛥𝑡 𝔵^T̂ 𝑉𝔣(𝔵),̂

                                                    𝔵^T̂ 𝑉𝔵^′ − 𝔵^T𝑉𝔵^′ = (𝔵 − 𝔵̂       )^T𝑉𝔵^′ = 𝛥𝑡 𝔣(𝔵^′)^T𝑉𝔵^′.

By identity (4), both terms vanish. Hence 𝑤̂ = 𝑤, which shows that 𝑤 is exactly preserved by the discrete map (3).                  

2.3.      Hamiltonian formulation

Suppose the original system (1) is Hamiltonian, i.e., it can be written as

𝑑𝑝𝑖                  𝜕𝐻         𝑑𝑞𝑖             𝜕𝐻                                                            = −         ,                =          ,         𝑖 = 1,…,𝑛.                                                    𝑑𝑡            𝜕𝑞𝑖               𝑑𝑡         𝜕𝑝𝑖 Then the doubled system takes the form                                                    𝑑𝑝𝑖                   𝜕𝐻′             𝑑𝑞𝑖             𝜕𝐻′                                                 ⎧⎪ 𝑑𝑡 = − 𝜕𝑞𝑖′ ,           𝑑𝑡 = 𝜕𝑝𝑖′ ,           𝑖 = 1,…,𝑛, ⎨⎪𝑑𝑝𝑖′ = −𝜕𝐻 ,     𝑑𝑞′𝑖 = 𝜕𝐻 ,                  𝑖 = 1,…,𝑛. ⎩ 𝑑𝑡      𝜕𝑞𝑖                    𝑑𝑡               𝜕𝑝𝑖 where 𝐻 = 𝐻(𝑝,𝑞) and 𝐻′ = 𝐻(𝑝′,𝑞′). We introduce
𝐻2 = 𝐻(𝑝,𝑞) + 𝐻(𝑝′,𝑞′) and rewrite the system in the following way	(5)
𝑑𝑝             𝜕𝐻                ′                                                 ⎧         𝑖 = −        ′𝑖2 ,  𝑑𝑞𝑑𝑡𝑖 = 𝜕𝐻𝜕𝑝𝑖2 ,    𝑖 = 1,…𝑛,                                                 ⎪ 𝑑𝑡              𝜕𝑞                                                 ⎨𝑑𝑝𝑖′            𝜕𝐻2        𝑑𝑞𝑖        𝜕𝐻2                                                 ⎪           = −           ,                =          ′ ,         𝑖 = 1,…𝑛.                                                 ⎩ 𝑑𝑡              𝜕𝑞𝑖                𝑑𝑡          𝜕𝑝𝑖	(6)

This combined Hamiltonian allows us to write the doubled system in a compact Hamiltonian form. It is now clear that this system is Hamiltonian, with the same Hamiltonian as in (5). Its symplectic structure is defined by the 1-form

𝑛

𝜔 = ∑(𝑝_𝑖 𝑑𝑞^′_𝑖 + 𝑝_𝑖^′ 𝑑𝑞_𝑖).

𝑖=1

Note the “cross” pairing of coordinates and momenta — each momentum is paired with the other set of coordinates.

The finite-difference scheme obtained by the doubling method reads

                                𝑝̂𝑖 − 𝑝𝑖             𝜕𝐻′       ′,𝑞′),               𝑞̂𝑖 − 𝑞𝑖 =          𝜕𝐻′ (𝑝′,𝑞′),          𝑖 = 1,…,𝑛,

                           ⎪⎧    𝛥𝑡         = − 𝜕𝑞′_𝑖 (𝑝                            𝛥𝑡                𝜕𝑝𝑖′

⎨𝑝̂𝑖′ − 𝑝𝑖′            𝜕𝐻                             ⎪                   = −         (𝑝,̂ 𝑞)̂ ,                             ⎩      𝛥𝑡                𝜕𝑞𝑖 Or, using (6), the scheme can be written as	𝑞̂′𝑖 − 𝑞′𝑖               𝜕𝐻 = (𝑝,̂ 𝑞)̂ , 𝛥𝑡            𝜕𝑝𝑖	𝑖 = 1,…,𝑛.
𝑝̂                                  ⎧       𝑖 − 𝑝𝑖 = −𝜕𝐻′2 ,                                  ⎪      𝛥𝑡                 𝜕𝑞𝑖                                  ⎨𝑝̂𝑖′ − 𝑝𝑖′            𝜕𝐻2 ||                                  ⎪                   = −                                        𝛥𝑡                   𝜕𝑞 |	𝑞̂′𝑖 − 𝑞′𝑖       𝜕𝐻2 ||                   , = 𝛥𝑡                    𝜕𝑝𝑖 |𝑝=𝑝,𝑞=̂ 𝑞̂ 𝑞̂      ,                      𝑖 − 𝑞𝑖 = 𝜕𝐻′2 ,                   𝛥𝑡             𝜕𝑝	𝑖 = 1,…,𝑛, 𝑖 = 1,…,𝑛.	(7)
⎩                                  𝑖     𝑝=𝑝,𝑞=̂ 𝑞̂                                          𝑖

Figure 1. Time evolution of the variables 𝑝(𝑡) and 𝑝^′(𝑡) for the Jacobi oscillator

3.     Results

3.1.      Jacobi oscillator

To understand the purpose of doubling the variables, consider a simple and well-studied example — the Jacobi oscillator:

𝑑𝑝

⎧⎪ 𝑑𝑡 = 𝑞𝑟,

⎪𝑑𝑞

⎪ 𝑑𝑡 = −𝑝𝑟,

⎪

𝑑𝑟 = −𝑘2𝑝𝑞,

⎨𝑑𝑡

^⎪𝑝(0) = 0, ⎪

⎪𝑞(0) = 1, ⎪

⎩𝑟(0) = 1. It is well known that the solution of this system is periodic and can be expressed in terms of elliptic functions. Figure 1 shows the time dependence of 𝑝(𝑡) and 𝑝^′(𝑡): the amplitude of oscillations of 𝑝 grows over time, while the amplitude of 𝑝^′ decreases. The exact solution lies somewhere between these two trajectories. The behaviour shown in Fig. 1 closely resembles the classical idea of upper and lower solutions for differential equations, first introduced by Chaplygin [11].

Looking at this figure, one can expect that the average of 𝑝 and 𝑝^′ approximates the exact solution, while half the absolute difference |𝑝 − 𝑝^′| provides an estimate of the error incurred by the method.

The graph of this error estimate is shown in Fig. 2.

Nevertheless, taking arithmetic averages does not produce the behavior typically expected from a geometric integrator. This is particularly clear from the phase portrait in the (𝑝,𝑞)-plane. Figure 3 displays the trajectories of the points (𝑝,𝑞) and (𝑝^′,𝑞^′). As anticipated, (𝑝,𝑞) spirals outward to infinity,

Figure 2. Error function 𝐸(𝑝) as a function of time 𝑡

 

Figure 3. Trajectories of the points (𝑝,𝑞) and (𝑝^′,𝑞^′)

q², a.u.

 

Figure 4. Trajectory of the point (𝑝,𝑞)

 

Figure 5. Phase portrait of the Jacobi oscillator in the (𝑞,𝑝)-plane obtained using the geometric mean

whereas (𝑝^′,𝑞^′) spirals inward toward the origin. Forming the arithmetic mean — namely, the point

𝑝 + 𝑝^′ 𝑞 + 𝑞^′

                                                                (𝑝,𝑞) = (                 ,               ),

                                                                                  2              2

partially corrects this behavior (see Fig. 4), but the improvement is incomplete. The averaged trajectory still fails to lie on a closed curve, in contrast to the exact periodic orbits or those obtained with Kahan’s method. This indicates that the discretization introduces a small amount of numerical dissipation into the system.

However, if we use geometric means in our calculations, the situation improves dramatically. As clearly seen in Fig. 5, the point

(𝑝,𝑞) = (√𝑝𝑝^′,√𝑞𝑞^′)

lies exactly on the circle

                                                                           2         2

                                                                         𝑝     + 𝑞     = 1,

which is precisely the first integral of the Jacobi oscillator,

𝑝² + 𝑞² = 1.

Thus, in this case, using geometric averages exactly preserves the original system’s motion integral.

As an example of the quadratic integral preservation, the Jacobi oscillator possesses the quadratic first integral

𝑣 = 𝑝² + 𝑞².

Consequently, the doubled system admits the invariant

𝑤 = 𝑝𝑝^′ + 𝑞𝑞^′,

which is exactly conserved under the numerical scheme (3). When expressed in terms of geometric

means,

this invariant takes the form                         

                                                                         𝑤 = 𝑝      + 𝑞 ,

which explains the exact preservation of the circular phase trajectories observed in Fig. 5.

3.2.     Vanhaecke system

To see how integrals of degree higher than two behave under the discrete scheme (3), we look at a well-known example of a completely integrable Hamiltonian system — the Vanhaecke system [12– 14], first introduced by P. Vanhaecke in the 1990s. This system has two degrees of freedom and

Hamiltonian

𝐻 = 1(𝑝12 + 𝑝22) + 12(𝑞21 + 𝑞22)2 + 𝛼𝑞21 + 𝛽𝑞22 − 𝛼 − 𝛽, 2

where 𝛼 and 𝛽 are parameters of the system. This system has two polynomial first integrals:

𝐹 = (𝑝₂𝑞₁ − 𝑝₁𝑞₂)² − (2𝑞⁴₁ + 2𝑞²₁𝑞²₂ + 2𝛼𝑞₁² + 𝑝₁²)(𝛼 − 𝛽),

and

𝐺 = (𝑝₂𝑞₁ − 𝑝₁𝑞₂)² + (2𝑞²₁𝑞²₂ + 2𝑞⁴₂ + 2𝛽𝑞₂² + 𝑝₂²)(𝛼 − 𝛽),

                      Figure 6. Plots of 𝐻 and 𝐻^′                                                                                                    Figure 7. Plot of the averaged Hamiltonian

F, a.u.

                      Figure 8. Plots of 𝐹 and 𝐹^′                                                                                                          Figure 9. Plots of averaged values of 𝐹

which are in involution. Moreover, their difference is proportional to the Hamiltonian 𝐻. We integrated the system using the scheme (3) with the initial conditions

𝑞1(0) ⎧ ⎪𝑞2(0) ⎨𝑝1(0) ⎪ ⎩𝑝2(0)

= 0, = 1, = 1, = 0,

parameters 𝛼 = 1, 𝛽 = 2, and time step 𝛥𝑡 = 0.05.

Figure 6 shows the values of 𝐻 = 𝐻(𝑝,𝑞) and 𝐻^′ = 𝐻(𝑝^′,𝑞^′). Even with such a small time step, the energy changes noticeably — already in the first decimal place. However, 𝐻 and 𝐻^′ oscillate out of phase around a value close to the exact one, and — most importantly — neither shows any long-term drift upward or downward.

Averaging in different ways reduces the size of these oscillations. We tried three approaches:

• the simple average (𝐻 + 𝐻^′) — shown as the solid line in Fig. 7,
• 𝐻(𝑝,𝑞), where 𝑝 = ^𝑝+𝑝′ and 𝑞 = ^𝑞+𝑞′ — dashed line,

                                            2                            2

• 𝐻(𝑝,𝑞), where 𝑝 = √𝑝𝑝^′ and 𝑞 = √𝑞𝑞^′ — doted line.

This keeps the energy with an error of order 10⁻³, and, most importantly, the energy shows no long-term drift upward or downward.

Figure 8 shows how the integral 𝐹 behaves over time. When computing the averages (Fig. 9), we used only the first two of the three methods described above. The geometric mean was not used here because the expression for 𝐹 involves more than just squares of coordinates and momenta, which caused technical difficulties.

The same pattern is clearly visible as in the energy plot: the two copies 𝐹 and 𝐹^′ oscillate out of phase around the true value, and averaging reduces the oscillation amplitude without introducing any drift.

3.3.      Symplecticity

The scheme (7) is a partitioned (or split) finite-difference method [15, 16]. The first 2𝑛 variables

(𝑝,𝑞^′) and the last 2𝑛 variables (𝑝^′,𝑞) are updated using different formulas, which is typical for such methods. Because of this structure, the scheme is expected to be symplectic.

Theorem 11. If the Hamiltonian splits as

𝐻 = 𝑇(𝑝) + 𝑈(𝑞),

then the discrete scheme (7) is symplectic, i.e.,

𝜔̂ − 𝜔 = 𝑑𝑆

for some function 𝑆.

Proof. Let us denote the partial derivatives of 𝑇 and 𝑈 by

                                                                        𝜕𝑇                           𝜕𝑈

                                                                𝑇_𝑖 =   ,            𝑈_𝑖 =   .

                                                                        𝜕𝑝_𝑖                                        𝜕𝑞_𝑖

Recall that the discrete symplectic form is

𝑛

𝜔̂ = ∑(𝑝̂_𝑖 𝑑𝑞̂^′_𝑖 + 𝑝̂_𝑖^′ 𝑑𝑞̂_𝑖).

𝑖=1

We treat the two parts separately. First, using the update rule from (7),

𝑞̂^′_𝑖 = 𝑞^′_𝑖 + 𝛥𝑡 𝑇_𝑖(𝑝)̂ ,

so

                                  ∑𝑝̂_𝑖 𝑑𝑞̂^′_𝑖                 + 𝛥𝑡 𝑇_𝑖(𝑝)̂) = ∑𝑝̂_𝑖 𝑑𝑞^′_𝑖 + 𝛥𝑡 ∑𝑝̂_𝑖 𝑑𝑇_𝑖(𝑝)̂ .

The second term is an exact differential. Indeed,

∑𝑝̂_𝑖 𝑑𝑇_𝑖(𝑝)̂ = 𝑑(∑𝑝̂_𝑖𝑇_𝑖(𝑝)̂ − 𝑇(𝑝)̂),

so we can write

∑𝑝̂_𝑖 𝑑𝑞̂^′_𝑖 = ∑𝑝̂_𝑖 𝑑𝑞^′_𝑖 + 𝑑𝑆₁(𝑝)̂ ,

where 𝑆₁(𝑝)̂ = 𝛥𝑡(∑𝑝̂_𝑖𝑇_𝑖(𝑝)̂ − 𝑇(𝑝)̂).

Now use the other half of the scheme: 𝑝̂_𝑖 = 𝑝_𝑖 − 𝛥𝑡 𝑈_𝑖(𝑞^′), so

∑𝑝̂_𝑖 𝑑𝑞^′_𝑖 = ∑𝑝_𝑖 𝑑𝑞^′_𝑖 𝑝_𝑖 𝑑𝑞^′_𝑖 − 𝑑𝑈(𝑞^′).

Putting everything together,

∑𝑝̂_𝑖 𝑑𝑞̂^′_𝑖 = ∑𝑝_𝑖 𝑑𝑞^′_𝑖 − 𝑑𝑈(𝑞^′) + 𝑑𝑆₁(𝑝)̂ .

A similar calculation for the second part gives

∑𝑝̂_𝑖^′ 𝑑𝑞̂_𝑖 = ∑𝑝_𝑖^′ 𝑑𝑞̂_𝑖 + 𝑑𝑈(𝑞)̂ = ∑𝑝_𝑖^′ 𝑑𝑞_𝑖 + 𝑑𝑆₂(𝑝^′) + 𝑑𝑈(𝑞)̂ ,

where 𝑆    ) − 𝑇(𝑝^′)).

Adding both contributions and subtracting the original form

𝜔 = ∑(𝑝_𝑖 𝑑𝑞^′_𝑖 + 𝑝_𝑖^′ 𝑑𝑞_𝑖),

we find that all non-exact terms cancel each other out, and

𝜔̂ − 𝜔 = 𝑑𝑆

with 𝑆 = 𝑆1(𝑝)̂ + 𝑆2(𝑝′) − 𝑈(𝑞′) + 𝑈(𝑞)̂. Hence the scheme is symplectic.                                                                                

4.     Discussion

In general, it is highly convenient when a numerical method inherently includes a built-in error estimator. In contrast, methods like Runge–Kutta require such error estimation techniques to be added separately [17].

Nevertheless, taking arithmetic averages does not produce the behavior typically expected from a geometric integrator. This is particularly clear from the phase portrait in the (𝑝,𝑞)-plane. Figure 3 displays the trajectories of the points (𝑝,𝑞) and (𝑝^′,𝑞^′). As anticipated, (𝑝,𝑞) spirals outward to infinity,

whereas (𝑝^′,𝑞^′) spirals inward toward the origin.

Forming the arithmetic mean — namely, the point

𝑝 + 𝑝^′ 𝑞 + 𝑞^′

                                                                (𝑝,𝑞) = (                 ,               ),

                                                                                  2              2

partially corrects this behavior (see Fig. 4), but the improvement is incomplete. The averaged trajectory still fails to lie on a closed curve, in contrast to the exact periodic orbits or those obtained with Kahan’s method. This indicates that the discretization introduces a small amount of numerical dissipation into the system.

However, if we use geometric means in our calculations, the situation improves dramatically. As clearly seen in Fig. 5, the point

(𝑝,𝑞) = (√𝑝𝑝^′,√𝑞𝑞^′)

lies exactly on the circle

                                                                           2         2

                                                                         𝑝     + 𝑞     = 1,

which is precisely the first integral of the Jacobi oscillator,

𝑝² + 𝑞² = 1.

Thus, in this case, using geometric averages exactly preserves the original system’s motion integral.

In this sense, scheme (3) actually works better than Kahan’s method. Unlike the doubling scheme, Kahan’s method does not preserve the original integral 𝑝² + 𝑞². Instead, it conserves a modified expression dependent on 𝛥𝑡. Finding this expression for the Kahan scheme was the main result of the work [18], in [5] this expression was found according to the Lagutinsky method.

Exact preservation is important in many cases. For example, one conservation law for the spinning top simply says that the sum of the squares of the direction cosines is always 1 [19]. Kahan’s method preserves a modified expression [20], i.e., adds tiny corrections that don’t have a clear geometric meaning. Because scheme (3) preserves this integral exactly, it looks very promising for building geometric integrators — methods that respect the underlying geometry of the problem.

Overall, we can say that the scheme (3) keeps all integrals of the system close to their exact values allowing only small oscillations around them. This is a very good property and suggests that the scheme may be reliable for long-time simulations. At the same time, it is important to note that the energy integral is not preserved exactly, no matter how we average it.

It should also be mentioned that the Vanhaecke system has cubic right-hand sides, so it cannot be directly discretized by Kahan’s method. Even if it could, there are no general results guaranteeing that Kahan’s method preserves arbitrary polynomial integrals. The only case that is well understood is the preservation of the energy (Hamiltonian) itself [2].

Energy conservation in symplectic integrators has been studied extensively. It is well known that the original Hamiltonian is not preserved exactly. However, one can construct a modified Hamiltonian 𝐻^(𝑚) that is conserved up to any desired order of accuracy [15, 21]. Because of this, numerical experiments typically show the energy oscillating around a constant value that is close to — but not exactly equal to — the true energy. This is precisely what we observed in theVanhaecke system. Looking at Fig. 7, one might even guess that the different averaging strategies we used correspond to modified Hamiltonians of different orders.

5.     Conclusion

In this work, we studied reversible finite-difference schemes for dynamical systems based on the doubling method introduced by V. N. Abrashin and S. N. Sytova. We carried out a detailed analysis of algebraic integral preservation for two benchmark systems: the classical Jacobi oscillator and the Vanhaecke system.

First of all, quadratic integrals are proven to preserve exactly (Th. 9). In particular, the quadratic integrals for the Jacobi oscillator are preserved exactly if we replace the squares of the quantities with the geometric mean. For comparison’s sake, Kahan’s method does not preserve the original quadratic integral 𝑝² + 𝑞², but conserves a modified expression dependent on 𝛥𝑡. Exact preservation is important for geometrical interpretation of expressions like the sum of the squares of the direction cosines.

For Hamiltonian systems, this method is symplectic in the sense of the theorem 11. Thus the original Hamiltonian is not preserved exactly, but one can construct a modified Hamiltonian 𝐻^(𝑚) that is conserved up to any desired order. It seems that the different averaging strategies we used correspond to modified Hamiltonians of different orders. In any way, the scheme keeps all integrals of the Vanhaecke system close to their exact values allowing only small oscillations around them. Thus the scheme may be reliable for long-time simulations.

About the authors

Lyubov O. Lapshenkova

RUDN University

Email: lapshenkova-lo@rudn.ru
ORCID iD: 0000-0002-1053-4925

Assistant of Department of Computational Mathematics and Artificial Intelligence

6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation

Kseniia S. Mashkovtseva

RUDN University

Email: 1132226438@rudn.ru
ORCID iD: 0009-0002-6927-4467

Student of Department of Computational Mathematics and Artificial Intelligence

6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation

Alina A. Trusova

RUDN University

Email: 1132246715@rudn.ru
ORCID iD: 0009-0008-5140-5629

Student of Department of Computational Mathematics and Artificial Intelligence

6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation

Mikhail D. Malykh

RUDN University; Joint Institute for Nuclear Research

Author for correspondence.
Email: malykh-md@rudn.ru
ORCID iD: 0000-0001-6541-6603
Scopus Author ID: 6602318510
ResearcherId: P-8123-2016

DSc., Head of Department of Computational Mathematics and Artificial Intelligence of RUDN University; Senior Researcher of Joint Institute for Nuclear Research

6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation; 6, Joliot-Curie St., Dubna, Moscow Region, 141980, Russian Federation

References

Supplementary files