Efficient Spectral-Galerkin Method

Gypsophila
  2025-01-13 20:38:26 2025-01-13 20:38:26 Created   2025-04-02 19:13:26 2025-04-02 19:13:26 Updated    

This note is based the series of seminal papers by Jie Shen:
1.Direct Solvers for the Second and Fourth Order Equations Using Legendre Polynomials;
2.Direct Solvers for the Second and Fourth Order Equations Using Chebyshev Polynomials;
3.Polar and Cylindrical Geometries;
4.Spherical Geometries.

The codes in this note can be found in MATLAB Online folder.

Second-Order Two-Point Boundary Value Problems

Consider solving the following second-order two-point boundary value ordinary differential equation:

in , where , with the general boundary conditions:

Without loss of generality, to ensure the well-posedness of the problem, we assume the following conditions to be satisfied:

  1. ;
  2. , ;
  3. , ;
  4. for all ;
  5. if ;
  6. if .

For convenience, we first reduce the problem to a problem with homogeneous boundary conditions by introducing a new simple function satisfying the boundary conditions:

  1. If and , set ;
  2. If , set ;

where and are constants to be determined by the boundary conditions. Then the original problem can be reduced to the following problem:

where , and the boundary conditions are homogeneous

To solve the problem by Galerkin method, we introduce the solution space

and construct the following standard weak formulation: find such that

where the bilinear form and the linear functional are defined as

and , where

To simplify the problem, we restrict ourselves to a special case where and , namely

in with the boundary conditions

The weak formulation of the problem is to find such that

where the bilinear form and the linear functional are defined as

and , where are defined as above.

Weight Galerkin Formulation

To use the Galerkin method, we need to specify the trial and test spaces (same in the case of the Galerkin method). A straightforward choice is the following polynomial space

however, this choice will lead to some shortcomings. To see this, we construct the weak formulation of the problem with a given weight function in the following way: find such that

There are two bad news about this formulation:

  1. The above formulation does not make sense if does not exist, except for the case of Dirichlet boundary conditions. Hence it can not be used for the Jacobi weight function with , including the Chebyshev weight functions.
  2. Even in the case of , the above formulation will not lead to a sparse or special linear system that can be solved efficiently.

A better choice is to use the following space

as the trial and test spaces by enforce the Robin boundary conditions exactly. Now the new weighted Galerkin formulation reads: find such that

where is usually taken to be the interpolation of associated with the nodes of the Gauss-type quadraturem points.

Given a set of basis function for , the solution can be represented as

and we denote

set , and

Then the weighted Galerkin formulation can be written as

Lagendre-Galerkin Method

The Lagendre-Galerkin method use the Legendre weight function and set to be the Legendre interpolation of at the Legendre-Gauss-Lobatto quadrature points. Now the weighted Galerkin formulation reads: find such that

Since the final linear system depends on the choice of the basis functions of , and just as in the finite element method, where neighboring points are used to form the basis functions so as to minimize their interations in the physical space, neighboring orthogonal polynomials should be used to form the basis functions in a spectral-Galerkin method so as to minimize their interactions in the frequency space. Now we look for basis functions as a compact combination of Legendre polynomials:

where are chosen to satisfy the boundary conditions. Such basis functions are called modal basis functions. The existence of such basis functions is guaranteed by the following lemma:

Lemma. For all , there exists a unique pair of real numbers and such that satisfies the boundary conditions

given by

where

Hint. To prove the lemma, it’s sufficient to notice and .

Especially, for Dirichlet and Neumann boundary conditions, we have

  • Dirichlet boundary conditions ( and ): ;
  • Neumann boundary conditions ( and ): .

Therefore, the stiffness matrix of the Lagendre-Galerkin method is a diagonal matrix with the diagonal entries

for . The mass matrix is symmetric and penta-diagonal with the entries

Hint. The diagonal entries of can be derived by using the orthogonality of Legendre polynomials

and the fact that , i.e.

The mass matrix can be obtained by the orthogonality of Legendre polynomials directly.

