Fourier 伪谱法

形式陈述 ​

周期网格上的函数值，怎样通过频率乘子计算导数？为避开 Nyquist 端点约定，先取奇数节点 $N=2K+1$，其中 $K\ge 0$ 为整数，$x_j=2\pi j/N$。离散 Fourier 系数与插值为

$${\hat{u}}_k=\frac{1}{N}\sum _{j=0}^{N-1}{u}_j{e}^{-ik{x}_j},p_N(x)=\sum _{k=-K}^K{\hat{u}}_k{e}^{ikx}.$$

按Fourier 模态公理库Fourier 级数Fourier series · 傅里叶级数把周期函数投影到整数频率的正交指数基上所得的离散频谱展开。逐项微分，导数系数为 $ik{\hat{u}}_k$，二阶导数系数为 $-k^2{\hat{u}}_k$。使用FFT公理库快速 Fourier 变换Fast Fourier transform · FFT利用单位根的偶奇分解在 $O(n\log n)$ 时间计算离散 Fourier 变换。在节点值与系数间转换，一次微分的复数运算成本为 $O(N\log (N+1))$；下面给出奇数长度如何调用二次幂 FFT 的具体转换。

Fourier 伪谱法用这一导数处理微分方程的空间项，非线性项通常转到节点空间逐点相乘，再变换回系数空间。空间离散后得到系数或节点值的常微分方程，还需另选时间积分器。

奇数长度如何调用二次幂 FFT ​

先按余数 $k=0,\dots ,N-1$ 存放变换输出，最后把 $k>K$ 重排为频率 $k-N$。由 $2jk=j^2+k^2-(j-k{)}^2$，令

$$a_j=u_j{e}^{-\pi i{j}^2/N},b_{\ell }=e^{\pi i{\ell }^2/N}(-(N-1)\le \ell \le N-1),$$

就有

$${\hat{u}}_k=\frac{e^{-\pi i{k}^2/N}}{N}\sum _{j=0}^{N-1}{a}_j{b}_{k-j}.$$

把 $b_{\ell }$ 存为数组 $B_{\ell +N-1}$，则右侧和是长度 $N$ 与 $2N-1$ 数组的线性卷积在下标 $k+N-1$ 处的值。选不小于 $3N-2$ 的最小二次幂 $M$，补零后用 FFT 的两次正变换、逐点相乘和一次逆变换计算；逆变换除以 $M$，所得卷积不发生绕回。这里 FFT 工具采用正号正变换，但它的正逆变换配对计算同一线性卷积，外层的负号和 $1/N$ 由上式明确给出。

反向合成 $u_j=\sum _{k=0}^{N-1}{\hat{u}}_k{e}^{2\pi ijk/N}$ 时，把这些 chirp 指数的符号全部反转，并去掉外层 $1/N$，即可使用同一卷积构造。由于 $M=O(N)$，转换成本仍为 $O(N\log (N+1))$。与 FFT 的运算计数一致，这里采用单位成本的精确复数加乘，并假设所需单位根已给定；chirp 因子可按相邻项比值逐次生成，只需 $O(N)$ 次运算。浮点实现还需另计根与蝶形的舍入误差。

直觉

每个 Fourier 模态都是一根已知振动频率的弦。求导不需要重新拟合整条曲线，只需把每根弦的振幅乘上对应频率与相位因子。周期首尾自然相接，因此无需另添两个互不相干的端点。

这里的导数是有限频带插值的精确导数。只有当采样充分解析目标函数时，它才同时是原函数导数的好近似。

例子与边界

一个非平凡而可完全核对的导数 ​

取 $u(x)=\sin \left(2x\right)+\frac{1}{2}\cos \left(3x\right)$，用 $N=9$ 个点，即 $K=4$。所有频率均在可表示范围内，因此离散系数精确保留这两个模态，谱微分得到

$$u^{\prime }(x)=2\cos \left(2x\right)-\frac{3}{2}\sin \left(3x\right).$$

在 $x=0$，导数为 $2$；在 $x=2\pi /3$，导数为 $-1$。把解析导数在各网格节点的值与 FFT 结果逐点核对，可以同时检查频率排序、正负号及变换归一化。有限差分可作为另一个离散方法比较收敛速度，解析导数才是本例的误差基准。

由热方程看时间误差 ​

对周期热方程 $u_t=\nu u_{xx}$，取 $\nu >0$、时间步长 $\mathrm{\Delta }t>0$，模态满足 ${\hat{u}}_k^{\prime }=-\nu k^2{\hat{u}}_k$，故精确时间推进是乘 $e^{-\nu k^2\mathrm{\Delta }t}$。若改用显式 Euler，放大因子为 $1-\nu k^2\mathrm{\Delta }t$，稳定性要求最高频率满足 $\nu K^2\mathrm{\Delta }t\le 2$。空间谱精度并未免除时间步长限制。

非周期数据也要注意周期延拓。比如 $u(x)=x$ 在 $[0,2\pi )$ 内很光滑，但首尾相接出现跳跃，Fourier 插值就会在接口产生振荡。应改用适合非周期区间的基，或先处理边界失配。

推论与应用

实现流程是采样、FFT、乘微分因子、逆 FFT，最后在方程中组合导数与非线性项。长度为 $L>0$ 的周期区间应使用频率 $2\pi k/L$；实值数据还应保持共轭对称，使结果只剩舍入量级的虚部。

偶数节点时，最高 Nyquist 模态对一阶导数需要一致约定，常把其导数系数置零。非线性乘积产生超出原频带的频率，必须检查混叠与去混叠公理库谱混叠与去混叠Spectral aliasing · Dealiasing从频率模网格数的折叠推导二次非线性的补零与截断规则，完整算出一个伪低频反例。，这是伪谱离散中特有的一步。

参考资料

• Trefethen, Spectral Methods in MATLAB, Chapters 2–4，周期谱微分与 FFT。
• KTH spectral-methods notes, Pseudo-spectral method and aliasing errors, Sections 7.1–7.2，系数空间与节点空间的切换。