 # On conservative difference schemes for the many-body problem

A new approach to the construction of difference schemes of any order for the many-body problem that preserves all its algebraic integrals is proposed. We introduced additional variables, namely, distances and reciprocal distances between bodies, and wrote down a system of differential equations with respect to coordinates, velocities, and the additional variables. In this case, the system lost its Hamiltonian form, but all the classical integrals of motion of the many-body problem under consideration, as well as new integrals describing the relationship between the coordinates of the bodies and the additional variables are described by linear or quadratic polynomials in these new variables. Therefore, any symplectic Runge-Kutta scheme preserves these integrals exactly. The evidence for the proposed approach is given. To illustrate the theory, the results of numerical experiments for the three-body problem on a plane are presented with the choice of initial data corresponding to the motion of the bodies along a figure of eight (choreographic test).

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## Introduction

One of the main continuous models is a dynamic system described by an autonomous system of ordinary differential equations, that is, a system of equations of the form

 dxidt=fi(x1,…,xn),i=1,2,…n, (1)

where is an independent variable commonly interpreted as time, are the coordinates of a moving point or several points. In practice, the right-hand sides are often rational or algebraic functions of coordinates or can be reduced to such form by a certain change of variables.

In the framework of the classical approach, the solutions of Hamiltonian dynamical systems are found in two stages: first, algebraic integrals of motion (Whittaker, § 164) are determined and then, based on Liouville’s theorem (Whittaker, § 148), the problem is reduced to the calculation of integrals of algebraic functions referred to as quadratures. However, construction of these integrals is usually not enough to reduce differential equations to quadratures.

For example, the classical problem of bodies Marshal consists in finding solutions to an autonomous system of ordinary differential equations

 mi¨→ri=n∑j=1γmimjr3ij(→rj−→ri),i=1,…,n (2)

Here

is the radius vector of the

-th body, is its masses, is the distance between the -th and -th bodies, and is the gravitational constant.

Let us for brevity denote the components of the velocity of the -the body as , , and and combine them in vector . This problem has 10 classical integrals of motion:

1. the momentum conservation (3 scalar integrals)

 n∑i=1mi→vi=const
2. the angular momentum conservation (3 scalar integrals)

 n∑i=1mi→vi×→ri=const
3. the center-of-mass inertial motion (3 scalar integrals)

 n∑i=1mi→ri−tn∑i=1mi→vi=const
4. the energy conservation (1 scalar integral)

 n∑i=1mi2|→vi|2−γ∑i,jmimjrij=const

Bruns (Whittaker, § 164) proved that every algebraic integral of motion in this system is expressed algebraically in terms of the above 10 classical integrals. Note, that these algebraic integrals do not provide reduction of the problem to quadratures.

When it is not possible to integrate a dynamic system analytically, it is done numerically. Among the methods of numerical analysis, the finite difference method is the most universal. The numerical solution of the systems of autonomous equations (1) is performed via a difference scheme. Every difference scheme describes a transition from the value of taken at some time to the value of

taken at the next moment in time

. Hereafter we denote these new values as . Usually this relationship is described by an algebraic equation similar, in its form, to a differential equation. For example, an explicit Euler scheme is described by a difference scheme

 ^x−x=f(x)Δt.

However, much more sophisticated schemes can also be used, including implicit ones, which have the form of algebraic equations that not solved with respect to , and multi-stage schemes which use additional variables. By this reason, we do not restrict ourselves to any essential conditions in regarding the form of difference scheme.

## 1 Finite-difference schemes preserving all the algebraic integrals of motion

Let us start from the basic definition. A difference scheme for autonomous system (1) is said to preserve the integral

 g(x1,…,xn)=C

of this system if the equality

 g(^x)=g(x)

is valid for any non-special choice of . We call a scheme conservative if it preserves all the algebraic integrals of motion of the considered autonomous system.