In short, the implementation of the Lagendre-Galerkin method is as follows:

  1. Pre-computation. Compute the Legendre-Gauss-Lobatto (LGL) quadrature nodes and weights , basis functions parameters , and the nonzero entries of and ;
  2. Forward Legendre Transform. Evaluate the Legendre coefficients of from function values at the LGL quadrature nodes and compute the right-hand side ;
  3. Solve Linear System. Solve the linear system for ;
  4. Backward Legendre Transform. Evaluate the function values at the LGL quadrature nodes from the Legendre coefficients .

Remark. Though Step 3 need only operations, the two discrete Legendre transforms in Steps 2 and 4 require operations. The overall complexity of the Lagendre-Galerkin method is .

Numerical Example

The following code is a simple implementation of the Lagendre-Galerkin method for solving the second-order two-point boundary value problem (with constant coefficients)

with the Robin boundary condition . The exact solution is .

The solver in code first computes the Legendre-Gauss-Lobatto quadrature nodes and weights, the basis functions parameters, and the stiffness and mass matrices. Then it evaluates the Legendre coefficients of the interpolation of the function at the LGL quadrature nodes and solves the linear system for . Finally, it evaluates the function values at the LGL quadrature nodes from the Legendre coefficients.

LegendreGalerkin.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
% Test of Legendre-Galerkin Method
% -u_xx + u = (1+pi^2)*sin(pi*x)-pi*x+pi with Robin boundary condition 
% u(-1)+u_x(-1)=0, u(1)=0
% Exact solution is u(x) = sin(pi*x)-pi*x+pi

%% Problem Setting
alpha = 1;                                  % -u'' + alpha*u = f

% Boundary coefficients
a1 = 1; b1 = 1;                             % a1*u(-1)+b1*u_x(-1)=0
a2 = 1; b2 = 0;                             % a2*u(+1)+b2*u_x(+1)=0

% Right hand function
f = @(x) (1+pi^2)*sin(pi*x)-pi*x+pi;

%% Error Plot and Solution Plot
%%%%%%%%%%%%%%%%%%%%%%%%%%% Error Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 30;                                     % use max N+1 points
E = zeros(N-3, 1);                          % Error

for n = 4:N
    [x,uv] = LegendreSolver(n,alpha,f,a1,a2,b1,b2); % Approximant
    u = sin(pi*x)-pi*x+pi;                  % Exact solution
    e = norm(u - uv, "inf");                % error
    E(n-3) = e;
end

figure(1)
semilogy(4:N,E, LineWidth=2, Marker="o")
grid on
yticks(10 .^ (-16:2:0))
title("\bfError of Legendre-Galerkin Method")
xlabel("$\textbf{Approx order}$ \boldmath{$n$}", "Interpreter","latex")
ylabel("\boldmath{$\log(|u_n-u|)$}", "Interpreter","latex")

%%%%%%%%%%%%%%% Approximant solution and exact solution %%%%%%%%%%%%%%%%%%%
[x,uv] = LegendreSolver(20,alpha,f,a1,a2,b1,b2);
u = sin(pi*x)-pi*x+pi;                      % Exact solution

figure(2)
plot(x, uv, LineWidth=2, Marker="o")
hold on
plot(x, u, LineWidth=1)
grid on
legend("Legendre-Galerkin Method with " + num2str(n+1) + " LGL points", ...
    "Exact solution $u(x) = \sin(\pi x)-\pi x+\pi$", ...
    Location="southwest", Interpreter="latex")
