Chebyshev Interpolation

Journal of Online Mathematics and Its Applications

Volume 6. August 2006. Article ID 1297

Chebyshev Interpolation: An Interactive Tour

Scott A. Sarra

Contents

  1. Introduction
  2. Chebyshev Polynomials
  3. Continuous Chebyshev Expansion
  4. Discrete Chebyshev Expansion
    1. Interpolating Partial Sum
    2. Aliasing
  5. Rates of Convergence
  6. Filters
  7. Current Research Areas
  8. Further Explorations
  9. Summary
  10. References

1. Introduction

Most areas of numerical analysis, as well as many other areas of mathematics as a whole, make use of the Chebyshev polynomials. In several areas, e.g. polynomial approximation, numerical integration, and pseudospectral methods for partial differential equations, the Chebyshev polynomials take a significant role. In fact, the following quote has been attributed to a number of distinguished mathematicians:

"The Chebyshev polynomials are everywhere dense in numerical analysis."

In this article we use Java applets to interactively explore some of the classical results on approximation using Chebyshev polynomials. We also discuss an active research area that uses the Chebyshev polynomials. Mason and Handscomb (2003) and Rivlin (1974) are devoted to the Chebyshev polynomials and may be consulted for more detailed information than we provide in this brief presentation. The Chebyshev polynomials are named for Pafnuty Chebyshev. You can read a brief biography of Chebyshev at Wikipedia

The article uses four applets:

The CP applet and the CA applet are used frequently and thus open in separate windows that you can keep open as you read the text. The CA applet window also gives instructions for using the applet and definitions of the functions used in the applet.

2. Chebyshev Polynomials

The Chebyshev Polynomials (of the first kind) are defined by as

(1) Image: Eqn1.gif

They are orthogonal with respect to the weight Image: Eqn2.gif on the interval Image: Eqn3.gif. Intervals Image: Eqn4.gif other than Image: Eqn3.gif are easily handled by the change of variables Image: Eqn5.gif. Although not immediately evident from definition (1), Tn is a polynomial of degree n. From definition (1) we have that Image: Eqn7.gif and Image: Eqn8.gif.

Exercise. Use basic trig identities to establish the triple recursion relation

(2) Image: Eqn10.gif

Using equation (2) we see

Image: Eqn11.gif

and that the Chebyshev polynomial Tn is indeed a polynomial of degree n.

Applet Activity

What do the Chebyshev polynomials look like? The Chebyshev polynomials of degree n = 0, 1, ..., 12 can be plotted in the CP applet. Move the slider to change the degree. Notice that Image: Eqn12.gif. Since Tn is a degree n polynomial we can observe as expected that it has n zeros, which in this case are real and distinct and located in .Image: Eqn3.gif.

Exercise. Show that the zeros of Tn are

(3) Image: Eqn13.gif

The zeros are known as the Chebyshev-Gauss (CG) points.

3. Continuous Chebyshev Expansion

The infinite continuous Chebyshev series expansion is

(4) Image: Eqn14.gif

where

(5) Image: Eqn15.gif

The single prime notation in the summation indicates that the first term is halved. Truncating the series after N + 1 terms, we get the truncated continuous Chebyshev expansion:

(6) Image: Eqn16.gif

There are several functions in which the integral for the coefficients Image: Eqn17.gif can be evaluated explicitly, but this is not possible in general. Examples included in the CA applet for which a continuous truncated expansion can be derived are the sign function f1, the square root function f4, and the absolute value function f5 (open the applet window to review the definitions of these functions).

The conditions which must be placed on f to ensure the convergence of the series (4) depend on the type of convergence to be established: pointwise, uniform, or L2. At the lowest level, the series (4) converges pointwise to f at points where f is continuous in Image: Eqn3.gif and converges to the left and right limiting values of f at any of a finite number of jump discontinuities in the interior of the interval.

Applet Activity

The sign function in the CA applet has a jump discontinuity at x0 = 0 and has the limiting values on each side of the discontinuity of Image: Eqn20.gif and Image: Eqn21.gif. Thus the series converges to zero at this point, i.e.

Image: Eqn22.gif

for sufficiently large N. In the applet select the sign function from the Functions menu and check the blue continuous, S option on the Approximation menu. Using the slider at the bottom of the applet, slowly adjust N from N = 7 to N = 128 and observe that the value of Image: Eqn23.gif is approximately zero.

