切比雪夫伪谱方法以切比雪夫多项式作为全局基函数，在有限区间上选取非等间距配点，用全体节点函数值的加权组合计算空间导数。它保留了伪谱方法的高精度，同时避免傅里叶基天然要求周期边界的问题，适合含有明确物理边界的有限区域。

函数近似与正交函数基

设函数 $f(x)$ 定义在区间 $[a,b]$ 上，用有限个基函数近似：

$$f(x)\approx g_N(x)=\sum_{n=0}^{N}a_n\Phi_n(x).$$

若基函数关于权函数 $w(x)>0$ 满足

$$\langle \Phi_m,\Phi_n\rangle_w =\int_a^b\Phi_m(x)\Phi_n(x)w(x)\,dx =\gamma_n\delta_{mn},$$

则称它们构成加权正交函数集。展开系数可以由投影得到：

$$a_n=\frac{\langle f,\Phi_n\rangle_w}{\gamma_n}.$$

正交性使不同阶基函数之间相互独立，能够减小求解展开系数时的耦合。傅里叶方法选取三角函数或复指数函数；切比雪夫方法则选取定义在 $[-1,1]$ 上的正交多项式。

切比雪夫多项式

第一类切比雪夫多项式

第一类切比雪夫多项式（Chebyshev polynomials of the first kind）定义为

$$T_n(x)=\cos(n\arccos x),\qquad x\in[-1,1].$$

令 $x=\cos\theta$，则

$$T_n(\cos\theta)=\cos(n\theta).$$

前几阶多项式为

$$T_0(x)=1,$$

$$T_1(x)=x,$$

$$T_2(x)=2x^2-1,$$

$$T_3(x)=4x^3-3x.$$

利用余弦的递推关系可以得到

$$T_{n+1}(x)=2xT_n(x)-T_{n-1}(x).$$

因此不需要显式展开高次幂，就可以稳定地逐阶计算 $T_n(x)$。

正交性

第一类切比雪夫多项式在 $[-1,1]$ 上关于权函数

$$w(x)=\frac{1}{\sqrt{1-x^2}}$$

正交：

$$\int_{-1}^{1} \frac{T_m(x)T_n(x)}{\sqrt{1-x^2}}\,dx = \begin{cases} 0, & m\ne n,\\ \pi, & m=n=0,\\ \pi/2, & m=n\ne0. \end{cases}$$

作代换 $x=\cos\theta$ 后，上式转化为余弦函数在 $[0,\pi]$ 上的正交积分，这说明切比雪夫展开与余弦展开具有紧密联系。

零点和极值点

$T_N(x)$ 的零点为

$$x_j=\cos\left(\frac{2j+1}{2N}\pi\right), \qquad j=0,1,\ldots,N-1.$$

$T_N(x)$ 的极值点为

$$x_j=\cos\left(\frac{j\pi}{N}\right), \qquad j=0,1,\ldots,N.$$

后一组节点同时包含 $x=\pm1$，称为 Chebyshev-Gauss-Lobatto 节点，最适合需要显式施加边界条件的伪谱计算。

切比雪夫配点与区间映射

非等间距节点

Chebyshev-Gauss-Lobatto 节点

$$x_j=\cos\left(\frac{j\pi}{N}\right)$$

在区间中心最稀疏，在两端最密集。相邻节点间距为

$$\Delta x_j=x_j-x_{j+1}.$$

中心附近最大间距约为 $O(N^{-1})$，边界附近最小间距约为 $O(N^{-2})$。端点聚集可以减弱等距高次多项式插值的 Runge 现象，并提高边界附近的分辨率，但最小网格间距也会使显式时间推进的稳定条件更加严格。

映射到一般区间

切比雪夫节点定义在标准区间 $\xi\in[-1,1]$。若物理区间为 $x\in[a,b]$，采用线性映射

$$x=\frac{a+b}{2}+\frac{b-a}{2}\xi,$$

$$\xi=\frac{2x-a-b}{b-a}.$$

由链式法则

$$\frac{\partial}{\partial x} =\frac{2}{b-a}\frac{\partial}{\partial\xi},$$

$$\frac{\partial^2}{\partial x^2} =\left(\frac{2}{b-a}\right)^2 \frac{\partial^2}{\partial\xi^2}.$$

