通过切比雪夫极值点的快速离散傅里叶变换进行插值

Question

我试图在任意点插值函数，并且我在切比雪夫极点处有函数值。我使用快速傅里叶变换的实数值来计算切比雪夫系数。然后我用 2/N 缩放它们，然后使用多项式库来评估一组点处的切比雪夫多项式系列。这会产生错误的函数近似。我哪里错了？

import numpy as np
import matplotlib.pyplot as plt

# Define the number of 
# Chebyshev extreme points
N = 10

# Define the function to be 
# approximated
def f(x):
   return x**2

# Evaluate the function at the 
# Chebyshev extreme points
x = np.cos(np.arange(N) * np.pi / N)
y = f(x)

# Compute the discrete Fourier 
# transform (DFT) of the function 
# values using the FFT algorithm
DFT = np.fft.fft(y).real

# Compute the correct scaling 
# factor
scaling_factor = 2/N

# Scale the DFT coefficients by 
# the correct scaling factor
chebyshev_coefficients = scaling_factor * DFT

# Use Chebval to 
# evaluate the approximated 
# polynomial at a set of points
x_eval = np.linspace(-1, 1, 100)
y_approx =  np.polynomial.chebyshev.chebval(x_eval, chebyshev_coefficients[::-1])

# Plot the original function 
# and the approximated function
plt.plot(x, y, 'o', 
label='Original function')
plt.plot(x_eval, y_approx, '-', 
label='Approximated function')
plt.legend()
plt.show()

Answer 1

你犯了几个错误。

1. 对于 $i \in [0, N - 1]$，切比雪夫 $N$ 点的位置由 $cos(i\pi / (N - 1))$ 给出。
2. 要获取切比雪夫系数，您需要使用函数值的填充向量，即 np.hstack((y, y[N - 2: 0: -1]))
3. 正确的比例变为 2 / (2 * N - 2)
4. 缩放后第一个系数必须乘以0.5
5. 如果切比雪夫点的数量是奇数，那么最后一个系数也必须乘以0.5

这是代码的更新版本，运行良好。

import numpy as np
import matplotlib.pyplot as plt

# Define the number of 
# Chebyshev extreme points
N = 11

# Define the function to be 
# approximated
def f(x):
   return x**2

# Evaluate the function at the 
# Chebyshev extreme points
x = np.cos(np.arange(N) * np.pi / (N - 1))
y = f(x)

# Compute the discrete Fourier 
# transform (DFT) of the function 
# values using the FFT algorithm
DFT = np.fft.fft(np.hstack((y, y[N - 2: 0: -1]))).real / (2 * N - 2)

# Scale the DFT coefficients by 
# the correct scaling factor
chebyshev_coefficients = DFT[:N] * 2
chebyshev_coefficients[0] /= 2

if N % 2 != 0:
    chebyshev_coefficients[-1] /= 2    


# Use Chebval to 
# evaluate the approximated 
# polynomial at a set of points
x_eval = np.linspace(-1, 1, 100)
y_approx =  np.polynomial.chebyshev.chebval(x_eval, chebyshev_coefficients)

# Plot the original function 
# and the approximated function
plt.plot(x, y, 'o', 
label='Original function')
plt.plot(x_eval, y_approx, '-', 
label='Approximated function')
plt.legend()
plt.show()