CC BY
13
Discrete & Continuous Models
#& Applied Computational Science 2022, 30 (2) 105-114
ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph
Research article
UDC 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2022-30-2-105-114
Numerical solution of Cauchy problems with multiple
poles of integer order
Aleksandr A. Belov1'2, Nikolay N. Kalitkin3
1 Lomonosov Moscow State University, 1, bld. 2, Leninskie Gory, Moscow, 119991, Russian Federation 2 Peoples' Friendship University of Russia (RUDN University), 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation 3 Keldysh Institute of Applied Mathematics RAS, 4 A, Miusskaia Sq., Moscow, 125047, Russian Federation
(received: March 12, 2022; revised: April 18, 2022; accepted: April 19, 2022)
Abstract. 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.
Key words and phrases: Cauchy problem, singularities, continuation through a pole, multiple poles
1. 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
© Belov A. A., Kalitkin N.N., 2022
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4-0/
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 k > 1, accuracy sharply deteriorated. The reason was as follows: the reciprocal function had a zero of order k > 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 k > 1 if the differential equation is autonomous.
2. Generalized reciprocal function 2.1. Method
Let us write down the Cauchy problem for an ordinary differential equation of the first order
du/dt = f(u,t), u(0) = u°. (1)
Its solution is assumed to have a sequence of poles at points t*m of integer orders km. 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 tn. 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 un. 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 U > 0 be introduced. If the condition lun| > U is met, then the calculation proceeds from the function u(t) to the reciprocal function v(t) = [u(t)]-1. It satisfies the following equation:
dv/dt = -v2 f(v-1 ,t). (2)
The initial condition at the transition point is assumed to be vn = (un)-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 u(t) of multiplicity k corresponds to the zero of the reciprocal function v(t) of the same multiplicity. For k = 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 k > 1, the zero v(t) 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 (k — 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 w(t). Suppose the multiplicity of the nearest pole k is known. Then for any k, we can put
w(t) = [v(t)]1/k. (3)
This expression has k complex branches. We choose the only real one from them. The generalized reciprocal function satisfies the following differential equation:
dw/dt = —k-1 w1+k f(w-k, t). (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 u(t) function.
2.2. Multiplicity determination
Sometimes, from a theoretical study of the Cauchy problem, it is possible to determine a priori the multiplicities of the poles km. In general case, one has to find km a posteriori in the course of calculation. To do this, we propose the following procedure.
Near the pole, the following relation holds: v(t) & A(t* — t)k. Then in two adjacent n°de^ vn & A(t* — tn)k, vn+i & A(t* — tn+i)k, fn =, fn+i. This is an over-determined system in unknowns A, t*, k. Excluding A and t*, one obtains
ln(fn f-h )
k
1—
lnK Vn+1).
(5)
If the resulting k 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):
vnvn+1 > 0, fn fn+1 > 0, Vn fn < 0, \Vn \ > lvn+11 (6)
2.3. Test problem Let us construct Cauchy problem with the following exact solution:
u(t) = sin t cos-k t. (7)
This function has poles at t*m = 0.5^ + nm. Its derivative equals
du/dt = cos1-k t + k sin2 t cos-fc-11. (8)
In the intervals between neighboring poles, the derivative preserves the sign. For odd k, the derivative is always positive, and for even k, 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 k > 1, the solution (7) and equation (8) can be converted to the form
u(t) = tan t(1 + tan2 t)(k-1)/2, (9)
du/dt = (1 + k tan2t)(1 + tan2 t)(k-1)/2. (10)
Let us consider (9) as equation in tant and express tant in terms of u. 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 k = 1 is trivial
u(t) = tan t, du/dt = 1 + u2. (11)
This example was used in [8] as a test for simple pole.
The second case with k = 3 is non-trivial
u(t) = tan t + tan31, du/dt =(l+i(u)2 )(l + 3i(u)2), (12)
^(u) = —2 • 3-0'5sign(u) sinh ^(u), <p(u) = 3-1 arsinh(0.5 • 315 |m|).
This test is used in the present work.
2.4. Numerical example
Calculation of the test (12) was performed on the segment 0 < t < 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 t = 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 p = 4. One can see that the error reaches round-off errors ~ 10-14 (which is only 100 times greater than the error of
10
0 5 10 15
t
Figure 1. Calculation of the test (12) with step t = 0.15 using the ERK4 scheme. Comparison of the numerical solution on the first grid (markers) with the exact one (solid
line)
a single rounding equal 10-16) at N « 105 of grid nodes. This indicates high accuracy and reliability of the method. The position of the poles is determined by interpolation of wn 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.
Figure 2. Dependence of the error of the solution and the fifth pole position (triangles)
on the step size for the test (12)
3. Non-autonomous problems 3.1. Difficulties
For illustration, we used a test in which the differential equation was autonomous. However, in applied problems we also have to deal with non-autonomous equations. For such problems, the reliability of numerical methods deteriorates. In the vicinity of the pole, the calculated profile vn 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 vn, 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 t is determined by the position of the exact pole. The value vn 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 vn 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.
3.2. Even multiplicity
For a pole of even multiplicity, the Cauchy problem can be non-autonomous only. In fact, near the pole u & A(t — t*)-k, and du/dt & —kA(t* — t)-k-1. For even k, du/dt has different signs on different sides of the pole. Therefore, it cannot be an unambiguous function of f(u). Thus, any problems for an even k face all the difficulties that are typical for non-autonomous problems. The ways to overcome them are also indicated above.
3.3. Example
Consider the following non-autonomous problem
du/dt = (0.5+ ^0.25 +u2 + 2u2) cost, «(0) = 0. (13)
The exact solution is as follows:
u(t) = sin t cos-2 t. (14)
It has poles of the order k = 2 at t*m = n/2 + nm.
Calculations were performed using the ERK4 scheme. Figure 3 shows the numerical solution for t = 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
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 p = 4, and very high accuracy is achieved on moderate meshes, close to unit rounding errors.
Figure 3. Calculation of the test (13) with step t = 0.15 using the ERK4 scheme. The notations correspond to figure 1
1 1.5 2 2.5 3 3.5 -lg 7
Figure 4. Dependence of the error of the solution and the fifth pole position on the step size for the test (13). The notations correspond to figure 2
3.4. 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
du/dt = cos-11 + 2 sin2 t cos-3 t. (15)
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.
Acknowledgments
This work was supported by grant MK-3630.2021.1.1.
References
[1] L. F. Janke E. Emde F., Taffeln horere Functionen. B.G. Teubbner Verlagsgesellschaft, Stuttgart, 1960.
[2] NIST digital library of mathematical functions, https://dlmf.nist.gov.
[3] C. F. Corliss, "Integrating ODE's in the complex plane - pole vaulting", Mathematics of Computation, vol. 35, pp. 1181-1189, 1980. DOI: 10. 1090/S0025-5718-1980-0583495-8.
[4] B. Fornberg and J. A. C. Weideman, "A numerical methodology for the Painleve equations", Journal of Computational Physics, vol. 230, pp. 5957-5973, 2011. DOI: 10.1016/j.jcp.2011.04.007.
[5] M. Fasondini, B. Fornberg, and J. A. C. Weideman, "Methods for the computation of the multivalued Painleve transcendents on their Riemann surfaces", Journal of Computational Physics, vol. 344, pp. 3650, 2017. DOI: 10.1016/j.jcp.2017.04.071.
[6] I. M. Willers, "A new integration algorithm for ordinary differential equations based on continued fraction approximations", Communications of the ACM, vol. 17, pp. 504-508, 1974. DOI: 10.1145/361147.361150.
[7] A. A. Abramov and L. F. Yukhno, "A method for calculating the Painleve transcendents", Applied Numerical Mathematics, vol. 93, pp. 262-267, 2015. DOI: 10.1016/j.apnum.2014.05.002.
[8] A. A. Belov and N. N. Kalitkin, "Reciprocal function method for Cauchy problems with first-order poles", Doklady Mathematics, vol. 101, no. 2, pp. 165-168, 2020. DOI: 10.1134/S1064562420020040.
[9] E. Hairer, G. Wanner, and S. P. N0rsett, "Solving ordinary differential equations: I. Nonstiff problems", in Springer Series in Computational Mathematics. Berlin: Springer, 1993, vol. 8. DOI: 10.1007/978-3-54078862-1.
[10] E. Hairer and G. Wanner, "Solving ordinary differential equations: II. Stiff and differential-algebraic problems", in Springer Series in Computational Mathematics. Berlin: Springer, 1996, vol. 14. DOI: 10.1007/9783-642-05221-7.
[11] A. A. Belov and N. N. Kalitkin, "Efficient numerical integration methods for the Cauchy problem for stiff systems of ordinary differential equations", Differential equations, vol. 55, no. 7, pp. 871-883, 2019. DOI: 10.1134/S0012266119070012.
For citation:
A. A. Belov, N. N. Kalitkin, Numerical solution of Cauchy problems with multiple poles of integer order, Discrete and Continuous Models and Applied Computational Science 30 (2) (2022) 105-114. DOI: 10.22363/2658-46702022-30-2-105-114.
Information about the authors:
Belov, Aleksandr A. — Candidate of Physical and Mathematical Sciences, Associate Professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University); Researcher of Faculty of Physics, Lomonosov Moscow State University (e-mail: aa.belov@physics.msu.ru, phone: +7(495)9393310, ORCID: https://orcid.org/0000-0002-0918-9263, ResearcherID: Q-5064-2016, Scopus Author ID: 57191950560)
Kalitkin, Nikolay N. — Doctor of Physical and Mathematical Sciences,
Professor, Corresponding member of the RAS, head of department, Keldysh Institute of Applied Mathematics RAS (e-mail: kalitkin@imamod.ru, phone: +7(499)2509726, ORCID: https://orcid.org/0000-0002-0861-1792)
УДК 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2022-30-2-105-114
Численное решение задач Коши со множественными полюсами целого порядка
А. А. Белов1,2, Н. Н. Калиткин3
1 Московский государственный университет им. М. В. Ломоносова, Ленинские горы, д. 1, стр. 2, Москва, 119991, Россия 2 Российский университет дружбы народов, ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия 3 Институт прикладной математики им. М. В. Келдыша РАН, Миусская пл., д. 4-А, Москва, 125047, Россия
Аннотация. Рассмотрена задачи Коши для обыкновенного дифференциального уравнения с решением, обладающим последовательностью кратных полюсов целого порядка. Предложен обобщённый метод обратной функции, который сводит вычисление кратного полюса к расчёту простого нуля соответственно выбранной функции. Преимущества такого подхода проиллюстрированы на численных примерах. Предложены сложные тестовые задачи, которые представляют интерес для проверки других численных методов для задач с полюсами.
Ключевые слова: задача Коши, сингулярности, продолжение за полюс, кратные полюсы
