Discrete and Continuous Models and Applied Computational ScienceDiscrete and Continuous Models and Applied Computational Science2658-46702658-7149Peoples' Friendship University of Russia2369410.22363/2658-4670-2020-28-1-17-34Research ArticleNumerical determination of the singularity order of a system of differential equationsBaddourAli<p>PhD student of Department of Applied Probability and Informatics</p>alibddour@gmail.comMalykhMikhail D.<p>Doctor of Physical and Mathematical Sciences, assistant professor of Department of Applied Probability and Informatics</p>malykhmd-md@rudn.ruPaninAlexander A.<p>Candidate of Physical and Mathematical Sciences, assistant professor of Faculty of Physics</p>a-panin@yandex.ruSevastianovLeonid A.<p>Doctor of Physical and Mathematical Sciences, professor of Department of Applied Probability and Informatics</p>sevastianov-la@rudn.ruPeoples’ Friendship University of Russia (RUDN University)M. V. Lomonosov Moscow State University15122020281173409052020Copyright © 2020, Baddour A., Malykh M.D., Panin A.A., Sevastianov L.A.2020<p>We consider moving singular points of systems of ordinary differential equations. A review of Painlevs 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 A. Alshina (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 Hamiltonian 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.</p>CROSFinite-difference methodssageCalogero systemPainlevé propertyCROSSageметод конечных разностейсистема Калоджеросвойство Пенлеве<p>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 Painlev in the very end of the 19-th century and was presented in his famous Stockholm Lectures [3]. Theorem 1 (Painlev, 1897; [2]). A moving singular point of a solution of an ordinary first-order differential equation (,̇ , ) = 0, ℚ[, , ], (1) is always algebraic, i.e., in the vicinity of the singular point such a solution can be expanded in Puiseux series = ( - ) + , ℚ. (2) Remark 1. Here it is assumed that the Puiseux series converges absolutely and uniformly in a certain vicinity of point = . Remark 2. In some papers on numerical solution of ordinary differential equations such a singularity is referred to as a pole of the order || even when ℤ. A generalization of this result to the case of a system having the form 1 1 = ⋯ = = 0 (3) requires some stipulations. Let 0, , be polynomials from ℚ[, 1, ]. Generally the hypersurfaces 0(1, , , ) = 0, , (1, , , ) = 0 (4) have a finite number of crossing points in the projective space ℙ+1. Fixed singular points of the systems solution can be calculated as projections of these points on the -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: 1 = 1( - ) + , ℚ. (5) Remark 3. All systems, in which 1, do not contain 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 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 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. Marchuks problem was solved by N. N. Kalitkin and E. A. Alshina 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. Alshina 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 (, ) can be expanded in an asymptotic series (, ) = () + (), where = () is the exact solution and is the order of approximation. 3. Therefore, at regular points the ratio (, ) - (, /2) 1 - 1 ≃ 2 = 2. - (, /2) - (, /22) 1 2 1 22 A remarkable fact is that after the moving singularity it equals (, ) - (, /2) (, /2) - (, /22) where is the order of algebraic singularity. = 2, The complex first-order Rosenbrock scheme for autonomous system of differential equations is written as = () (6) = + Re , ⎨ 1 {⎧ ̂ (7) + { ( - ) = (). ⎩ 2 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 () (0 - )-. Then the behavior of the derivative has the form 1 () (0 - )--1 = 1+ . Let us denote this effective right-hand side as (). In this case = + 1 1 1 = ( + 1) . The Rosenbrock schemes have the form (1 - ()) = (), {̂ - = Re , where is the time step of the grid. In application to our right-hand side we get 1 1 1 1+ , (1 - ( + 1) ) = 1+ = 1+ 1 1 1 - ( + 1) ̂ - = Re 1 , 1 - ( + 1) or, for the complex Rosenbrock scheme, which is of interest for us, i.e., for = (1 + )/2, 1 1+ 1+ 1 1 = 1 - 1+ 1 2 ( + 1) = 1 2 1 - 1 ( + 1) 1 = 2 - ( + 1) 1 1 1 1 - 1 ( + 1) + ( + 1) 2 = 1+ 2 2 1 1 2 , 2 (1 - 1 ( + 1) ) 2 + ( 1 ( + 1) ) 1 - 1 ( + 1) 1 1 Re = 1+ 2 . 1 2 1 2 Thus, 2 (1 - 1 ( + 1) ) 2 + ( 1 ( + 1) ) 1 - 1 ( + 1) 1 1 ̂ - = 1+ 2 , 1 2 1 2 2 (1 - 1 ( + 1) ) 2 + ( 1 ( + 1) ) from which it follows that the fixed point of the scheme is 2 0 = ( ( + 1)) . (8) Now let us apply the Richardson estimate of the effective accuracy order, which has the form eff = log (/) - () , (9) (/2) - (/) where is the grid densening factor (we used 3 grids with the step , / and /2). From (8) and (9) we get ( ) 2 (+1) - ( 2 (+1) ) 1 2 ( 2 eff = log (+1) ) - ( ) 2 (+1) = log = -, which demonstrates that the Richardson estimate of the effective accuracy order using the Rosenbrock complex scheme yields the order of the algebraic singularity. 1. 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. 1. CROS Function cros(F,ics) for the Cauchy problem described in terms of its arguments ̇ = (), ∣ = = 0, (10) considered on segment ⩽ ⩽ , calculates the values of column at + 1 points of the grid, covering segment [, ] uniformly with the step = ( - )/(). The function returns a list of values of column , calculated at points + - , = 0, 1, . Natural number indicates by how many times the grid formed by division of the initial segment into parts becomes denser. Necessary arguments: is a list of right-hand sides of the ODE, its element type is symbolic expression, 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: is the initial value of variable , by default = 0 is the finite value of variable , by default = 1, is the number of grid nodes before densening, by default = 10, is the densening index; ultimately the scheme uses the step = ( - )/(), by default = 1. As initial conditions, it is recommended to take the numbers from ℝ, otherwise the calculations become extremely time-consuming, which is a specific feature of the CROS scheme. Example 1. Consider the Cauchy problem = 2, ∣ =0 = 1, ∣ = 1 =0 on the segment 0 3. Let us rewrite it in the form (10) ⎧ {̇ = 1, { {̇ = ⎨ { ̇ = 2, { { ⎩ = 0, = = 1, t=0. When specifying the initial conditions let us use decimal fractions, which will be perceived by the system as numbers from ℝ. 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([1,z,y^2],[x==0, y==1.0, z==1.0], b=3.0)]) Graphics object consisting of 1 graphics primitive y 20 15 10 5 x 0 0.5 1 1.5 2 2.5 3 Figure 1. Solution of the initial problem from example 1 using the CROS scheme at = 10, = 1 (points) and = 4 (line) 2. The order Let = (, , ) and let , , ‴ be three solutions to problem (10) calculated by the CROS scheme. The first solution is obtained with the division of the segment by points, the second by 2, and the third by 22 points. Then the expression - = log2 ∣ ‴ ∣ - before a singular point asymptotically in 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 -th node of the grid. The function eff_order(F,ics) computes the effective order for the solutions of the Cauchy problem ̇ = (), ∣ = = 0, using the CROS scheme on the segment ⩽ ⩽ . The function returns a list of points ( , ). 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: is the initial value of variable , by default = 0, is the final value of variable , by default = 1, is the number of grid nodes, by default = 10 . As initial conditions it is recommended to take the numbers from ℝ, 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 = 2, ∣ =0 = 1, ∣ = 1 =0 from example 10 on the segment 0 3. sage: eff_order([1,z,y^2],[x==0, y==1.0, z==1.0], b=3.0, N=30) [[0.100000000000000, 2.00989109022581], [0.200000000000000, 2.01268790858706], [0.300000000000000, 2.01483966814792], [0.400000000000000, 2.01644106174018], [0.500000000000000, 2.01754741718014], [0.600000000000000, 2.01818124385354], [0.700000000000000, 2.01833386808129], [0.800000000000000, 2.01796281442251], [0.900000000000000, 2.01698449712759], [1.00000000000000, 2.01526061186748], [1.10000000000000, 2.01257491160427], [1.20000000000000, 2.00859410681839], [1.30000000000000, 2.00280107081003], [1.40000000000000, 1.99437737792691], [1.50000000000000, 1.98198860002381], [1.60000000000000, 1.96337299865009], [1.70000000000000, 1.93450889341422], [1.80000000000000, 1.88781839346059], [1.90000000000000, 1.80800548867119], [2.00000000000000, 1.66166416882944], [2.10000000000000, 1.36965695426710], [2.20000000000000, 0.734022444594345], [2.30000000000000, -0.678328975943153], [2.40000000000000, -2.12961258492029], [2.50000000000000, -2.03173073133296], [2.60000000000000, -1.98022752373820], [2.70000000000000, -1.99096481400289], [2.80000000000000, -1.99866158256879], [2.90000000000000, -2.00066044407114]] It is well seen that the singularity is somewhere near = 2.2 and has an order of -2, see also Figure 2. y 2 1.5 1 0.5 0 -0.5 x 0.5 1 1.5 2 2.5 -1 -1.5 -2 Figure 2. Effective order of the solution to the initial problem from example 1 using the CROS scheme 2. Testing the package 1. Examples with integer orders An ODE is said to possess Painlev property if the order 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. = 62 + , { (0) = (1) = 1. (11) Applying the CROS package we see that has a pole of the second order. This result agrees well with V. V. Golubevs theorem (1912) [2, p. 199]. Example 4. = 23 + + 1, { (0) = (1) = 1. (12) We see that has a pole of the first order. This contradicts V. V. Golubevs theorem (1912) [2, p. 200]. However, at present a number of solutions of the second Painlev 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. 2 ⎧{ = ( ) + (2 + 1) + 2 (3 + 1 ) , (13) {⎨ ⎩(0) = (1) = 1. We see that has a pole of the first order. 2. Calogero system Among mechanical systems, the Calogero system [15] is the most suitable for testing. Let us consider 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 be the position of the -th point, then where ̈ = - , (14) = ( - ), () = This system is Hamiltonian with ||2 . = 1 2 + () 2 and, what is most essential for the subsequent analysis, it can be written in the form where ̂ = [,̂ ], (15) 1 - and (, ) = diag(1, 2, , ) + ( 1 - ) - 1 - (, ) = diag(1, 2, , ) - ( ( - )2 ) , = . ( - )2 This representation was guessed and can be checked by direct substitution. Equation (15) means that the eigenvalues of matrix ((), ()) are independent of . Indeed, let us introduce the matrix ̂ () as a solution to the Cauchy problem ̂ = ̂ ̂ , ̂(0) = . A conjugate matrix satisfies the equation ̂ and, therefore, ̂ ̂ = -̂̂ = ̂ ̂̂ - ̂ ̂ ̂ = -[,̂ (̂ ̂ )], ̂ ̂ (0) = . The only solution to this Cauchy problem is a unitary matrix. It remains to note that ̂ ̂̂ ̂ ̂ = , therefore, ̂ is so that = -̂̂ ̂̂ + ̂ [, ]̂ ̂ + ̂ ̂̂ ̂= 0, ̂ ()∣ ̂ () = ∣ , 0 i.e., the eigenvalues of the matrices ̂ ((), ()) and ̂ ((0), (0)) on any solution of (14) coincide, which was to be proved. The proved statement means that the eigenvalues of matrix (, ) are integrals of motion for system (14). It is convenient to use them for constructing symmetric functions (, ) = Sp (, ), ( = 1, ), 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 + will spread (i.e., | - | ). Therefore = + and (, ) 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., (, ) = 0. Hence, the Calogero system is completely integrable and has rational integrals of motion. The description of solution used above can be derived as follows: for the matrix it is valid that () = diag(1, ) ̂ ̂̂ = -̂̂ ̂̂ + ̂ diag(1 , )̂ + ̂ ̂̂ ̂ = 1 = ̂ ([,̂ ] + diag( , ))̂ = ̂ ̂̂ = ̂∣ , =0 from where using the Taylor formula we immediately have ̂ ̂̂ = ∣ =0 + ∣ . =0 Therefore, general solution 1 = 1(), is a set of eigenvalues of the matrix ̂∣ =0 + ̂∣ =0 and is given by algebraic functions of and initial data. Thus, the Calogero system is an example of a system with rational Hamiltonian (, ), which has rational integrals in involution, and whose general solution is provided by algebraic, but not rational functions of and initial data. Two bodies attract, the order of singularity 1/2 (Figure 3): sage: var('t, q1,q2,p1,p2') (t, q1, q2, p1, p2) sage: L=cros([1, p1,p2, diff(1/(q1-q2)^2,q1), diff(1/(q1-q2)^2,q2)], [t==0, q1==0.0, q2==1.0, p1==0.0, p2==0.0], N=10^2) q 1 0.8 0.6 0.4 0.2 t 0.2 0.4 0.6 0.8 1 r 2 1.5 1 0.5 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, q1,q2, q3, p1, p2, p3') (t, q1, q2, q3, p1, p2, p3) sage: V=(1/(q1-q2)^2 + 1/(q1-q3)^2 + 1/(q2-q3)^2) sage: L=cros([1, p1,p2,p3, diff(V,q1),diff(V,q2),diff(V,q3)], [t==0, q1==0.0, q2==1.0, q3==3.0, p1==0.0, p2==0.0, p3==0], b=1, N=10^2) q 3 2.5 2 1.5 1 0.5 r t 0.2 0.4 0.6 0.8 1 6 5 4 3 2 1 t 0.2 0.4 0.6 0.8 1 Figure 4. Collision in a three-body problem 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 ⎧ 3 { = - , ⎨ {∣ 6 = ∣ = 0 ⎩ =0 =0 on the segment 0 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. y 2.4 2.2 2 1.8 1.6 1.4 1.2 1 x 0 0.2 0.4 0.6 0.8 1 Figure 5. Solutions found using the CROS scheme (line) and rk4 (dots) p 3 2 1 x 0.2 0.4 0.6 0.8 1 -1 -2 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 - () = 2 + (3) is used assuming without any justification that 0. In the present case we deal with a single node, in which = 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 yh/2 2e-3 1.5e-3 1e-3 5e-4 x 0.2 0.4 0.6 0.8 1 Figure 7. Plot of difference between two approximate solutions calculated on grids with the step = 0.01 and 0.005 Remark 4. Note that the point = 0 is always such a point of superconvergence, 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. 1. 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 Painlev 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.</p>[W. A. Stein, Sage Mathematics Software (Version 6.7), The Sage Development Team, 2015.][W. W. Golubev, Vorlesungen über Differentialgleichungen im Komplexen. Berlin: VEB Deutscher Verlag der Wissenschaften, 1958.][P. Painlevé, “Leçons sur la theorie analytique des equations differentielles,” in Œuvres de Paul Painlevé. 1973, vol. 1.][C. L. Siegel and J. Moser, Lectures on Celestial Mechanics. Springer, 1995.][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.][A. A. Belov, “Numerical detection and study of singularities in solutions of differential equations,” Doklady Mathematics, vol. 93, no. 3, pp. 334- 338, 2016. DOI: 10.1134/S1064562416020010.][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.][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.][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.][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.][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.][A. Baddour and M. D. Malykh. (2019). Cros for sage, RUDN, [Online]. Available: https://malykhmd.neocities.org.][H. Airault, “Rational solutions of Painlevé equations,” Stud. Appl. Math., vol. 61, no. 1, pp. 31-53, 1979. DOI: 10.1002/sapm197961131.][(2019). Nist digital library of mathematical functions. version 1.0.25, The National Institute of Standards and Technology, [Online]. Available: http://dlmf.nist.gov.][J. Moser, Integrable Hamiltonian Systems and Spectral Theory. Edizioni della Normale, 1983.]