Clenshaw–Curtis & Gauss Quadrature

Clenshaw–Curtis Quadrature

给定一个在区间上定义的函数，Clenshaw–Curtis积分方法通过在Chebyshev点

上对进行插值，并利用插值得到的多项式的积分来近似原函数的积分。具体来说，设插值多项式为

其中是Chebyshev多项式，系数可以通过离散余弦变换（DCT）高效计算得到:

Clenshaw–Curtis积分的近似值由以下公式给出:

其中权重可以通过Chebyshev多项式的积分性质计算得到:

Code: MATLAB code for Clenshaw–Curtis quadrature via Chebfun.

clenshaw_curtis.m

function I = clenshaw_curtis(f,n)           % (n+1)-pt C-C quadrature of f 
x = cos(pi*(0:n)'/n);                       % Chebyshev points 
fx = feval(f,x)/(2*n);                      % f evaluated at these points 
g = real(fft(fx([1:n+1, n:-1:2])));         % Fast Fourier Transform 
a = [g(1); g(2:n)+g(2*n:-1:n+2); g(n+1)];   % Chebyshev coefficients 
w = 0*a'; w(1:2:end) = 2./(1-(0:2:n).^2);   % weight vector 
I = w*a;                                    % the integral

Integral Matrix

现在我们考虑计算以下不定积分：

首先与之前类似地在Chebyshev基下对进行近似

注意到当时Chebyshev多项式满足

且

于是

因此，积分可以表示为

其中

此处的为Chebyshev基下的积分矩阵，其Matlab实现见下方代码块。

Code: MATLAB code for Chebyshev based integral matrix.

cheb_intmat.m

function V = cheb_intmat(N)
% CHEB_INT_MATRIX 生成Chebyshev基下的谱积分矩阵
%
% V = cheb_int_matrix(N)
%
% 输入:
%   N : 输入多项式的最高阶数 (输入系数向量 a 的长度应为 N+1)
%
% 输出:
%   V : 大小为 (N+2) x (N+1) 的积分矩阵
%
% 用法:
%   如果 f(x) = sum(a_k * T_k(x))，令 a = [a_0; ...; a_N]
%   则 F(x) = int_{-1}^x f(t) dt 的系数 b 由 b = V * a 给出
%   其中 b 对应 F(x) = sum(b_k * T_k(x))，长度为 N+2

    % 初始化稀疏矩阵，大小为 (N+2) x (N+1)
    % 行索引对应 b_0 到 b_{N+1} (Matlab索引 1 到 N+2)
    % 列索引对应 a_0 到 a_N (Matlab索引 1 到 N+1)
    V = sparse(N+2, N+1);
    
    % 处理 k=0 (a_0) -> 贡献给 b_1 (T_1)
    % int T_0 = T_1
    V(2, 1) = 1;
    
    if N > 0
        % 处理 k=1 (a_1) -> 贡献给 b_2 (T_2)
        % int T_1 = 1/4 * T_2 (+ const term handled by BC later)
        V(3, 2) = 1/4;
    end
    
    % 处理 k >= 2 的通用情况
    % int T_k = 1/(2(k+1)) * T_{k+1} - 1/(2(k-1)) * T_{k-1}
    if N >= 2
        k = 2:N; % 输入的阶数 k
        
        % 对应的列索引 (a_k) 是 k + 1
        col_idx = k + 1;
        
        % 贡献给 T_{k+1} (即 b_{k+1}, 行索引 k+2)
        % 系数: 1 / (2*(k+1))
        row_idx_lower = k + 2;
        vals_lower = 1 ./ (2 * (k + 1));
        
        % 贡献给 T_{k-1} (即 b_{k-1}, 行索引 k)
        % 系数: -1 / (2*(k-1))
        row_idx_upper = k;
        vals_upper = -1 ./ (2 * (k - 1));
        
        % 填充矩阵
        ind_lower = sub2ind(size(V), row_idx_lower, col_idx);
        ind_upper = sub2ind(size(V), row_idx_upper, col_idx);
        
        V(ind_lower) = vals_lower;
        V(ind_upper) = vals_upper;
    end
    
    % 计算第一行 (b_0) 以满足 F(-1) = 0
    % F(-1) = b_0 - b_1 + b_2 - b_3 + ... = 0
    % b_0 = b_1 - b_2 + b_3 - ...
    % 也就是说，第一行是其余各行的交替和
    
    % 生成符号向量 [1, -1, 1, -1, ...] 对应行 2, 3, 4, ...
    row_indices = 2:(N+2);
    signs = (-1).^(row_indices); % b_1对应+1 (k=1, (-1)^2), b_2对应-1
    
    % 注意：b_k 的 k 从 1 开始对应矩阵的第 2 行。
    % b_1 (row 2): T_1(-1)=-1. Wait.
    % expansion: b_0*T_0(-1) + b_1*T_1(-1) + b_2*T_2(-1)... = 0
    % b_0(1) + b_1(-1) + b_2(1) + b_3(-1)... = 0
    % b_0 = b_1 - b_2 + b_3 - ...
    % Row 2 是 b_1，系数应为 +1
    % Row 3 是 b_2，系数应为 -1
    
    scale_factors = (-1).^(0:N); % [1, -1, 1, -1 ...] 长度为 N+1 (对应行2到end)
    
    % 将行向量 scale_factors 作用于 V 的剩余部分
    V(1, :) = scale_factors * V(2:end, :);
    
end

Gauss Quadrature

给定函数在个点处的插值信息，可以构造出唯一一个次数不超过的插值多项式作为的近似，再进行积分来作为的积分的近似。特别地，如果是一个次数不超过阶多项式，则在代数上精确成立，于是此类基于多项式的积分方法在使用个点时，能够至少精确积分所有不超过阶的多项式函数。进一步地，通过选择合适的插值节点，Gauss发现通过在一组特别选取的个节点—积分权重对应的正交多项式的零点—上进行插值，可以使得该积分方法精确积分所有不超过阶的多项式函数，即具有最高阶的代数精度。当固定使用个节点时，可以证明插值积分方法的代数精度最多为阶，因此Gauss积分方法是所有基于插值多项式的积分方法中代数精度最高的，进而在代数精度的意义下，Gauss积分方法是最优的。

对于积分权重的情形，Gauss积分方法对应的正交多项式为Legendre多项式，其节点称为Gauss-Legendre节点。具体地来讲，Gauss-Legendre积分实际上对应于首先将在Gauss-Legendre节点上进行插值，得到以下插值多项式

然后计算该插值多项式的积分来近似的积分，即

其中为Legendre多项式的第个零点，为对应的Lagrange基函数，权重由下式给出:

对于带有更一般的权重函数的积分

Gauss积分方法同样适用，只需将对应的正交多项式替换为与权重函数相关的正交多项式即可。例如，对于权重函数，其中，对应的正交多项式为Jacobi多项式，其节点称为Gauss-Jacobi节点。Gauss-Jacobi积分的形式与Gauss-Legendre积分类似:

其中为Jacobi多项式的第个零点，权重同样由对应Lagrange基函数的积分给出。

Gauss-Jacobi Quadrature Nodes and Weights

使用Gauss积分方法的关键在于如何高效且准确地计算出Gauss节点和权重。解决这一问题的经典方法是通过求解与正交多项式相关的三对角矩阵的特征值问题来获得节点和权重，即Golub-Welsch算法。

给定归一化的正交多项式序列，它们满足三项递推关系:

其中为递推系数，于是我们有

其中为大小为的Jacobi矩阵:

我们将上式简记为

对Jacobi矩阵进行特征值分解得

其中特征向量都是归一化的，即。基于以上分解，Golub和Welsch证明了以下事实：

1. Jacobi矩阵的特征值即为正交多项式的零点；
2. 权重满足
   
   其中为特征向量的第一个分量。

由于Jacobi矩阵为对称三对角矩阵，利用上这一结构可以在复杂度内完成特征分解，因此传统的Golub-Welsch算法可以在的计算成本下稳定地计算出Gauss节点和权重。

Code: MATLAB code for computing Gauss-Legendre quadrature by solving the tridiagonal eigenvalue problem, with a cost of .

gauss.m

function I = gauss(f,n)             % (n+1)-pt Gauss quadrature of f 
beta = .5./sqrt(1-(2*(1:n)).ˆ(-2)); % 3-term recurrence coeffs 
T = diag(beta,1) + diag(beta,-1);   % Jacobi matrix 
[V,D] = eig(T);                     % eigenvalue decomposition 
x = diag(D); [x,i] = sort(x);       % nodes (= Legendre points) 
w = 2*V(1,i).ˆ2;                    % weights 
I = w*feval(f,x);                   % the integral

事实上，近年来已经有多种方法被提出以进一步提升Gauss节点和权重的计算效率和精度。例如，Hale和Townsend提出了一种基于正交多项式的渐进表示以及Newton迭代的方法，可以在的复杂度下高效地计算Gauss-Legendre节点和权重[1]；Bogaert同样基于渐进展开进一步提出了一种无需迭代的直接计算方法，同样实现了的复杂度[2]。这些方法在处理大规模Gauss积分节点和权重的计算时在速度与精度上均远优于传统的Golub-Welsch算法。

References

1. Hale N, Townsend A. Fast and accurate computation of Gauss–Legendre and Gauss–Jacobi quadrature nodes and weights. SISC, 2013, 35(2): A652-A674.
2. Bogaert I. Iteration-free computation of Gauss–Legendre quadrature nodes and weights. SISC, 2014, 36(3): A1008-A1026.
3. Trefethen L N. Is gauss quadrature better than Clenshaw–Curtis?. SIREV, 2008, 50(1): 67-87.
4. Trefethen L N. Exactness of quadrature formulas. SIREV, 2022, 64(1): 132-150.