因此在标准区间构造的微分矩阵必须乘以相应尺度因子，才能用于物理坐标。

切比雪夫多项式插值

给定 $N+1$ 个 Lobatto 节点及函数值 $f_j=f(x_j)$，构造次数不超过 $N$ 的插值多项式

$$p_N(x)=\sum_{j=0}^{N}f_j\ell_j(x),$$

其中 $\ell_j(x)$ 是 Lagrange 基函数，满足

$$\ell_j(x_k)=\delta_{jk}.$$

也可以写成截断切比雪夫展开

$$p_N(x)=\sum_{n=0}^{N}{}' a_nT_n(x),$$

其中撇号表示端点系数采用相应的半权约定。由于 $x_j=\cos(j\pi/N)$，节点值和切比雪夫系数之间的变换本质上是离散余弦变换，可以借助 FFT 高效计算。

对于解析且光滑的函数，切比雪夫系数随阶数近似指数衰减，插值误差也快速下降，这通常称为谱收敛。若函数只有有限阶光滑性，收敛速度退化为代数型；若存在间断，则会出现 Gibbs 振荡。

注

所谓“谱精度”并不意味着有限节点下误差严格为零，而是指对足够光滑的函数，误差随节点数增长得比任意固定阶幂律更快。函数的光滑性决定实际收敛速度。

切比雪夫微分矩阵

从插值多项式得到导数

对插值表达式求导：

$$p_N'(x_i)=\sum_{j=0}^{N}f_j\ell_j'(x_i).$$

定义一阶微分矩阵 $\mathbf D$：

$$D_{ij}=\ell_j'(x_i),$$

则所有节点处的导数可以写为

$$\mathbf f'=\mathbf D\mathbf f.$$

这表明任一节点处的导数由全体节点函数值共同决定，体现了切比雪夫伪谱法的全局性。

一阶微分矩阵元素

令

$$c_j= \begin{cases} 2, & j=0\ \text{或}\ j=N,\\ 1, & 1\leq j\leq N-1, \end{cases}$$

对于 $i\ne j$，一阶微分矩阵可写为

$$D_{ij}=\frac{c_i}{c_j} \frac{(-1)^{i+j}}{x_i-x_j}.$$

内部对角元素为

$$D_{ii}=-\frac{x_i}{2(1-x_i^2)}, \qquad i=1,\ldots,N-1.$$

两个端点的对角元素为

$$D_{00}=\frac{2N^2+1}{6}, \qquad D_{NN}=-\frac{2N^2+1}{6}.$$

若节点按 $x_0=1$ 到 $x_N=-1$ 的顺序排列，上述符号成立；若程序把节点升序排列，需要同步重排矩阵。数值实现中还可以利用

$$D_{ii}=-\sum_{j\ne i}D_{ij}$$

计算对角元素，使常数函数的数值导数严格为零并减小舍入误差。

二阶导数矩阵可以取

$$\mathbf D^{(2)}=\mathbf D^2,$$

更高阶导数同理。但直接矩阵乘方在高阶和大规模问题中可能放大舍入误差，需要根据算法规模选择递推或变换形式。

物理区间上的微分矩阵

若 $\mathbf D_\xi$ 是标准区间上的微分矩阵，则

$$\mathbf D_x=\frac{2}{b-a}\mathbf D_\xi,$$

$$\mathbf D_{xx}=\left(\frac{2}{b-a}\right)^2 \mathbf D_\xi^2.$$

二维张量积网格中，沿两个坐标方向的求导可以分别表示为

$$\frac{\partial\mathbf U}{\partial x} =\mathbf D_x\mathbf U,$$

$$\frac{\partial\mathbf U}{\partial z} =\mathbf U\mathbf D_z^{T},$$

具体左右乘关系取决于数组中两个坐标轴的排列。

切比雪夫伪谱法求解波动方程

空间与时间离散

以二维声波方程为例：

$$\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).$$

若两个方向均采用切比雪夫配点，空间 Laplace 算子可离散为

$$\nabla^2\mathbf P \approx \mathbf D_{xx}\mathbf P +\mathbf P\mathbf D_{zz}^{T}.$$

