Numerical determination of the singularity order of a system of differential equations

Full Text

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 system’s 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. 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 (, Δ) 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. Golubev’s 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. Golubev’s 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.

About the authors

Ali Baddour

Peoples’ Friendship University of Russia (RUDN University)

Author for correspondence.
Email: alibddour@gmail.com

PhD student of Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation

Mikhail D. Malykh

Peoples’ Friendship University of Russia (RUDN University)

Email: malykhmd-md@rudn.ru

Doctor of Physical and Mathematical Sciences, assistant professor of Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation

Alexander A. Panin

M. V. Lomonosov Moscow State University

Email: a-panin@yandex.ru

Candidate of Physical and Mathematical Sciences, assistant professor of Faculty of Physics

Leninskie Gory, Moscow, 119991, Russian Federation

Leonid A. Sevastianov

Peoples’ Friendship University of Russia (RUDN University)

Email: sevastianov-la@rudn.ru

Doctor of Physical and Mathematical Sciences, professor of Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation

References