Exercise. Show that if f is an even function then Image: Eqn24.gif for k = 1, 3, 5, ... If f is an odd function then Image: Eqn24.gif for k = 0, 2, 4, ....

Applet Activity

The result in the last exercise can be observed in the truncated continuous expansion of Image: Eqn91.gif and Image: Eqn92.gif (even) and f1(x) = sign(x) (odd) in the CA applet. For example, select the even function f4 which is labeled as sqrt on the Functions menu and select the blue continuous, S option on the Approximation menu. Then on the Options menu check plot coefficients and using the slider slowly adjust N from N = 7 to N = 21. In the right window observe that Image: Eqn24.gif for k = 1, 3, 5, .... The magnitude of the coefficients can also be viewed with the y-axis scaled logarithmically (semiLogY on the Options menu). However, in this case the coefficients which are zero are not plotted as log(0) is undefined.

4. Discrete Chebyshev Expansions

When the integral in (5) can not be evaluated exactly, we can introduce a discrete grid and use a numerical quadrature (integration) formula. Several possible grids, and related quadrature formulas exist. The Chebyshev-Gauss-Lobatto (CGL) points

(7) Image: Eqn28.gif

are a popular choice of quadrature points. The CGL points are where the Image: Eqn29.gif extrema of Image: Eqn6.gif occur plus the endpoints of the interval Image: Eqn3.gif.

Applet Activity

Using the CP applet, observe how the extrema of the Chebyshev polynomials are not evenly distributed and how they cluster around the boundary. In the CA applet, the CGL points may be plotted by checking plot CGL points on the Options menu. Try this with the sign function starting with N = 9 and then with increasing N.

Exercise. Show that Image: Eqn30.gif at the Image: Eqn29.gif CGL points.

The corresponding CGL quadrature formula is

(8) Image: Eqn31.gif

The double prime notation in the summation indicates that the first and last terms are halved. If f is a polynomial of degree less than or equal to Image: Eqn32.gif, the CGL quadrature formula is exact. This is remarkable accuracy considering that the values of the integrand are only known at the N +1 CGL points. Using the CGL quadrature formula to evaluate the integral in (5), the discrete Chebyshev coefficients Image: Eqn33.gif are defined to be

(9) Image: Eqn34.gif

and the discrete truncated partial sum is

(10) Image: Eqn35.gif

Using definition (9) takes Image: Eqn39.gif floating point operations (flops) to evaluate the discrete Chebyshev coefficients. For large N, a better choice is the fast cosine transform (FCT) (Briggs and Henson, 1995) which takes Image: Eqn36.gif flops. For example, if N = 1000, Image: Eqn37.gif while Image: Eqn38.gif. The extreme efficiency of the FCT is one reason for the popularity of Chebyshev approximations in applications.

5.1 Interpolating Partial Sum

Requiring that the approximation be interpolating, i.e., requiring that it satisfy

(11) Image: Eqn40.gif

we get the interpolating partial sum

(12) Image: Eqn41.gif

The interpolating partial sum would be equal to the truncated series with the coefficients approximated via CGL quadrature except the last coefficient is halved. This is due to the choice of quadrature points. If Gaussian quadrature, which uses the Chebyshev-Gauss (CG) points, had been used instead of CGL quadrature, the interpolating and discrete truncated partial sum would be identical. The CG points are the zeros of Tn and do not include Image: Eqn42.gif. Chebyshev pseudospectral methods for solving PDEs usually incorporate the CGL points and not the CG points. The reason for this is that the discrete grid must include the boundary points so that the boundary conditions of the PDE can be incorporated into the numerical approximation.

Applet Activity

Using the CA applet, we can observe the difference between SN, PN, and IN. For example if we use the sign function (select sign from the Functions menu) with N = 11 (set N using the slider at the bottom of the applet) and plot the CGL points (check plot CGL points on the Options menu) we see that IN goes through the interpolation sites while SN and PN do not (On the Approximations menu, select the blue interpolation, I and then the red discrete, P. Then select the red continuous, S to make the next comparison).

