Вычисление модифицированного гамильтониана в Sage

Полный текст

Introduction Explicit Runge-Kutta difference schemes do not preserve the most important properties of the dynamical systems they approximate. The most painful case is Hamiltonian dynamical systems, where the application of explicit Runge-Kutta schemes leads to the violation of the conservation law of the total mechanical energy [1]. In the 1990s, the concept of geometric integrators emerged, i.e. difference schemes that inherit the symplectic structure of Hamiltonian systems [1-4]. These difference schemes preserve exactly linear and quadratic integrals of motion [1, Cooper’s theorem], but do not preserve more complicated integrals of motion of the system [1, 5]. On the other hand, we can construct difference schemes that precisely preserve motion integrals but do not preserve a symplectic structure. The first such scheme for the many-body problem was proposed by Greenspan [6-9], there are such schemes of arbitrarily large order [10]. © 2026 M. D. Malykh, M. A. Konyaeva This work is licensed under a Creative Commons “Attribution-NonCommercial 4.0 International” license. M. D. Malykh, M. A. Konyaeva Calculation of modified Hamiltonian in Sage 215 For symplectic integrator, the change in the integral of motion on an approximate solution approximating a dynamic system with order also begins with a term of order or higher. Using Richardson’s method [11-17], it is easy to verify that this term has a non-zero coefficient in the overwhelming majority of nonlinear problems. Therefore, in particular, the change in energy turns out to be quite noticeable in the plots obtained from approximate solutions. We will illustrate this with an example in Section 2. Nevertheless, the first computer experiments were carried out with Hamiltonian systems having one degree of freedom. In these systems, it was found every time that the points of the approximate solution on the phase plane are aligned along some line that is close to the line of constant energy level, but slightly different from it. This allows us to hope that the points of the approximate solution lie on the level lines of the modified Hamiltonian, that is, some function ̃( , , ), depending on the step and transforming into the Hamiltonian ( , ) of the original system at → 0. In Ref. [18, n. 10.1.2] it was proved that for any symplectic difference scheme approximating a given Hamiltonian system of order , and any order > , one can compute a modified Hamiltonian ( , , ) = ( , ) + ( )( , ) + ⋯ + ( -1) -1, which is preserved on the approximate solution of order . The problem of finding is solved symbolically. We formulate it as a computer algebra problem in Section 5 and then present our solution as software for Sage. Automating the calculation of the modified Hamiltonian allows us to clarify, using specific examples, what exactly we see when we think that the points of the approximate solution line up along some curves. In Section 6 we will see that the points line up along the level lines of the modified Hamiltonian +1 and it takes considerable effort to discern the change in +1 on the approximate solution. Conservation of mechanical energy on symplectic schemes Consider a Hamiltonian system = , = - (1) on a fixed interval 0 < < . Let the approximate solution ( 0, 0), ( 1, 1), … , ( , ) be found by some difference scheme of order with step = / , and let ( ), ( ) be the exact solution of this system corresponding to the initial conditions Then (0) = 0, (0) = 0. - ( ) = ( ), - ( ) = ( ). The law of energy conservation is satisfied exactly on the exact solution, therefore the change in energy in the interval 0 < < is equal to zero: ( ( ), ( )) - ( 0, 0) = 0. On the approximate solution, the change in energy can be represented as follows: 216 Modeling and Simulation DCM&ACS. 2026, 34 (2), 214-225 1.0 E(∆H) y = 2.01 x - 1.03 10 2 0.5 p 10 3 0.0 0.5 10 4 1.0 2.0 1.5 1.0 0.5 0.0 0.5 1.0 q 10 5 ∆t 10 2 10 1 Figure 1. Trajectory in the phase plane for example 1: exact orbit (solid line) and approximate orbit (dots). Figure 2. Richardson diagram for energy change in example1 that is, the change in energy begins with a term of order or higher = + … . We emphasize that here we consider to be a fixed value, and = / to be a variable. The symplecticity of the difference scheme does not contribute to increasing the order of approximation in this equality. Example 1. For example, let us take the Hamiltonian = 2 2 + 2 3 2 + and let us assume for definiteness that = 0.166. With this choice of parameter and initial conditions = 0, = 1, we obtain a closed trajectory, which is noticeably different from an ellipse. Therefore, we will use these values further. In Figure 1, both the exact orbit , the energy level line, and the approximate orbit found using the midpoint scheme in our fdm for sage [19, 20] system are clearly visible. This scheme is the simplest symplectic scheme of the second order [1]. From the Richardson diagram for the energy change (Figure 2), which is a standard tool in our system for estimating the errors of numerical methods [19], it is clearly seen that = 2 + … and that lg ≃ 1.03. Thus, the energy change is a second-order quantity as for any other second-order scheme. However, Figure 1 shows that the points of the approximate solution in the phase plane line up along some line different from the level line for . The situation described in Example 1 is typical and this allows us to hope that some expression ( , ), depending on the step , is preserved on the approximate solution found using the symplectic scheme: ( , , ) = ( 0, 0, ). The search for such a modified Hamiltonian must be preceded by a more general study of modified equations. M. D. Malykh, M. A. Konyaeva Calculation of modified Hamiltonian in Sage 217 Modified differential equation Let us consider the system of ordinary differential equations = ( 1, … , ), = 1, … , , which for brevity will be written as = ( ). (2) Let be a positive constant, which we will interpret as a time step, and let ̂ be the value of variable at time + . The solution of Eq. (1) can be expanded in a Taylor series ̂ = + ̇ + 1 ̈ 2 + 1 ⃛ 3 + … 2! 3! (3) = + + 1 ( ) 2 + 1 2( ) 3 + … , 2! 3! where the differentiation is expressed as =1 = ∑ . By a difference scheme we mean any system of algebraic equations that relate and ̂: ( 1, … , , ̂1, … , ̂ , ), = 1, … , , which for brevity will be written in the form ( , ̂, ) = 0. (4) Definition 1. A difference scheme (4) is said to approximate the differential equation (1) with order of approximation if substituting in (4) for ̂ its Taylor series expression (3) yields ( , ̂, ) = ( +1). Recall that what is said in this definition is enough to prove that the error in the determination of in one step is of order +1, and in steps is of order [1, 12]. We can save the definition 1 for the case when the right-hand side (1) depends on . Definition 2. A difference scheme (4) is said to approximate the differential equation = ( , ) (5) with order of approximation if substituting for ̂ in (4) its Taylor series expression ̂ = + + 1 ( ) 2 + 1 2( ) 3 + … , 2! 3! gives Here and below ( , ̂, ) = ( +1). = ∑ , . =1 218 Modeling and Simulation DCM&ACS. 2026, 34 (2), 214-225 Statement 1. A difference scheme of order and a natural number > are given. It is required to find a differential equation approximated by this difference scheme with order . Analytical methods for studying finite-difference equations are much less developed than methods for studying differential equations. Having learned to solve the problem 1 symbolically, we can reduce the study of the difference scheme to the study of a “modified” differential equation. Definition 3. Let the difference scheme (4) approximate the differential equation (1) with order , and the differential equation (5) with order > . Then we will call the differential equation (5) a modified differential equation of order , indicating, if necessary, that the original equation is being modified. Algorithm for finding a modified equation Let the difference scheme approximate the differential equation (1) with order . Let us find the modified equation +1 = ( , ) (6) of order + 1. Let us find its right-hand side in the form +1( , ) = ( ) + ( ) . The solution of the modified equation is given by the Taylor series 2! ̂ = + ( + ( ) ) + 1 +1( ( ) + ( ) ) 2 + … , where +1 = ∑( + ) . After removing the brackets, we have =1 ̂ = + ( ) + 1 ( ) 2 + ⋯ + 1 ( ) 2! ! + ( 1 ( + 1)! +1( ) + ( )) +1 + ( +2) This expression up to terms of the order + 1 coincides with the Taylor series for the solution of the original equation (1), so its substitution into the difference scheme (4) for an arbitrary yields series in powers of that start with +1. The coefficients at +1 are linear functions of , so the condition ( , ̂, ) = ( +2) gives exactly linear equations for finding 1, … , . From this system, is uniquely determined as a set of rational functions 1, … , . Thus, if the difference scheme approximates the original differential equation (1) with order , then this scheme approximates with order + 1 the equation (6), the right-hand side of which is a polynomial of degree with respect to . To find the modified equation, it is sufficient to compose and solve the appropriate system of linear algebraic equations. Now it is easy to find the modified equation +2 = ( , ) M. D. Malykh, M. A. Konyaeva Calculation of modified Hamiltonian in Sage 219 of order + 2. We will seek for its right-hand side in the form +2 = +1( , ) + ( ) +1, where is the next correction to the right-hand side of the modified equation, subject to determination. The solution of the modified equation is given by the Taylor series ̂ = + +1( +1) + ⋯ + 1 ( + 1)! +1( +1) +1 +1 + ( 1 ( + 2)! +2( +1) + ( )) +2 + ( +3). This expression up to ( + 2)-th order terms coincides with the Taylor series for solving the modified equation (6), so substituting it into the difference scheme (4) for an arbitrary yields series in powers of that start with +2. The coefficients of +2 are linear functions of , so the condition ( , ̂, ) = ( +3) gives exactly linear equations for finding 1, … , . Thus, we find +2 as a polynomial of degree + 1 in . Proceeding in this way, we can reach any given order . Theorem 1. Let the difference scheme approximate the differential equation (1) with order . Then this scheme approximates with order > the equation = ( , ) the right-hand side of which is a polynomial of degree - 1 in . The problem of finding a modified Hamiltonian We apply the developed theory to the Hamiltonian system (1). Any difference scheme for this system can be written in the form ( ̂, ̂, , , ) = 0, ( ̂, ̂, , , ) = 0. (7) A difference scheme is called symplectic if there exists a function such that the equality ̂ ̂ = + is satisfied on the manifold (7). Using the method described above, for any difference scheme we can construct a modified system = ( , , ), = ( , , ), (8) approximated by the difference scheme (7) with order + 1. If the difference scheme is symplectic [1], then the modified system is Hamiltonian [18, n. 10.1.2]. We denote its Hamiltonian as . The energy conservation law for the system (8) is satisfied exactly: ( ̃( ), ̃( )) = ( 0, 0), 220 Modeling and Simulation DCM&ACS. 2026, 34 (2), 214-225 where ̃( ), ̃( ) is the exact solution of the modified system (8). Therefore, on the approximate solution, the change in energy over the interval 0 < < is equal to = ( , , ) - ( 0, 0, ) = ( , , ) - ( ̃( ), ̃( ), ) = ( ), i.e., the change in energy starts with a term of order or higher = ( ). Thus, the algorithm for finding a modified equation of a given order (Section 4) allows us to find a function ( , , ) that, although not exactly preserved on the approximate solution, is preserved with any given order. Then the energy conservation law of the modified system must be satisfied on the considered scheme with the order . Therefore, it must be true = ( , , ) - ( 0, 0, ) = ( ). Thus, the Hamiltonian of the modified system is the function that is preserved on the difference scheme with the given order . We will call such a function the modified Hamiltonian of the -th order. Statement 2. A symbolic expression , containing two independent variables , , a symplectic difference scheme of order , and an integer > are given. It is required to find a modified Hamiltonian of order , that is, a symbolic expression of the form ( , ) = ( , ) + ( )( , ) + ⋯ + ( -1) -1, which is preserved with order . In Ref. [18, n. 10.1.2], expressions for modified Hamiltonians of order +1 were found for a separate scheme. We implemented an algorithm for solving a more general problem 1 in the Sage computer algebra system [21] and tested its operation on several examples; the program is available in the public repository [22]. Example 2. To verify our program, consider the one-stage split Runge-Kutta scheme ̂ - = - ′( ̂) , ̂ - = ′( ) for a system with Hamiltonian of the form = ( ) + ( ). This scheme is symplectic and has the first order of approximation [23, p. 3.1]. We define it in our system as follows: var('p,q,dt,pp,qq') T=function('T')(p) U=function('U')(q) H=T+U scheme=[pp-p + diff(H,q).subs(q=qq)*dt, qq-q-diff(H,p)*dt] Our function mod_ham has three arguments: the Hamiltonian, the scheme, and the order. In the case under consideration, the function mod_ham(H,scheme,2) returns the modified Hamiltonian 2 = ( ) + ( ) + ′( ) ′( ) 2 M. D. Malykh, M. A. Konyaeva Calculation of modified Hamiltonian in Sage 221 of the second order, and mod_ham(H2,scheme,3) returns the modified Hamiltonian ′ ′ 2 ″ ′ 2 ′ 2 ″ 3 = ( ) + ( ) + 2 ( ) ( ) + 12 ( ( ) ( ) + ( ) ( )) of the third order. This coincides with the expressions specified in [18, n. 10.1.1]. Example 3. Midpoint scheme - = ( , ) , - = - ( , ) , where = + , = + , 2 2 is a symplectic scheme and has 2nd order approximation [1, 5]. We define the scheme in Sage as follows: scheme=[pp-p + diff(H,q).subs(q=(q+qq)/2)*dt, \ qq-q-diff(H,p).subs(p=(p+pp)/2)*dt] The modified 3rd order Hamiltonian is 2 ′ 2 ″ ′ 2 ″ 3 3 = ( ) + ( ) - 24 ( ( ) ( ) + ( ) and coincides with the modified 4th order Hamiltonian. Remark 1. It is worth noting that for a harmonic oscillator ( )) + ( ) 2 2 = 2 + 2 it is true ′ 2 ″ ′ 2 ″ 2 2 2 ( ) ( ) + ( ) ( ) = 2 - 2 = . Thus, the original Hamiltonian is preserved. This is as it should be according to Cooper’s theorem [1]. Potentially, our program allows searching for Hamiltonians of any order. However, for the midpoint scheme, it is not possible to calculate a Hamiltonian of the 5th order, since the Sage system cannot calculate the appropriate integrals in the general form. However, our program supports the specification of the Hamiltonian as any symbolic expression, so when considering particular cases, a specific expression for the energy can be used instead of the general expression. In this case, the integrals can be calculated and we get the opportunity to trace the role of the various terms in (8), not restriciting ourselves to the first one, as is usually done. System with cubic Hamiltonian Consider the midpoint scheme for the system with cubic Hamiltonian = 2+ 2 + 3, the 2 consideration of which we began in the Example 1. We define this scheme in the usual way: var('p,q,dt,pp,qq,a') H=p^2 / 2 + q^2 / 2 + a * q^3 scheme=[pp-p + diff(H,q).subs(q=(q+qq)/2)*dt, \ qq-q-diff(H,p).subs(p=(p+pp)/2)*dt] 222 Modeling and Simulation DCM&ACS. 2026, 34 (2), 214-225 H 1.0 0.665 0.5 0.660 p 0.0 0.655 0.650 0.5 0.645 1.0 2.0 1.5 1.0 0.5 0.0 0.5 1.0 q Figure 3. Phase trajectory for = 0.166 t 0 20 40 60 80 100 Figure 4. Time dependence of (dashed), 3 (dotdash) and 5 (solid) for = 0.166 First, our function allows us to find 3 = - 1 (9 2 4 + 6 2 + 6 3 + 2 + 2)dt2 24 + 1 2 + 1 2 + 3, 2 2 then 4 = 3, and then 5 = 1 160 (54 3 5 + 45 2 4 + 10 2 + 12 3 + (30 2 2 + 1) 2 + 2)dt4 - 1 (9 2 4 + 6 2 + 6 3 + 2 + 2)dt2 24 + 1 2 + 1 2 + 3. 2 2 Figure 3 shows the phase trajectory passing through the point ( , ) = (1, 0). The dots represent the approximate solution found using the midpoint scheme. The true trajectory, i.e. the level line of the Hamiltonian , is shown by the red solid line, and the level line 3 is shown by the green dotted line. It is clearly seen that with graphical accuracy the points of the approximate solution lie on the curve 3( , ) = const, and not on ( , ) = const. Figure 4 shows the change in the Hamiltonians , 3, and 5 on the approximate solution. It is clearly seen that the amplitude of the change in is several hundredths, 3 is several thousandths, and 5 is not visible in the plot at all. Results The presented program allows finding a symbolic expression for a modified Hamiltonian of a given order > for a given symplectic difference scheme of order . We tested our program on a split Runge-Kutta scheme, which was previously investigated without involving computer algebra systems. Then we conducted a series of numerical experiments with a midpoint scheme for a cubic Hamiltonian. In these experiments, it turned out every time that with graphical accuracy the points of the approximate orbit lie on the level lines of the modified Hamiltonian of the 3rd order. It is possible to see that this Hamiltonian is not preserved exactly only on a specially plotted change in energy over time. The change in the modified Hamiltonian of the 5th order over the entire observation period lies within the limits of the calculation error.

Об авторах

М. Д. Малых

Российский университет дружбы народов; Объединённый институт ядерных исследований

Автор, ответственный за переписку.
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, head of the department of Mathematical Modeling and Artificial Intelligence, RUDN University, research fellow of MLIT JINR (Dubna)

ул. Миклухо-Маклая, д. 6, Москва, 117198, Российская Федерация; ул. Жолио-Кюри, д. 6, Дубна, Московская обл., 141980, Российская Федерация

М. А. Коняева

Российский университет дружбы народов

Email: 1032217044@pfur.ru
ул. Миклухо-Маклая, д. 6, Москва, 117198, Российская Федерация

Список литературы

Дополнительные файлы