Robust formula for $N$-point Pad\'e approximant calculation based on Wynn identity

TL;DR

This paper introduces a new criterion based on Wynn's identity for selecting the optimal N-point Padé approximant, enhancing the accuracy and reliability of rational approximations in physics and applied mathematics.

Contribution

The work presents a novel formula and criterion for optimal Padé approximation using Wynn's identity, with practical applications demonstrated in physics and differential equations.

Findings

Wynn's identity provides a criterion for minimal |η| as the optimal Padé approximant

The method improves multipoint Padé approximation accuracy

Application to physics problems like solar corona heating series summation

Abstract

The performed numerical analysis reveals that Wynn's identity for the compass $1/(N-C)+1/(S-C)=1/(W-C)+1/(E-C)=1/\eta$ (here C stands for center, the other letters correspond to the four directions of the compass) gives the long sought criterion, the minimal $|\eta|$, for the choice of the optimal Pad\'e approximant. The work of this method is illustrated by calculation of multipoint Pad\'e approximation by a new formula for calculation of this best rational approximation. The work of the criterion for the calculation of optimal Pad\'e approximant is illustrated by the frequently seen in the theoretical physics problems of calculation of series summation and multipoint Pad\'e approximation used as a predictor for solution of differential equations motivated by the magneto-hydrodynamic problem of heating of solar corona by Alv\'en waves. In such a way, an efficient and valuable control mechanism for $N$-point Pad\'e approximant calculation is proposed. We believe that the suggested method and criterion can be useful for many applied problems in numerous areas not only in physics but in any scientific application where differential equations are solved. The solution of the Cauchy-Jacobi problem is illustrated by a Fortran program. The algorithm is generalized for the case of the first $K$ derivatives at $N$ nodal points.

Equations

$$f(z)\approx r_{l,m}(z)=\frac{\sum _{i=0}^l{a}_i{z}^i}{\sum _{k=0}^m{b}_k{z}^k},b_0=1,$$

$$\frac{1}{{\eta }_{lm}}\equiv \frac{1}{r_{l,m+1}-r_{l,m}}+\frac{1}{r_{l,m-1}-r_{l,m}}=\frac{1}{r_{l+1,m}-r_{l,m}}+\frac{1}{r_{l-1,m}-r_{l,m}},$$

$$\frac{1}{\eta }\equiv \frac{1}{S-C}+\frac{1}{N-C}=\frac{1}{E-C}+\frac{1}{W-C},$$

$$C\equiv r_{l,m},N\equiv r_{l,m-1},E\equiv r_{l+1,m},W\equiv r_{l-1,m},S\equiv r_{l,m+1}.$$

$$f(z)=\underset{l,m\to \infty }{\lim }{r}_{l,m}(z)$$

$$\underset{l,m}{\lim }{\eta }_{lm}\to 0.$$

$${\eta }_{\min }\equiv |{\eta }_{L,M}|={\min }_{l,m}|{\eta }_{l,m}|$$

$$r_{l,0}=\sum _{i=0}^l{a}_i{z}^i=S_l.$$

$$S=C+\frac{1}{\frac{1}{E-C}+\frac{1}{W-C}-\frac{1}{N-C}}$$

$$r_{l,m+1}=r_{l,m}+\frac{1}{\frac{1}{r_{l+1,m}-r_{l,m}}+\frac{1}{r_{l-1,m}-r_{l,m}}-\frac{1}{r_{l,m-1}-r_{l,m}}}$$

$${\eta }_{lm}=\frac{1}{\frac{1}{r_{l+1,m}-r_{l,m}}+\frac{1}{r_{l-1,m}-r_{l,m}}}\equiv \eta =\frac{1}{\frac{1}{E-C}+\frac{1}{W-C}}.$$

$$\ln (1+x)=\underset{n\to \infty }{\lim }{S}_n,S_n(x)\equiv \sum _{k=1}^n(-1{)}^{(k+1)}\frac{x^k}{k}.$$