Since (12) is a polynomial of at most degree N that satisfies the interpolation condition (11) at N + 1 distinct points, a standard result from numerical analysis tells us that IN is the unique interpolating polynomial (see Burden and Faires (1995), p. 106). The interpolating polynomial may be written in several equivalent forms: Lagrange, Newton, and Barycentric. For information on the merits of each form, see Berrut and Trefethen (2004). The Lagrange form of the interpolating polynomial is

Image: Eqn46.gif

where

(13) Image: Eqn47.gif

are cardinal polynomials that satisfy

(14) Image: Eqn48.gif

The Lagrange form gives an error term of the form

(15) Image: Eqn49.gif

where

(16) Image: Eqn50.gif

The underlying function f(x) is often unknown and the number Image: Eqn51.gif is only known in simple examples. Thus, Image: Eqn52.gif is the only part of the error term which can be controlled. By using the CG or CGL points as interpolation cites, Image: Eqn52.gif is made nearly as small as possible (see Burden and Faires (2005), p. 507). On the other hand, it is well known that polynomial interpolation in equally spaced points can be troublesome. The classic example provided by Runge is the function

(17) Image: Eqn53.gif

For the function (17), equidistant polynomial interpolation diverges for Image: Eqn54.gif. By using the CGL points (7), which cluster densely around the endpoints of the interval, as interpolation sites the nonuniform convergence (the Runge Phenomenon) associated with equally spaced polynomial interpolation is avoided.

Applet Activity

The RP applet below illustrates equidistant and Chebyshev interpolation for the Runge example (17). The applet starts with N = 15 and equidistant interpolation. Use the slider to increase N and observe that the oscillations near the boundary become larger and that the approximation is good for |x| < 3.63. Select the CGL button at the top of the applet and observe that the oscillations near the boundary disappear.

The RP Applet
The RP applet

5.2 Aliasing

Introducing a discrete grid leads to aliasing. The discrete coefficients can be expressed in terms of the continuous coefficients as

(18) Image: Eqn55.gif

As an example consider the sign function with N = 9. The difference between the discrete coefficient a5 and the continuous coefficient Image: Eqn56.gif can be quantified by the aliasing relation (18) as

Image: Eqn57.gif

This relation is a result of the fact that on the discrete grid, T5 is identical to T23, T41, T59, ... and also to T13; T31; T49, ... as is illustrated in Figure 1.

Figure 1. On the CGL grid (open black circles) for N = 9, T5 is identical to T13 (green)
and T23 (cyan). Points of intersection on the CGL grid are marked with red *'s.
convergence plot

The image was produced with the following Matlab script:

   N = 9;   M = 600;
   x = -cos(pi*(1:9)./N);      % extrema and endpoints of T10
  xp = linspace(-1,1,M);

   T23 = cos(23*acos(xp));     %  cyan
   T13 = cos(13*acos(xp));     %  green
    T5 = cos( 5*acos(xp));     %  blue
   T5g = cos( 5*acos(x));      %  T_5(x) (red *)

  XGL10 = zeros(1,length(x));  %  CGL pts (open black circles)

  plot(xp,T5,'b',xp,T13,'g',x,T5g,'r*',x,XGL10,'k-o',xp,T23,'c')
  xlabel 'x', ylabel 'T'

Applet Activity

In the CA applet, observe the difference between the odd numbered coefficients of the S9, P9 and I9 approximations of the sign function (select sign from the Functions menu and set N = 9 using the slider at the bottom of the applet). On the Approximations menu, select the blue interpolation, I and then select the red continuous, S. On the Approximations menu select plot coefficients. There is no difference in the even numbered coefficients, as the sign function is odd. Thus the continuous even coefficients that are involved in the aliasing relation are all zero. The difference in the odd coefficients is due to aliasing. Make similar comparison with the truncated discrete series by selecting the blue discrete, P from the approximations. Again there is a difference in the odd coefficients that is due to aliasing.

Now compare the two discrete approximations, I9 (blue interpolation, I) and P9 (red discrete, P). The coefficients are identical, but the approximations are different due to Image: Eqn58.gif being halved in the interpolating approximation but not in the truncated series.

5. Rates of Convergence

Repeatedly integrating equation (5) by parts we get

(19) Image: Eqn59.gif

Thus, if f is m-times (Image: Eqn60.gif) continuously differentiable in Image: Eqn3.gif the above integral will exist and we can conclude that