xlabel("\boldmath{$x$}", "Interpreter","latex")
ylabel("\boldmath{$u$}", "Interpreter","latex")
title("\boldmath{$-u'' + u = (1+\pi^2) \sin(\pi x)-\pi x+\pi$} " + ...
    "with \boldmath{$u(-1)+u'(-1)=0, u(1)=0$}", ...
    "Interpreter","LaTex")
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%% LEGENDRE SOVLER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,uv] = LegendreSolver(n,alpha,f,a1,a2,b1,b2)
%% Step 1. Precomputation
% Compute the Legendre-Gauss-Lobatto points via eigenvalue method
% See sec 3.3.2 in Spectral method by Shen (P98).
av = zeros(1,n-1); 
j  = (1:n-2)'; 
bv = j.*(j+2)./((2*j+1).*(2*j+3));
A  = diag(av)+diag(sqrt(bv),1)+diag(sqrt(bv),-1);   % form the Jacobi matrix
x  = sort(eig(sparse(A)));                          % find the eigenvalues
x  = [-1;x;1];                                      % n+1 LGL points

% LGL Weight
L = zeros(n+1);                                     % Legendre value matrix
L(1, :) = ones(1, n+1);                             % L0(x) = 1
L(2, :) = x;                                        % L1(x) = x
for j = 3:n+1
    % Lj(x) = (2j-1)/j x*L(j-1)(x) - (j-1)/j L(j-2)(x)
    L(j, :) = (2*j-3)/(j-1) * x' .* L(j-1, :) - (j-2)/(j-1) * L(j-2, :);  
end
w = 2/(n*(n+1)) ./ L(end, :).^2;

% Basis functions parameters
k     = (0:n-2)';
Delta = 2*a1*a2 + a1*b2*(k+2).^2 - a2*b1*(k+2).^2 - b1*b2/2*(k+1) .* (k+2).^2 .* (k+3);
a     = (a2*b1+a1*b2) * (2*k+3) ./ Delta;
b     = (-2*a1*a2 + (a2*b1-a1*b2)*(k+1).^2 + b1*b2/2* k.* (k+1).^2 .* (k+2)) ./ Delta;

% Stiff Matrix
sv = -(4*k+6) .* b;
S = diag(sv);

% Mass Matrix
m0 = 2 ./ (2*k+1) + 2*a.^2 ./ (2*k+3) + 2*b.^2 ./ (2*k+5);
m1 = 2*a(1:end-1) ./ (2*k(1:end-1)+3) + 2*a(2:end) .* b(1:end-1) ./ (2*k(1:end-1)+5);
m2 = 2*b(1:end-2) ./ (2*k(1:end-2) + 5);
M  = diag(m0) + diag(m1, 1) + diag(m1, -1) + diag(m2, 2) + diag(m2, -2);

%% Step 2. Forward Legendre Transform: fhat = G*L*W*fv
% See sec 3.3.3 in Spectral method by Shen (P100).
fv = f(x);                  % function value on LGL points

W = diag(w);                % LGL weight matrix

