On a set of tests for numerical methods of integrating differential equations, based on the Calogero system

Full Text

1. Introduction We are now developing a system for integrating ordinary differential equations in the Sage computer algebra system called fdm for sage [1] (https: //github.com/malykhmd/fdm). The problem that arises when developing and implementing numerical methods for integrating dynamic systems is the limited number of test examples, most of them taken from classical monographs [2, 3]. Non-integrability itself is an important property of a dynamic system, although very difficult to formalize. Therefore, tests based on integrable systems are obviously doomed to be somewhat one-sided, which, at least at the present stage, cannot be corrected. The second property of widely used tests is their small dimension. Concerning the many-body problem, for tests the two-body problem and special cases in which the three-body problem, limited or complete, has elementary or at least obviously periodic solutions are chosen. Although many such solutions have been found [4, 5], all of them are rather specific and, in particular, lack singular points - or rather - the collisions of the bodies. A very small number of papers [6] are devoted to the development of difference schemes that inherit periodicity. On the contrary, the question of numerically determining the position of moving singular points of the solution has been well elaborated [7-11]. However, this method has not actually been tested on mechanical problems of many bodies, although in these problems the bodies can approach each other at arbitrarily small distances and cases of false positive singularity tests can be expected. It is equally interesting to explore the question of whether these tests can handle multiple collisions of bodies. On the contrary, the attention has always been focused at the preservation of symplectic structure and integrals of motion. In classical many-body problem the questions of integrals of motion are poor. All of these algebraic integrals, except the energy integral, are linear or quadratic and are preserved by any symplectic Runge-Kutta scheme, and the question of conservation of all integrals is reduced to the question of conservation of energy [12, 13]. Therefore, testing conservative difference schemes requires the development of tests based on Hamiltonian systems that have a large supply of algebraic integrals, and therefore a large dimension. Among the integrable high-dimensional Hamiltonian systems, the most promising for creating this kind of tests is the Calogero system [14]. This system describes the motion of n particles of the same mass in one dimension, their interaction potential being inversely proportional to the square of the distance between the bodies. Its solution is described by algebraic functions of time, and the integrals of motion are rational functions. In fact, this system is the only one among many body problems that can be integrated for any number of bodies, and, moreover, the integration does not require any transcendental functions. For this reason, when developing a set of tests in our fdm system, we gave a prominent place to tests based on the Calogero system. In this paper, we present a set of tools for testing integration methods based on the Calogero problem. This system is very convenient for implementation in computer algebra systems, since it has a purely algebraic properties of solutions and integrals of motion. 2. Tools for specifying the Calogero system The Calogero system describes the motion of n material points of unit mass in one dimension, repelling or attracting each other with a force inversely proportional to the third power of the distance. Let us start numbering the points from zero. Let qi be the position of the i-th point, then ∂U i q¨i = - ∂q , i = 0, . . . , n - 1, (1) M. D. Malykh et al., On a set of tests for numerical methods … 389 where b x2 U = '\" V (qi - qj ), V (x) = . i<j In fdm, the process of specifying the initial problem is separated from the application of the numerical method. To describe the initial problem, a special class Initial_problem is used, described in [1]. To set the initial problem for the Calogero system, the calogero_problem function has been added to fdm, which has 4 optional arguments: o ics is the list of initial values, with the positions of the bodies coming first, and then their momenta, o n is the number of bodies, o T is the final time; we always take t = 0 for the initial time, o b is the value of parameter b, default b = -1. As an illustration, let us take the problem of 5 bodies, which at the initial moment of time occupy the positions qi(0) = i, i = 0, . . . , n - 1. Let us take the velocities close to zero, and change t in the interval t ∈ [0, 0.5]. The time unit in this article is seconds. This problem in our system is set as follows: n=5 ics=list(range(n))+[0.1,0.2,-0.1,0,0] problem_calogero=calogero_problem(ics=ics,n=n, T=0.5, b=-1) It can be solved using standard fdm tools, for example, using the 4th order Runge-Kutta method and plotting the dependence of q0 on t: sol=erk(problem_calogero,N=100) This function is described in [1]. Figure 1 presents the result of the calculations. q0 0.35 0.30 0.25 0.20 0.15 0.10 0.05 t 0.1 0.2 0.3 0.4 0.5 Figure 1. The problem of 5 bodies 390 DCM&ACS. 2023, 31 (4) 387-398 3. Tools for the analytical solution of the Calogero system The exact solution is described by means of the Lax pair [14]. The solution q0(t), . . . , qn-1(t) is a set of eigenvalues of the matrix M = Q|t=0 + tL|t=0, where and Q = diag(q0, . . . qn-1) √ ( 1 - δjk \ L = diag(p0, p2, . . . , pn-1) + j k -b q - q . Therefore, calculating the coordinates of bodies at time t is reduced to calculating the eigenvalues of the matrix M, i.e., they are the zeros of the polynomial of symbolic variable q. F = det(M - qE) (2) For a fixed value of t, these calculations are implemented as a function calogero_q, which has two required arguments ics and n and two optional arguments b and t. Eigenvalues are calculated in the algebraic number field. Unfortunately, the order of the eigenvalues may not coincide with the numbering order of the bodies. Therefore, this function returns the coordinates of the bodies up to some permutation. For example, for the problem of 5 bodies considered above at time t = 0.1, we can find the positions of the bodies as follows: calogero_q(ics,n,b=-1,t=0.1) [0.02174568377617522?, 1.022008478256518?, 1.989567507785786?, 2.998508164438269?, 3.988170165743252?] This can be compared with the values of qi found using the Runge-Kutta method: Q=[symbolic_expression('q'+str(i)) for i in range(n)] [sol.value(q,0.1) for q in Q] [0.0217456837809339, 1.02200847823443, 1.98956750781722, 2.99850816441585, 3.98817016575156] The main disadvantage of this tool is the calculation of the determinant of the matrix M - qE. For large n, this operation becomes very costly. Function calogero_curve returns the polynomial (2) itself. For example, calogero_curve(ics,n,b=-1) -q^5 + 1/5*q^4*t 18089/3600*q^3*t^2 + 23053/36000*q^2*t^3 526651/129600*q*t^4 + 104483/1296000*t^5 + 10*q^4 2*q^3*t M. D. Malykh et al., On a set of tests for numerical methods … 391 + 18101/600*q^2*t^2 5347/2000*q*t^3 + 266021/32400*t^4 35*q^3 + 34/5*q^2*t 196199/3600*q*t^2 + 6829/4000*t^3 + 50*q^2 43/5*q*t + 51523/1800*t^2 24*q + 12/5*t Standard Sage tools allow plotting qi versus time and compare the plots with the results of calculations using a difference scheme. For example, for the considered example of the 5-body problem, the plots of q0 and q1 versus t, obtained numerically (solid line) and analytically (dashed line) can be displayed in one figure as follows: F=calogero_curve(ics,n) sol.plot(t,q0)+sol.plot(t,q1)+\ implicit_plot(F,(t,0,0.7),(q,0,1.2), color='red', linestyle='-', axes_labels=['$t$','$q$'], aspect_ratio=1/2) The result is presented in figure 2. It is clearly seen that the bodies collide at t "' 0.6. The function implicit_plot is used when drawing graphics, which can quickly draw contour maps, so different types of curves in the figure represent different contours rather than different particle trajectories. 3.0 2.5 2.0 q 1.5 1.0 0.5 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 t Figure 2. Problem of 5 bodies 4. Collision of bodies Collisions of bodies occur when the polynomial (2) has multiple roots. Therefore, the first collision of bodies can be calculated as the smallest positive value of t for which the equations ∂F F = 0, = 0 ∂q 392 DCM&ACS. 2023, 31 (4) 387-398 are compatible. Function calogero_solution_crash allows finding it exactly in algebraic numbers. For example, calogero_solution_crash(ics,n) 0.618840603733536? Unless making special efforts, at such a point ∂2F ∂q2 /= 0. From the implicit function theorem it follows that at the collision point the coordinates of the bodies have an algebraic singularity of the order of 1/2. Known the moment of impact, it is possible to observe the behavior of a particular numerical method near the point of impact. For example, let us increase the final time in the example considered to t = 0.7 and make the plots: problem_calogero=calogero_problem(ics=ics,n=n, T=0.7, b=-1) sol=erk(problem_calogero,N=100) pl=sol.plot(t,q0)+sol.plot(t,q1, axes_labels=['$t$','$q$']) pl.show(ymax=1, ymin=0) The result of the calculations is presented in figure 3. q 3.0 2.5 2.0 1.5 1.0 0.5 t 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Figure 3. Problem of 5 bodies. Collision of bodies q1 и q2 In fdm, there is a tool for numerically determining the position and order of the solution singularity using the Alshina method [7-10]. Previously, we tested it on the two-body problem [11]; testing on problems with a large number of bodies indicated the inefficiency of our implementation and forced us to significantly update it. The syntax remains the same. Function eff_order(problem_calogero,q0,N=100) M. D. Malykh et al., On a set of tests for numerical methods … 393 returns the plot 4. The dotted line in the figure represents the theoretical order, which does not change with time. order 2 1 0.2 0.4 0.6 0.8 1.0 1 r= 0.485661585050870 Figure 4. Problem of 5 bodies, effective order according to Alshina 2 Let us recall how to read this plot [11]. The CROS scheme on which Alshina’s method is based is a second-order scheme. Up to the point of collision, the scheme order, calculated approximately, coincides with the theoretical one, that is, equals 2. Near the collision point, this rule ceases to apply and, which is a specific property of the CROS circuit, the order jumps sharply to a new constant value, which is equal to the order of this special point [7]. In figure 4 we actually see that at the point of impact the order changes sharply from a value of 2 to some value close to 1 . However, in addition to this jump, there are two more ’spikes’, at the very beginning of the graph and in the region of t = 0.8. Near these two points, the order changes sharply, but eventually returns to its original value. The theory does not explain the appearance of these artifacts, which require further research. We plan to implement and test other numerical methods in our system for identifying singularities [15, 16]. 5. Preservation of integrals of motion As noted above, the question of preserving all integrals of motion by a difference scheme has been fruitfully discussed for a long time. The Calogero system has n rational integrals of motion being in involution and, therefore, is completely Liouville integrable [14]. These integrals can be described as traces of the matrix L powers: Fk = Sp Lk , (k = 1, . . . n). 394 DCM&ACS. 2023, 31 (4) 387-398 They are, of course, symmetric functions with respect to the group of permutations of bodies. The calculation of these integrals is implemented as the function calogero_integral(k,n). For example, calogero_integral(1,n) p0 + p1 + p2 + p3 + p4 It should be noted that for b > 0 the matrix L is complex, and the traces of its degree are real. Our function performs the appropriate simplification and returns a rational function with real coefficients. The integral F1 is linear, so all Runge-Kutta methods preserve it exactly. The energy integral F2 is not quadratic for the Calogero system, so even symplectic Runge-Kutta schemes do not preserve it [3]. Figure 5 shows the dependences of F2 and F3 for our 5-body problem; it is clearly seen that these functions monotonically increase in absolute value on approaching the collision point. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 1e1 t 0.0070 0.0075 0.0080 0.0085 0.0090 F2 0.0 0.1 0.2 0.3 0.4 0.5 0.6 t 1.086 1.087 1.088 1.089 1.090 1.091 F3 Figure 5. 5-body problem, integrals F2 and F3 in the explicit Runge-Kutta scheme M. D. Malykh et al., On a set of tests for numerical methods … 395 It should also be noted that sharp jumps in the integrals, which are often observed when applying explicit schemes to the classical, three-dimensional many-body problem, did not appear in these plots. This suggests that the Calogero problem, whose solutions are algebraic functions t, is much simpler than the classical many-body problem and, therefore, the test based on it does not detect this feature of explicit schemes. 6. Conclusion The completely integrable dynamical Calogero system, which describes the one-dimensional many-body problem, allows creating a convenient tool for testing difference schemes by means that do not go beyond the algebraic framework. We have implemented these tools in the new version of the fdm for sage software. The first application showed that the implementation of the method for numerical detection of moving singularities is complicated by artifacts, the false spikes in the order plot, not previously described in theory. This allows an idea that the developed tools provide a fairly representative set of tests that will allow revealing previously unnoticed difficulties. On the other hand, it should be noted that the trajectories of bodies in this dynamic system are arranged quite simply, without any ’loops’ and ’fine structures’, which makes this test rather rough. In particular, we associate this with the absence of jumps of integrals of motion in the plots calculated using explicit schemes. This offers a prospect for the development of new tests.

About the authors

Mikhail D. Malykh

RUDN University; Joint Institute for Nuclear Research

Author for correspondence.
Email: malykh_md@rudn.ru
ORCID iD: 0000-0001-6541-6603
Scopus Author ID: 6602318510
ResearcherId: P-8123-2016

Doctor of Physical and Mathematical Sciences, Assistant Professor

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

Wang Shiwei

RUDN University

Email: 1995wsw@gmail.com
ORCID iD: 0009-0007-6504-8370

Ph.D. student

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

Yu Ying

Kaili University

Email: 45384377@qq.com
ORCID iD: 0000-0002-4105-2566

Assistant Professor of Department of Mathematics and Applied Mathematics

3 Kaiyuan Rd., Kaili, 556011, People’s Republic of China

References