Accuracy Enhancement for Higher Derivatives using Chebyshev Collocation and a Mapping Technique.
A new method is investigated to reduce the roundoff error in computing derivatives using Chebyshev collocation methods. By using a grid mapping derived by Kosloff and Tal-Ezer, and the proper choice of the parameter $\alpha$, the roundoff error of the $k$th derivative can be reduced from $O(N^{2k})$ to $O((N \lge)^k)$, where $\epsilon$ is the machine precision and $N$ is the number of collocation points.  This drastic reduction of roundoff error makes mapped Chebyshev methods competitive with any other algorithm in computing second or higher derivatives with large $N$. Several other aspects of the mapped Chebyshev differentiation matrix are also studied, revealing that    the mapped Chebyshev methods require much less than $\pi$ points to resolve a wave; the eigenvalues are less sensitive to perturbation by roundoff error; and  larger time steps can be used for solving PDEs.   All these advantages of the mapped Chebyshev methods can be achieved while maintaining spectral accuracy.
Introduction
In [5], we addressed the issue of roundoff error in computing the derivative using the Chebyshev
collocation methods. For details on these methods, see Canuto, et. al. [1], or Gottlieb and Orszag
[2]. We showed how to construct the Chebyshev collocation derivative with only O(fflN 2 ) roundoff
error, where ffl is the machine precision and N is the number of collocation points, by carefully
constructing the entries of the derivative matrix.
There are PDEs that involve higher derivatives than the first. For example, the viscosity term in
the Navier-Stokes Equation and the fourth derivative term in the Kuramoto-Sivashinsky Equation.
For these problems, the roundoff error in the k-th derivative will be O(fflN 2k ). This can ruin the
computed solution, even if k and N are not that large. This limits the applicability of Chebyshev
collocation methods to certain types of PDEs and values of N .
In this paper we investigate a way of modifying the Chebyshev collocation derivative which
reduces the roundoff error, and also improves its accuracy. As stated in [5], the roundoff error
cannot be reduced further for basis functions based on polynomials, since the roundoff error of
the Chebyshev collocation method already achieves the theoretical minimum. Therefore, to have
any hope of reducing the roundoff error, the polynomial basis functions must be replaced with
something else. One way to do this is to apply a coordinate transformation. More specifically,
the Chebyshev collocation points - are mapped to a new set of points
with a parameter ff. The transformation function g(-; ff) is one-to-one
and onto. The mapping is also applied to the polynomial basis functions, and so they are changed
as well.
In this paper, we will concentrate on a mapping of the form
In [4], Kosloff and Tal-Ezer showed that this mapping increases the minimum spacing \Deltax between
collocation points from O(N \Gamma2 ) to O(N \Gamma1 ), and argued that this reduces the spectral radius of the
differentiation matrix from O(N 2 ) to O(N ). Here, we will show that it has a similar effect on the
roundoff error. Moreover, the mapping not only reduces the roundoff error but also requires less
than - points per wave number as with standard Chebyshev methods.
In section 2, we introduce the Chebyshev collocation method and its derivative matrix. Section 3
illustrates the problem of large roundoff error when computing higher derivatives using the standard
Chebyshev collocation methods. We discuss the transformation of the Chebyshev collocation points
and its properties in section 4. The minimum roundoff error is estimated for the standard and
mapped differentiation matrix in section 5. In section 6, some numerical results of the standard vs.
mapped Chebyshev collocation methods are demonstrated. We study the effect of the mapping on
the eigenvalue spectrum of the matrix with Dirichlet boundary conditions in section 7. Section 8
gives a heuristic estimate for the resolution power of the mapped Chebyshev methods. We discuss
in section 9 some pitfalls and proper procedures in computing the mapped Chebyshev derivatives.
In this section, we present the Chebyshev collocation method. Consider the Chebyshev-Gauss-
Lobatto collocation points,
which are the extrema of the Nth order Chebyshev polynomial
Let v(x) be a smooth function in the domain x 2 [\Gamma1; 1]. Then v(x) is interpolated by constructing
the N order interpolation polynomial g j (x) such that g j
where u(x) is the polynomial of degree N and u . It can be shown that
where
Differentiation
The derivative of u(x) at the collocation points x j can be computed in many different ways. The
most obvious way to compute the derivative is the matrix-vector multiplication. The entries of
the Chebyshev derivative matrix D, are computed by taking the analytical derivative of g j (x) and
evaluating it at the collocation points x k for
Then the entries
of the matrix are, for
sin -
sin
and
The derivative of u(x i ) becomes ( 0 denotes derivative unless specified otherwise)
There are also O(N log 2 N) algorithms for computing the derivative that involve Fast Discrete
Fourier and/or Cosine Transforms. Their roundoff error properties are mostly the same as the
matrix-vector multiplication. We will concentrate on the matrix-vector multiplication algorithm
and only mention them occasionally.

