CC BY
3
Discrete & Continuous Models & Applied Computational Science 2022, 30 (4) 364-373
ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph
Research article
UDC 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2022-30-4-364-373
Conservative finite difference schemes for dynamical
systems
Yu Ying, Zhen Lu
Kaili University, Kaiyuan Road 3, Kaili, 556011, China
(received: November 6, 2022; revised: December 8, 2022; accepted: December 19, 2022)
Abstract. The article presents the implementation of one of the approaches to the integration of dynamical systems, which preserves algebraic integrals in the original fdm for Sage system . This approach, which goes back to the paper by del Buono and Mastroserio, makes it possible, based on any two explicit difference schemes, including any two explicit Runge-Kutta schemes, to construct a new numerical algorithm for integrating a dynamical system that preserves the given integral. This approach has been implemented and tested in the original fdm for Sage system. Details and implementation difficulties are discussed. For testing, two Runge-Kutta schemes were taken having the same order, but different Butcher tables, which does not complicate the method due to paralleling. Two examples are considered — a linear oscillator and a Jacobi oscillator with two quadratic integrals. The second example shows that the preservation of one integral of motion does not lead to the conservation of the other. Moreover, this method allows us to propose a practical application of the well-known ambiguity in the definition of Butcher tables.
Key words and phrases: finite difference method, dynamical systems, explicit Runge-Kutta methods
1. Introduction
Many dynamical systems have algebraic integrals of motion [1], but standard numerical methods do not allow preserving these integrals exactly on the approximate solution [2]. This means that the approximate solution satisfies such fundamental laws of nature as the law of conservation of energy also approximately, and, in view of the importance of this law itself, this circumstance is always striking.
Consider the dynamical system
— = /(x), xe\Rm, (1)
© YingY., LuZ., 2022
This work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by-nc/4-0/legalcode
whose right side is a rational function with rational coefficients. Since difference schemes are described by algebraic equations, there are no obstacles to the case that among these schemes there be those that preserve all the algebraic integrals of this system. Linear integrals of motion are preserved by almost all schemes used. In the late 1980s, among the Runge-Kutta schemes, a subclass of schemes was discovered that preserves not only linear, but also quadratic integrals of motion. This class is called Runge-Kutta symplectic schemes [2]. These schemes make it possible, for example, to preserve all algebraic integrals in the top rotation problem, with the exception of the Ko-valevskaya case [3]. In the 1990s, Greenspan constructed the first difference scheme for the many-body problem that preserves all algebraic integrals of this problem [4-7], using the principle of energy quadratization; such schemes can be constructed for the many-body problem based on any Runge-Kutta symplectic scheme [8, 9].
The main disadvantage of symplectic Runge-Kutta schemes is their implicit nature: in calculations using these schemes, a system of nonlinear equations has to be solved at each step, which significantly complicates the calculations compared to the commonly used explicit schemes. Unfortunately, in any case, there are no such schemes among Runge-Kutta schemes, which became clear at the dawn of the theory of symplectic Runge-Kutta schemes [2].
In the paper by del Buono and Mastroserio [10] an approach to constructing conservative difference schemes was proposed, which can be described as follows. Let the dynamical system (1) have an integral g and let there be two difference schemes. Let the step of the first scheme be described as x = fi(x,dt), and the step of the second scheme as x = ^(x,dt). Consider the composite scheme
x = 4>(x, dt) + ^(x, dt).
We use the parameter ^ in such a way that the given integral g be preserved on the approximate solution found exactly by the composite scheme.
To do this, we will describe the transition from n to n+1 as follows: xn is given, and the next value xn+1 is found from the system of equations
xn+1 = &(xn, dt) + vn+Mxn, dt), g(xn+i) = 9(xn), whose solution reduces to solving one equation
g(^(xn, dt) + ^n+i4)(xn, dt)) = g(xn)
with respect to ^n+1.
Such an approach, of course, does not avoid the main difficulty that arises when using implicit difference schemes: at each step, you still have to solve a nonlinear algebraic equation. But this is only one equation and its degree coincides with the degree of the integral g. Therefore, the pioneers called it the explicit conservative Runge-Kutta method.
In fact, this method certainly does not belong to the Runge-Kutta family of methods. In terms of difference schemes, it can be written as follows. One more variable ^ is added to the variables x and the transition from (x,^) to
(x, p) is described by a system of algebraic equations
x = 4>(x, dt) + dt), g(x) = g(x).
It should be noted that this system is not a difference scheme for any dynamical system in the variables (x,^). Here we are dealing with a certain generalization of the very concept of a difference scheme.
In Ref. [11] it was proposed to use as 0 the Runge-Kutta scheme nested in the 0 scheme and having a smaller order. Such an approach was fully justified then, since it promised to reduce the cost of calculating 0. However, now that multi-core processors have come into general use, it seems superfluous to impose such a restriction: it is easier to calculate 4>(xn,dt) and ^(xn,dt) at each step in parallel.
In fact, the described approach, which we will call del Buono and Mastroserio approach, makes it possible on the basis of any two explicit difference schemes, including any two explicit Runge-Kutta difference schemes, to construct a new numerical algorithm for integrating a dynamical system that preserves the integral g. There is some arbitrariness [12] when calculating explicit Butcher tables, which give the highest possible order of approximation for a fixed number of stages. Therefore, two Runge-Kutta schemes of the same order, but with different Butcher tables, can be taken as two schemes.
We have implemented del Buono and Mastroserio's approach in the Sage computer algebra system, and in a series of computer experiments we have seen that it does indeed preserve the integrals of motion [13]. Since then, a number of our developments related to the solution of ordinary differential equations using the finite difference method was compiled into the fdm for Sage package presented at ITTMM'2022 [14]. In this article, we want to discuss the details of the implementation of del Buono and Mastroserio's approach in our package. Since it has built-in tools for working with Butcher tables and estimating the approximation error, we will further evaluate the possibility of using two Runge-Kutta schemes of the same order, but with different Butcher tables. The implementation itself is available at https://github.com/malykhmd/fdm.
2. Implementation of del Buono and Mastroserio's
approach
Let a dynamical system (1) be given, its integral of motion g, which will be preserved in the numerical solution, and two explicit Runge-Kutta difference schemes x = fi(x,dt) and x = ^(x,dt), which are given using two Butcher tables. Our implementation of the del Buono and Mastroserio approach in fdm for Sage supports any pair of Butcher tables.
With two stages, the highest order of approximation that Runge-Kutta schemes can have is 2. In this case, one parameter of the Butcher table remains undefined. Therefore, in our experiments, schemes with Butcher
tables were taken as two initial schemes
and
1 1
2 2
3 1
4 4
To calculate an approximate solution of the Cauchy problem
x = f^^ x\t=0 = x0
on the segment 0 < t < T, we divide this segment into N equal parts with the step dt = T/N. Then, in a loop over n = 1,2,...,N, we calculate X1,X2, ..•. At each step in n, we first calculate the values x'n+1 = <fi(xn, dt) and x"i+1 = ^(xn, dt) by these schemes in parallel, then we calculate the root of the equation
9(xn+1 + tlxn+1) = 9(xn) (2)
relative to Then we find
Xn+1 Xn+1 + l1Xn+1.
In this case, a certain arbitrariness arises in the choice of the method for solving the algebraic equation (2). When applying implicit Runge-Kutta methods, iterative methods are used, since there is a natural approximation to the desired root [7]. In this case, however, we do not know a good approximation to ^, and the equation itself has a small degree and, at least in the examples considered below, admits a solution in radicals. We tried different methods for solving this equation, in particular, explicitly expresses the solution in radicals, but this did not give any noticeable increase in comparison with the standard method of finding the roots of polynomials \R[v], implemented in Sage.
3. Test examples
As a first test case, we took a linear oscillator
d d
—x = y, —
dt dt
x(0) = 0, y(0) = 1
dtx У, dty x, (3)
in the segment 0 < t < 10.
In figure 1, it is clearly seen that the quadratic integral of motion x2 + y2 is preserved on the approximate solution found by the composite scheme, which favorably distinguishes this scheme from the usual explicit Runge-Kutta scheme. The slope of the Richardson diagram (figure 2) [15] shows that the order of approximation of the composite circuit is equal to the order of the original ones, that is, 2, as expected from theoretical considerations. It should also be noted that the value of ^ did not change from step to step.
x2+y2
1.0025 H
1.0020
1.0015
1.0010
1.0005
1.0000 -K
10
Figure 1. Variation of integral x2 + y2 on the approximate solutions of the problem (4): the explicit Runge-Kutta scheme (dotted line) and the composite scheme (solid line)
\E(x 10° -
101
10-2
10"3
10^ 1
— 2/=2.00 x +0.130
10
101
10°
Figure 2. Richardson diagram for the value of x at t = 9 in the example (4)
As a second test case, we took the Jacobi oscillator:
d d d -p = qr, —q = —pr,
dt dt j>(0) = 0, <7(0) = 1, r(0) = l
dtr = ~ VKL
(4)
in the segment 0 < t < 10. This system admits two independent quadratic integrals, we took p2 + q2 as g.
In figure 3, it is clearly seen that this integral of motion p2 + q2 is preserved on the approximate solution found by the composite scheme. However, the second integral p2/A + r2 is not preserved in this case (see figure 4). According
to the Richardson diagram (figure 5), it can be seen that in this case, too, the order of approximation of the composite circuit is equal to that of the original ones. Unlike the linear case, the value of ^ changed from step to step, and very smoothly (figure 6).
p2 + q2
1.0025 1.0020 1.0015 1.0010 1.0005
Figure 3. Variation of the integral x2 + y2 on approximate solutions of the problem (3): explicit Runge-Kutta scheme (dotted line) and composite scheme
\p2+r2
10
Figure 4. Variation of the integral x2 + y2 on approximate solutions of the problem (3): explicit Runge-Kutta scheme (dotted line) and composite scheme (solid line)
t
o
2
6
S
10
t
0
2
6
8
4. Conclusion
The performed experiments confirm that the approach, which goes back to the article by del Buono and Mastroserio [10], can also be used with a non-standard choice of initial schemes. Moreover, this method allows us to offer a practical application of the well-known ambiguity in the definition of Butcher tables.
№)l
Figure 5. Richardson diagram for value x at t = 9 for example (3)
0 2 4 6 8 10
Figure 6. Dependence of parameter ¡i on time t for example (3)
However, we should admit that the main problem of all known conservative schemes, the necessity to solve nonlinear equations at each step, has not been completely eliminated. The gain that the described method gives can be described as follows. If only one integral of motion of a dynamical system with n unknowns needs to be preserved, the Runge-Kutta symplectic method leads to the solution of a system of n algebraic equations at each step, and the method under discussion only one equation with one unknown. It is well known that numerical methods for solving one equation with one unknown are much simpler and more reliable than methods for solving systems [16]. The smoothness of changing the values of the parameter fi noted above suggests that iterative algorithms, for example, Newton's method, can be used to calculate it.
Unfortunately, as was seen in the second test case, the inheritance of one integral of motion does not entail the inheritance of the second. The del Buono and Mastroserio approach has already been generalized to the case of several integrals of motion [11] and in this case leads to the solution of
a system of nonlinear equations at each step with respect to several auxiliary
parameters. Thus, we return to the difficulty typical of the symplectic Runge-Kutta methods. The study of the gain that the del Buono and Mastroserio approach gives, when applied to systems with several algebraic integrals of motion, in comparison with classical implicit methods requires further study.
It is no less interesting how the del Buono and Mastroserio approach will show itself in solving the many-body problem, in which the approach of bodies significantly worsens the convergence of iterative methods used to solve nonlinear systems that arise when using the implicit Runge-Kutta methods.
References
[1] A. Goriely, Integrability and Nonintegrability of Dynamical Systems. Singapore; River Edge, NJ: World Scientific, 2001.
[2] E. Hairer, G. Wanner, and S. P. N0rsett, Solving Ordinary Differential Equations I, Nonstiff Problems, 3rd ed. Springer, 2008. DOI: 10.1007/ 978-3-540-78862-1.
[3] V. V. Golubev, Vorlesungen über Differentialgleichungen im Komplexen. VEB Deutscher Verlag der Wissenschaften, 1958.
[4] D. Greenspan. "Completely Conservative and Covariant Numerical Methodology for N-Body Problems With Distance-Dependent Potentials. Technical Report no. 285." (1992), [Online]. Available: http://hdl. handle.net/10106/2267.
[5] 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.
[6] 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.
[7] D. Greenspan, N-Body Problems and Models. World Scientific, 2004.
[8] Y. Ying, A. Baddour, V. P. Gerdt, M. Malykh, and L. Sevastianov, "On the quadratization of the integrals for the many-body problem," Mathematics, vol. 9, no. 24, 2021. DOI: 10.3390/math9243208.
[9] A. Baddour and M. Malykh, "On difference schemes for the many-body problem preserving all algebraic integrals," Phys. Part. Nuclei Lett., vol. 19, pp. 77-80, 2022. DOI: 10.1134/S1547477122010022.
[10] N. Del Buono and C. Mastroserio, "Explicit methods based on a class of four stage fourth order Runge-Kutta methods for preserving quadratic laws," Journal of Computational and Applied Mathematics, vol. 140, pp. 231-243, 2002.
[11] M. Calvo, D. Hernandez-Abreu, J. I. Montijano, and L. Randez, "On the preservation of invariants by explicit Runge-Kutta methods," SIAM Journal on Scientific Computing, vol. 28, no. 3, pp. 868-885, 2006.
[12] Y. Ying, "The symbolic problems associated with Runge-Kutta methods and their solving in Sage," Discrete and Continuous Models and Applied Computational Science, vol. 27, no. 1, pp. 33-41, 2019. DOI: 10.22363/ 2658-4670-2019-27-1-33-41.
[13] Y. Ying and M. Malykh, "On the realization of explicit Runge-Kutta schemes preserving quadratic invariants of dynamical systems," Discrete and Continuous Models and Applied Computational Science, vol. 28, no. 4, pp. 313-331, 2020. DOI: 10.22363/2658-4670-2020-28-4-313-331.
[14] L. González and M. D. Malykh, "On a new package for the numerical solution of ordinary differential equations in Sage," in Information and telecommunication technologies and mathematical modeling of high-tech systems. Materials of the All-Russian Conference with international participation, In Russian, Moscow: RUDN, 2022.
[15] A. Baddour and M. Malykh, "Richardson-Kalitkin method in abstract description," Discrete and Continuous Models and Applied Computational Science, vol. 29, no. 3, pp. 271-284, 2021.
[16] W. H. Press. "Numerical Recipes Home Page." (2019), [Online]. Available: http://numerical.recipes.
For citation:
Y. Ying, Z. Lu, Conservative finite difference schemes for dynamical systems, Discrete and Continuous Models and Applied Computational Science 30 (4) (2022) 364-373. DOI: 10.22363/2658-4670-2022-30-4-364-373.
Information about the authors:
Ying, Yu — Assistant Professor of Department of Algebra and Geometry, Kaili University, China (e-mail: 45384377@qq.com, ORCID: https://orcid.org/0000-0002-4105-2566)
Lu, Zhen — Associate Professor, Department of Fine art, Kaili University, China (e-mail: 157739594@qq.com, ORCID: https://orcid.org/0000-0002-7526-9026)
УДК 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2022-30-4-364-373
Консервативные конечно-разностные схемы для динамических систем
Юй Ин, Чжэнь Лу
Университет Кайли, 3, Кайюань Роуд, Кайли, 556011, Китай
Аннотация. В статье представлена реализация одного из подходов к интегрированию динамических систем, при котором сохраняются алгебраические интегралы в оригинальной системе fdm for sage. Этот подход, восходящий к статье дель Буоно и Мастросерио, позволяет на основе двух любых явных разностных схем, в том числе любых двух явных схем Рунге—Кутты, сконструировать новый численный алгоритм интегрирования динамической системы, сохраняющий заданный интеграл. Этот подход реализован и протестирован в оригинальной системе fdm for sage. Обсуждены детали и трудности реализации. Для тестирования в качестве двух схем взяты две схемы Рунге—Кутты одного порядка, но с разными таблицами Бутчера, что не приводит к усложнению метода благодаря распараллеливанию. Рассмотрено два примера — линейный осциллятор и осциллятор Якоби, имеющий два квадратичных интеграла. На втором примере показано, что сохранение одного интеграла движения не приводит к сохранению другого. Проделанные эксперименты подтверждают, что данный подход может быть использован и при нестандартном выборе исходных схем. Более того, этот метод позволяет предложить практическое применение хорошо известной неоднозначности в определении таблиц Бутчера.
Ключевые слова: метод конечных разностей, динамические системы, явные методы Рунге-Кутты