Standard finite-difference schemes, and first of all explicit Euler and Runge-Kutta schemes, do not preserve the algebraic integrals of motion, but yield an approximate solution that is close to the exact solution at sufficiently small values of the time step in the considered time interval Wanner. However, after its discretization, the model can be changed qualitatively, e.g., because of the violation of energy conservation or geometric relationships and sometimes the appearance of dissipation (cf.Lubich; Casc-2019).

The idea to construct difference schemes that precisely preserve the integrals of motion in dynamical systems arose in the late 1980s. In the works of Cooper Cooper and Yu.B. Suris Suris-1988; Suris-1990 a large family of Runge-Kutta schemes was discovered that preserve all the polynomial second-order integrals of motion that do not depend on explicitly in any dynamical system and also on the symplectic structure in Hamiltonian systems Lubich; Sanz-Serna.

Strictly speaking, these studies are related only to integrals that are not explicitly dependent on time, and when they are presented in the literature Wanner, these results are often limited to the case of quadratic forms. In reality, these reservations are not needed. In fact, one can always introduce a new independent variable and rewrite the system (1) in the form

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩dxidτ=fi(x1,…,xn),i=1,2,…n,dtdτ=1. (3)

Any integral of motion of the original system, quadratic in , is also a quadratic integral of the system (3) that does not contain

. If necessary, a linear transformation can be converted to a quadratic form. By Cooper’s theorem, this integral is preserved by any symplectic Runge-Kutta scheme applied to (

3). Since the right-hand side of the last equation of the system (3) is identically equal to unity, in any Runge-Kutta scheme, the projections of the slopes on the - axis will be equal to , and, hence, a step made with respect to is equal to the step along . As a result, the transition from layer to layer will be carried out in accordance to the same formulas as govern the transition done according to the Runge-Kutta scheme with the same Butcher matrix, but already written for the original system (1). Therefore, the quadratic integral is preserved in the calculations according to the symplectic Runge-Kutta scheme applies the original system (1).

If a dynamical system admits only transcendental integrals of motion, one can try to construct difference schemes preserving these integrals. For example, recently Wan et al. (Wan-2017, § 5.2.1) proposed a finite difference scheme for the Lotka-Volterra system, preserving its transcendental integral. The schemes obtained in this way cannot be expressed in terms of arithmetic operations only; in particular, the scheme proposed by Wan et al. involves logarithms. Since we intend to apply computer algebra methods to analyze of difference schemes in future research, it is of fundamental importance for us to deal with algebraic equations only. For this reason, we avoid discussion of the use of transcendental integrals.

In the many-body problem (2), 9 out of 10 classical integrals are quadratic in the coordinates of the bodies, their velocities and time . Thereby, symplectic schemes exactly preserve 9 from 10 integrals, but they do not preserve the energy integral.Thus, the discretization introduces dissipation in the system.

The three-body problem describes the motion of a planet and its satellite around the Sun, and in this situation two characteristic parameters appear that have the dimension of time: the (approximate) period of revolution of the planet around the Sun and the period of revolution of the satellite around the planet. One can hope to obtain reasonable results by taking a time step substantially smaller than these two periods. At the same time, the positions of bodies after many periods can be of interest. In this case, we need in solving the problem of many bodies in a large time scale. From general considerations, one could hope that symplectic integrators would allow greater accuracy in calculations for large times. The advantages and drawbacks of symplectic integrals in application to the many-body system long-term behavior have been discussed many times Hernandez-2019.

The first finite-difference scheme for the many-body problem, preserving all classical integrals of motion, was proposed in 1992 by Greenspan Greenspan-1992; Greenspan-1995-1; Greenspan-1995-2; Greenspan-2004 and independently in somewhat different form by Simo and González Simo-1993; Graham. The Greespan’s scheme is a kind of combination of the midpoint method and discrete gradient method McLachlan; Christiansen preserving a number of symmetries of the initial problem Greenspan-1992.

