Chebyshev Interpolation: an interactive tour

Scott A. Sarra
Marshall University

* NOTE: For the proper typesetting of the mathematical symbols in this document, it must be viewed with Internet Explorer. In order for the applets to function, the Java 1.5 (5.0) runtime environment must be installed on your computer.

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 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 manuscript we make use of Java applets to interactively explore some of the classical results on approximation using Chebyshev polynomials. The three applets used are the Chebyshev approximation (CA) applet, the Chebyshev polynomial (CP) applet, and the exponential filter (EF) applet. We also discuss an active research areas that uses the Chebyshev polynomials. The following books [7,8] devoted to the Chebyshev polynomials may be consulted for more detailed information than we provide in this brief presentation.

2  Chebyshev Polynomials

The Chebyshev Polynomials (of the first kind) are defined as
Tn(x) = cos[ n arccos(x) ].
(1)
They are orthogonal with respect to the weight w(x) = (1-x2)-1/2 on the interval [-1,1]. Intervals [a,b] other than [-1,1] are easily handled by the change of variable x® [1/2][ (b-a)x + a + b ]. Although not immediately evident from definition (1), each Tn(x) is a polynomial of degree n. From definition (1) we have that T0( x ) = cos(0 ) = 1 and T1( x ) = cos(arccosx) = x. For n ³ 1 we have the triple recursion relation
Tn+1(x) = 2 x Tn(x) - Tn-1(x)
(2)
which can be shown to be true using basic trig identities. Using equation (2) we see
T2( x )
=
2x (x) - 1 = 2x2 - 1
T3( x )
=
2x (2x2 - 1) - x = 4x3 - 3x
T4( x )
=
2x (4x3 - 3x) - (2x2 - 1) = 8x4 - 8x2 + 1
:
and that the Chebyshev polynomials are indeed polynomials of degree n. What do the Chebyshev polynomials look like? The Chebyshev polynomials of degree k=0,1,¼,12 can be plotted in the CP applet below. Notice that |Tn(x)| £ 1.

3  Continuous Chebyshev Expansion

The infinite continuous Chebyshev series expansion is
f(x) » ¥
å
n=0 
 ¢an Tn(x)
(3)
where
an = 2
p
 ó
õ 		1

-1 
 (1-x2)-1/2  f(x)  Tn(x) dx.
(4)
Truncating the series after N+1 terms, we get the truncated continuous Chebyshev expansion
SN (x)= N
å
n=0 
 ¢an Tn(x).
(5)
The single prime notation in the summation indicates that the first term is halved. There are several functions in which the integral for the coefficients an 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 (26), f(x)=Ö{1-x2}, and f(x) = |x|.
The conditions which must be placed on f to ensure the convergence of the series (3) depend on the type of convergence to be established: pointwise, uniform, L2, or Cesaro mean. At the lowest level, the series (3) converges pointwise to f at points where f is continuous in [-1,1] and to the left and right limiting values of f at any of a finite number of jump discontinuities in the interior of the interval. For example, 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 f(x0+) = 1 and f(x0-) = -1. Thus the series converges to zero at this point, i.e.
SN ( x0 ) » 1
2
 [ f(x0+) +f(x0-) ].
If f(x) is an even function then ak = 0 for k = 1,3,5,.... If f(x) is an odd function then ak = 0 for k = 0,2,4, .... This can be observed in the truncated continuous expansion of f(x) = Ö{1-x2} and f(x)=|x| (even) and f(x) = sign(x) (odd) in the CA applet.

4  Discrete Chebyshev Expansions

When the integral in (4) 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
xk = -cos æ
è 		k p
N
 ö
ø 		    k=0,1,¼,N
(6)
are a popular choice. The CGL points are the n-1 extrema of Tn(x) plus the endpoints of the interval [-1,1]. 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 increasing N.
The corresponding CGL quadrature formula is
ó
õ 		1

-1 
 f(x) dx

Ö
1-x2
 » p
N
 N
å
j=0 
 ¢¢ f(xk).
(7)
Using the CGL quadrature formula, the discrete Chebyshev coefficients an are defined to be
an » an = 2
N
 N
å
k=0 
 ¢¢ f(xk) Tn(xk)
