On Fast Implementation of Higher Order Hermite-Fejer Interpolation
aa r X i v : . [ m a t h . NA ] J un ON FAST IMPLEMENTATION OF HIGHER ORDER HERMITE-FEJ´ERINTERPOLATION ∗ SHUHUANG XIANG † AND GUO HE ‡† Abstract.
The problem of barycentric Hermite interpolation is highly susceptible to overflows or under-flows. In this paper, based on Sturm-Liouville equations for Jacobi orthogonal polynomials, we consider the fastimplementation on the second barycentric formula for higher order Hermite-Fej´er interpolation at Gauss-Jacobior Jacobi-Gauss-Lobatto pointsystems, where the barycentric weights can be efficiently evaluated and cost linearoperations corresponding to the number of grids totally. Furthermore, due to the division of the second barycen-tric form, the exponentially increasing common factor in the barycentric weights can be canceled, which yieldsa superiorly stable method for computing the simplified barycentric weights, and leads to a fast implementationof the higher order Hermite-Fej´er interpolation with linear operations on the number of grids. In addition, theconvergence rates are derived for Hermite-Fej´er interpolation at Gauss-Jacobi pointsystems.
Key words.
Hermite-Fej´er interpolation, barycentric, Jacobi polynomial, Gauss-Jacobi point, Lobatto-Gauss-Jacobi point, Chebyshev point.
AMS subject classifications.
1. Introduction.
There are many investigations for the behavior of continuous functionsapproximated by polynomials. Weierstrass [58] in 1885 proved the well known result that everycontinuous function f ( x ) in [ − ,
1] can be uniformly approximated as closely as desired bya polynomial function. This result has both practical and theoretical relevance, especially inpolynomial interpolation.Polynomial interpolation is a fundamental tool in many areas of numerical analysis. La-grange interpolation is a well known, classical technique for approximation of continuous func-tions. Let us denote by x ( n )1 , x ( n )2 , . . . , x ( n ) n (1.1)the n distinct points in the interval [ − ,
1] and let f ( x ) be a function defined in the same interval.The n th Lagrange interpolation polynomial of f ( x ) is uniquely defined by the formula L n [ f ] = n X k =1 f ( x ( n ) k ) ℓ ( n ) k ( x ) , ℓ ( n ) k ( x ) = ω n ( x ) ω ′ n ( x ( n ) k )( x − x ( n ) k ) , (1.2)where ω n ( x ) = ( x − x ( n )1 )( x − x ( n )2 ) · · · ( x − x ( n ) n ). However, for an arbitrarily given system ofpoints { x ( n )1 , x ( n )2 , . . . , x ( n ) n } ∞ n =1 , Bernstein [2] and Faber [13], in 1914, respectively, showed thatthere exists a continuous function f ( x ) in [ − ,
1] for which the sequence L n [ f ] ( n = 1 , , . . . ) isnot uniformly convergent to f in [ − , . Additionally, Bernstein [3] proved that there exists acontinuous function f ( x ) also for which the sequence L n [ f ] is divergent. Particularly, Gr¨unwald[23] in 1935 and Marcinkiewicz [32] in 1937, independently, showed that even for the Chebyshevpoints of first kind x ( n ) k = cos (cid:18) k − n π (cid:19) , k = 1 , , . . . , n, n = 1 , , . . . , (1.3)there is a continuous function f ( x ) in [ − ,
1] for which the sequence L n [ f ] is divergent everywherein [ − , ∗ This work was supported by National Science Foundation of China (No. 11371376). † Department of Applied Mathematics and Software, Central South University, Changsha, Hunan 410083, P.R. China. ‡ Corresponding author A very simple proof was given by Fej´er [15] in 1930.1
One of the proofs of Weierstrassapproximation theorem using interpolation polynomials was presented by Fej´er [14] in 1916 basedon the above Chebyshev pointsystem (1.3): If f ∈ C [ − , H n − ( f, x ) of degree at most 2 n − n →∞ k H n − ( f ) − f k ∞ = 0, where H n − ( f, x )is determined by H n − ( f, x ( n ) k ) = f ( x ( n ) k ) , H ′ n − ( f, x ( n ) k ) = 0 , k = 1 , , . . . , n. (1.4)This polynomial is known as the Hermite-Fej´er interpolation polynomial.The convergence result has been extended to general Hermite-Fej´er interpolation of f ( x ) atnodes (1.1), upon strongly normal pointsystems introduced by Fej´er [16]: Given, respectively,the function values f ( x ( n )1 ), f ( x ( n )2 ), . . . , f ( x ( n ) n ) and derivatives d ( n )1 , d ( n )2 , . . . , d ( n ) n at these grids,the Hermite-Fej´er interpolation polynomial H n − ( f ) has the form of H n − ( f, x ) = n X k =1 f ( x ( n ) k ) h ( n ) k ( x ) + n X k =1 d ( n ) k b ( n ) k ( x ) , (1.5)where h ( n ) k ( x ) = v ( n ) k ( x ) (cid:16) ℓ ( n ) k ( x ) (cid:17) , b ( n ) k ( x ) = ( x − x ( n ) k ) (cid:16) ℓ ( n ) k ( x ) (cid:17) and v ( n ) k ( x ) = 1 − ( x − x ( n ) k ) ω ′′ n ( x ( n ) k ) ω ′ n ( x ( n ) k ) (see Fej´er [17]).The pointsystem (1.1) is called strongly normal if for all nv ( n ) k ( x ) ≥ c > , k = 1 , , . . . , n, x ∈ [ − , c . The pointsystem (1.1) is called normal if for all nv ( n ) k ( x ) ≥ , k = 1 , , . . . , n, x ∈ [ − , . (1.7)Fej´er [16] (also see Szeg¨o [45, pp 339]) showed that for the zeros of Jacobi polynomial P ( α,β ) n ( x )of degree n ( α > − β > − v ( n ) k ( x ) ≥ min {− α, − β } for − < α ≤ − < β ≤ k = 1 , , . . . , n and x ∈ [ − , . While for the Legendre-Gauss-Lobatto pointsystem (the roots of (1 − x ) P (1 , n − ( x ) = 0 ), v ( n ) k ( x ) ≥ , k = 1 , , . . . , n, x ∈ [ − , . This result is extended to Jacobi-Gauss-Lobatto pointsystem (the roots of (1 − x ) P ( α,β ) n − = 0)and Jacobi-Gauss-Radau pointsystem (the roots of (1 − x ) P ( α,β ) n − = 0 or (1 + x ) P ( α,β ) n − = 0) byV´ertesi [52, 53]: for all k and x ∈ [ − , v ( n ) k ( x ) ≥ min { − α, − β } for { x ( n − k } S {− , } with 1 ≤ α ≤ ≤ β ≤ v ( n ) k ( x ) ≥ min { − α, − β } for { x ( n − k } S { } with 1 ≤ α ≤ − < β ≤ v ( n ) k ( x ) ≥ min {− α, − β } for { x ( n − k } S {− } with − < α ≤ ≤ β ≤ . Based upon the (strongly) normal pointsystem, Gr¨unwald [24] in 1942 showed that for every f ∈ C [ − , n →∞ k H n − ( f ) − f k ∞ = 0 if { x ( n ) k } is strongly normal satisfying (1.6) and { d ( n ) k } satisfies | d ( n ) k | < n c − δ for some given positive number δ, k = 1 , , . . . , n = 1 , , . . . , while lim n →∞ k H n − ( f ) − f k ∞ = 0 in [ − ǫ, − ǫ ] for each fixed 0 < ǫ < { x ( n ) k } is normaland { d ( n ) k } is uniformly bounded for n = 1 , , . . . .To get fast convergence on suitable smooth functions, higher order Hermite-Fej´er interpola-tion polynomials were considered in Goodenough and Mills [21], Sharma and Tzimbalario [39],Szabados [44], V´ertesi [55], etc.: for p = 0 , , . . . , m − q = 1 , , . . . , n , H mn − ( f, x ) = n X k =1 m − X j =0 f ( j ) ( x ( n ) k ) A jk ( x ) , H ( p ) mn − ( f, x ( n ) q ) = f ( p ) ( x ( n ) q ) , (1.8)where the polynomial A jk ( x ) of degree at most mn − A ( p ) jk ( x ( n ) q ) = δ jp δ kq , j = 0 , , . . . , m − , k = 1 , , . . . , n, (1.9)and δ is the Kronecker delta function. For simplicity, in the following we abbreviate x ( n ) k as x k , ℓ ( n ) k ( x ) as ℓ k ( x ), h ( n ) k ( x ) as h k ( x ), and b ( n ) k ( x ) as b k ( x ).The convergences of the higher order Hermite-Fej´er interpolation polynomials have beenextensively studied (see e.g. Byrne et al. [9], Goodenough and Mills [22], Locher [30], Mathurand Saxena [31], Moldovan [33], Nevai and V´ertesi [34], Popoviciua [35], [39], Shi [40, 41], Shishaet al. [42], Sun [43], Szili [46], Vecchia et al. [12], V´ertesi [52, 54] etc.). The convergence ratesare achieved most on Gauss-Jacobi or Jacobi-Gauss-Lobatto pointsystems. As is well known inapproximation theory, the right approach is to use point sets that are clustered at the endpointsof the interval with an asymptotic density proportional to (1 − x ) − / as n → ∞ emphasised,for example, in Berrut and Trefethen [5] and Trefethen [48]. Hence, in this paper we confineourselves to Gauss-Jacobi or Jacobi-Gauss-Lobatto pointsystems. In general, the Hermite interpolation is to find a polynomial H N − ( f, x ) of degree at most N − d r dx r H N − ( f, x ) (cid:12)(cid:12)(cid:12)(cid:12) x = x k = f k,r for r = 0 , , . . . , n k − , where f k, , f k, , . . . , f k,n k − denote the function value and its first n k − x k ( k = 1 , , . . . , n ), respectively, and N = n + n + · · · + n n .The polynomial H N − ( f, x ) can be represented in either the Newton form or the barycentricform. In the Newton form, the grid points x k must be ordered in a special way (see Schneiderand Warner [38]). If the grid points are not carefully ordered, the Newton form is susceptibleto catastrophic numerical instability. For more details, see Fischer and Reichel [18], Tal-Ezer[47], Berrut and Trefethen [5], Butcher et al. [8] and Sadiq and Viswanath [36]. In contrast, thebarycentric form does not depend on the order in which the nodes are arranged, which treats allthe grid points equally. Barycentric interpolation is arguably the method of choice for numericalpolynomial interpolation.The first barycentric formula for the Hermite interpolation is of the form of H N − ( f, x ) = H ∗ ( f, x ) P nk =1 f k,nk − ( n k − (cid:16) w k, x − x k (cid:17) + f k,nk − ( n k − (cid:16) w k, ( x − x k ) + w k, x − x k (cid:17) + · · · + f k, (cid:16) w k, ( x − x k ) nk + · · · + w k,nk − x − x k (cid:17) , (1.10)where H ∗ N ( f, x ) = Q nk =1 ( x − x k ) n k and w k,r is called the barycentric weights. Applying 1 ≡ H ∗ N ( f, x ) P nk =1 P n k − r =0 w k,r ( x − x k ) r − n k derives the second barycentric form H N − ( f, x ) = P nk =1 P n k − s =0 f k,s s ! P n k − s − r =0 w k,r ( x − x k ) r + s − n k P nk =1 P n k − r =0 w k,r ( x − x k ) r − n k (1.11) Two typos occur in (1.4) [36]: ( x − x k ) n k − r − s should be ( x − x k ) r + s − n k and ( x − x k ) n k − r be ( x − x k ) r − n k . (see e.g. [8, 36]).The second barycentric form is more robust in the presence of rounding errors in the weights w k,r . It is obvious from inspection that either of the two forms can be used to evaluate theinterpolant H N − ( f, x ) at a given point x using O ( P nk =1 n k ) arithmetic operations once thebarycentric weights w k,r are known.Schneider and Werner [38] used divided differences to evaluate the barycentric weights. Thismethod requires n ( n − / N − P nk =1 n k ) / Y j = k ( x + x k − x j ) − n j = Y j = k ( x k − x j ) − n j Y j = k (cid:18) − x − x k x k − x j (cid:19) − n j which can be obtained by the following recursion I k,r = r X s =1 P k,s I k,r − s /r, I k, = 1 , P k,r = X j = k n j ( x j − x k ) − r , r = 1 , , . . . , n k − , and then w k,r = C k I k,r with C k = Q j = k ( x k − x j ) − n j . It costs O ( n P nk =1 n k + P nk =1 n k )multiplications. Roughly half these operations are additions or subtractions and roughly half aremultiplications. Furthermore, if an additional derivative is prescribed at one of the interpolationpoints, update the barycentric coefficients use only O ( N ) operations [36].Notice that the barycentric Hermite interpolation problem is highly susceptible to overflowsor underflows. The weights w k,r in (1.10) and (1.11) usually vary by exponentially large factors.Figures 1.1-1.2 illustrate the magnitudes of the barycentric weights computed by the methodof Sadiq and Viswanath [36] at the Chebyshev pointsystem (1.3) or Legendre pointsystem withdifferent multiple number m of derivatives ( n = n = · · · = n n = m ). We can see from thesetwo figures that the barycentric weights become extremely large while the number of points andthe multiple number of derivatives are not so large, which will lead to overflows for larger n or m . Table 1.1 shows the threshold S that the algorithm suffers overflows for computation of theweights if n ≥ S with different m , respectively. l og | w k , r | m=2 w k,0 w k,1 l og | w k , r | m=3 w k,0 w k,1 w k,2 l og | w k , r | m=4 w k,0 w k,1 w k,2 w k,3 Fig. 1.1 . Magnitudes of the barycentric weights w k,r by method of Sadiq and Viswanath interpolating atChebyshev pointsystem (1.3). In Matlab , the largest positive normalized floating-point number in IEEE double precision is (1 + (1 − − ))2 ≈ . × , the smallest positive normalized floating-point number in IEEE double precisionis 2 − ≈ . × − . l og | w k , r | m=2 w k,0 w k,1 l og | w k , r | m=3 w k,0 w k,1 w k,2 l og | w k , r | m=4 w k,0 w k,1 w k,2 w k,3 Fig. 1.2 . Magnitudes of the barycentric weights w k,r by method of Sadiq and Viswanath interpolating atLegendre pointsystem. Table 1.1
The threshold S that the algorithm [36] collapses for computation of the barycentric weights if n ≥ S withdifferent m at the Chebyshev pointsystem (1.3) and Legendre pointsystem Chebyshev pointsystem (1.3) Legendre pointsystem m = 2 m = 3 m = 4 m = 2 m = 3 m = 4S=524 S=347 S=263 S=523 S=346 S=262 Thus, in [36], Sadiq and Viswanath used 2 cos((2 k − π/ (2 n )) instead of cos((2 k − π/ (2 n ))and considered Leja reordering of the points to get more stable computation of the weights anddecrease the chance of overflows or underflows. However, reordering of the points does notchange the magnitudes of the barycentric weights. Furthermore, Leja reordering needs O ( n )operations [10].Fortunately, from the second barycentric form (1.11), we see that the weights w k,r appear inthe denominator exactly as in the numerator. Due to the division, the barycentric weights canbe simplified by cancelling the common factors without altering the result (see Figures 2.1-2.4and Tables 2.1-2.4 below).In this paper, we are concerned with fast implementation of the higher order Hermite-Fej´erinterpolation polynomial (1.8), based on the second barycentric form (1.11), at Gauss-Jacobi orJacobi-Gauss-Lobatto pointsystems.Recently, a new algorithm on the evaluation of the nodes and weights for the Gauss quadra-ture was given by Glaser, Liu and Rokhlin [20] with O ( n ) operations, which has been extendedby both Bogaert, Michiels and Fostier [7], and Hale and Townsend [25]. A Matlab routine forcomputation of these nodes and weights can be found in
Chebfun system [51].As a result of these developments, in Section 2, we will discuss in details on calculation of thesecond barycentric weights of the higher order Hermite-Fej´er interpolation polynomial (1.11) atGauss-Jacobi or Jacobi-Gauss-Lobatto pointsystems, and present two algorithms with O ( nm )operations. Particularly, due to division, the common factor in the barycentric weights canbe canceled, which yields a superiorly stable method for computing the simplified barycentricweights. In Section 3, we will consider the stabilities of these implementations for the secondbarycentric formula on higher order Hermite-Fej´er interpolation and present numerical examplesillustrating the efficiency and accuracy. A final remark on the convergence rate of Hermite-Fej´erinterpolation (1.5) is included in Section 4.All the numerical results in this paper are carried out by using Matlab
R2012a on a desktop(2.8 GB RAM, 2 Core2 (32 bit) processors at 2.80 GHz) with Windows XP operating system.
2. Fast computation of the barycentric weights on higher order Hermite-Fej´erinterpolation.
In this section, we will introduce two methods for fast computation of thebarycentric weights on higher order Hermite-Fej´er interpolation at Gauss-Jacobi or Jacobi-Gauss-Lobatto pointsystems, which both share O ( nm ) operations and lead to fast and stablecalculation of the barycentric weights due to the division in (1.11), where the exponentiallyincreasing common factor is cancelled. Algorithm 1 : Following the barycentric Hermite interpolation formula in [8, 36, 38], werewrite the higher order Hermite-Fej´er interpolation (1.8) as H mn − ( f, x ) = H ∗ mn ( x ) n X k =1 f k,m − ( m − (cid:18) w k, x − x k (cid:19) + f k,m − ( m − (cid:18) w k, ( x − x k ) + w k, x − x k (cid:19) + · · · + f k, (cid:16) w k, ( x − x k ) m + · · · + w k,m − x − x k (cid:17) (2.1)with d r H mn − ( f, x ) dx r (cid:12)(cid:12)(cid:12)(cid:12) x = x k = f k,r for r = 0 , . . . , m − k = 1 , · · · , n ,(2.2)where H ∗ mn ( x ) = Q nk =1 ( x − x k ) m and f k,j = f ( j ) ( x k ). Denote by ¯ h k ( x ) = ω n ( x ) x − x k = ω ′ n ( x k ) ℓ k ( x ),then, expression (2.1) can be represented as H mn − ( f, x ) = n X k =1 f k, (cid:0) w k, ¯ h mk ( x ) + w k, ¯ h mk ( x )( x − x k ) + · · · + w k,m − ¯ h mk ( x )( x − x k ) m − (cid:1) + · · · + f k,m − ( m − w k, ¯ h mk ( x )( x − x k ) m − . (2.3)Furthermore, from (1.9) and (2.1)-(2.3), we have w k, ¯ h mk ( x ) (cid:12)(cid:12)(cid:12)(cid:12) x = x k = 1 ,w k, (¯ h mk ( x )) ( j ) + w k, (¯ h mk ( x )( x − x k )) ( j ) + · · · + w k,j (¯ h mk ( x )( x − x k ) j ) ( j ) (cid:12)(cid:12)(cid:12)(cid:12) x = x k = 0 (2.4)for j = 1 , . . . , m −
1. By using the Taylor expansion of ¯ h mk ( x ) at x k , it leads to (cid:0) ¯ h mk ( x )( x − x k ) i (cid:1) ( j ) j ! (cid:12)(cid:12)(cid:12)(cid:12) x = x k = (¯ h mk ( x )) ( j − i ) ( j − i )! (cid:12)(cid:12)(cid:12)(cid:12) x = x k , i = 0 , , . . . , j − . Thus, from (2.4), we get the following formulas w k, = 1( ω ′ n ( x k )) m , w k,j = − j − X i =0 w k,i ( j − i )! (cid:18) ℓ mk ( x ) (cid:19) ( j − i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x k , j = 1 , , . . . , m − . (2.5)In addition, from the definition of ℓ k ( x ) = ω n ( x ) ω ′ n ( x k )( x − x k ) , it is not difficult to deduce byapplying the Taylor expansion of ω n ( x ) at x = x k that( ℓ k ( x )) ( r ) (cid:12)(cid:12) x = x k r ! = ω ( r +1) n ( x k )( r + 1)! ω ′ n ( x k ) := M k,r , r = 0 , · · · , m − . (2.6)Set b k,j = (cid:2) ℓ mk ( x ) (cid:3) ( j ) x = xk j ! ( j = 0 , , . . . , m −
1) and a k,i = m ( i − (cid:2) ℓ ′ k ( x ) ℓ k ( x ) (cid:3) ( i − x = x k ( i = 1 , , . . . , m − b k,j = mj ! (cid:2) ℓ m − k ( x ) ℓ ′ k ( x ) (cid:3) ( j − x = x k = mj ! (cid:2) ℓ mk ( x ) ℓ ′ k ( x ) ℓ k ( x ) (cid:3) ( j − x = x k , it follows by Leibniz formulathat b k,j = 1 j j X i =1 a k,i b k,j − i , j = 1 , , . . . , m − b k, = 1 . (2.7)In the following, we shall show that a k,i can be computed from M k,r , then b k,j can beevaluated from the recursion (2.7). In fact, a k,i can be calculated from the coefficient of theTaylor expansion of ℓ ′ k ( x ) /ℓ k ( x ) at x = x k by the following lemma. Lemma 2.1.
Let F ( x ) = A ( x ) B ( x ) = P sj =1 F j ( x − x k ) j − + O (( x − x k ) s ) , where A ( x ) = P sj =1 A j ( x − x k ) j − + O (( x − x k ) s ) and B ( x ) = P sj =1 B j ( x − x k ) j − + O (( x − x k ) s ) , then it obtains B F j = A j − j X i =2 B i F j − i +1 , j = 1 , · · · , s. (2.8) Proof . The proof is trivial.Since ℓ k ( x ) = m − X r =1 M k,r − ( x − x k ) r − + O (( x − x k ) m − ) , (2.9)and ℓ ′ k ( x ) = m − X r =1 rM k,r ( x − x k ) r − + O (( x − x k ) m − ) , (2.10)by Lemma 2.1 and using M k, = 1 and ℓ ′ k ( x ) ℓ k ( x ) = m P m − r =1 a k,r ( x − x k ) r − + O (( x − x k ) m − ), wehave a k, = mM k, , a k,i = imM k,i − i X j =2 M k,j − a k,i − j +1 , i = 2 , . . . , m − , (2.11)and then by (2.5) we get w k, = 1( ω ′ ( x k )) m , w k,j = − j − X i =0 w k,i b k,j − i , j = 1 , . . . , m − . (2.12)Thus, if { M k,r } m − r =0 and w k, are known, from (2.11) the total computation of { a k,j } m − j =0 costs O ( m ) operations, the same as those for { b k,j } m − j =0 and { w k,j } m − j =0 by (2.7) and (2.12), respec-tively, and then the implementation of (1.11) with n = · · · = n n = m costs O ( nm ) operations. Algorithm 2 : The barycentric weights w k,r of the Hermite-Fej´er interpolation (2.1) canalso be calculated by another fast way based on the formula given in Szabados [44], H mn − ( f, x ) = n X k =1 m − X j =0 f ( j ) ( x k ) ℓ k ( x ) m j ! m − j − X i =0 [ ℓ k ( x ) − m ] ( i ) x = x k i ! ( x − x k ) i + j . (2.13)We rewrite the interpolation (2.13) in the first barycentric interpolation form H mn − ( f, x ) = H ∗ mn ( x ) n X k =1 ω ′ n ( x k )) m m − X j =0 f k,j j ! m − j − X i =0 (cid:2) ℓ k ( x ) − m (cid:3) ( i ) x = x k i ! ( x − x k ) i + j − m , (2.14)and second barycentric interpolation form H mn − ( f, x ) = P nk =1 w k, P m − j =0 f k,j j ! P m − j − i =0 (cid:2) ℓ k ( x ) − m (cid:3) ( i ) x = xk i ! ( x − x k ) i + j − m P nk =1 w k, P m − i =0 (cid:2) ℓ k ( x ) − m (cid:3) ( i ) x = xk i ! ( x − x k ) i − m (2.15)respectively. Next we concentrate on the computation of ˜ b k,j := (cid:2) ℓ k ( x ) − m (cid:3) ( j ) x = xk j ! ( j = 0 , , . . . , m − w k,j = ˜ b k,j w k, , j = 0 , , . . . , m − . (2.16) From [44], it follows˜ b k,j = 1 j j X i =1 ˜ a k,i ˜ b k,j − i , j = 1 , , . . . , m −
1; ˜ b k, = 1 , (2.17)where ˜ a k,i = m ( i − (cid:2) x − x k − ω ′ n ( x ) ω n ( x ) (cid:3) ( i − x = x k ( i = 1 , · · · , m − ω n ( x ) = ω ′ n ( x k ) ℓ k ( x )( x − x k ), it is easy to see that ˜ a k,i = m ( i − (cid:2) ℓ k ( x ) ω ′ n ( x k ) − ω ′ n ( x ) ω n ( x ) (cid:3) ( i − x = x k = m ( i − (cid:2) − ℓ ′ k ( x ) ℓ k ( x ) (cid:3) ( i − x = x k . Similarly,by Lemma 2.1 and (2.9)-(2.10), we get˜ a k,j = − mM k, , ˜ a k,j = − jmM k,j − j X i =2 M k,j − ˜ a k,j − i +1 , j = 2 , . . . , m − . (2.18)Thus, from (2.18) and (2.17), the total computation of barycentric weights { e b k,r } m − r =0 costs also O ( m ) operations if { M k,r } m − r =0 and w k, are known.From the above illustrations, we see that fast computation of M k,r and w k, leads to fastimplementation of higher order barycentric Hermite-Fejer interpolation. We shall show thatfor each k , { M k,r } m − r =0 can be rapidly calculated with O ( m ) operations for Gauss-Jacobi orJacobi-Gauss-Lobatto pointsystems. Let { x k } nk =1 be the zeros of the Jacobi polynomial P ( α,β ) n ( x ). Thus ω n ( x ) = P ( α,β ) n ( x ) K n , where K n is the leading coefficient of P ( α,β ) n ( x ). From (2.6),we get M k,r − = ω ( r ) n ( x k ) r ! ω ′ n ( x k ) = d r P ( α,β ) n dx r ( x k ) r ! P ′ ( α,β ) n ( x k ) , r = 1 , . . . , m. It is known that P ( α,β ) n ( x ) is the unique solution of the second order linear homogeneous Sturm-Liouville differential equation(1 − x ) y ′′ + ( β − α − ( α + β + 2) x ) y ′ + n ( n + α + β + 1) y = 0 , (2.19)from which it is not difficult to deduce that(1 − x ) y ( r +2) + [ β − α − ( α + β + 2( r + 1)) x ] y ( r +1) + [ n ( n + α + β + 1) − r ( α + β + r + 1)] y ( r ) = 0 , (2.20)for r = 0 , , . . . . Thus, we have M k,r +1 = ( α + β +2( r +1)) x k + α − β − x k r +2) M k,r + r ( α + β + r +1) − n ( n + α + β +1)1 − x k r +2)( r +1) M k,r − , (2.21)with M k, = 1, M k, = α − β +( α + β +2) x k − x k ) and then { M k,r } m − r =0 can be computed in O ( m ). Con-sequently, from (2.6)-(2.7), (2.11)-(2.12), (2.14)-(2.15) and (2.18), the barycentric form (1.11)with n = · · · = n n = m and (2.15) can be achieved in O ( nm ) operations if w k, is known.From Wang et al. [57], w k, has the explicit form w k, = (cid:20) C ( α,β ) n ( − k +1 q (1 − x k ) w k (cid:21) m , k = 1 , , . . . , n, (2.22)where w k is the Gaussian quadrature weight corresponding to x k for the Jacobi weight, C ( α,β ) n = σ n Γ(2 n + α + β + 1)2 n + α + β +12 p n !Γ( n + α + β + 1)Γ( n + α + 1)Γ( n + β + 1)and σ n = +1 for n odd and σ = − n even. Moreover, both { x k } nk =1 and { w k } nk =1 can beefficiently calculated by routine jacpts in Chebfun [51] with O ( n ) operations.Additionally, it is worth noting that due to the division of the second barycentric form (1.11)or formula (2.15), the common factor ( C ( α,β ) n ) m of weights w k,r from (2.12) and (2.16) can becancelled without affecting the value of H mn − ( f, x ) (we call these new weights as simplifiedbarycentric weights). Then the barycentric weight w k, can be simplified as w k, = ( − m ( k +1) (cid:20)q (1 − x k ) w k (cid:21) m . (2.23)Comparing these two algorithms with Sadiq and Viswanath’s [36], we find that the newalgorithms cost O ( nm ) operations much less than that O ( n m + nm ) given by Sadiq andViswanath [36] if n ≫ m . Moreover, due to the cancellation of the common factor ( C ( α,β ) n ) m ,the computation of { w k,j } m − j =0 is quite efficient and stable (see Figures 2.1-2.2 and Tables 2.3-2.4). Tables 2.1-2.2 show the values of the common factor (cid:16) C ( α,β ) n (cid:17) m with respect to Chebyshevpointsystem (1.3) and Gauss-Legendre pointsystem, respectively. Table 2.1 (cid:16) C ( α,β ) n (cid:17) m with respect to Chebyshev pointsystem (1.3): α = β = − . m = 2 m = 3 m = 4 m = 10 n = 100 1 . ∗ − . ∗ . ∗ . ∗ n = 200 1 . ∗ − . ∗ . ∗ . ∗ n = 500 1 . ∗ − . ∗ . ∗ . ∗ n = 1000 9 . ∗ − . ∗ . ∗ . ∗ Table 2.2 (cid:16) C ( α,β ) n (cid:17) m with respect to Gauss-Legendre pointsystem: α = β = 0 m = 2 m = 3 m = 4 m = 10 n = 100 2 . ∗ − . ∗ . ∗ . ∗ n = 200 2 . ∗ − . ∗ . ∗ . ∗ n = 500 3 . ∗ − . ∗ . ∗ . ∗ n = 1000 1 . ∗ − . ∗ . ∗ . ∗ −20−15−10−50 k l og | w k , r | m=2 w k,0 wk,1 −30−20−10010 k l og | w k , r | m=3 k=0k=1k=2 −40−30−20−10010 k l og | w k , r | m=4 k=0k=1k=2k=3 Fig. 2.1 . Magnitude of the simplified barycentric weights w k,r by the Algorithm 2 interpolating at theChebyshev pointsystem (1.3): k = 1 : 10 and r = 0 : m − . Remark 2.2.
The simplified barycentric weights (2.23), in the case m = 1 ( w k, =( − k +1 p (1 − x k ) w k ), is exact the barycentric weight for the barycentric formula of Lagrangeinterpolation at the Jacobi pointsystem { x k } nk =1 , which was derived for Chebyshev points of sec-ond kind in Henrici [27], for Legendre points in Wang and Xiang [56], and extended to Jacobipoints in Hale and Trefethen [26]. −25−20−15−10−5 k l og | w k , r | m=2 w k,0 wk,1 0 2 4 6 8 10x 10 −40−30−20−10010 k l og | w k , r | m=3 k=0k=1k=2 −50−40−30−20−10010 k l og | w k , r | m=4 k=0k=1k=2k=3 Fig. 2.2 . Magnitude of the simplified barycentric weights w k,r by the Algorithm 2 interpolating at Leg-endre pointsystem: k = 1 : 10 and r = 0 : m − . Table 2.3
The CPU time for computation of the simplified barycentric weights w k,r by the Algorithm 2 at theChebyshev pointsystem (1.3): k = 1 : n and r = 0 : m − . m = 2 m = 3 m = 4 m = 10 n = 10 n = 10 n = 10 Table 2.4
The CPU time for computation of the simplified barycentric weights w k,r by the Algorithm 2 at theLegendre pointsystem: k = 1 : n and r = 0 : m − . m = 2 m = 3 m = 4 m = 10 n = 10 n = 10 n = 10 Suppose − x < x < · · · < x n − 1, it follows M ,r − = M r − ( − 1) = (cid:0) P ( α,β ) n − (cid:1) ( r − ( − r − P ( α,β ) n − ( − − (cid:0) P ( α,β ) n − (cid:1) ( r − ( − r − P ( α,β ) n − ( − , r = 2 , . . . , m, (2.25) M n,r − = M r − (1) = (cid:0) P ( α,β ) n − (cid:1) ( r − (1)( r − P ( α,β ) n − (1) + (cid:0) P ( α,β ) n − (cid:1) ( r − (1)( r − P ( α,β ) n − (1) , r = 2 , . . . , m, (2.26)and for the case x k = ± r = 2 , . . . , m , M k,r − = M r − ( x k ) = (cid:0) P ( α,β ) n − (cid:1) ( r ) ( x k ) r ! (cid:0) P ( α,β ) n − (cid:1) ′ ( x k ) + x k x k − (cid:0) P ( α,β ) n − (cid:1) ( r − ( x k )( r − (cid:0) P ( α,β ) n − (cid:1) ′ ( x k ) + x k − (cid:0) P ( α,β ) n − (cid:1) ( r − ( x k )( r − (cid:0) P ( α,β ) n − (cid:1) ′ ( x k ) , (2.27)where (cid:0) P ( α,β ) n − (cid:1) ( r ) ( x k ) r ! (cid:0) P ( α,β ) n − (cid:1) ′ ( x k ) can be evaluated by (2.21) with n − n for x k = ± 1, and for x k = ± (cid:16) P ( α,β ) n − (cid:17) ( r ) ( ± r ! (cid:16) P ( α,β ) n − (cid:17) ( ± 1) = ( n − α + β + n − − ( r − α + β + r ) r ( ± ( α + β + 2 r ) + α − β ) (cid:16) P ( α,β ) n − (cid:17) ( r − ( ± r − (cid:16) P ( α,β ) n − (cid:17) ( ± r = 1 , , . . . , m − w k, has the explicit form of w k, = h C ( α,β ) n − ( − k +1 p η k b w k i m , η k = β, k = 1 ,α, k = n, , otherwise , (2.28)where b w k = w k − x k for k = 2 , , . . . , n − b w = 2 α + β − Γ( β )Γ( β + 1)Γ( n + α + 1) n !Γ( n + β + 1)Γ( n + α + β + 1) , b w n = 2 α + β − Γ( α )Γ( α + 1)Γ( n + β + 1) n !Γ( n + α + 1)Γ( n + α + β + 1)and { w k } n − k =2 is the Gaussian quadrature weight corresponding to { x k } n − k =2 . Due to the divisionin (2.15), the barycentric weight w k, can be simplified as w k, = ( − m ( k +1) ( η k b w k ) m/ . (2.29)Figures 2.3-2.4 show magnitude of the simplified barycentric weights w k,r by the Algorithm 2 interpolating at the Jacobi-Gauss-Lobatto pointsystem with α = β = 1 . n = 10 and m = 2 , , −20−15−10−50 k l og | w k , r | m=2 w k,0 w k,1 −30−20−10010 k l og | w k , r | m=3 w k,0 w k,1 w k,2 −40−30−20−10010 k l og | w k , r | m=4 w k,0 w k,1 w k,2 w k,3 Fig. 2.3 . Magnitude of the simplified barycentric weights w k,r by the Algorithm 2 interpolating at theJacobi-Gauss-Lobatto pointsystem with α = β = 1 . : k = 1 : 10 and r = 0 : m − . Fig. 2.4 . Magnitude of the simplified barycentric weights w k,r by the Algorithm 2 interpolating atLegendre-Lobatto pointsystem: k = 1 : 10 and r = 0 : m − . w k, = 0 except the first and the last term. Remark 2.3. In the case m = 1 and α = β = 0 . , the nodes { x k } nk =1 are the Chebyshevpoints of second kinds x k = cos (cid:18) k − n − π (cid:19) , k = 1 , , . . . , n (2.30) and w k, = ( − k +1 η k with η k = (cid:26) , k = 1 , n , otherwise derived by Salzer [37]. Algorithm 1 Barycentric weights at Gauss-Jacobi or Jacobi-Gauss-Lobatto pointsystems Input parameters n , m , α , β . Compute the nodes { x k } nk =1 and simplified barycentric weights { w k, } nk =1 by jacpts . Compute M k,r = ω ( r +1) n ( x k )( r +1)! ω ′ n ( x k ) ( k = 1 , · · · , n, r = 1 , · · · , m − 1) by recursion. Compute a k,i = imM k,i − P i − j =1 a k,j M k,i − j ( k = 1 , · · · , n, i = 1 , · · · , m − Compute b k,i = i P iv =1 a k,v b k,i − v ( k = 1 , · · · , n, i = 1 , · · · , m − 1) with b k, = 1. Let c k, = 1, compute c k,i = − P i − j =1 c k,j b k,i − j ( k = 1 , · · · , n, i = 1 , · · · , m − Return w k,i = w k, c k,i ( k = 1 , · · · , n, i = 0 , · · · , m − Algorithm 2 Barycentric weights at Gauss-Jacobi or Jacobi-Gauss-Lobatto pointsystems Input parameters n , m , α , β . Compute the nodes { x k } nk =1 and simplified barycentric weights { w k, } nk =1 by jacpts . Compute M k,r = ω ( r +1) n ( x k )( r +1)! ω ′ n ( x k ) ( k = 1 , · · · , n, r = 1 , · · · , m ) by recursion. Compute ˜ a k,i = − imM k,i − P i − j =1 ˜ a k,j M k,i − j ( k = 1 , · · · , n, i = 1 , · · · , m − Compute ˜ b k,i = i P iv =1 ˜ a k,v ˜ b k,i − v ( k = 1 , · · · , n, i = 1 , · · · , m − 1) with ˜ b k, = 1. Return w k,i = w k, ˜ b k,i ( k = 1 , · · · , n, i = 0 , · · · , m − In particular, from (2.5), the barycen-tric weights for lower order Hermite-Fej´er barycentric interpolation can be given in the explicitforms. • m = 2: w k, = − ω ′′ n ( x k ) ω ′ n ( x k ) w k, . Moreover, for the Gauss-Jacobi pointsystem, w k, = ( β − α − ( α + β + 2) x k ) w k , , k = 1 , , . . . , n, while for the Jacobi-Gauss-Lobatto pointsystem, w , = (cid:18) n − n + α + β − β + 1 (cid:19) β b w , w n, = (cid:18) − − ( n − n + α + β − α + 1 (cid:19) α b w n and w k, = ( β − α − ( α + β − x k ) b w k − x k , k = 2 , . . . , n − . • m = 3: w k, = − ω ′′ n ( x k )2 ω ′ n ( x k ) w k, , w k, = − ω (3) n ( x k )2 ω ′ n ( x k ) + 32 (cid:18) ω ′′ n ( x k ) ω ′ n ( x k ) (cid:19) ! w k, . • m = 4: w k, = − ω ′′ n ( x k ) ω ′ n ( x k ) w k, , w k, = − ω (3) n ( x k )3 ω ′ n ( x k ) + 52 (cid:18) ω ′′ n ( x k ) ω ′ n ( x k ) (cid:19) ! w k, and w k, = − ω (4) n ( x k )6 ω ′ n ( x k ) + 5 ω ′′ n ( x k ) ω (3) n ( x k )3 ω ′ n ( x k ) ω ′ n ( x k ) − (cid:18) ω ′′ n ( x k ) ω ′ n ( x k ) (cid:19) ! w k, . 3. Illustration of numerical stability and numerical examples. The stability for thesecond barycentric formulas for Lagrange interpolation has been extensively studied by Henrici[27], Berrut and Trefethen [5]. Rigorous arguments that make this intuitive idea precise areprovided by Higham [28, 29]. For more details, see [5, 28, 29].3These arguments can be directly applied to the second barycentric higher order Hermite-Fej´er formula (2.15) if the barycentric weights can be evaluated well and do not suffer fromoverflows or underflows, which makes barycentric interpolation entirely reliable in practice forsmall values of m at the pointsystems, such as the Chebyshev pointsystem (1.3) and Jacobi-Gauss-Lobatto pointsystem for α = β = 1 . 5, i.e. the roots of (1 − x ) P ( , ) n − ( x ), studied in thispaper.To be pointed out especially, although both the Algorithm 1 and Algorithm 2 enjoy thefast implementation at the cost of O ( nm ) operations for pointsystems discussed in Section 2,the performances are not always same. More specifically, they own the same high accuracy forsmall m and different outcomes for larger m . The Algorithm 2 manifests better stability forlarger m in our numerical experiments. Form the descriptions of the two algorithms, we can seethat the first algorithm needs one more step than the second algorithm. This extra step leadsto a great loss of significance due to the fast growth of the entries in these two algorithms forlarge m . So, the Algorithm 2 is recommended in practice.Here, we use four functions to test the accuracy and stability for the second Hermite-Fej´er barycentric interpolation form with Algorithm 2 at the Gauss-Jacobi and Jacobi-Gauss-Lobatto pointsystems with f ( x ) = 1 / (1+ x ), which is analytic in a neighborhood of [ − , C ∞ function f ( x ) = e − /x , and nonsmooth functions f ( x ) = 1 − | x | . Figures 3.1-3.3 illustrate theperformance of the second barycentric interpolant (2.15) for the above first three functions atChebyshev pointsystem (1.3) and Jacobi-Gauss-Lobatto pointsystem with α = β = 1 . n = 10 : 20 : 2000 grid points and m = 2 j − ( j = 1 , , . . . , 6) derivatives under the ∞ -norm forthe vector k f ( x ) − H mn − ( f, x ) k ∞ at x = − . 02 : 1. −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=1 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=2 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=410 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=8 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=16 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=32Jacobi: α = β =−0.5Lobatto: α = β =1.5 Fig. 3.1 . k f ( x ) − H mn − ( f, x ) k ∞ at x = − . 02 : 1 with n = 10 : 20 : 2000 and m = 2 j − ( j = 1 , , . . . , )for f ( x ) = 1 / (1 + x ) at the Chebyshev pointsystem (1.3) and Jacobi-Gauss-Lobatto pointsystem with α = β =1 . , respectively. From these examples, we can see that the barycentric interpolation is quite stable for smallvalues of m . However, when m is too large, the simplified barycentric weights will suffer fromoverflows or underflows too. Figure 3.5 shows the maximum m for a fixed n in the computation4 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=1 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=2 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=410 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=8 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=16 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=32Jacobi: α = β =−0.5Lobatto: α = β =1.5 Fig. 3.2 . k f ( x ) − H mn − ( f, x ) k ∞ at x = − . 02 : 1 with n = 10 : 20 : 2000 and m = 2 j − ( j = 1 , , . . . , )for f ( x ) = e − /x at the Chebyshev pointsystem (1.3) and Jacobi-Gauss-Lobatto pointsystem with α = β = 1 . ,respectively. −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=1 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=2 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=410 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=8 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=16 10 −20 −15 −10 −5 n=10(2k−1), k=1:100 || f ( x ) − H m n − ( f, x ) || ∞ a t x =− : . : m=32Jacobi: α = β =−0.5Lobatto: α = β =1.5 Fig. 3.3 . k f ( x ) − H mn − ( f, x ) k ∞ at x = − . 02 : 1 with n = 10 : 20 : 2000 and m = 2 j − ( j = 1 , , . . . , )for f ( x ) = 1 − | x | at the Chebyshev pointsystem (1.3) and Jacobi-Gauss-Lobatto pointsystem with α = β = 1 . ,respectively. w k,r ( k = 1 : n, r = 0 : m − 1) by Algorithm 2 before theoverflows or underflows occurre. m m Fig. 3.4 . The maximum number m for a fixed n before the overflows or underflows occurred for Gauss-Jacobi pointsystem α = β = − . (left) and for Jacobi-Gauss-Lobatto pointsystem α = β = 1 . (right): n = 10 :10 . 4. Final remarks. It is remarkable that Chebyshev pointsystems (1.3) and (2.32) arefairly nice in Lagrange polynomial approximation (see [49, 50, 59]). However, for (higher or-der) Hermite-Fej´er interpolation, the Chebyshev pointsystem (2.32) completely fails (see Figure4.1). The good choice is Chebyshev pointsystem (1.3) or the roots of (1 − x ) P ( , ) n − ( x ), sincepointsystem (2.32) is not normal and the latter two pointsystems are strongly normal for m = 2. −1 −0.5 0 0.5 10123456 x 10 −11 m=2, n=10 x=−1:0.002:1 A b s o l u t e e rr o r s −1 −0.5 0 0.5 100.20.40.60.81 x 10 −8 m=2, n=10 x=−1:0.002:1 A b s o l u t e e rr o r s −1 −0.5 0 0.5 100.511.522.533.5 x 10 −7 m=2, n=10 x=−1:0.002:1 A b s o l u t e e rr o r s −1 −0.5 0 0.5 10123456 x 10 −7 m=3, n=10 x=−1:0.002:1 A b s o l u t e e rr o r s −1 −0.5 0 0.5 100.511.522.533.54 x 10 −3 m=4, n=10 x=−1:0.002:1 A b s o l u t e e rr o r s −1 −0.5 0 0.5 100.10.20.30.40.50.60.7 m=8, n=10 x=−1:0.002:1 A b s o l u t e e rr o r s Fig. 4.1 . The absolute errors of H mn − ( f, x ) − f ( x ) at x = − . 002 : 1 with different m and n by usingChebyshev pointsystem (2.32) for f ( x ) = x . For strongly normal pointsystem satisfying (1.6), V´ertesi [52] proved that for each f ∈ C [ − , k E ( f ) k ∞ = max x ∈ [ − , | H n − ( f, x ) − f ( x ) | ≤ (cid:18) c (cid:19) min q j ∈P n − k f ′ − q j k ∞ where P n − denotes the set of all polynomials of degree at most 2 n − f is analytic or of finite limited regularity, the convergence rate on Hermite-Fej´er inter-polation H n − ( f, x ) at Gauss-Jacobi pointsystem can be improved and given explicitely basedon the asymptotics of the coefficients of Chebyshev series for f .6 Suppose f ( x ) satisfies a Dini-Lipschitz condition on [ − , f ( x ) = ∞ X j =0 ′ c j T j ( x ) , c j = 2 π Z − f ( x ) T j ( x ) √ − x dx, j = 0 , , . . . . (4.1)where the prime denotes summation whose first term is halved, T j ( x ) = cos( j cos − x ) denotesthe Chebyshev polynomial of degree j . Lemma 4.1. (i) (Bernstein [4]) If f is analytic with | f ( z ) | ≤ M in the region bounded bythe ellipse E ρ with foci ± and major and minor semiaxis lengths summing to ρ > , then foreach j ≥ , | c j | ≤ Mρ j . (4.2) (ii) (Trefethen [49, 50]) For an integer k ≥ , if f ( x ) has an absolutely continuous ( k − stderivative f ( k − on [ − , and a k th derivative f ( k ) of bounded variation V k = Var( f ( k ) ) < ∞ ,then for each j ≥ k + 1 , | c j | ≤ V k πj ( j − · · · ( j − k ) . (4.3) Lemma 4.2. Suppose { x j } nj =1 are the roots of P ( α,β ) n ( x ) ( α, β > − ), then it follows ( x − x j ) ℓ j ( x ) = σ n q (1 − x j ) w j ( α + β +1) / s n !Γ( n + α + β + 1)Γ( n + α + 1)Γ( n + β + 1) P ( α,β ) n ( x ) , j = 1 , , . . . , n. (4.4) Proof . Let z n = R − (1 − x ) α (1+ x ) β [ P ( α,β ) n ( x )] dx and K n the leading coefficient of P ( α,β ) n ( x ).From Abramowitz and Stegun [1], we have z n = 2 α + β +1 n + α + β + 1 · Γ( n + α + 1)Γ( n + β + 1) n !Γ( n + α + β + 1) , K n = 12 n Γ(2 n + α + β + 1) n !Γ( n + α + β + 1) . Furthermore, by Hale and Townsend [25] and Wang et al. [57], we obtain( x − x j ) ℓ j ( x ) = w n ( x ) w ′ n ( x j ) = σ n ( − j r K n n (1 − x j ) w j n (2 n + α + β +1) z n w n ( x )= σ n ( − j r (1 − x j ) w j (2 n + α + β +1) z n P ( α,β ) n ( x ) , which leads to the desired result (4.4). Theorem 4.3. Suppose { x j } nj =1 are the roots of P ( α,β ) n ( x ) ( − < α, β ≤ ), then theHermite-Fej´er interpolation (1.5) at { x j } nj =1 has the convergence rate k E ( f ) k ∞ ≤ τ n M [2 nρ + (1 − n ) ρ ]( ρ − ρ n ( n ≥ , if f analytic in E ρ with | f ( z ) | ≤ M τ n V k ( k − π (2 n − n − · · · (2 n − k + 1) , if f, . . . , f ( k − absolutely continuousand V k < ∞ , n ≥ k/ , k ≥ , (4.5) where E ( f, x ) = f ( x ) − H n − ( f, x ) , and τ n = O ( n − . − min { α,β } log n ) , if − < min { α, β } ≤ max { α, β } ≤ − O ( n { α,β }− min { α,β }− ) , if − < min { α, β } ≤ − < max { α, β } ≤ O ( n { α,β } ) , if − < min { α, β } ≤ max { α, β } ≤ . (4.6)7 Proof . Since the Chebyshev series expansion of f ( x ) is uniformly convergent under theassumptions of Theorem 4.3, and the error of Hermite-Fej´er interpolation (1.5) on Chebyshevpolynomials satisfies | E ( T j , x ) | = | T j ( x ) − H n − ( T j , x ) | = 0 for j = 0 , , . . . , n − 1, then ityields | E ( f, x ) | = | f ( x ) − H n − ( f, x ) | = | ∞ X j =0 | c j || E ( T j , x ) | ≤ ∞ X j =2 n | c j || E ( T j , x ) | . (4.7)Furthermore, | E ( T j , x ) | = | T j ( x ) − P ni =1 T j ( x i ) h i ( x ) − P ni =1 T ′ j ( x i ) b i ( x ) | . In the following, wewill fucus on estimates on | E ( T j , x ) | for j ≥ n .Notice that the pointsystem is normal which implies h i ( x ) ≥ i = 1 , , . . . , n and x ∈ [ − , ≡ n X i =1 h i ( x ) = n X i =1 v i ( x ) ℓ i ( x )(see [16]) and then | n X i =1 T j ( x i ) h i ( x ) | ≤ n X i =1 h i ( x ) = 1 , j = 0 , , . . . . (4.8)Additionally, by Lemma 4.2, it obtains for j = 2 n, n + 1 , . . . that | P ni =1 T ′ j ( x i ) b i ( x ) | = j | P ni =1 U j − ( x i )( x − x i ) ℓ i ( x ) | = j ( α + β +1) / q n !Γ( n + α + β +1)Γ( n + α +1)Γ( n + β +1) | P ( α,β ) n ( x ) P ni =1 U j − ( x i ) p (1 − x i ) w i ℓ i ( x ) | = j ( α + β +1) / q n !Γ( n + α + β +1)Γ( n + α +1)Γ( n + β +1) | P ( α,β ) n ( x ) P ni =1 sin(( j − 1) arccos( x i )) √ w i ℓ i ( x ) | = jO (cid:16) | P ( α,β ) n ( x ) | p k{ w i } ni =1 k ∞ Λ n (cid:17) since q n !Γ( n + α + β +1)Γ( n + α +1)Γ( n + β +1) is decreasing as n increases and then uniformly bounded on n for − < α, β ≤ 0, where Λ n = max x ∈ [ − , P ni =1 | ℓ i ( x ) | is the Lebesgue constant, which, togetherwith P ( α,β ) n ( x ) = (cid:26) O ( n − ) , if max { α, β } ≤ − O ( n max { α,β } ) , if max { α, β } > − , w i = (cid:26) O ( n − − { α,β } ) , if min { α, β } ≤ − O ( n − ) , if min { α, β } > − (see Szeg¨o [45, pp 168, 354]) andΛ n = (cid:26) O (log n ) , if max { α, β } ≤ − O ( n max { α,β } + ) , if max { α, β } > − ([45, pp 338]) , yields | n X i =1 T ′ j ( x i ) b i ( x ) | = jτ n . (4.9)Thus, by (4.8), (4.9) and (1.5), we find | E ( T j , x ) | ≤ jτ n for j ≥ n , and then the errorof Hermite-Fej´er interpolation (4.7) satisfies | E ( f, x ) | = | f ( x ) − H n − ( f, x ) | ≤ ∞ X j =2 n | c j || E ( T j , x ) | = 2 τ n ∞ X j =2 n j | c j | , which, following [59], leads to the desired result.8 From the definition of τ n (4.6), we see that when α = β = − the convergence order on n is the lowest. In addition, from Szabados [44] (also see Sadiq and Viswanath [36]), we seethat the convergence of the higher order Hermite-Fej´er interpolation (2.15) at the Chebyshevpointsystem (1.3) satisfies k f − H mn − ( f ) k ∞ = (cid:26) O (log n ) k f − p ∗ k C m − [ − , , if m is odd O (1) k f − p ∗ k C m − [ − , , if m is even(4.10)where p ∗ is the best approximation polynomial of f with degree at most mn − k f − p ∗k C m − [ − , = max ≤ j ≤ m − k f ( j ) − ( p ∗ ) ( j ) k ∞ .Numerical examples also illustrate that the roots of (1 − x ) P ( , ) n − ( x ) are appropriate tohigher order Hermite-Fej´er interpolation. In the future work, we will consider the convergencerates on this pointsystem.It is worth noting that the new methods for Hermite barycentric weights at Gauss-Jacobipointsystems or Jacobi-Gauss-Lobatto pointsystems can be extended to Jacobi-Gauss-Radaupointsystems or the roots of other kinds of orthogonal polynomials, such as Laguerre polynomi-als, Hermite polynomials, etc., based on the works of [20], [57] and Chebfun [51]. REFERENCES[1] M. Abramowitz and I.A. Stegun , Handbook of Mathematical Functions , National Bureau of Standards,Washington, D.C., 1964.[2] S. Bernstein , Quelques remarques sur l’interpolation , Comm. Soc. Math. Charkow, 14 (1914).[3] S. Bernstein , Sur la limitation des valeurs d’un polynome etc. , Bull. de l’Acad. des Science de I’U.R.S.S.,(1931), 1025-1050.[4] S. Bernstein , Sur l’ordre de la meilleure approximation des fonctions continues par les polynˆomes de degr´edonn´e , Mem. Cl. Sci. Acad. Roy. Belg., 4 (1912), 1-103.[5] J. P. Berrut and L.N. Trefethen , Barycentric Lagrange interpolation , SIAM Rev., 46 (2004), 501-517.[6] R. Bojanic , A note on the precision of interpolation by Hermite-Fej´er polynomials , in Proceedings, Confer-ence on Constructive Theory of Functions, Budapest 1969 (G. Alexits et al., Eds.) pp. 69-76, AkademiaiKiado, Budapest, 1972.[7] I. Bogaert, B. Michiels and J. Fostier , O (1) Computation of Legendre Polynomials and Gauss-LegendreNodes and Weights for Parallel Computing , SIAM J. Sci. Comput., 34 (2012), C83-C101.[8] J. C. Butcher, R. M. Corless, L. Gonzalez-Vega, and A. Shakoori , Polynomial algebra for Birkhoffinterpolants , Numer. Alg., 56 (2011), 319-347.[9] G. J. Byrne, T. M. Mills and S.J. Smith , On Hermite-Fej´er type interpolation on the Chebyshev nodes ,Bull. Austral. Math. Soc., 47(1993), 13-24.[10] D. Calvetti and L. Reichel , On the evaluation of polynomial coefficients , Numer. Alg., 33(2003), 153C161.[11] E. W. Cheney , Introduction to Approximation Theory , McGraw-Hill, New York, 1966.[12] B. Della Vecchia, G. Mastroianni, P. V´ertesi , One-sided convergence conditions for Hermite-Fej´erinterpolation of higher order of Lagrange type , Results in Math. 34 (1988), 294-309.[13] G. Faber , ¨Uber die interpolatorische Darstellung stetiger Funktionen , Jahresber. Deut. Math. Verein. 23(1914), 192-210.[14] L. Fej´er , ¨Uber Interpolation , Nachrichten der Gesellschaft der Wissenschaften zu G¨ottingen Mathematisch-physikalische Klasse, 1916, 66-91.[15] L. Fej´er , Die Absch¨atzung eines Polynoms in einem Intervalle, wenn Schranken f¨ur seine Werte understen Ableitungswerte in einzelnen Punkten des Intervalles gegeben sind, und ihre Anwendung auf dieKonvergenzfrage Hermitescher Interpolationsreihen , Math. Zeitschrift, 32(1930), 425-457.[16] L. Fej´er , Lagrangesche interpolation und die zugeh¨origen konjugierten Punkte , Math. Ann., 106(1932),1-55.[17] L. Fej´er , Bestimmung derjenigen Abszissen eines Intervalles, f¨ur welche die Quadratsumme der Grund-funktionen der Lagrangeschen Interpolation im Intervalle ein M¨oglichst kleines Maximum Besitzt , An-nali della Scuola Norm sup. di Pisa, 1 (1932), 263-276.[18] B. Fischer and L. Reichel , Newton interpolation in Fej´er and Chebyshev points , Math. Comp., 53 (1989),265-278.[19] W. Gautschi , High-order Gauss-Lobatto formulae , Numer. Alg., 25 (2000), 213-222.[20] A. Glaser, X. Liu and V. Rokhlin , A fast algorithm for the calculation of the roots of special functions ,SIAM J. Sci. Comput., 29 (2007), 1420-1438.[21] S. J. Goodenough and T.M. Mills , A new estimate for the approximation of functions by Hermite-Fej´erinterpolation polynomials , J. Approx. Theory, 31(1981), 253-260.[22] S. J. Goodenough and T.M. Mills , On interpolation polynomials of the Hermite-Fej´er type II , Bull.Austral. Math. Soc., 23(1981), 283-291.[23] G. Gr¨unwald , ¨Uber Divergenzerscheinungen tier Lagrangeschen Interpolationspolynome , Acta Szeged, 7(1935), 207-211. [24] G. Gr¨unwald , On the theory of interpolation , Acta Math., 75(1942), 219-245.[25] N. Hale and A. Townsend , Fast and accurate computation of Gauss-Legendre and Gauss- Jacobi quadra-ture nodes and weights , SIAM J. Sci. Comput., 35(2013), A652-A674.[26] N. Hale and L. N. Trefethen , Chebfun and numerical quadrature , Science in China, 55 (2012), 1749-1760.[27] P. Henrici , Essentials of Numerical Analysis , Wiley, New York, 1982.[28] N. J. Higham , Accuracy and Stability of Numerical Algorithms , 2nd ed., SIAM, Philadelphian, 2002.[29] N. J. Higham , The numerical stability of barycentric Lagrange interpolation , IMA J. Numer. Anal., 24(2004), 547-556.[30] F. Locher , On Hermite-Fejer interpolation at Jacobi zeros , J. Approx. Theory, 44(1985), 154-166.[31] K. K. Mathur and R. B. Saxena , On the convergence of quasi-Hermite-Fej´er interpolation , Pacific J.Math., 20(1967), 197-392.[32] J. Marcinkiewicz , Sur la divergence des polynomes d’interpolation , Acta Szeged, 8 (1937), 131-135.[33] E. Moldovan , Observatii asupra unor procede de interpolare generalizate , Acad. Repub. Pop. Rom. Bul.Stinte Sect. Stinte Mat. Fiz. 6 (1954), 472-482.[34] P. Nevai and P. V´ertesi , HermiteCFej´er interpolation at zeros of generalized Jacobi polynomials , Approx-imation Theory IV. Academic Press, 1983, 629-630.[35] T. Popoviciua , Asupra demonstratiei teoremei lui Weierstrass cu ajutorul polynoamelor de interpolare ,Acad. Rep. &. Pop. Rom. Lucrarile sessiunii generala stintifice din 2-12 iunie 1950 (1951) 1664-1667.[36] B. Sadiq and D. Viswanath , Barycentric Hermite interpolation , SIAM J. Sci. Comput., 35(2013), 1254-1270.[37] H. E. Salzer , Lagrangian interpolation at the Chebyshev points x n,v = cos( vπ/n ) , v = O (1) n ; someunnoted advantages , Comput. J., 15 (1972), 156-159.[38] C. Schneider and W. Werner , Hermite interpolation: The barycentric approach , Computing, 46 (1991),35-51.[39] A. Sharma and J. Tzimbalario , Quasi-Hermite-Fej´er interpolation of higher order , J. Approx. Theory,13(1975), 431-442.[40] Y. Shi , A theorem of Gr¨unwald-type for Hermite-Fej´er interpolation of higher order , Constr. Approx.,10(1994), 439-450.[41] Y. Shi , On Hermite interpolation , J. Approx. Theory, 105(2000), 49-86.[42] O. Shisha, C. Sternin and M. Fekete , On the accuracy of approximation to given functions by certaininterpolartory polynomials of given degree , Riveon Lematematika 8 (1954), 59-64.[43] X. Sun , Approximation of continuous functions by Hermite-Fej´er type interpolation polynomials , J. Math.Research Expo., 3(1983), 45-50.[44] J. Szabados , On the order of magnitude of fundamental polynomials of Hermite interpolation , Acta Math.Hungar., 61 (1993), 357-368.[45] G. Szeg¨o , Orthogonal Polynomials , Colloquium Publications 23, A, Providence, Rhode Island, 1939.[46] L. Szili , Uniformly weighted convergence of Gr¨unwald interpolation process on the roots of Jacobi polyno-mials , Annales Univ. Sci. Budapest., Sect. Comp., 29(2008), 245-261.[47] H. Tal-Ezer , High degree polynomial interpolation in Newton form , SIAM J. Sci. Statisti. Comput., 12(1991), 648-667.[48] L. N. Trefethen , Spectral Methods in MATLAB , SIAM, Philadelphia, 2000.[49] L. N. Trefethen , Is Gauss quadrature better than Clenshaw-Curtis? , SIAM Rev., 50(2008), 67-87.[50] L. N. Trefethen , Approximation Theory and Approximation in Practice , SIAM, Philadelphia, 2012.[51] L. N. Trefethen and others , Chebfun Version 4.0 P. V´ertesi , Hermite-Fej´er type interpolations. III , Acta Math. Acad. Sci. Hung., 34 (1979), 67-84.[53] P. V´ertesi , ρ -normal point systems , Acta Math. Acad. Sci. Hung., 34 (1979), 267-277.[54] P. V´ertesi , Convergence criteria for Hermite-Fejer interpolation based on Jacobi abscissas , in Foundations,Series and Operators, Proc. of Int. Conf. in Budapest, 1980, v.II, pp.1253-1258, North Holland, 1983.[55] P. V´ertesi , Hermite-Fej´er interpolation of higher order , Acta Math. Hung., 54 (1989), 135-152.[56] H. Wang and S. Xiang , On the convergence rates of Legendre approximation , Math. Comp., 81 (2012),861-877.[57] H. Wang, D. Huybrechs and S. Vandewalle , Explicit barycentric weights for polynomial interpolationin the roots or extrema of classical orthogonal polynomials , arXiv: 1202.0154, 2013, Math. Comp., toappear.[58]