(20) Image: Eqn61.gif

If we make a careful choice of which definition of the integral to use, the same result can be shown to be true if f is (Image: Eqn62.gif)-times differentiable a.e. (almost everywhere) in Image: Eqn3.gif with its (Image: Eqn62.gif)'th derivative of bounded variation in Image: Eqn3.gif.

Since the absolute value of each Tk is bounded above by 1 on Image: Eqn3.gif, it follows that the truncation error for the continuous expansion is bounded by the sum of the absolute value of the neglected coefficients:

(21) Image: Eqn63.gif

A similar bound, with an additional factor of two, holds for the interpolating partial sum:

(22) Image: Eqn64.gif

From (20), (21), and (22) we conclude that

(23) Image: Eqn65.gif

and

(24) Image: Eqn66.gif

If f is infinitely differentiable the convergence is faster than Image: Eqn67.gif no matter how large we take m. This is commonly termed spectral accuracy or exponential accuracy. If f can be extended to an analytic function in a suitable region of the complex plane, the pointwise error on Image: Eqn3.gif can be shown to be

(25) Image: Eqn68.gif

for some r > 1 (Mason and Handscomb (2003)). In Figure 2 the rather slow decay rate of the error with increasing N is illustrated for the absolute value function f5 for which m = 1. This can be contrasted with the rapid spectral convergence of the infinitely smooth function f2. Notice that the decay of error for the smooth function ceases at about N = 140. This is due to the accuracy of the representation of floating point numbers on the computer which limits accuracy to about 14 or 15 decimal places.

Figure 2. Convergence of an infinitely differentiable function Image: Eqn90.gif
versus convergence of a continuous function Image: Eqn70.gif.
Convergence plot

No matter what rate of decay the coefficients have, the convergence rate is only observed for n > n0. Using an approximation with fewer than n0 terms may result in a very bad approximation. For example, the decay rate of the coefficients of the infinitely smooth function in the applet is not yet evident for N = 17 and the approximation is very poor.

Applet Activity

Equation (19) allows us to conclude that if f is a polynomial of degree N, then Image: Eqn71.gif for all n > N since Image: Eqn72.gif for n > N. In the CA applet select the 7th degree polynomial from the Functions menu. Use the slider at the bottom of the applet to set N to 9. From the Options menus check plot coefficients and semiLogY. Observe that Image: Eqn73.gif (to within machine precision) for n > 7.

If m = 0, i.e., f is discontinuous, the accuracy of the Chebyshev approximation methods reduces to O(1) near the discontinuity. Sufficiently far away from the discontinuity, the convergence will be slowed to Image: Eqn74.gif. Oscillations will be present near the discontinuity and they will not diminish as Image: Eqn75.gif. Additionally, the oscillations will not even be localized near a discontinuity. This situation is referred to as the Gibbs phenomenon. A nice history of the Gibbs phenomenon can be found in Hewitt and Hewiit (1979).

Applet Activity

In the CA applet, select the sign function from the Functions Menu. From the Options menus uncheck plot coefficients and check semiLogY. Use the slider at the bottom of the applet to slowly change N from 10 to 256. Observe that the maximum amplitude of the overshoot at the discontinuity does not decrease with increasing N. Observe that sufficiently far away from the discontinuity that the oscillations are slowly decaying. Now check plot coefficients on the Options menu and again use the slider at the bottom of the applet to slowly change N from 10 to 256. Notice that the coefficients are decaying, but at a very slow rate. Spectral convergence has been lost due to the discontinuity. Select the smooth function from the Functions menu and compare how fast the coefficients of this function decay compared to the sign function.

6. Filters

Spectral filters can be used to enhance the decay rate of the Chebyshev coefficients (Vandeven (1992)) and to lessen the effects of the Gibbs phenomenon. The filtered Chebyshev approximation is

(26) Image: Eqn76.gif

where Image: Eqn77.gif is a spectral filter. A pth (p > 1) order spectral filter is defined as a sufficiently smooth function satisfying

(27) Image: Eqn78.gif
(28) Image: Eqn79.gif
(29) Image: Eqn80.gif

The convergence rate of the filtered approximation is determined solely by the order of the filter and the regularity of the function away from the point of discontinuity.