In the 2000s, several general approaches to construction of difference schemes were proposed that preserve all algebraic integrals of motion. To our opinion, these methods can be divided into two large groups: the methods based on modification of the Runge-Kutta schemes, and the methods used certain transformations of the original differential system. The first group includes the parametric symplectic partitioned Runge-Kutta method Brugnano-2012; Wang-2012, in which the Runge-Kutta method is considered with the order of approximation chosen such that there is a free parameter in the procedure of determining the Butcher matrix coefficients. This parameter can be chosen to preserve a chosen quadratic integral.

Naturally, for practical implementation, just this step is the most difficult part of the procedure. If the scheme is chosen to be simplectic, then all quadratic integrals are preserved automatically. For example, for the many-body problem (2) one can construct an implicit finite difference scheme of any order. On the contrary, it is possible to take an explicit Runge-Kutta scheme with one or several parameters and define these parameters so that given integrals of motion become preserved. The time step itself is often used as a parameter. In this way a few explicit schemes have been constructed preserving both quadratic and non-quadratic integrals Buono-2002; Zhang-2020-2. The Hamiltonian boundary value methods Brugnano-2010; Brugnano-2012-2; Brugnano-2012-3, can be considered as belonging to the same group. In these methods, the initial system is approximated according to the Galerkin procedure and then the scalar products are approximated by finite sums. As a result, the obtained method resembles the Runge-Kutta method, which preserves the integral of energy not exactly, but with any predetermined order of accuracy (Brugnano-2012-2, Th. 3).

The invariant energy quadratization method (IEQ, proposed recently Yang et al. Yang-2016; Yang-2017) can be considered to be a member of the second group. Based on this method, for a number of Hamiltonian systems including the two-body problem (Kepler problem), Hong Zhang et al. Zhang-2020 constructed finite-difference schemes conserving the energy. The IEQ method suggests finding such a change of variables, after which the energy becomes a quadratic form, so that the standard simplectic Runge-Kutta schemes, written for this transformed system, conserve its energy. The procedure admits increasing the number of unknowns and the loss of non-Hamiltonian form of equations.

In the past, every effort was made to perform replacements that preserve the symplectic structure of the many-body problem. In particular, when studying the simple collision of two bodies, the aim was to find not just a regularizing transformation, but a canonical transformation. Due to this, the famous Weierstrass theorem on simple collisions was considered as an obstacle Siegel. There is, however, a simple way to find a regularizing transformation, which was proposed by Burdet Burdet and Heggie Heggie and is thoroughly described in the book by Marchal (Marshal, ch.  6). For this purpose, the simplectic structure of the problem should be abandoned and additional variables can be introduced. Therefore, the construction proposed by Hong Zhang et al. Zhang-2020 within the IEQ method for the two-body problem seems to be quite natural for the considered class of problems.

One of the delicate features of IEQ method is the preservation of constraints. The fact is that when additional variables are introduced, equations describing the relation between the coordinates and velocities of bodies, on the one hand, and the auxiliary variables, on the other hand. In the IEQ method these equations are not necessarily quadratic, so that these constraints are not preserved exactly. We argue that it is possible if one looks at such a change of variables that makes all integrals of motion and all constraint equations to be described by quadratic functions. Therefore, in this paper we want to present the solution of the following problem.

###### Problem 1.

For the -body problem (2) construct another system containing additional variables and having the following properties:

1. It possesses a sufficient number of algebraic integrals of motion in order to express additional variables in terms of the coordinates of the bodies,

2. With some choice of constant values in these integrals, its solutions coincide with those of the original system,

3. It has integrals of motion, which, if one takes into account the relationship between additional variables and the coordinates of the bodies, are transformed into 10 classical integrals of the many-body problem,

4. All integrals of motion of the system are quadratic in the coordinates and velocities of bodies, as well as in additional variables.

## 2 Rationalization of the n-body problem

First of all, let us eliminate irrationality by introducing new variables , related to the coordinates by equation

 r2ij−(xi−xj)2−(yi−yj)2−(zi−zj)2=0

Differentiating this relation with respect to , we get

 rij˙rij=(xi−xj)(ui−uj)+(yi−yj)(vi−vj)+(zi−zj)(wi−wj)

