The Ultraspherical Spectral Method | ALL YOU NEED IS LOVE

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:

Representing derivatives in terms of ultraspherical polynomials. This results in diagonal differentiation matrices.
Representing conversion from Chebyshev polynomials to ultraspherical polynomials by a banded operator.
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.
Imposing the boundary conditions using boundary bordering, that is, rows of linear system are used to impose boundary conditions.
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,

Numerical Example

Now we consider the following first order differential equation

with boundary condition , where . The exact solution is given by

We can use the ultraspherical spectral method to solve this problem.

Code: An implementation of the ultraspherical spectral method for the first order differential equation.

% Solve the First order BVP u' + a*u = f on [-1,1]
% with boundary condition u(-1) = c by Ultraspherical Spectral Method

%% Problem setting
% Solving u'(x) + 1/(\alpha x^2+1) u(x) = 0
% u(x) = \exp(-(\tan^{-1}(\sqrt{\alpha}x) + \tan^{-1}(\sqrt{\alpha}))/\sqrt{\alpha})

alpha = 5*10^4;
a = @(x) 1./(alpha*x.^2+1);             % variable coefficient term
f = @(x) 0*x;                           % source term
c = 1;

n = 1000;                               % use n+1 CGL points

%% Step 1. Approx a and f by Chebyshev interpolation
% Compute the Chebyshev-Gauss-Lobatto points
j = (0:n)';
x = -cos(pi*j/n);                                   % n+1 CGL points
% w = pi/n * [1/2; ones(n-1,1); 1/2];               % CGL Weight

% function values at CGL points
av = a(x);
fv = f(x);

% Forward Chebyshev Transform via FFT
ahat = sparse(chebtransf(av));                      % Chebyshev coefficients of a
fhat = chebtransf(fv);                              % Chebyshev coefficients of f

index = find(abs(ahat) < 2^(-52));
ahat(index) = 0;

%% Step 2. Construction Differential, Multiplication and Conversion Matrix
% First Order Differential Matrix
D = zeros(n+1);
D(1:end-1, 2:end) = D(1:end-1, 2:end) + diag(1:n);

% Multiplication Matrix with respect to function a
T = toeplitz([2*ahat(1); ahat(2:end)]);                 % Toeplitz matrix
H = [zeros(1, n+1); [hankel(ahat(2:end)) zeros(n,1)]];  % zero-plus-hankel
M = 0.5 * (T + H);                                      % Multiplication Matrix

% Conversion Matrix, transform Chebyshev coefficients to C^1 coefficients
S = sparse(eye(n+1) / 2);
S(1, 1) = 1;
S(1:end-2, 3:end) = S(1:end-2, 3:end) - eye(n-1) / 2;

%% Step 3. Boundary Bording
FirstRow = (-1) .^ (0:n);
A = [FirstRow; D+S*M];
F = [c; S*fhat];

%% Step 3. Solve linear system: A*uhat = F
uhat = A \ F;                                           % To be improved by QR decomposition

%% Step 4. Backward Chebyshev Transform
uv = ichebtransf(uhat);                                 % approximant values on CGL points

%% MAG
plot(x,uv,LineStyle=":", LineWidth=2, Marker="o")
hold on

t = -1:0.0001:1;
u = exp(-1/sqrt(alpha)* ( atan(sqrt(alpha)*t) + atan(sqrt(alpha)) ));
plot(t, u, LineWidth=2)
grid on
legend("Ultraspherical Approximant with order " + num2str(n+1), ...
    "Exact solution $\exp(-(\tan^{-1}(\sqrt{\alpha}x) + \tan^{-1}(\sqrt{\alpha}))/\sqrt{\alpha})$", ...
    Interpreter="latex")
hold off

%% Fast Chebyshev Transform
% CHEBTRANSF Chebyshev Transformation via FFT.
% Given the value of f on CGL points, use FFT to compute the Chebyshev coefficients 
% of f efficiently.
function fhat = chebtransf(fk)
    n = length(fk) - 1;
    % fk are the function values of f at CGL points 
    % -1 = x1 < x2 < ... < xn < x(n+1) = 1
    % filp the vector s.t. the values are on CGL in descending order
    fk = fk(end:-1:1);              
    % Now fk(1) = f(1=x(n+1)), fk(2) = f(xn), ..., fk(end) = f(-1=x1)
    gk = [fk; fk(end-1:-1:2)];      % g(theta) = f(cos(theta)), even and periodic
    
    fhat = ifft(gk);                % Discrete Chebyshev Transform
    
    fhat = fhat(1:n+1);             % delete repeat elements
    fhat(2:end) = fhat(2:end) * 2;  % scaling
