Skip to main content
Internet Archive's 25th Anniversary Logo

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