$$|x-x_1|\le |x-x_2|\le |x-x_2|\le \cdots \le |x-x_N|.$$

$${\tilde{y}}_i=y_i+(x-x_i)y_i^{\prime }+\frac{1}{2}(x-x_i{)}^2{y}_i^{\prime \prime }+\cdots +\frac{1}{K!}(x-x_i{)}^K{y}_i^{(K)}.$$

$$x_0-x$$

$$x_1-x$$

$$x_2-x$$

$$\dots \dots \dots$$

$$x_n-x$$

$$(f_0,f_k)=\frac{(x_0-x)f_k-(x_k-x)f_0}{(x_0-x)-(x_k-x)}.$$

$$(f_0,f_1,f_k)=\frac{(x_1-x)(f_0,f_k)-(x_k-x)(f_0,f_1)}{(x_1-x)-(x_k-x)}.$$

$$|x_0-x|<|x_1-x|<|x_2-x|<\dots |x_n-x|.$$

$$S_0\equiv f_0,S_1\equiv (f_0,f_1),S_2\equiv (f_0,f_1,f_2),\dots ,S_n\equiv (f_0,f_1,f_2,\dots ,f_n).$$

$$S=\underset{n\to \infty }{\lim }{S}_n$$

$$\phi (x;x_0,\tilde{x})=f(\tilde{x})+\frac{f_0-f(\tilde{x})}{x_0-\tilde{x}}(x-\tilde{x})=\frac{(x_0-\tilde{x})f(\tilde{x})-(\tilde{x}-x)f_0}{(x_0-x)-(\tilde{x}-x)},$$

$$(f_0,f_1)=\phi (x;x_0,x_1),(f_0,f_2)=\phi (x;x_0,x_2),(f_0,f_3)=\phi (x;x_0,x_3),,\dots ,(f_0,f_n)=\phi (x;x_0,x_n).$$

$$f_1\leftarrow (f_0,f_1),f_2\leftarrow (f_0,f_2),f_3\leftarrow (f_0,f_3),,\dots ,f_k\leftarrow (f_0,f_n).$$

$$f_k\leftarrow \frac{(x_1-x)f_k-(x_k-x)f_1}{(x_1-x)-(x_k-x)},\text{mbox}fork=2,3,\dots ,n.$$

$$f_{k}\leftarrow\frac{(x_{l}-x)f_{k}-(x_{k}-x)f_{l}}{(x_{l}-x)-(x_{k}-x)}=\frac{\left|\begin{array}[]{cc}(x_{l}-x)&f_{l}\\ (x_{k}-x)&f_{k}\end{array}\right|}{x_{l}-x_{k}},\qquad\mbox{for}\qquad k=l+1,\dots,n;$$

$$(f_0,f_1,f_2,\dots f_n)=(S_0,S_1,S_2,\dots S_n)$$

$$x_0,(x_0,x_1),(x_0,x_1,x_2),(x_0,x_1,x_2,\dots ,x_n)$$

$$\varphi (x;x_0)=f_0+f_0^{\prime }(x-x_0)+\frac{1}{2}{f}_0^{\prime \prime }(x-x_0{)}^2+\cdots +\frac{1}{K!}{f}_0^{(K)}(x-x_0{)}^K,$$

$$f_0=\varphi (x;x_0),f_1=\varphi (x;x_1),f_2=\varphi (x;x_2),,\dots ,f_n=\varphi (x;x_n),$$

Full text

I Introduction and motivation

As the analytical structures take central place in the theoretical and mathematical physics, the Padé approximants of functions become one of the most important problems of applied mathematics. For a general introduction in the problem of Padé approximants, see the well-known monographs Refs. Pade ; Brezinski:91 ; Brezinski:96 . Since the problem was first systematically studied by Frobenius Brezinski:H ; Frobenius:81 and the Padé table introduced by Padé Brezinski:H ; Pade:92 in the 19th century, the contemporary literature is enormous. We will mention only two popular problems, which can be of broad interest for physicists:

