Padé Approximants, Their Properties, and Applications to Hydrodynamic Problems

: This paper is devoted to an overview of the basic properties of the Padé transformation and its generalizations. The merits and limitations of the described approaches are discussed. Par-ticular attention is paid to the application of Padé approximants in the mechanics of liquids and gases. One of the disadvantages of asymptotic methods is that the standard ansatz in the form of a power series in some parameter usually does not reflect the symmetry of the original problem. The search for asymptotic ansatzes that adequately take into account this symmetry has become one of the most important problems of asymptotic analysis. The most developed technique from this point of view is the Padé approximation.


Introduction
"Taylor series per se are almost useless. The real power of power series is in combination with other ideas. … Power series in isolation are almost never useful to solve boundary-value problems because it is only occasionally that the radius of convergence is sufficiently large to embrace both boundaries. Nevertheless, when combined with another idea, Padé approximants (PAs), power series become an effective tool for ordinary differential equation boundary-value problems" [1].
The use of PAs allows us in many cases to overcome the local nature of the obtained solutions inherent in the methods of series and the method of perturbations [1][2][3][4][5][6].
In this paper, we first provide an overview of PAs, and then describe the advantages of using PAs for solving specific hydrodynamic problems.
Problems of the flow of incompressible viscous liquid in rotating channels are considered. More exactly, we investigate thin domains near horizontal walls perpendicular to the axis of rotation. At high speeds of channel rotation, the so-called Ekman layers are formed in these domains. As a result, the flow region in the rotating channel is divided into boundary layers with an inhomogeneous distribution of velocity components and a core with a uniform velocity distribution.
The appearance of Ekman layers leads to the formation of two flows in the rotating channel. The main flow with a velocity w is formed in the direction of the Oz axis, and the secondary flow u is formed in the direction of the Ox axis.
In this paper, we also aim to study the structure of Ekman boundary layers by asymptotic methods using PAs. This paper is organized as follows. We present the definition of PAs and describe their general properties in Section 2. Sections 3 and 4 describe the merits and drawbacks of PAs and their generalizations. Section 5 is devoted to the application of PAs in nonlinear mechanics. The use of PAs in boundary layer problems is described in Section 6. The application of the PA to the problem of rotating fluid of the Ekman layer type is analyzed in Section 7. Finally, Section 8 presents our concluding remarks.

Definition and Properties
The PA is the "best" approximation of a function by a rational function of a given order [7]. PAs are rational functions, with a denominator that does not vanish at zero, and whose series expansion matches a given series as far as possible. For power series Usually, one normalizes by a coefficient condition such as 0 1 = b (4) Another normalization is also possible [8] 1 = b (5) where ... is the Euclidean norm of the vector.
Under normalization (4), equating the coefficients at the same degrees of ε , we obtain the system of linear algebraic equations where The terms of the first row of the Padé For a description of different types of continued fractions, algorithms for their construction, and application areas, see [9]. The advantage of continued fractions is the significant reduction in the number of "long" arithmetic operations (multiplications and divisions).
Many methods, such as the Aitken, Shanks, eta, and epsilon methods, are mathematically equivalent to the evaluation of certain Padé approximants but have traditionally been formulated algorithmically by fast recurrence relations related to continued fractions in the interests of speed [10,11].
Let us note some properties of PAs [6]: 1. If a PA exists for the chosen m and n , it is unique.
2. If the sequence PA converges to a given function, the roots of its denominator tend to the poles of the function. 3. A PA realizes a meromorphic continuation of a function given by a power series. 4. The PA of the inverse function is equal to the PA inversion of the function itself. This property is called duality and is more strictly formulated as follows. If provided that one of these approximations exists. 5. Diagonal PAs are invariant under linear fractional transformations of the argument.
Let the function be given by its expansion (1). Consider a linear fractional transfor- ε ε that preserves the origin and a function provided that one of these approximations exists. 6. Diagonal PAs are invariant under linear fractional transformations of functions. Let function (1) be given. We put