or

 rij˙rij=(→ri−→rj)⋅(→vi−→vj)

Let us rewrite the differential system (2) of second-order equations into the system of the first-order equations in the variables

 x1,…,zn,u1,…,wn,r12,…,rn−1,n

This system consists of three coupled subsystems:

 ˙→ri=→vi,i=1,…,n (4) mi˙→vi=n∑j=1γmimjr3ij(→rj−→ri),i=1,…,n (5) ˙rij=1rij(→ri−→rj)⋅(→vi−→vj),i,j=1,…,n;i≠j. (6)

Every solution to the many-body problem (2) satisfies this system, but, in general, the converse is not true. Not every solution to (4)-(6) satisfies the relation

 r2ij=(xi−xj)2+(yi−yj)2+(zi−zj)2

and, thus, the systems (2) and (4)-(6) have distinct solution sets. Therefore, generally speaking, the new system may have fewer integrals than the original. Howeverm, in the given case, we can prove that the new system inherits the set of classical integrals of the many-body problem (2).

###### Theorem 1.

System (4)-(6) has 10 classical integrals of the problem (2) and additionally the integrals

 r2ij−(xi−xj)2−(yi−yj)2−(zi−zj)2=const.
###### Proof.

To check the conservation of expressions indicated in the theorem we will calculate the derivations of these expressions.

• The momentum conservation

 ddtn∑i=1mi→vi=n∑i=1mi˙→vi=n∑i=1n∑j=1γmimjr3ij(→rj−→ri)=0. (7)
• The angular momentum conservation

 ddtn∑i=1mi→vi×→ri=n∑i=1mi˙→vi×→ri=n∑i=1n∑j=1γmimjr3ij(→rj−→ri)×→ri=0. (8)
• The center-of-mass inertial motion

 ddtn∑i=1mi→ri−tn∑i=1mi→vi=n∑i=1mi→vi−n∑i=1mi→vi=0. (9)
• The energy conservation

 ddtn∑i=1mi2(u2i+v2i+w2i) =n∑i=1mi→vi⋅˙→vi=n∑i=1n∑j=1γmimjr3ij→vi⋅(→rj−→ri) =γ∑ijmimjr3ij(→vi⋅(→rj−→ri)+→vj⋅(→ri−→rj)) =γ∑ijmimjr3ij(→vi−→vj)⋅(→rj−→ri)

and, due to (6),

 ddt∑i,jmimjrij=−∑i,jmimjr2ij˙rij=−∑i,jmimjr3ij(→ri−→rj)⋅(→vi−→vj)

 r2ij−(xi−xj)2−(yi−yj)2−(zi−zj)2=const (10)

follow from equations (6), since the derivative

 ddt(r2ij−(xi−xj)2−(yi−yj)2−(zi−zj)2)

is equal to

 2rij˙rij−2(→ri−→rj)(→vi−→j)=0.

Now we have an autonomous system with rational right-hand side, all integrals of which are quadratic polynomials, except the following rational expression

 n∑i=1mi2(u2i+v2i+w2i)−γ∑i,jmimjrij=const

which corresponds to the energy conservation.

## 3 System with quadratic polynomial integrals

To get rid of the denominators in the energy conservation law, we introduce new additional variables such that

 rijρij=1.

Note that this relation again is a quadratic polynomial what allow us to obtain a quadratic polynomial integral instead of the rational one.

Differentiating this relation with respect to , we obtain

 rij˙ρij+ρij˙rij=0

or

 ˙ρij=−ρijr2ij(→ri−→rj)⋅(→vi−→vj).

Now we rewrite the system (4) - (6) as the system of first-order equations with respect to the extended set of unknowns

 x1,…,zn,u1,…,wn,r12,…,rn−1,n,ρ12,…,ρn−1,n

Now we obtain the system of four coupled subsystems

 ˙→ri=→vi,i=1,…,n (11) mi˙→vi=n∑j=1γmimjρijr2ij(→rj−→ri),i=1,…,n (12) ˙rij=1rij(→ri−→rj)⋅(→vi−→vj),i,j=1,…,n;i≠j (13) ˙ρij=−ρijr2ij(→ri−→rj)⋅(→vi−→vj),i,j=1,…,n;i≠j (14)