1. series summation and 2) extrapolation of functions related to predictor-corrector methods for solution of ordinary differential equations. Searching for numerical recipes in this field, a beginner can read that the downside of the Padé approximation is that it is uncontrolled and the dangers of extrapolation cannot be overemphasized NumRep . And after that our beginner will most probably start searching for a commercial software, where the appropriate analytical formulae are programmed. By the very famous law: whatever can go wrong, will go wrong, in very interesting for the physics problems, the commercial software stops working and the research would continue only in co-authorship with the source code authors.

In 1947 Hannes Alvfén Alfven:47 concluded that the solar corona heating is due to the absorption of magneto-hydrodynamic (now called Alfvén) waves. In order to describe the almost 100-time increase of the plasma temperature in a very narrow region of the solar atmosphere. The programming of the magneto-hydrodynamic equations is relatively easy, but after 30-time temperature increase all available for us commercial computation software stops working. Facing such a problem, we had to develop an own method for a frequently encountered numerical problem (four first order linear differential equations coupled to four first order non-linear differential equations) MHD ; MHD:19 . Our research goal was to develop this method based on calculation of Padé approximants with a mechanism to select the most accurate approximant. This mechanism is precisely the main development of the current research and allows the usage of the described method for solution of differential equations. That is why we suppose that many colleagues are in our situation and we share know-how in calculation of Padé approximants.

Before starting let us recall the definition of multipoint Padé approximationPade . A rational function $f(z)$ which fits given function values $y_{i}$ at various points $x_{i}$, $i=1,\dots,\,N$, i.e. $f(x_{i})=y_{i}$ is called multipoint Padé approximant. The associate problem of interpolation by rational functions is called Cauchy-Jacobi problem. Multipoint Padé approximants are also called: rational interpolants, N-point Padé approximants, Newton-Padé approximants, etc, depending on the context. Addressing to experts in the present paper we present a general robust solution of the Cauchy-Jacobi problem which can attract interest and become a convenient tool in the applied numerical analysis.

II Announcement of the results

For a practical implementation of an analytical result, we have to use a computer with finite accuracy of digital representation of real numbers. The Padé approximants give fast convergence and actually the best rational approximation but they are sensitive with respect to the noise of the discrete representation of real numbers. Roughly speaking every value representing a real number suffers from error of discretization as it is multiplied by a random factor $1+\varepsilon_{\mathrm{mach}}(\mathrm{rnd}-1/2)$, where $\varepsilon_{\mathrm{mach}}$ is the machine epsilon (the biggest positive value for which $1+\varepsilon_{\mathrm{mach}}=1$), and “$\mathrm{rnd}$” is a programming operator generating a random number homogeneously distributed between 0 and 1.

As formulas for Padé approximants picturally speaking are very sensitive to the random noise of truncation, it is necessary to study the influence of the machine error $\varepsilon_{\mathrm{mach}}$ on the final result. In such a way beyond the analytical problem we have also a statistical one which already belongs to another branch of mathematics. Probably for some special cases it is possible to calculate the probability distribution functions (PDF) of the errors of the final result of calculation of Padé approximants, but it deserves to start with descriptive statistics of some well-known examples, which in our opinion illustrate how to calculate and choose the optimal Padé approximant in physical applications. Here the determination of the optimal Padé approximant is crucial, as our beginner has already learned from the numerical recipes book. To experts we propose a robust method which at least empirically works well.

Let $r_{l,m}(z)$ are different Padé approximants of the complex function $f(z)$

[TABLE]

where $l$ is the maximal power of the numerator and $m$ is the maximal power in the polynomial in the denominator. For these approximants in 1966 Wynn discovered the remarkable relation (Wynn:66, , Eqs. (15) and (16))

[TABLE]

later on baptized by Gragg (Gragg:72, , Theorem 5.5) as missing identity of Frobenius. The criterion for choosing of optimal Padé approximant is based on this Wynn relation which can be rewritten also as

[TABLE]

