The Ultraspherical Spectral Method

This note is based on the paper “A Fast and Well-Conditioned Spectral Method” by Sheehan Olver and Alex Townsend.

Consider the family of linear differential equations on :

where is an th order linear differential operator

and denotes boundary conditions, , and are sufficiently smooth functions on . Assume this problem is not singular; i.e., the leading variable coefficient does not vanish on the interval .

The Ultraspherical Spectral Method solves this problem by the following steps:

1. Representing derivatives in terms of ultraspherical polynomials. This results in diagonal differentiation matrices.
2. Representing conversion from Chebyshev polynomials to ultraspherical polynomials by a banded operator.
3. Representing multiplication for variable coefficients by banded operators in coefficient space. This is achieved by approximating the variable coefficients by a truncated Chebyshev (and thence ultraspherical) series.
4. Imposing the boundary conditions using boundary bordering, that is, rows of linear system are used to impose boundary conditions.
5. Using an adaptive QR decomposition to determine the optimal truncation denoted by . Deverlop a sparse representation for the resulting dense matrix, allowing for efficient back substitution.

First order differential equations

First, we begin by considering the first order differential equation

where and are continuous functions with bounded variation on . This continuity assumption ensures there is a unique continuously differentiable solution on , and the bounded variation assumption ensures that the solution can be represented by a uniformly convergent Chebyshev series: an exact representation of a continuous function with bounded variation on by a Chebyshev series is given by

where is the th Chebyshev polynomial of the first kind.

Differentiation Operator

The derivative of Chebyshev polynomials satisfy

where is the th Chebyshev polynomial of the second kind, also known as the ultraspherical polynomial of order 1.

Suppose is represented by a Chebyshev series

then differentiating scales the coefficients and change the basis to ultraspherical polynomials:

Denote as the (infinite) vector of Chebyshev coefficients of , then the (infinite) vector of first order ultraspherical coefficients of is given by , where

Note this differentiation operator is sparse, contrast to the classic differentiation matrix in spectral collocation methods.

Multiplication Operator

Given two Chebyshev series

the product can be represented as a new Chebyshev series:

where satisfies

As a result, we have

where is a Toeplitz matrix plus an almost Hankel operator given by

Though above matrix seems to be dense, however, since is continuous with bounded variation, we are able to uniformly approximate with a finite number of Chebyshev coefficients to any desired accuracy. If we truncate the series with terms, that is, set for , then the matrix becomes banded with bandwidth .

Conversion Operator

Now we need an operator that maps coefficients in a Chebyshev series to those in a series to cooperate the differentiation operator and the multiplication operator. Recall the Chebyshev polynomials can be written in terms of the polynomials by the recurrence relation

Hence

therefore, the conversion operator is given by

which ensures to be the vector of coefficients of .

Discretization of the system

The first order differential operator can be discretized as

and the discretized system (without its boundary condition) is given by

both sides of the equation are infinite vectors whose entries are the coefficients of the corresponding functions, and are the vectors of Chebyshev coefficients of and respectively.

The operators above are all infinite dimensional, we need to truncate them to derive a practical numerical scheme, which corresponds to applying the projection operator given by

Truncating the differentiation operator to results in an martix with a zero last row, this observation motivates us to impose the boundary condition by replacing the last row of . We take the convention of permuting this boundary row to the first row so that the linear system is close to upper triangular. That is, in order to obtain an approximate solution to the differential equation, we solve the following linear system

with

The solution is then approximated by the computed solution,

High order differential equations

Now we consider the general th order linear differential equation of the form

on with general boundary conditions . Assume that the boundary operator is given in terms of the Chebyshev coefficients of , for example, the operators corresponding to Dirichlet and Neumann boundary condition have the following forms:

1. Dirichlet boundary conditions:
   
2. Neumann boundary conditions:
   

To generalize the spectral method to high order differential equations we use similar relations, inculding differentiation, multiplication, and conversion, in terms of higher order ultraspherical polynomials.

The ultraspheical polynomials are a family of polynomials orthogonal with respect to the weight function on , where is a real parameter. Here we only need the ultraspherical polynomials for , defined uniquely by normalizing the leading coefficient so that

where is the Pochhammer symbol. In particular, the ultraspherical polynomials of order 1 are the Chebyshev polynomials of the second kind, that is, .

Ultraspherical polynomials satisfy the following important differentiation formula

where . They also have similar recurrence relation

provided .