gamma = 2 ./ (2*(0:n)'+1);
gamma(end) = 2/n;           % special case of LGL points
G = diag(1 ./ gamma);       % Scaling matrix

fhat = G*(L*(W*fv));        % Forward Legendre Transform

% Basis coefficients: phik = Lk + ak*L(k+1) + bk*L(k+2)
gamma(end) = 2/(2*n+1);
F = gamma(1:end-2) .*     fhat(1:end-2) +...
    gamma(2:end-1) .* a.* fhat(2:end-1) +...
    gamma(3:end)   .* b.* fhat(3:end);

%% Step 3. Solve Linear System: (S+alpha*M)*U = F
U = (S+alpha*M) \ F;

%% Step 4. Backward Legendre Transform: uv = Phi*U
% Basis value matrix
Phi = L(1:end-2, :)' + L(2:end-1, :)' * diag(a) + L(3:end, :)' * diag(b);      

uv = Phi*U;                 % values of solution function at LGL points
end

Output of this code is the following error graph, in which we can see the error decreases in an exponential rate as the number of LGL points increases. Hence the Lagendre-Galerkin method is spectrally accurate.

Something is Wrong
Error of solving the BVP with the Robin boundary condition by the Lagendre-Galerkin method with Legendre-Gauss-Lobatto points.

By calling the function solver one additional time, we can also plot the approximate solution given by the Lagendre-Galerkin method with 21 Legendre-Gauss-Lobatto points and the exact solution . The graph is shown below.

Something is Wrong
The approximate solution given by the Lagendre-Galerkin method with 21 Legendre-Gauss-Lobatto points and the exact solution .

There are several things to note in the code:

  1. The Legendre-Gauss-Lobatto quadrature nodes and weights are computed by the eigenvalue method. In general, the interior nodes of Legendre-Gauss-type nodes are eigenvalues of the following Jacobian matrix:

    where
    • Legendre-Gauss: , , and ;
    • Legendre-Gauss-Radau: , , and ;
    • Legendre-Gauss-Lobatto: , , and .
  2. The forward Legendre transform is defined by

    where for , except for the LGL nodes case, where , and the weight is given by

    In matrix form, the forward Legendre transform can be written as
  3. The backward Legendre transform is defined by

    In matrix form, the backward Legendre transform can be written as

    where

    in which we use the MATLAB notation.
  4. The left hand of the linear equation is given by F in the code, which is compute by

Chebyshev-Galerkin Method

The main cost of the Lagendre-Galerkin method is the two discrete Legendre transforms in Steps 2 and 4, which require operations. By using the Chebyshev-Galerkin method, we can reduce the discrete transform cost to operations with the help of the Fast Fourier Transform (FFT), however, at the cost of the matrices and being no longer as symmetric and sparse as in the Lagendre-Galerkin method. Therefore, iterative methods are more suitable for solving the linear system in the Chebyshev-Galerkin method than direct methods.

Now we set and , which is the Chebyshev interpolant of at the Chebyshev-Gauss-Lobatto (CGL) quadrature points. The weak formulation of the BVP deduced to the Chebyshev-Galerkin method reads: find such that

which is equivalent to

Similar to the Lagendre-Galerkin method, the Chebyshev-Galerkin method uses neighboring Chebyshev polynomials to form the basis functions. The basis functions have the form

where is the Chebyshev polynomial of the first kind. And the basis functions parameters are chosen to satisfy the boundary conditions. The existence of such basis functions is guaranteed by the following lemma:

Lemma. For all , there exists a unique pair of real numbers and such that satisfies the boundary conditions

given by

where

Hint. To prove the lemma, it’s sufficient to notice and .

In the Chebyshev-Galerkin method, we define the solution space to be the space spanned by above basis functions

Recall and . The mass matrix of the Chebyshev-Galerkin method can be obtained by the orthogonality of Chebyshev polynomials directly:

where and for . Unfortunately, the stiffness matrix of the Chebyshev-Galerkin method is much harder to compute. We give the formula of the stiffness matrix when the boundary conditions are Dirichlet or Neumann:

  • Dirichlet boundary conditions ( and ): , and
  • Neumann boundary conditions ( and ): , and

In these two case, the stiffness matrix is an upper triangle matrix, and by the odd-even property of Chebyshev polynomials, we have for is odd. For general Robin boundary conditions, the stiffness matrix will become a upper full matrix and the computation of is more complicated. Since we seldomly use the Chebyshev-Galerkin method to deal with Problem with Robin boundary conditions (we may use the following Chebyshev-Legendre Galerkin method for example), we will not discuss the computation of in this case.

As a summary, the implementation of the Chebyshev-Galerkin method is as follows:

  1. Pre-computation. Compute the Chebyshev-Gauss-Lobatto (CGL) quadrature nodes and weights , basis functions parameters , and the nonzero entries of and ;
  2. Forward Chebyshev Transform. Evaluate the Chebyshev coefficients of from function values at the CGL quadrature nodes and compute the right-hand side ;
  3. Solve Linear System. Solve the linear system for ;
  4. Backward Chebyshev Transform. Evaluate the function values at the CGL quadrature nodes from the Chebyshev coefficients.

Remark. The Chebyshev-Galerkin method is more efficient than the Lagendre-Galerkin method in step 2 and 4, but the linear system in step 3 is more difficult to solve. The overall complexity of the Chebyshev-Galerkin method is still in general if we use direct methods to solve the linear system.

Numerical Example

The following code is a simple implementation of the Chebyshev-Galerkin method for solving the second-order two-point boundary value problem (with constant coefficients)

with the Dirichlet boundary condition . The exact solution is .

ChebyshevGalerkin.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
% Test of Chebyshev-Galerkin Method
% -u_xx + u = (1+pi^2)*sin(pi*x)+x^2-3 with Dirichlet boundary condition 
% u(-1)=0, u(1)=0
% Exact solution is u(x) = sin(pi*x)+x^2-1

%% Problem Setting
alpha = 1;                                  % -u'' + alpha*u = f

% Boundary coefficients
a1 = 1; b1 = 0;                             % a1*u(-1)+b1*u_x(-1)=0
a2 = 1; b2 = 0;                             % a2*u(+1)+b2*u_x(+1)=0

% Right hand function
f = @(x) (1+pi^2)*sin(pi*x)+x.^2-3;

%% Error Plot and Solution Plot
%%%%%%%%%%%%%%%%%%%%%%%%%%% Error Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 30;                                                 % use max N+1 points
E1 = zeros(N-3, 1); E2 = zeros(N-3, 1);                 % Error

for n = 4:N
    [x1,uv1] = ChebyshevSolver(n,alpha,f,a1,a2,b1,b2);  % Approximant by Chebyshev-Galerkin
    [x2,uv2] = LegendreSolver(n,alpha,f,a1,a2,b1,b2);   % Approximant by Chebyshev-Galerkin
    u1 = sin(pi*x1)+x1.^2-1;
    u2 = sin(pi*x2)+x2.^2-1;
    e1 = norm(u1 - uv1, "inf");                         % Chebyshev error
    e2 = norm(u2 - uv2, "inf");                         % Legendre error
    E1(n-3) = e1; E2(n-3) = e2;
end

figure(1)
semilogy(4:N,E1, LineWidth=2, Marker="o")
hold on
semilogy(4:N,E2, LineWidth=2, Marker="x")
legend("Chebyshev-Galerkin", "Legendre-Galerkin")
grid on
yticks(10 .^ (-16:2:0))
title("\bfError of Chebyshev/Legendre-Galerkin Method")
xlabel("$\textbf{Approx order}$ \boldmath{$n$}", "Interpreter","latex")
ylabel("\boldmath{$\log(|u_n-u|)$}", "Interpreter","latex")
hold off

%%%%%%%%%%%%%%% Approximant solution and exact solution %%%%%%%%%%%%%%%%%%%
n = 20;
[x1,uv1] = ChebyshevSolver(n,alpha,f,a1,a2,b1,b2);
[x2,uv2] = LegendreSolver(n,alpha,f,a1,a2,b1,b2);
t = -1:0.001:1;
u = sin(pi*t)+t.^2-1;                      % Exact solution

figure(2)
plot(x1, uv1, LineWidth=2, Marker="o", LineStyle=":")
hold on
plot(x2, uv2, LineWidth=2, Marker="x", LineStyle="--")
plot(t, u, LineWidth=1)
grid on
legend("Chebyshev-Galerkin Method with " + num2str(n+1) + " CGL points", ...
    "Legendre-Galerkin Method with " + num2str(n+1) + " LGL points", ...
    "Exact solution $u(x) = \sin(\pi x)+x^2-1$", ...
    Location="southeast", Interpreter="latex")
xlabel("\boldmath{$x$}", "Interpreter","latex")
ylabel("\boldmath{$u$}", "Interpreter","latex")
title("\boldmath{$-u'' + u = (1+\pi^2) \sin(\pi x)+x^2-3$} " + ...
    "with \boldmath{$u(-1)=0, u(1)=0$}", ...
    "Interpreter","LaTex")
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHEBYSHEV SOVLER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,uv] = ChebyshevSolver(n,alpha,f,a1,a2,b1,b2)
%% Step 1. Precomputation
% 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