In the analysis of algorithms is is convenient to represent different approximants $r_{l,m}$ in a Padé table. If $C$ is in the Center, the notations of the other elements $N,\,E,\,W,\,S$ follow the direction of the compass: North, East, Wests, and South, also found in the original paper of Wynn Wynn:66 (i.e. we use Central NEWS algorithm). Probably this representation was an idea of Gragg, who has also described it Gragg:72 , because of the footnote in Wynn’s paper on page 266 (if that’s the case, Gragg was a referee of Wynn’s paper). For convergent Padé approximants

[TABLE]

in the Padé table $C,\,N,\,E,\,W,\,S\rightarrow f(z)$, $1/\eta\rightarrow\infty$, and

[TABLE]

In numerical implementations this limit leads to minimal-$|\eta|$ criterion at CNEWS algorithm.

Let us repeat in other words. Often the difference between sequential Padé approximants gives a reasonable evaluation of the error, see Ref. Pade . In case of convergence, we have vanishing differences between the values of different cells of the Padé table $(r_{l,m+1}-r_{l,m})\rightarrow 0$, for $l,m\rightarrow\infty$ and $l/m=\mathrm{const}$. In this case also $\eta_{lm}\rightarrow 0$ and the minimal value of $\eta_{lm}$ gives a reasonable criterion for the minimal error and the choice of the optimal Padé approximant. This theoretical hint we prove on many examples of calculation of Padé approximation and arrive at the conclusion that the long sought criterion for the choice of the optimal Padé approximant is simply a search for the minimal value $|\eta_{l,m}|$, i.e.

[TABLE]

and $f(z)\approx r_{L,M}$. It is technologically and aesthetically attracting that the Wynn identity gives simultaneously

1. explicit method for the calculation of Padé approximants in $\varepsilon$-table and
2. method for the evaluation of the error; $\lim\eta_{\mathrm{min}}=0$ is actually a criterion for convergence of Padé approximation if we consider real numbers. Our criterion is an alternative of the singular value decomposition (SVD) method described in great detail in Refs. Gonnet:13 ; Beckermann:13 , see also Refs. Karlsson:76 ; Gonchar:78 ; Gonchar:89 ; Levis:18 ; Yattselev:18 . Our minimal $\eta$ criterion and the tolerance level of SVD are perhaps different implementations of one and the same idea. For both methods the rounding error is something external for the theory of Padé approximation. There are no theoretical justifications which method is better and the numerical experiment and the descriptive statistics are just the first steps in the practical realizations on the robust Padé approximants.

The purpose of the present work is to demonstrate how this criterion works and can be used in applications. In the next section we illustrate in detail this criterion on the problem of of calculation of limes of a numerical series.

III The new algorithm for determination of the optimal Padé approximant

Let us have a numerical sequence $\{S_{0},S_{1},S_{2},\dots,S_{N}\}$ and we need to calculate the limit $S=\lim\limits_{l\rightarrow\infty}S_{l}$. The well-known method is to initialize $r_{l,0}=S_{l}$, for $l=0,\dots,N$, and also $r_{l,-1}=\infty$ for $l=0,\dots,N+1$. In other words we treat every term of the initial series as polynomial approximation of some function

[TABLE]

The method can work even if this sequence is convergent, and those are the most interesting cases for applications. Then we calculate in the “south direction” the corresponding Padé approximants

[TABLE]

or

[TABLE]

and simultaneously calculate the empirical error

[TABLE]

The minimal absolute value in the $\eta$-table $\eta_{\mathrm{min}}$ is our criterion for the determination of the optimal Padé approximant. According to the best we know, this minimal-$|\eta|$ criterion has never been implemented in the numerical recipes so far (authors will very much appreciate any information about this).

For the programming task, we have to calculate all the values, for which division is possible. In the next section we thoroughly describe 2 technical examples.

IV Two easy pieces

IV.1 Calculation of divergent series $\ln(1+x)$

Perhaps the most famous example summation of divergent series by Padé approximants is the calculation of the divergent Taylor series $\ln(1+x)$ beyond the radius of convergence $|x|<1$