(8)
and the discrete truncated partial sum is
PN (x)= N
å
n=0 
 ¢an Tn(x).
(9)
Using definition (8) takes O(N2) operations to evaluate the discrete Chebyshev coefficients. For large N, a better choice is the fast cosine transform which takes O(N logN) operations. If f(x) is a polynomial of at most degree 2N-1, the CGL quadrature formula is exact, and PN (x) and SN (x) coincide.

4.1  Interpolating Partial Sum

Requiring that the approximation be interpolating, i.e., requiring that it satisfy
PN (xi) = f(xi)     i = 0,1,¼,N
(10)
we get the interpolating partial sum
IN (x)= N
å
n=0 
 ¢¢ an Tn(x).
(11)
The double prime notation in the summation indicates that the first and last terms are halved. 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(x) and do not include x ±1. 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. Using the applet we can observe the difference between SN, PN, and IN. For example if we use the sign function with N=11 and plot the CGL points we see that IN goes through the interpolation sites while SN and PN do not.
By using the CGL points (6), which cluster densely around the endpoints of [-1,1], as interpolation sites, the nonuniform convergence (the Runge Phenomena) associated with equally spaced polynomial interpolation is avoided.

4.2  Aliasing

Introducing a discrete grid leads to aliasing. The discrete coefficients can be expressed in terms of the continuous coefficients as
an = an + ¥
å
j=1 
 ( an+2jN + a-n+2jN ).
(12)
As an example consider the sign function with N=9. The difference between the discrete coefficient a5 and the continuous coefficient a5 can be quantified by the aliasing relation (12) as
a5 - a5
=
a23 + a41 + a59 + ¼
+
a13 + a31 + a49 + ¼.
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. In the CA applet, observe the difference between the odd numbered coefficients of S9, P9 and I9. 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.

convergence plot

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.

5  Rates of Convergence

Repeatedly integrating equation (4) by parts we get
an = 1
nm
 2
p
 ó
õ 		1

-1 
 (1-x2)-1/2  f(m)(x)  Tn(x) dx.
(13)
Thus, if f(x) is m-times (m ³ 1) continuously differentiable in [-1,1] the above integral will exist and we can conclude that
an = O(n-m),        n = 1,2,¼.
(14)
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(x) is (m-1)-times differentiable a.e. (almost everywhere) in [-1,1] with its (m-1)th derivative of bounded variation in [-1,1].
Since the absolute value of each Tk is bounded above by one on [-1,1], it follows that the truncation error for the continuous expansion is bounded by the sum of the absolute value of the neglected coefficients:
|f(x) - SN (x) | £ ¥
å
n=N+1 
 |an|.
(15)
A similar bound, with an additional factor of two, holds for the interpolating partial sum:
|f(x) - IN (x) | £ 2 ¥
å
n=N+1 
 |an|
(16)
for x Î [-1,1]. From (14), (15), and (16) we conclude that
|f(x) - SN (x) | = O(N-m)
(17)
and
|f(x) - IN (x) | = O(N-m).
(18)
If f is infinitely differentiable the convergence is faster than O(N-m) 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 [-1,1] can be shown to be
O(r-N)
(19)
for some r > 1. In figure 2 the rather slow decay rate of the error with increasing N is illustrated for the absolute value function for which m=1. This can be contrasted with the rapid spectral convergence of the infinitely smooth function which is also illustrated in figure 2.
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.

convergence plot

Figure 2: Convergence of an infinitely differentiable function f2(x) = ecos(8x3+1) vs. convergence of a continuous function f5(x)=|x|.
Equation (13) allows us to conclude that if f is a polynomial of degree N then an = 0 for all n > N. For example, observe that the seventh degree polynomial in the CP approximation applet has ak = 0 (to within machine precision of 2.2204e-16) 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 O(N-1). Oscillations will be present near the discontinuity and will not diminish as N® ¥. Additionally, the oscillations will not even be localized near a discontinuity. This situation is referred to as the Gibbs Phenomenon. The Gibbs' oscillations are very visible in the Chebyshev approximations of the sign function in the CP applet.

6  Filters

Spectral filters can be used to enhance the decay rate of the Chebyshev coefficients [12]. The filtered Chebyshev approximation is defined as
FN ( x) = N
å
n=0 
 s æ
