On the realization of explicit Runge–Kutta schemes preserving quadratic invariants of dynamical systems

We implement several explicit Runge–Kutta schemes that preserve quadratic invariants of autonomous dynamical systems in Sage. In this paper, we want to present our package ex.sage and the results of our numerical experiments. In the package, the functions rrk_solve, idt_solve and project_1 are constructed for the case when only one given quadratic invariant will be exactly preserved. The function phi_solve_1 allows us to preserve two specified quadratic invariants simultaneously. To solve the equations with respect to parameters determined by the conservation law we use the elimination technique based on Gröbner basis implemented in Sage. An elliptic oscillator is used as a test example of the presented package. This dynamical system has two quadratic invariants. Numerical results of the comparing of standard explicit Runge–Kutta method RK(4,4) with rrk_solve are presented. In addition, for the functions rrk_solve and idt_solve, that preserve only one given invariant, we investigated the change of the second quadratic invariant of the elliptic oscillator. In conclusion, the drawbacks of using these schemes are discussed.


Quadratic invariant and conservative RK scheme
One of most widespread mathematical models is an autonomous system of ordinary differential equations, i.e., the system of the form where: is an independent variable, commonly interpreted as time; is a vector ( 1 , … , ); is a vector function ( 1 , 2 , … , ), when in applications its element ( = 1, 2, … , ) is often taken as a rational or an algebraic function of the coordinates 1 , … , or one can be reduced to this form by some transformation of variables.
Definition 1 (Goriely [1]). If there exists a function of , such that, for any solution ( ) of system (1), the condition holds, then is called the first integral or invariant of the system (1). If ( ) is any exact solution of system (1) then ( ( )) is independent of . If is a polynomial of degree 2 with respect to then it is called a quadratic invariant.
We will interpret { } as an approximation to the exact solution ( ) at time 0 + Δ , i.e.
For the system (1) Runge-Kutta scheme (RK scheme) with stages can be written as and Below the parameters and ( = 1, 2, ..., , = 1, 2, ..., ) will be arranged in an array 1 11 12 , is known as the Butcher table [2], [3] and will be called the coefficients of the RK scheme. The RK method is explicit if = 0 for ≤ , otherwise is implicit.
In this case the RK method will be called S-conservative RK scheme. Ref.
[4]- [6] indicated that the RK method preserves the quadratic first integrals of system (1) iff the coefficients of such RK method satisfy Such Runge-Kutta methods are called symplectic.
Obviously, no explicit RK schemes satisfies the symplectic condition [6], [7]. Unfortunately, during using the implicit schemes, we must solve a system of non-linear algebraic equations at each step. This is very complex problem, so implicit schemes require more resources than explicit RK schemes [8]. Furthermore numerical solutions of nonlinear system (for ex., by the Newton method) introduce new errors that sometimes we cannot estimate effectively [9]. Thus, the integrals could not be preserved exactly. For this reason, many authors try to construct numerical methods for solving the system of differential equations (1) with the preservation of algebraic integrals without the need to solve nonlinear algebraic equations.
To overcome these difficulties Buono and Mastroserio [10] suggested a method that uses explicit RK schemes for the construction of new finitedifference schemes which exactly preserve invariants. Below we will call it the Buono method for shorthand. Of course, these new schemes are not standard RK schemes, but they are usually called an explicit RK scheme preserving invariants [8]. These schemes preserve only one specified invariant. We implemented several such schemes in Sage and investigated what happens to other invariants. Next, we investigated the method from th article [11] by Calvo et al. which is an extension of the Buono method and can be used as a conservation one or more invariants. Below we will call it the Calvo method for shorthand.

The Buono method
To make the explicit RK scheme conservative, we follow to Buono and Mastroserio [10]. We scale the weights by a parameter ∈ ℝ at the step , i.e. use instead of +1 obtained by the formula (4) numerical solution after one time step.
Using the shorthand the parameter could be estimated by the conservative condition, i.e.
Thus, we preserve the invariant ⟨ , ⟩, if we take For Runge-Kutta schemes of order this expression is close to 1, i.e this is called the increment direction technique (IDT). Note that the value of the scalar parameter at each step depends on the quadratic invariant that appears in rrk or idt method. (3) and (4) has order , then the method defined by (3) and (7) has: -(RRK method) If the solution +1 is interpreted as an approximation to ( + Δ ), the method has order . -(IDT method) If the solution +1 is interpreted as an approximation to ( + Δ ), the method has order − 1.

