伪谱方法

伪谱方法（pseudospectral method）在基函数系数与物理节点值之间转换，以谱表示完成微分，再在节点上计算乘积等非线性运算。周期问题通常结合 Fourier 变换使用。它能避开直接计算大量系数卷积，但节点数有限时，高频会被误认成低频，产生混叠误差。

同一组节点辨认不出的两条曲线

取八个等距节点 $x_{j} = 2 π j / 8$ ， $j = 0, \dots, 7$ 。在这些节点上 $\cos 6 x_{j} = \cos (8 x_{j} - 2 x_{j}) = \cos 2 x_{j} .$ 因此节点数据无法区分频率六和频率二。这不是浮点误差，而是有限采样本身的等价关系。

现有函数 $u (x) = \sin 3 x$ 。它的平方是 $u^{2} (x) = \frac{1}{2} - \frac{1}{2} \cos 6 x .$ 如果仅在八个节点上平方，再用这些数据重建低频插值多项式，便得到 $\frac{1}{2} - \frac{1}{2} \cos 2 x$ 。频率六被折成频率二。后者会错误地进入低频演化，所以只检查节点处两条曲线相等并不能验收非线性项。

真实平方函数与混叠后的低频函数在八个采样点完全重合但点间不同 — 实线为 sin²(3x)，虚线为 1/2−cos(2x)/2。标记处的八个采样值相同，高频信息在节点之间丢失。

系数、节点值与谱微分

设截断函数为 $u_{K} (x) = \sum_{k = - K}^{K} {\hat{u}}_{k} e^{i k x}$ 。其导数仍是有限和： ${u_{K}}^{'} (x) = \sum_{k = - K}^{K} i k {\hat{u}}_{k} e^{i k x} .$ 在频域将每个系数乘以 $i k$ 即可微分。若只保存节点值，先作离散 Fourier 变换得到系数，乘以 $i k$ 后再逆变换，便得到该插值多项式在节点处的导数。真实函数若含未解析的模式，这个结果仍只是插值导数，而非真导数的自动保证。

对于非线性方程，例如黏性 Burgers 方程 $u_{t} + (u^{2} / 2)_{x} = ν u_{x x}$ ，可将系数变到节点计算 $u^{2} / 2$ ，变回后乘以 $i k$ ，同时用 $- ν k^{2} {\hat{u}}_{k}$ 表示扩散。这种流程通常称为 Fourier 伪谱离散。时间积分和边界条件仍需单独选择。Fourier 系数与节点值转换的基本形式可参见 Driscoll 与 Braun 的三角插值章节；下文对乘积混叠直接作有限求和推导。^[1]

混叠关系的推导

在 $M$ 个等距节点上， $e^{i (k + ℓ M) x_{j}} = e^{i k x_{j}}$ ，因为 $e^{2 π i ℓ j} = 1$ 。设乘积的真实 Fourier 系数为 $c_{r}$ ，节点离散变换得到的第 $k$ 个系数为 ${\tilde{c}}_{k} = \sum_{ℓ \in ℤ} c_{k + ℓ M} .$ 可以直接证明：把 Fourier 和代入离散系数公式，内层出现 $\frac{1}{M} \sum_{j = 0}^{M - 1} e^{2 π i (r - k) j / M} .$ 若 $r - k$ 是 $M$ 的整数倍，每项均为一；否则它是首项为一、公比不等于一的有限几何级数，分子 $1 - e^{2 π i (r - k)}$ 为零，和为零。因此只有相差整数倍采样数的频率混在一起。

这一式子也给出消除二次混叠的条件。两个最高频率均为 $K$ 的函数相乘，只产生 $| r | \leq 2 K$ 。如果只保留 $| k | \leq K$ ，要避免其他模式满足 $r = k \pm M$ ，只需 $M > 3 K$ 。例如本例 $K = 3$ ，用十二点计算乘积再保留前三个频率，频率六不会再落入保留范围。

补零与三分之二截断

补零法把已知低频系数放入更长的频率数组，在较密节点上求乘积，再截回目标频率。二次非线性的常见实现称为“三分之二规则”或“二分之三补零”：保留的有效带宽与节点数之间留下去混叠余量。具体程序必须处理正负频率、偶数节点的 Nyquist 模式和变换归一化，名称不能替代下标检查。

本例在十二点上算出的平方仍含真实的常数和频率六。投影回 $| k | \leq 3$ 后只剩常数 $1 / 2$ 。这一步是正确的低频 Galerkin 投影，但并不等于完整的平方函数；高频模式是有意截去的。去混叠消除了错误折回，不能恢复被目标空间排除的所有信息。

有解变式取 $u = \cos 2 x$ 、 $K = 2$ 。平方为 $1 / 2 + \cos 4 x / 2$ ，投影回 $| k | \leq 2$ 仍为 $1 / 2$ 。若用五点直接相乘，频率四会折成频率负一，在实函数中形成错误的 $\cos x$ 项；至少七点满足 $M > 6$ ，可正确完成这个低频投影。

可复算的去混叠平方

下面用直接 Fourier 求和展示下标和归一化。输入字典保存 $- K$ 到 $K$ 的复系数，缺失项视为零；输出平方投影的系数字典，不修改输入。为了便于逐式核对，没有调用 FFT，计算量为 $O (M K)$ ，当 $M$ 与 $K$ 成比例时为 $O (K^{2})$ ，存储为 $O (M + K)$ 。实际 FFT 实现可将变换成本降至 $O (M \log M)$ 。

Python 3

import cmath
import math


def dealiased_square(coefficients, K, M):
    if not isinstance(K, int) or K < 0 or not isinstance(M, int) or M <= 3*K:
        raise ValueError("integer K and M with K nonnegative and M greater than 3K required")
    modes = {}
    for k, value in coefficients.items():
        if not isinstance(k, int) or abs(k) > K:
            raise ValueError("input mode outside retained band")
        value = complex(value)
        if not math.isfinite(value.real) or not math.isfinite(value.imag):
            raise ValueError("finite coefficients required")
        modes[k] = value
    xs = [2*math.pi*j/M for j in range(M)]
    squared = [sum(c * cmath.exp(1j*k*x) for k, c in modes.items())**2 for x in xs]
    result = {k: sum(v * cmath.exp(-1j*k*x) for x, v in zip(xs, squared))/M
              for k in range(-K, K+1)}
    if not all(math.isfinite(v.real) and math.isfinite(v.imag) for v in result.values()):
        raise OverflowError("nonfinite transformed product")
    return result


answer = dealiased_square({3: -0.5j, -3: 0.5j}, 3, 12)
print(round(answer[0].real, 12))  # 0.5
print(max(abs(v) for k, v in answer.items() if k != 0))  # near roundoff

这里的三倍带宽条件针对二次乘积。三次非线性有更宽的卷积支撑，需要相应增加节点或分步控制混叠；非周期多项式谱法也有不同的去混叠积分条件。对含激波的解，去混叠只能修正采样产生的频率错误，不能自行解决 Gibbs 振荡或代替熵稳定离散。

参考资料

↑ Tobin A. Driscoll, Richard J. Braun. Fundamentals of Numerical Computation: Trigonometric interpolation，复三角表示、离散 Fourier 矩阵及 FFT。

[fnc-1] Tobin A. Driscoll, Richard J. Braun. Fundamentals of Numerical Computation: Trigonometric interpolation，复三角表示、离散 Fourier 矩阵及 FFT。

[1]