Idea:

% 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);

% error


積分權重:

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

• $[\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$ 是最後得到的數值積分值

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));

% error


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))