Clenshaw–Curtis quadrature

前情提要: 數值積分初探

前情提要: 以內插多項式來做數值積分

前情提要: 高斯積分

我們想要利用 Chebyshev polynomial 來算積分

Goal: 任意給定一可積分函數 $f(x)$, $x\in[-1, 1]$, 我們想要算 $\int^1_{-1} f(x) dx$.

Idea:

我們簡單做個變數變換 $x = \cos\theta$ 可以得到 $$ \int^1_{-1} f(x) dx = \int^{\pi}_0 f(\cos\theta)\sin\theta d\theta. $$ 所以如果我能將 $f(x)$ 這個函數以 first kind Chebyshev polynomial 來做展開, i.e., $$ f(x) = \sum^N_{k=0} a_k T_k(x) = \sum^N_{k=0} a_k \cos(k\cos^{-1}x). $$ 那代入後我們就有 $$ \int^1_{-1} f(x) dx = \int^{\pi}_0 \sum^N_{k=0} a_k \cos(k\theta)\sin\theta d\theta = \sum^{N/2}_{k=0}\frac{2 a_{2k}}{1-(2k)^2}. $$

因此, 若我們知道如何求一個函數的 Chebyshev coefficients, $a_k$, 那我們很容易就可以計算 $\int^1_{-1}f(x)dx$.

簡單的 Matlab 程式如下:

% generate grid points
N = 25;
x = -sin((-N:2:N)'*pi/(2*N));

% evaluate at grid points
f = cos(x).*exp(sin(x));

% extend the function and take scaled FFT
f_extend = [f; f(N: -1:2)];
f_hat = fft(f_extend)/(2*N);

% obtain the Chebyshev coefficients
a = zeros(N+1,1);
a(1) = real(f_hat(1));
a(2:N) = 2*real(f_hat(2:N));
a(N+1) = real(f_hat(N+1));

% Clenshaw–Curtis weights
w = zeros(1,N+1);
if(mod(N,2)==0)
    M = N/2;
else
    M = (N-1)/2;
end
w(1:2:end) = 2./(1-4*(0:M).^2);

% Clenshaw–Curtis quadrature
cc_quad = w*a;

% error
exact_quad = exp(sin(1)) - exp(sin(-1));
abs(cc_quad - exact_quad)

以上做法是先將函數取值, 算其係數, 再求其積分. 不過實際上我們應該能將積分權重直接算出來, 如此一來我們就能直接由函數值求積分了. 也就是說我們希望能寫成以下這式子: $$ \int^1_{-1} f(x)dx \approx \sum^N_{k=0} w_k f(x_k), \quad x_k = \cos\left(\frac{k\pi}{N}\right). $$

積分權重:

要求得積分權重我們首先回顧一下上面的作法：

1. 首先將點做 even extension
2. 做 FFT
3. 由 Fourier coefficients 求得 Chebyshev coefficients
4. 再將係數乘上積分權重即得到積分值

以上步驟可以寫成矩陣形式: $$ q = [\pmb{w}][\pmb{I}]^T[\pmb{F}][\pmb{I}][\pmb{f}] $$ 其中

• $[\pmb{f}]$ 是函數值, 是個 $(N+1)\times 1$ 的矩陣
• $[\pmb{I}]$ 是 even extension 是個 $2N\times (N+1)$ 的矩陣
• $[\pmb{F}]$ 是 Fourier transform 是個 $2N\times 2N$ 的矩陣
• $[\pmb{w}]$ 是 Clenshaw–Curtis 對於係數的 weights 是個 $1\times (N+1)$ 的矩陣
• $q$ 是最後得到的數值積分值

所以這樣看就很清楚, Clenshaw–Curtis 對於函數值的 weights 事實上就是 $[\pmb{w}][\pmb{I}]^T[\pmb{F}][\pmb{I}]$ 這串矩陣乘法. 也就是說如果我定義 $$ [CCW]=[\pmb{w}][\pmb{I}]^T[\pmb{F}][\pmb{I}], $$ 則我們有 $$ q = [CCW][\pmb{f}] $$ 而由於 $[\pmb{F}]$ 是對稱矩陣, 我們可以把這整個做矩陣轉置, 這樣就很清楚要怎樣求 weight 了 $$ [CCW]^T = [\pmb{I}]^T[\pmb{F}]^T[\pmb{I}][\pmb{w}]^T, $$ 其中 $[\pmb{F}]^T$ 這個矩陣乘法則是以 iFFT 來加速計算.

簡單的 Matlab 程式如下:

N=25;

w = zeros(1,N+1);
if(mod(N,2)==0)
    M=N/2;
else
    M = (N-1)/2;
end
w(1:2:end) = 2./(1-4*(0:M).^2);

w2 = [w, w(N: -1:2)];
CCW = ifft(w2');
CCW = [CCW(1), 2*CCW(2:N)', CCW(N+1)];

有了權重後即可由函數值直接求得積分值:

% generate grid points
N = 25;
x = -sin((-N:2:N)'*pi/(2*N));

% evaluate at grid points
f = cos(x).*exp(sin(x));

% Clenshaw–Curtis quadrature
cc_quad2 = CCW*f;

% error
exact_quad = exp(sin(1)) - exp(sin(-1));
abs(cc_quad2 - exact_quad)

最後附上 python 程式:

import numpy as np
from numpy.fft import fft,ifft,fftshift

N = 25

# Clenshaw–Curtis weights 
if N % 2 ==0:
    M= int(N/2)
else:
    M = int((N-1)/2)

w = np.zeros(N+1)

w[0:N+1:2] = 2./(1-4*np.linspace(0,M,M+1)**2);
w2 = np.append(w, w[N-1:0:-1])
CCW2 = np.real(ifft(w2))

CCW = 2*CCW2[0:N+1]
CCW[0] = CCW2[0]
CCW[N] = CCW2[N]

# Chebyshev points of the second kind
x = -np.sin(np.arange(-N, N+1, 2)*np.pi/(2*N))

#f = np.exp(np.sin(x));
f = np.cos(x)*np.exp(np.sin(x))
cc_quad2 = np.inner(CCW, f)

# error
exact_quad = np.exp(np.sin(1)) - np.exp(np.sin(-1));
np.abs(cc_quad2 - exact_quad)