Theorem 1 (Zhang et al. [8]). Let the original RK scheme be defined by
In our package ex.sage for Sage [12] the function rrk_solve(P1,F,ics) returns the numeric points (0, 0 ), ( 0 Δ , 1 ), ⋯ with the parameters: -P1 is a quadratic invariant; -F is the right sides of system (1) ics is the initial condition -the default Δ = 0.1, = 10 -4-stage explicit RK scheme with the Butcher table Function idt_solve() returns the numerical points (0, 0 ), (Δ , 1 ), ⋯ with the parameters: -P1 is a quadratic invariant; -F is the right sides of the system (1) ics is the initial condition -the default Δ = 0.1, = 10 -4-stage explicit RK scheme with coefficients the same as the previous table. Users can redefine these variables in both functions, for ex., by adding dt=0.01 or new explicit RK method.

Elliptic function test
To test this routine, we investigate a nonlinear oscillator. By the definition of Jacobi functions [13], = sn , = cn , = dn is a particular solution to a nonlinear autonomous system of differential equationṡ with the initial conditions = 0, = = 1 at = 0.
This autonomous system has two quadratic integrals of motion 2 + 2 = 1 and We can solve the autonomous system (10) by rrk or idt methods that preserve only the first or second integral. For certainly, we take = 1/2 and indicated above initial condition.
In figure 1 we can see a graph of the solution found by rrk method with exact conservation of 2 + 2 = 1. Rrk give a condensation of the greed points in those arches of the graph where the curvature has a maximum. The second integral 2 2 + 2 − 1 is not exactly preserved, but its value fluctuates with a small amplitude of 10 − 7 (figure 2). We also use rrk with exact conservation of 2 2 + 2 = 1. In this case first integral grows by leaps bounds and quickly becomes larger than 10 −3 ( figure 3). Thus, the preservation of one integral does not preserve others.
Let's compare the Buono method with the standard rk4 at the same step size (which denote Δ in the rrk method).   The rrk method with the exact conservation of the first invariant 2 + 2 = 1 preserves both integrals better than rk4 (figure 5) but the rrk with the exact conservation of 2 2 + 2 = 1 preserves the first integral worse than rk4 (figure 4).  From the numerical experiments, we came to the conclusion that preserving multiple integrals requires a different approach.

The Calvo method
Let us take two popular explicit RK methods for examples: the RK4 and Euler methods. RK4 method has 4 stages and 4th approximation order. We calculate four axillary quantities at th step and then the quantity which is used in standard way as +1 . Similarly, by the Euler method, we calculate in the step one axillary quantity 1 = and the quantity which is used in standard way as +1 . We can describe this scheme by Butcher table with additional row: Calvo et al. [11] extend the Buono method by coupling these two schemes: in the step we take where ∈ ℝ is the scalar parameter that can be determined by the conservation law, i.e.
Using (12) and (13)  Thus, (15) gives an algebraic equation to determinating the parameter at each step. Calvo et al [11] proved that approximation order of this new scheme is equal to 3. They investigated more general case when coupling any two explicit RK scheme.
If is a quadratic invariant then we have a quadratic equation for , one of the roots of which goes to 0 at Δ → 0 and the other goes to ∞. In numerical experiments we choose the parameter as real number which is close to the value 0. In CAS sagemath, we use symbolic calculation to solve this equation (15) with respect to , and use the function roots to get this value, since our calculation is performed in the ring ℝ (Real Field with 53 bits of precision). So, there is a small error, since in ℚ will get the exact value, but it is very time-consuming.

