谱方法

有限差分方法在每个网格点附近使用局部多项式近似导数；谱方法（Spectral Method）则选取一组定义在整个计算区域上的全局基函数，用有限项展开逼近未知函数。对于足够光滑的函数，谱展开系数通常随阶数快速衰减，因此谱方法能够以较少的空间节点获得很高的精度。

考虑偏微分方程及边界条件

$$Lu(\mathbf x)=s(\mathbf x),\quad \mathbf x\in U,$$

$$Bu(\mathbf x)=0,\quad \mathbf x\in\partial U,$$

其中 $L$ 是微分算子，$B$ 是边界算子。选取有限维近似空间 $P_N$ 及其一组试探函数（trial functions）$\{\phi_n\}_{n=0}^{N}$，用

$$u_N(\mathbf x)=\sum_{n=0}^{N}\tilde u_n\phi_n(\mathbf x)$$

逼近精确解。代入控制方程后产生残差

$$R(\mathbf x)=Lu_N(\mathbf x)-s(\mathbf x).$$

数值方法的任务是确定展开系数 $\tilde u_n$，使残差在某种意义下尽可能小。

强形式、弱形式与加权残差

原偏微分方程要求在区域内每一点成立，称为强形式（strong form）。若取权函数 $w$ 并要求

$$\int_U\left(Lu-s\right)w\,dU=0,$$

则得到相应的积分形式。对于近似解，通常不能使残差处处为零，但可以令它与一组检验函数（test functions）$\{\chi_n\}$ 正交：

$$(\chi_n,R)=\int_U\chi_nR\,dU=0, \qquad n=0,1,\ldots,N.$$

代入展开式可得

$$\sum_{k=0}^{N}\tilde u_k(\chi_n,L\phi_k) =(\chi_n,s).$$

记

$$L_{nk}=(\chi_n,L\phi_k),\qquad s_n=(\chi_n,s),$$

便得到线性代数方程组

$$\mathbf L\tilde{\mathbf u}=\mathbf s.$$

根据检验函数或残差约束方式的不同，可以得到多种方法：

• Galerkin 方法：取 $\chi_n=\phi_n$，即检验函数与试探函数相同；
• Petrov-Galerkin 方法：检验空间与试探空间不同；
• 最小二乘法：使残差范数 $\|R\|_2$ 最小；
• 配点法：选取节点 $x_j$，直接要求 $R(x_j)=0$。

傅里叶伪谱法采用傅里叶函数作为全局基，并利用配点处的函数值完成空间微分。

伪谱方法

经典谱方法直接求解谱展开系数；伪谱方法（Pseudospectral Method, PSM）则在物理空间节点和谱空间系数之间来回变换：

$$\{u(x_j)\} \xleftrightarrow{\mathrm{DFT/IDFT}} \{\hat u_k\}.$$

空间导数在谱空间中计算，而介质参数乘法、震源加入和部分边界处理通常在物理空间中完成。这种“变换—运算—逆变换”的方式避免了显式构造和求解稠密谱矩阵。

谱方法和伪谱方法都是全局方法：某个节点处的导数一般依赖整个计算区间上的所有节点。其优点是高精度，代价是周期性、边界处理和并行通信更复杂。

傅里叶级数与傅里叶变换

周期函数的谱展开

设 $f(x)$ 是定义在长度为 $L$ 的区间上的周期函数，满足 $f(x+L)=f(x)$，可以写成复指数傅里叶级数

$$f(x)=\sum_{m=-\infty}^{\infty}\hat f_m e^{ik_mx}, \qquad k_m=\frac{2\pi m}{L},$$

其中

$$\hat f_m=\frac{1}{L}\int_0^L f(x)e^{-ik_mx}\,dx.$$

复指数函数满足正交关系

$$\int_0^L e^{ik_mx}e^{-ik_nx}\,dx =L\delta_{mn},$$

因此不同波数分量可以独立提取。对于实函数，谱系数满足共轭对称性

$$\hat f_{-m}=\hat f_m^*.$$

离散傅里叶变换

在区间 $[0,L)$ 上取 $N$ 个等间距节点

$$x_j=j\Delta x,\qquad \Delta x=\frac{L}{N},qquad j=0,1,\ldots,N-1.$$

一种常用的离散傅里叶变换（DFT）约定为