% Basis functions parameters
k     = (0:n-2)';
Delta = 2*a1*a2 + ((k+1).^2 + (k+2).^2)*(a1*b2-a2*b1) - 2*b1*b2*(k+1).^2 .* (k+2).^2;
a     = 4*(a2*b1+a1*b2) * (k+1) ./ Delta;
b     = (-2*a1*a2 + (k.^2 + (k+1).^2)*(a2*b1-a1*b2) + 2*b1*b2* k.^2 .* (k+1).^2) ./ Delta;

% Stiff Matrix
if a1 == 1 && a2 == 1 && b1 == 0 && b1 == 0         % Dirichlet Boundary
    S = diag(2*pi*(k+1).*(k+2));
    t = 2;
    while t <= n-1
        S = S + diag(4*pi*(k(1:end-t)+1), t);
        t = t+2;
    end
elseif a1 == 0 && a2 == 0 && b1 == 1 && b2 == 1     % Neumann Boundary
    S = diag(2*pi*(k+1) .* k.^2 ./ (k+2));
    t = 2;
    while t <= n-1
        j = k + t;
        S = S + diag(4*pi* j(1:end-t).^2 .* (k(1:end-t)+1) ./ (k(1:end-t)+1).^2, t);
        t = t+2;
    end
else                                                % Robin Boundary
    fprintf("The algorithm for Robin BC has not been implemented yet.\n")
    return