[TABLE]

The application of Eq. (10) and Eq. (11) gives the optimal Padé approximant for every positive $x$ shown in Fig. 1.

Even for $x\approx 20$ for the calculation of the series $\ln(1+x)$ a pixel accuracy is present. In many articles results of numerical calculations are illustrated by the figures presented by computer graphics with million pixels. For such a figures pixel accuracy means maximal error in which the approximation is undistinguished from the exact result when graphically presented on screen. For many applied studies this is an acceptable beginning. For larger values of $x$ the dots representing the calculation evaporate from the line representing the exact value. The most important detail is how and when the method stops working and when the resources of the numerical accuracy are exhausted. The real $\varepsilon_{\mathrm{real}}$ and empirical $\varepsilon_{\mathrm{emp}}\equiv\eta_{\mathrm{min}}$ of the calculation are shown in Fig. 2.

The empirical error $\varepsilon_{\mathrm{emp}}\equiv\eta_{\mathrm{min}}$ shows saturation and is a reliable indicator for the calculation accuracy. Statistically analysis of the correlation between both errors from Fig. 2 shows very strong dependence between them in Fig. 3.

Now we can safely tell how accurate our Padé approximation is and even how far out in $x$ it can be extended. This dependence reveals that the long sought criterion for empirical evaluation of the accuracy of the Padé approximants calculated by $\varepsilon$-algorithm has already been found.

IV.2 Extrapolation of a sine arch from a preceding one

The second technical example we present is the problem of extrapolation of function which we illustrate in the case of the $\sin(x)$ function. We take $N$ equidistant interpolation points from one arch of the $\sin(x)$ function and extrapolate the next arch and even beyond. We use the well-known Aitken interpolation method Aitken (1932); Bronshtein & Semedyayev (1955); Abramowitz & Stegun (1972); Korn2 to order the interpolation points within the arch and the point to be extrapolated. In more details: we calculate $y_{i}=\sin(x_{i})$ for all $N$ points $i=1,\dots,N.$ The digital noise is minimized if the points are ordered with the proximity with the point of interpolation, i.e.

[TABLE]

Then using sequentially $l$ interpolation points in the Aitken scheme we calculate the corresponding polynomial approximation $S_{l}$ for $l=1,\dots,N.$ Next we apply the CNEWS method and use the minimal-$|\eta|$ for the best Padé approximant. In this manner we have a numerical sequence that we give to the $\varepsilon$-algorithm to calculate its limit. The described calculation for the $\sin(x)$ function is shown in Fig. 4.

The preceding arch is in the interval $(-\pi,0)$ and contains 21 interpolation points. Using these points, we extrapolate one arch with 2000 points in the interval $(0,\pi)$ and continue the extrapolation in an attempt to obtain a second arch with the same number of points in the next interval $(\pi,2\pi)$. Repeating, every extrapolated point is calculated independently using the Aitken method for interpolation followed by calculation of Padé approximants by the Wynn identity and the optimal Padé approximant is chosen by the advocated in the present paper $\eta_{\mathrm{min}}$ criterion. The deviation of the points from the real function represented with the line in Fig. 4 shows the limit of applicability of the Aitken-Wynn extrapolation algorithm. Detailed error estimates of the extrapolation are shown in Fig. 5.

The $\eta_{\mathrm{min}}$ criterion shown in Fig. 5 gives the reliable order estimation of the error. The similar behaviour of both errors shows that our criterion is a reliable method for error estimation. Furthermore, a statistical analysis of the correlation between the real $\varepsilon_{\mathrm{real}}$ and the empirical error $\varepsilon_{\mathrm{emp}}$ shows very strong dependence in Fig. 6.