$$\hat f_m=\sum_{j=0}^{N-1}f_j e^{-i2\pi jm/N},$$

逆变换为

$$f_j=\frac{1}{N}\sum_{m=0}^{N-1}\hat f_m e^{i2\pi jm/N}.$$

不同软件可能把归一化系数放在正变换、逆变换或两者中。推导和编程时必须使用一致的约定。

对于偶数 $N$，按 FFT 数组顺序排列的波数通常写为

$$k=\frac{2\pi}{L} \left(0,1,\ldots,\frac N2, -\frac N2+1,\ldots,-1\right).$$

Nyquist 波数为

$$k_N=\frac{\pi}{\Delta x}.$$

它对应网格能够区分的最高空间频率。任何高于 Nyquist 波数的分量都会折叠到低波数区间，形成混叠（aliasing）。

傅里叶伪谱求导

根据傅里叶变换的微分性质，若

$$u(x)=\sum_m\hat u_m e^{ik_mx},$$

则

$$\frac{\partial u}{\partial x} =\sum_m ik_m\hat u_m e^{ik_mx},$$

$$\frac{\partial^2u}{\partial x^2} =\sum_m-k_m^2\hat u_m e^{ik_mx}.$$

因此离散形式可以写为

$$D_xu=\mathcal F^{-1}\left[ik_x\mathcal F(u)\right],$$

$$D_{xx}u=\mathcal F^{-1}\left[-k_x^2\mathcal F(u)\right].$$

二维情况下，分别沿两个方向变换：

$$\frac{\partial u}{\partial x} =\mathcal F^{-1}_{2D}\left[ik_x\hat u(k_x,k_z)\right],$$

$$\nabla^2u =\mathcal F^{-1}_{2D} \left[-(k_x^2+k_z^2)\hat u(k_x,k_z)\right].$$

计算步骤

以一阶空间导数为例：

1. 在等间距网格上取得函数值 $u_j$；
2. 计算 $\hat u=\operatorname{FFT}(u)$；
3. 对每个波数分量计算 $ik\hat u$；
4. 计算 $u_x=\operatorname{IFFT}(ik\hat u)$；
5. 对实值问题，舍去舍入误差产生的极小虚部。

伪代码为：

k  = build_wavenumber_array(N, dx)
uh = FFT(u)
ux = real(IFFT(i * k * uh))

Nyquist 分量的处理

对偶数节点，Nyquist 模式在离散网格上表现为 $(-1)^j$。其一阶导数在节点上的表示存在符号歧义。为了保证实函数求导结果仍为实数，实际程序常把一阶导数乘子中的 Nyquist 分量设为零。二阶导数乘子 $-k^2$ 没有同样的奇对称问题，但仍需保持波数数组与 FFT 库的排列一致。

快速傅里叶变换

直接计算 DFT 需要 $O(N^2)$ 次运算。快速傅里叶变换（Fast Fourier Transform, FFT）利用复指数的周期性和对称性递归分解计算，将复杂度降低为

$$O(N\log_2N).$$

对二维 $N_x\times N_z$ 网格，二维 FFT 可以分解为先沿一个方向、再沿另一个方向的一维 FFT，总体计算量约为

$$O\left(N_xN_z\log(N_xN_z)\right).$$

FFT 并不是新的数学变换，而是计算 DFT 的高效算法。节点数取具有较小质因子的整数时通常效率较高，但现代 FFT 库也可以处理一般长度。

傅里叶伪谱法求解声波方程

考虑常密度二维声波方程

$$\frac{\partial^2p}{\partial t^2} =c^2(x,z)\left( \frac{\partial^2p}{\partial x^2} +\frac{\partial^2p}{\partial z^2} \right)+s(x,z,t).$$

空间二阶导数采用傅里叶伪谱算子，时间二阶导数采用中心差分：

$$\frac{p^{n+1}-2p^n+p^{n-1}}{\Delta t^2} =c^2\nabla^2_{PS}p^n+s^n.$$

得到时间推进格式

$$p^{n+1} =2p^n-p^{n-1} +\Delta t^2 \left[c^2\nabla^2_{PS}p^n+s^n\right],$$

其中

$$\nabla^2_{PS}p^n =\mathcal F^{-1}_{2D}\left[ -(k_x^2+k_z^2)\mathcal F_{2D}(p^n) \right].$$

