CC BY
6
Research article
UDC 519.872, 519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2021-29-1-63-72
On conjugate difference schemes: the midpoint scheme and the trapezoidal scheme
Yu Ying1, Mikhail D. Malykh2
1 Kaili University 3, Kaiyuan Road, Kaili, 556011, China 2 Peoples' Friendship University of Russia (RUDN University) 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation
(received: February 10, 2021; accepted: March 12, 2021)
The preservation of quadratic integrals on approximate solutions of autonomous systems of ordinary differential equations x = f(x), found by the trapezoidal scheme, is investigated. For this purpose, a relation has been established between the trapezoidal scheme and the midpoint scheme, which preserves all quadratic integrals of motion by virtue of Cooper's theorem. This relation allows considering the trapezoidal scheme as dual to the midpoint scheme and to find a dual analogue for Cooper's theorem by analogy with the duality principle in projective geometry. It is proved that on the approximate solution found by the trapezoidal scheme, not the quadratic integral itself is preserved, but a more complicated expression, which turns into an integral in the limit as At ^ 0. Thus the concept of conjugate difference schemes is investigated in pure algebraic way. The results are illustrated by examples of linear and elliptic oscillators. In both cases, expressions preserved by the trapezoidal scheme are presented explicitly.
Key words and phrases: dynamical systems, quadratic integrals, difference schemes, conservation laws, midpoint scheme, trapezoidal scheme
1. Introduction
Dynamical systems are the most important mathematical models in mechanics and physics. Only a few of these models are integrated in a closed form [1], therefore, they have to be investigated using numerical methods, of which the most important is the finite difference method.
Let x be a point in an m-dimensional affine space. Any difference scheme that approximates differential equation
Ah**) w
© YingY., Malykh M.D., 2021
/T\ This work is licensed under a Creative Commons Attribution 4.0 International License
®
— http://creativecommons.org/licenses/by/4-0/
describes a transition from the value x at some initial moment of time to the value of x at the moment of time shifted from the initial value by the quantity At, called the step. We will consider algebraic schemes, i.e., those in which the above correspondence is specified using a system of algebraic equations
F(x,x,At) = 0. (2)
If the original equation has an algebraic integral g(x) = C, and it follows from the equations (2) that
g(x) = g(x),
then this difference scheme is said to preserve this integral.
If we use explicit difference schemes for integrating dynamical systems, then the values of the integrals of motion will change monotonically step by step. At the turn of the 1980s and 1990s, the first difference schemes were constructed that preserve exactly the algebraic integrals of dynamical systems. For example, the scheme constructed by D. Greenspan preserves all classical integrals of N-body problems [2]-[5], the symplectic Runge-Kutta schemes, including the simplest of them, the midpoint scheme
£-x = f(^-)At, (3)
preserve linear and quadratic integrals in virtue Cooper's theorem [6]-[9]. This circuit has a whole bunch of wonderful properties inherited from the original differential equation [10].
This seems to be a simple consequence of the ¿-symmetry of the midpoint circuit: the equation (3) is invariant under the transformation
At — -At, x —y x, x —y x.
The trapezoidal scheme has the same property
£-x = (f(£) + f(x))Af, (4)
however, in experiments with an elliptic oscillator performed by Yu. A. Blinkov for PCA'2019 [11], the quadratic integrals oscillated, although they did not increase monotonically. The absence of monotonicity in the variation of the values of the integrals of motion on approximate solutions is extremely important from the physical point of view, since, on average, all fundamental conservation laws are satisfied on solutions of this type.
The noted behavior of the approximate solutions found by the trapezoidal scheme can be explained by the fact that it is conjugated to the midpoint scheme and therefore some more complex expression is retained on it [9, §VI.8.1-2].
The very concept of conjugate difference schemes [9, def 8.1] is formulated locally in terms of power series. The implicit function theorem can be applied to the system of algebraic equations (2) and, under certain conditions, we can assert that
x = x + f(x)At + ... = $At (x),
where $ is a series in powers of At, the coefficients of which are rational functions of x. Difference schemes
x = $At(x) and x = ^At(x) are referred to as mutually conjugate, if there exists a change of coordinates
x = XAt(y) = V +
such that
$ At = Xa1 0 ^ At ° XAt.
It is clear from this definition that the exact preservation of the expression of one of these schemes entails the preservation of some expression by the other scheme. Say, if the scheme ^At preserves the integral g(x) exactly, then the scheme $At preserves the expression
9(XAt(x)) = 9(x) + 9i(x)At + . ,
depending on At [9, §VI.8.2].
In this article, we will clarify the geometric meaning of the conjugacy of the two above schemes and write down explicitly the expression that preserves the trapezoidal scheme.
2. Relationship between trapezoidal and midpoint
schemes
The approximate solution of the system (1), found by the scheme (2) with constant step At, is a finite or infinite sequence of points
x0 ,xl,x2, . (5)
the first element of which is chosen in an arbitrary way, and all the others are defined recursively: xn+1 is the root x of the equation
F(xn, x, At) = 0,
tending to xn at At ^ 0.
Theorem 1. Let x0,x1 ,x2,..■ be an approximate solutions of equation (1), calculated using the midpoint scheme (3). Then coordinates x'0,x1 ,x'2,... of middles of links of broken line x0x1 x2 ... (figure 1) yield another approximate solution of the same equation, calculated by the trapezoidal scheme (4).
Proof. The middle of link xnxn+1 of the solution found by the midpoint scheme is given by the formula
r _ xn+1 + xn
Xn O .
1
•En— 1
Figure 1. The solutions and ... found by the midpoint schemes
and by the trapezoidal scheme
In this case, the ends of the link are unambiguously reconstructed from its given midpoint:
_ xn+1 + xn + xn+1 xn _ ' + f( '
xn+i 2 2 xn ' J\xn) 2
and
_ _ xn+1 + Xn _ Xn+1 _ Xn _ / _ r// ) A
xn 2 2 Ln ^ (xn) 2 '
Since xn+1 belongs to two links xnxn+1 and xn+1 xn+2, we have
xn+1 xn + f(xn) 2 Xn+1 f(xn+l) 2
from where it immediately follows that
At
x
n+1 Xn = (f(xn+l) + f(xn)) 0 •
Thus, the midpoints can be calculated using the trapezoidal difference scheme (4). □
By virtue of the theorem 1, the solutions x0x1 x2 ... and X q X i X ' .... found by the midpoint schemes (3) and by the trapezoidal scheme (4) turn out to be coupled with each other. By analogy with the duality principle in projective geometry [12] it is hoped that any statement about the midpoint scheme should have a 'twin' in the trapezoidal scheme.
One of the most interesting properties of the midpoint scheme is Cooper's theorem [9, th. 2.2], according to which this scheme preserves any quadratic integral of motion.
Equality
g(xn) = g(xn+1)
is easily rewritten by expressing xn though x'n, and xn+1 though x'n+1:
9 (xn — f(xn)_2_) = ^ {xn+1 — f(xn+1 .
Therefore, for dual scheme (4) the conservation law takes the form
g(x-f(x)A) = g(x-f(x)A). (6)
Thus, a quadratic integral is also inherited by the trapezoidal scheme (4), but the expression for the conserved quantity coincides with g only in the limit At ^ 0. This circumstance made complicated finding it.
Definition 1. Let us say that a difference scheme inherits the integral g(x) = C if there exists a rational function G(x, At) such that
1) from the equations that specify the scheme it follows that
G(x,At) = G(x,At),
2) in the limit At ^ 0 expression G(x, At) turns into g(x)
Function G itself will be referred to as the difference analog of the integral g.
Theorem 2 (Cooper's dual theorem). The trapezoidal scheme inherits all linear and quadratic integrals of motion, and the difference analogue of the integral g will be
9(x- f(x)A).
3. Examples
Consider several examples.
3.1. Linear oscillator
In the case of a linear dynamical system, the midpoint scheme and the trapezoidal scheme are the same, so the midpoint scheme becomes self-conjugate. This circumstance greatly simplifies the study of the midpoint scheme for a linear oscillator.
Consider a dynamic system
x = -y, y = x (7)
which has a quadratic integral
x2 + y2 = C.
The midpoint scheme will give points (x0, y0),(x1 ,y1),... lying on the circle
X2 + y2 _ R2, the radius of which is determined by the initial point
R _ ^x2 + yQ'
The midpoints of the links lying on the circle
x2 + y2 _ r2,
the radius of which can be determined by the first link
r _ J(xi+xq )2 + (vi+y0Y _ R +(—>
Thus, the trajectory on the phase plane turns out to be a broken line, the vertices of which lie on a circle of radius R, and the links touch a concentric
circle, the radius of which is Jl +At2/4 times less than R. In particular,
the trajectory will be closed, and the solution will be periodic if R and r are the radii of the circumscribed and inscribed circle in the A-gon, that is, if
w
r/R _ cos —'
This immediately gives the formula for choosing a step
^l + At2/4 _ cos ^' This formula was previously obtained by us analytically [10].
3.2. Elliptic oscillator
By the definition of Jacobi functions [13],
p_ sn t, q_ cn t, r_ dn t is a particular solution of the autonomous system of differential equations
p _ qr, q _ _pr, r _ _k2pq (8)
with the initial conditions
p _ 0, q _ r _ 1 at t _ 0' The midpoint scheme preserves both intergals
p2 + q2 _ const and k2p2 + r2 _ const (9)
of this system. Now the trapezoidal scheme
,At
p-p=(qr + qr) —,...
does not coincide with the midpoint scheme and, therefore, its invariants are more complicated.
Nevertheless, the integral
p2 + q2 = const
corresponds to the integral
At\2 ( At\ 2 2 2 2 2 At2 2 2 At2 p — qr— \ + \q + pr—\ = p2 + q2 + q2 r2 ——+ p2
or
(p2 + q2 ) (1 + r~AA^
The integral correspods to
k2 p2 + r2 = const
( At\ ( At\2 At2 At2
k2 ( p — qr—J + (r + k2Pq^J = k2p2 + r2 + k2q2r2——+ k4p2q:
or
2
At2
(k2p2 + r2 ) (1 + k2 q2 —
Thus, in the space pqr the vertices of the trajectory lie on the elliptic curve (9), and the midpoints of the links of the broken line lie on a more complex curve
(p2 + q2)\l + r~^ )=Ci, (k2p2 + r2)(l + k2q2 ^ \=C2. (10)
This means that the trapezoidal scheme for an elliptic oscillator inherits both quadratic integrals, and their difference counterparts are the expressions (10).
If we follow the change in p2 + q2 on the approximate solution found by the trapezoidal scheme, then we will see a deviation from a constant value equal to
. 2 2. r2At2 (P2 + Q2
The exact solution is periodic, so in the plots these deviations appear as periodic fluctuations.
4. Conclusion
If you do not use specially developed difference schemes, discretization of continuous models by the method of finite differences introduces completely new properties into these models: calculations lead to a monotonic change in quantities, which, from physical considerations, must remain constant. For example, in computer experiments, dissipation appears even in those cases when energy was conserved in the original continuous model. In calculations for sufficiently long time intervals, this dissipation becomes very noticeable, and the parameters of the dynamical system are significantly distorted.
The number of schemes that preserve algebraic integrals of motion exactly is small and their drawbacks are well known. Difference schemes, in which the integrals of motion fluctuate around their initial values, significantly expand this set. However, the noted property is usually accepted without explanation and even more rigorous proof. Theorem 2, which is dual to Cooper's theorem, allows us to fill in this gap for the trapezoidal scheme (4) by explicitly specifying expressions that coincide in the limit At ^ 0 with exact integrals and at the same time are preserved on approximate solutions exactly.
It would be very interesting to generalize this result to other schemes, the use of which does not lead to a monotonic increment of the integrals of motion. For such a generalization, in our opinion, it is necessary to investigate in more detail the question of schemes that are, in a sense, dual to the symplectic Runge-Kutta schemes.
Acknowledgments
The publication was supported by the RUDN University Strategic Academic Leadership Program.
References
[1] A. Goriely, "Integrability and nonintegrability of dynamical systems," in Advanced Series in Nonlinear Dynamics. Singapore; River Edge, NJ: World Scientific, 2001, vol. 19. doi: 10.1142/3846.
[2] D. Greenspan, "Completely conservative, covariant numerical methodology," Computers & Mathematics with Applications, vol. 29, no. 4, pp. 3743, 1995. DOI: 10.1016/0898-1221(94)00236-E.
[3] D. Greenspan, "Completely conservative, covariant numerical solution of systems of ordinary differential equations with applications," Rendiconti del Seminario Matematico e Fisico di Milano, vol. 65, pp. 63-87, 1995. DOI: 10.1007/BF02925253.
[4] J. C. Simo and M. A. González, "Assessment of Energy-momentum and Symplectic Schemes for Stiff Dynamical Systems," in American Society of Mechanical Engineers. ASME Winter Annual Meeting, New Orleans, Louisiana, 1993.
[5] E. Graham, G. Jelenic, and M. A. Crisfield, "A note on the equivalence of two recent time-integration schemes for N-body problems," Communications in Numerical Methods in Engineering, vol. 18, pp. 615-620, 2002. DOI: 10.1002/cnm.520.
[6] G. J. Cooper, "Stability of Runge-Kutta methods for trajectory problems," IMA Journal of Numerical Analysis, vol. 7, pp. 1-13, 1 1987. DOI: 10.1093/imanum/7.1.1.
[7] Y. B. Suris, "Hamiltonian methods of Runge-Kutta type and their variational interpretation [Gamil'tonovy metody tipa Runge-Kutty i ikh variatsionnaya traktovka]," Matematicheskoe modelirovaniye, vol. 2, no. 4, pp. 78-87, 1990, in Russian.
[8] J. M. Sanz-Serna, "Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more," SlAM review, vol. 58, pp. 3-33, 2016. DOI: 10.1137/151002769.
[9] E. Hairer, G. Wanner, and C. Lubich, Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations. Berlin Heidelberg New York: Springer, 2000.
[10] V. P. Gerdt, M. D. Malykh, L. A. Sevastianov, and Yu Ying, "On the properties of numerical solutions of dynamical systems obtained using the midpoint method," Discrete & Continuous Models & Applied Computational Science, vol. 27, no. 3, pp. 242-262, 2019. DOI: 10.22363/ 2658-4670-2019-27-3-242-262.
[11] Y. A. Blinkov and V. P. Gerdt, "On computer algebra aided numerical solution of ODE by finite difference method," in International Conference Polynomial Computer Algebra'2019; St. Petersburg, April 15-20, 2019, N. N. Vassiliev, Ed., SPb: VVM Publishing, 2019, pp. 29-31.
[12] F. Klein, Vorlesungen über Nicht-Euklidische Geometrie. Springer, 1967. DOI: 10.1007/978-3-642-95026-1.
[13] P. F. Byrd and M. D. Friedman, Handbook of Elliptic Integrals for Engineers and Scientists. Springer, 1971. DOI: 10.1007/978-3-64265138-0.
For citation:
Y. Ying, M. D. Malykh, On conjugate difference schemes: the midpoint scheme and the trapezoidal scheme, Discrete and Continuous Models and Applied Computational Science 29 (1) (2021) 63-72. DOI: 10.22363/2658-4670-202129-1-63-72.
Information about the authors:
Ying, Yu — Candidate of Physical and Mathematical Sciences, Assistant Professor of Department of Algebra and Geometry, Kaili University (e-mail: 45384377@qq.com, ORCID: https://orcid.org/0000-0002-4105-2566) Malykh, Mikhail D. — Doctor 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@pfur.ru, phone: +7(495)9550927, ORCID: https://orcid.org/0000-0001-6541-6603, ResearcherlD: P-8123-2016, Scopus Author ID: 6602318510)
УДК 519.872, 519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2021-29-1-63-72
О** о
сопряженных разностных схемах: схема средней точки и схема трапеций
Юй Ин1, М. Д. Малых2
1 Университет Кайли Kaiyuan Road 3, Кайли, 556011, Китай 2 Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия
В статье исследован вопрос о сохранении квадратичных интегралов на приближённых решениях автономных систем обыкновенных дифференциальных уравнений х = f(x), найденных по схеме трапеций. Установлена связь между схемой трапеции и схемой средней точки, которая сохраняет все квадратичные интегралы движения в силу теоремы Купера. Эта связь позволяет рассматривать схему трапеций как двойственную к схеме средней точки и отыскать двойственный аналог для теоремы Купера. Доказано, что на приближённом решении, найденном по симметрической схеме, сохраняется не сам квадратичный интеграл, а более сложное выражение, которое переходит в интеграл в пределе при At ^ 0. Результаты проиллюстрированы примерами — линейным и эллиптическим осцилляторами. В обоих случаях в явном виде выписаны выражения, которые сохраняет схема трапеций.
Ключевые слова: динамические системы, квадратичные интегралы, разностные схемы, законы сохранения, схема средней точки, схема трапеций
