On the properties of numerical solutions of dynamical systems obtained using the midpoint method Текст научной статьи по специальности «Математика»

Computational modeling and simulation

Research article

UDC 517.9

DOI: 10.22363/2658-4670-2019-27-3-242-262

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

Vladimir P. Gerdt1, Mikhail D. Malykh2, Leonid A. Sevastianov1,2, Yu Ying2,3

1 Joint Institute for Nuclear Research Joliot-Curie St. 6, Dubna, Moscow Region 141980, Russian Federation 2 Department of Applied Probability and Informatics Peoples' Friendship University of Russia (RUDN University) Miklukho-Maklaya St. 6, Moscow 117198, Russian Federation 3 Department of Algebra and Geometry

Kaili University Kaiyuan Road 3, Kaili 556011, China

(received: November 20, 2019; accepted: December 23, 2019)

The article considers the midpoint scheme as a finite-difference scheme for a dynamical system of the form x = f(x). This scheme is remarkable because according to Cooper's theorem, it preserves all quadratic integrals of motion, moreover, it is the simplest scheme among symplectic Runge-Kutta schemes possessing this property.

The properties of approximate solutions were studied in the framework of numerical experiments with linear and nonlinear oscillators, as well as with a system of several coupled oscillators. It is shown that in addition to the conservation of all integrals of motion, approximate solutions inherit the periodicity of motion. At the same time, attention is paid to the discussion of introducing the concept of periodicity of an approximate solution found by the difference scheme.

In the case of a nonlinear oscillator, each step requires solving a system of nonlinear algebraic equations. The issues of organizing computations using such schemes are discussed. Comparison with other schemes, including those symmetric with respect to permutation of x and x.

Key words and phrases: conservative finite-difference schemes, dynamical systems, Sage, Maple

© Gerdt V. P., Malykh M.D., Sevastianov L.A., YingYu., 2019

This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/

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 n

X = f(x). (1)

Here x = (xx,..., xn) is a point in the affine space, f = (fx, ■■■, fn) is a set of rational functions belonging to x). Except the cases where it leads to ambiguity, we will use the vector notation x = f(x). By the algebraic integral of motion of this system, we mean the algebraic function g, constant on any particular solution of the system (1), i.e., satisfying the equation