Table

1: Absolute maximum error for the k-th derivative of
3 Accuracy in computing higher derivatives
In this section we show some numerical experiments which illustrate the large roundoff error of
the standard Chebyshev collocation method. We computed the first four derivatives of the test
function any mapping. The kth derivative is computed by
multiplying the vector \Gamma! by the derivative matrix D k times. Table 1 shows
the absolute maximum error for the first four derivatives of different numbers of
Chebyshev collocation points. It clearly illustrates the problem when one wishes to obtain higher
derivatives using the Chebyshev collocation methods. The absolute maximum error is fine for the
first derivative and acceptable for the second derivative if N ! 128. For the third and the fourth
derivatives, the error becomes unacceptable for N ? 64 and N ? 32, respectively.
The rate of growth of the roundoff error shown in table 1 is approximately O(N 2k
(See section 5 for more details). Similar conclusions can be reached for other functions, since this
situation is quite generic to the Chebyshev differentiation matrix D. Since the matrix D has
already achieved the minimum roundoff error possible, any further manipulation of the elements
of the matrix is useless. The situation seems hopeless unless other basis functions are found.
For example, we know that the Fourier collocation methods only have roundoff error growth like
for the k-th derivative. So if we modify the Chebyshev basis functions to resemble the
Fourier methods in some way, then maybe we can reduce the roundoff error. However, this has to
be done in such a way as to preserve the spectral accuracy that the standard Chebyshev method
has. One way to achieve this is to map the Chebyshev collocation points to some other set of
interpolation points. We discuss this in the next section.
4 Grid Transformation
In [4], Kosloff and Tal-Ezer derived a grid transformation that mapped the Chebyshev collocation
points to a new set of interpolation points. The grid transformation was defined as
sin
are the Chebyshev collocation points and x j are a new set of
interpolation points with the parameter ff 2 [0; 1].
Given a function u(x) and the derivative of u(x) can be evaluated as
du
d-
d-
dx (7)
where
d-
dx
and
sin
So computing at the grid points x or (7) can be rewritten as
\Gamma! u= MD \Gamma! u (10)
where the diagonal matrix M has elements
We define our new differentiation matrix.
We will now give a few important properties of the mapping g(-; ff).
ffl the mapping derivative g 0 (-; ff) is singular if
ffl For \Deltax one can show that
ffl if, for some positive number c [4],
then
-N
The coordinate transform results in some approximation error due to the fact that the transformation
function g is not entirely smooth (g(-) has singularities at In [4, section 4.8]
Kosloff and Tal-Ezer showed that if ff is chosen to be
then the approximation error is roughly ffl. By choosing ffl to be the machine precision, the error of
the coordinate transformation is essentially guaranteed to be harmless. With this choice for ff, in
the limit as N !1,
If we substitute this into (12) and assume that jln fflj AE -, this gives a minimum grid point spacing
of
\Deltax min (0)

Table

2: ff; \Deltax min ; \Deltax min (ff); \Deltax min (0) and their ratio as a function of N with
In table 2, ff is shown as a function of N for a fixed accuracy which is the
machine precision of the Cray C-90. It also lists several values for equation (15), the minimum
spacing with and without mapping and their relative size for each ff(N ).
Equation (15) implies that \Deltat = O(\Deltax min ) is inversely proportional to jln fflj. Hence, larger \Deltat
can only be achieved by reducing the accuracy ffl of the interpolation. This situation might not
be acceptable for those problems demanding high accuracy in the solution. However, Table 2 also
clearly shows that one can still get a slightly larger time step without any degradation of accuracy
of the approximation. In this paper, we are interested not on the issue of stability but on the
accuracy. This allows us to fix the accuracy requirement ffl to be the machine precision and choose
ff to be a function of N alone.
5 Estimates of the roundoff error
In this section we discuss two estimates of the roundoff error. In [5, section 4] we constructed
the following estimate of the roundoff error: We ignored all of the rounding errors except those
that occurred when the function vector u is stored in single precision. Further, we assumed that
those errors occurred randomly. Specifically, we assumed that they were uncorrelated, zero-mean,
and that their variance (i.e., their rms average value) was the machine precision. Under these
assumptions, the average roundoff error in the i-th component of the derivative is
ik
The average L 2 norm of the error is
and a reasonable estimate of the maximum error is
We can see that if e cheb; i is the average roundoff error for the Chebyshev collocation derivative,
then
is the corresponding error for the derivative when using the mapping. Since M ii is small near the
edges of the grid where e cheb; i is large, we expect the mapping to reduce the roundoff error.
Usually the error will be larger than this because of the other rounding errors that happen
during the computation. Occasionally it will be smaller because of a lucky cancelation of errors.
This estimate is useful as a lower bound. If an algorithm for computing Du is able to come
close to it, then there is no use in trying to improve the algorithm. Many algorithms for computing
the derivative are able to approach this bound.
The other estimate for the roundoff error is a heuristic estimate based on the minimum grid
point spacing. In finite difference methods with equally spaced points, the formula for the k-th
derivative is
where the p jk are the numbers in Pascal's triangle and are O(2 k ) in size. We see from the size
of the coefficients that rounding errors can be multiplied by roughly (2=\Deltax) k . We claim that the
situation is the same for spectral approximations of the derivative. Therefore a rough estimate of
norm of the error is
Taking the value of \Deltax min from equation (15), this gives
The corresponding error for the unmapped Chebyshev derivative is
6 Numerical Results on the mapped Chebyshev methods
In this section we compute the higher derivatives using the Chebyshev collocation with and without
the mapping technique. We use the test function sin(mx) for this purpose. We also
tested other functions, e.g. results. The relative
and absolute error between the exact and numerical results are computed by the L1 (pointwise
norm.
Here we note that there are two ways to computing the k-th derivative using Chebyshev collocation
methods with a mapping. One way is the multiply the matrix by the vector \Gamma! u
times, i.e. \Gamma! u (k)
\Gamma! u or \Gamma! u (k)
:). The other way is to form an expression
involving D l u(-) and g l (-; by differentiating u(g(x)) k times. For example,
\Gamma! u= \Gamma! g
where \Gamma! In numerical experiments, both techniques gave exactly the same results. In
this paper, we used the first method because of its simplicity in programming.
Mapping with Mapping
MV FFT CFT MV FFT CFT
First Derivative
N Second Derivative
Third Derivative
N Fourth Derivative

Table

3: Absolute maximum error for the first four derivative of
In

Table

3, we computed the maximum absolute error of the 1st, 2nd, 3rd and 4th derivative of
the function sin(2x). The computations used a Cray C90, whose machine precision is about 6:5 \Theta
\Gamma15 in single precision. We used three different algorithms, namely, Matrix-Vector Multiplication,
FFT-Recursion and CFT-Recursion with and without the mapping g(-; ff). The Cosine transform-
recursion algorithm is computed using the forward and inverse cosine transform subroutines FCR
and FCRINV from LARCLIB. This is a library local to NASA Langley Research Center, but
the source code is available. Both CFT and FFT algorithm are based on the same RFFTMLT
subroutine in the Cray Scientific Library Routine. The Matrix-Vector Multiplication algorithm
(MV) used the MXM subroutine in the same library package.
Based on Table 3, We can see that the solutions using the mapped Chebyshev methods are
uniformly better than standard one for all four differentiations. Another observation is that the
CFT algorithm is uniformly worse than both the MV and FFT algorithm in either the standard
or mapped category, especially for third and fourth derivative and large N . Moreover, the Matrix-Vector
Multiply (MV) algorithm and FFT Recursion algorithm (FFT) had a similar order of error.
We will concentrate our discussion on the data obtained by MV algorithm alone.
In

Figure

1 we compare the roundoff error of the mapped Chebyshev derivative with the error
estimates in equations (18) and (22) for the test function sin(4x). The vertical
axis of the graph is (Max error)=ffl. We see that the actual error comes fairly close to the lower
bound. We can also see that the heuristic error estimate is fairly close to real error and the lower
bound, although the dependence on k is not exactly right.10001000001e+071e+091e+11100 1000
error/epsilon
lower bound
error
error estimate

Figure

1: Maximum error vs. number of grid points for k-th derivative. Lower bound is equation
(18) and error estimate is equation (22). Error has been scaled by machine precision, i.e., the
vertical axis is (max err)=ffl, with
Mapping
With Mapping

Table

4: Absolute maximum error for the k-th derivative of
6.1 Accuracy for Functions with Boundary Layers
It is well known that the Chebyshev method has good resolution power for functions with boundary
layers. How well does the mapped Chebyshev method do on such functions? We tested functions
of the type
The gradient of the boundary layer can be adjusted by varying the parameter ffi ? 0. The smaller
the ffi becomes, the steeper the gradient will be near the boundary. The other parameters k and m
were set to be 2. We computed the derivatives for two different

Table

4 showed the absolute maximum error for case of 0:9. The gradient in this case is not
very large. The data showed that both algorithms performed equally well for small N . However, for
large N ? 32 and higher derivatives, the mapped Chebyshev method outperforms the unmapped
one.

Table

5 shows the absolute and relative maximum error for case of As seen in Table
5, the unmapped Chebyshev method performs better for small N - 32. We speculate that this
is because the mapped Chebyshev method will need more points to resolve such high gradient
near the boundary. It is quite self-evident that once the function is well resolved at N - 64, the
mapped Chebyshev methods again perform better. For higher derivatives, since the error is masked
by the large functional value, we examined the error in the relative sense for the third and fourth
derivatives. The relative accuracy using the mapping is rather good. The fourth derivative with
only :15% error with the mapping rather than 430% error of the standard Chebyshev
methods.
Absolute
Mapping
With Mapping

Table

5: Absolute and relative maximum error for the k-th derivative of
7 Eigenvalue Spectrum of the Differentiation Matrix
In this section, we analyze the eigenvalue spectrum of both the standard and mapped Chebyshev
Matrices D and D respectively, subject to the Dirichlet boundary conditions It is
equivalent to deleting the first row and column of the matrix D and D. Their eigenvalues are
computed by the IMSL subroutine EVLRG.
It is an important issue because the (time-) stability condition (\Deltat fixed; t ! 1) of any
explicit time marching scheme, for example, Runge-Kutta or Adams-Bashford scheme, is controlled
by its largest eigenvalues.
The eigenvalues of D and D can be separated into two groups; Those with large imaginary
parts, i.e. larger than O(N) which we will call outliers, and all the others. The outliers are
stable to perturbations of the matrix. In [3], Trefethen and Trummer found that the non-outlying
eigenvalues of D were extremely sensitive to perturbations of the matrix elements, and simply
storing the elements of the matrix in single or double precision would completely change the
spectrum of the matrix. The changed eigenvalues would have large negative real parts of order
O(N 2 ), as N !1. However, they showed that for any given precision ffl, then for non-outlying
eigenvalues - with negative real part Re- of those eigenvalues are
lost. In [3], they simulated the perturbation by lowering the precision of the arithmetic. Here the
perturbation comes from representing the matrix elements in finite precision.
For figures (2.a) and (2.b), show the eigenvalue spectrum of the standard Chebyshev
differentiation matrix D in 32 and 64 bits precision arithmetic respectively. Both cases exhibited
spurious eigenvalues with large negative real part in the form of an arc to the left of 1
the largest negative real part computed by bits precision is much larger than those computed
with 64 bits precision. For any fixed N , the only way to eliminate those spurious eigenvalues is
to use much higher precision arithmetic. This is undesirable and impractical for the numerical
solution of PDEs.
However, the eigenvalues of the mapped Chebyshev differentiation matrix D with the Dirichlet
boundary conditions behave very differently. For figures (2.c) and (2.d), showed the
eigenvalue spectrum of the standard Chebyshev differentiation matrix D in 32 and 64 bits precision
arithmetic respectively. For the 32-bit precision matrix, there are a few eigenvalues of D to the
left of the line
ffl. in the form of a small arc. For the 64-bits precision, there are none.
There are a few (less than 10) eigenvalues that have negative real part - that are less than 1ln ffl.
Hence we can conclude that the stability of the numerical scheme based on the mapped Chebyshev
collocation methods is much less sensitive to the precision of the arithmetic than the unmapped
method.
An important observation is that even for 64, the largest eigenvalues are considerably
smaller in the mapped case than in the unmapped case. Also the largest eigenvalue of the single-
precision-mapped matrix is smaller than the largest eigenvalue for the double-precision-mapped
matrix. The largest eigenvalue of D is (\Gamma91:9; \Sigma351:977) and its absolute value is e
The largest eigenvalue of D is (\Gamma52:1; \Sigma207:987) and its absolute value is e 214:426. The ratio
which is very close to the value 1:732 in Table 2 for
-100200(a) bits precision
(b) 64 bits precision
With Mapping
(c) bits precision
With Mapping
(d) 64 bits precision

Figure

2: Eigenvalues spectrum of the standard and mapped Chebyshev differentiation matrix D
(a,b) and D (c,d) respectively, in 32 (a,c) and 64 (b,d) bits precision arithmetic. the ff of
the mapping computed using the different values of ffl for the 32-bit and 64-bit precision matrices.
As seen in table 6, the most negative eigenvalue of D lying on the real axis grows like O(N 1:25 ),
at least for the values of N in the table. This is substantially smaller than the O(N 2 ) growth seen
when using the standard differentiation matrix D.
This result of O(N 1:25 ) is actually quite surprising, for the following reason: recall that the trace
Detail of spectrum of single precision mapped derivative matrix

Figure

3: Trajectories of some eigenvalues of the mapped Chebyshev derivative matrix as the size
of perturbations of the matrix increases from ffl double precision to ffl single precision The ff is computed
according to a 32-bit precision ffl, computations are in IEEE double precision.
r 1.27 1.33 1.20 1.18

Table

The most negative eigenvalue of D lying on the real axis, - 0 and the local exponent of its
growth, i.e., O(N r ).
of a matrix is defined as
are the eigenvalues
of M .
Both D and D are anti-centrosymmetric (ACS) matrices, i.e., they satisfy equation (4). The
trace of an ACS matrix can easily be shown to be zero. Suppose we modify these matrices by
throwing away the first row and column, and call these matrices ~
D and ~
D. Then
and
ff
Since M 00 is supposed to be small, the mapping forces the eigenvalues of D to be closer to the
imaginary axis than those of D. Letting taking the limit as N !1,
The trace of ~
D is real, so
Also, the eigenvalues of ~
are experimentally always found to be in the left half-plane. If some of
them were not, it would make D useless for solving parabolic PDEs, so for the rest of the section,
we will assume that
Now we can use our knowledge of tr( ~
D) to say some things about the spectrum of ~
D.
Let us consider two extreme arrangements of eigenvalues which satisfy equations (30) and (31).
First, suppose all of the eigenvalues but one lay on the imaginary axis. Then the real part of
the remaining one would be tr( ~
D). This maximizes suppose that all of the
eigenvalues have the same real part, so Re-
D)=N; 8i. This minimizes
reality falls between these extremes, so
D) -
D)
-jln fflj: (33)
If the O(N 1:25 ) growth rate for the eigenvalues of ~
D lying on the real axis shown in table 6 holds
for all N then for some large enough N , it will exceed the limit given by equation (32). Therefore
the O(N 1:25 ) growth rate cannot be the asymptotic rate as N ! 1. The asymptotic growth rate
of the real part of any eigenvalue of D can be at most O(N ). In table 6, we see that r is decreasing.
Presumably, as N!1, r approaches 1.
In the case of the Legendre derivative matrix, Trefethen and Trummer found that the sensitivity
of the eigenvalues of that matrix effectively raised its spectral radius from O(N) to O(N 2 ). In
infinite precision, the spectrum of the Legendre matrix is similar to that of the Chebyshev matrix,
except that it doesn't have the O(N 2 ) outliers. Heuristically, we expect that the mapping reduces
the spectral radius of ~
D to O(N jln fflj). A reasonable thing to fear is that eigenvalue sensitivity might
create eigenvalues with large negative real parts, increasing the spectral radius above O(N jln fflj).
This cannot happen, at least not by generating eigenvalues with large negative real parts. Equation
prevents any eigenvalue from having real part bigger than O(N jln fflj).
?From equation (33), at least some of the eigenvalues have to be to the left of the line
1jln fflj. It is not clear if this means that some of the eigenvalues have to be changed by roundoff
error, since the eigenvalues farthest to the left appear to be the outliers, which are stable anyway.
In

