4.5 谱方法

4.5 谱方法

四种方法的精度天花板不同：差分与有限体积每加密一倍误差缩四倍，有限元高阶单元可以更高，但都停在代数收敛——误差按 n 的负幂下降。谱方法换一口气：用全局光滑基函数（正弦、切比雪夫多项式）展开解，只要真解足够光滑（解析），误差随 n 指数下降——每加一个基函数，精度多出一个数量级。它是 3.1 分离变量的计算化直系后裔，也是第一章傅里叶论战最远的回声。本节用周期热传导问题走通全套流程，并标出它的两条能力边界。

精度有没有天花板

谱方法的假设直接写在名字里：解本身用一组全局基函数表示。周期问题选傅里叶基 u(x) = Σ û_k·e^(ikx)，系数 û_k 由变换（离散情形是快速傅里叶变换，1965 年的库利–图基算法让每步变换代价从 n² 降到 n·log n，这是谱方法实用化的前提）算出。方法的核心动作只有一个：在系数空间里做微积分。对 e^(ikx) 求导恰是乘 ik：

∂x → ik，∂xx → −k²。

于是热传导方程 u_t = a²·u_xx 在系数空间里解耦成一族独立常微分方程：

dû_k/dt = −a²k²·û_k。

每个模式各自指数衰减（3.1 的特征值三件套在此重现：特征值 a²k²、特征函数 e^(ikx)），解为 û_k(t) = û_k(0)·e^(−a²k²t)。整个求解流程三步：变换到系数空间、每个系数独立演化、变换回物理空间。没有网格差商、没有相邻点耦合——空间离散被代数恒等式替换。

指数收敛的来源也在这里：差分方法的误差来自"差商近似导数"的局部近似，而谱方法的"求导"对基函数是精确的，唯一误差是"用 n 个基函数截断无穷级数"的截断误差。解析函数的傅里叶系数按指数速度衰减（复分析里的柯西不等式给出衰减率），截断误差自然指数下降。函数越光滑，谱方法越占便宜；解只有有限光滑度时，优势骤减到代数收敛——光滑度是它的全部本钱。

数值实验：与二阶差分同台对表

周期区间上比较谱方法与二阶中心差分解热传导问题（初始为一个光滑波包，与第三章同款初值族）。对每个分辨率 n 报告两者在 t = 0.1 的最大误差，见证指数与二阶的分野。

import numpy as np

L, a, T_end = 2*np.pi, 1.0, 0.1
x_full = lambda n: np.arange(n) * (2*np.pi/n)
u0 = lambda x: np.exp(-2*(np.sin(x/2))**2)        # 周期光滑初值

def spectral_solve(n):
    x = x_full(n)
    u = u0(x)
    k = np.fft.fftfreq(n, d=1.0/n)                 # 整数波数
    uh = np.fft.fft(u)
    uh_t = uh * np.exp(-a**2 * k**2 * T_end)       # 系数空间精确演化
    return np.fft.ifft(uh_t).real

def fd_solve(n):
    x = x_full(n); dx = 2*np.pi/n; dt = 0.4*dx**2
    u = u0(x)
    for _ in range(int(round(T_end/dt))):
        u[1:-1] = u[1:-1] + 0.4*(u[2:]-2*u[1:-1]+u[:-2])
        u[0], u[-1] = u[-2], u[1]                  # 周期边界
    return u

ref = spectral_solve(512)                          # 高分辨率参照解
for n in [8, 16, 32, 64]:
    es = np.abs(spectral_solve(n) - ref[:n]).max()
    ef = np.abs(fd_solve(n) - ref[:n]).max()
    print(f"n={n:3d}  谱方法误差 {es:.2e}   二阶差分误差 {ef:.2e}")
​

典型输出：n = 8 时谱方法误差约 10⁻³，n = 16 时跳到 10⁻⁷，n = 32 时逼近机器精度——每翻一倍，多个数量级；二阶差分同表按四倍下降，追谱方法需要上万格点。这份对比表是"光滑问题选谱方法"的全部理由；反过来，解带棱角时两条曲线的差距会迅速缩窄，理由见下。

