∞

1

a2n T B2n (y 1/3 ) (17.55)

=

2 )1/3

(1 + y n=0

which is equivalent to mapping the problem from y to the trigonometric coordinate t using

y = L cot(t3 ).

4 Theindexing convention for the SB has been adjusted from the ¬rst edition of this book so that SB0 is the

lowest basis function and the even degree SB are symmetric with respect to y = 0.

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

366

Similarly, a change-of-coordinate will redeem Hermite and sinc series if the mapping

transforms algebraic decay into exponential decay as noted by Eggert, Jarratt, and Lund

(1987). An example on the semi-¬nite interval:

y = exp(x) ’ 1 y & x ∈ [0, ∞] (17.56)

Decay as O(1/y), for example, becomes decay as O(exp[’x]).

Nevertheless, if one is forced to accept the burden of a mapping, it is far simpler to solve

the transformed problem using a Fourier series rather than something more complicated

like Hermite or sinc expansions. Consequently, the earlier assertion remains true: the best

way to approximate algebraically decaying functions is to use a rational Chebyshev series,

which is a Fourier series with a change-of-coordinate. However, the possibility of applying

sinc and Hermite functions even to “bad” cases like algebraic decay is again a reminder of

the tremendous power of mappings.

17.10 Numerical Examples for the Rational Chebyshev

Functions

In oceanography, the steady north-south current of the “Yoshida jet” solves

vyy ’ y 2 v = y (17.57)

For large y,

1

v(y) ∼ ’ (17.58)

y

This slow, algebraic decay with |y| destroys both Hermite series and sinc expansions. Be-

cause v(y) is an antisymmetric function with odd powers of 1/y, Table 17.6 implies that the

T Bn (y) will fail, too. However, the mapping from y to t will still succeed provided that

we use a sine series instead of the usual cosine series to solve the transformed differential

equation. Fig. 17.7 shows the exact solution and the 2-point, 3-point, and 4-point colloca-

tion approximations. The N = 6 series listed in the table is indistinguishable from the exact

solution to within the thickness of the curve.

Table 17.7 lists the coef¬cients for various numbers of interpolation points. As always,

the error is the sum of two contributions. The truncation error is the neglect of all coef-

¬cients an with n > N ; the missing coef¬cients are represented by dashes. In addition,

those coef¬cients that are calculated differ from those of the exact solution. As we read

a row from left to right, however, we can see this “discretization error” decrease as each

coef¬cient converges to its exact value as N increases.

Note that because v(y) is antisymmetric, only odd basis functions (sin(t[y]), sin(3t[y]),

sin(5t[y]), etc.) were used in the expansion. All the collocation points were on the interior

of y ∈ [0, ∞] (instead of [’∞, ∞]) for the same reason.

17.10. NUMERICAL EXAMPLES: RATIONAL CHEBYSHEV FUNCTIONS 367

Figure 17.7: The latitudinal velocity v(y) for the “Yoshida jet”. The map parameter L = 3.

SBn were used instead of the more common T Bn to capture the very slow 1/y decay of

v(y) with spectral accuracy.

Table 17.7: The coef¬cients of the spectral series for the Yoshida jet v(y) as computed using

various numbers of collocation points with a basis of odd SB functions.

n 2 pts. 3 pts. 4 pts. 7 pts. 21 pts.

1 -0.47769 -0.388976 -0.394428 -0.395329 -0.395324

3 0.199475 0.179639 0.176129 0.179237 0.179224

5 ”” -0.050409 -0.04471 -0.051163 -0.051147

7 ”” ”” 0.001399 0.002620 0.002651

9 ”” ”” ”” 0.003856 0.003740

11 ”” ”” ”” -0.000541 -0.000592

13 ”” ”” ”” -0.000452 -0.000446

15 ”” ”” ”” ”” 0.000060

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

368

The second example is

uyy + n(n + 1) sech2 (y) + 1 u = Pn (tanh[y]) (17.59)

where n is an integer ≥ 0 and Pn (x) is the n-th Legendre polynomial.

Eq. (17.59) is merely Legendre™s differential equation transformed from x ∈ [’1, 1]

to y ∈ [’∞, ∞] via y = arctanh(x). On the ¬nite interval, the differential equation is

singular at both endpoints. The boundary conditions are that the solution must be analytic

at x = ±1 in spite of the singularities. The Legendre polynomials are the only solutions

that satisfy this requirement, and only when n is an integer.

The solutions of the transformed equation (17.59) asymptote to the constant 1 as y ’ ∞

for all n and to either 1 (even n) or -1 (odd n) as y ’ ’∞. The coef¬cients of the Hermite

series for such a function would decay as 1/n1/4 at the fastest. However, a rational Cheby-

shev expansion is extremely effective. The boundary conditions are natural (behavioral)

both before and after the arctanh-mapping, so it is not necessary to impose any constraints

on the T Bn (y).

Fig. 17.8 compares the exact solution (solid line) with the sixteen-point pseudospec-

tral approximation for the case n = 12. Because P13 (tanh[y]) is symmetric about y = 0,

only symmetric basis functions were used, and the sixteen-point approximation included

T30 (y). The approximation for N = 38 (20 collocation points, all with y ≥ 0), was indistin-

guishable from the exact solution to within the thickness of the curve.

Figure 17.8: Solid curve: P12 (tanh[y]) where P12 (x) is the Legendre polynomial. Dotted

curve: numerical approximation by solving the transformed Legendre differential equation

using 16 symmetric basis functions, i. e. T B0 (y), T B2 (y), T B4 (y), . . . , T B30 (y).

17.11. SEMI-INFINITE INTERVAL: RATIONAL CHEBYSHEV T LN 369

17.11 Rational Chebyshev Functions on y ∈ [0, ∞]: the T Ln (y)

These basis functions are de¬ned by

T Ln (y) ≡ Tn (x) = cos(nt) (17.60)

where the argument of the rational Chebyshev functions, y, is related to the argument of

the Chebyshev polynomials, x, and the argument of the trigonometric functions, t, via

y’L

L(1 + x)

” (17.61)

y= x=

1’x y+L

t y

”

L cot2 = 2 arccot (17.62)

y= t

2 L

The ¬rst few functions are given explicitly in Table 17.8. The pseudospectral interpolation

points are

(2i ’ 1) π

ti

” ti ≡

y = L cot2 (17.63)

i = 1, . . . , N + 1

2 2N + 2

Table 17.8: Rational Chebyshev functions for the semi-in¬nite interval: T Ln (y).

n T Ln (y; L)

0 1

(y ’ L)/(y + L)

1

(y 2 ’ 6yL + L2 )/(y + L)2

2

(y ’ L)(y 2 ’ 14yL + L2 )/(y + L)3

3

(y 4 ’ 28Ly 3 + 70L2 y 2 ’ 28L3 y + L4 )/(y + L)4

4

We will discuss these functions only brie¬‚y because their characteristics are very similar

to those of the rational Chebyshev functions on y ∈ [’∞, ∞] and Boyd (1987b) gives a

thorough treatment. If f (y) is rational with no singularity at ∞, then the T Ln (y) series will

have geometric convergence. Usually, however, the solution will be singular at in¬nity so

that convergence will be subgeometric.

If the differential equation has polynomial or rational coef¬cients (in y), the Galerkin

matrix will be banded.

If f (y) asymptotes to a constant or decays algebraically with y as y ’ ∞, then the

T Ln (y) series may still have an exponential rate of convergence. A series of Laguerre

functions, however, will converge very poorly unless f (y) decays exponentially for y 1.

In a modest difference from the T Bn (y) case, the map (17.62) shows that a cosine series in

t will give rapid convergence if f (y) has an asymptotic series as y ’ ∞ in inverse negative

powers of y. The sines are needed only to represent half-integral powers of y.

Boyd (1982b) offers guidelines for optimizing the map parameter L, but some trial-

and-error is inevitable. (The simplest criterion is: choose L to roughly equal the scale of

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

370

variation of the solution.) Fortunately, the accuracy is usually quite insensitive to L so long

as it is of the same order-of-magnitude as the optimum value. In the next section, we offer

an illustration.

17.12 Numerical Examples for Mapped Chebyshev

Methods on a Semi-In¬nite Interval

Example #1: Laguerre Eigenvalue Problem

11»

’’ (17.64a)

uyy + u=0

y4y

u(y) analytic and bounded at both y = 0 and y = ∞ (17.64b)

where » is the eigenvalue. The exact solution is

y

exp ’ y L1 (y) (17.65a)

u = n

2

≡ n, n an integer ≥ 0 (17.65b)

»

where L1 (y) denotes the ¬rst order Laguerre polynomial of degree n. Since the boundary

n

conditions are “behavioral” at both endpoints, unmodi¬ed T Ln (y) are a good basis. The

(N + 1) — (N + 1) pseudospectral matrix always has (N +1) eigenvalues, but the (N +1)-st

eigenfunction of the original differential equation is always oscillating too rapidly to be

resolved by interpolation using only (N +1) points. Consequently, we face the dif¬culty

” inevitable for any eigenvalue problem, whether on a ¬nite or unbounded interval ”

that the lowest few eigenvalues of the matrix are good approximations to those of the

differential equation whereas the larger matrix eigenvalues are numerical artifacts.

Fig. 17.9 shows the number of “good” eigenvalues computed with 41 collocation points

as a function of the map parameter L where a “good” eigenvalue is arbitarily de¬ned as

one within 0.05 of the exact eigenvalue (17.65b). Because N is large and the eigenfunctions

are entire, the theory of Boyd (1982a) argues that L should be fairly large, and the graph

con¬rms it. (Loptimum = 32.) It is reassuring, however, that taking L too large or small by

a factor of 2 still yields at least N/3 “good” eigenvalues. It is striking that the lowest two

eigenvalues are accurately calculated for L anywhere between 1 and 512, a range of three

orders-of-magnitude!

Example #2: Global Expansion of the K-Bessel Function, i. e.

f (r) ≡ r K1 (r) r ∈ [0, ∞]

on (17.66)

where K1 (r) is the usual imaginary Bessel function. Library software for evaluating this

function normally employs (i) a power series representation with logarithmic terms for

small r and (ii) an asymptotic series in inverse powers of r for r 1. Our goal is much

more ambitious: to use mapping to obtain a single expansion that is accurate for all r ≥ 0.

The Bessel function K1 has a simple pole at the origin, but the factor of r in (17.66) re-

moves this so that f (r) is ¬nite. However, no multiplicative factor can remove the log(r),

which itself multiplies an in¬nite series in r. Fortunately, since f (r) is bounded at the ori-

gin, we can use the exponential mapping trick introduced by Stenger. Since K1 (r) decays

as exp(’r) for large r, it is convenient to use a mapping which is linear in r for large r

since we do not have any unusual problem at ∞. Therefore, we ¬rst apply the change of

coordinate

r = arcsinh(ey ) (17.67)

17.12. NUMERICAL EXAMPLES: CHEBYSHEV FOR SEMI-INFINITE INTERVAL 371

Figure 17.9: The number of “good” eigenvalues of the Laguerre equation as a function

of the map parameter L. A “good” eigenvalue is arbitrarily de¬ned as one whose abso-

lute error is less than 0.05. All runs used 41 collocation points and the basis functions

T L0 (y), . . . , T L40 (y) on the interval y ∈ [0, ∞].

to transform to y ∈ [’∞, ∞] and then apply the T Bn (y). Although it is a little confusing

to use the Chebyshev rational functions on the in¬nite interval in solving a problem which

was originally de¬ned on a semi-in¬nite domain, this example does illustrate the power of

mappings.

The end result, after conversion back to the original coordinate r, is an approximation

of the form

0.534 ’ 0.6 T B1 (log[sinh(r)]) ’ 0.068 T B2 (log[sinh(r)]) (17.68)

r K1 (r) =

+ 0.125 T B3 (log[sinh(r)]) + 0.032 T B4 (log[sinh(r)])

’ 0.032 T B5 (log[sinh(r)]) + terms smaller than 0.01

Because we are resolving a logarithmic singularity at r = 0, a good choice of L must be

much smaller than for the previous example. Fig. 17.10 shows the result with map param-

eter L = 4: the logarithm of the error to base 10 is graphed against N , the truncation.

The approximation is almost too good. If the series converged geometrically, then the

curve would approximate a straight line. The actual graph is a fairly good approximation

to a straight line ” but one can rigorously prove that the rate of convergence must be subge-

ometric because r K1 (r) is singular at both endpoints. A careful look shows that the curve is

beginning to ¬‚atten out a bit near the right edge of the graph; the slope would tend to zero

from above if we extended the graph to suf¬ciently large N . This would be rather foolish,

however, because N = 24 already gives six decimal place accuracy, and within the range

illustrated, the subgeometrically convergent series does a remarkably good job of imitating