As well as the system (4)-(6) the dynamical system system (11)-(14) inherits the conservation laws of (2).

###### Theorem 2.

The system (11)-(14) has 10 classical integrals of the many-body problem (2), namely, Eqs. (7)-(9), the energy conservation in the form

 n∑i=1mi2(u2i+v2i+w2i)−γ∑i,jmimjρij=const,

 rijρi,j=const,i≠j. (15)
###### Proof.

In perfect analogy to the proof of Theorem 1 we perform explicit differentiation and simplify the obtained expressions. For Eqs. (7)-(10) it is done exactly as the proof of Theorem 1. Since now the energy conservation has a slightly different form, we present the relevant computation

 ddtn∑i=1mi2(u2i+v2i+w2i) =n∑i=1mi→vi⋅˙→vi=n∑i=1n∑j=1γmimjρijr2ij→vi⋅(→rj−→ri) =γ∑ijmimjρijr2ij(→vi⋅(→rj−→ri)+→vj⋅(→ri−→rj)) =γ∑ijmimjρijr2ij(→vi−→vj)⋅(→rj−→ri)

and, in view of (14),

 ddt∑i,jmimjρij=−∑i,jmimjρijr2ij(→ri−→rj)⋅(→vi−→vj).

The validity of the conservation law (15) follows from equations (13)-(14):

 ddtrijρij=˙rijρij+rij˙ρij=ρijrij(→ri−→rj)⋅(→vi−→vj)−rijρijr2ij(→ri−→rj)⋅(→vi−→vj)=0.

Since the differential equations of the considered system were derived by differentiating relations (10) and (15), the appearance of the above additional integrals is obvious. ∎

In general, it is not possible to state that the system (11)-(14) has no other algebraic integrals, functionally independent of those listed above. However, a solution to the many-body problem satisfies to the extended system (11)-(14), so that any algebraic integral of motion of the extended system remains constant on solutions of the many-body problem. By the Bruns theorem (Whittaker, § 164) on a manifold

 r2ij−(xi−xj)2−(yi−yj)2−(zi−zj)2=0,rijρij=1,i≠j,

such integral is expressed algebraically in terms of the classical integrals.

According to Theorem 2, the autonomous differential system (11) - (14) containing additional variables and , has the following properties:

• It has the quadratic integrals of motion

 r2ij−(xi−xj)2−(yi−yj)2−(zi−zj)2=const

and

 rijρij=const,

which allow to express the additional variables and in terms of the coordinates of the bodies.

• If the constants in these integrals are chosen in such a way that

 r2ij−(xi−xj)2−(yi−yj)2−(zi−zj)2=0

and

 rijρij=1,

then solutions to system coincide with those to the initial one (2).

• The new system has quadratic integrals of motion, which, taking into account the relationship between additional variables and the coordinates of the bodies, turn into 10 classical integrals of the many-body problem.

Thus, the constructed system posses all the properties listed in the Problem 1.

## 4 Conservative schemes for N body problem

Since all classical integrals of the many-body problem, as well as the additional integrals, are quadratic in their variables, any symplectic Runge-Kutta difference scheme, including the simplest midpoint one, that is

 ^x−x=f(^x+x2)dt, (16)

preserves all these integrals. In particular, the midpoint scheme written for the system (11)-(14), preserves all its algebraic integrals exactly and is invariant under permutations of bodies and time reversal. It is not difficult to create high-order schemes which preserve all integrals of motion in the many-body problem.

The autonomous system of differential equations (11) - (14) preserves the symmetry of the original problem with respect to permutations of bodies and time reversal. As noted above, the midpoint scheme is invariant under these symmetries.