Figure

3 we show how the eigenvalues of the mapped derivative matrix change as the matrix
perturbations increase in size. We computed the eigenvalues of the matrix
~
where the entries of the matrix E are random numbers, uniformly distributed in the interval
the curves shown in Figure 3 are the paths taken by the eigenvalues of ~
D t in the complex
plane as t went from from 10 \Gamma15 to 5 \Theta 10 \Gamma8 . The mapping used a value of ff consistent with
. All of the computations were done in IEEE double precision using the LAPACK
subroutine DGEEVX. The ends of the curves with the diamond are the small-t ends. Only the
non-outlying eigenvalues with non-negative imaginary part are shown. It can be clearly seen that
at least a few of the eigenvalues are effected by roundoff error. The diamonds that don't appear
to have a curve attached to them are eigenvalues which are not affected by perturbations.
This contradicts the statement in [4], that the eigenvalues of the mapped Chebyshev matrix
will not be affected by the roundoff error. They are, but much less than for the unmapped matrix.
Resolution power of the mapped Chebyshev methods
The standard Chebyshev method requires at least - points per wave before it begins to resolve an
oscillatory function like sin(m-x). How many points does the mapped Chebyshev method require?
In the limit as ff ! 1, the grid points become equally spaced, the Chebyshev polynomials become
cosine functions, and resolution requires two points/wave. However ff, as a function of N , never
reaches one, so its not clear what really happens.
We claim (without proof) that the grid begins to resolve a function when the density of points
reaches two points per wave in the center of the domain, where the density of points is lowest. This
implies that
if the function cos(mx) is to be resolved. In the limit as N !1, this means that
per wave are required.
We computed the normalized L 2 Error of the function cos(mx) for many m and N in figure
4. The x axis represents the inverse of the wave number m=- which is normalized by the factor
?From the previous section, we have
Examining the data carefully, for any fixed wave number m=-, the error decays exponentially
fast to ffl as N increases and as soon as N=(mg 0 (0; ff)) ? 1. This implies that a minimum of g 0 (0; ff)-
points is needed to resolve a wave. Table 7 shows some typical values of g 0 (0; ff)- for various sample
N . ffl is fixed at 6:5 \Theta 10 \Gamma15 . For large N , only about two points per wave are needed. This is close
to the performance of Fourier methods.

