第 4 章：离散化基础

离散化的意义

PDE 描述的是连续空间和连续时间中的流场，而计算机只能存储有限数量的数据，并用有限个网格点上的值近似连续函数。

一旦离散，连续 PDE 的很多性质都必须重新检查：

  • 是否仍然守恒；
  • 是否稳定；
  • 是否收敛；
  • 是否能够正确传播波；
  • 是否产生额外耗散或色散。

差分公式近似导数

差分格式利用邻近网格点的信息来近似导数。

例如中心差分：

ux≈uj+1−uj−12Δxu_x \approx \frac{u_{j+1}-u_{j-1}}{2\Delta x}

任何差分近似都必然存在截断误差。

所谓截断误差，截断的是 Taylor 展开的高阶项。

例如前向差分是一阶精度的：

(∂u∂x)i,j=ui+1,j−ui,jΔx+O(Δx)\left(\frac{\partial u}{\partial x}\right)_{i,j} = \frac{u_{i+1,j} - u_{i,j}}{\Delta x} + O(\Delta x)

这里的 O(Δx)O(\Delta x) 表示余下的项是 Δx\Delta x 的 同阶无穷小。

在数学上，小 (o) 和大 (O) 有严格的区别，Anderson 书上这里的表述“o(Δx)o(\Delta x) 表示与Δx\Delta x同量级”是不严谨的。

稳定性问题

真实数值计算中一定存在误差：

ujn=uexact,jn+ϵjnu_j^n = u_{\mathrm{exact},j}^n + \epsilon_j^n

误差可能来自：

  • 截断误差；
  • 舍入误差；
  • 初值误差；
  • 边界误差；
  • 迭代误差。

稳定性分析关心的核心问题是：

ϵn⟶ϵn+1\boxed{ \epsilon^n \longrightarrow \epsilon^{n+1} }

这一过程中误差会不会被差分格式自身不断放大。

Fourier 分析

一般的空间误差

ϵ(x)\epsilon(x)

可能非常复杂。

直接研究任意误差形状几乎不可能。

Fourier 分析的思想是：

复杂误差=很多不同波数的简单波形叠加\boxed{ \text{复杂误差} = \text{很多不同波数的简单波形叠加} }

写成复指数形式：

ϵ(x)=∑mAmeikmx\epsilon(x) = \sum_m A_m e^{ik_mx}

每一个

eikmxe^{ik_mx}

称为一个 Fourier 模态。

因此可以把一个复杂的误差传播问题转化成：

分别研究每一种波数的误差在数值格式中如何演化。

复指数表示的意义

Euler 公式给出

eikx=cos⁡(kx)+isin⁡(kx)e^{ikx} = \cos(kx)+i\sin(kx)

因此复指数同时包含了正弦和余弦信息。

复指数形式可以显著简化差分运算。

例如：

eik(x+Δx)=eikxeikΔxe^{ik(x+\Delta x)} = e^{ikx}e^{ik\Delta x}

因此：

网格上向右移动一个点，在 Fourier 空间里只是乘一个 eikΔxe^{ik\Delta x} 。

这使 Fourier 模态成为差分算子的天然分析基底。

标准复 Fourier 级数写成

f(x)=∑m=−∞∞cmeikmxf(x) = \sum_{m=-\infty}^{\infty} c_m e^{ik_mx}

对于实函数，有

c−m=cm∗c_{-m}=c_m^*

因此正波数和负波数不是相互独立的自由度。

可以只研究正波数，同时把负波数信息吸收到复振幅中。

因此 Anderson 的教材中写成只对

m=1,2,…m=1,2,\dots

求和。

波数与波长

波数定义为

k=2πλk=\frac{2\pi}{\lambda}

因此：

  • 小 kk 对应长波
  • 大 kk 对应短波

数值计算中特别需要关注高波数：

高波数意味着波长接近网格尺度，也是数值误差最难正确解析的区域。