At each step of the scheme, new values will be determined not only for the coordinates and velocities, but also for the auxiliary quantities and . If at the initial moment only the coordinates and velocities of the bodies were specified and the auxiliary variables were defined by the equalities

 rij=√(xi−xj)2+(yi−yj)2+(zi−zj)2,ρij=1rij,

then these equalities are preserved exactly (up to the signs of the radicals) due to the fact that the auxiliary integrals (10) and (15) are quadratic and exactly preserved under the action of midpoint scheme. Therefore, the quantities and remain the distances between bodies and the inverse distances between bodies.

## 5 Numerical experiments with planar three-body problem

We developed a numerical algorithm based on the midpoint scheme YY-2019 and realized it in SageMath sage; SageMath. Discussion of this algorithm and its program implementation is beyond the scope of this paper and will be presented in other paper. Now we illustrate the above considerations by a numerical example. For simplicity, we consider the planar dimensionless problem of the motion of three bodies of equal mass with .

Before starting the numerical experiments, it is necessary to discuss the issue of organizing the calculations by means of the midpoint scheme. The problem is that this scheme is implicit, so that at each step it is necessary to solve the system of algebraic equations (16).

We are going to use the method of simple iterations to solve the system of nonlinear algebraic equations (16) with respect to : starting from the sequence

 ^x(n+1)=x+f(^x(n)+x2)dt,n=0,1,2,… (17)

is constructed. If this sequence has a limit , then it is a root of the system of equations (16

). In practice it is inconvenient to use the known estimates for

, since they are rather cumbersome, particularly, for systems of equations with a large number of unknowns. However, we use the midpoint scheme in order to preserve the integrals of motion, therefore, we will monitor that the increment of the integrals of motion does not go beyond the given boundaries ( in our examples) rather than that is small. When conducting numerical experiments, we often encountered divergence of the method. In this case we reduce the step .

### 5.1 Test on Lagrange solutions Figure 1: Lagrange test, step dt=1/10. Trajectories of three bodies at 0 Figure 2: Lagrange test, step dt=1/10. Solid line shows the exact solution, dots represent the approximate one.

As is well known since Euler Marshal; Montgomery-2001, the planar three-body problem has two families of particular solutions that can be described in elementary functions. The first of them, traditionally associated with the name of Lagrange, can be described as follows: the bodies form a regular triangle with sides of a fixed length , which rotates around the center of mass with constant speed. We took and a quite noticeable step , nevertheless, the nature of the movement remained the same, which is clearly visible along the trajectories (Fig. 1). The difference from the exact solution is barely noticeable in Fig. 2.

At the same time, the distances between the points remain constant at the level of rounding error. The error in calculating the distance between the center of gravity and the bodies, as well as in determining the energy and angular momentum, is much larger, at the level of the specified . Thus, throughout the entire considered time interval, the bodies fall at the vertices of a regular triangle, which rotates around the center of mass.

### 5.2 Choreographic test Figure 3: Choreographic test, step dt=1/10. Trajectory of the first body for 0 Figure 4: Choreographic test, step dt=1/10. Dependence of energy H on time for approximate solutions found using our algorithm (blue) and the rk4 scheme (red).

In 1993, Chris Moore discovered Moor, as part of a numerical experiment, a new periodic solution of the three-body problem, in which three bodies write out the eight; later Chenciner and Montgomery Montgomery-2000; Montgomery-2001 gave a justification to this fact. The initial conditions are known only approximately. We will use the values found numerically by Carles Simó and given in Montgomery-2000.

First of all, we note that the figure eight really turns out (Fig. 3). However, at a coarse step it is noticeable that this eight moves in time. It is not clear whether this movement can be stopped by the refinement of initial data.

To compare our algorithm and the standard algorithms based on explicit schemes, we took the fourth-order standard Runge-Kutta scheme (rk4) with the same coarse step . Of course, the calculation according to this scheme is faster, and the approximation is better by two orders in magnitude. Nevertheless, the error in calculating the energy from rk4 grows almost linearly like and becomes observable soon, and according to our scheme it does not leave the error corridor (Fig. 4).

