CC BY
15
Research article
UDC 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2020-28-4-327-345
On the realization of explicit Runge—Kutta schemes preserving quadratic invariants of dynamical systems
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: October 23, 2020; accepted: November 12, 2020)
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 Grobner 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.
Key words and phrases: Explicit Runge-Kutta method, quadratic invariant, dynamical system, Sage
1. 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
dx
Tt=f(x)' i>0' (1)
,x(0) = x0,
where: t is an independent variable, commonly interpreted as time; x is a vector (x1,..., xn); f is a vector function (f1, f2,..., fn), when in applications
© YingY., Malykh M.D., 2020
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4-0/
its element fi (i _ 1,2,...,n) is often taken as a rational or an algebraic function of the coordinates x1,... ,xn or one can be reduced to this form by some transformation of variables.
Definition 1 (Goriely [1]). If there exists a function I of x, such that, for any solution x(t) of system (1), the condition
VI(x)f(x) _ £
n dxk dl
k=i dt dxk k=i
Y,fk (Mt))
dl
dxu
_0,
holds, then I is called the first integral or invariant of the system (1). If x(t) is any exact solution of system (1) then I(x(t)) is independent of t. If I is a polynomial of degree 2 with respect to x then it is called a quadratic invariant.
Any quadratic invariant after a linear transformation can be rewritten in the form
I(x(t)) _ x(t)TSx(t) _ (Sx(t), x(t)) _ const, (2)
where (■ , ■) denote the Euclidean inner product on and S £ is
a symmetric, constant matrix.
To determine an uniform grid (with a step At) of the time interval [0, T] we take
tn _ nAt (n_0,...,N).
We will interpret {xn} as an approximation to the exact solution x(t) at time t0 + n At, i.e.
x(t0 + nAt) & xn.
For the system (1) Runge-Kutta scheme (RK scheme) with s stages can be written as
kt = xn + At£aZJf(kJ), i _ 1,2,... ,s
3=1
and
+ At£bz f(kz).
(3)
(4)
¿=i
Below the parameters airj and (i _ l,2,...,s,j _ \,2,...,s) will be arranged in an array
C1 a11 aV2 * ' a1s
«21 a22 ' ' a2s
cs as1 as2 ' ass
h * • K
where
Ci ^ y aij ,
3=1
is known as the Butcher table [2], [3] and will be called the coefficients of the RK scheme. The RK method is explicit if a^ = 0 for i < j, otherwise is implicit.
We want to construct numerical solutions x°,x1 ,...,xN such that the quadratic invariant I(x) is preserved numerically, i.e.
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 finite-difference 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.
2. Explicit RK scheme of preserving one quadratic
invariant
To make the explicit RK scheme conservative, we follow to Buono and Mastroserio [10]. We scale the weights bi by a parameter G \R at the step
n, i.e. use
(Sxn,xn) = (Sx°,x°) i = I,., N.
(5)
bzaij + biaii - bib1 =0, i,j = l,...,s.
(6)
2.1. The Buono method
instead of xn+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.
Vn+1 =^ AtO./S>rn A ) + jnAt(An,An)). (8)
For Runge-Kutta schemes of order p this expression is close to 1, i.e
as At ^ 0, see Buono at al. [10, prop. 4] for p = 4 and Zhang at al. [8, lemma 3.3]. Thus, the new numerical solution x"^1 can be considered as an approximation either to x(tn + At) with the RK weights scaled by 7n, or to x(tn + jnAt) with the time At scaled by 7n. Zhang at al. [8] denote the method defined by (3) and (7) with the interpretation
as the relaxation Runge-Kutta method (RRK), while the method using
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.
Theorem 1 (Zhang et al. [8]). Let the original RK scheme be defined by (3) and (4) has order p, then the method defined by (3) and (7) has:
— (RRK method) If the solution xrl'+1 is interpreted as an approximation to x(tn + 7n At), the method has order p.
— (IDT method) If the solution xrl'+1 is interpreted as an approximation to x(tn + At), the method has order p — 1.
In our package ex.sage for Sage [12] the function rrk_solve(P1,F,ics) returns the numeric points (0,x°), (j°At,x1),-- with the parameters:
— P1 is a quadratic invariant;
— F is the right sides of system (1)
— ics is the initial condition
— the default At = 0.1, T =10
— 4-stage explicit RK scheme with the Butcher table
A. = ^ f(kz)
(9)
= 1 + 0(AtP-1
x^1 «x(tn + jn At)
x1l+1 « x(tn + At)
0 0 0 0 0
1 2 1 2 0 0 0
1 2 0 1 2 0 0
1 0 0 1 0
1 1 1 1
6 3 3 6
Function idt_solve() returns the numerical points (0, x°), (At, x1), ••• with the parameters:
— P1 is a quadratic invariant;
— F is the right sides of the system (1)
— ics is the initial condition
— the default At = 0.1, T =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.
2.2. Elliptic function test
To test this routine, we investigate a nonlinear oscillator. By the definition of Jacobi functions [13], p = sn t, q = cn t, r = dn t is a particular solution to a nonlinear autonomous system of differential equations
p = qr, q = —pr, r = —k2pq (10)
with the initial conditions
p = 0, q = r = 1 at t = 0.
This autonomous system has two quadratic integrals of motion
p2 + q2 = 1 and k2p2 + r2 = 1 (11)
We can solve the autonomous system (10) by rrk or idt methods that preserve only the first or second integral. For certainly, we take k =1/2 and indicated above initial condition.
sage: var( 'p,q,r')
sage: load( 'ex.sage')
sage: k=1/2
sage: s=4
sage: F=[r*q,-p*r,-k~2*p*q]
sage: list_of_integral=[p~2+q~2,k~2*p~2+r~2]
sage: ics=[p==0,q==1,r==1]
sage: idt_solve(P1=list_of_integral[0],ics=ics,F=F,dt=0.2,T=1)
sage: rrk_solve(P1=list_of_integral[0],ics=ics,F=F,dt=0.2,T=1) sage:
^ B1=rrk_solve(P1=list_of_integral[0],ics=ics,F=F,dt=0.1,T=40) sage:
^ B2=rrk_solve(P1=list_of_integral[1],ics=ics,F=F,dt=0.1,T=40)
sage: G=line([[t, p] for [t,p,q,r] in ^ Bl] ,color= 'red',axes_labels=[ '£ ', 'p '] )+point( [ [t,p] for ^ [t,p,q,r] in B2],frame=true) sage: max(abs(k~2*p~2+r~2-l) for [t,p,q,r] in B2 ) sage: max(abs(p"2+q"2-l) for [t,p,q,r] in B2) sage: max(abs(k"2*p"2+r"2-l) for [t,p,q,r] in Bl ) sage: max(abs(p"2+q"2-l) for [t,p,q,r] in Bl) sage: Gl=line([[t,k~2*p~2+r~2-l] for [t,p,q,r] in - Bl] ,axes_labels=[ '$t$ ', ' $k~2p~2+r~2-^ 1$ '] ,tick_formatter=[None,RR(10e-^ 18).n(digits=l)],frame=true) sage: G2=line([[t,p"2+q"2-l] for [t,p,q,r] in ^ B2],axes_labels=[ '$t$ ', '$p~2+q~2-^ 1$ '] ,tick_formatter=[None,RR(10e-^ 18).n(digits=l)],frame=true)
In figure 1 we can see a graph of the solution found by rrk method with exact conservation of p2 + q2 = 1. Rrk give a condensation of the greed points in those arches of the graph where the curvature has a maximum.
t
Figure 1. Graph of p(t), rrk, dt = 0.1
The second integral k2p2 + r2 — 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 k2p2 + r2 = 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 dt (which denote At in the rrk method).
sage: var( 'p, q,r, t ')
sage: QRK=desolve_system_rk4(F,[p,q,r], ics= [0,0,1,1,], ivar=t, ^ end_points=20)
sage: G3=line([[t,p~2+q~2] for [t,p,q,r] in ^ B2],color= 'red')+point([[t,p~2+q~2] for [t,p,q,r] in ^ QRK],axes_labels=[ '$t$', '$p~2+q~2$'], tick_formatter=[None, ^ RR(10e-18).n(digits=l)],frame=true) sage: G4=line([[t,k"2*p"2+r"2-1] for [t,p,q,r] in B1],axes_labels=[ '$t$', '$k~2p~2+r~2-1$'],tick_formatter=[None, ^ RR(10e-18).n(digits=1)], color= 'red') +
^ point([[t,k~2*p~2+r~2-1] for [t,p,q,r] in QRK],frame=true)
5.0 x 10"!
0.00
-5.0 x 10_!
i-H -1.0 x 10"7
-1.5 x 10"7
+
-2.0 x 10"7
-2.5 x 10"7
-3.0 x 10"7
-3.5 x 10"7
V V
ll II II 11 H H V
10 15 20 25 30 35 t
Figure 2. Second invariant for rrk method with exact conservation of first invariant
q2 = 1
o
5
40
0.00099-
0.00080-
i-H
J 0.00060 -
OI +
®< 0.000400.00020 -
0.00 4 ........................................
0 5 10 15 20 25 30 35 40 t
Figure 3. First integral for rrk with exact conservation of k2p2 + r2 = 1
The rrk method with the exact conservation of the first invariant p2 + q2 = 1 preserves both integrals better than rk4 (figure 5) but the rrk with the exact conservation of k2p2 + r2 = 1 preserves the first integral worse than rk4 (figure 4).
i.o -
cs
&< 1.0 -+ -
CN ft,
i.o -
i.o -
1.0 -
Figure 4. First invariant for rrk method with exact conservation of second invariant k2p2 + r2 = 1 (red) and for standard rk4
0.00 -1.0 x HT7 i-H -2.0 x HT7
I
^ -3.0 xHT7 -4.0 x HT7 -5.0 x HT7 -6.0 x HT7
Figure 5. Second invariant for rrk method with exact conservation of first invariant p2 + q2 = 1 (red) and for standard rk4 (blue)
From the numerical experiments, we came to the conclusion that preserving multiple integrals requires a different approach.
V.
v
0 5 10 15 20 25 30 35 40
t
3. The Calvo method for one invariant
3.1. 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 nth step
¿-> -ry»^
rvi —JL ,
h —xn + 1 f(h )At, k3 —xn + 1 f(k2 )At, k4—xn + f(k3 )At
and then the quantity
t(xn) = X- + At (i/fo) + 3f(k2) + 1/(^3) + 6/(*4)) (12)
which is used in standard way as xn+1. Similarly, by the Euler method, we calculate in the step n one axillary quantity k1 = xn and the quantity
(xn)=xn + f(k1 )At (13)
which is used in standard way as xn+1. We can describe this scheme by Butcher table with additional row:
0 0 0 0 0
1 2 1 2 0 0 0
1 2 0 1 2 0 0
0 0 0
1 1 1 1 1 1
9 6 3 3 6
01 1 0 0 0
Calvo et al. [11] extend the Buono method by coupling these two schemes: in the step n we take
XI+1 = ^(xn) - Xn(<f>(xn) - (xn)), (14)
where \n G \R is the scalar parameter that can be determined by the conservation law, i.e.
I(xl+1) = l(xn). (15)
Using (12) and (13) we have
4+1 = X- + (6 - Ik) №1 )At + (3 + ^) fih)At+
+ (1 + ) f(.h)At + (6 + 1xn) f(k4)At. (16)
Thus, (15) gives an algebraic equation to determinating the parameter \n 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 I is a quadratic invariant then we have a quadratic equation for \n, one of the roots of which goes to 0 at At ^ 0 and the other goes to to. In numerical experiments we choose the parameter \n 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 \n, and use the function roots to get this value, since our calculation is performed in the ring R (Real Field with 53 bits of precision). So, there is a small error, since in Q will get the exact value, but it is very time-consuming.
3.2. Elliptic function test
We implement the described scheme in Sage as the function project_1().
Function project_1() returns the numerical points (0, x°), (At, x1), ••• 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
— dt is the step size,T is the end point of time t. Let us take the elliptic function for example.
sage: load( 'ex.sa.ge') sage: var( 'p,q,r,k') sage: k=1/2
sage: list_of_integral=[p~2+q~2,k~2*p~2+r~2] sage: F=[r*q,-p*r,-k~2*p*q] sage: ics=[p==0,q==1,r==1]
sage: B1=project_1(list_of_integral[0],F,ics,dt=0.1,T=20) sage: max([abs(p~2+q~2-1) for [t,p,q,r] in B1]) sage: P=line([[t,p] for [t,p,q,r] in ^ B1],axes_labels=[ '$t$', '$p$'],tick_formatter=[None,RR(10e-^ 15).n(digits=1)],frame=true, color='red')+point([[t,p] for ^ [t,p,q,r] in B1]) sage: P1=line([[t,p~2+q~2-1] for [t,p,q,r] in ^ B1],axes_labels=[ '$t$', '$p~2+q~2-^ 1$'],tick_formatter=[None,RR(10e-^ 15).n(digits=1)],frame=true, ^ color='red') sage: P2=line([[t,k"2*p"2+r"2-1] for [t,p,q,r] in ^ B1],axes_labels=[ '$t$', '$k~2*p~2+r~2-^ 1$'],tick_formatter=[None,RR(10e-^ 15).n(digits=1)],frame=true) sage: B2=project_1(list_of_integral[1],F,ics,dt=0.1,T=20) sage: max([abs(p~2+q~2-1) for [t,p,q,r] in B2]) sage: P0=line([[t,p] for [t,p,q,r] in ^ B2],axes_labels=[ '$t$', '$p$'],tick_formatter=[None,RR(10e-^ 15).n(digits=1)],frame=true, color='red')+point([[t,p] for ^ [t,p,q,r] in B2])
sage: P01=line([[t,p~2+q~2-l] for [t,p,q,r] in ^ B2] , axes_labels= [ ' $t$ ', '$p~2+q~2-^ 1$'1 ,tick_formatter=[None,RR(10e-^ 15).n(digits=l)],frame=true) sage: P02=line([[t,k's2*p~2+r's2-l] for [t,p,q,r] in ^ B2] ,axes_labels=[ '$t$ ', 1 $k~2*p~2+r~2-^ 1$'1 ,tick_formatter=[None,RR(10e-^ 16).n(digits=l)],frame=true)
In figure 6 we can see that the solution founded by the Calvo method with exact conservation of the first invariant p2 + q2 = 1 and the second invariant k2p2 + r2 = 1 gives a condensation of greed points in those arches of the graph where the curvature has a maximum.
t
t
Figure 6. Graph of p(t), for the Calvo method with exact conservation of p2 + q2 = 1, dt = 0.1(up one) and that of k2p2 + r2 = 1 (down one)
The second integral k2p2 + r2 = 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 p2 + q2 = 1. The Calvo method with exact conservation of the second invariant k2p2 + r2 = 1 shows that the first invariant p2 + q2 = 1 grows with an error of no more than 10-5 (figure 8). Thus, preserving one integral does not preserve the other.
t
t
Figure 7. First and second integrals for the Calvo method with exact conservation of
p2 + q2 = 1
0.000025
0.000020
T—i
I
^ 0.000015 +
N ft,
0.000010 5.0 x 10-6
0.00
0 5 10 15 20
t
6.0 x 10~16
4.0 x 10~16 f 2.0 x 10~16
I
0.00
+
"a, -2.o x io"16 *
cs
-4.0 x 10~16 -6.0 x 10~16 -8.0 x 10~16
Figure 8. First and second integrals for the Calvo method with exact conservation of
k2p2 +r2 = 1
4. 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:
t
0 0 0 0 0 0 0 0 0 0
1 2 i 2 0 0 0 1 2 1 2 0 0 0
1 2 0 1 2 0 0 and 1 2 0 1 2 0 0
1 0 0 1 0 0 1 0 1 0
fa 1 1 1 1 fa 1 1 1 1
6 3 3 6 6 3 3 6
fa 1 0 0 0 0 1 0 0
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:
x^1 = <f>(xn) - (<f>(xn) - 0i(xn)) - /3n(<f>(xn) - 02(xn)), (17)
where an,fin <G R are two scalar parameters, which can be determined by using the conservation laws, i.e.
Thus, (18) gives us a system of two equations with respect to two unknowns an and . In numerical experiment, we choose the parameters an,fin 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 Grobner basis [14] to solve (18). Namely, in each step n, after constructing multivariate polynomial ideal in variables an, generated by (18), we use Sage built-in function elimination_ideal to obtain an univariate equation in variable . The function roots over ring R is used to solve this equation. Substituting one of the value of which is close to 0 to (18), the value of another parameter an 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, x0), (At, x1), ••• 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;
— dt is the step size,T is the end point of time t.
(18)
and
<f>(xn) - 01 (xn ) = (~5f( ki) + 3f(k2 ) + 1/(^3 ) + 1/(*4)) At 0(x-) - 02(xn) = (1ki - 2^2 + 1*3 + 1*4) At.
Let us take the elliptic function as example.
sage: var( 'p,q,r,t') sage: k=1/2 sage: s=4
sage: F=[r*q,-p*r,-k~2*p*q]
sage: list_of_integral=[p~2+q~2,k~2*p~2+r~2] sage: ics=[p==0,q==1,r==1]
sage: L=phi_solve_1(list_of_integral=list_of_integral, F=F, ^ ics=ics, dt=0.1, T=20) sage: max([abs(k~2*p~2+r~2-1) for [t,p,q,r] in L]) sage: G4=line([[t,k~2*p~2+r~2-1] for [t,p,q,r] in ^ L],axes_labels=[ '$t$', '$k~2*p~2+r~2-^ i$'],tick_formatter=[None,RR(10e-^ 20).n(digits=1)],frame=true) sage: max([abs(p~2+q~2-1) for [t,p,q,r] in L]) sage: G5=line([[t,p~2+q~2-1] for [t,p,q,r] in ^ L],axes_labels=[ '$t$', '$p~2+q~2-^ i$'],tick_formatter=[None,RR(10e-^ 20).n(digits=1)],frame=true)
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 k2p2 + r2 — 1 remains constant in size at 10-16, while figure 10 shows the error p2 + q2 — 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 R instead of the algebraic closed field Q 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.
-5.0 x 10
+
CN
-1.0x10
-1.5 x 10
10
t
15
20
0
5
Figure 9. The first integral k2p2 + r2 - 1 for T = 20, dt = 0.1
9.9 x 1CT13
8.0 x 1013 6.0 x 10-13
I
CN
^ 4.0 x 10-13
CN
2.0 x 1013 0.00
-2.0 x 10~13
Figure 10. The second integral p2 + q2 — 1 for T = 20, dt = 0.1
5. 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.
10
t
15
20
0
5
Acknowledgments
The authors are thanks Prof. L. A. Sevastianov for his discuss and suggestions with this paper. This paper was supported by the foundation of KaiLi university (Z1602) and by the Russian Foundation for Basic Research (RFBR) according to the research project No 18-07-00567.
References
[1] A. Goriely, Integrability and nonintegrability of dynamical systems. Singapore; River Edge, NJ: World Scientific, 2001. DOI: 10.1142/3846.
[2] J. C. Butcher, "On Runge-Kutta processes of high order," Journal of the Australian Mathematical Society, vol. 4, no. 2, pp. 179-194, 1964. DOI: 10.1017/S1446788700023387.
[3] E. Hairer, G. Wanner, and S. P. N0rsett, Solving ordinary differential equations I: Nonstiff Problems, 3rd ed. New York: Springer, 1993, vol. 1. DOI: 10.1007/978-3-540-78862-1.
[4] K. Burrage and J. C. Butcher, "Stability criteria for implicit Runge-Kutta methods," SIAM Journal on Numerical Analysis, vol. 16, no. 1, pp. 46-57, 1979. DOI: 10.1137/0716004.
[5] G. Cooper, "Stability of Runge-Kutta methods for trajectory problems," IMA journal of numerical analysis, vol. 7, no. 1, pp. 1-13, 1987. DOI: 10.1093/imanum/7.1.1.
[6] E. Hairer, G. Wanner, and C. Lubich, Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations. Berlin Heidelberg New York: Springer, 2000.
[7] J.-M. Sanz-Serna and M.-P. Calvo, Numerical hamiltonian problems. Courier Dover Publications, 2018.
[8] H. Zhang, X. Qian, J. Yan, and S. Song, "Highly efficient invariant-conserving explicit Runge-Kutta schemes for nonlinear Hamiltonian differential equations," Journal of Computational Physics, p. 109 598, 2020.
[9] W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, et al., Numerical recipes. Cambridge University Press Cambridge, 1989.
[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. DOI: 10.1016/S0377-0427(01)00398-3.
[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. DOI: 10.1137/04061979X.
[12] Y. Ying. (2020). "Package ex.sage," [Online]. Available: https : / / malykhmd.neocities.org.
[13] Y. S. Sikorsky, Elements of the theory of elliptic functions with applications to mechanics [Elementy teoriy ellipticheskikh funktsiy s prilozheniyami v mekhanike]. Moscow-Leningrad: ONTI, 1936, In Russian.
[14] D. Cox, J. Little, D. O'Shea, and M. Sweedler, "Ideals, varieties, and algorithms," American Mathematical Monthly, vol. 101, no. 6, pp. 582586, 1994.
[15] The Sage Developers, Sagemath, the Sage Mathematics Software System (Version 7.4), https://www.sagemath.org, 2016.
For citation:
Y. Ying, M. D. Malykh, On the realization of explicit Runge-Kutta schemes
preserving quadratic invariants of dynamical systems, Discrete and Continuous
Models and Applied Computational Science 28 (4) (2020) 327-345. DOI: 10.22363/2658-4670-2020-28-4-327-345.
Information about the authors:
Ying, Yu — Assistant Professor of Department of Algebra and Geometry, Kaili University (e-mail: yingy6165@gmail.com, phone: +7(495)9550927)
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-2020-28-4-327-345
О реализации явных схем Рунге—Кутты с сохранением квадратичных инвариантов динамических систем
Юй Ин1, М. Д. Малых2
1 Университет Кайли Кайли, 556011, Китай 2 Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия
Авторами реализовано несколько явных схем Рунге—Кутты, которые сохраняют квадратичные инварианты автономных динамических систем в Sage. В статье представлен пакет ex.sage и результаты численных экспериментов.
В пакете функции rrk_solve, idt_solve и project_1 построены для случая, когда только один заданный квадратичный инвариант будет сохранён точно. Функция phi_solve_1 позволяет сохранить одновременно два указанных квадратичных инварианта. Для решения уравнений относительно параметров, определяемых законом сохранения, использована методика исключения на основе базисов Грёб-нера, реализованная в Sage. В качестве тестового примера представленного пакета используется эллиптический осциллятор. Эта динамическая система имеет два квадратичных инварианта. Представлены численные результаты сравнения стандартного явного метода Рунге—Кутты RK(4,4) с rrk_solve. Кроме того, для функций rrk_solve и idt_solve, сохраняющих только один инвариант, исследовано изменение второго квадратичного инварианта эллиптического осциллятора. В заключение рассматриваются недостатки использования этих схем.
Ключевые слова: явный метод Рунге—Кутты, квадратичный инвариант, динамическая система, Sage