Merits and Drawbacks of PAs
The mathematical theory of PAs has been developed quite deeply [6,11,12]. However, rigorous mathematical theorems are focused on the widest possible classes of functions, including those that are pathological from the point of view of practical applications. Therefore, the general assessments obtained along this path are often too pessimistic. For example, the very general theorem of Gonchar [12] states: if none of the diagonal PAs and for any domain that is the union of these circles. A significant drawback for practice is the need to check all diagonal PAs. If in a disc of radius R only some sequence of diagonal PAs has no poles, then its uniform convergence to the original function holomorphic in the given disc is guaranteed only for [13].
For practical applications, the most acceptable means of studying PA properties are numerical experiments. In particular, on their basis, it can be argued with a high degree of reliability that in the overwhelming majority of practically important cases The epsilon method, as shown in [14], is the most promising from the point of view of summing diverging or poorly converging series.
As a rule, the rate of convergence of rational approximations significantly exceeds the rate of convergence of polynomial approximations. For example, the function z e in the circle of convergence is approximated by a diagonal PA of the n-th degree 4n times better than by an algebraic polynomial of degree 2n. This is even more noticeable for functions of limited smoothness. E.g., the function ε on the interval [ 1,1] − cannot be approximated by an algebraic n-polynomial so that the order of approximation is better than 1 n. The PA gives the rate of convergence ~exp( 2 ) − n [15].
As a result of numerical experiments, it was found that diagonal and close to diagonal PA sequences often have the autocorrection property [16][17][18][19]. Its essence is as follows.
To determine the coefficients of the numerator and denominator, one has to solve systems of linear algebraic equations. This is an ill-posed operation, so the PA coefficients will be determined with large errors. However, these errors are in a sense self-compensated; therefore, the PA's approximation of the original function is high. This is the radical difference between PA and the Taylor series, in the calculation of which the error only grows with an increase in the number of terms.
The occurrence of the autocorrection property has been numerically tested for a number of special functions. At the same time, for elliptic functions, the so-called spurious pole-zero pairs or "Froissart doublets" appear, formed by closely spaced (but different and irreducible) zeros and PA poles. This suggests that the PA in the general case is illposed. The reason the PA is ill-posed is that it is related to an analytic continuation, since the aim is typically to gain information about a function in a region of the complex plane based on information at a single point [8].

Generalizations of the PA
Among the methods for removing doublets, one can be mentioned that is based on explicit computation of zeros and poles [20,21].
For eliminating spurious poles, in [22] the smoothed PA was proposed. Its essence is to use the averaged (in a certain sense) expression instead of the usual diagonal PA The authors of [23] propose another variant of a weighted average in which the spurious zero-pole pairs give a strongly reduced contribution.
The so-called robust PA [8] is based on slightly reducing the accuracy of the solution to a certain matrix problem; the benefit is that the resulting PAs will be pole-free in most regions where the original function itself is pole-free.
The PA described above is sometimes called the PA of the first kind [5,24]. The PA of the second kind, or multipoint PA (MPA), is a rational function [ / ] ( ) n m f ε of the form (2) such that its values for the 1 + + n m values of the argument coincide with the values of ( ) f ε at the same points. This gives a system of linear algebraic equations for determining the coefficients of the polynomials of the numerator and denominator [24]. The MPA in many cases provides effective interpolation formulas and allows for extrapolating the values of a function specified on a limited interval beyond its limits [25,26].
The PA allows you to successfully work with functions that have poles. However, it is often necessary to investigate functions that have branch points and build all their branches. In this case, Hermite-Padé [27] approximations can be used. Let it be a question of a function having the following expansion If function (11) has a branch point, then to search for branches, one can go to an implicit function and construct a polynomial whose coefficients are determined from the condition unknowns and condition (12) yields N linear algebraic equations; therefore, there should be After the polynomial p F has been found, one can easily find the p branches of the solution to the equation 0 = p F (16) There also exist many other generalizations of the PA. E.g., the Padé-Fourier and Padé-Chebyshev approximants [7,11] use trigonometric functions or orthogonal polynomials. A PA in two variables is called a Chisholm approximant; a PA in multiple variables is called a Canterbury approximant [7].

Using PAs in Nonlinear Mechanics
Initially, PAs were applied in physics for the analysis of divergent and asymptotic series [5,6] and in the theory of critical phenomena. They are also widely used in the mechanics of solids, especially in such areas as wave dynamics (improving convergence, overcoming the Wilbraham-Gibbs phenomenon, and building concentrated nonlinear waves [3]), discrete lattices theory (construction of improved continuous models), and the theory of composites [28].
It is worth noting the geometric series method for the construction of nonlinear concentrated waves [29][30][31]. In this approach, the solution is found in the form of a series in powers of the exponential function and, using the PA, the conditions under which the series becomes geometric are revealed.
Using PAs in asymptotic approaches with unconventional perturbation parameters [3,32] can greatly extend the range of applicability of perturbation series. Especially noteworthy is the method of the homotopy ("artificial small") parameter, which made it possible to propose a new technique for solving mixed problems of mathematical physics [3]. Two-and multipoint PAs (the TPPA, the MPA) proved to be very effective tools for constructing global approximations from asymptotic expansions [2,3,25,26]. A natural generalization of this approach is the two-point quasifractional approximant method [26], which is also called the asymptotically equivalent functions method [2,3] or the asymptotic approximants approach [33].