If p is chosen increasing with N, the filtered expansion recovers exponential accuracy away from a discontinuity. Assuming that f has a discontinuity at x0 and setting Image: Eqn81.gif, the estimate

(30) Image: Eqn82.gif

holds where K is a constant. If p is sufficiently large, and d(x) not too small, the error goes to zero faster than any finite power of N, i.e. spectral accuracy is recovered. When x is close to a discontinuity the error increases. If d(x) = O(1/N) then the error estimate is O(1).

Many different filter functions are available, but perhaps the most versatile and widely used filter is the exponential filter

(31) Image: Eqn83.gif

of order 2p. In order for condition (29) to be satisfied, the parameter Image: Eqn84.gif is taken as Image: Eqn85.gif where Image: Eqn86.gif is defined as machine zero. On a 32-bit machine using double precision floating point operations, Image: Eqn87.gif and Image: Eqn88.gif.

Applet Activity

The exponential filter is implemented in the CA applet. The default order of the filter is 4 (p = 2). Select the sign function from the Functions menu. From the Approximations menu select the blue interpolation and red filter options. From the Options menu check semiLogY and uncheck connect. Use the slider to increase N and observe the rapid decrease in the error of the filtered approximation away from the discontinuity. The filter has restored spectral accuracy at points sufficiently far away from the discontinuity. Next, check plot coefficients on the Options menu and compare the filtered and unfiltered coefficients. Now, display the parameters dialog from the Options menu and enter 1 in the filter order box to change the order of the filter to 2. Repeat the above experiments. Observe how the sharp front at the discontinuity is rounded or smeared in the filtered approximation by the low order filter. Enter 4 in the filter order box to change the order of the filter to 8 and repeat. What do you observe?

Applet Activity

In the CA applet, select the absolute value function from the Functions menu and repeat the previous applet activity.

Applet Activity

The EF applet illustrates the strength of the damping applied in equation (26) to the coefficients ak from k = 0, 1, ..., N for filters of order 2 to 32. The slider at the bottom of the applet can be used to change the order of the filter. Observe that as the order of the filter increases that less damping is applied to the coefficients with small k.

The EF Applet

7. Current Research Area

Chebyshev approximation is an old and rich subject. However, many areas that employ Chebyshev polynomials have open questions that have attracted the attention of current researchers. One example is pseudospectral methods for the numerical solution of partial differential equations (PDEs). Chebyshev pseudospectral methods, which are based on the interpolating Chebyshev approximation (12), are well established as powerful methods for the numerical solution of PDEs with sufficiently smooth solutions. Interpolation means that f, the function that is approximated, is a known function. The terms collocation and pseudospectral are applied to global polynomial interpolatory methods for solving differential equations for an unknown function f. Detailed information on pseudospectral methods may be found in the standard references: Boyd (2000), Canuto, et al. (1988), Funaro (1992), Gottlieb, et al. (1984), Gottlieb and Orszag (1977), and Trefethen (2000).

Many important PDEs have discontinuous (or nearly discontinuous) solutions. See the article Sarra (2003) for a discussion of one such class of PDEs, nonlinear hyperbolic conservation laws. In these cases, the Chebyshev pseudospectral method produces approximations that are contaminated with Gibbs oscillations and suffer from the corresponding loss of spectral accuracy, just like the Chebyshev interpolation methods that the pseudospectral methods are based on.

An active research area is the development of postprocessing methods to remove the Gibbs oscillations from PDE solutions and to restore spectral accuracy. Spectral filters may be used but they perform poorly in the neighborhood of discontinuities. More sophisticated methods that do better in the area of discontinuities, but they may need to know the exact location of the discontinuities. The methods include Spectral Mollification, Gegenbauer Reconstruction Gottlieb (1997), Padé Filtering, and Digital Total Variation Filtering.

Several postprocessing methods with applications are discussed in Sarra (2003) with supporting web material at the Matlab Postprocessing Toolbox. The ultimate goal is a "black box" postprocessing algorithm, which can be given an oscillatory PDE solution and return a postprocessed solution with spectral accuracy restored.

8. Further Explorations