### 5.3 Test with small distances Figure 5: Test with collision, step dt=6⋅10−3. Trajectories of the three bodies for 0 Figure 6: Test with collision, step dt=6⋅10−3. Trajectories of the three bodies for 0

In all the tests considered above, the bodies remained at a considerable distance from each other, so there was no computational singularity rather typical for the many-body problem, i.e., the situation when bodies pass close to each other.

At the initial moment of time, the following conditions are taken: the bodies lie on the line , the first body at the point , the second at , the third at . The first body is at rest, and the initial velocities of the other two bodies are directed along the -axis and are equal to and , respectively. Then in the process of computation, it turns out that the bodies several times pass close to each other (Fig. 5). Moreover, the applicability conditions of the iteration method make it necessary to reduce the step. Of course, we cannot assume that at this moment a true collision occurs, since, according to the Weierstrass theorem, the collision in the many-body problem with arbitrary initial data is improbable Siegel.

The question of the stability of the many-body problem with respect to the initial data is very complex and we would not want to invade this area here. However, it seemed important to us to note that under small perturbations of the initial data our scheme yields trajectories close to the original ones. We checked this in several numerical experiments, e.g., in the test with possible collision, we increased the velocity of the third body by and compared the trajectories (Fig. 6). It is clearly seen that the new trajectories are close to the old ones, and the point of possible collision turns out to be almost at the initial position.

## Conclusion

The major result of our work is Theorem 2, which yields a solution to Problem 1 (Section 1) and can be used to construct finite-difference schemes exactly preserving all its algebraic integrals and invariant under the permutations of bodies and time reversal. The midpoint scheme is an example of such scheme when it applied (11)-(14). It should be emphasized that in this way we can construct high-order schemes that preserve all integrals of motion in the many-body problem.

The key problem with the application of these schemes in practice, of course, is their implicitness, which can be overcome by both numerical and numerical-analytical methods. Obviously, errors in solving a nonlinear system of algebraic equations describing one step in a scheme can completely level out all the advantages of the proposed scheme Zhang-2020-2. Therefore, in our opinion, this issue requires further study, including numerical experiments. The numerical example given above is intended only to illustrate theoretical results.

However we believe that the conservative difference schemes constructed above for the many-body problem will not only make it possible to carry out calculations for large times, but will also make it possible to qualitatively investigate the properties of solutions of the many-body problem by the finite difference method.

It should also be noted that when comparing the schemes, the focus of the researchers’ attention was always on the quantitative proximity of the exact and approximate solutions, while the question of preserving the qualitative properties of the exact solution has remained insufficiently studied, although the tempting possibility to ‘‘determine the nature of the dynamic process using only rough calculations with a large grid pitch’’ according to the symplectic scheme was noted in Gevorkyan-2013. Meanwhile, for a linear oscillator, the simplest Runge-Kutta symplectic scheme, the midpoint scheme, allows not only accurate preservation of the energy integral, but also obtaining an exact periodic approximate solution and closed polygons as an analogue of closed phase trajectories of the continuous model YY-2019.

In other words, after discretization, a model is obtained that inherits almost all the algebraic and qualitative properties of a continuous model. In the theory of partial differential equations, discretizations inheriting certain properties of differential equations and operations on the dependent variables included in them, are called

mimetic (i.e., imitative) or compatible Shashkov-1996; Arnold-2006; Castillo-2013; Vega-2014. Therefore, we find it appropriate to consider the construction of conservative difference schemes in the context of the more general and more important issue of constructing mimetic difference schemes for the many-body problem, of course, clarifying the concept itself.

Nevertheless, on the basis of the midpoint scheme and idea of introducing auxiliary variables for the many-body problem, a difference scheme is explicitly presented that preserves all the classical integrals of motion precisely and is invariant under permutations of bodies and time reversal. In further studies we plan to explore numerical and algebraic properties of this system. We believe that this scheme inherits at least a part of the algebraic and group-theoretical properties of exact solutions of the many-body problem.

This work is supported by the Russian Science Foundation (grant no. 20-11-20257).