高阶N下傅里叶级数近似不收敛反而出现不稳定问题

你的现象是两个因素共同作用的结果：

1. 数值积分的高频振荡误差 #

当N增大到1000时，被积函数是round(x)*cos(nx)或round(x)*sin(nx)，其中n=1000，这是高频振荡的分段常数函数。scipy.integrate.quad的默认积分精度无法准确捕捉这种高频振荡的积分值，导致计算出的傅里叶系数An、Bn出现严重误差。这些错误的系数叠加后，就会让近似曲线出现随机扭曲，而非正常收敛。

2. 吉布斯现象（次要因素） #

round(x)是分段不连续函数（在整数点有跳变），傅里叶级数对这类函数的逼近本身会存在吉布斯现象：在跳变点附近，即使N很大，也会始终存在固定幅度的过冲，但这只会在跳变点附近出现局部异常，不会导致全局的随机扭曲——你看到的怪异结果主要还是数值积分误差导致的。

针对数值积分的问题，有两种可靠的解决途径：

方法一：提高数值积分的精度 #

修改my_fourier_coef函数，给quad设置更严格的误差容忍参数，让它更精确地计算高频振荡函数的积分：

def my_fourier_coef(f, n):
    integrand_a = lambda x: f(x) * np.cos(n * x)
    # 设置更高的精度要求，降低积分误差
    an_integral, _ = quad(integrand_a, -np.pi, np.pi, epsabs=1e-12, epsrel=1e-12)
    An = an_integral / np.pi
    
    integrand_b = lambda x: f(x) * np.sin(n * x)
    bn_integral, _ = quad(integrand_b, -np.pi, np.pi, epsabs=1e-12, epsrel=1e-12)
    Bn = bn_integral / np.pi

    return [An, Bn]

这种方法能缓解高n时的积分误差，但对于极端大的N（如10^4），仍可能存在性能和精度瓶颈。

方法二：手动推导解析傅里叶系数 #

因为round(x)是分段常数函数，我们可以直接推导傅里叶系数的解析表达式，完全避免数值积分的误差。

在区间[-π, π]内，round(x)的分段区间和对应值明确，可拆分区间计算积分：

def analytic_fourier_coef(n):
    # 定义round(x)在[-π, π]的分段区间与对应值
    intervals = [
        (-np.pi, -3), -3,
        (-3, -2), -3,
        (-2, -1), -2,
        (-1, 0), -1,
        (0, 1), 0,
        (1, 2), 1,
        (2, 3), 2,
        (3, np.pi), 3
    ]
    an = 0.0
    bn = 0.0
    for i in range(0, len(intervals), 2):
        a, b = intervals[i]
        val = intervals[i+1]
        # 计算cos(nx)的积分
        cos_int = (np.sin(n*b) - np.sin(n*a)) / n if n !=0 else (b - a)
        an += val * cos_int
        # 计算sin(nx)的积分
        sin_int = (np.cos(n*a) - np.cos(n*b)) / n if n !=0 else 0.0
        bn += val * sin_int
    An = an / np.pi
    Bn = bn / np.pi
    return [An, Bn]

修改plot_results函数使用解析系数：

def plot_results(N):
    import matplotlib.pyplot as plt
    x = np.linspace(-np.pi, np.pi, 10000)
    # 计算A0系数
    [A0, _] = analytic_fourier_coef(0)
    y = A0 * np.ones(len(x)) / 2
    for n in range(1, N):
        [An, Bn] = analytic_fourier_coef(n)
        y += An * np.cos(n * x) + Bn * np.sin(n * x)
    # 生成目标函数值
    f_vals = np.array([round(xi) for xi in x])
    plt.figure(figsize=(10, 6))
    plt.plot(x, f_vals, label="analytic")
    plt.plot(x, y, label="approximate")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.grid()
    plt.legend()
    plt.title(f"{N}th Order Fourier Approximation")
    plt.show()

# 测试N=1000
plot_results(1000)

这种方法计算出的系数是精确的，即使N=1000，近似曲线也只会在跳变点附近出现吉布斯现象的过冲，不会出现全局的随机扭曲。

内容的提问来源于stack exchange，提问作者sai prabhav