Polynomial series versus sinc expansions for functions with corner or endpoint singularities

JOURNAL

OF COMPUTATIONAL

PHYSICS

64,

266-270 (1986)

Note Polynomial Series versus Sine Expansions for Functions with Corner or Endpoint Singularities In a review of methods that use “Whittaker cardinal” or “sine” functions, Stenger [l] shows that these basis functions-in combination with a change-ofvariable-are a powerful tool for approximating a function with weak singularities at the ends of the interval. Although Stenger himself is careful to note that the same optimal convergence rate can be obtained with other basis functions, he does not elaborate or give examples. In this note, we show that the change-of-variable-not the use of sine functions-is the key to successin coping with endpoint singularities. We explicitly construct approximations using mapped orthogonal polynomials’ which have the property of “exponential” or “infinite order” convergence for f(x) which has weak singularities at the endpoints but is regular on the interior of the interval. (We define what we mean by a “weak” singularity in Eq (1) below.) For simplicity, most of the analysis is limited to functions of one variable, but corner singularities for two-dimensional boundary value problems are a very important application which will be briefly discussed at the end. We shall standardize the interval to x E [ - 1, 1] and assume f(x) has singularities at the endpoints of the form f(x) = (1 - X2Y g(x)

(1)

where g(x) is free of singularities on [ - 1, l] and where CI> 0 so that the f(x) is bounded even at the branch points. (This boundedness off(x) is what we mean by a “weak” singularity.) If we expand f(x) in a Chebyshev series, Elliott [4] has shown that the coefficients (for LX> 0) are a,-0(1/n’+*‘).

This poor convergence*oeflicients

(2)

decreasing only as an algebraic function of

’ Because of the mapping, the terms of the series will be transcendental functions of the original coordinate x. This is what makes it possible to evade the theorem that polynomials in x could not give better than algebraic convergence because of the singularities.

266 OOZl-9991/86$3.00 Copyright 0 1986 by Academic Press, Inc All rights of reproductmn in any iorm reserved.

POLYNOMIALS

VERSUS

SING

267

SERIES

n-is the best possible for polynomials, and has prompted a search for alternatives. The successfuloptions use a series of terms which are transcendental functions of x. Stenger [l] showed that it is possible to create an approximation through a three-step procedure whose error decreasesexponentially rather than algebraically. The first step is to transform the interval x E [ - 1, 11 to y E [-co, cc] through the mapping x = tanh (ky)

(3)

where k is an arbitrary scale factor. The second step is to choose a grid spacing h and then approximate f( y[x]) through the “sine expansion” or “Whittaker cardinal” approximation, which is f(y)=

f f(3) sine(CY -31/h); ,=--co

YE

[--co, 001

where the “sin? function is defined by sine (z) 5 sin (7rz)/(rcz).

(5)

The final step is to truncate the infinite series in (4) so that we sum over a finite number of grid points N. Stenger [l] goes on to show that the sine expansion can be used to solve differential and integral equations, but for simplicity, we will assumef(x[ y]) is a known function. There are two sources of error in (4): a “grid-spacing” error because the interpolation points are a finite distance h apart and a “grid-span” error because we must truncate the infinite sum in (4), which implicitly restricts the grid points to some finite span of the interval in y. Stenger shows that minimizing the combined effects of these two sources of error requires taking hwO( l/N”*) where N is the total number of grid points. The error is then O(exp [ -PN”~]) for some constant p provided that c?> 0,

(6)

which is the condition that f(x) remain bounded at x = f. 1. However, the key to defusing the branchpoints is the tanh-mapping. If we substitute (3) into (l), we find

f(xCYI)=sech*“(k~)g(xCyl)

(7)

where the lack of singularities in g(x) on [ - 1, l] implies that g(x[y]) is bounded for all real y. Once the problem has been converted to the interval y E [ - co, co], one can use any of a wide variety of alternatives to approximate the function, not just the sine series (4). 581/64/l-18

268

JOHN P. BOYD

Boyd [3] shows that for a function like (7), the expansion in Hermite functions, $n(Yh will converge for all real y. As for the sine series, the error is O(exp [ -qN’/‘]) for some constant q. In terms of the original coordinate x, we are escaping the ineffectiveness of a polynomial series by using the transcendental basis functions $J y[x]) instead. A Hermite pseudospectral method would use a computational grid composed of the roots of the (N + 1)-degree Hermite function, which are separated by an average grid spacing of O(NP1”). (For a given N, the distance between neighboring grid points is variable, but does not differ much from the average.) Stenger [l] shows that h- O(N- l”) is usually optimum for the sine expansion as well. For either series, the mapping is crucial: the tanh function transforms a grid of points which is evenly (or almost evenly) spaced in y into a grid in x in which the gridpoints are clustered near the endpoints at x = + 1. The spacing of adjacent points in x decreasesexponentially with N near the branch-points, which is precisely what is needed to defeat these singularities. In contrast, the Chebyshev pseudospectral method (without a mapping of the coordinate) would use a grid with a nearestneighbor separation no smaller than 0( 1/N2); for any positive a, Eq. (2) shows that the highest calculated Chebyshev coefficient, aN, is proportional to a negative power of N also. Thus, a grid spacing which decreasesalgebraically fast with N gives an error that decays algebraically with N, too. In contrast, decreasing the grid spacing in x exponentially fast near x = f 1 gives an error which decreasesexponentially with N also in spite of the bounded endpoint singularities. Estimating the constants p and q in error formulas like (7) is difficult in general, but it is instructive to consider the simple example

f(y) = se& (C7@1”2~)

(8)

which is (7) for the special case a = l/2, k = (n/2)‘/*, and g(y) = 1. This choice of the scale factor k optimizes the convergence of the Hermite series and we find from [3] N) - O(exp [ - 1.772N”*] ). E Hermite(

(9)

Using Theorem 2.1 of Lund and Riley [2] to optimize the sine grid gives h = (27~/N)“~ where N is the total number of interpolation points and Esinc(N) w O(exp [ - 1.5708N”2]).

(10)

Thus, the series of orthogonal polynomials is superior (at least asymptotically as N + 03) to employing the sine expansion. It would be wrong to conclude, however, that the sine expansion is an inferior method. The difference between (9) and (10) is quite small, and sine functions are simpler and easier to program than Hermite series.The point is rather that the sine expansion is not the only choice.

POLYNOMIALS VERSUS SINC SERIES

269

Indeed, Boyd [S] is the basis for a third alternative for expanding functions with branch-points at the ends of the approximation interval. That earlier article discusses how problems on y E [ - co, co] can be efficiently solved by mapping the infinite interval into [ - 1, 1] and then applying Chebyshev polynomials. The map x = tanh (ky) is poor for this purpose because it will convert a function that decays exponentially in y into one with singularities at x = f 1 which spoil the convergence of the Chebyshev series as shown by (2).* It is amusing to see that the tanh (ky) map is useful when applied in the opposite direction to convert a function that really is singular at x = &-1 into one with exponentially convergent Hermite and sine series. This suggests the third alternative for expanding this class of functions: employing the tanh (ku) map to transform from [ - 1, l] to [-co, cc] and then using an algebraic map of the kind described in [S] to transform back to [ - 1, 11; the doubly-mapped Chebyshev polynomials will give exponential convergence. (Without the change of coordinates, of course, the Chebyshev series would converge algebraically with N.) The class of functions with weak endpoint singularities may seem rather special, but Lund and Riley [2] give a number of one-dimensional examples. Partial differential equations are a very rich source of examples since boundary value problems may have singular solutions even if the equation is linear and constant coefficient. The classic illustration is v*u= -1

(11)

on the square [ - 1, l] x [ - 1, l] with u E 0 on all four sides of the domain. Strang and Fix [6] show that U(X) has singularities of the form 24= Y*log (r) + less singular terms

(12)

in each corner where r is the radius of a local polar coordinate centered on the corner. Bowers and Lund [7] have reported success in solving problems like (10) using mapped sine expansions, but the transformed Hermite series and the doublemapped Chebyshev polynomials would probably be equally effective. The moral of this note is that there are many ways to deal with endpoint and corner singularities. The sine series of Stenger is simple and exponentially convergent, but there are alternative basis functions that are just as good. For all these pseudospectral methods, however, a change of coordinates which gives an exponential clustering of gridpoints near the singularities is essential.

ACKNOWLEDGMENTS This work was supported by NSF Grant OCE8305648. I thank Professors Frank Stenger and John Lund for helpful comments, and also the two anonymous reviewers. ZThere are some tricks that allow the tanh-mapping to be effective, at least occasionally, but the author is still investigating.

270

JOHN P. BOYD REFERENCES

F. STENGER,SIAM Rev. 23, 165 (1981). J. R. LUND AND B. V. RILEY, IMA J. Numer. Anal. 4, 83 (1984). J. P. BOYD, J. Comput. Phys. 54, 382 (1984). D. E. ELLIOTT, Math. Comput. 18, 274 (1964). 5. J. P. BOYD, J. Comput. Phys. 45, 43 (1982). 6. G. STRANG AND G. J. FIX, An Analysis of the Finire Element Method (Prentice-Hall, New York, 1. 2. 3. 4.

1971), Chapter 8. 7. K. BOWERS AND J. LUND, RECEIVED:

submitted for publication.

February 28, 1985; REVISED August 8, 1985 JOHN P. BOYD Department of Atmospheric and Oceanic Science University of Michigan 2455 Hayward Avenue, Ann Arbor, Michigan 48109