On the properties of numerical solutions of dynamical systems obtained using the midpoint method

Full Text

Introduction One of the most common and popular mathematical models is the Cauchy problem for an autonomous system of ordinary differential equations. Analytical methods make it possible to find the algebraic integrals of motion [1] for such systems, and numerical methods allow approximate plotting the particular solutions [2]. Consider an autonomous system of differential equations in an affine space of dimension ̇ = (). (1) Here = (1, … , ) is a point in the affine space, = (1, … , ) is a set of rational functions belonging to ℚ(). Except the cases where it leads to ambiguity, we will use the vector notation ̇ = (). By the algebraic integral of motion of this system, we mean the algebraic function , constant on any particular solution of the system (1), i.e., satisfying the equation 1 1 + ⋯ + = 0. (2) It can be shown that the existence of an algebraic integral implies the existence of a rational integral, therefore only rational integrals are considered below [3]. The finite difference method is a standard numerical method for solving systems of ordinary differential equations [2]. The finite-difference scheme for solving the system of equations (1) describes the transition from the value of taken at some instant of time to the value of taken at the next instant of time + Δ. This new value will be denoted below by .̂ Of course, by the difference scheme for the system (1) we understand a correspondence in some sense approximating the system of differential equations, rather an arbitrary correspondence between the variables and .̂ Usually, by this approximation we mean that the system of equations defining the difference scheme tends to the original system as Δ → 0. Definition 1. By a particular solution of the system (1), found using the finite-difference scheme, we mean a finite or infinite sequence of points (0), (1), (2), … , (), … of the -dimensional space (or subset considered in it), the first element of which is taken arbitrarily, and each next element is obtained from the previous one according to the difference scheme: (+1) = (), = 0, 1, 2, … This approximate particular solution will be associated with the exact solution of the Cauchy problem ̇ = (), |=0 = (0). Analytical and numerical methods cannot always be reconciled. It often turns out that the algebraic integral of motion is known, but a difference scheme is used, which does not preserve this integral. Therefore, in numerical experiments, it is often nothing to do but to observe with regret how the quantity which remains constant on the exact solution, changes step by step. Very frequently, the integrals of motion express fundamental laws of nature, the violation of which introduces new properties into the mathematical model under consideration, or trivial geometric relationships, the violation of which makes it difficult to interpret the results of integration. The idea of contructing finite-difference schemes exactly preserving the integral of motion of dynamical systems was proposed in late 1980s in the papers by Yu. B. Suris [4] and Cooper [5], approaching the problem from different sides, namely, Yu. B. Suris from composing difference schemes for Hamiltonian systems that preserve a symplectic structure, and Cooper from preserving integrals. As a result, a large family of Runge-Kutta schemes was discovered that preserve all quadratic integrals of any dynamical system and the symplectic structure of the Hamiltonian system [6]. The simplest of this class of schemes is the midpoint scheme. By construction, the approximate solutions found using this scheme retain all quadratic integrals. Traditionally, the question of ‘improvement’ of convergence due to the conservation of integrals remained in the focus of attention. In this article, we intend to clarify what other qualitative properties of the exact solution are inherited by the approximate one. For completeness, we give an elementary proof of Cooper’s theorem. 1. Conservation of quadratic integrals Let be an integral of the system of Eqs. 1. According to Lagrange theorem about the mean value, (1̂ , … , ̂ ) - (1, … , ) = ∑ (1, … ) ⋅ Δ, =1 where the derivatives are calculated at the point , lying somewhere in the segment connecting the points and ,̂ i.e., = + Δ, ∈ (0, 1). The Lagrange theorem ensures the existence of a suitable function (, )̂ , taking the values between 0 and 1 at the real values of , .̂ However, it provides no method to find this function. If we knew such a function, we could easily construct a difference scheme that preserves the manifold = 0 exactly. Indeed, let Then Δ Δ = ( + Δ) , = 1, 2, … , . (3) (1̂ , … , ̂ ) - (1, … , ) = ∑ (1, … ) ⋅ (1, … , ), =1 and this expression is exactly equal to zero due to (2). For some classes of functions, the Lagrange theorem allows a constructive formulation. Lemma 1. If is a polynomial from ℂ[1, … ], the degree of which does not exceed 2 inclusively, then (1̂ , … , ̂ ) - (1, … , ) = ∑ (1, … ) ⋅ Δ, where = (̂ + ) /2. =1 Proof. Let be a new auxiliary variable and () = (1 + Δ1, … , + Δ). According to the Lagrange theorem, there is such value ∈ (0, 1), that (1) - (0) = ′(), (4) or, in more detail, (1̂ , … , ̂ ) - (1, … , ) = ∑ =1 (1, … ) ⋅ Δ, where = + Δ. We have to prove that = 1/2. By the hypothesis of the lemma, is a polynomial whose degree does not exceed 2, therefore ″ is a polynomial with respect to whose degree does not exceed 2, i.e. () = 2 + + , ′() = 2 + . But then the Lagrange formula (4) reduces to + = 2 + , from which it immediately follows that = 1/2. D The finite-difference scheme (3) at = 1/2, i.e., the scheme Δ Δ = 1̂ ( + 1 , … , 2 ̂ + ) , = 1, 2, … , , (5) 2 is called the midpoint scheme. Lemma 1 immediately leads to the following theorem. Theorem 1 (Cooper, 1983; [5]). The midpoint finite-difference scheme (5) not only approximates the system (1), but also preserves all linear and quadratic integrals of this system. 2. Harmonic oscillator Consider the simplest dynamical system ⎧{̇ = -, ⎨{⎩̇ = , (6) describing a harmonic oscillator. This system is Hamiltonian, the energy conservation law for it yields the integral 2 + 2 = . Usually it is immediately said without further ado that the solutions of this system describe periodic rotations along concentric circles in the phase plane . The standard discretization does not allow to preserve this description. Consider for definiteness, the explicit Euler scheme ⎧{̂ - = -Δ, ⎨{⎩̂ - = Δ. Figure 1 presents the curve in the phase plane, obtained using the Euler method instead of a unit circle. Thus it is easy to describe the difference between the exact and approximate solutions. Instead of closed curves, a spiral in the phase plane appears, or, in terms of classical mechanics, the time discretization leads to gradual increase of the system energy. y 1.5 1 0.5 x -1.5 -1 -0.5 0.5 1 -0.5 -1 Figure 1. Solution of the initial value problem for Eqs. (6) using the Euler method (solid) and the meanpoint method (dashed), Δ = 0.1, 100 steps are maid Now let us consider the meanpoint scheme ̂ - = -(̂ + ) Δ , 2 ̂ - = (̂ + ) Δ 2 (7) and try to understand the qualitative difference between the approximate solution and the exact one. According to Cooper’s theorem, this difference scheme exaclty preserves this integral. Therefore, in the phase plane we get ovals that seemingly do not differ from a circle (see Figure 1). However, these curves are still not closed, and the motion along them cannot be considered perodic, since the values of , do not repeat. Therefore, it could be executed that the absence of closedness and periodicity completely allows one to distinguish between the approximate solution and the exact one. In fact, it is possible to achieve the solution relativity by the appropriate step choice. Theorem 2. If we take for the minimal positive root of the equation (1 + )2 = (1 + 2), (8) which in terms of trigonometric functions can be expressed as = tan , then the calculation according to the meanpoint finite-difference scheme (7) with the step Δ = 2 = 2 tan , in steps leads to the initial values of , . Remark 1. For us it is convenient to use the transcendent formula = tan , however, one could do without it using purely algebraic means. To emphasize this fact, in the statement of the theorem we have preserved the algebraic equation, for which this number is a root. For proof, let us express ,̂ ̂ in terms of , , denoting for brevity Δ Then = . 2 or, in the matrix form ̂ + ̂ = - , { - ̂ + ̂ = + 1 ̂ ( ) ( ) = ( 1 - ) ( ) . - 1 Inverting the matrix, we get ̂ 1 ) = 1 2 (1 - (9) ̂ 1 + 2 ) ( ) 1 now the proof of theorem 2 is reduced to the proof of the following purely algebraic lemma. Lemma 2. If we take for the minimal positive root of equation (8), which in terms of trigonometric functions can be written as = tan , then 2 is a unit matrix. Proof. Set for brevity 1 = √ (1 -) . 1 + 2 1 Assume that we performed steps of the finite-difference scheme, having started from the values , at a certain value = 0, and finished with , at = Δ. Then ( -1 ) = 2 ( ) ) = ⋯ = 2 ( . -1 The eigenvalues of the matrix are 1 + 1 - 1 = √ , 2 = √ . 1 + 2 1 + 2 When a matrix is raised to a power, its eigenvalues are raised to this power, too. The equation can be satisfied assuming that 2 1 = 1 (10) 1 + √ = cos 1 + 2 or + sin = tan . The equality (10) can be rewritten in the form of the algebraic equation (8); in this way the minimal positive root of this equation has been found. 2 At real Δ from the equality (10) the equality 2 = 1 follows. In this case both eigenvalues of the matrix 2 are equal to 1. The eigenvectors of matrix are also eigenvectors of any power of this matrix, therefore, in this case two linearly independent eigenvectors correspond to the unit eigenvalue. This is possible only if the matrix 2 coincides with the unit matrix. D The proved theorem gives rise to the following definition. Definition 2. The particular solution , ′, (2), … found using a certain finite-difference scheme approximating the system (1), at a certain numerical value of the step Δ will be called periodic, if for some ∈ ℕ () = . The number Δ will be called a period of this particular solution. Theorem 2 means that the midpoint scheme (7) yields periodic solutions at a number of step values that form a descending sequence: Δ3 = 2 tan 3 √ = 2 3 = 3.46, Δ4 = 2 tan = 2, 4 Δ5 = 2 tan 5 = 2 √-2 2 √ √ 5 + 5 = 1.45, Δ6 = 2 tan = 6 3 … 3 = 1.15, √ √ Δ20 = 2 tan = 2 20 … 5 - 2 √2 5 + 5 + 2 = 0.31, converging to zero. The corresponding sequence of periods Δ = 2 tan = 2 + 23 3 1 2 + … converges to the exact solution period 2. Thus, the second most important qualitative property of the exact solution, i.e., its periodicity, can be preserved, too. Moreover, even for small values of , in this way it is possible to get a solution that is qualitatively similar to the exact one. For example, for = 10 and the step value Δ = 2 tan 10 2 √ = √-10 5 + 25 = 0.64 5 we arrive at a periodic scheme with the period = 10Δ, that differs from the period of the exact solution by a noticeable value of - 2 = 0.21. In the phase plane a regular decagon is obtained that obviously differs from a circle, too. However, even in this case the trajectory in the phase plane is exactly closed and the motion is periodic. In other words, the approximate solution found using the midpoint scheme turns out to be qualitatively similar to the exact one, being rather different quantitatively. In our opinion, when using the harmonic oscillator model, there is not enough reason to think that the characteristic time scale, taken as physically small, is really small from the point of view of computational mathematics. Therefore, the real reason in favor of using the continuous form of the Newton equations rather than the discrete one is that the traditional method of discretization leads to a violation of the fundamental laws of mechanics (the law of energy conservation) and to a solution, whose properties differ from the expected periodicity. The use of periodic difference schemes removes this difficulty. 3. A system of coupled oscillators A system of coupled oscillators can be described as a Hamiltonian system with the Hamiltonian = ∑ ∑ - ∑ ∑ . =1 =1 =1 =1 It is assumed that the kinetic and potential energy ∑ , =1 ∑ =1 are described by positively defined quadratic forms. In the matrix form the system is written as ̈⃗ = -.⃗ Since the matrices and are symmetric and positively defined, the generalized eigenvalue problem ⃗ = ⃗ yields eigenvectors ⃗ , … , ⃗ 1 that form a basis ℝ, orthonormalized with the weight , corresponding to 1 positive eigenvalues 1 = 2, … , = 2 . In this basis ⃗ = ∑ () ⃗ , where = ⃗ .⃗ =1 The differential equation in these variables is separated into independent equations ⃗ =1 =1 ⃗ or, due to the basis orthogonality, ̈ + 2 = 0, = 1, … , . Each of these equations separately describes the vibration of a harmonic oscillator with the frequency . Therefore, the positive numbers 1, … , are called eigenfrequencies of the system of coupled oscillators. Now from the formula ⃗ = ∑ () ⃗ =1 it is seen that the oscillations of the coupled system will be a superposition of harmonic oscillations at the eigenfrequencies. Now it is important for us that this system has algebraic integrals 2̇ + 22 = , which in old variables can be written as (⃗ )⃗̇ 2 + 2(⃗ )⃗ 2 = , = 1, 2, … , . (11) Thus, the initial system has quadratic integrals. These integrals will be called partial integrals. Lemma 3. If the eigenfrequencies are incommensurable, then any algebraic integral of a system of coupled oscillators is expressed algebraically in terms of partial integrals. The proof of this lemma is somewhat lengthy and is not given here. In essence, it relies on constructions from the proof of Bruns theorem proposed by Painlevé [7]. Remark 2. The condition of incommensurability of frequencies is essential. If 1 = 2, then the expression 1̇ 2 - 12̇ = is an additional algebraic integral. In fact, the general solution of the system is 1 = 1 cos 1 + 2 sin 1, 2 = 3 cos 1 + 4 sin 1, … therefore 1̇ 2 - 12̇ = (23 - 14)1, i.e., is really a constant. Thus, without a separate study, it cannot be ruled out that in degenerate cases the system has additional integrals of motion that are not preserved by the midpoint scheme. According to the theorem 1, the midpoint scheme preserves all partial integrals, and, therefore, all algebraic integrals of the system. The step Δ can always be chosen so that the normal oscillation at a certain natural frequency is described by a periodic particular solution. However, this step depends on the harmonic number. Therefore, it is impossible to choose a time step at which the oscillation, in the framework of the continuous model, which is a superposition of periodic normal oscillations (mixed oscillation), would be the sum of periodic particular solutions. Thus, for the case of linear systems of ordinary differential equations, the midpoint scheme is a difference scheme that preserves all algebraic integrals. Moreover, for periodic particular solutions, it is possible to select a step in a purely algebraic way so that the approximate solution becomes periodic in the sense of the definition 2. 4. Elliptic oscillator We now proceed to the nonlinear case. Solid state dynamics provides many excellent examples of autonomous systems with periodic solutions. The simplest of them are integrated in Jacobi elliptic functions [8]. By the definition of Jacobi functions [9], = sn , = cn , = dn is a particular solution of the nonlinear autonomous system of differential equations ̇ = , ̇ = -, ̇ = -2 (12) with the initial conditions = 0, = = 1 at = 0. This autonomous system has two quadratic integrals of motion 2 + 2 = const and 22 + 2 = const. (13) This system is convenient for us, because its properties are not only well studied, but the solutions themselves are easily accessible in any computer algebra system, e.g., in Sage [10]. The system (12) is a convenient test for studying the conservatism of finite-difference schemes used in the common numerical solvers of ordinary differential equations. We, in cooperation with Yu.A. Blinkov, conducted tests on the following systems: 1. lsoda: Real-valued Variable-coefficient Ordinary Differential Equation solver, with fixed-leading-coefficient implementation. It provides automatic method switching between implicit Adams method (for non-stiff problems) and a method based on backward differentiation formulas (BDF) (for stiff problems). Source: http://www.netlib.org/odepack; 2. vode: Real-valued Variable-coefficient Ordinary Differential Equation solver, with fixed-leading-coefficient implementation. It provides implicit Adams method (for non-stiff problems) and a method based on backward differentiation formulas (BDF) (for stiff problems). Source: http://www. netlib.org/ode/vode.f; 3. dopri5, dop853: This is an explicit Runge-Kutta method of order (4)5 due to Dormand & Prince (with stepsize control and dense output). In Figures 4-7 it is well seen that in all cases the value of 2 + 2 increases or decreases almost linearly. Only in the first solver (see Figures 2, 3) “random” fluctuations are observed, but with a trend towards linear growth. 8e-5 +1 kp2 +q2 -1k∞ = 8.1e-05 numeric exact 6e-5 p2 +q 2 4e-5 2e-5 0 0 200 400 600 800 1000 t Figure 2. The value of 2 + 2 - 1, calculated with the solver Lsoda using the method based on backward differentiation formulas From an analytical point of view, this system is remarkable because any particular solution of it is representable as the ratio of two everywhere convergent series in powers of [8]. From the point of view of the finite difference method, this system is remarkable because it can be approximated by a difference scheme, namely, the midpoint scheme (5), which preserves its integrals exactly. 0.000175 +1 kp2 +q2 -1k∞ = 1.9e-04 numeric exact 0.000150 0.000125 0.000100 0.000075 0.000050 0.000025 0.000000 0 200 400 600 800 1000 t Figure 3. The value of 2 + 2 - 1, calculated with the solver Lsoda using the method based on the default settings 1.0014 1.0012 numeric exact kp2 +q2 -1k∞ = 1.4e-03 1.0010 p2 +q 2 1.0008 1.0006 1.0004 1.0002 1.0000 0 200 400 600 800 1000 t Figure 4. The value of 2 + 2 - 1, calculated with the solver Vode using the method based on backward differentiation formulas (BDF) with the choice of Jacobian However, for nonlinear differential equations, the organization of the transition from layer to layer according to the midpoint scheme requires solving a system of nonlinear algebraic equations. The organization of calculations according to such a scheme is a subject of discussion. 0.0007 0.0006 +1 kp2 +q2 -1k∞ = 7.3e-04 numeric exact 0.0005 p2 +q 2 0.0004 0.0003 0.0002 0.0001 0.0000 0 200 400 600 800 1000 t Figure 5. The value of 2 + 2 - 1, calculated with the solver Vode using the method based on backward differentiation formulas (BDF) without the choice of Jacobian 0.00010 +9.999×10-1 kp2 +q2 -1k∞ = 9.2e-05 0.00008 p2 +q 2 0.00006 0.00004 0.00002 numeric exact 0 200 400 600 800 1000 t Figure 6. The value of 2 + 2 - 1, calculated with the solver Vode using the default settings The midpoint scheme for the system (12) is written as follows: ⎧̂ - { ̂ + ̂ + - = 0, { Δ 2 2 { ̂ - + ̂ + ̂ + (14) { ⎨ Δ 2 = 0, 2 { { ̂ - 2 ̂ + ̂ + + ⎩ Δ 2 = 0. 2 0.0000100 +9.9999×10-1 kp2 +q2 -1k∞ = 3.6e-06 numeric exact 0.0000095 0.0000090 q 0.0000085 p 0.0000080 0.0000075 0.0000070 0.0000065 0 200 400 600 800 1000 t Figure 7. The value of 2 + 2 - 1, calculated with the solver Dopri5 To pass from one layer to another it is necessary from the given numerical values of Δ and , , to find ,̂ ,̂ ,̂ having solved the system of nonlinear equations. In this case, in addition to the required root, close to the values of , , at small Δ, extraneous roots are also obtained. Let us investigate this system in Sage, replacing the hat notation with doubling the letter, e.g., using the notation pp instead of .̂ The system of equations (14) generates an ideal in the ring ℚ[, , , ̂, ,̂ ,̂ Δ, ]. In Sage this ideal can be specified by the following code: sage: vars=var('p,q,r,pp,qq,rr,dt,k') sage: K=QQ[vars] sage: eqs=[4*(pp-p)-(q+qq)*(r+rr)*dt,\ 4*(qq-q)+(p+pp)*(r+rr)*dt,\ 4*( rr-r)+k^2*(p+pp)*(q+qq)*dt] sage: J=K.ideal(eqs) Now we can compose an equation connecting ̂ with , , by means of elimination ideals [11]: sage: J.elimination_ideal([K(qq),K(rr)]).gens()[0] as 5dt44 + 34p̂ dt44 + 23p̂ 2dt44 - 22p̂ 3dt44 - 3p̂ 4dt44 - p̂ 5dt44 +162dt32 + 32p̂ dt32 + 16p̂ 2dt32 - 323dt22 - 642dt22 -322p̂ dt22 - 642p̂ dt22 + 32p̂ 2dt22 + 32p̂ 3dt22 - 642dt2 - 642p̂ dt2 +256dt + 256 - 256p̂ = 0. Thus, one value of , , on the previous layer generally corresponds to 5 different values of ̂ rather than one value. Since the higher degrees of ̂ enter this equation with the factor Δ, for Δ → 0 only one of these roots tends to a finite value, which is obviously ̂ = , and the four others go to infinity. Thus, applying the midpoint scheme to a nonlinear system requires solving a nonlinear equation and then choosing from its roots one root that is ’close’ to the values in the previous step. Traditionally, such equations are solved numerically by iterative methods, and the values of , , are used as the first approximation for ,̂ ,̂ .̂ If the step Δ is sufficiently small, we can expect the fast convergence of the described method. At the same time, it is not possible to control the error of the iterative method due to the extreme cumbersomeness of the known estimates, and instead the number of iterations is simply fixed. However, as correctly noted in Numerical Recipies [12], there are no universal numerical methods for solving systems of algebraic equations. Symbolic methods, primarily the Gröbner basis technique, allow us to propose a different approach, circumvent this difficulty, and implement calculations according to the midpoint scheme in a different way: · the first stage (symbolic analysis of the finite-difference scheme): find equations having the form (̂, , , , Δ) = 0, (,̂ , , , Δ) = 0, (,̂ , , , Δ) = 0, that follow from Eq. (14); the first of these equations was written out above; · the second stage (computations with floating point): to pass from layer to layer solve three uncoupled equations to find ̂, ,̂ ̂ and select the roots close to the values of , , at the previous layer. Thus, from the numerical solution of the system of nonlinear equations, we proceed to the solution of one algebraic equation of the 5th degree with numerical coefficients. The methods of numerical solution are well established and in practice give errors close to the errors in calculating a radical from a number. Remark 3. Unfortunately, the algorithms implemented in Sage do not give us, together with the roots, an estimate of this error. The question of the precision solution of algebraic equations with real coefficients is currently discussed at all conferences on numerical methods and computer algebra (MMA’2018, CASC’2019). We hope that in the near future this issue will be completely closed. We implemented both approaches to the implementation of the midpoint method and made sure that the first approach gives the same result faster. Figures 8-9 show the elliptic sine plot calculated in Sage using the built-in algorithm for calculating elliptic functions and directly using the midpoint method. It is clearly seen that even at extremely large values of , the oscillation amplitude does not decrease. Thus, when using the midpoint scheme: · both algebraic integrals of motion are preserved exactly, · rounding error does not accumulate in a way that is noticeable in numerical experiments, · the periodic nature of the motion is preserved. A theoretical study of the last two issues faces significant difficulties, which, in our opinion, make us search not for any conservative schemes, but for schemes in which the transition from layer to layer requires solving linear equations. One such scheme for the elliptic oscillator (12) was specified by us in [13]; in the general case of the system (1) there are fundamental obstacles to the existence of such schemes [14]. p p 1 1 0.5 t 5 10 15 20 0.5 0 t 990 992 994 996 998 1000 -0.5 -0.5 -1 -1 2 Figure 8. sn (, 1 ) , = 200 2 Figure 9. sn (, 1 ) , = 200 5. Comparing midpoint scheme with other symmetric schemes One could think that the conservation of integrals is the result of the scheme symmetry. Therefore, for completeness, we compare the midpoint scheme with another popular symmetric second-order scheme: Δ Δ ()̂ + () = 2 . (15) The results of our numerical experiments, presented in the Figure 10, show that this scheme does not preserve the integrals, but gives them values that oscillate with a constant period around some fixed value. This can be accepted for preserving the integrals of the motion on average. 0.0006 +1 kp2 +q2 -1k∞ = 6.2e-04 numeric exact 0.0005 0.0004 p2 +q 2 0.0003 0.0002 0.0001 0.0000 990 992 994 996 998 1000 t Figure 10. The value of 2 + 2 - 1, calculated using the scheme (15) This effect can be explained as follows. The scheme (15) for the Jacobi system yields Δ = ̂ ̂ + 2 Δ, Δ = - ̂ ̂ + 2 Δ, Δ = -2 ̂ ̂ + 2 Δ. (16) Therefore Δ(2 + 2) = (̂ + )Δ + ( ̂+ )Δ = = [(̂ + )(̂ ̂ + ) - ( ̂+ )(̂̂ + )] Δ Δ = 2 = [̂ ̂ + ̂ - ̂ - ̂]̂ 2 Δ = Δ = ( ̂- ̂)Δ 2 = (Δ - Δ)Δ = 2 Δ2 = -[(̂̂ + ) + (̂ ̂ + )]Δ = 4 2 3 = [(̂ + )̂ ̂ + (2 + 2)](̂ ̂+ ) or, expanding in series in powers of Δ Δ 8 (17) 2 Δ(2 + 2) = (2 + 2) Δ3 2 + … . For small Δ in the plot of 2 + 2 versus the step number (or, which is similar, versus ) we will see periodic oscillations of the first term 2 (2 + 2) Δ3 . 2 Thus, the oscillation of the value of the integral 2 +2 near its exact value 1 observed in the numerical experiment is not related to the conservativeness of the scheme, but is due to the fact that the main term in the expansion of the increment of this integral in power of depends on the values of , , , approximating periodic functions. The considered representation scheme is a perfect example of the scheme symmetric with respect to permutation of and .̂ Conclusion For the study of dynamical systems with quadratic integrals, the midpoint scheme is perfect. This scheme preserves all the integrals of these systems precisely, that is, discretization does not introduce any dissipativity or antidissipativity into the model. Linear problems have just quadratic integrals, and the calculation according to the midpoint scheme requires solving linear equations, so the issue of constructing conservative difference schemes for linear differential equations can be considered closed. It worth noting that the Kepler problem after passing to proper time and introducing the Runge-Lenz-Laplace vector turns into a linear dynamical system, therefore, it is possible to construct for it a conservative finite-difference scheme [15]. Conservation of integrals leads to the preservation of a number of qualitative properties of the model, e.g., the periodicity of solutions and the closedness of phase trajectories. The calculation according to the midpoint scheme for nonlinear systems, even if all their integrals are quadratic, makes it necessary to solve a nonlinear system of algebraic equations at each step, which significantly complicates both the calculation and the study of the properties of approximate solutions. Overcoming this difficulty is the main challenge of the theory of difference schemes for ordinary differential equations.

About the authors

Vladimir P. Gerdt

Joint Institute for Nuclear Research

Author for correspondence.
Email: gerdt@jinr.ru

Doctor of Physical and Mathematical Sciences, Full Professor at the Joint Institute for Nuclear Research (JINR) where he is the head of the Group of Algebraic and Quantum Computations

Joliot-Curie St. 6, Dubna, Moscow Region 141980, Russian Federation

Mikhail D. Malykh

Peoples’ Friendship University of Russia (RUDN University)

Email: malykh-md@rudn.ru

Candidate of Physical and Mathematical Sciences, assistant professor of Department of Applied Probability and Informatics of Peoples’ Friendship University of Russia

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

Leonid A. Sevastianov

Joint Institute for Nuclear Research; Peoples’ Friendship University of Russia (RUDN University)

Email: sevastianov-la@rudn.ru

professor, Doctor of Physical and Mathematical Sciences, professor of Department of Applied Probability and Informatics of Peoples’ Friendship University of Russia (RUDN University)

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

Yu Ying

Peoples’ Friendship University of Russia (RUDN University); Kaili University

Email: yingy6165@gmail.com

postgraduate student of Department of Applied Probability and Informatics of Peoples’ Friendship University of Russia (RUDN University); assistant professor of Department of Algebra and Geometry, Kaili University

Miklukho-Maklaya St. 6, Moscow 117198, Russian Federation; Kaiyuan Road 3, Kaili 556011, China

References