Using PAs in Boundary Layer Problems
The pioneering papers by Van Dyke [34][35][36] are devoted to the application of PA in hydrodynamics. An overview of the work done in this area up to 1994 is contained in the monograph by Pozzi [37]. Since that time, PAs have taken their rightful place in the mechanics of liquids and gases [38,39]. The results obtained after 1994 are partially highlighted in our review [39]. In this section, we focus on applying PAs to boundary layer problems.
The history of the emergence and development of the concept of a boundary layer is described in detail in [40][41][42]. Recall that the key (for singular asymptotics) concept of an edge effect occurred in the works of H. Lamb and A.B. Basset in 1890, while the concept of a boundary layer occurred in fluid mechanics in 1904 [43]. An extremely interesting and useful visual illustration of this concept is the album of photographs of liquid and gas flows [44]. The works of Vishik and Lyusternik [45,46] played a significant role in the mathematical formulation of the boundary layer idea. They considered boundary layers for ODEs and linear and quasilinear PDEs of elliptic and parabolic type. The initial layers [47,48] were also studied.
In various fields of science, the terms 'stiff equation', 'contrast structure', 'edge effect', and 'skin effect' are used to describe the phenomena of the boundary layer. In any case, such an asymptotic phenomenon is based on the principle of spatial or temporal localization [49].
Before the computer revolution, methods for solving boundary layer equations relied on phenomenological or asymptotic approaches. Kármán and Polgauzen [50,51] proposed that we limit ourselves to satisfying these equations only on average over the thickness of the boundary layer. For this purpose, one should use the momentum theorem and replace the differential equations with integral relations obtained from the equation of motion by integrating it over the thickness of the boundary layer.
The Kochin-Loitsyanskii [50,51] method uses the idea of an "intermediate" asymptotic. An arbitrary distribution of velocities on the outer boundary of the boundary layer is replaced by a linear distribution in the first approximation, a quadratic distribution in the second, etc., representations that perform smooth matching with the solution in the main zone. Later, this idea was developed by Loitsyanskii in the method of "generalized similarity" [51].
The computer revolution has led to the creation of effective numerical methods for solving the boundary layer equations [52,53]. This devalued a number of phenomenological approaches, but only increased the role of asymptotic methods [54][55][56].
A feature of boundary layer problems is the nonuniformity of the asymptotic expansion. The boundary layer is the region of singularity, and its solution is called the inner asymptotic. The solution in the outer region is called the regular (outer) asymptotic. If the conditions of the Kaplun theorem [57] are satisfied, the regular asymptotic can be uniformly extended in such a way that an overlap region arises where it is possible to compare the inner and outer asymptotics. The method of matched asymptotic expansions [54,58] is based on this. However, it is not always easy to justify the applicability of the Kaplun theorem [57]. This implies the relevance of methods for the synthesis of various limiting asymptotics. Below, we consider the matching of inner and outer asymptotics using the PA. Let us consider as an example the laminar boundary layer near a flat plate [50] described by the Blasius equation with BCs where ( ) = x ψ φ ζ The inner asymptotic (17) The construction of the outer asymptotic is nontrivial due to the presence of logarithmic terms. In [38,39], the mechanism for obtaining and estimating both the main and minor terms of the outer asymptotic is described in detail. The expression for the outer as- The constants 2 , , a c D are determined by the integral relations The procedure for obtaining the PA consists in fulfilling local equalities for relations (20) and (21)

