Fourier Method :

Since the Fourier differentiation and smoothing matrix are Toeplitz matrices, i.e. D(k,j) = D(k-j), only the first 2N-1 elements need to be stored. The Matrix-Matrix Multiply algorithm reduces to SAPXY (dot product).

Chebyshev and Legendre Methods :

Multiplies a matrix by each data vector, using BLAS calls. The matrix is very structured, and for large N, special algorithms are usually faster. But for small enough N, matrix multiplication is sometimes the fastest choice.

This is an O(N²) Algorithm.

This algorithm takes advantage of the (anti) centro-antisymmetry property of the differentiation and filtering (smoothing) matrices.

1. The data vectors are decomposed into the even and odd parts.

2. Two smaller matrices modified from the original differentiation matrix multiply the even and odd parts of the data.

3. An Even-Odd reconstruction procedure is then employed to obtain the desired final result.

This has roughly half as many operations as the matrix-multiply algorithm.

This is an O(N²)/2 Algorithm.

For small to moderate N, this is usually the fastest choice.

Fourier Transform Algorithm (FFT) :

The Transform Algorithm consisted of a Fast Fourier Transform, an O(N) operations to modify the Fourier coefficients and an Inverse Fast Fourier Transform.

Depending on the computational platform, a native optimized FFT (recommended) or one of the generic FFT library from NetLib can be used.

This is O(N log N) operations. Roundoff error usually is not an issue for the Fourier method.

Chebyshev Transform-Recursion Algorithm (CFT) :

The Transform Algorithm consisted of a Fast Cosine Transform (CFT), an O(N) operations recursive operation on the transformed data and another Fast Cosine Transform.

Depending on the computational platform, a native optimized CFT (recommended) or one of the generic CFT library from NetLib can be used.

This is O(N log N) operations, while the other two algorithms (MXM, EOD) are both O(N²).

For large enough N, this is always the fastest.

Roundoff Error :

Roundoff error of Transform-Recursion Algorithm is usually somewhat worse than for the other two algorithms.

The roundoff error of Cosine Transform from IBM ESSL v2 is good, as good as the Matrix-Matrix-Multiply Algorithms.

VFFTPack from NetLib was 2 or 3 digits worse accuracy for large N.

Speed :

Speed strongly depends on the quality and availability of optimized library software.

The critical value of N, where the Transform-Recursion Algorithm becomes faster than the Even-Odd Algorithm, can vary a great deal from one computer to another. On vector supercomputers it is usually between 50 and 150. On workstations it might be less. Or not.

For functions that are either even or odd, appropriate fast algorithms are employed to take advantages of such facts.

For Full Matrix-Multiply Algorithm, it is deferred to an Even-Odd Decomposition Algorithm immediately.

For the Even-Odd Decompositions Algorithm, only the even or odd part of the matrices is used for computations.

Since in most cases, both even and odd functions exist together in the computation, the even and odd parts of the differentiation/smoothing matrices are both computed and kept.

For the Transform Algorithm, it is a bit more complicated.

To compute derivative for the Fourier Method,

1. Fast Cosine (Sine) Transform is used to compute the real (imaginary) part of the Fourier coefficients from the Even (Odd) data.

2. The resulting Fourier coefficients are multiplied by ik.

3. Fast Sine (Cosine) Transform returns the derivative of the function with Odd (Even) symmetry.

To compute derivative for the Chebyshev Method,

1. Fast Cosine (Quarter-Wave) Transform is used to compute the Even (Odd) number of the Chebyshev coefficients from the Even (Odd) data.

2. The resulting Chebyshev coefficients for the derivative are computed recursively as in the standard case.

3. Quarter-Wave (Cosine) Transform returns the derivative of the function with Odd (Even) symmetry.

Similar Algorithms are also used for filtering.

Chebyshev and Legendre Methods :

In the classical Chebyshev and Legendre collocation methods, the grid points are much denser near the ends of the domain than in the middle. This is necessary for the polynomial interpolation process to be well behaved, but leads to problems with time-stepping and condition numbers. The minimum spacing between grid points is O(N^-2) for N grid points.

For complete description of the mapping capability of the library, please consult the postscript file Mapping.ps or the Acrobat PDF file Mapping.pdf.

Angular mapping for the Fourier Method :

This is a coordinate transformation that concentrates points at a specific Angle for the Fourier Method.

x = ArcTan((1-a²)*Sin(y-b)/((1+a²)*Cos(y-b)-2*a)

where a=alpha, |a|<=1 and b=Angle of concentration of grid points. The larger alpha(a) is, the more points will be concentrated near Angle(b).

Kosloff-Tal-Ezer mapping for Chebyshev and Legendre Methods :

This is a coordinate transformation that reduces the variation in grid point spacing for the Chebyshev and Legendre Methods.

x = ArcSin(alpha*y)/ArcSin(alpha)

The "strength" of the mapping depends on a parameter alpha :

    alpha =  0 :  No mapping, original Chebyshev/Legendre grid.
    alpha -> 1 :  Equally spaced points. 

The mapping creates some approximation error, which can be reduced by reducing the strength of the mapping. In particular, the mapping error can be set to the precision of the arithmetic, which is the default here.

Additional mapping for Chebyshev and Legendre Methods :

Several Built-in Coordinate Transformations are included in this release. A subroutine is also provided to allow user to add any mapping function of their own choice. (See the postscript file Mapping.ps or the Acrobat PDF file Mapping.pdf for detail.)

In summary, numbers of mapping function for Finite, Semi-Infinite and Infinite domain are included.

For Finite Domain, Kosloff-Tal-Ezer and Tangent Mapping are available.

For Semi-Infinite Domain, it has Algebraic, Exponential and Logarithm Mapping.

For Infinite Domain, Tangent, ArcTanh and Algebric Mapping are included.

For non-linear PDE and/or non-smooth solutions, numerical scheme employing spectral methods may exhibit instability.

Filtering (Smoothing) allows one to achieve numerical stability.
Four different filters are incorporated into the package :

1. Exponential Filter (Variable Order)

   For k = 0 ,..., M_N
   

       Sigma(k) = Exp(- Omega ((k-M_C)/(M_N-M_C))^Order)

   where
   

     MN    = Last resolvable mode by the methods.
     MC    = Cutoff mode < MN
     Omega = -ln(epsilon), epsilon=machine zero 
     Order = Order of filter, typically in the range between 4-20.
   

2. Second Order Lanczos Filter

3. Second Order Raised Cosine Filter

4. Eighth Order Sharpened Raised Cosine Filter

Caution : The combined effect of mapping and filtering has not yet been fully understood mathematically.