# Full text of "Padé-type rational and barycentric interpolation"

## See other formats

o O (N < (N X Pade-type rational and barycentric interpolation Claude Brezinski* Michela Redivo-Zaglia^ November 1, 2011 O ■ Abstract (N In this paper, we consider the particular case of the general rational Hermite interpolation problem where only the value of the function is interpolated at some points, and where the function and its first derivatives agree at the origin. Thus, 0^ I the interpolants constructed in this way possess a Pade-type property at 0. Nu- merical examples show the interest of the procedure. The interpolation procedure can be easily modified to introduce a partial knowledge on the poles and the zeros of the function to approximated. A strategy for removing the spurious poles is explained. Applications are given. ,S^ ■ Keywords: Rational interpolation, Pade-type interpolation, rational representa- C^ , tion, barycentric representation, piecewise rational interpolation. Key words and phrases. Rational interpolation, Pade-type approximation, barycentric formula. > ; 2010 Mathematics Subject Classification. Primary 65D05, 65D15, 41A20, 41A21. vn 00 ^^ ; 1 Pade-type approximation and rational interpolation O I For representing a function /, rational functions are usually more powerful than poly- nomials. The information on the function / can consist either of the first coefficients of its Taylor series expansion around zero, or in its values at some points of the complex plane. In the first case, Pade-type, Pade, or partial Pade approximants can be used. These are rational functions whose series expansion around zero (obtained by EucUdean di- vision in ascending powers of the numerator by the denominator) coincides with the series / as far as possible. In Pade-type approximation, the denominator can be arbi- trarily chosen and, then, the coefficients of the numerator are obtained by imposing the preceding approximation-through-order conditions. In Pade approximation, both the *Laboratoire Paul Painleve, UMR CNRS 8524. UFR de Mathematiques Pures et Appliquees, Universite des Sciences et Technologies de Lille, 59655-Villeneuve d'Ascq cedex, France, E-mail: Claude .BrezinskiSuniv-lillel .f r. TUniversita degli Studi di Padova, Dipartiniento di Matematica Pura ed Applicata, Via Trieste 63, 35121-Padova, Italy. E-mail: Michela. RedivoZagliaSunipd. it. denominator and the numerator are fully determined by these conditions. For partial Fade approximants, a part of the denominator and/or a part of the numerator can be arbitrarily chosen, and their remaining parts are given by the approximation-through- order conditions. On these topics, see |l][7|[T ^ fT5 ] . In the second case, an interpolating rational function can be built using Thiele's formula, which comes out from continued fractions (see, for example, [T3| pp. 102ff.] or |T6| Sec. III.3-4]). It achieves the maximum number of interpolation conditions, and, so, no choice is left for its construction jT^. The same is true for Hermite rational interpolants, a subject treated in many publications (see, for example, [25]) which is related to Newton-Fade approximants [TB| p. 157]. On the other hand, when the degrees of the denominator and of the numerator are the same, writing the rational interpolant in a barycentric form allows to freely choose the weights appearing in this formula. These weights can be chosen by imposing various additional conditions such as monotonicity or the absence of poles pi-H|[TB]. For an interesting discussion between the coefficients of the interpolating rational function and the weights of its barycentric representation, see |5]. For the important problems of the ill-conditioning of rational interpolation, and of the numerical stability of the algorithms for its solution, consult |5|fT9]. In this paper, we will construct rational functions possessing both properties, that is interpolating / at some points of the complex plane, and whose series expansion around zero coincides with the Taylor series / as far as possible. Of course, this case is a particular instance of the general rational Hermite interpolation problem treated in its full generality in [22], for example. Then, using a different number of conditions than required, we are able to construct rational interpolants in the least squares sense. We will also show how informations on the poles and the zeros of / could be included into these interpolants in a style similar to the definition of the partial Fade approximants [8] . 2 Problem statement We consider two different arguments. • Let / be a function whose Taylor series expansion around zero is known. A Fade- type approximant of / is a rational function with an arbitrarily chosen denominator of degree /c, and whose numerator, also of degree /c, is determined such that the power series expansion of the approximant around zero coincides with the devel- opment of / as far as possible, that is up to the term of degree k inclusively [H]- By choosing the denominator appropriately, this rational function has a series ex- pansion which agrees with that of / up to the term of degree 2k inclusively. It is then called a Fade approximant, and there is no freedom in the choice of the coefficients of the numerator and the denominator of the rational approximant. On this topic, see, for example [HE]. • Let / be a function whose values at A; + 1 distinct points in the complex plane are known. It is possible to construct a rational function, with a numerator and a denominator both of degree k, which interpolates / at these points. If this rational function is written in barycentric form, it depends on k nonzero weights which can be arbitrarily chosen. But, by Thiele's interpolation formula, it is also possible to obtain a rational function, with a numerator and a denominator both of degree k which interpolates / at 2A; + 1 distinct points in the complex plane. In that case, there is no freedom in the construction of the rational interpolant. We now consider these two themes together and work in both directions in a different way. Each of these choices leads to a different rational function whose series expansion agrees with that of / as far as possible, and which interpolates / at distinct points in the complex plane. • We determine the denominator of the Pade-type approximant so that it also in- terpolates / at as many distinct points in the complex plane as possible, which is k points. Thus we obtain a rational function interpolating f ai k points and with order of approximation fc + 1 at 0. Such a rational function will be called a Pade-type rational interpolant. • We determine the weights of the barycentric formula for the rational interpolant so that its power series expansion coincides with that of / as far as possible, that is up to the term of degree k — 1 inclusively. This approach produces a rational function with order of approximation k at 0, and interpolating f at k + 1 points. Such a rational function will be called a Pade-type barycentric interpolant. In each case, different interpolation or approximation conditions can be considered, and the rational function can be computed in the least squares sense. Rational inter- polants with arbitrary degrees in the numerator and in the denominator of the inter- polant could also be defined similarly. Let us mention that it is also possible to work with the reciprocal function g of /, and its reciprocal series which is defined by the algebraic relation f{t)g{t) = 1. In the sequel, the formal power series / will be written as f{t) = Co + cit + C2t^ + ■■■ 3 Pade— type rational interpolants We will begin by treating the case of a formal power series and, then, we will consider a series in Chebyshev polynomials. 3.1 Power series Let Rk be written as Rk{t) Nk{t) _ao + ait-\ h Okt'' Dk{t) ~ bo + bit + --- + bkt''' If the coefficients bi of the denominator are arbitrarily chosen (with bk 7^ 0), and if the coefficients a, of the numerator are computed by the relations ak cobo cibo + cobi Ckbo + Ck-ibi + cobk (1) then Rk is the Pade-type approximant {k/k)f of / which satisfies the approximation- through-order conditions f(t) — Rk{t) = 0{t''~^^). Let us remind that this condition means that /, R^., and their derivatives up to the kth inclusively take the same values at the point t = 0. Replacing Qq, . . . ,ak by their expressions ([T]) in A'^^, and gathering the terms corresponding to each bi, we obtain Nk{t) = boSkit) + bitSk-iit) + ■■■ + bkt'Soit), with n 0,1, (2) (3) Suit) = co + cit + --- + cX, Let us now determine bo, . . . ,bk such that R{Ti) = f{Ti){=: fi) ioi i = 1, . . . ,1, that is such that iVfc(r.) - /.A(r,) = 0, t = l,...,l, where Ti, . . . ,Ti are distinct points in the complex plane (none of them being 0) . We obtain the system {Sk{r^) - Idbo + r,(5fe„i(r,) - j-)b^ + ■■■ + rf (5'o(r.) - /.)6, = 0, ^ = 1, ...,/. (4) Since a rational function is defined up to a multiplying factor, we set bo = 1 (imposing another normalization condition could lead to feo = and, so, Oq = 0, thus reducing the degree), and we obtain a system of / linear equations in the k unknowns bi, . . . ,bk. li I ^ k, we consider its least squares solution (/ > k, overdetermined system) or its minimum norm solution (/ < k, underdetermined system). Therefore, the 6j's are first determined by the interpolation conditions and, then, the a^'s are calculated by formulae O- Multiplying each equation in (jlj) by t[^'' (the reason will be made clear later) and using ([2]), we obtain the following Property. Property 1 When I = k, it holds Rk{t) Skit) r^'iSkin) - -h) tSk-i{t) r^'-^'iSk^^in) - -h) ■ ■ t'Soit) ■ So{n)-f, rk\Sk{Tk) - -fk) r^'^\Sk-i{rk) - -fk) ■ ■ SoiTk)-fk 1 n'iSkin) - -h) t n'-^'iSk^iin) - -/i) ■■ t^ ■ So{T,)-h Tk'iSkiTk) - -fk) Tk'^'iSk^liTk) - -fk) ■■ ■ SoiTk) - fk Let us take t = Tj in this formula, and multiply the first row of the numerator and of the denominator by t~^. Then, subtract the row z + 1 of the numerator from the first one. This first row becomes t^^ fi^ T~''~^^fi, . . . , fi, and we obtain Rk{Ti) = fi, for i k, since the first row of the denominator is r,,- -fc+i 1. Thus the interpolation property of Rk has been recovered from its determinantal expression. Let us now define the linear functionals Li acting on the vector space of polynomials by (this is the reason for multipying each equation in (jlj) by r^^ ) Ht^ \S,{n)-f,), J = 0,1, 1,2, The polynomial Pk{t) = Dk ifc-i rriSk{ri)-fi) -k+l (S.^iin) - h) r^'{S,{n)-U) r,'^'iSk-i{r,)-h) Soin) - /i Soin) - fk where Dk is any nonzero normalization factor, satisfies the so-called hiorthogonality conditions L,{Pkit))=0, z = l,...,/c. Such a polynomial is the kth member of the family of formal biorthogonal polynomials with respect to the linear functionals {Li} [HI pp. 104ff.], and we see that the denomi- nator of Rk is equal to Pkit) = t''Pk(t~^). Let now c be the linear functional acting on the vector space of polynomials and defined by c(x*) = Cj for i = 0,1, . . ., let Qk be the polynomial of degree A; — 1 in t Qkit) C X Pk{x)-Pk{t) X and set Qkif) = t^ ^Qk(t ^). From the definitions of Pk, Qk, and the determinantal formula of Rk given in Property [H we have the following Property. Property 2 Rkit) = Co + t Qkjt) Pk{t) ' when I = k. This Property shows that Rk is exactly the generalization of the Pade-type approx- imants defined in j2l pp. 97ff.], and, thus, it holds Rk{t) — f{t) = 0{t'^~^^) as required by our approximation-through-order conditions. It is possible to construct Pade-type rational interpolants {p/q)f with an arbitrary degree p in the numerator and q in the denominator, and then to determine its denom- inator in order to satisfy q (or even / ^ q) interpolation conditions [TitlT5]. Let us set Np{t) = flo + Ciit + ■ ■ ■ + apf^, and Dq{t) = bo + bit + ■ ■ ■ + bqf^. The coefficients of the denominator are first computed as the solution of tlie system (j4]) with / = g (or even I ^ q). Then, the coefficients of the numerator are given by Co = = Cofco ai = = cibo + cobi Op = = Cpbo + Cp„ 161 r Cp^qOq, with the convention that Q = for i < 0, and the partial sums ([3]) computed accordingly. Such an interpolant satisfies {p/q)f{Ti) = fi ior i = 1, . . . ,q and {p/q)f{t) — f{t) = If some poles and some zeros of / are known, this information could be included into the construction of the rational interpolant. Let pi, . . . ,Pm and Zi, . . . ,Zn be these poles and zeros, respectively. Setting Pm(t) = (t — Pi) ■ ■ ■ {t — Pm) and Z„(t) = (t — Zi) ■ ■ ■ (t — Zn), we are looking for the rational function _ N,{t)Z^{t) """^'^ - D,{t)PUt) again such that Rkiji) = f{Ti)(= fi) for i = 1, . . . ,k, and such that f(t) — Rk(t) = 0{t''~^^). Such a rational function is called a partial Pade-type rational interpolant since it is similar to the partial Pade approximants introduced in jB], but with a lower order of approximation. We must have Nk{n)Zn{Ti) - fiDk{Ti)Pm{ri) = Nk{T.i) - f^^^Duir,) = 0, z = 1, . . . , fc. Setting Nk and Dk as above, the coefficients of Dk are first determined as the preceding ones with fi replaced by fiPm{Ti)/Zn{Ti) in the system (j4]), and then the coefficients of Nk are obtained by the same relations as before where, now, the coefficients q have to be replaced by those of the series expansion of f{t)Pm{t)/Zn{t) in (^. Thus, we first compute the coefficients of h(t) = f{t)/Zn{t) by identification in the relation f[t) = h{t)Zn,{t). Then the coefficients of f{t)Pm{t)/Zn{t) = h{t)Pm{t) are obtained by a simple product. These coefficients replace the q's in the definition of the partial sums ([3]) . Let us mention that this division and the following multiplication can be performed monomial by monomial in order to avoid the computation of the coefficients of the polynomials Zn and Pm- Indeed, we can begin by computing the coefficients of f{t)/{t — Zi), then, from these coefficients, we compute those of {f (t) / {t — Zi)) / {t — Z2) , and so on until the division by (t — Zn)- Thus, we obtain the coefficients of h. Then, we formally multiply h{t) by (t — pi), the result by (t — P2), and so on until (t — Pm) which gives the coefficients of h{t)P^{t) = /(t)P^(t)/Z„(t). 3.2 Fourier and Chebyshev series Fourier series can be approximated similarly by a procedure introduced in [30] and de- veloped in [TT]. It consists in adding to the Fourier series its conjugate series, thus transforming it, by a change of variable, into a power series, then computing the inter- polants as described above, and finally keeping only their real part. The approximation of parametric curves is another topic which could be explored in the continuation of | Let us consider the case of a series in Chebyshev polynomials oo i=l where Tj(t) = cos(zarccost). The rational interpolant Rk is defined as ho/2 + hiTi{t) + --- + hkTk{t) Rkit) eo/2 + eiTi{t) + --- + ekTk{t) Adapting to our case a general approach due to Hornecker |22 [ [23 ] and particularized by Paszkowski Pl] using the multiplication law Ti(t)Tj{t) = {T\i_j\{t) + Tj+j(t))/2 of the Chebyshev polynomials, we have Rk{t) — f{t) = 0{Tk^i{t)) for any choice of the coefficients Cj of the denominator, if the coefficients hi of the numerator are computed by k ho = coeo/2 + ^ CiCi i=l k hn = {cneo + ^{c\n-j\ + Cn+j)ej)/2, n = l,...,k. i=i Let us now choose cq, . . . , e^ such that Rk{Ti) = fi for i = 1, . . . ,k. Similarly to the procedure followed for a power series, these coefficients must satisfy k k / k \ /'^X coeo/2+^Cjej+^ I c„,eo + ^{c\n-j\ + Cn+j)ej J Tn{Ti)-fi I cq + 2^ejTj{Ti) | = 0, j=l n=l V J=l / V j=l / for i = 1, . . . ,k, thus leading to the system k \ k / k Co/2 + ^ CnTnin) - /i I eo + ^ I Cj + ^(C|„_j| + Cn+j)Tn{Ti) - 2fiTj{T, n=l / i=l V n=l for 2 = 1, . . . , fc. Since a rational function is defined apart a multiplying factor, we set Co = 1 for solving it. This approach can be extended to a numerator of degree n + /c, /c > 1 [TUl pp. 161ff.], [3 pp. 220ff.]. Moreover, since a Chebyshev series is a cosine series, its conjugate series could be added to it, as indicated above for Fourier series, then a rational Fade- type interpolant could be constructed, keeping only its real part. 7 4 Pade— type barycentric interpolants We consider the following rational function, in barycentric form, k Wi ^-^ t — Ti Rtit)-- "° tj^ i=0 t -Ti where fi = f{Ti). This rational function interpolates / at the k + 1 points Ti,i = 0, . . . ,k, whatever the Wi ^ are. It is well-known that, by the Lagrangian interpolation formula for the denominator of Rk, Wi = qi/v'iji) with v{t) = n7=o('^ ~ '^i)' ^"^^ '^''i'^i) — Y[i=o j^ii'^i ~ '''3)1 where gj is the value of the denominator of R^ at the point Tj. This remark shows that, as in the case of Pade-type rational interpolation, the rational interpolant Rk is fully determined by its denominator as mentioned in [5]. Let us remind that, for the choice Wi = v'{Ti), Rk becomes a polynomial and that, for several choices of the points Ti closed expressions of the weights Wi are known. Let us now determine wq, . . . ,Wk such that fit)~Rkit) = Oit'). In that case, Rk is a Pade-type approximant {k/k)f of /, but with a lower order k of approximation instead of k+1. This condition means that / and Rk and their derivatives up to the {k — l)th inclusively take the same values at the point t = 0. Let us mention that it is not possible to improve the order of approximation for obtaining an exact Pade-type approximant. The preceding approximation-through-order condition reads — t — Ti ^ — ' t — Ti Dividing each fractional term by the corresponding Tj (obviously all the Tj have to be different from zero, which is not a restriction since our Pade-type barycentric interpolant will interpolate / at t = 0), changing the signs, and using the formal identity 1 t t^ l + - + ^ + --- , t/Ti Ti T- 2 we have fc / , ,9 s k ,_n '« \ '^ 'i / i_n '* \ '« 'j / i=0 - ^ 'I / j=o Identifying on both sides the coefficients of the identical powers of t leads to E?/.-«E i=0 Ti 1=0 Wi 1 v-^ W - + Cl > , r. EWi ± Y^ Wi 1 Y^ Wi —fi- = Co > + Cl X — Ti Ti ^-^ Ti Ti ^ — ^ '^• i=0 « « j=0 * * j=0 fc fc fc A: EWi \ sr^Wil sr^Wil sr^VJi j=0 i=0 i=0 i=0 and so on up to the term of degree k — \ inclusively. Thus, the Wi must be the solution of the linear system ^{fi - Co' Wi Ti j=0 k >; j=0 a _ Co _ \ \Ti Ti ^) Wi -0 Wi Ti k i=0 ^ /. Co Cl = 0. (5) Since a rational fraction is defined apart a multiplying factor in its numerator and in its denominator, we will set wq = 1 and, thus, we obtain a system of k equations in the k unknowns Wi, . . . ,Wk- This approach needs the knowledge of the values of / at A; + 1 points, and that of the coefficients Cq, . . . , Ck-i- Let us write the system ([5]) as y^^ajiWj = 0, j = l,...,k. i=0 Then, we obtain two determinantal expressions for Rk, the ffist one in a barycentric form, and the second one in a Lagrangian-type basis. Property 3 Rk{t) /o/(t-ro) ■ ■ fklit-Tk) aio aik a/co akk l/(t-ro) ■ ■ l/it-Tk) aio aik OfcO akk foLoit) O-IO OfcO Lo{t) aio fkLk{t) a-ik dkk Lk{t) O-lk dkk with, for i = 0, . . . ,k, an = ifi-co) — _ Qj-i,i - Cj.iWj Ujji , J Z,,...,fv, Ti and H-t) = \{it-rj j=0 The second formula comes out from Li{t) = L{t)/(t — Ti) with L(t) = (t — tq) ■ ■ ■ (t — Tfc). Since Li{Tm) = for m ^ i and Li{Ti) 7^ 0, we immediately recover, from the second expression, the interpolation property Rk{Ti) = fi for i = 0, . . . , /c. For recovering the approximation-through-order property, the expressions l/(t — Tj) in the numerator and in the denominator of Rk have to be replaced by — l/(rj(l —t/ri)) = — (1 + t/rj + t^/rf + ■ ■ ■)/Ti, and the coefficient of each power of t has to be separately identified up to the degree k — 1 inclusively. Assume now that only cq, . . . , Q_i are known, with / < k. We can choose wq, . . . ,Wk such that f{t) — Rk{t) = (9(t') by considering only the first / equations of the preceding system, and replacing the last ones by the equations = 0, j = l,...,/c-/, which is equivalent to considering that the coefficients Q, . . . , Ck-i are zero in the system (^ . The rational function R^ interpolates f in k + 1 points and its expansion coincides with that of / up to the term of degree / — 1 inclusively. It is also possible to consider the case where / > k coefficients of the series of / are known. Adding to the preceding system the equations k / sr^ I fi Co Ci \ Wi Cj] — = 0, 3 = k,--- ,1-1, ^ \t^ t^ t^'^ V Ti i=0 \'t 'i 'i / * and solving it in the least squares sense leads to an approximation i?^ whose series expansion agrees with that of / only in the least squares sense, and which interpolates / at fc + 1 points. Let us again consider the case where some poles and some zeros of / are known. The rational function h Wi XI J^fiPm{Ti)/Zn{Ti Pmit) ^ E Wi 1=0 10 interpolates / at the k + 1 points Ti, i = 0, . . . ,k, whatever the Wi ^ are, and it has the poles Pi, . . . ,Pm and the zeros Zi, . . . , Zn- Thus is can be constructed as above after replacing everywhere fi by fiPm{ri)/Zn{Ti), and we obtain f(t) - Rk{t)Zn(t)/Pm(t) = When the poles of / are known, an explicit expression for the weights of the near- best rational interpolants in a Chebyshev sense can be obtained [29]. As mentioned in this paper, the knowledge of the poles dramatically improves the interpolation process as can be seen from the numerical examples given there, and also in [5|I2B]. 5 Numerical examples We will now show some numerical examples which gather several interesting properties that will allow us to exemplify the effectiveness of our procedures. But, before, let us give the following consistency property Property 4 If f is a rational function with a numerator and a denominator both of degree smaller or equal to k, then, our two procedures produce a rational function R^ which is identical to f when I = k. This property comes out from the fact that R^ is defined by a set of linear equations which is the same as the set of equations which defines /, and the result follows from the uniqueness of R^ . Let us remind that the solution of a rectangular system of equations Ax = 6 of maximal rank r = min(/, k) with A G C'^'^ is x = A%, where A^ is the Moore-Penrose pseudo-inverse of A defined by A^ = {A*A)~^A* ii r = k < I (overdetermined system) and A^ = A*{AA*)"^ ii r = I < k (undertermined system). If the rank r is not maximal, then A^ = V'E'^U* where A = UT,V* is the singular value decomposition of A. The Matlab instruction pinv(A)*b gives the least squares solution when the system is overdetermined (that is the unique solution minimizing the 2-norm of the residual if the matrix is full rank, and the vector of minimal 2-norm among those minimizing the 2-norm of the residual, if not), and the minimal 2-norm solution when the system is underdetermined. In all cases, the computation is based on the singular value decomposition of A. Example 1: a function with poles We consider the following function, and its series expansion ^ tanM ^ , ^ 1^.,. ^ 1^4,4 ^ iL^6,6 + _6L^8,8 + . . . ■'^ ^ ut 3 15 315 2835 This function has poles at odd multiples of 7r/(2a;), and zeros at odd multiples of ir/u, except at 0. We considered two sets of interpolations points: equidistant points in the interval [—1, +1], and the {k+1) roots of unity on the unit circle, defined by Tj = exp[2ij/{k+l)) 11 tan(4*x)./(4*x): Pade-type rat. interp.: complex pts (dashed), real pts (solid) Figure 1: Pade-type rational interpolant with A; = 8 for tan(4t)/(4t) for j = 0, . . . ,k, where i = \/—l. This complex choice was discussed in J2T]. Let us mention that none of the interpolation points Tj should be 0, since, because it is the point where the Pade-type approximants are computed, this point will appear twice in the list of interpolation points, which leads to a redundancy and a loss of information. Pade-type rational interpolants The errors obtained with the Pade-type rational interpolants are given in Figure [T] for w = 4 and k = 8. The solid line corresponds to the case of the real interpolation points, while the dashed one refers to the points on the unit circle (in that case, we only keep the real part of Rk{t)). The poles of / are, in the interval considered in Figure [H at the points ±0.39269908 ... and ±1.1780972 . . ., and the zero at ±0.78539816 . . . Since the poles and the zeros are known, we also took Z2{t) = {t — 7r/4)(t + 7r/4) and P2{t) = (t — 7r/8)(t + 7r/8). The errors obtained with the partial Pade-type rational interpolants are displayed in Figure [2j The improvement brought by taking partially into account the knowledge on the poles and on the zeros is clear. Choosing the zeros of the Chebyshev polynomials as the real interpolation points in [—1, ±1] does not change much the quality of the results for such small values of k. Pade-type barycentric interpolants Let us now consider the same example but with the Pade-type barycentric inter- polants. The results are given in Figure [3l With the partial Pade-type barycentric interpolants, we obtain the results of Figure m In both figures, the solid line corresponds to the real interpolation points, and the dashed one to the interpolation points on the unit circle. Let us mention that, with the Shepard's weights Wi = l/(t — Tj) [27], the interpolants 12 tan(4*x)./(4*x): Partial Pade-type rat. interp.: complex pts (dashed), real pts (solid) 10"^ Figure 2: Partial Pade-type rational interpolant with A; = 8 for tan(4t)/(4t) have poles around —0.4 and +0.4. Example 2: a function with a cut We consider the series fit) log(l + 1) t 1- t t' which converges in the unit disk and on the unit circle except at the point —1 since there is a cut from —1 to — oo. Pade-type rational interpolants For a Pade-type interpolant of degree 7, we consider equidistant real interpolation points in the interval [—0.9, +1.2]. For 7 points (solid line), the system to be solved is square. For 3 points (dashed line) and 14 points (dash-dotted line), the system is solved in the least squares sense as explained above. See Figure O We see that quite good results are obtained for values of t far outside the convergence interval. Pade-type barycentric interpolants The interpolation points Tj are equidistant in [—0.9, +4], and k = 7. In FigureEl three types of weights Wi are considered: those corresponding to the Pade-type barycentric interpolants are the same as explained above (solid line), the weights Wi = (—1)* of Berrut |i2j (dashed line), and the weights Wi = l/{t — Ti) suggested by Shepard |27] (dash-dotted line), these last two choices ensuring pole-free interpolants on the real line. 13 Pade-type baryc. interp.: real pts. (solid), complex pts. (dashed) -1.5 -1 -0.5 0.5 1 1.5 Figure 3: Pade-type barycentic interpolant with k = 8 for tan(4t)/(4t) Partial Pade-type baryc. interp.: real pts. (solid), complex pts. (dashed) -1.5 -1 -0.5 0.5 1 1.5 Figure 4: Partial Pade-type barycentric interpolant with k = 8 for tan(4t)/(4t) Example 3: a continuous function We consider the exponential function t e /(t) = e* = i + Yy + ^ + 14 Degree k=7, 1=3 (dashed), 1=7 (solid), 1=14 (dash-dot) Figure 5: Pade-type rational interpolants with A; = 7 for log(l + t)/t, and 3, 7 and 14 int. pts Pade-type baryc. interp. (plain), Berrut weights (dashed), Shepard weights (dash-dotted) 10° 10" 10" 10" 10" 10" 10" 10" 10" -10 12 3 4 5 Figure 6: Pade-type barycentric interpolants with k = 7 ioi log(l -|- t)/t Let us now compare, for the degree k = 4, the Pade-type rational interpolant, the Pade-type barycentric interpolant, and the Pade approxiniant [4/4] which is given by ^3 , .4n [4/4]/(t) = (1680 + 840t + ISOr + 20r + r)/(1680 - 840t + 180r - 20t^ + T). 15 error rational (solid), barycentric (dash-dott) and Pade (dashed) Figure 7: Pade-type rational and barycentric interpolants, and Pade approximant for Let us remind that [4/4] j(t) — e* = O(t^), and that its construction makes use of the first 8 coefficients of the power series. The results are given in Figure [71 where the solid line represents the error of the Pade-type rational interpolant, the dashed line corresponds to the Pade approximant, and the dash-dotted line to the Pade-type barycentric interpolant. The interpolation points were chosen equidistant in the interval [0.1,0.8]. Notice that, for the interpolants, the error is smaller around the interpolation points, while the errors of the Pade approximant is more symmetric around the origin. Example 4: spurious pole removal Let us now give an example showing that the rational interpolant can have poles even if the function is continuous. In fact, it is known [23] that if, after cancelation of common factors between the numerator and the denominator and ordering the interpolation points, two consecutive weights Wi and Wj+i in the barycentric formula have the same sign, then the reduced interpolant has an odd number of poles in [tj, Tj+i). We consider the Pade-type rational interpolant of the cosine function with 5 equidis- tant interpolation points in the interval [— 7r/2, +7r/8]. As may be seen in Figure [HI the interpolant (dashed line) has one real pole at t = —2.8636 . . . (its other poles are com- plex). When t goes to infinity, the interpolant tends to 25.269 . . .. However, the results are quite good (the cosine function is the solid line in Figure [8]) from the right of the pole up to almost vr (see Figure [9]). It is possible to remove a spurious pole p by forcing the Pade-type interpolant to go through the point {p,f{p))- In Figure [HI the first of the equidistant interpolation 16 Pade-type rational interpolant for cosine function 4 3 2 1 -1 -2 -3 -4 -5 -4-202468 Figure 8: Pade-type rational interpolation of the cosine function points is replaced by the pole, a procedure which removes it and leads to a better result (dashed line). Error before (solid) and after (dasfied) pole removal for cosine function -202468 Figure 9: Error before (solid) and after (dashed) the pole removal for cost If the interpolant exhibits several poles, they can be eliminated successively. If a new pole is introduced during the procedure, then it can be removed similarly. 17 In our case, the location of the pole was directly computed from the coefficients bi of the denominator of the interpolant since they were available. It is also possible to locate approximately a pole when the absolute value of the interpolant becomes larger than a fixed threshold, or when the interpolant has a sudden change of sign, and then to impose it as an interpolation point. This procedure was tried on the Pade-type barycentric interpolant for cost in the same interval as before, but with k = 13. The interpolant was computed at 500 points in [— TT, +27r]. A sudden change of the sign was observed in the interval [5.3010, 5.3199]. We had i?i3(5.3010) = 15.068 and i?i3(5.3199) = -25.101. Replacing the first interpolation point To = —7t/2 by tq = 5.3050, the spurious pole was removed, and no other pole appeared. The advantage of this procedure is that it can also be be used for Pade-type barycen- tric interpolants where the coefficients of the denominator are not explicitly known. The same techniques can be applied in the case of partial Pade-type rational and barycentric interpolation. 6 Applications Let us now briefiy discuss some possible applications to numerical analysis problems. 6.1 Convergence acceleration We consider the sequence {Sn = /(Tn)) where (r„) is a sequence of parameters such that lim„_5.oo T„ = Too 7^ 0, ±oo, and where / is a function whose first coefficients of the series expansion around are known. We set S = lim„_j,oo Sn- The convergence of the sequence (Sn) can be accelerated by computing the Pade-type rational interpolant or the Pade-type barycentric interpolant Rj^ satisfying Rj^ (ri) = Si for i = 0, . . . ,k — 1 (or k in the second case), and setting T^ = i?^" (roo)-This is the essence of an extrapolation method for accelerating the convergence of a sequence |13] . Under certain assumptions, the sequences {Tjf' ) converge to 5* faster than (Sn) either when k is fixed and n goes to infinity or vice versa. 6.2 Inversion of the Laplace transform We consider the Laplace transform POO F{p)= / e-^'f{s)ds. Jo Assume that F is known at some points p„ for n = 0, ...,/;; — 1 (or k), and also the first coefficients of its series expansion around 0. F can be approximated by a Pade-type rational interpolant or by a Pade-type barycentric interpolant Rk, and the interpolant then inverted, thus leading to an approximation of /. Let us remark that, since limp_^oo F{p) = 0, the degree of the numerator of the interpolant must be smaller 18 than the degree of its denominator. The inversion can be performed without decompos- ing F into its partial fractions by a procedure due to Longman and Sharir [SO]- Let F have the form F{p) = with m < n. They showed that oo with Vi = Ui+rn + aiMi+m-l H h OmMi, Z = 0, 1 5-^5 where Ui = 0, Z = 0, ... ,77, — 2, Un-1 = 1, Mi = -{PlUi-i-\ h/3nMj-n), 2=n,n + l,... Usually, the series furnishing / is quickly converging. Let us take the example considered in [T3l p. 350] F{p) = log(l + aV/), /(s) = 2(1 - cosas)/s. We make the change of variable t = a^/p^, and we set F{p) = Git) = log(l + t)=t-^ + ^-... The Pade-type (rational and barycentric) interpolants will be approximations of G. Replacing t by a^/p^ in a Pade-type rational interpolant of degree k in t produces an interpolant of degree 2k in p, and we obtain an approximation of F of the form feoP^'^ + hia?p'^'^~'^ + ■ ■ ■ + hkO^^ Notice that, since cq = in the series expansion of G'(t), the relations ([T]) lead to ao = 0, and, thus, linip^oo R2k{p) = which is consistant with the asymptotic property of the Laplace transform. Thus, this approximant can be written as R2k{P) = A- p2k + p^p2k-2 + . . . + /3. 2k with A = a^tti/bo, a2i = a^*(aj+i/ai), for z = 1, . . . ,k — 1, f32i = a^*(6j/6o), for i = l,...,/c, the a's and the /3's with an odd index being zero. We see that the series expansion of R2k{p) only contains even powers of 1/p as the series F{p) itself. Inverting i?2fc by the procedure of Longman and Sharir [20] (after replacing m hy 2k — 2 and 19 n by 2k in the formulae for the Wj's and the Mj's), or performing its partial fraction decomposition, gives an approximation of /. For A; = 5, a = 1, Tj = 1/pf with pi = 0.1 +ih ioi i = 0, . . . ,k — l, and h = 2/{k — 1), the Pade-type rational interpolant leads to the results of Figure [101 using 12 terms in the series expansion of /. Although /(O) =0 and the series expansion by the method of Longman and Sharir is also at s = (since vq = 0), there is a loss of accuracy around this point due to the indeterminacy. These results have to be compared with those given in [T3| p. 350] which were obtained by constructing a rational interpolant with a numerator of degree 7 and a denominator of degree 8, that is using 16 interpolation points. We see that our Pade-type rational interpolant provides a much better precision. Moreover, the precision can be even improved by taking more terms in the series for / at almost no additional price. -1 -0.5 0.5 1 Figure 10: Error for the inversion of the Laplace transform Another treatment of this example could consist in making the change of variable t = a/p, thus leading to F{p) = G{t) = log(l + 1^) = t^ - 1^/2 + t^/3 . 6.3 Piecewise rational interpolation Our approach can be used for constructing piecewise rational interpolants. Let a < a' < < b' < b. We construct a first Pade-type rational or barycentric interpolant in [a, a'], and then a second one in [b',b]. Due to the Pade-type property of these interpolants and the fact that, for all i, Ci = f^^\0)/i\, the two interpolants and their first derivatives will have the same values at the point t = 0. Obviously, by a change of variable, the same construction holds at a point different from the origin, and it can be repeated. 20 One of the advantages of such a construction is to obtain a good accuracy with a low degree in the interpolants, thus avoiding the usually bad conditioning when using more interpolation points and a rational interpolant with a higher degree. We interpolate the function f{t) = log(l +t)/t on the intervals [—0.9,-0.1] and [+0.1, +1] with k = 2, which means that the first rational function interpolates / only at the points —0.9 and —0.1, and the second ones interpolates it at +0.1 and +1. These two interpolants and their first and second derivatives agree with that of / at t = 0. The solid line in Figure [11] corresponds to the curve formed by these two interpolants. The two systems have a condition number of 3.25 x 10^ and 2.86 x 10^, respectively. Then, we construct the Pade-type rational interpolant interpolating / at the 4 points —0.9, —0.1, +0.1 and +1, and with a 0{t^) error at the origin. The system is overdetermined since / > k, its condition number is 1.90 x 10^, and the error is given by the dashed line. Finally, with the same 4 interpolation points, we construct the interpolant of degree k = 3. The system is also overdetermined, its condition number is 6.37 X 10^^, and we obtain the results given by the dash-dotted line. 10" 10" 10" 10" 10" 10" 10" 10" 10 -1 -0.5 0.5 1 1.5 Figure 11: Piecewise Pade-type rational interpolant with /c = 2 for log(l + t)/t. 7 Conclusions In this paper, we presented in details the particular case of the general rational Hermite interpolation problem where only the value of the function is interpolated at some points, and where the function and its first derivatives agree at the origin. Thus, the interpolants constructed in this way possess a Pade-type property at 0. Numerical examples show the interest of the procedure. The interpolation procedure can be easily modified to 21 introduce a partial knowledge on the poles and the zeros of the function to approximated. We also showed how spurious poles can be eliminated. The ideas developed in this paper need additional investigations. An expression for the error will be useful. Another important open problem which remains to be studied is the convergence of the interpolants when the degree tends to infinity. In this direction, the results of Eiermann [T7] could be useful. Acknowledgment: We would like to thank Jean-Paul Berrut for interesting discus- sions and comments. This work was partially supported by MIUR, PRIN grant no. 20083KLJEZ-003, and by University of Padova, Project 2010 no. CPDA104492. References [4 [5 [6 [7; [8 [9 [10 ill G.A. Baker Jr., P.R. Graves-Morris, Pade Approximants, 2nd edition, Cambridge University Press, Cambridge, 1996. J. -P. Berrut, Rational functions for guaranteed and experimentally well- conditioned global interpolation, Comput. Math. AppL, 15 (1988) 1-16. J. -P. Berrut, R. Baltensperger, H.D. Mittelmann, Recent developments in barycen- tric rational interpolation, in Trends and Applications in Constructive Approxima- tion, M.G. de Bruin, D.H. Mache, J. Szabados eds., ISNM vol. 151, Birkhauser Verlag, Basel, 2005, pp. 27-51. J. -P. Berrut, H.D. Mittelmann, Lebesgue constant minimizing linear rational inter- polation of continuous functions over the interval, Comput. Math. AppL, 33 (1997) 77-86. J. -P. Berrut, H.D. Mittelmann, Matrices for the direct determination of the barycentric weights of rational interpolation, J. Comput. Appl. Math., 78 (1997) 355-370. C. Brezinski, Rational approximation to formal power series, J. Approx. Theory, 25 (1979) 295-317. C. Brezinski, Pade-Type Approximation and General Orthogonal Polynomials, ISNM, vol. 50, Birkhauser-Verlag, Basel, 1980. C. Brezinski, Partial Pade approximation, J. Approx. Theory, 54 (1988) 210-233. C. Brezinski, Biorthogonality and its Applications to Numerical Analysis, Marcel Dekker, New York, 1992. C. Brezinski, Computational Aspects of Linear Control, Kluwer, Dordrecht, 2002. C. Brezinski, Extrapolation algorithms for filtering series of functions, and treating the Gibbs phenomenon, Nunier. Algorithms, 36 (2004) 309-329. 22 [12] C Brezinski, Interpolation and implicitization of parametric curves, unpublished. [13] C Brezinski, M. Redivo-Zaglia, Extrapolation Methods. Theory and Practice, North-Holland, Amsterdam, 1991. [14] C Brezinski, J. Van Iseghem, Pade approximations, in Handbook of Numerical Analysis, vol. Ill, P.G. Ciarlet and J.L. Lions eds., North-Holland, Amsterdam, 1994, pp. 47-222. [15] C Brezinski, J. Van Iseghem, A taste of Fade approximation, in Acta Numerica 1995, A. Iserles ed., Cambridge University Press, Cambridge, 1995, pp. 53-103. [16] A. Cuyt, L. Wuytack, Nonlinear Methods in Numerical Analysis, North-Holland, Amsterdam, 1987. [17] M. Eiermann, On the convergence of Fade-type approximants to analytic functions, J. Comput. Appl. Math., 10 (1984) 219-227. [18] M.S. Floater, K. Hormann, Barycentric rational interpolation with no poles and high rates of approximation, Nunier. Math., 107 (2007) 315-331. [19] F.R. Graves-Morris, Efficient reliable rational interpolation, in Pade Approximation and its Applications. Amsterdam 1980, M.G. de Bruin and H. van Rossum eds., LNM vol. 888, Springer-Verlag, Berlin, Heidelberg, New York, 1981, pp. 28-63. [20] I. M. Longman, M. Sharir. Laplace transform inversion of rational functions, Geo- phys. J. R. Astron. Soc, 25 (1971) 299-305. [21] F. Gonnet, R. Fachon, L.N. Trefethen, Robust rational interpolation and least- squares, ETNA, to appear. [22] G. Hornecker, Approximations rationnelles voisines de la meilleur approximation au sens de Tchebyscheff, C. R. Acad. Sci. Faris, 249 (1959) 939-941. [23] G. Hornecker, Determination des meilleures approximations rationnelles (au sens de Tchebycheff), des fonctions reelles d'une variable sur un segment fini et des bornes d'erreur correspondantes, C. R. Acad. Sci. Faris, 249 (1959) 2265-2267. [24] S. Faszkowski, Approximation uniforme des fonctions continues par les fonctions rationnelles, Zastos. Mat., 6 (1963) 441-458. [25] C. Schneider, W. Werner, Some new aspects of rational interpolation. Math. Com- put., 47 (1986) 285-299. [26] C. Schneider, W. Werner, Hermite interpolation: the barycentric approach, Com- puting, 46 (1991) 35-51. [27] D. Shepard, A two-dimensional interpolation function for irregularly-spaced data, in Proceedings of the 23rd ACM National Conference, ACM Fress, New York, 1968, pp. 517—524. 23 [28] J. Van Deun, Eigenvalue problems to compute almost optimal points for rational interpolation with prescribed poles, Numer. Algorithms, 45 (2007) 89-99. [29] J. Van Deun, Computing near-best fixed pole rational interpolants, J. Comput. Appl. Math., 235 (2010) 1077—1084. [30] P. Wynn, Transformations to accelerate the convergence of Fourier series, in Gertrude Blanch Anniversary Volume, Wright Patterson Air Force Base, Aerospace Research Laboratories, Office of Aerospace Research, United States Air Force, 1967, pp. 339 379. 24