end

% Mass Matrix
c  = [2; ones(n-2, 1)];
m0 = pi/2 * (c + a.^2 + b.^2);
m1 = pi/2 * (a(1:end-1) + a(2:end) .* b(1:end-1));
m2 = pi/2 * b(1:end-2);
M  = diag(m0) + diag(m1, 1) + diag(m1, -1) + diag(m2, 2) + diag(m2, -2);

%% Step 2. Forward Chebyshev Transform via iFFT
fv = f(x);                  % function value on LGL points

fhat = chebtransf(fv);      % Forward Chebyshev Transform

% Basis coefficients: phik = Lk + ak*L(k+1) + bk*L(k+2)
gamma = [2; ones(n, 1)] * pi/2;
F = gamma(1:end-2) .*     fhat(1:end-2) +...
    gamma(2:end-1) .* a.* fhat(2:end-1) +...
    gamma(3:end)   .* b.* fhat(3:end);

%% Step 3. Solve Linear System: (S+alpha*M)*U = F
U = (S+alpha*M) \ F;

%% Step 4. Backward Chebyshev Transform via FFT
uv = ichebtransf([U;0;0]) + ichebtransf([0; diag(a)*U; 0]) + ichebtransf([0;0;diag(b)*U]);
end

% 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 output of this code is the following error graph, in which we can see the error decreases in an exponential rate as the number of CGL points increases. Hence the Chebyshev-Galerkin method is spectrally accurate.

Something is Wrong
Error of solving the BVP with the Dirichlet boundary condition by the Chebyshev-Galerkin method with CGL points and the Legendre-Galerkin method with LGL points.

By calling the function ChebyshevSolver and LegendreSovler one additional time, we can also plot the approximate solution given by the Chebyshev/Lagendre-Galerkin method with 21 CGL/LGL points and the exact solution . The graph is shown below.

Something is Wrong
The approximate solution given by the Chebyshev/Lagendre-Galerkin method with 21 CGL/LGL points and the exact solution .

Remark. The theory and detailed implementation of the discrete Chebyshev transform via FFT can be found in Fast Chebyshev Transform.

  • Title: Efficient Spectral-Galerkin Method
  • Author: Gypsophila
  • Created at : 2025-01-13 20:38:26
  • Updated at : 2025-04-02 19:13:26
  • Link: https://chenx.space/2025/01/13/SpectralGalerkinByShen/
  • License: This work is licensed under CC BY-NC-SA 4.0.
Comments
On this page
Efficient Spectral-Galerkin Method
  1. Second-Order Two-Point Boundary Value Problems
    1. Weight Galerkin Formulation
    2. Lagendre-Galerkin Method
      1. Numerical Example
    3. Chebyshev-Galerkin Method
      1. Numerical Example