最短可分辨波长是 2Δx2\Delta x

在网格上：

xj=jΔxx_j=j\Delta x

Fourier 模态变为

eikjΔxe^{ikj\Delta x}

如果将波数增加

2πΔx\frac{2\pi}{\Delta x}

则

ei(k+2π/Δx)jΔx=eikjΔxei2πje^{i(k+2\pi/\Delta x)j\Delta x} = e^{ikj\Delta x}e^{i2\pi j}

而由于 jj 为整数：

ei2πj=1e^{i2\pi j}=1

所以这两个波数在所有网格点上的取值完全一样。

因此离散网格无法区分

kk

和

k+2πΔxk+\frac{2\pi}{\Delta x}

于是独立可分辨波数只需取：

−πΔx≤k≤πΔx-\frac{\pi}{\Delta x} \le k\le \frac{\pi}{\Delta x}

最大可分辨波数为

kmax⁡=πΔx\boxed{ k_{\max}=\frac{\pi}{\Delta x} }

对应最短波长

λmin⁡=2Δx\boxed{ \lambda_{\min}=2\Delta x }

这就是 Nyquist 极限。

如果 k=πΔxk=\frac{\pi}{\Delta x}，也就是最短波长。

那么 eiπj=(−1)je^{i\pi j} =(-1)^j

所以网格上看到的是：

j:       0    1    2    3    4    5
         ●         ●         ●
         |         |         |
---------+----+----+----+----+----
              |         |
              ●         ●

value:  +1   -1   +1   -1   +1   -1

每两个网格间距重复一次：

λ=2Δx.\boxed{\lambda=2\Delta x}.

这已经是网格能表示的最快的交替变化。

你不可能让离散点做到比正负峰值交替变化得更快了。

若计算区域长度为

L=NΔxL=N\Delta x

Fourier 波数满足

km=2πmLk_m=\frac{2\pi m}{L}

最大可分辨波数为

πΔx\frac{\pi}{\Delta x}

令二者相等：

2πmmax⁡NΔx=πΔx\frac{2\pi m_{\max}}{N\Delta x} = \frac{\pi}{\Delta x}

得到

mmax⁡=N2\boxed{ m_{\max}=\frac{N}{2} }

于是有最高阶谐波为 N/2N/2

然而“可以表示”并不等于“可以准确解析”

当

λ=2Δx\lambda=2\Delta x

时，网格上的取值只能呈现：

+1,−1,+1,−1,…+1,-1,+1,-1,\dots

这虽然是合法的离散模态，但并不意味着差分格式能够准确计算它的导数。

例如中心差分

uj+1−uj−12Δx\frac{u_{j+1}-u_{j-1}}{2\Delta x}

作用在

uj=(−1)ju_j=(-1)^j

上时，会得到零。

但连续 Fourier 波的真实导数并不是零。

Nyquist 极限只是可表示极限，而不是高精度解析极限

实际 CFD 中，一个物理波长通常需要多个网格点才能被准确解析。

von Neumann 稳定性分析

假设某个单一 Fourier 误差模态为

ϵjn=AneikjΔx\epsilon_j^n = A_n e^{ikj\Delta x}

时间推进一步以后：

ϵjn+1=An+1eikjΔx\epsilon_j^{n+1} = A_{n+1} e^{ikj\Delta x}

定义放大因子

G(k)=An+1An\boxed{ G(k)= \frac{A_{n+1}}{A_n} }

于是经过 nn 步：

An=GnA0A_n = G^n A_0

因此：

  • ∣G∣<1|G|<1：模态衰减；
  • ∣G∣=1|G|=1：模态振幅不变；
  • ∣G∣>1|G|>1：模态增长。

稳定性通常要求：

∣G(k)∣≤1\boxed{ |G(k)|\le1 }

并且必须对所有允许波数都成立,一个不稳定模态就足以毁掉整个计算.