In addition to the exponential filter, other postprocessing methods for lessening the effects of the Gibbs phenomenon exist. Explore some of them which include:

  1. Other spectral filters. See Boyd (1996).
  2. Reprojection methods. Reprojection methods work by projecting the slowly converging Chebyshev approximation onto a Gibbs complementary basis in which the convergence is faster. See Gelb (2005), Gottlieb (1997), and Sarra (2003).
  3. Padé based reconstruction. Padé methods reconstruct the Chebyshev polynomial approximation as a rational approximation (Mace, 2003).
  4. Digital Total Variation (DTV) filtering. DTV methods which were developed in image processing have been used to postprocess Chebyshev approximations. See Sarra (2006).

9. Summary

Chebyshev approximation and its relation to polynomial interpolation at equidistant nodes has been discussed. We have illustrated how the Chebyshev methods approximate with spectral accuracy for sufficiently smooth functions and how less smoothness slows down convergence. We have illustrated how the presence of a discontinuity leads to lack of convergence at the discontinuity and leads to slowed convergence away from the discontinuity. We have described the Gibbs phenomenon which is characterized by a lack of or slow convergence as well as non-physical oscillations. Spectral filtering was discussed as a method used to lessen the effects of the Gibbs phenomenon and to restore spectral accuracy sufficiently far away from a discontinuity. Postprocessing methods to lessen the effects of the Gibbs oscillations are an active research area which would be an excellent topic for undergraduate research or as the topic of a Masters thesis.

10. References

  1. J. Berrut and L. N. Trefethen. Barycentric Lagrange interpolation. SIAM Review, 46(3): 501-517, 2004.
  2. John P. Boyd. The Erfc-Log Filter and the asymptotics of the Vandeven and Euler sequence accelerations. Houston Journal of Mathematics, 267--275, 1996.
  3. John P. Boyd. Chebyshev and Fourier Spectral Methods. Dover Publications, Inc, New York, second edition, 2000.
  4. W. Briggs and V. Henson. The DFT: An Owner's Manual for the Discrete Fourier Transform. SIAM, 1995.
  5. R. Burden and J. Faires. Numerical Analysis. Brooks Cole, eighth edition, 2005.
  6. Claudio Canuto, M. Y. Hussaini, Alfio Quarteroni, and Thomas A. Zang. Spectral Methods for Fluid Dynamics. Springer-Verlag, New York, 1988.
  7. D. Funaro. Polynomial Approximation of Differential Equations. Springer-Verlag, New York, 1992.
  8. A. Gelb and J. Tanner. Robust reprojection methods for the resolution of the Gibbs phenomenon. To appear in Applied and Computational Harmonic Analysis, 2005.
  9. David Gottlieb, M. Y. Hussaini, and Steven A. Orszag. Theory and application of spectral methods. In R. G. Voigt, D. Gottlieb, and M. Y. Hussaini, editors, Spectral Methods for Partial Differential Equations, 1--54. SIAM, Philadelphia, 1984.
  10. David Gottlieb and Steven A. Orszag. Numerical Analysis of Spectral Methods. SIAM, Philadelphia, PA, 1977.
  11. David Gottlieb and Chi-Wang Shu. On the Gibbs phenomenon and its resolution. SIAM Review, 39(4): 644--668, 1997.
  12. E. Hewitt and R.E. Hewitt. The Gibbs-Wilbraham phenomenon: an episode in Fourier analysis. History of Exact Science, 21: 129--160, 1979.
  13. R. Mace. Reduction of the Gibbs Phenomenon in Chebyshev approximations via Chebyshev-Padé filtering. Master's thesis, Marshall University, 2005.
  14. J. Mason and D. Handscomb. Chebyshev Polynomials. CRC, 2003.
  15. T. Rivlin. The Chebyshev Polynomials. Wiley, 1974.
  16. S. A. Sarra. The method of characteristics with applications to conservation laws. Journal of Online Mathematics and its Applications, 3, 2003. (accessed December 16, 2005).
  17. S. A. Sarra. The spectral signal processing suite. ACM Transactions on Mathematical Software, 29(2): 1--23, 2003.
  18. S. A. Sarra. Digital total variation filtering as postprocessing for Chebyshev pseudospectral methods for conservation laws, Numerical Algorithms, 41:17-33, 2006.
  19. L. N. Trefethen. Spectral Methods in Matlab. SIAM, Philadelphia, 2000. 9
  20. H. Vandeven. Family of spectral filters for discontinuous problems. SIAM Journal of Scientific Computing, 6:159--192, 1991.