£ + + £ =° (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 x taken at some instant of time t to the value of x taken at the next instant of time t + At. This new value will be denoted below by x. 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 x and x. Usually, by this approximation we mean that the system of equations defining the difference scheme tends to the original system as At ^ 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

x(0), x(1), x(2) ,.._,XM ,.._

of the n-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:

x

(rn+l) = £(m), m = 0,1,2,

This approximate particular solution will be associated with the exact solution of the Cauchy problem

x = f(x), x\t=0 = x(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 g be an integral of the system of Eqs. 1. According to Lagrange theorem about the mean value,

n

g(Xi,..., xn) 9(%i, •••, %n) ^ y ) ' ,

i=1 i

where the derivatives are calculated at the point c, lying somewhere in the segment connecting the points x and x, i.e., ci = xi + dAxi, d £ (0,1).

The Lagrange theorem ensures the existence of a suitable function d(x, x), taking the values between 0 and 1 at the real values of x,x. 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 g = 0 exactly. Indeed, let

AxA = fz (Xi + 9Axi), i = 1,2,...,n. (3)

Then

n

g(Xi,..., xn) g(Xi,..., xn) ^ ^ ~dx^(ci, "' C'a) ' ^(ci, "', C'a^

i=1 i

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 g is a polynomial from C^,...xn], the degree of which does not exceed 2 inclusively, then

n

g(Xi,..., x^a) g(Xi,..., xn) ^ ^ (Ci,... cn) • ax.,

i=1 i

where ci = (xi + xi) /2.

Proof. Let u be a new auxiliary variable and

G(u) = g(xi + uAxx,..., xn + uAxn). According to the Lagrange theorem, there is such value d £ (0,1), that

G(1)-G(0) = G' (9), (4)

or, in more detail,

n

g(Xi,..., x^a) g(Xi,..., xn) ^ ^ (c-i, ..• cn) • AXi,

i=1 i

where ci = xi + 9 Ax..

We have to prove that d = 1/2.

By the hypothesis of the lemma, g is a polynomial whose degree does not exceed 2, therefore G" is a polynomial with respect to u whose degree does not exceed 2, i.e.

G(u) = au2 + bu + c, G'(u) = 2au + b.

But then the Lagrange formula (4) reduces to a + b = 2ad + b, from which it immediately follows that 0=1/2. □

The finite-difference scheme (3) at d = 1/2, i.e., the scheme

AXi =.U i = 1,2,...,n, (5)

At Ji\ 2 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

X = -y, (6)

.y = x,

describing a harmonic oscillator. This system is Hamiltonian, the energy conservation law for it yields the integral x2 + y2 = C.

Usually it is immediately said without further ado that the solutions of this system describe periodic rotations along concentric circles in the phase plane xy.

The standard discretization does not allow to preserve this description. Consider for definiteness, the explicit Euler scheme

{x — x = —yAt, y y = xAt.

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.

Figure 1. Solution of the initial value problem for Eqs. (6) using the Euler method (solid) and the meanpoint method (dashed), At = 0.1, 100 steps are maid

Now let us consider the meanpoint scheme

•At ,At

x — x = —{y + y) — , y — y = (x + x)— (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 x, y 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 a the minimal positive root of the equation

(1 + ia)2N = (1 + a2 )N, (8)

which in terms of trigonometric functions can be expressed as

a = tan

then the calculation according to the meanpoint finite-difference scheme (7) with the step

At = 2a = 2 tan ^, in N steps leads to the initial values of x,y.

Remark 1. For us it is convenient to use the transcendent formula a = tan jj, 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 x, y in terms of x,y, denoting for brevity

At

Then

or, in the matrix form

— _ a. 2

x + ay _ x — ay, — ax + y _ ax + y

1 a|[X|_/l —a if x -a W \yJ \a 1 ) \yy

Inverting the matrix, we get

2

X \ 1 I 1 I X

y I 1 + OL2 \ a 1 I \y

(9)

now the proof of theorem 2 is reduced to the proof of the following purely algebraic lemma.

Lemma 2. If we take for a the minimal positive root of equation (8), which in terms of trigonometric functions can be written as

n

a = tan —,

then A2N is a unit matrix.

Proof. Set for brevity

A =

\f\+a2

Assume that we performed N steps of the finite-difference scheme, having started from the values x,y at a certain value t = 0, and finished with xN, yN at t = NAt. Then

) = A2 ( XN-1

W \yN-1

The eigenvalues of the matrix A are

1 + ia

Ai =

= - = A2N

V1+

a9 =

1 — ia

a2

V1+

a2

When a matrix is raised to a power, its eigenvalues are raised to this power, too.

The equation

X2N = 1 (10)

can be satisfied assuming that

1 + ia n n

= cos — + I sin —

V1+

a2

N

N

or

■k

a = tan —.

N

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.

At real At from the equality (10) the equality \22N = 1 follows. In this case both eigenvalues of the matrix A2N are equal to 1. The eigenvectors of matrix A 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 A2N coincides with the unit matrix.

□

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 At will be called periodic, if for some

N

(N) X == X.

The number NAt 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:

At3 ж = 2 tan — 3 = 2\13 = 3.46,

At 4 ж = 2 tan — 4 = 2,

A t5 n = 2 tan — 5 2 л/5+ 5 = 1.45

A t6 n = 2 tan — 6 =2^ = 1.15,

At20 = 2 tan — = 2^5-2 ^2^5 + 5 + 2 = 0.31,

converging to zero. The corresponding sequence of periods

n 2n3 1

NAtN = 2 N tan - = 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 N, in this way it is possible to get a solution that is qualitatively similar to the exact one.

For example, for N = 10 and the step value

At = 2 tan — = 2 J-10 л/б + 25 = 0.64 10 5 v

we arrive at a periodic scheme with the period T = 10At, that differs from the period of the exact solution by a noticeable value of

T- 2ж = 0.21.

In the phase plane xy 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

n n n n

H = mz3y% y3 - ky xtxj.

i=1 j=1 i=1 j=1

It is assumed that the kinetic and potential energy

n n

^ y mijViyj, ^ y

i=1 i=1

are described by positively defined quadratic forms.

In the matrix form the system is written as Mx = —Kx. Since the matrices K and M are symmetric and positively defined, the generalized eigenvalue problem Kf = \M.£ yields n eigenvectors that form a basis , orthonormalized with the weight M, corresponding to n positive eigenvalues \1 = u2,..., Xn = u2a. In this basis

n i=1

where zi = tjMx.

The differential equation in these variables is separated into n independent equations

n n

i=1 i=1 or, due to the basis orthogonality, Zj + u2zrj = 0, j = 1,... ,n.

Each of these equations separately describes the vibration of a harmonic oscillator with the frequency u>j. Therefore, the positive numbers u1,... ,un are called eigenfrequencies of the system of coupled oscillators. Now from the formula

n

1 =

i=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 n algebraic integrals

4 + ^ zj = G3,

which in old variables can be written as

ÇijM!)2 + u2 (ijMx)2 = Cj, j=1,2,...,n. (11)

Thus, the initial system has n 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 n 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 Painleve [7].

Remark 2. The condition of incommensurability of frequencies is essential. If u1 = u2, then the expression

¿1^2 — Z1 ¿2 = C

is an additional algebraic integral. In fact, the general solution of the system is

z1 = C1 cosu1t + C2 sinu11, z2 = C3 cosu11 + C4 sinu11,...

therefore z1z2 — z1 z2 = (C2C3 — C1C4)u1, 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 At 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],

p = sn t, q = cn t, r = dn t

is a particular solution of the nonlinear autonomous system of differential equations

p = qr, q = —pr, r = —k2pq (12)

with the initial conditions p = 0, q = r = 1 at t = 0.

This autonomous system has two quadratic integrals of motion

p2 + q2 = const and k2 p2 + r2 = 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 p2 + q2 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

6e-5 + 4e-5 2e-5

0

Figure 2. The value of p2 + q2 — 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 t [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.

pf + cf - 1||„ =8.1e-05

0 200 400 600 800 1000

t

0.000175 0.000150 0.000125 0.000100 0.000075 0.000050 0.000025 0.000000

Figure 3. The value of p2 + q2 — 1, calculated with the solver Lsoda using the method

based on the default settings

1.0014 1.0012 1.0010

2 1.0008

+

2 1.0006 1.0004 1.0002 1.0000

Figure 4. The value of p2 + q2 — 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.

+1

11^+ q - 1||„ = 1.9e-04

numeric ....... exact

^ jMihi r

f/Nf^

K i A j|

/ wy fW

M

NN

0 200 400 600 800 1000

t

|p* + q - 1||„ = 14e-03

—:—i

---- exact

0 200 400 600 800 1000

t

0.0007 0.0006 0.0005

2 0.0004

+

(N

^ 0.0003 0.0002 0.0001 0.0000

Figure 5. The value of p2 + q2 — 1, calculated with the solver Vode using the method based on backward differentiation formulas (BDF) without the choice of Jacobian

0.00010 0.00008

2 0.00006

+

(N

a

0.00004 0.00002

Figure 6. The value of p2 + q2 — 1, calculated with the solver Vode using the default settings

+1

||p2 + q2 - 1||„ =7.3e-04

- numeric

—- exact

200

400

600

800

1000

0

t

+9.999X10"

kp2 + q - 1||„ =9.2e-05

T* U Wv 4 i «* J

v \

y

- numeric ....... exact —i- V V

200

400 600

t

800

1000

0

The midpoint scheme for the system (12) is written as follows:

'p—p q + qr + r

At 2 2

q — q p + pr + r

At

+

= 0, = 0,

r — r + ,2 p + pq + q = o At + 2 2

0.0000100 0.0000095 0.0000090 0.0000085 , 0.0000080 0.0000075 0.0000070 0.0000065

Figure 7. The value of p2 + q2 — 1, calculated with the solver Dopri5

To pass from one layer to another it is necessary from the given numerical values of At and p, q, r to find p, q, r, having solved the system of nonlinear equations. In this case, in addition to the required root, close to the values of p, q, r at small At, 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 p. The system of equations (14) generates an ideal J in the ring ^[p, q, r,p, q, r, At, k]. 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 p with p, q, r by means of elimination ideals [11]:

sage: J.elimination_ideal([K(qq),K(rr)]).gens()[0]

as

p5 dt4k4 + 3p4pdt4k4 + 2p3p2 dt4k4 - 2p2p3 dt4k4 - 3pp4 dt4k4 - p5 dt4k4 + 16p2qrdt3k2 + 32pqrpdt3k2 + !6qrp2dt3k2 - 32p3dt2k2 - 64pq2dt2k2 -32p2pdt2k2 - 64q2pdt2k2 + 32pp2dt2k2 + 32p3dt2k2 - 64pr2d,t2 - 64r2pdt2 +256qrdt + 256p - 256p = 0.

Thus, one value of p, q, r on the previous layer generally corresponds to 5 different values of p rather than one value. Since the higher degrees of p enter this equation with the factor At, for At ^ 0 only one of these roots tends

to a finite value, which is obviously p = p, 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 p, q, r are used as the first approximation for p, q, r. If the step At 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 Grobner 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 P(p,p,q,r,At) = 0, Q(q,p,q,r, At) = 0, R(r,p, q, r, At) = 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 p,q,r and select the roots close to the values of p, q, r 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 t, 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].

Figure 9. sn (t, ■§) , N = 200

Figure 8. sn (t, , N = 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:

At' 2 • ( ^

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

0.0005

0.0004

+ 0.0003 a.

0.0002 0.0001 0.0000

Figure 10. The value of p2 + q2 — 1, calculated using the scheme (15)

This effect can be explained as follows. The scheme (15) for the Jacobi system yields

qr + qr . pr + pr . 10PQ + PQ * /„„x

Ap = H 2 At, Aq = ~ 2 At, Ar = —k2 At. (16)

Therefore

A(p2 + q2) = (p + p)Ap + (q + q)Aq =

= [(p + p)(qr + qr) — (q + q)(pr + pr)]A =

r......,At

= [pqr + pqr — qpr — qpr]— =

= (PQ — pq)ata^ = (pAq — qAp)Ar=

At2

= —[p(pr + pr) + q(qr + qr)]Ar=

k,2 At3

= [(pp + qq)r+ (p2 + q2 )r](pq + pq)—-— (17)

8

or, expanding in series in powers of At

k2 At3

A(p2 + q2) = (p2 + q2)pqr—^--+ - .

For small At in the plot of p2 + q2 versus the step number (or, which is similar, versus t) we will see periodic oscillations of the first term

(p2 + q2 )pqr^ ^ .

Thus, the oscillation of the value of the integral p2 + q2 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 t depends on the values of p, q, r, approximating periodic functions. The considered representation scheme is a perfect example of the scheme symmetric with respect to permutation of x and x.

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.

Acknowledgments

The studies of difference schemes that inherit the properties of a continuous model carried out by V. P. Gerdt and L. A. Sevastianov were supported by the RUDN 5-100 program. The calculation results presented in Figures 2-7 and 10 were obtained by V. P. Gerdt together with Yu. A. Blinkov in the SymPy system and were first presented at the PCA'2019 conference [16]. The authors are grateful to Yu. A. Blinkov for the materials provided. Calculations using the midpoint scheme were performed in the Sage system by M. D. Malykh and Yu Ying as part of studies supported by RFBR grant No. 18-51-18005.

References

[1] A. Goriely, Integrability and nonintegrability of dynamical systems. Singapore; River Edge, NJ: World Scientific, 2001.

[2] E. Hairer, G. Wanner, and C. Lubich, Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations. Berlin Heidelberg New York: Springer, 2000.

[3] L. Konigsberger, Die Principien der Mechanik. Leipzig: Teubner, 1901.

[4] Y. B. Suris, "Hamiltonian methods of Runge-Kutta type and their variational interpretation [Gamil'tonovy metody tipa Runge-Kutty i ikh variatsionnaya traktovka]," Matematicheskoye modelirovaniye, vol. 2, no. 4, pp. 78-87, 1990, in Russian.

[5] G. J. Cooper, "Stability of Runge-Kutta methods for trajectory problems," IMA Journal of Numerical Analysis, vol. 7, pp. 1-13, 1987. DOI: 10.1093/imanum/7.1.1.

[6] J. M. Sanz-Serna, "Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more," SIAm REVIEW, vol. 58, pp. 3-33, 1 2016. DOI: 10.1137/151002769.

[7] P. Painlevé, "Mémore sur les intégrales du problème des n corps," in Œuvres de Paul Painlevé. 1975.

[8] V. V. Golubev, Vorlesungen über Differentialgleichungen im Komplexen. Leipzig: VEB Deutscher Verlag der Wissenschaften, 1958.

[9] P. F. Byrd and M. D. Friedman, Handbook of elliptic integrals for engineers and scientists. Springer, 1971. DOI: 10.1007/978-3-64265138-0.

[10] W. A. Stein. (2015). Sage Mathematics Software (Version 6.7), The Sage Development Team, [Online]. Available: http://www.sagemath.org.

[11] D. A. Cox, J. Little, and D. O'Shea, Ideals, varieties, and algorithms: an introduction to computational algebraic geometry and commutative algebra, 4th ed. Springer, 2015. DOI: 10.1007/978-3-319-16721-3.

[12] Numerical Recipes: The Art of Scientific Computing, 3rd ed. Cambridge University Press, 2007, 1256 pp.

[13] E. A. Ayryan, M. D. Malykh, L. A. Sevastianov, and Yu Ying, "Finite difference schemes and classical transcendental functions," LNCS, vol. 11189, pp. 3-33, 2019. DOI: 10.1007/978-3-030-10692-8_26.

[14] E. A. Ayryan, M. D. Malykh, L. A. Sevastianov, and Yu Ying, "On explicit difference schemes for autonomous systems of differential equations on manifolds," LNCS, vol. 11661, pp. 343-361, 2019. DOI: 10.1007/978-3-030-26831-2_23.

[15] R. Kozlov, "Conservative discretizations of the Kepler motion," Journal of Physics A, vol. 40, no. 17, pp. 4529-4539, 2007. DOI: 10.1088/17518113/40/17/009.

[16] Y. A. Blinkov and V. P. Gerdt, "On computer algebra aided numerical solution of ODE by finite difference method," in Polynomial Computer Algebra '2019. April 15-20, 2019. Euler International Mathematical Institute, St. Petersburg, RUSSIA. 2019.

For citation:

V. P. Gerdt, M. D. Malykh, L. A. Sevastianov, Yu. Ying, On the properties of numerical solutions of dynamical systems obtained using the midpoint method, Discrete and Continuous Models and Applied Computational Science 27 (3) (2019) 242-262. DOI: 10.22363/2658-4670-2019-27-3-242-262.

Information about the authors:

Vladimir P. Gerdt — 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 (e-mail: gerdt@jinr.ru, phone: + 7 (49621) 63437, ORCID: https://orcid.org/0000-0002-0825-1811, ResearcherlD: T-5179-2019, Scopus Author ID: 57209076002)

Mikhail D. Malykh (Russian Federation) — Candidate of Physical and Mathematical Sciences, assistant professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University) (e-mail: malykh-md@rudn.ru, phone: +7(495)9550927, ORCID: https://orcid.org/0000-0001-6541-6603, ResearcherID: P-8123-2016, Scopus Author ID: 6602318510)

Yu Ying (China) — 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 (e-mail: yingy6165@gmail.com, phone: +7(495)9550927, ORCID: https://orcid.org/0000-0002-4105-2566, ResearcherlD: AAC-8344-2019, Scopus Author ID: 57208127921) Leonid A. Sevastianov (Russian Federation) — professor, Doctor of Physical and Mathematical Sciences, professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University) (e-mail: sevastianov-la@rudn.ru, phone: +7(495)955-07-83, ORCID: https://orcid.org/0000-0002-1856-4643, ResearcherID: B-8497-2016, Scopus Author ID: 8783969400)

УДК 517.9

DOI: 10.22363/2658-4670-2019-27-3-242-262

О о о

свойствах численных решении динамических систем, полученных по схеме средней точки

В. П. Гердт1, М. Д. Малых2, Л. А. Севастьянов12, Юй Ин2 3

1 Объединенный институт ядерных исследований ул. Жолио Кюри, д. 6, г. Дубна, Московская область, 141980, Россия 2 Кафедра прикладной информатики и теории вероятностей Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, г. Москва, 117198, Россия 3 Кафедра алгебры и геометрии Университет Каили Kaiyuan Road 3, Каили 556011, Китай

В статье рассматривается схема средней точки как разностная схема для динамической системы вида х = f(x). Эта схема замечательна тем, что в силу теоремы Купера сохраняет все квадратичные интегралы движения, более того, это — простейшая схема из числа симплектических схем Рунге-Кутты, обладающих названным свойством.

Свойства приближённых решений изучены в рамках численных экспериментов с линейным и нелинейным осцилляторами, а также с системой нескольких связанных осцилляторов. Показано, что помимо сохранения всех интегралов движения, приближённые решения наследуют периодичность движения. При этом уделено внимание обсуждению введения понятие периодичности приближённого решения, найденного по разностной схеме.

В случае нелинейного осциллятора выполнение каждого шага требует решения системы нелинейных алгебраических уравнений. Обсуждены вопросы организации вычислений по таким схемам. Дано сравнение с другими схемами, в том числе симметрическими относительно перестановки х и х.

Ключевые слова: консервативные конечно-разностные схемы, динамические системы, Sage, Maple