The literature on the problem of extrapolation is enormous, however we have not found a discussion of the problem of the mathematical statistics. If we know the values of the function with some accuracy, how to choose the best Padé approximant among the many, some of which have serious noise of rounding and others are even uncontrolled NumRep . There is a huge literature about the formulae, but a convenient criterion for the choice of optimal Padé approximant is not discussed. We give our illustrative example not as an exercise of programming, but as an illustration that we have a good working criterion for many practical cases. As a rule, articles on Padé approximation are illustrated by examples describing how some method is good working, but the instructive analysis reveals where the method stops working and what is most important whether some criterion warns that calculation resource has been exhausted. From a practical point of view, exhaustion of resources for extrapolation is demonstrated by evaporation of the extrapolated points from the analytical curve, and one of the achievements of the present work is that order of the magnitude of this deviation is reliably indicated by our criterion based on the Wynn identity.

V Discussion and conclusions

In spite that the literature for Padé approximants is enormous we were unable to find alternative formulas, methods, algorithms and programs to compare the work of our method. One can put in the agenda for the development of the numerical analysis a simple test: which method for extrapolation of functions from tabulated values in $N$ points gives the best extrapolation of a function far from the interpolation interval. This juxtaposition will give the final verdict which method is appropriate to be implemented in commercial software like Mathematica and Maple.

The new result of the implementation of calculation of Padé approximants and their application is the modulus minimization $\eta_{\mathrm{min}}$ of Wynn’s $\eta_{lm}$ as reliable empirical criterion of the error.

The preformed analysis of several simple examples has revealed that for practical implementation of Padé approximants we can use the empirical $\eta_{\mathrm{min}}$ error extracted from the Wynn’s identity. In the agenda the statistical problem of calculation of PDF of the Padé approximants has already been set. The comparison of descriptive statistics data for the PDF of errors of calculation of Padé approximants by different criteria will give what the answer what general recommendation as a numerical recipe have to be given to users not willing to understand how.

In short, the practical implementation of Padé approximants can reach one order of magnitude more applications in theoretical physics and applied mathematics. More than half a century after its discovery, the $\varepsilon$-algorithm has not been included for calculation of divergent series with convergent Padé approximants and for extrapolation of functions in commercial software. Now the time for this inclusion has come, the herein implemented control mechanism has rendered this mission possible.

Last but not least, the suggested criterion Eq. (11) is applicable in solution of differential equations, numerical analytical continuation, perturbation, series summation and other analogous problems of theoretical physics.

Let us summarize our results:

1. we have suggested a method for empirical choice of optimal Padé approximant which woks for all examples we have met in the literature and seems to have universal applicability,
2. we derived a new formula for the $N$-point Padé approximant, for a function tabulated in $N$-points. This formula based on Eitken interpolation formula and Wynn identity is maximally robust with respect of truncation errors in the numerical calculations and can be used even for extrapolation. We have developed this method in order to obtain one predictor method for predictor-corrector methods for solving of one irrelevant for the present paper magnetohydrodynamic problem, the programs where our formulae are implemented are given as appendices of the arXiv version of the present paper arxiv and early and more detailed version has already been published in conference proceedings AIP_Pade . The figures in the paper are calculated by those programs. Reproducing the figures is the criterion that our formula can be used directly in numerical calculations. Up to now a non-specialist knows that extrapolation and summation of divergent series is a forbidden procedure used only by elite mathematicians. We have made a popularization addressed to user of applied mathematics working for the industry.
3. The described method can be used as a predictor method for solving ordinary differential equations $\mathrm{d}y/\mathrm{d}x=\mathcal{F}(x,y)$. The described N-point Padé approximation can be used as a predictor using former calculated points and their first derivatives as interpolation points to calculate the predicted value of the function $y_{\mathrm{pred}}$ at the new point $x_{\mathrm{new}}$ of the argument. Then we calculate the derivative at the new point $\mathcal{F}(x_{\mathrm{new}},y_{\mathrm{pred}})$. We consider this possibility as a new method for solution of ordinary differential equations. The method can be additionally improved if in the ($x$-$y$) plane we perform a rotation, and the new abscissa of the interpolation to be taken along the local tangent of the curve. We started with the problem of calculation of the temperature profile of the solar corona and velocity of solar windMHD:19 but this method can be applied for many other problems for which standard software is not (good) working.
4. If we know the first $K$ derivatives of the function using of Taylor expansion we can calculate

