Numerical solution of Cauchy problems with multiple poles of integer order

. We consider Cauchy problem for ordinary differential equation with solution possessing a sequence of multiple poles. We propose the generalized reciprocal function method. It reduces calculation of a multiple pole to retrieval of a simple zero of accordingly chosen function. Advantages of this approach are illustrated by numerical examples. We propose two representative test problems which constitute interest for verification of other numerical methods for problems with poles.


Introduction
There are a number of important applied problems in which the solution has multiple singularities. In such problems, it is required to find a chain of sequentially located singularities. Similar problems are often found in the theory of special functions (elliptic functions, gamma function, etc.).
Numerical methods are widely used to compile tables of special functions [1] and for standard direct calculation programs [2]. Standard schemes (for example, Runge-Kutta schemes) allow one to calculate smooth sections of the solution with good accuracy. However, near the singularity, the error of such schemes increases catastrophically. Direct continuation of the solution beyond the pole, as a rule, is impossible. Therefore, the solution is continued beyond the pole with some artificial techniques. Continuation through a number of poles is an even greater problem and requires the development of special procedures.
The literature describes methods based on the Pade approximation [3]- [5] and on the approximation of the solution by chain fractions [6]. Abramov and Yukhno proposed a special replacement for an unknown function that translates the solution into a non-singular one, see [7] and the bibliography there. However, these methods are applicable only for calculating the transcendental Painleves, for which there is a lot of a priori information. In addition, the coefficients of the Pade approximation are calculated from the coefficients of the Taylor series, and to find the latter, you need to solve the original problem with some difference scheme. The problems that arise along this path are described above.
In [8], we have constructed the reciprocal function method which for the first time allowed to perform highly accurate calculations through a sequence of first-order poles. However, for poles of order > 1, accuracy sharply deteriorated. The reason was as follows: the reciprocal function had a zero of order > 1. Calculation of such zero is an ill-conditioned problem conjuncted with considerable loss of accuracy.
In the present work, we propose the generalized reciprocal function method which overcomes the mentioned difficulty. It provides high accuracy in computation of a sequence of poles with multiplicity > 1 if the differential equation is autonomous.

Method
Let us write down the Cauchy problem for an ordinary differential equation of the first order / = ( , ), (0) = 0 .
(1) Its solution is assumed to have a sequence of poles at points * of integer orders . The orders of the different poles may not be the same. At the same time, we assume that the solution does not have special points of other types.
Let us introduce some fine enough mesh . Let us choose some one-step method of numerical integration. A large number of such methods is given in the monographs [9], [10]. One can detect approach to the nearest pole by rapid increase of the numerical solution . However, this does not allow us to determine the position of the pole with sufficient accuracy, calculate the solution in its vicinity, and continue the solution beyond the pole.
To overcome this difficulty in the case of first-order poles, we proposed the reciprocal function method [8]. Let adjusting parameter > 0 be introduced. If the condition | | > is met, then the calculation proceeds from the function ( ) to the reciprocal function ( ) = [ ( )] −1 . It satisfies the following equation: The initial condition at the transition point is assumed to be = ( ) −1 . Note that such a transition at any mesh node is possible only when using one-step schemes (for example, explicit Runge-Kutta methods).
The pole of the original function ( ) of multiplicity corresponds to the zero of the reciprocal function ( ) of the same multiplicity. For = 1, this is a simple zero, in which the solution of the equation (2) does not present any problem. This is illustrated by examples of numerical calculations in [8]. In this case, the solution is calculated with good accuracy in the vicinity of the pole and continues beyond it. This makes it possible to perform through calculations of the sequence of poles of the first order with good accuracy.
However, for multiplicity > 1, the zero ( ) turns out to be a special point of the equation (2). At this point, the reciprocal function itself and all its derivatives up to the ( − 1)th inclusive turn to zero. Numerical solution through this feature leads to a strong decrease in accuracy and even failure of the calculation. The solution cannot be confidently continued even beyond the first pole.
To overcome this difficulty, we propose to introduce a generalized reciprocal function ( ). Suppose the multiplicity of the nearest pole is known. Then for any , we can put This expression has complex branches. We choose the only real one from them. The generalized reciprocal function satisfies the following differential equation: (4) For it, this zero turns out to be simple, and its calculation does not cause fundamental difficulties. After passing this zero, one can return to calculation of the ( ) function.

Multiplicity determination
Sometimes, from a theoretical study of the Cauchy problem, it is possible to determine a priori the multiplicities of the poles . In general case, one has to find a posteriori in the course of calculation. To do this, we propose the following procedure.
Near the pole, the following relation holds: ( ) ≈ ( * − ) . Then in two adjacent nodes, . This is an over-determined system in unknowns , * , . Excluding and * , one obtains If the resulting is close enough to some integer on several sequential mesh steps, this integer number can be taken as the pole multiplicity. Note that in order to apply the formula (5), the following conditions are necessary (although not sufficient):