è 		n
N
 ö
ø 		an Tn(x)
(20)
where s is a spectral filter. A pth (p > 1) order spectral filter is defined as a sufficiently smooth function satisfying
s(0)
=
1
(21)
s(m)(0)
=
0        m = 1,2,¼,p-1
(22)
s(m)(1)
=
0        m = 0,1,¼,p-1.
(23)
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(x) has a discontinuity at x0 and setting d(x) = x -x0, the estimate
|f(x) - FN(x)| £ K
d(x)p-1Np-1
(24)
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
s(w)=exp(-aw2p)     p = 1,2,¼
(25)
of order 2p. In order for condition (23) to be satisfied, the parameter a is taken a = -lnem where em is defined as machine zero. On a 32-bit machine using double precision floating point operations, em=2-52 » 2.2204e-16 and ln(e) » -36.0437. The exponential filter is implemented in the CA applet. To apply the filter in the applet select the filter option on the approximation menu. The filter will be applied to the approximation that is marked in blue on the approximation menu. The order (the default is a fourth order filter) of the filter may be changed in the parameter dialog box which is available from the options menu. Select the sign functions and apply the filter. Observe the improved accuracy away from the discontinuity as well as the effect it has on the magnitude of the coefficients. Experiment with filters of various orders. Notice how a lower order filter smears or rounds off the discontinuity while a higher order filter does not completely remove the oscillations around the discontinuity.
The EP applet illustrates the strength of the damping applied in equation (20) to the coefficients ak from k=0,1,¼,N for filters of order 2 to 32.

7  Applet

Use the CA applet to further experiment with the concepts we have discussed. The example functions included in the applet:

8  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 (11), are well established as powerful methods for the numerical solution of PDEs with sufficiently smooth solutions.
Interpolation means that f(x), the function that is approximated, is already a known function. The terms collocation and pseudospectral are applied to global polynomial interpolatory methods for solving differential equations for an unknown function f(x). Detailed information on pseudospectral methods may be found in the standard references [1,2,3,4,5,11].
Many important PDEs have discontinuous (or nearly discontinuous) solutions. See the article [9] 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 may need to know the exact location of the discontinuities, include Spectral Mollification, Gegenbauer Reconstruction [6], Padé Filtering, and Digital Total Variation Filtering. Several postprocessing methods with applications are discussed in [10] with supporting web material at http://www.scottsarra.org/signal/signal.html. 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.

References

[1]
John P. Boyd. Chebyshev and Fourier Spectral Methods. Dover Publications, Inc, New York, second edition, 2000.
[2]
Claudio Canuto, M. Y. Hussaini, Alfio Quarteroni, and Thomas A. Zang. Spectral Methods for Fluid Dynamics. Springer-Verlag, New York, 1988.
[3]
D. Funaro. Polynomial Approximation of Differential Equations. Springer-Verlag, New York, 1992.
[4]
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, pages 1-54. SIAM, Philadelphia, 1984.
[5]
David Gottlieb and Steven A. Orszag. Numerical Analysis of Spectral Methods. SIAM, Philadelphia, PA, 1977.
[6]
David Gottlieb and Chi-Wang Shu. On the Gibbs phenomenon and its resolution. SIAM Review, 39(4):644-668, 1997.
[7]
J. Mason and D. Handscomb. Chebyshev Polynomials. CRC, 2003.
[8]
T. Rivlin. The Chebyshev Polynomials. Wiley, 1974.
[9]
S. A. Sarra. The method of characteristics with applications to conservation laws. Journal of Online Mathematics and its Applications, 3, 2003. http://www.joma.org/mathDL/4/?pa=content&sa=viewDocument&nodeId=389 (accessed March 1, 2005).
[10]
S. A. Sarra. The spectral signal processing suite. ACM Transactions on Mathematical Software, 29(2):1-23, 2003.
[11]
L. N. Trefethen. Spectral Methods in Matlab. SIAM, Philadelphia, 2000.
[12]
H. Vandeven. Family of spectral filters for discontinuous problems. SIAM Journal of Scientific Computing, 6:159-192, 1991.


File translated from TEX by TTH, version 3.66.
On 01 Mar 2005, 09:38.