每一时间步的基本流程为：

1. FFT2(p_n)
2. 乘以 -(kx^2 + kz^2)
3. IFFT2 得到 Laplacian(p_n)
4. 加入介质参数和震源
5. 用中心差分推进到 p_(n+1)
6. 施加吸收边界或阻尼层
7. 记录接收点数据并滚动时间层

虽然空间离散具有谱精度，时间推进仍只有二阶精度。整体误差由空间误差、时间误差、边界误差和浮点误差共同决定。

频散与稳定性

空间频散

对可解析的平面波 $e^{i(kx-\omega t)}$，傅里叶伪谱求导直接使用真实波数 $k$，因此理想的空间微分本身不引入有限差分型的修正波数误差。与二阶或高阶有限差分相比，傅里叶伪谱法可以用更少的每波长节点数描述光滑波场。

但是，“空间无有限差分频散”并不等于整个算法完全无频散。二阶时间差分仍给出

$$\frac{4}{\Delta t^2} \sin^2\left(\frac{\omega\Delta t}{2}\right) =c^2k^2,$$

因此时间离散会造成相速度误差。

稳定性条件

由上式要求右端不超过正弦函数的最大值，可得

$$c_{\max}\Delta t\,k_{\max}\leq 2.$$

二维等间距网格中

$$k_{\max}=\sqrt{k_{x,N}^2+k_{z,N}^2} =\frac{\pi\sqrt2}{\Delta x},$$

因而有保守条件

$$\Delta t\leq \frac{2\Delta x}{\pi\sqrt2\,c_{\max}}.$$

若两个方向步长不同，应使用

$$\Delta t\leq \frac{2}{\pi c_{\max} \sqrt{\Delta x^{-2}+\Delta z^{-2}}}.$$

稳定性条件只保证误差不发散。为了限制时间频散，实际时间步长通常还要取得更小。

周期性、边界和人为噪声

隐含周期边界

DFT 把有限数组视为一个周期信号。如果波场从计算区域右边界离开，它会从左边界重新进入；上下边界同理。因此实际波场模拟必须在物理区域外设置足够宽的吸收层或阻尼层，使波在到达周期拼接处之前明显衰减。

Gibbs 现象

当函数或介质参数存在间断时，截断傅里叶级数会在间断附近产生振荡。增加节点数会使振荡集中到更窄范围，但最大过冲不会简单消失。常见缓解方式包括：

• 平滑过陡的模型界面；
• 对高波数分量施加谱滤波；
• 对非线性乘积采用去混叠；
• 在间断附近改用局部方法或混合方法。

奇偶解耦与 artefacts

在某些一阶方程的同位网格傅里叶求导中，最高波数分量、介质间断和变量配置会共同产生棋盘状或条纹状人为噪声。把不同变量布置在交错网格上，相应的谱导数乘子变为带半网格相移的形式，例如

$$D_x^{+}u =\mathcal F^{-1} \left[ik\,e^{ik\Delta x/2}\hat u\right],$$

$$D_x^{-}u =\mathcal F^{-1} \left[ik\,e^{-ik\Delta x/2}\hat u\right].$$

正负半网格相移使导数与变量实际位置一致，可以减弱奇偶节点解耦和相关人为噪声。

与有限差分方法的比较

特征	有限差分法	傅里叶伪谱法
空间近似	局部差分模板	全局傅里叶展开
节点	通常等间距	等间距
空间精度	由差分阶数决定	对光滑周期函数具有谱精度
每波长节点数	低阶格式要求较多	通常较少
边界	可局部构造	隐含周期延拓，需额外处理
介质间断	相对灵活	容易出现 Gibbs 振荡
计算	局部运算，易于区域分解	FFT 高效，但涉及全局数据交换
存储	只需少量时间层和模型	与差分法相近

两类方法的时间推进可以完全相同，主要差别在于空间导数的计算方式。对于规则区域、光滑介质和高精度波场模拟，伪谱法优势明显；对于复杂边界、强间断介质和大规模分布式并行，局部高阶差分往往更加灵活。

版权所有

版权归属：Staaaaaaaaar

许可证：署名-非商业性-相同方式共享 4.0 国际 (CC-BY-NC-SA-4.0)

编辑此页

贡献者: Staaaaaaaaar

上一页高阶差分算子

下一页切比雪夫伪谱方法