Figure

4: The normalized L 2 Error of the function cos(mx) for many m and N

Table

7: Number of points/wave vs. N
9 How not to compute the derivative
Since computing the k-th derivative with k and N large can result in very large roundoff errors, it
is worthwhile to work extra-hard to minimize them. In this section we discuss several non-obvious
sources of roundoff error and show how to avoid them.
9.1 Ill-conditioning of the computations associated with the mapping
Recall that the grid transformation function is
ff. and then the mapped grid points are x are the
unmapped Chebyshev grid points. The derivative of g(-) is
When jff-j is near 1, jg 0 (-)j is large. This can be a problem. The roundoff errors associated with
storing the - i in single precision will result in errors in x i of size O(fflg 0 (- i )). The largest values of
1. The computation of the transformation derivative M ii is also
ill-conditioned. To avoid excessive errors resulting from these problems, both x i and M ii need to
be computed in extended precision before being stored in single precision. Failure to do this can
result in loss of about 1 digit of precision for N larger than about 100.
9.2 Ill-conditioning of the matrix multiply used to form D k
The most obvious way to obtain the k-th derivative matrix D k is by constructing the first derivative
matrix D and doing multiplications. This is a bad idea; the matrix multiplication is
ill-conditioned. The sum that must be evaluated to obtain the entries of the 2nd derivative matrix
is
D (2)
(i), and ignoring the terms in the sum. The terms in the sum
change sign twice; when j. This is a problem. It means that the partial
sums can be much larger than the final sum. The roundoff error is basically proportional to the
size of the partial sums, see [6]. This can be partially fixed by adding up the terms in the sum in
random order, rather than in order of increasing k. But a better idea is simply to use closed-form
expressions for the entries of the matrix. Of course, there are no closed-form expressions for the
entries of the matrix (MD) k . This is a strong argument in favor of using the long form of the
mapped derivative matrix, i.e., according to equation (24).
9.3 Cosine transform Algorithms
There are other algorithms for computing this matrix-vector product. One involves doing a discrete
cosine transform on u, doing O(N) operations on the transformed data, and then an inverse cosine
transform. One way of computing the cosine transform is to symmetrically extend u into a vector
of length 2n and do an FFT on the longer vector.
Two faster alternative are the algorithms described in Dollimore [8], and in Cooley, Lewis, and
Welch [7]. They both involves an O(N) preprocessing step, an FFT on a vector of length n, and an
O(N) post-processing step. Unfortunately the pre- and post-processing steps are ill-conditioned,
and these algorithms have O(N ffl) roundoff error, compared to O(ffl) roundoff error for the symmetric
extension algorithm. When used to compute the Chebyshev derivative, this gives O(N 3 ffl) roundoff
error rather than O(N 2 ffl).
Swarztrauber [9] describes a cosine transform which appears to be free of roundoff error prob-
lems, and at least as fast as the Dollimore algorithm. However, there are many publicly-available
subroutines which encode the Dollimore algorithm, and we have not been able to find any examples
of Swarztrauber's algorithm.

