<<

. 76
( 136 .)



>>

the geometrically converging series

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

<<

. 76
( 136 .)



>>