Elliptic function test
We implement the described scheme in Sage as the function project_1().
Function project_1() returns the numerical points (0, 0 ), (Δ , 1 ), ⋯ with the parameters: -list_of_integral is an invariant required to be conserved; -F is the right sides of the system (1) ics is the initial condition is the step size, is the end point of time . Let us take the elliptic function for example. In figure 6 we can see that the solution founded by the Calvo method with exact conservation of the first invariant 2 + 2 = 1 and the second invariant 2 2 + 2 = 1 gives a condensation of greed points in those arches of the graph where the curvature has a maximum. The second integral 2 2 + 2 = 1 in not be exactly preserved, but its value fluctuates with a small amplitude 10 −7 (figure 7) by the Calvo method with exact conservation of the first invariant 2 + 2 = 1. The Calvo method with exact conservation of the second invariant 2 2 + 2 = 1 shows that the first invariant 2 + 2 = 1 grows with an error of no more than 10 −5 ( figure 8). Thus, preserving one integral does not preserve the other.

Scheme for preserving two invariants
Theoretically, the Calvo method allows to construct schemes that preserve several invariants. We take two pairs of RK methods defined by the following two extended Butcher tables: The first embedded RK method corresponding to the first table is constructed by combining the standard RK4 method, which has order 4, with the Euler method, which has order 1. The second embedded RK method corresponding to second table is constructed by combining the standard RK4 method, which has order 4, with a second order explicit RK scheme.
We are trying to find an explicit RK scheme of the type: where , ∈ ℝ are two scalar parameters, which can be determined by using the conservation laws, i.e.
where 1 and 2 are two invariants of system (1). By definition, we have Thus, (18) gives us a system of two equations with respect to two unknowns and . In numerical experiment, we choose the parameters , as real numbers close to the value 0.
In this way we have a system of algebraic equations for calculating parameters and we must solve this system at each step. Thus, we lose the main advantage of the exact methods. Sage has numerous tools for applying operations on the field of ideals. In our numerical experiments we use the elimination technique based on Gröbner basis [14] to solve (18). Namely, in each step , after constructing multivariate polynomial ideal in variables , generated by (18), we use Sage built-in function elimination_ideal to obtain an univariate equation in variable . The function roots over ring ℝ is used to solve this equation. Substituting one of the value of which is close to 0 to (18), the value of another parameter can be obtained by using the function roots again.
Consider elliptic function test. We construct the function phi_solve_1() to implement the routine described above.
We implement the described scheme in Sage as the function phi_solve_1(). Function phi_solve_1() returns the numerical points (0, 0 ), (Δ , 1 ), ⋯ with the parameters: -list_of_integral is two invariants that are required to be conserved; -F is the right sides of the system (1); ics is the initial condition; -is the step size, is the end point of time .
Let us take the elliptic function as example. The Calvo method allows to preserve both invariants and thus significantly surpasses the other methods presented in the previous sections. Figure 9 shows that the error 2 2 + 2 − 1 remains constant in size at 10 −16 , while figure 10 shows the error 2 + 2 − 1 remains constant in size at 10 −13 . We can say that these errors are due to the implementation of the calculation in solving equation (18) over the ring ℝ instead of the algebraic closed field ℚ in CAS Sage [15]. From this point of view, we can conclude that this method can be considered as a method that exactly preserves exactly both quadratic invariants in the elliptic function test.

Conclusion
We have investigated several implementation of explicit RK schemes that preserve invariants. To preserve one invariant the Buono and Calvo methods require solving one algebraic equation with one unknown at each step. In our example, this equation is quadratic, so we can find its numerical solution without any difficulties. From the numerical experiments, we concluded that the exact conservation of one invariant is not an obstacle for changing of the other invariants. Thus, the conservation of multiple integrals requires a different approach.
The Calvo method which preserves several invariants has a drawback: it requires solving a system of algebraic equations with several unknown variables at each step, i.e. has the same drawback as standard implicit RK methods have. Fortunately, in our tests the system we obtained is much simpler than the system described by the midpoint scheme. Авторами реализовано несколько явных схем Рунге-Кутты, которые сохраняют квадратичные инварианты автономных динамических систем в Sage. В статье представлен пакет ex.sage и результаты численных экспериментов.