Acknowledgments

The authors would like to thank Prof. Alvin Bayliss for bringing this subject to the authors'
attention. Computing time on the Cray YMP was provided by the US Army Corps of Civil
Engineering, Waterway Experiment Station, High Performance Computer Center. The first author
would like to acknowledge the support of this research by AFOSR grant 93-1-0090, DARPA grant
N00014-91-J-4016 and NSF grant DMS-92-11820.



--R

Spectral Methods in Fluid Mechanics
Numerical Analysis of Spectral Methods: Theory and Applications
An Instability Phenomenon in Spectral Methods
Modified Chebyshev Pseudospectral Methods With O(N
Accuracy and Speed in Computing the Chebyshev Collocation Derivative
"The Accuracy of Floating Point Summation"
"The Fast Fourier Transform Algorithm: Programming considerations in the Calculation of Sine, Cosine, and Laplace Transforms"
"Some Algorithms for use with the Fast Fourier Transform"
"Symmetric FFTs"
--TR

--CTR
J. Mead , B. Zubik-Kowal, An iterated pseudospectral method for delay partial differential equations, Applied Numerical Mathematics, v.55 n.2, p.227-250, October 2005
R. Barrio , J. M. Pea, Numerical evaluation of the pth derivative of Jacobi series, Applied Numerical Mathematics, v.43 n.4, p.335-357, December 2002
Elsayed M. E. Elbarbary , Salah M. El-Sayed, Higher order pseudospectral differentiation matrices, Applied Numerical Mathematics, v.55 n.4, p.425-438, December 2005
W. Lyons , H. D. Ceniceros , S. Chandrasekaran , M. Gu, Fast algorithms for spectral collocation with non-periodic boundary conditions, Journal of Computational Physics, v.207 n.1, p.173-191, 20 July 2005
differentiation matrix;tal-ezer mapping;chebyshev collocation;roundoff error