初始误差可能包含很多 Fourier 模态。

即使某个不稳定模态最初只有机器舍入量级：10−1510^{-15}

只要∣G∣>1|G|>1，误差是指数增长的。

放大因子在复数域上的意义

考虑热传导方程

ut=αuxxu_t=\alpha u_{xx}

连续 Fourier 模态

u=A(t)eikxu=A(t)e^{ikx}

代入以后：

A′(t)=−αk2A(t)A'(t) = -\alpha k^2A(t)

所以

A(t)=A(0)e−αk2tA(t) = A(0)e^{-\alpha k^2t}

放大因子始终为实数正数。

这意味着：

时间的推进主要改变 Fourier 模态的振幅，而不会把波形向左或向右搬运。

这里也体现了扩散的性质。

一般情况下：

G=∣G∣eiϕG=|G|e^{i\phi}

其中：∣G∣|G| 控制振幅变化；ϕ\phi 控制相位变化。

因此：

  • ∣G∣<1|G|<1 表示数值耗散；
  • 相位误差表示数值色散。

所以 Fourier 分析不仅可以研究稳定性，还可以研究：

  • 数值耗散
  • 数值色散
  • 不同波长的解析能力

但实放大因子不一定意味着完全没有相位变化

如果 G>0G>0 相位不变。

但如果 G<0G<0 则

G=∣G∣eiπG=|G|e^{i\pi}

相当于产生了 π\pi 的相位翻转。

因此更准确地说：

实放大因子意味着不存在连续传播型相位偏移，但可能出现 00 或 π\pi 的离散相位变化。

迎风格式

考虑线性对流方程：

ut+aux=0u_t+a u_x=0

如果

a>0a>0

信息向右传播。

因此某一点未来状态的信息来自左侧，也就是上游。

于是空间导数使用：

uj−uj−1Δx\frac{u_j-u_{j-1}}{\Delta x}

如果

a<0a<0

信息向左传播，则应使用右侧数据。

对于双曲型 PDE，信息沿特征线传播。

因此“上游”不是一个纯几何概念，而是由特征速度决定。

一维线性对流中：

dxdt=a\frac{dx}{dt}=a

所以 aa 的符号决定上游方向。

而高阶 CFD 中的迎风实际上是：

将系统分解到特征空间，再分别根据各个特征速度的符号选择信息来源。

CFL 条件

例如

ut+aux=0u_t+a u_x=0

定义 Courant 数

C=aΔtΔx\boxed{ C= \frac{a\Delta t}{\Delta x} }

一阶迎风显式格式通常要求

0≤C≤1\boxed{ 0\le C\le1 }

即：

aΔt≤Δxa\Delta t\le\Delta x

一个时间步内，物理信息传播距离为

aΔta\Delta t

而最简单的一阶差分格式一步只能利用一个网格间距范围内的信息。

这就是 CFL 条件的实际意义。

更一般地说，CFL 条件要求：

数值依赖域必须覆盖物理依赖域

如果真实 PDE 的特征线从某个位置进入当前网格点，但数值模板完全没有使用那个方向的数据，那么该数值格式就不可能正确传播信息。

离散信息传播速度必须能够跟上连续 PDE 的信息传播速度。

稳定性、相容性和收敛性

相容性

当

Δx,Δt→0\Delta x,\Delta t\to0

时，差分方程能够逼近原 PDE。

稳定性

数值误差不会在时间推进过程中被无限放大。

收敛性

当网格不断加密时：

unumerical→uexactu_{\mathrm{numerical}} \rightarrow u_{\mathrm{exact}}

对于适定线性问题，Lax 等价定理给出：

相容性+稳定性  ⟺  收敛性\boxed{ \text{相容性} + \text{稳定性} \iff \text{收敛性} }

因此仅仅拥有高阶截断误差并不够。

一个高阶但不稳定的格式没有实际价值。