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

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 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.


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.

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].
is always algebraic, i.e., in the vicinity of the singular point such a solution can be expanded in Puiseux series

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 A remarkable fact is that after the moving singularity it equals where is the order of algebraic singularity. The complex first-order Rosenbrock scheme for autonomous system of differential equations is written as 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 Then the behavior of the derivative has the form Let us denote this "effective" right-hand side as ( ). In this case The Rosenbrock schemes have the form where is the time step of the grid.
In application to our right-hand side we get or, for the complex Rosenbrock scheme, which is of interest for us, i.e., for = (1 + )/2, from which it follows that the fixed point of the scheme is ) .
Now let us apply the Richardson estimate of the effective accuracy order, which has the form where is the grid densening factor (we used 3 grids with the step , / and / 2 ). From (8) and (9) we get eff = log which demonstrates that the Richardson estimate of the effective accuracy order using the Rosenbrock complex scheme yields the order of the algebraic singularity.

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.

CROS
Function cros(F,ics) for the Cauchy problem described in terms of its argumentṡ= 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 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
on the segment 0 < < 3. Let us rewrite it in the form (10) When specifying the initial conditions let us use decimal fractions, which will be perceived by the system as numbers from ℝ. For plotting (Figure 1) it is possible to use a standard procedure:

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 2 2 points. Then the expression 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 probleṁ 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 from example 10 on the segment 0 < < 3. It is well seen that the singularity is somewhere near = 2.2 and has an order of −2, see also Figure 2.

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.
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.
{ ″ = 2 3 + + 1, 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.
We see that has a pole of the first order.

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 This system is Hamiltonian with and, what is most essential for the subsequent analysis, it can be written in the form̂= wherê( This representation was guessed and can be checked by direct substitution. Equation (15)  The only solution to this Cauchy problem iŝ̂ * = , therefore,̂is a unitary matrix. It remains to note that * ̂̂= −̂ * ̂̂̂+̂ * [,]̂+̂ * ̂̂̂= 0, i.e., the eigenvalues of the matriceŝ( ( ), ( )) 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̂( from where using the Taylor formula we immediately havê * ̂̂=̂| =0 +| =0 . Therefore, general solution 1 = 1 ( ), … is a set of eigenvalues of the matrix 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.

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 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. 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 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. 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.