时间上使用二阶中心差分：

$$\mathbf P^{n+1} =2\mathbf P^n-\mathbf P^{n-1} +\Delta t^2\left[ \mathbf C^2\odot \left(\mathbf D_{xx}\mathbf P^n +\mathbf P^n\mathbf D_{zz}^{T}\right) +\mathbf S^n \right],$$

其中 $\odot$ 表示逐点乘法。

基本计算流程

1. 在标准区间生成 Chebyshev-Gauss-Lobatto 节点；
2. 将节点映射到物理区域；
3. 构造一阶或二阶微分矩阵并进行尺度变换；
4. 在非均匀节点上赋值介质参数和初始条件；
5. 计算空间导数，加入震源；
6. 按时间差分格式推进；
7. 在每个时间步施加边界条件；
8. 输出波场快照和接收点记录。

边界条件处理

Lobatto 节点包含计算区间端点，因此边界条件可以直接在端点行上施加。

Dirichlet 边界

若边界上要求

$$u(a,t)=g_a(t),\qquad u(b,t)=g_b(t),$$

可在每个时间步直接覆盖端点值，或在方程组中把对应矩阵行替换为单位行。

Neumann 边界

若边界上要求

$$u_x(a,t)=q_a(t),$$

则利用微分矩阵边界行：

$$\sum_{j=0}^{N}D_{Nj}u_j=q_a$$

或相应端点行，具体取决于节点顺序。弹性波自由表面要求法向应力和切向应力为零，也可通过边界处的谱导数关系修正速度或应力更新。

人工边界

切比雪夫方法解决了有限区间上的配点问题，但不会自动吸收向外传播的波。在截断的开放区域中仍需设置吸收层、阻尼函数或其他无反射边界条件。边界处节点高度聚集，吸收参数应按实际非均匀间距设计。

震源处理

点震源的空间尺度通常小于网格尺度。若直接把源加在单一节点上，会激发接近网格极限的高波数成分并产生振荡。因此实际计算常采用具有有限宽度的平滑空间分布，例如 Gaussian 函数：

$$S(x,z,t)=A(t) \exp\left[-\frac{(x-x_s)^2+(z-z_s)^2}{\sigma^2}\right].$$

常用 Ricker 子波为

$$A(t)= \left[1-2\pi^2f_0^2(t-t_0)^2\right] \exp\left[-\pi^2f_0^2(t-t_0)^2\right].$$

震源中心一般不恰好落在配点上，需要按插值权重或平滑核把源分配给邻近节点，并注意非均匀节点对应的积分权重。

稳定性与计算代价

切比雪夫节点的最小间距位于边界附近，约随 $N^{-2}$ 缩小。显式时间推进的稳定步长受最小间距或微分矩阵最大特征值控制，通常比等距网格更加严格。对二阶空间导数问题，最大特征值会随 $N$ 快速增大，因此高分辨率下时间步长可能成为主要瓶颈。

直接微分矩阵乘法的单次成本为 $O(N^2)$；二维张量积计算的成本也明显高于局部有限差分。利用切比雪夫系数与离散余弦变换可以把部分运算加速到 $O(N\log N)$，但边界条件和非均匀介质仍需谨慎实现。

与傅里叶伪谱方法的比较

特征	傅里叶伪谱法	切比雪夫伪谱法
基函数	复指数或三角函数	切比雪夫多项式
节点	等间距	两端密集、中心稀疏
自然区间	周期区间	有限区间 $[-1,1]$
边界	周期边界最自然	Lobatto 节点直接包含边界
求导	FFT 后乘以 $ik$	微分矩阵或余弦变换
时间步长	由均匀步长控制	易受边界最小间距限制
适用问题	光滑周期模型、规则区域	明确有限边界、非周期问题

两者都属于全局方法，对光滑函数具有高精度，对不连续函数都可能出现 Gibbs 现象。傅里叶方法计算效率通常更高；切比雪夫方法的有限区间边界处理更自然，但节点聚集和稳定步长限制更突出。

版权所有

版权归属：Staaaaaaaaar

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

编辑此页

贡献者: Staaaaaaaaar

上一页傅里叶伪谱方法

下一页二维地球介质中地震波传播计算