end
% ICHEBTRANSF Inverse Chebyshev Transformation via FFT.
% Given the Chebyshev coefficients of a function f, use FFT to compute the
% function values of f at CGL efficiently.
function fk = ichebtransf(fhat)
    % fhat are the Chebyshev coefficients of f obtained by interpolating at CGL points.
    n = length(fhat) - 1;
    
    fhat(2:end) = fhat(2:end) / 2;  % scaling
    ghat = [fhat; fhat(end-1:-1:2)];
    
    fk = fft(ghat);                 % Inverse Chebyshev Transform
    
    fk = fk(1:n+1);                 % delete repeat elements
    fk = real(fk);                  % delete imaginary part
    fk = fk(end:-1:1);              % filp back
end

The above code solves the first order differential equation with variable coefficient and boundary condition. The result is shown in the figure below.

The numerical solution of the first order differential equation

by ultraspherical spectral method and the exact solution.

And the sparsity of the matrix is shown below.

The sparsity of the matrix

given by the ultraspherical spectral method (in the case of

and

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:

Dirichlet boundary conditions:
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 .

Differentiation operator

Now suppose is represented by a Chebyshev series

then differentiating with respect to for times we have

Using the properties of ultraspherical polynomials above we can express the differentiation operator in terms of ultraspherical polynomials:

Therefore, the -order differentiation operator can be expressed as

where the first columns are all zeros.

Multiplication operator

In order to handle the variable coefficients in the differential operator, we need to represent the multiplication of two ultraspherical series in coefficient space. Given two ultraspherical series of same order

directly multiplying them gives

To compute the coefficients , we use the linearization formula by Carlitz

where are the coefficients of the linearization formula, given by

To avoid the overflow of the factorials, we can use the following recurrence relation to compute :

Hence can be deduced to

and

As a result, the multiplication operator is given by

for . In practice, will be approximated by a truncated Chebyshev series, i.e. truncation of its series

by applying times conversion operator of order in next subsection, and with this approximation the matrix is banded with bandwidth for .

Recurrence Construction

Another more practical way to compute the multiplication operator is to use the recurrence relation of ultraspherical polynomials. Now we write the product of two ultraspherical polynomials in terms of quasimatrics as follows:

where is the vector of coefficients of in series. In this method, we define the multiplication operator as the operator such that

the th column of is the coefficient vector of :

The key observation is that the operator is linear: if we represent in basis

then we have

Since ultraspherical polynomials satisfy the following recurrence relation

we have

Therefore to construct , we only need the two former matrix , and an additional easily computed matrix :

by using

and recall that , . It’s trivial to verify that

Hence we can use and the initial condition

to construct by .

Conversion operator

To construct the final equation in coefficient space, we need to convert the lower order ultraspherical coefficients to higher order ultraspherical coefficients in both sides.

The conversion operator is defined as the operator that transforms the coefficients in a series to coefficients in a series. By using the recurrence relation of ultraspherical polynomials of different orders, we have

Discretization

By using the operators defined in previous sections, we can discretize the -th order differential operator as

which takes coefficients in a Chebyshev series to those in a series. Therefore, to obtain the discrete system, the right hand side must be also expressed in terms of coefficients, which can be done by applying the conversion operator to the Chebyshev coefficients of .

Moreover, the boundary conditions can be imposed by replacing the last rows of the matrix with the boundary operator , and then permuting the rows to the top of the matrix.

Finally, we have the following linear system to solve for the Chebyshev coefficients of :

where

and is again the vector of Chebyshev coefficients of . In this end, the solution is approximated by the term Chebyshev series

Numerical Example

In this section we give two examples to illustrate the performance of the ultraspherical spectral method. The first example is a second order differential equation with variable coefficients, and the second example is ten order differential equation with variable coefficients. The numerical results and the sparsity of the matrix are shown, and the codes are provided.

A Second Order Example
First we consider the following variable coefficient second order differential equation

with boundary conditions , , set , and the solution is given by

where is the Airy function. The numerical solution with order and the exact solution are shown in the figure below.

The numerical solution of the second order differential equation

by ultraspherical spectral method and the exact solution.

And the sparsity of the matrix in the final discretization is shown below.

Code: An implementation of the ultraspherical spectral method for the second order example.

% Solve the Second order BVP a2*u'' + a1*u' + a0*u = f on [-1,1]
% with boundary condition [u(-1);u(1)] = c by Ultraspherical Spectral Method

%% Problem setting
% Solving Airy equation eps*u''(x) - x*u(x) = 0 with boudary condition
% u(-1) = airy(-eps^(-1/3)), u(1) = airy(eps^(-1/3))
% Exact solution: u(x) = airy(eps^(-1/3) * x)

eps = 1e-6;                             % parameter
N = 2;                                  % order of the equation

a2 = @(x) 0*x + eps;
a1 = @(x) 0*x;
a0 = @(x) -x;                           % variable coefficient term
varcoefs = {a0, a1, a2};

f = @(x) 0*x;                           % source term

% Boundary condition
c = [airy(-eps^(-1/3)); airy(eps^(-1/3))];
K = length(c);

n = 10000;                              % use n+1 CGL points

%% Step 1. Approx a and f by Chebyshev interpolation
% Compute the Chebyshev-Gauss-Lobatto points
j = (0:n)';
x = -cos(pi*j/n);                                   % n+1 CGL points
% w = pi/n * [1/2; ones(n-1,1); 1/2];               % CGL Weight 

% function values at CGL points
av = zeros(n+1, N+1);
for i = 1:N+1
    av(:, i) = varcoefs{i}(x);
end
fv = f(x);

% Forward Chebyshev Transform via FFT
% Chebyshev coefficients of a's
tol = 2^(-48);
ahat = zeros(n+1, N+1);
alength = zeros(N+1);
for k = 1:N+1
    ahat(:, k) = chebtransf(av(:, k));
    index = find(abs(ahat(:, k)) < tol);
    ahat(index, k) = 0;
    alength(k) = n+1 - length(index);
end

fhat = chebtransf(fv);                              % Chebyshev coefficients of f
%% Step 2. Construction Differential, Multiplication and Conversion Matrix
D = cell(1, N);
for k = 1:N
    D{k} = 2^(k-1)*factorial(k-1)* spdiags((0:n)', k, n+1, n+1);                  % kth order differential matrix
end

% Conversion Matrix, transform Chebyshev coefficients to C^1 coefficients
S = cell(1, N);
b0 = [1; ones(n, 1)/2];
b2 = -ones(n+1, 1)/2;
d = [0; 2];
S{1} = spdiags([b0, b2], d, n+1, n+1);
% S{1} = diag([1; ones(n, 1)/2]) - diag(ones(n-1, 1)/2, 2);

for k = 2:N
    b0 = [1, (k-1) ./ (k:k+n-1)]';
    b2 = [zeros(1,2), -(k-1) ./ (k+1:k+n-1)]';
    S{k} = spdiags([b0, b2], d, n+1, n+1);
    % S{k} = diag([1, (k-1) ./ (k:k+n-1)]) - diag((k-1) ./ (k+1:k+n-1), 2);
end

% Multiplication Matrix with respect to function a
M = cell(1, N+1);       % N order equation, N+1 terms

% first we transform the coefficients of ahat to appropriate order
for lambda = 1:N
    for k = 1:lambda
        ahat(:, lambda+1) = S{k} * ahat(:, lambda+1);
    end
end

% 0 order term: M(:, :, 1)
T = toeplitz([2*ahat(1, 1); ahat(2:end, 1)]);               % Toeplitz matrix
H = [zeros(1, n+1); [hankel(ahat(2:end, 1)) zeros(n,1)]];   % zero-plus-hankel
M{1} = 0.5 * (T + H);                                 % Multiplication Matrix

for lambda = 1:N
    M0 = speye(n+1);    % M[C_0^lambda]
    
    % construct M[x]
    bn1 = [((0:n-1) + 1) ./ (2*(0:n-1) + 2*lambda), 0]';
    bp1 = [0, ((1:n)+2*lambda-1) ./ (2*(1:n) + 2*lambda)]';
    Mx = spdiags([bn1, bp1], [-1, 1], n+1, n+1);
    % Mx = diag( ((0:n-1) + 1) ./ (2*(0:n-1) + 2*lambda), -1) +...
    %      diag( ((1:n)+2*lambda-1) ./ (2*(1:n) + 2*lambda), 1);
    
    M1 = 2*lambda*Mx;   % M[C_1^lambda]

    % Construct next multiplication operator by three term recurrence
    M{lambda+1} = ahat(1, lambda+1) * M0 + ahat(2, lambda+1) * M1;
    for k = 2 : (alength(lambda)+1)

        tempM0 = M0;

        M0 = M1;
        M1 = 2*(k-1+lambda)/k * Mx * M1 - (k-1+2*lambda-1)/k * tempM0;

        M{lambda+1} = M{lambda+1} + ahat(k+1, lambda+1) * M1;
    end
end

%% Step 3. Discrete Matrix
% Left hand matrix
L = zeros(n+1);
for k = 1:N+1
    if k == 1
        Lk = M{1};
        for j = 1:N
            Lk = S{j}*Lk;
        end
        L = L + Lk;
    elseif k == N+1
        L = L + M{N+1}*D{N};
    else
        Lk = M{k};
        for j = k:N
            Lk = S{j}*Lk;
        end
        Lk = Lk * D{k-1};
        L = L + Lk;
    end
end
% Right hand vector
F = fhat;
for k = 1:N
    F = S{k}*F;
end

%% Step 4. Boundary Bording
B = [(-1).^(0:n); ones(1, n+1)];
L = [B; L(1:end-K, :)];

F = [c; F(1:end-K)];

%% Step 5. Solve linear system: L*uhat = F
% uhat = L \ F;
uhat = AlmostBandedSolver(L, F, K);

%% Step 6. Backward Chebyshev Transform
uv = ichebtransf(uhat);                                 % approximant values on CGL points

%% MAG
figure(1)
plot(x,uv,LineStyle=":", LineWidth=2, Marker="o")
hold on

t = -1:0.0001:1;
u = airy(eps^(-1/3) * t);
plot(t, u, LineWidth=2)
grid on
legend("Ultraspherical Approximant with order n=" + num2str(n), ...
    "Exact solution $u(x) = \mathrm{Airy}(\epsilon^{-1/3} x), \epsilon=$" + num2str(eps), ...
    Interpreter="latex")
title("\bfUltraspherical Spectral Method")
xlabel("\boldmath{$x$}", "Interpreter","latex")
ylabel("\boldmath{$u_n$}", "Interpreter","latex")
hold off

figure(2)
spy(L, 'k')
title("\bfSparsity of Linear System")
xlabel("\boldmath{$j$}", "Interpreter","latex")
ylabel("\boldmath{$i$}", "Interpreter","latex")


%% Fast Chebyshev Transform
% CHEBTRANSF Chebyshev Transformation via FFT.
% Given the value of f on CGL points, use FFT to compute the Chebyshev coefficients 
% of f efficiently.
function fhat = chebtransf(fk)
    n = length(fk) - 1;
    % fk are the function values of f at CGL points 
    % -1 = x1 < x2 < ... < xn < x(n+1) = 1
    % filp the vector s.t. the values are on CGL in descending order
    fk = fk(end:-1:1);              
    % Now fk(1) = f(1=x(n+1)), fk(2) = f(xn), ..., fk(end) = f(-1=x1)
    gk = [fk; fk(end-1:-1:2)];      % g(theta) = f(cos(theta)), even and periodic
    
    fhat = ifft(gk);                % Discrete Chebyshev Transform
    
    fhat = fhat(1:n+1);             % delete repeat elements
    fhat(2:end) = fhat(2:end) * 2;  % scaling
end
% ICHEBTRANSF Inverse Chebyshev Transformation via FFT.
% Given the Chebyshev coefficients of a function f, use FFT to compute the
% function values of f at CGL efficiently.
function fk = ichebtransf(fhat)
    % fhat are the Chebyshev coefficients of f obtained by interpolating at CGL points.
    n = length(fhat) - 1;
    
    fhat(2:end) = fhat(2:end) / 2;  % scaling
    ghat = [fhat; fhat(end-1:-1:2)];
    
    fk = fft(ghat);                 % Inverse Chebyshev Transform
    
    fk = fk(1:n+1);                 % delete repeat elements
    fk = real(fk);                  % delete imaginary part
    fk = fk(end:-1:1);              % filp back
end

A Ten Order Example
The second example is a ten order differential equation with variable coefficients

with boundary conditions , and . The real solution is hard to obtain, but we can still use the ultraspherical spectral method to solve it numerically. The numerical solution with order are shown in the figure below.

The numerical solution of the ten order differential equation by ultraspherical spectral method.

And the sparsity of the matrix in the final discretization is shown below.

The code for this example is similar to the previous one, and the only difference is that we need to change the variable coefficients and the boundary conditions. The code is provided below.

Code: An implementation of the ultraspherical spectral method for the ten order example.

% Solve the higher order BVP an*u^(n) + ... + a1*u' + a0*u = f on [-1,1]
% with some boundary condition \mathcal{B}u = c by Ultraspherical Spectral Method

%% Problem setting
% Solving u^(10) + cosh(x) u^(8) + x^2 u^(6) + x^4 u^(4) + cos(x) u^(2) +
% x^2 u = 0 with the following boundary condition:
% u(-1) = u(1) = 0;
% u'(-1) = u'(1) = 1;
% u^(i)(-1) = u^(i)(1) = 0 for i = 2,3,4.

N = 10;                                 % order of the equation

a10 = @(x) 0*x + 1;
a9 = @(x) 0*x;
a8 = @(x) cosh(x);
a7 = @(x) 0*x;
a6 = @(x) x.^2;
a5 = @(x) 0*x;
a4 = @(x) x.^4;
a3 = @(x) 0*x;
a2 = @(x) cos(x);
a1 = @(x) 0*x;
a0 = @(x) x.^2;
% variable coefficient term
varcoefs = {a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10};

f = @(x) 0*x;                           % source term

% Boundary condition
c = [zeros(2, 1); ones(2, 1); zeros(6, 1)];
K = length(c);

n = 100;                               % use n+1 CGL points

%% Step 1. Approx a and f by Chebyshev interpolation
% Compute the Chebyshev-Gauss-Lobatto points
j = (0:n)';
x = -cos(pi*j/n);                                   % n+1 CGL points
% w = pi/n * [1/2; ones(n-1,1); 1/2];               % CGL Weight 

% function values at CGL points
av = zeros(n+1, N+1);
for i = 1:N+1
    av(:, i) = varcoefs{i}(x);
end
fv = f(x);

% Forward Chebyshev Transform via FFT
% Chebyshev coefficients of a's
tol = 2^(-48);
ahat = zeros(n+1, N+1);
alength = zeros(N+1);
for k = 1:N+1
    ahat(:, k) = chebtransf(av(:, k));
    index = find(abs(ahat(:, k)) < tol);
    ahat(index, k) = 0;
    alength(k) = n+1 - length(index);
end

fhat = chebtransf(fv);                              % Chebyshev coefficients of f
%% Step 2. Construction Differential, Multiplication and Conversion Matrix
D = cell(1, N);
for k = 1:N
    D{k} = 2^(k-1)*factorial(k-1)* spdiags((0:n)', k, n+1, n+1);                  % kth order differential matrix
end

% Conversion Matrix, transform Chebyshev coefficients to C^1 coefficients
S = cell(1, N);
b0 = [1; ones(n, 1)/2];
b2 = -ones(n+1, 1)/2;
d = [0; 2];
S{1} = spdiags([b0, b2], d, n+1, n+1);
% S{1} = diag([1; ones(n, 1)/2]) - diag(ones(n-1, 1)/2, 2);

for k = 2:N
    b0 = [1, (k-1) ./ (k:k+n-1)]';
    b2 = [zeros(1,2), -(k-1) ./ (k+1:k+n-1)]';
    S{k} = spdiags([b0, b2], d, n+1, n+1);
    % S{k} = diag([1, (k-1) ./ (k:k+n-1)]) - diag((k-1) ./ (k+1:k+n-1), 2);
end

% Multiplication Matrix with respect to function a
M = cell(1, N+1);       % N order equation, N+1 terms

for lambda = 1:N
    for k = 1:lambda
        ahat(:, lambda+1) = S{k} * ahat(:, lambda+1);
    end
end

% 0 order term: M(:, :, 1)
T = toeplitz([2*ahat(1, 1); ahat(2:end, 1)]);               % Toeplitz matrix
H = [zeros(1, n+1); [hankel(ahat(2:end, 1)) zeros(n,1)]];   % zero-plus-hankel
M{1} = 0.5 * (T + H);                                 % Multiplication Matrix

for lambda = 1:N
    M0 = speye(n+1);
    
    bn1 = [((0:n-1) + 1) ./ (2*(0:n-1) + 2*lambda), 0]';
    bp1 = [0, ((1:n)+2*lambda-1) ./ (2*(1:n) + 2*lambda)]';
    Mx = spdiags([bn1, bp1], [-1, 1], n+1, n+1);
    % Mx = diag( ((0:n-1) + 1) ./ (2*(0:n-1) + 2*lambda), -1) +...
    %      diag( ((1:n)+2*lambda-1) ./ (2*(1:n) + 2*lambda), 1);
    
    M1 = 2*lambda*Mx;

    M{lambda+1} = ahat(1, lambda+1) * M0 + ahat(2, lambda+1) * M1;
    for k = 2 : (alength(lambda)+1)

        tempM0 = M0;

        M0 = M1;
        M1 = 2*(k-1+lambda)/k * Mx * M1 - (k-1+2*lambda-1)/k * tempM0;

        M{lambda+1} = M{lambda+1} + ahat(k+1, lambda+1) * M1;
    end
end

%% Step 3. Discrete System
% Left hand matrix
L = zeros(n+1);
for k = 1:N+1
    if k == 1
        Lk = M{1};
        for j = 1:N
            Lk = S{j}*Lk;
        end
        L = L + Lk;
    elseif k == N+1
        L = L + M{N+1}*D{N};
    else
        Lk = M{k};
        for j = k:N
            Lk = S{j}*Lk;
        end
        Lk = Lk * D{k-1};
        L = L + Lk;
    end
end

% Right hand vector
F = fhat;
for k = 1:N
    F = S{k}*F;
end

%% Step 4. Boundary Bording
% use T_n^{(k)}(\pm 1) = (\pm 1)^(n+k) \prod_{r=0}^{k-1} (n-r)(n+r)/(2r+1)
B = zeros(K, n+1);
B(1:2, :) = [(-1).^(0:n); ones(1, n+1)];
for j = 3:K
    if mod(j, 2) == 1
        k = (j-1)/2;
        for nn = 1:n+1
            r = 0 : k-1;
            B(j, nn) = (-1)^(k+nn-1) * prod( (nn-1-r) .* (nn-1+r) ./ (2*r+1) );
        end
    else
        B(j, :) = abs(B(j-1, :));
    end
end

L = [B; L(1:end-K, :)];
F = [c; F(1:end-K)];

%% Step 5. Solve linear system: L*uhat = F
% uhat = L \ F;
uhat = AlmostBandedSolver(L, F, K);

%% Step 6. Backward Chebyshev Transform
uv = ichebtransf(uhat);                                 % approximant values on CGL points

%% MAG
figure(1)
plot(x,uv,LineStyle=":", LineWidth=2, Marker="o")
grid on
legend("Ultraspherical Approximant with order n=" + num2str(n), ...
    Interpreter="latex")
title("\bfUltraspherical Spectral Method")
xlabel("\boldmath{$x$}", "Interpreter","latex")
ylabel("\boldmath{$u_n$}", "Interpreter","latex")

figure(2)
spy(L, 'k')
title("\bfSparsity of Linear System")
xlabel("\boldmath{$j$}", "Interpreter","latex")
ylabel("\boldmath{$i$}", "Interpreter","latex")


%% Fast Chebyshev Transform
% CHEBTRANSF Chebyshev Transformation via FFT.
% Given the value of f on CGL points, use FFT to compute the Chebyshev coefficients 
% of f efficiently.
function fhat = chebtransf(fk)
    n = length(fk) - 1;
    % fk are the function values of f at CGL points 
    % -1 = x1 < x2 < ... < xn < x(n+1) = 1
    % filp the vector s.t. the values are on CGL in descending order
    fk = fk(end:-1:1);              
    % Now fk(1) = f(1=x(n+1)), fk(2) = f(xn), ..., fk(end) = f(-1=x1)
    gk = [fk; fk(end-1:-1:2)];      % g(theta) = f(cos(theta)), even and periodic
    
    fhat = ifft(gk);                % Discrete Chebyshev Transform
    
    fhat = fhat(1:n+1);             % delete repeat elements
    fhat(2:end) = fhat(2:end) * 2;  % scaling
end
% ICHEBTRANSF Inverse Chebyshev Transformation via FFT.
% Given the Chebyshev coefficients of a function f, use FFT to compute the
% function values of f at CGL efficiently.
function fk = ichebtransf(fhat)
    % fhat are the Chebyshev coefficients of f obtained by interpolating at CGL points.
    n = length(fhat) - 1;
    
    fhat(2:end) = fhat(2:end) / 2;  % scaling
    ghat = [fhat; fhat(end-1:-1:2)];
    
    fk = fft(ghat);                 % Inverse Chebyshev Transform
    
    fk = fk(1:n+1);                 % delete repeat elements
    fk = real(fk);                  % delete imaginary part
    fk = fk(end:-1:1);              % filp back
end

Efficient Almost Banded System Solver

Without any preconditioning grows proportionally with , which is significantly better than the typical growth of in the condition number for the standard tau and collocation methods. In fact this condition number can be bounded by using a trivial, diagonal preconditioner that scales the columns. The preconditioner used here is

where the first entries in brackets are all .

Now we aim to solve the almost banded system of equations

in which the matrix A is banded, with non-zero superdiagonals and non-zero subdiagonals (so that ), except for the first dense boundary rows. The following algorithm solves the system in operations and storage, which is based on computing a QR factorization using Given’s rotations. However, the resulting upper triangular part will be dense because of fill-in caused by the boundary rows. We will show that these dense rows can still be represented sparsely, and that the resulting upper triangular matrix can be solved in operations.

(Figure 5.1 in paper A Fast and Well-Conditioned Spectral Method by Olver & Townsend) Structure of operators during QR factorization. Left: Depicts the structure of the original differential operator. Right: Depicts a

filled-in matrix obtained after upper triangularizing the first

columns.

The key observation to see the low rank structure of the filled-in part is to note that only the first rows of the matrix are dense and the Givens rotations only affect two rows at a time, so that the filled-in elements are only dependent on the first rows as long as they are enough far away from the diagonal.

This algorithm are basically consisted of two steps:

Reduce the almost banded matrix to upper triangular form using Givens rotations, in this process, we need to store and update the linear combination coefficients of the filled-in rows in a -row matrix as elimination proceeds.
Solve the resulting upper triangular system by back substitution using the sparse representation of the filled-in rows.

Recall that the final upper triangular matrix will be -banded if the original matrix is -banded, in fact we can show that the elements outside of the th superdiagonal are all linear combinations of the first rows of the original matrix.

Lemma. Let is a filled-in matrix if, for , the th row of has the form

where is the th standard basis vector, is a -row matrix and is part of the -row matrix of the first rows of the original matrix. Furthermore, every row of has the form

The remaining rows have the form

In each step, we have the following relation for the low rank part of the filled-in rows:

where is the partly filled-in matrix and is the dense first rows of the original matrix. If at some stage the Givens rotation is applied to the th and th rows of the matrix, then we update the th and th columns of the matrix by

to preserve , i.e., holds.

Code: An implementation of the efficient almost banded system solver.

AlmostBandedSolver.m

% ALMOSTBANDEDSOLVER. An efficient solver for almost banded system by QR
% decompostion with Givens Rotation.
%-------------------------------------------------------------------------------
% Problem: solve linear system L*u = F
% Input:    L - almost banded matrix
%           F - right hand vector
%           K - number of dense rows
% Output:   u - solution vector
%-------------------------------------------------------------------------------
% Remark.
% Almost banded - Only first K rows of L are dense, L(K+1:end, :) is banded.
%-------------------------------------------------------------------------------
% Ref.  A Fast and Well-Conditioned Spectral Method, Sheehan Olver, Alex
%        Townsend. SIAM Rev. 55, 2013
% Note. https://chenx.space/2024/12/21/UltrasphericalSpectral/
function u = AlmostBandedSolver(L, F, K)
    [n, ~] = size(L);

    % Boundary Bording rows
    B = L(1:K, :);                          % first K rows are dense
    
    % check band length
    indexL = find(L(n, :) ~= 0);
    indexR = find(L(K+1:end, n) ~= 0);
    
    % (mL, mR)-banded
    mL = n - indexL(1);
    mR = n - indexR(1) - K;
    m = mL + mR + 1;
    
    %------------------ triu(L, m) = triu(b'*B, m) -----------------------%
    % initialize a vector to store the linear combination coefficients for
    % filled-in rows
    % kth column stores the combination coefficients for triu(L, m)'s kth row
    % in the beginning, triu(L, m) = triu((I_K; 0)* L(1:K, :), m)
    b = zeros(K, n);      
    b(1:K, 1:K) = eye(K); 
    %---------------------------------------------------------------------%
    
    %============================ QR decompostion ========================%
    for j = 1 : n-1
        for i = min(n, mL+j) : -1 : 1+j % eliminate all the lower banded part
            %------------ Construct Givens Rotation ----------------------%
            [c, s, nu] = rotator(L(i-1:i, j));
            Q = [c, -s; s, c];
            %-------------------------------------------------------------%
    
            %-------- Eliminate elements by Givens rotator ---------------%
            L(i-1:i, j) = [nu; 0];

            if i <= K + j
                L(i-1:i, j+1:end) = Q' * L(i-1:i, j+1:end);               % -----------------
            else                                                          %                 |
                L(i-1:i, j+1:min(j+m-1, n)) = Q' *L(i-1:i, j+1:min(j+m-1, n));%             |
                % using min(j+m-1, n) to replace 'end' since filled-in elements             |
                % can be computed by vector b and matrix B.                                 |
            end                                                           %                 |
                                                                          %                 |
            F(i-1:i) = Q'*F(i-1:i);                     % update right hand vector          |
            %-------------------------------------------------------------%                 |
                                                                          %                 |
            % ---------Update the combination coeficients ----------------%                 |
            % recall triu(L, m) = triu(b'*B, m)                                             |
            % (b'*B)(i-1:i, j+1:end) = Q' * (b'*B)(i-1:i, j+1:end)  <------------------------
            % hence b(:, i-1:i)' = Q' * b(:, i-1:i)' by associative law
            % update the linear combination coefficients as QR continuing
            b(:, i-1:i) = b(:, i-1:i) * Q; 
            %-------------------------------------------------------------%
        end
    end
    %======================================================================%
    % figure() spy(L, 'k')
   

    % Sparse Representation for the upper triangular matrix obtianed by QR decompostion
    % 
    %                               ******-----
    %                                *    * L2|
    %                                 *    *  |
    %                                  * L1 * |
    %           matrix L:               *    *|
    %                                    *    *
    %                                     *   *
    %                                 0    *  *
    %                                       * *
    %                                         *
    %
    % mL+mR+1 = m
    % L = tril(L, mL+mR) + triu(L, mL+mR+1)
    %   = L1             + L2
    % Notice: triu(b'*B, m+1) = triu(L, m+1), therefore
    %   L1 = tril(L, m-1);    % (0, m)-banded upper triangular matrix
    %   L2 = b'*B;            % strictly m-upper triangular matrix
    % since L2 is semi-separable, the back-subsitution can be done fast
    
    p = zeros(K, 1);
    u = zeros(n, 1);
    %======================= Back Substitution ============================%
    for k = n:-1:1
        if k == n                                                          % compute the last entry
            u(k) = F(k) / L(k,k);
        elseif k >= n-m                                                    % classical substitution for last (n-m:n) rows
            u(k) = (F(k) - L(k, k+1:n)*u(k+1:n)) / L(k,k);
        else  % 1 <= k <= n-m-1                                            % special substitution for filled-in part
            p = u(k+m+1)*B(:, k+m+1) + p;                                  % update p
            u(k) = (F(k) - L(k, k+1:k+m)*u(k+1:k+m) - b(:, k)'*p) / L(k,k);% rows of L2 are linear combinations of rows of B 
            %                           |                     |            %
            %------------------------ L1 part ------------- L2 part -------%
        end
    end
    %======================================================================%

    % u = L \ F;
end


% ROTATOR. construct Givens Rotator
% Input:    vector x = [x1; x2]
% Output:   parameters c, s, nu
% s.t.  Q = [c, -s; s, c]
%       Q'x = [nu; 0], nu = norm(x, 2)
% ----------------------------------------------------------------------
% ref: Fundamentals of Matrix Computations (3rd edition) by D.S. Watkins
%       sec 3.2, P195.
function [c, s, nu] = rotator(x)
    beta = norm(x, "inf");              % check whether x=0
    if beta < 2^(-52)                   % if x=0, no need to rotate x
        c = 1; s = 0;                   % Q=I
        nu = 0;
    else
        % ------------ x!=0 ------------%
        x(1) = x(1) / beta;             % normalization
        x(2) = x(2) / beta;             % scaling to aovid overflow or harmful underflow
        nu = norm(x);                   % Q is unitary, norm(Qx) = tau = norm(x)
        c = x(1) / nu;
        s = x(2) / nu;
        nu = beta * nu;                 % re-scraling
    end
end

The Matlab codes in this note are saved in the Matlab Drive folder.