Test problem
Let us construct Cauchy problem with the following exact solution: This function has poles at * = 0.5 + . Its derivative equals In the intervals between neighboring poles, the derivative preserves the sign. For odd , the derivative is always positive, and for even , its signs are opposite in neighboring intervals separated by a pole. Therefore, in both cases, the solution (7) and has no special points other than poles.
The equation (8) is of no interest to be considered as Cauchy problem, since the solution is reduced to quadrature calculation. However, in the case of an odd ⩾ 1, the solution (7) and equation (8) can be converted to the form Let us consider (9) as equation in tan and express tan in terms of . Next, we substitute the obtained expression into (10) and obtain autonomous form of the equation.
Practically, explicit relations expressed in elementary functions can be derived only in two cases. The first one corresponding to = 1 is trivial This example was used in [8] as a test for simple pole. The second case with = 3 is non-trivial This test is used in the present work.

Numerical example
Calculation of the test (12) was performed on the segment 0 ⩽ ⩽ 15 containing 5 poles of the third order. The calculation was performed on a sequence of uniform meshes using an explicit Runge-Kutta scheme of the fourth order of accuracy (ERK4). The first grid had a step of = 0.15, the remaining grids were obtained by successive decreasing of all steps by the factor of 2 from mesh to mesh. Figure 1 compares the numerical solution on the first grid (markers) with the exact one (solid line). The vertical lines show the asymptotes of the exact solution. Even with such a large step, one can see good agreement between the numerical solution and the exact one. Figure 2 shows the solution error in mean-squared analogue of the Hausdorff metrics [11] as function of the mesh step. The plot is given in double logarithmic scale. The calculated points lie on a straight line with a slope of −4. This corresponds to the power-law nature of convergence with the theoretical order of accuracy = 4. One can see that the error reaches round-off errors ∼ 10 −14 (which is only 100 times greater than the error of a single rounding equal 10 −16 ) at ≈ 10 5 of grid nodes. This indicates high accuracy and reliability of the method. The position of the poles is determined by interpolation of at two points to the right and left of zero. This procedure is described in [8]. The error of the fifth pole position is shown in figure 2 with triangles.

Difficulties
For illustration, we used a test in which the differential equation was autonomous. However, in applied problems we also have to deal with nonautonomous equations. For such problems, the reliability of numerical methods deteriorates. In the vicinity of the pole, the calculated profile may look like an alternating "saw", which does not allow to determine the position of zero. Let's explain the reason of this phenomenon.
Take, for example, the non-autonomous equation (8). The zero of the grid solution , understood in the sense of different signs of this value in neighboring nodes, does not coincide with the exact pole. At the same time, the sign of the right side (8) that depends only on is determined by the position of the exact pole. The value changes sign when passing through the "mesh" pole, and the right part does so when passing through the exact pole. This lack of synchronization can lead to an unpredictable sign of increment of the value at the next step. The higher is the pole multiplicity the stronger is this effect.
These effects usually reveal on insufficiently fine meshes. To overcome these difficulties, we recommend to choose fine enough mesh. Increasing digit capacity is also a helpful strategy.

Even multiplicity
For a pole of even multiplicity, the Cauchy problem can be non-autonomous only. In fact, near the pole ≈ ( − * ) − , and / ≈ − ( * − ) − −1 . For even , / has different signs on different sides of the pole. Therefore, it cannot be an unambiguous function of ( ). Thus, any problems for an even face all the difficulties that are typical for non-autonomous problems. The ways to overcome them are also indicated above.
The exact solution is as follows: It has poles of the order = 2 at * = /2 + . Calculations were performed using the ERK4 scheme. Figure 3 shows the numerical solution for = and the exact solution (the notation corresponds to figure 1). One can see that the numerical calculation through 5 poles is successful, although the visual difference at the end of the calculation is somewhat greater than for the autonomous problem in figure 1. Figure 4 shows the dependence of the error of the solution itself and the one of the fifth pole position on the mesh step. The notations correspond to A. A. Belov, N. N. Kalitkin, Numerical solution of Cauchy problems … 111 figure 2. One can see that the calculated points are slightly scattered around the average line. This is a manifestation of the difficulties associated with solving non-autonomous problems. However, the average slope of the straight line corresponds to the theoretical order of accuracy = 4, and very high accuracy is achieved on moderate meshes, close to unit rounding errors.

Note
The same exact solution may correspond to different non-autonomous problem formulations. For example, the function (14) is the exact solution of a differential equation However, all attempts to calculate this equation using various quadrature formulas were unsuccessful due to "blow up" of the calculation. Therefore, the tests (12) and (13) constructed here are of value themselves. Solutions have sequences of poles of the specified orders, special points of other types are absent, and the influence of non-autonomy is minimized. These problems are recommended for validation of other methods of calculation through poles.