[TABLE]

Then we can use the Taylor approximants $\tilde{y}_{i}$ as interpolation points $x_{i}$ in the already described method. At the basic points the interpolated rational function will have exact coincidence not only for the function but also for its first $K$ derivatives. The main idea of the algorithm is the same: using a node point $x_{i}$ we calculate the interpolated value $\tilde{y}(x)$ at the argument of interpolation, and then use this calculated value $\tilde{y}(x;x_{i})$ at this point as a new value in the node point $y_{i}=f(x_{i})\leftarrow\tilde{y}(x;x_{i})$. If we need a formal proof of convergence it is necessary in the right side of Eq. (14) to put insert the evaluation of error $h^{(K+1)}M_{(K+1)}/(K+1)!$, where $h$ is the maximal distance between interpolating an node argument and $M_{(K+1)}$ is the supremum of the modulus of the $(K+1)$ derivative. This estimation of the error reveal why the analytical approximation is not applicable to function $|x|$ and the trajectory of a Brownial particle.

Finally we have to remark. There is no general criterion for the convergence of the Padé approximation, and even to choose the optimal approximant when the method is convergent, that is why there are no theorems which could be cited. As all known existing methods are given without any empirical prescription for the choice of the optimal approximant the considered minimal-$|\eta|$ criterion of the CNEWS method can not be compared with other methods. The area of the applicability of the suggested method is however well defined. For many applied problems it is known that solution exists, while for the analytical problems the solution is analytical. In many problems in the physics, for example, there are alternative method to check whether the searched solution is numerically obtained in spite that there will be no alternative methods. Contra-examples of mathematical series for which the Padé method does not work can be easily constructed, but we are unaware of a real problem for which the solution is experimentally measured but the series can not be summarized by the exciting combinations of Euler-McLaurin summation and Padé approximation. That is why we believe that suggested method, criterion and software can be useful for many applied problems in numerous areas not only in physics but in any scientific application where differential equations are solved.

VI Acknowledgments

The authors are thankful to Evgeni Penev for the collaboration in early stages in the present research, and for writing the $\varepsilon$-algorithm in Fortran-90, Michail Mishonov gave a version in Java and Nedeltcho Zahariev and Zlatan Dimitrov in C-language. The authors are also grateful to Claude Brezinski for pointing out and recommendation of Refs. Gonnet:13 ; Beckermann:13 , to Ognyan Kounchev for Refs. Karlsson:76 ; Gonchar:78 ; Gonchar:89 ; Levis:18 ; Yattselev:18 , to Emil Horozov, Radostina Kamburova, Dantchi Koulova, Peter Kenderov and Hassan Chamati for considerations, comments and stimulating discussions. This work is partially supported by grant KP-06-N38/6 from 5.12.2019 of the Bulgarian National Science Fund and National program “Young scientists and postdoctoral researchers” approved by DCM 577, 17.08.2018.

Appendix A $\varepsilon$-algorithm source code in Fortran-90

This section includes the $\varepsilon$-algorthm written in Fortran-90 subroutine and used for the calculations presented in the figures. This Fortran program is our understaning of Cauhi-Jakobi problem an if the time of floating point operation $\tau_{{}_{\mathrm{FPO}}}\rightarrow 0$ and also the machine epsilon $\epsilon_{\mathrm{mach}}\rightarrow 0$ we obtain the analytical result. The program can be converted in a constructive proof.

Appendix B Aitken algorithm

In order to explain the novelty of the present achievement in the beginning we will recall the well known Aitken algorithm. Later on using this basis we will explain the general polynomial exponential. We follow the well-known reference book by Bronstein and SemendyayevBronshtein & Semedyayev (1955). If we need to calculate value of the function $\varphi_{n}(x)$ at some fixed argument $x$, one can use following schematics (“crest by crest”) which is especially convenient for computer programming

[TABLE]

