On the Lebesgue constant of Berrut's rational ... - Semantic Scholar

Report 2 Downloads 83 Views
On the Lebesgue constant of Berrut’s rational interpolant at equidistant nodes Len Bos∗, Stefano De Marchi†and Kai Hormann‡ October 25, 2010

Abstract We study the Lebesgue constant of the rational interpolant of Berrut (cf. [1]) when the interpolation points are equally distributed. In the more general case of the rational interpolant of Floater and Hormann (cf. [6]), we show by several numerical results, that the behavior of the Lebesgue constant on equally distributed points is consistent with that of Berrut’s interpolant.

1

Introduction

Suppose we want to approximate a function f : [a, b] → R by some function g that interpolates f at the n + 1 distinct interpolation points a = x0 < x1 < · · · < xn = b. Given a set of basis functions bi which satisfy the Lagrange property ( 1, if i = j, bi (xj ) = δij = 0, if i 6= j, the interpolant g can be written as g(x) =

n X

bi (x)f (xi ).

i=0

The Lebesgue constant of this interpolation operator is Λ = max Λ(x), a≤x≤b

where Λ(x) is the associated Lebesgue function Λ(x) =

n X

|bi (x)| .

i=0 ∗

University of Verona (Italy), [email protected] University of Padua (Italy), [email protected] ‡ University of Lugano (Switzerland), [email protected]

1

(1.1)

While the Lebesgue constant has been studied intensively in the case of polynomial interpolation, that is when bi are the Lagrange basis polynomials (cf. e.g. [3, 4, 7] and references therein), our aim is to investigate the rational interpolant that was introduced by Berrut (cf. [1]) with basis functions (−1)i bi (x) = x − xi

X n j=0

(−1)j , x − xj

i = 0, . . . , n

and to derive an upper bound for the associate Lebesgue function  X n n n j X X (−1) 1 Λ(x) = |bi (x)| = . |x − xi | x − xj j=0 i=0 i=0

(1.2)

(1.3)

It is worth mentioning, that a lower bound for linear rational interpolation, was already studied by Berrut and Mittelmann in [2], and more recently Carnicer, in [5], has proposed a weighted Lagrange interpolation on equidistant or almost equidistant points with small Lebesgue constants. But in the stability analysis of interpolation processes one is more interested in finding an upper bound for the Lebesgue constant (1.1). Indeed, the knowledge of such a bound will give information on the stability of the process and, moreover, in the conditioning of the problem. In section 2, we prove the main results of the paper that is, Berrut’s rational interpolant on equally distributed points has a Lebesgue constant bounded above by 2 + ln(n) and below by π2 ln(n + 2). Floater and Hormann in [6], provided a general family of barycentric rational interpolants with arbitrarily high approximation orders which includes Berrut’s interpolant as a special case. In section 3 we outline the study of the Lebesgue constant growth on equispaced points in the general setting of Floater and Hormann’s rational interpolant. Supported by several numerical experiments, in section 4, we will show that also in this case the growth is most likely logarithmic in n.

2

Main result

In this section we prove our main results, that is that the Lebesgue constant of Berrut’s interpolant on equidistant points grows as ln n, n being the number of points. For simplicity (and without loss of generalization) we assume that the interpolation interval is [0, 1]. Let us also assume that the interpolation points are equally spaced with distance 1 h = , that is, n i xi = ih = , i = 0, . . . , n, n Our first result concerns the upper bound of the Lebesgue constant Theorem 2.1. Let the basis function bi (x) be given as in (1.2), the points xi as above, then Λ ≤ 2 + ln(n).

2

Proof. If x = xk for any k, then it follows from the interpolation property of the basis functions that Λ(x) = 1. So let xk < x < xk+1 for some k and consider the function n X

1 |x − xj | j=0 . Λk (x) = n X j (−1) (x − xk )(xk+1 − x) x − xj (x − xk )(xk+1 − x)

(2.4)

j=0

We first focus on the numerator in (2.4), N = (x − xk )(xk+1 − x)

n X j=0

= (x − xk )(xk+1 − x)

1 |x − xj |

k−1 X j=0

n X

1 1 1 1 + + + x − xj x − xk xk+1 − x j=k+2 xj − x k−1 X

!

n X 1 1 = (xk+1 − x) + (x − xk ) + (x − xk )(xk+1 − x) + x − xj j=k+2 xj − x j=0 ! n k−1 X X 1 1 + . = (xk+1 − xk ) + (x − xk )(xk+1 − x) x − xj j=k+2 xj − x j=0

As the xi are equally spaced with distance h = 1/n we have 1 n 1 = = xi − xj h(i − j) i−j for any i 6= j, and  2 1 h = 2 (x − xk )(xk+1 − x) ≤ 2 4n for xk < x < xk+1 . Therefore, ! n k−1 X X 1 1 1 1 + N≤ + 2 n 4n x − xj j=k+2 xj − xk+1 j=0 k ! k−1 n X X 1 n 1 n + = + 2 n 4n k − j j−k−1 j=0 j=k+2   1 1 1 1 1 1 1 1 = + + + ··· + + + + ··· + n 4n k k − 1 1 1 2 n−k−1  1 1 ≤ + ln(2k + 1) + ln(2n − 2k − 1) n 4n  1 1 = + ln (2k + 1)(2n − (2k + 1)) n 4n  1 1 ≤ + ln (2n/2)2 n 4n 1 1 = + ln(n). n 2n 3

!

We now turn to the denominator in (2.4), ignoring the absolute value and assuming both k and n to be even for the moment, so that D = (x − xk )(xk+1 − x)

n X (−1)j x − xj j=0

n k−1 X X 1 1 (−1)j (−1)j + + − = (x − xk )(xk+1 − x) x − xj x − xk xk+1 − x j=k+2 xj − x j=0 ! k−1 n j j X X (−1) (−1) = h + (x − xk )(xk+1 − x) − . x − xj j=k+2 xj − x j=0

!

Pairing the positive and negative terms in the rightmost factor adequately then gives k−1 n X X (−1)j (−1)j S= − x − xj j=k+2 xj − x j=0     1 1 1 1 1 1 + − + ··· + − − = x − x0 x − x2 x − x1 x − xk−2 x − xk−3 x − xk−1     1 1 1 1 1 + − − − + ··· + xk+2 − x xk+3 − x xk+4 − x xn−1 − x xn − x

(2.5)

Since both the leading term and all paired terms are positive, we have S>−

1 1 1 1 ≥− − = −2n − x − xk−1 xk+2 − x xk − xk−1 xk+2 − xk+1

and further D = h + (x − xk )(xk+1 − x)S ≥

1 1 1 1 1 + 2 (−2n) = − = . n 4n n 2n 2n

This bound also holds if n is odd as this only adds a single positive term 1/(xn − x) to S in (2.5), and if k is odd then a similar reasoning shows that D ≤ −1/(2n). Therefore, we have |D| ≥ 1/(2n) regardless of the parity of k and n, and combining the bounds for numerator and denominator in (2.4) yields   1 + 1 ln(n) = 2 + ln(n). Λ = max max Λk (x) ≤ n 2n1 k=0,...,n

xk <x<xk+1

2n

The next theorem concerns the lower bound of the Lebesgue constant for Berrut’s rational interpolant. The result is indeed an improvement of the corresponding result given by Berrut and Mittelmann in [2, Th. 3.1], that is Λ ≥ 1/(2n2 ). Theorem 2.2. Let the basis function bi (x) be given as in (1.2) and the points xi equally spaced in [0, 1], then Λ ≥ cn ln(n + 2). for constants cn > 0 such that limn→∞ cn = 2/π. 4

Proof. Let us consider the case of n = 2s + 1 and Pn 1

N j=0 |x−xj | Λ(x) = P . := n (−1)j D j=0 x−xj

Take x = 1/2. Then,

D

=

=

=

= |{z}

j 0 =j+s+1

= |{z}

j 0 =s−j

=

=

2s+1 2s+1 (−1)j X X (−1)j = 2 (2s + 1) j 1 (2s + 1) − 2j − 2 2s+1 j=0 j=0 s 2s+1 X j j X (−1) (−1) + 2 (2s + 1) (2s + 1) − 2j j=s+1 (2s + 1) − 2j j=0 s s X X (−1)j+s+1 (−1)j + 2 (2s + 1) (2s + 1) − 2j j=0 (2s + 1) − 2(j + s + 1) j=0 s s X j j X (−1) (−1) 2 + (−1)s+1 (2s + 1) (2s + 1) − 2j −2j − 1 j=0 j=0 s s X j X (−1)s−j (−1) 2 + (−1)s (2s + 1) (2s + 1) − 2(s − j) 2j + 1 j=0 j=0 s s j j X X (−1) (−1) 2 (−1)s + (−1)s (2s + 1) 2j + 1 2j + 1 j=0 j=0 s X (−1)j 4 (2s + 1) → π(2s + 1) . 2j + 1 j=0

Since it is an alternating series s X 4 (−1)j 4 π− ≤4 ≤π+ 2s + 3 2j + 1 2s + 3 j=0

that is 

4 π− 2s + 3



 (2s + 1) ≤ D ≤ π +

5

4 2s + 3

 (2s + 1)

(2.6)

For the numerator N at x = 1/2 we get N =

=

= =

2s+1 X

2s+1

X 1 1 1 = 2(2s + 1) − j |2s + 1 − 2j| 2s+1 j=0 j=0 2 ( s ) 2s+1 X X 1 1 2(2s + 1) + |2s + 1 − 2j| |2s + 1 − 2j| j=0 j=s+1 ) ( s 2s+1 X X 1 1 + 2(2s + 1) 2s + 1 − 2j j=s+1 2j − (2s + 1) j=0 Z s+1 s X 1 1 4(2s + 1) ≥ 4(2s + 1) dx 2j + 1 2x + 1 0 j=0

= 2(2s + 1) ln(2s + 3) . Hence   1 max Λ(x) ≥ Λ x∈[0,1] 2 2(2s + 1) ln(2s + 3)  ≥ 4 π + 2s+3 (2s + 1) 2 ln(2s + 3) 2 = ≈ ln(n + 2) . 4 π π + 2s+3 The same reasoning can be applied for the even case, i.e. n = 2s, which gives the same asymptotic limit. A comparison plot of the computed Lebesgue constant, the upper and lower bounds and the ln(n) is visible in Figure 1.

3

The case of Floater-Hormann rational interpolant

The rational interpolant introduced by Floater and Hormann, which generalizes Berrut’s rational interpolant, has basis functions (−1)i βi bi (x) = x − xi

X n j=0

(−1)j βj , x − xj

for certain positive weights β0 , . . . , βn given by Pj  d   k=0 k , βj = 2d ,   βn−j ,

i = 0, . . . , n

(3.7)

d-th row of the Bernoulli triangle if j ≤ d, if d ≤ j ≤ n − d, if j ≥ n − d.

(3.8)

This rational interpolant is based on blending local interpolating polynomials of degree d, with 2d ≤ n, hence its approximation order is O(hd+1 ). 6

Figure 1: Lebesgue constant for n ≤ 200 on uniform points for Berrut’s interpolant and the upper and lower bounds The study of the Lebesgue constant growth of such interpolant is more intriguing. Indeed, while the bound for the numerator trivially extends to   1 1 d N ≤2 + ln(n) n 2n by simply bounding all βk by 2d and pulling out this common factor, the denominator is more difficult to handle. All numerical experiments show that |D| ≥

1 ∀d>0 n

(3.9)

leading to the bound Λ ≤ 2d + 2d−1 ln(n)

(3.10)

and in particular, for d = 1, Λ ≤ 2 + ln(n) i.e. the same as for d = 0 which shows that the Lebesgue constants for d = 0 and d = 1 are almost identical (as shown below in Fig. 6). Hence, all we need to verify is that the inequality (3.9) holds. Once we shall be able to prove the above inequality, we could conclude that the following claim is true Conjecture 3.1. The Lebesgue function associated to the basis (3.7) is bounded from above as max Λ(x) ≤ 2d−1 (2 + ln(n)) . x∈[a,b]

At the present we are not able to prove this Conjecture. Nevertheless, in the next section we shall give numerical evidence of the bound above, by several examples on different sets of interpolation points.

7

4

Numerical experiments

4.1

Equally spaced points

All experiments have been done in the interval [−1, 1]. The plots of the Lebesgue function for d = 0, suggest that the maximum is attained near   1/n, if n = 4s,  2/n, if n = 4s + 1, ∗ x =  3/n, if n = 4s + 2,    0, if n = 4s + 3, for some k ∈ N, so that Λ = Λ(x∗ ). Notice that due to the symmetry of the Lebesgue function, the points x∗ could be taken in the left part of the interval. Numerically this can also be seen by taking more and more target points (as discretization of the interval [−1, 1]) because, as one can expect, finer is the grid the points where the maxima are taken, say x˜∗ , tend to the points x∗ (the real maxima). Consider the case n = 4s + 3. We have  X  X  X n 2s 2s 2s 2s n n i i i X X X X (−1) 1 (−1) 1 (−1) 1 Λ(x∗ ) = |bi (0)| = = = |x | x |x | x n − 2i n − 2i i i i i i=0 i=0 i=0 i=0 i=0 i=0 i=0 ! !  X  X 2s 2s 2s 2s X X 1 1 (−1)i 1 1 (−1)i 1 1 = − − − − = , n − 2i 2n n − 2i 2n 1 + 2i 2n 1 + 2i 2n i=0 i=0 i=0 i=0 (4.11) and since m X 1 ln(m) ∼ as m → ∞ 2i + 1 2 i=0 and

∞ X π (−1)i = , 2i + 1 4 i=0

we get the asymptotic estimate Λ∼

2 ln(n) as n → ∞. π

The same is true for the other three cases. In fact, Equation (4.11) can be generalized to !  2s ! 2s i X X 1 a (−1) b n n Λ(x∗ ) = , − − 1 + 2i 4 1 + 2i 4 i=0 i=0 where  1 +  n−1   2, an = n 1  −    n−1 1 − n−2

3 , n+1

if if 1 , if n+1 2 1 − n+2 , if n

n = 4s, n = 4s + 1, n = 4s + 2, n = 4s + 3,

and

8

 1 3 − n−1 + n+1 ,    2, bn = n 1 1  + n+1 ,  n−1   1 2 1 + n − n+2 , n−2

if if if if

n = 4s, n = 4s + 1, n = 4s + 2, n = 4s + 3.

Figure 2: Lebesgue function for 8 uniform points for Berrut’s interpolant

Figure 3: Lebesgue function for 8 uniform points for Floater-Hormann’s interpolant with d = 1. To support these considerations, we provide Figures 2, 3, 4. The case d = 0, i.e. Berrut’s interpolant, shows a Lebesgue constant with the same behavior of the corresponding constant at the Chebyshev points for polynomial interpolation. Indeed, by looking at Figure 5, the constant growth is close to π2 ln(n + 1) + 0.89 which is similar to the corresponding upper bound π2 ln(n + 1) + 1 for Chebyshev points (cf. [4]). For d > 1, the maximum seems to be taken closer to the end of the interval, but this requires further numerical investigations. However, even then the growth seems to be logarithmic as shown in Figure 6. In conclusion, for d fixed and equidistant points, the bound is of the form α + β ln(n), i.e. it 9

Figure 4: Lebesgue function for 8 uniform points for Floater-Hormann’s interpolant with d = 3.

Figure 5: Lebesgue constant growth up to n = 50 for d = 0, compared with the logarithm and the optimal growth at Chebyshev points for polynomial interpolation. grows logarithmically with n, in contrast to the exponential growth that happens for d = n (polynomial interpolation).

10

Figure 6: Lebesgue constant for uniform points.

4.2

Other point distributions

Let us now turn to another important case of interpolation points, namely the Chebyshev points iπ xi = − cos , i = 0, . . . , n. n Figure 7 shows the numerically computed Lebesgue constants for small degree (d = 0, 1, 2, 3) and 10 ≤ n ≤ 200. If we use the same weights wi as in the uniform setting, then the Lebesgue constant behaves very nicely (upper left), but if we use instead the weights that were proposed by Floater and Hormann (the ones that guarantee approximation order d + 1), then the Lebesgue constant grows logarithmically only for d = 0, 1 (other three pictures). In another experiment we used the logarithmically distributed interpolation points xi = ln (1/e + i/n(e − 1/e)) ,

i = 0, . . . , n,

and got a slightly more favorable behavior; see Figure 8. Here, the Lebesgue constant seems to grow logarithmically both for uniform weights (upper left) and for the ones that guarantee the optimal approximation order (other three figure), but still the uniform weights give much better bounds. Finally, we used randomly distributed points, but here the Lebesgue constant behaves badly for both choices of weights. Acknowledgments. This work has been done with support of the 60% funds, year 2010 of the University of Padua.

11

Figure 7: Lebesgue constant for Chebyshev points.

References [1] J.-P. Berrut, Linear Rational Interpolation of Continuous Functions over an Interval, in Mathematics of Computation 1943-1993: A Half-Century of Computational Mathematics, Proceedings of Symposia in Applied Mathematics, AMS, Providence, RI, 1994 (Ed. W. Gautschi), 261–264. [2] J.-P. Berrut and H. D. Mittelmann, Lebesgue Constant Minimizing Linear Rational Interpolation of Continuous Functions over the Interval, Computers Math. Appl., 33(6) (1997), 77–86. [3] L. Brutman, On the Lebesgue functions for polynomial interpolation, SIAM J. Numer. Anal., 15(4) (1978), 694–704. 12

Figure 8: Lebesgue constant for logarithmically distributed points. [4] L. Brutman, Lebesgue functions for polynomial interpolation-a survey, Annals Numer. Math., 4 (1997), 111–127. [5] J. M. Carnicer, Weighted interpolation for equidistant points, Numer. Algorithms 55(23) (2010), 223–232. [6] Michael S. Floater and Kai Hormann, Barycentric Rational Interpolation with no Poles and High Rates of Approximation, Numer. Math., 107(2) (2007), 315–331. [7] Simon J. Smith, Lebesgue constants in polynomial interpolation, Annales Math. Inf., 33 (2006), 109–123

13

[8] Q. Wang, P. Moin and G. Iaccarino, A rational interpolation scheme with superpolynomial rate of convergence, Annual Research Brief 2008, Centre for Turbulence Research, 31–54.

14