CC BY
7
Mathematical modeling
Research article
UDC 519.6
PACS 07.05.Tp, 02.70.Bf
DOI: 10.22363/2658-4670-2020-28-1-17-34
Numerical determination of the singularity order of a system of differential equations
Ali Baddour1, Mikhail D. Malykh1, Alexander A. Panin2, Leonid A. Sevastianov1
1 Department of Applied Probability and Informatics Peoples' Friendship University of Russia (RUDN University)
6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation
2 Department of Mathematics of Faculty of Physics
M. V. Lomonosov Moscow State University Leninskie Gory, Moscow, 119991, Russian Federation
(received: March 16, 2020; accepted: April 10, 2020)
We consider moving singular points of systems of ordinary differential equations. A review of Painlevé's results on the algebraicity of these points and their relation to the Marchuk problem of determining the position and order of moving singularities by means of finite difference method is carried out. We present an implementation of a numerical method for solving this problem, proposed by N. N. Kalitkin and E. A. Al'shina (2005) based on the Rosenbrock complex scheme in the Sage computer algebra system, the package CROS for Sage. The main functions of this package are described and numerical examples of usage are presented for each of them. To verify the method, computer experiments are executed (1) with equations possessing the Painlevé property, for which the orders are expected to be integer; (2) dynamic Calogero system.
This system, well-known as a nontrivial example of a completely integrable Hamil-tonian system, in the present context is interesting due to the fact that coordinates and momenta are algebraic functions of time, and the orders of moving branching points can be calculated explicitly. Numerical experiments revealed that the applicability conditions of the method require additional stipulations related to the elimination of superconvergence points.
Key words and phrases: CROS, Finite-difference methods, sage, Calogero system, Painlevé property
© Baddour A., Malykh M.D., Panin A. A., Sevastianov L.A., 2020
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/
1. Introduction
One of the main problems that arise in numerical analysis of systems of nonlinear ordinary differential equations is the appearance of moving singular points. It should be recalled that linear systems do not have such features and, therefore, the region of existence of the solution is always known in advance. In nonlinear problems, on the contrary, it is never clear beforehand whether a solution is defined for all considered values of the independent variable or not. Euler also noted that when approaching moving singular points of the solution of the Riccati equation, the approximate solution is increasingly deflecting from the exact one, and it has only recently been shown that, despite this, the finite difference method allows searching for the position and orders of moving singularities. In the present paper we report the implementation of this method in the computer algebra system [1], preceded by the necessary theoretical introduction.
2. Moving singularities of the solutions of ordinary
differential equations
If a singular point of the general solution to an ordinary differential equation or a system of such equations depends on the integration constant, it is called a moving singular point [2]. The behavior of the solution in the vicinity of a moving singular point was studied by Painleve in the very end of the 19-th century and was presented in his famous 'Stockholm Lectures' [3].
Theorem 1 (Painleve, 1897; [2]). A moving singular point of a solution of an ordinary first-order differential equation
F(x,x,t) = 0, FeQ[v,x,t], (1)
is always algebraic, i.e., in the vicinity of the singular point such a solution can be expanded in Puiseux series
x = C(t — a)p + ..., (2)
Remark 1. Here it is assumed that the Puiseux series converges absolutely and uniformly in a certain vicinity of point t = a.
Remark 2. In some papers on numerical solution of ordinary differential
equations such a singularity is referred to as a pole of the order IpI even when p£Z.
A generalization of this result to the case of a system having the form
A fn A
requires some stipulations. Let f0,..., fn be polynomials from ^[t, x1,... xn]. Generally the hypersurfaces
fo(x1, . , Xn, ty = 0, . , fn (X1, . , xn, t) = 0 (4)
have a finite number of crossing points in the projective space Pra+1. Fixed singular points of the system's solution can be calculated as projections of these points on the t-axis.
Definition 1. System (3) is called singular, if the system of algebraic equations (4) has an infinite number of solutions.
Theorem 2 (Painlevé, 1897; [3]). A moving singular point of a nonsingular system (3) is always algebraic, i.e., in its vicinity the solution can be expanded in a Puiseux series:
xi = Ci(t-a)p + ..., (5)
Remark 3. All systems, in which f1,... do not contain t explicitly, i.e., autonomous systems, are singular in the sense of Definition 1. However, a similar theorem can be proved for the most important autonomous systems, e.g., the problem of n bodies [4].
3. Evaluation of the order of an algebraic singularity using the finite difference method
If the solution is unknown, it can be found using the finite difference method. In this regard, the following problem naturally arises
Problem [G.I. Marchuk, 2003]. For a given Cauchy problem and interval a < x < b determine the position of moving singular points in this interval and their orders by analyzing one or several approximate solutions of the Cauchy problem.
Probably, many authors believed that this problem has no solution, since, as Euler noted, the finite difference method describes the solution the worse, the closer we get to a singular point. Nevertheless, G. I. Marchuk's problem was solved by N. N. Kalitkin and E. A. Al'shina in 2005 [5]. The proposed method was then approved and developed at Moscow State University [6]-[11]; the paper by A. A. Belov [6] contains historical information of great interest never published earlier.
The solution proposed by E. A. Al'shina can be described as follows.
1. There exist such schemes, e.g., the complex Rosenbrock scheme of the first order (CROS), for which the approximate solution tends to a finite value, when the exact solution has a pole.
2. At regular points the approximate solution x(tn, At) can be expanded in an asymptotic series
x(tn,At) = <p(tn) + r(tn)Atr,
where x = 4>(t) is the exact solution and r is the order of approximation.
3. Therefore, at regular points the ratio
X(xn, A) - x(tn, At/2) ^ 1 — -1 = T x(tn, At/2) — x(tn, At/22) 1 - _1_ •
A remarkable fact is that after the moving singularity it equals
x(tn,At) — x(tn, At/2)
x(tn ,At/2)—x(tn, At/22)
= 2^,
where p is the order of algebraic singularity.
The complex first-order Rosenbrock scheme for autonomous system of differential equations
f = F(-) <6>
is written as
'u = u + rRe k,
l + i dF\ , . (7)
E —
Not perfectly rigorous, but convincing substantiation of the present approach for the CROS scheme: let in the vicinity of a singular point the behavior of the function be described by expression
u(t)~(too —t)-P.
Then the behavior of the derivative has the form
u'(t) - p(t0 — t)-^-1 = fiu1+?.
Let us denote this "effective" right-hand side as F(u). In this case
Fu = $ = (J3+l)u $.
The Rosenbrock schemes have the form
(l — arFu (u))k = F(u), u — u = rRe k,
where r is the time step of the grid.
In application to our right-hand side we get
(l — ar(p+l)u ? )k = pu1+? ^ k =
I3u1+?
pu1+H
l — ar(ft+ l)u1
u — u = rRe
l — ar(ft+ l)u1,
or, for the complex Rosenbrock scheme, which is of interest for us, i.e.,
for a = (l + i)/2,
k =_, =Hu1+i. 1
l — ^ t(P + l)u ? l — 1 + l)u ? — 2 + l)u ?
„ 1+ 1 1 - 1 r(fi + 1)u n + 2 r(fi + 1)u n
= pu p-2-1
(l- 2t(I3+1)U?) + (±T(I3+1)U?)
1 - 2r(fl+1)u
Re k = flu 2
(l - ±T(fi+1)u? ) + (2r(fl+1)u? )
Thus,
~ , i+1 1- 2t(P + Ï)U?
u - U = Tpu P ■
2 •
(l- 2t(/3+1)U-p ) + (2T(/3+1)U-P ) from which it follows that the fixed point of the scheme is
»»^(^r (8)
Now let us apply the Richardson estimate of the effective accuracy order, which has the form
- l u{T/r) -u{T) (q\
Peff = logr u(j/r2) — u(j/r) ' (9)
where r is the grid densening factor (we used 3 grids with the step r, r/r and t/r2). From (8) and (9) we get
2 r
f = ^ f^Pf+% = log, ^ = ^
r(fi+1)) \r(fi+l)
which demonstrates that the Richardson estimate of the effective accuracy order using the Rosenbrock complex scheme yields the order of the algebraic singularity.
4. The package CROS for Sage
The computer algebra system Sage [1] is perfectly suitable for computer experiments with both symbolic and numerical methods. Therefore, we decided to implement the above method of determining the singularity order in this system by means of a small package CROS [12], which allows the calculation of the position and orders of the moving singular points.
4.1. CROS
Function cros(F,ics) for the Cauchy problem described in terms of its arguments
u = F(u), u\ = Uq , (10)
considered on segment a < t < b, calculates the values of column u at NS + l points of the grid, covering segment [a, b] uniformly with the step At = (b — a)/(NS). The function returns a list of values of column u, calculated at points
b — a
a + ^^3, j = 0,l,...N.
Natural number S indicates by how many times the grid formed by division of the initial segment into N parts becomes denser. Necessary arguments:
— F is a list of right-hand sides of the ODE, its element type is symbolic expression,
— ics is a list of initial data, its element type is equality of the form variable == value. The value can be any number or expression, arithmetic manipulations are supported by Sage.
Optional arguments:
— a is the initial value of variable t, by default a = 0
— b is the finite value of variable t, by default b = l,
— N is the number of grid nodes before densening, by default N = l0,
— S is the densening index; ultimately the scheme uses the step At = (b — a)/(NS), by default S = l.
As initial conditions, it is recommended to take the numbers from R, otherwise the calculations become extremely time-consuming, which is a specific feature of the CROS scheme.
Example 1. Consider the Cauchy problem
y" = y2, y\ n = l, y' \ =l
on the segment 0 < x < 3. Let us rewrite it in the form (10)
x = 1, y = z Z = y2,
^x = 0, y = z = 1, t=0.
When specifying the initial conditions let us use decimal fractions, which will be perceived by the system as numbers from R. sage: load('cros.sage') None
sage: var('x,y,z') (x, y, z)
sage: cros([1,z,y"2],[x==0, y==1.0, z==1.0], b=3.0) [[0, 1.00000000000000, 1.00000000000000], [0.300000000000000, 1.33821049499058, 1.37883146513243], [0.600000000000000, 1.81581265489380, 2.05236796585663], [0.900000000000000, 2.53640347402030, 3.28908311525028], [1.20000000000000, 3.68660561742990, 5.67422154586290], [1.50000000000000,
5.58820279611495, 10.4804480613621], [1.80000000000000,
8.65508112242668, 20.0486307394803], [2.10000000000000,
12.8620917982321, 36.1835389093226], [2.40000000000000,
16.9996351976843, 54.5720610523068], [2.70000000000000,
19.8371868069049, 67.2778475556264], [3.00000000000000,
21.3360253982491, 72.9430480385885]]
For plotting (Figure 1) it is possible to use a standard procedure:
sage: sage: point([ [xx,yy] for [xx,yy,zz] in
cros( [l,z,y"2],[x==0, y==l.0, z==l.0], b=3.0)]) Graphics object consisting of 1 graphics primitive
y
Figure 1. Solution of the initial problem from example 1 using the CROS scheme at N = 10,
S = 1 (points) and S = 4 (line)
4.2. The order
Let u = (x,y,...) and let it', it", it"' be three solutions to problem (10) calculated by the CROS scheme. The first solution is obtained with the division of the segment by Appoints, the second by IN, and the third by 22N points. Then the expression
Vi -Vi
Vj -v'j
before a singular point asymptotically in At equals 2, and after the singular point it equals the order of the passed singularity. This expression is referred to as effective order of the scheme at the j-th node of the grid.
The function eff _order (F, ics) computes the effective order for the solutions of the Cauchy problem
it = F(ii), = un,
using the CROS scheme on the segment a <t <b. The function returns a list of points (Xj, Sj).
Necessary arguments:
— F is a list of right-hand sides of the ODE, its element type is symbolic expression,
— ics is a list of initial data, its element type is equality of the form variable == value. The value may be any number or expression, arithmetic manipulations are supported by Sage.
Optional arguments:
— a is the initial value of variable t, by default a = 0,
— b is the final value of variable t, by default b = l,
— N is the number of grid nodes, by default N = l0 .
As initial conditions it is recommended to take the numbers from R, otherwise the computations appear to be extremely time-consuming (a specific feature of the CROS scheme).
Example 2. Let us calculate the effective order for the solution of the Cauchy problem
y" = y2, y\ n = l, y' \ =l
'x=0 'x=0
from example 10 on the segment 0 < x < 3.
sage: eff_order([1, [[0.100000000000000 2.01268790858706], [0.400000000000000, 2.01754741718014], [0.700000000000000, 2.01796281442251], [1.00000000000000, 2.01257491160427], [1.30000000000000, 1.99437737792691], [1.60000000000000, 1.93450889341422], [1.90000000000000, 1.66166416882944], [2.20000000000000, -0.678328975943153] [2.50000000000000, -1.98022752373820], [2.80000000000000, -2.00066044407114]]
z,y-2],[x==0, y==1. , 2.00989109022581] [0.300000000000000, 2.01644106174018], [0.600000000000000, 2.01833386808129], [0.900000000000000, 2.01526061186748], [1.20000000000000, 2.00280107081003], [1.50000000000000, 1.96337299865009], [1.80000000000000, 1.80800548867119], [2.10000000000000, 0.734022444594345], , [2.40000000000000 -2.03173073133296], [2.70000000000000, -1.99866158256879],
0, z==1.0], b=3.0, N=30) , [0.200000000000000, 2.01483966814792], [0.500000000000000, 2.01818124385354], [0.800000000000000, 2.01698449712759], [1.10000000000000, 2.00859410681839], [1.40000000000000, 1.98198860002381], [1.70000000000000, 1.88781839346059], [2.00000000000000, 1.36965695426710], [2.30000000000000, -2.12961258492029], [2.60000000000000, -1.99096481400289], [2.90000000000000,
It is well seen that the singularity is somewhere near x = 2.2 and has an order of —2, see also Figure 2.
y
2 -
0.5 0
-0.5
0.5 1
X
Figure 2. Effective order of the solution to the initial problem from example 1 using the
CROS scheme
i
5. Testing the package
5.1. Examples with integer orders
An ODE is said to possess Painleve property if the order p of all moving singularities is a negative integer [2]. We know all ODEs of the second order possessing this property, and this provides us with appropriate material for testing the developed package.
Example 3.
jr = %2 + i,
U(0)=Sf(1) = 1. ( )
Applying the CROS package we see that y has a pole of the second order. This result agrees well with V. V. Golubev's theorem (1912) [2, p. 199].
Example 4.
y = 2y3 + ty + 1,
J/(0) = y(1) = 1. ( )
We see that y has a pole of the first order. This contradicts V. V. Golubev's theorem (1912) [2, p. 200]. However, at present a number of solutions of the second Painleve equation are known in rational functions, which have poles of the first order [13], [14, §32.8]. Therefore, the statement of the above theorem should be considered erroneous.
Example 5.
y =
"- ^ + (z2 + 1) + e2t (y3 + -),
y V yj (13)
^y(0)=y(1) = 1.
We see that y has a pole of the first order.
5.2. Calogero system
xl2
Among mechanical systems, the Calogero system [15] is the most suitable for testing. Let us consider N material points of unit mass on a straight line, attracting or repulsing each other with the force inversely proportional to the cube of the distance. Let qn be the position of the n-th point, then
= -wn, (14)
where
U = Y.V^ -Qj), V(x) b
i<j
This system is Hamiltonian with
H=\T,Pn + U(q)
and, what is most essential for the subsequent analysis, it can be written in the form
^ = [A,L], (15)
where
L(p,q) = diag(p1 ,P2,... ,pN) + i (-—
\Qj
and
A(p,q) = idiagK,d2,...,dN)-i), 4 = Y;^—j^.
This representation was guessed and can be checked by direct substitution.
Equation (15) means that the eigenvalues of matrix L(p(t),q(t)) are independent of t. Indeed, let us introduce the matrix U(t) as a solution to the Cauchy problem
<^ = AU, U(0) = E.
A conjugate matrix satisfies the equation
dU * ~ -j. = -U* A
dt
and, therefore,
dUU
dt
= AUU* -UU*A = -[A,(UU*)], UU*(0) = E.
The only solution to this Cauchy problem is UU* = E, therefore, U is a unitary matrix. It remains to note that
dU * LU - -
= -U* ALU+U* [L,A]U+U*LAÜ = 0,
so that
U* (t)L\f U(t) = L\Q,
i.e., the eigenvalues of the matrices L(p(t),q(t)) and L (p(0),q(0)) on any solution of (14) coincide, which was to be proved.
The proved statement means that the eigenvalues of matrix L(p, q) are integrals of motion for system (14). It is convenient to use them for constructing symmetric functions
Fk(p,q) = Sp L(p,q)k, (k=l,...N),
which will be rational integrals of motion. These integrals are in involution. It is most simply seen from the fact that the repulsing particles at t ^ will spread (i.e., Iqrj — qk\ ^ to). Therefore
and (Fk, Fr) tends to zero on any trajectory. On the other hand, the Poisson bracket is an integral of motion. Therefore, this expression must be identically zero, i.e., (Fk,Fr) = 0. Hence, the Calogero system is completely integrable and has N rational integrals of motion.
The description of solution used above can be derived as follows: for the matrix
Q(q) = diag(g1 ,...qN)
it is valid that
^U*QU = —U*AQU + U* diag(Pl,... )U + U*QAU = dt
= U * ([A, Q] + diag(pi,... ))U = U * LU = L\±=q , from where using the Taylor formula we immediately have
U * QU = Q\ +tL\ .
Therefore, general solution q1 = q1 (t),... is a set of eigenvalues of the matrix
Q\ +tL\
and is given by algebraic functions of t and initial data.
Thus, the Calogero system is an example of a system with rational Hamil-tonian H(p, q), which has N rational integrals in involution, and whose general
solution is provided by algebraic, but not rational functions of t and initial data.
Two bodies attract, the order of singularity 1/2 (Figure 3):
sage: var('t, ql,q2,pl,p2') (t, ql, q2, pi, p2)
sage: L=cros([l, pl,p2, diff(1/(ql-q2)~2,ql), diff(i/(ql-q2)~2,q2)],
[t==0, ql==0.0, q2==l.0, pl==0.0, p2==0.0], N=10~2) q
r
.....................t
0 0.2 0.4 0.6 0.8 1
Figure 3. Collision in the two-body problem
Three bodies attract, the order of singularity at the collision points 1/2 (Figure 4):
sage: var('t, ql,q2, q3, pi, p2, p3') (t, ql, q2, q3, pi, P2, p3)
sage: V=(l/(ql-q2)"2 + l/(ql-q3)"2 + l/(q2-q3)"2) sage: L=cros([l, pl,p2,p3, diff(V,ql),diff(V,q2),diff(V,q3)], [t==0, ql==0.0, q2==l.0, q3==3.0, pl==0.0, p2==0.0, p3==0], b=l, N=10-2)
1.5
y
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
Figure 4. Collision in a three-body problem
5.3. Points of superconvergence
Experimentally, it was found that in some cases the use of the proposed technique leads to false operation. Thus, e.g. the Cauchy problem
y" = y
r 6' " =0
on the segment 0 < x < 1 using the CROS method is solved facing no singularities. For sure, in Figure 5 we compare solutions found using the CROS and rk4 schemes. However, the plot of effective orders (Figure 6) demonstrates strange peculiarities.
Figure 6. Effective order
The efforts of reducing the step and minimizing the rounding error did not lead to any changes, which suggests that this artifact is not a defect of the program. We explain it as follows. In theory the formula
Vn ~ 000 = 9rAx'2 + 0(Ax3)
is used assuming without any justification that gn ^ 0. In the present case we deal with a single node, in which gn = 0. This fact is well seen in Figure 7. Thus, the diagram of effective order indicates not only moving singular points, but also "removable" singular points, at which a superconvergence of the scheme takes place.
Yh ~ Yhi 2
2e-3
1.5e-3
le-3
5e-4
0.2
..............'
0.4
0.6
0.8
Figure 7. Plot of difference between two approximate solutions calculated on grids with the
step Ax = 0.01 and 0.005
Remark 4. Note that the point x = 0 is always such a point of supercon-vergence, due to which at the beginning of the plot strange fluctuations are always observed.
The revealed effect means that a rigorous substantiation of the method for determining the order of singularities requires elimination of a certain number of special cases, including superconvergence.
6. Conclusion
Numerical experiments convincingly verify the numerical method for determining the position and order of moving singular points of ordinary differential equations, based on the Rosenbrock difference scheme. Moreover, this method allows easy correction of errors in the order of singular points determined in the course of theoretical studies of equations that possess the Painleve property. On the other hand, we revealed a phenomenon of false operation of the algorithm of moving singular point recognition at the points of superconvergence.
Acknowledgments
The authors are grateful to A. A. Belov (Moscow State University) for discussing the topics covered in this article. All calculations presented in the paper are executed in the computer algebra system Sage [1].
References
[1] W. A. Stein, Sage Mathematics Software (Version 6.7), The Sage Development Team, 2015.
[2] W. W. Golubev, Vorlesungen über Differentialgleichungen im Komplexen. Berlin: VEB Deutscher Verlag der Wissenschaften, 1958.
[3] P. Painlevé, "Leçons sur la theorie analytique des equations differentielles," in Œuvres de Paul Painlevé. 1973, vol. 1.
[4] C. L. Siegel and J. Moser, Lectures on Celestial Mechanics. Springer, 1995.
[5] E. A. Al'shina, N. N. Kalitkin, and P. V. Koryakin, "Diagnostics of singularities of exact solutions in computations with error control,"
Computational Mathematics and Mathematical Physics, vol. 45, no. 10, pp. 1769-1779, 2005.
[6] A. A. Belov, "Numerical detection and study of singularities in solutions of differential equations," Doklady Mathematics, vol. 93, no. 3, pp. 334338, 2016. DOI: 10.1134/S1064562416020010.
[7] A. A. Belov, "Numerical diagnostics of solution blowup in differential equations," Computational Mathematics and Mathematical Physics, vol. 57, no. 1, pp. 122-132, 2017. DOI: 10.31857/S004446690002533-7.
[8] M. O. Korpusov, D. V. Lukyanenko, A. A. Panin, and E. V. Yushkov, "Blow-up for one Sobolev problem: theoretical approach and numerical analysis," Journal of Mathematical Analysis and Applications, vol. 442, no. 2, pp. 451-468, 2016. DOI: 10.26089/NumMet.v20r328.
[9] M. O. Korpusov, D. V. Lukyanenko, A. A. Panin, and E. V. Yushkov, "Blow-up phenomena in the model of a space charge stratification in semiconductors: analytical and numerical analysis," Mathematical Methods in the Applied Sciences, vol. 40, no. 7, pp. 2336-2346, 2017. DOI: 10.1002/mma.4142.
[10] M. O. Korpusov and D. V. Lukyanenko, "Instantaneous blow-up versus local solvability for one problem of propagation of nonlinear waves in semiconductors," Journal of Mathematical Analysis and Applications, vol. 459, no. 1, pp. 159-181, 2018. DOI: 10.1016/j.jmaa.2017.10.062.
[11] M. O. Korpusov, D. V. Lukyanenko, A. A. Panin, and G. I. Shlyapugin, "On the blow-up phenomena for a one-dimensional equation of ion-sound waves in a plasma: analytical and numerical investigation," Mathematical Methods in the Applied Sciences, vol. 41, no. 8, pp. 2906-2929, 2018. DOI: 10.1002/mma.4142.
[12] A. Baddour and M. D. Malykh. (2019). Cros for sage, RUDN, [Online]. Available: https://malykhmd.neocities.org.
[13] H. Airault, "Rational solutions of Painlevé equations," Stud. Appl. Math., vol. 61, no. 1, pp. 31-53, 1979. DOI: 10.1002/sapm197961131.
[14] (2019). Nist digital library of mathematical functions. version 1.0.25, The National Institute of Standards and Technology, [Online]. Available: http://dlmf.nist.gov.
[15] J. Moser, Integrable Hamiltonian Systems and Spectral Theory. Edizioni della Normale, 1983.
For citation:
A. Baddour, M. D. Malykh, A. A. Panin, L. A. Sevastianov, Numerical determination of the singularity order of a system of differential equations, Discrete and Continuous Models and Applied Computational Science 28 (1) (2020) 17-34. DOI: 10.22363/2658-4670-2020-28-1-17-34.
Information about the authors:
Baddour, Ali — PhD student of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University) (e-mail: [email protected], phone: +7(495)9550927, ORCID: https://orcid.org/0000-0001-8950-1781)
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: [email protected], phone: +7(495)9550927, ORCID: https://orcid.org/0000-0001-6541-6603)
Panin, Alexander A. — Candidate of Physical and Mathematical Sciences, assistant professor of Faculty of Physics, M. V. Lomonosov Moscow State University (MSU) (e-mail: [email protected], phone: +7(495)9391033, ORCID: https://orcid.org/0000-0003-1185-4717)
Sevastianov, Leonid A. — Doctor of Physical and Mathematical Sciences, professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University) (e-mail: [email protected], phone: +7(495)9522572, ORCID: https://orcid.org/0000-0002-1856-4643)
УДК 519.6
PACS 07.05.Tp, 02.70.Bf
DOI: 10.22363/2658-4670-2020-28-1-17-34
Численное определение порядка особенности системы дифференциальных уравнений
Али Баддур1, М. Д. Малых1, А. А. Панин2, Л. А. Севастьянов1
1 Кафедра прикладной информатики и теории вероятностей Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия 2 Кафедра математики физического факультета Московский государственный университет им. М. В. Ломоносова Ленинские горы, Москва, 119991, Россия
В статье рассматриваются подвижные особые точки систем обыкновенных дифференциальных уравнений. Дан обзор результатов Пенлеве об алгебраично-сти этих точек и их связи с задачей Г. И. Марчука об определении положения и порядка подвижных особых точек по методу конечных разностей. Представлена реализация численного метода решения этой задачи, предложенная Н. Н. Ка-литкиным и Е. А. Альшиной (2005) на основе комплексной схемы Розенброка, в системе компьютерной алгебры Sage — пакет CROS for Sage. Описаны основные функции этого пакета, приведены численные примеры использования каждой из них. В целях верификации метода проведены компьютерные эксперименты: (1) с уравнениями, обладающими свойством Пенлеве, для которых порядки должны получаться целыми числами; (2) с динамической системой Калоджеро.
Эта система, хорошо известная как нетривиальный пример вполне интегрируемой гамильтоновой системы, в данном контексте интересна тем, что координаты и импульсы являются алгебраическими функциями времени, причём порядки подвижных точек ветвления можно вычислить явно. В рамках численных экспериментов обнаружено, что условия применимости метода требуют дополнительных оговорок, связанных с исключением точек суперсходимости.
Ключевые слова: CROS, метод конечных разностей, Sage, система Калоджеро, свойство Пенлеве