Application of the PA to the Problem of Rotating Fluid of the Ekman Layer Type
We consider a flow established along the z-axis direction ( Figure 2) in a rotating rectangular channel with a side ratio of the order of unity. The distribution of the longitudinal component of the velocity is inhomogeneous along the transverse direction x; therefore, the influence of the convective terms cannot be neglected. However, at a high speed of rotation, the flow is still divided into a core and thin shear layers such as the Ekman layer on surfaces perpendicular to the axis of rotation [59]. In this case, in contrast to the slotted channel, the layers develop under conditions of an inhomogeneous external flow with a velocity shift in the transverse direction.
Problems of this type are considered in papers [60][61][62][63][64][65][66][67]. Papers [63,64] use the standard reception of decomposition of the sought functions into Fourier series with the further determination of parameters of decompositions by the method of successive approximations. This technique allows us to determine with sufficient accuracy the average characteristics of the flow in the channels. Numerical methods for the solution were used in works [61,62] and [65][66][67]. It should be noted that in order to correctly apply numerical methods, it is necessary to take into account the region of singularity of flow near the walls of the channel. In practice, this means that the finite difference grid (or cells in the finite element method) in these regions must be sufficiently shallow (thickened). The thickening methods can be determined analytically or a priori from the flow structure. Methods for such numerical solutions are considered, for example, in works [68,69]. The most stringent of these is Direct Numerical Simulation (DNS). This refers to the numerical integration of three-dimensional non-stationary Navier-Stokes equations on grids, which provide the resolution of all vortex flow structures with minimum scale dimensions. Thus, the DNS method is based on the basic principles of hydro-aerodynamics and is free from empiricism. However, this approach requires a lot of machine time and, given the nonlinear nature of the equations, it is extremely difficult to use parallel algorithms. The second-strictest approach is the so-called Large-Eddy Simulation (LES) method, first proposed in Smagorinsky's work [70]. Within the framework of this method, not the original, but prefiltered Navier-Stokes equations are solved. As a result, using LES, only the main "energycarrying" or "large" vortices are accurately resolved, the sizes of which depend on the size of the filter (its role is usually played by the computational grid used), and smaller vortex structures are approximately modeled using semi-empirical models. Finally, the third relatively new class of approaches, Hybrid RANS-LES Methods (HRLMs), includes a large group of methods that use one or another combination of RANS and LES methods.
The application of these methods in the calculation of internal currents in rotating channels requires experience with special and expensive application packages. The main difficulty of using numerical methods taking into account the problem of isolating the thin boundary layer area is the lack of a priori information about the structure of the calculation area. In this case, the method we propose can be used to evaluate the geometric parameters of the region. To characterize the thickness of the Ekman boundary layer, you Let us introduce the stream function In accordance with Praudman's theorem [59], the influence of rotation leads to the equalization of the flow characteristics in the direction of the axis of rotation. Strong viscosity-based rotational effects contribute to the formation of thin (relative to the flow region in the channel) Ekman boundary layers. Following the basic concept of boundary layer theory [50,51], we accept the following assumptions regarding the diffusion of velocity components within the Ekman layer: This assumption is based on the analysis of the dimensionality of the specified values for the boundary layer on a flat surface.

Rejecting the terms
x y x y y ψ ψ ψ ω (30) On the wall we suppose adhesion conditions We choose the channel half-height h and the local velocity in the median plane ( ) W x as the length and velocity scales in the section = x const . We also introduce the h dimensionless quantities ( ) ( ) ( ) Here, * P is the modified pressure, ω is the angular velocity of rotation, v is the coefficient of kinematic viscosity, and ρ is the density.
We introduce two dimensionless complexes Equations (32) The choice of parameter f in our example was made on the basis of speed measurement data in the rotating channel [60].
In this domain, the dependence of the flow characteristics on the local Reynolds number c Re is weak. This is a consequence of the fact that, due to the relative thinness of the Ekman layers, the scale of the velocity variation in the transverse direction significantly exceeds the linear scale of the velocity variation along the direction normal to the surface. Therefore, in Equations (32) and (33), the right-hand sides can be neglected. Then, in the local approximation, the inhomogeneity of the external flow is provided by the value characterizing the shear in the transverse direction [61]. It should be noted that a similar situation is observed when modeling the nonlinear Ekman layer in geophysical currents caused by the rotation of the Earth [59]. Then, Equations (32) and (33) can be written as follows ( ) The BCs are 0: In what follows, we will consider the channel half-height, counting the coordinate from the lower surface (bottom) of the channel. The presence of side walls of the channel assumes that the liquid flow rate through any section of the channel parallel to = x const is 2 0 0 =  u dη (38) Condition (36) Integrating Equation (34), we obtain To determine the asymptotics of Equations (34) and (35) We transform Equations (38) and (39), omitting the dashes above the dimensionless velocity components Further, we suppose − → ∞ Re ω and consider the flow formed under the strong influence of rotation on the flow. Then, the system of equations of motion has the form ( ) with BCs 0: : The last two conditions in (45) We substitute the expansions (49) and (50) into Equations (42) and (42), taking into account that from the boundary conditions (44) it follows that 0 Equating the coefficients at the same degrees ζ , we obtain a system for determining the parameters of the expansions (49)-(50) From this system, one obtains To determine the parameters 1 2 , a b , we use the integral relations [39]. Integrating Equations (42) and (43) with weight 1, we obtain The inner asymptotics for the velocity components depend on f and have the form ( ) Note that at 0.95 = f , the pressure side (the vertical channel wall) experiences the greatest influence of the secondary flow u. In addition, the choice of 0.95 = f is due to the experimental data [60,62].
For 0 = f , one obtains the well-known Ekman solution [59], which can be found analytically from (42) and (43)  The interest in solving the Ekman boundary layer problem on surfaces perpendicular to the axis of rotation is due to practical applications in more complex problems of calculating internal currents in rotating systems. These may be the regions of fluid flow between the blades of turbogenerators. The results of solving the problem can be used to estimate the intensity of secondary currents in the flow core, as well as to estimate the effects of rotation on the coefficient of resistance. The intensity of the secondary currents is determined by integrating the velocity component  [59,62] and are in qualitative agreement with the experimental data [60].