Every symbol $(f_{0},\,f_{1},\dots,f_{k})$ denotes the value in the point $x$ of the interpolating polynomial function, build on nodes $x_{0},\,x_{1},\dots,x_{k}$, i.e. $f_{0}=f(x_{0})$, $f_{1}=f(x_{1})$, …, $f_{n}=f(x_{n}).$ Those values calculates column by column, in the following manner. Numbers of the column $(f_{0},\,f_{k})$ are obtained by

[TABLE]

Every following column is calculating by the former one at the same schematics, for example

[TABLE]

The rounding error is minimized if we reorder the points as

[TABLE]

In such a manner we obtain the series

[TABLE]

In order to obtain an approximation of limit of this sequence

[TABLE]

we use the Wynn algorithm CNEWS algorithm with minimal-$|\eta|$ criterion, which works even if the sequence $(S_{0},\,S_{1},\,S_{2},\dots)$ is divergent. It is only necessary the problem which we solve to be analytical one born from a real physical problem.

Let us repeat the Aitken formula as an algorithm. Having fixed node point $x_{0}$ we can calculate the the linear approximation using a second node point $\tilde{x}$

[TABLE]

where node arguments $x_{0}$ and $\tilde{x}$ are parameters of the linear function. In such a way omittin in the notations the argument $x$ we can calculate the sequence

[TABLE]

In order to rewrite this conveniently for programming algorithm, we can use the same notations for the changing new values

[TABLE]

Then for the next column of the Aitken table according Eq. (24) using so re-denoted values we calculate

[TABLE]

The formula for a general column $l$ is the same

[TABLE]

i.e. linearly interpolated value of the function at point $x$ is ascribed as the new function value at the node point $x_{k}$. In other words polynomial interpolation is reduces to sequential linear interpolations This procedure repeats until $l=n-1$. In such a way we program the series identically by the so defined polynomial interpolations

[TABLE]

using sequentially

[TABLE]

set of nodes. According to the general theory of Padé approximants the Wynn CNEWS algorithm converts polynomial approximation into Padé ones, and for $N=n-1$ we obtain $N$-point Padé approximant.

If the first $K$ points coincide $x_{0}\rightarrow x_{1},\dots,\rightarrow x_{K}$ the corresponding polynomial approximant is just the Taylor series

[TABLE]

where $f_{0}^{\prime},\,f_{0}^{\prime\prime},\,\dots,\,f_{0}^{(K)}$ are the first $K$ derivatives of the approximated function, and $x_{0}$ is the parameter of this series expansion. Re-denoting

[TABLE]

we can apply the regular Aitken procedure using the new Taylor expanded functional values. In such a way finally we have solved the generalized Cauchy-Jacobi problem: to obtain a rational approximation of a function passing through the N-points $x_{1},\,x_{2},\dots,x_{N}$ with fixed values of its first derivatives $f_{i}^{\prime},\,f_{i}^{\prime\prime},\,\dots,\,f_{i}^{(K)}$ at every at those points $i=1,\dots,N.$ We gave only an idea for proof and a Fortran program illustrate our understanding.

We will use this achievement for many applied problems. For example, for solution of ordinary differential equation

[TABLE]

Using former N-points with first derivative Taylor approximation $K=1$ we can calculate the functional value $y(x)$ using the explained method of Padé extrapolation and then calculate the first derivative at this point explicitly $f^{\prime}(x)=\mathcal{F}(x,y(x))$. Renumbering later $x_{n+1}\leftarrow x_{n}$, $x_{n-1}\leftarrow x_{n},$…$x_{2}\leftarrow x_{1},$ $x_{1}=x_{0}$ and $x_{0}\leftarrow x$ we can continue the procedure. The number of former points is determined by the indices of the optimal Padé approximants according to the minimal-$|\eta|$ criterion. For the step $h$ for the next point $x=x_{0}+h$ one can apply many different criteria. To our knowledge such combination of methods for solution of ordinary differential equation without fixed order is a new one. We consider that content of this Appendix is implicitly given in the body of the article we give this clarification only for the beginner users.