两条边界：吉布斯与混淆

边界一：吉布斯现象。谱方法的截断对不连续函数失效：方波的傅里叶部分和在跳点附近顽固地过冲约 9% 的跳变高度，加密不会消除过冲，只把它压得更窄（第一章争论里的老朋友——有限和的极限可以不光滑，但逼近方式有病态）。推论：解一旦含间断或陡峭边界层，全局基就不是合适的语言，激波问题回到有限体积。

边界二：混淆误差。非线性问题里要在物理空间算乘积（如 u·u_x），两个截断到 n 个模式的函数相乘，乘积的波数超过奈奎斯特上限的部分会被"折叠"到低波数上冒充合法模式——这个冒名顶替叫混淆误差。处理办法要么定期把高波数系数清零（所谓 2/3 去混淆规则，谱方法湍流模拟的标准动作），要么换用无混淆的变换组合。第五章孤子方程的谱方法实验里会实际执行这条纪律。

💡 关键直觉：谱方法把"离散化"换成了"截断"，把"差商近似"换成了"基函数代数"。精度因此与解的光滑度绑定——光滑是燃料，间断是天花板。选它之前先估计解的正则性，2.2 节的三族性格档案是现成的判断依据。

非周期怎么办：切比雪夫与谱系数诊断

傅里叶基只服务周期问题，非周期的区间问题换切比雪夫多项式上岗：把区间映射到 [−1, 1]，用切比雪夫基展开，配点取切比雪夫聚集点（cos 均布——在两端天然加密，正好对付边界层）。求导仍在系数空间精确执行（三项递推），代数收敛换指数收敛的红利不变，只是代数比傅里叶版稍繁琐。非周期热方程、区别于周期湍流的槽道流计算，用的都是这套。

谱方法还有一个独特的"体检仪"：谱系数衰减曲线。把任何函数变换到谱系数空间，看系数随波数的下降速度——指数下降说明解析光滑，幂律下降说明有限光滑度，卡在某个水平不动则说明有噪声或间断。这条曲线比任何"看起来光不光滑"的主观判断都可靠，工程上用它决定"这个问题配不配用谱方法"，也用它诊断噪声水平（不动的那层平台高度就是噪声幅度）。

import numpy as np

# 谱系数诊断：光滑函数 vs 含噪函数 的系数衰减曲线
n = 256; x = np.linspace(-1, 1, n, endpoint=False)
smooth = np.exp(np.sin(2*np.pi*x))                    # 解析函数
noisy = smooth + 0.02*np.random.default_rng(0).standard_normal(n)
for name, f in [("光滑", smooth), ("含噪", noisy)]:
    c = np.abs(np.fft.rfft(f))
    print(f"{name}: 系数从 {c[1]:.2e} 衰减到 {c[-1]:.2e}，"
          f"尾部平台 {'无（继续下降）' if c[-1] < c[10]*1e-8 else '有（噪声水平约 %.0e）' % c[-1]}")
​

光滑函数的系数一路下探十个量级以上；含噪函数则先指数下降、随后压在噪声幅度构成的水平平台上不再动——平台高度就是噪声的"谱指纹"。这张诊断图是谱方法工程实践里最常用的一张表，成本只有一次 FFT。

本节要点

• 全局基与系数空间：微积分变成系数乘法，求导精确、误差只剩截断；
• 指数收敛：解析解的系数指数衰减，n 每翻倍误差多降一个数量级；
• FFT 地基：快速变换让每步代价 n·log n，谱方法由此实用；
• 吉布斯边界：间断函数的部分和过冲不灭，带棱角的解退回局部方法；
• 混淆边界：非线性乘积的高波数折叠需去混淆纪律，光滑度再次成为分水岭。

第四章四大方法各就各位。但它们默认的世界是线性的：解不会自己变陡、方程可以叠加、模式互不打扰。下一章抽掉线性这根拐杖——叠加失效之后，方程自己会长出激波、孤子这些原型里没有的角色。