2. 离散化方法选型框架
非线性PDE的离散化不是"挑一个顺手的格式"，而是由方程结构、解的 regularity、守恒性质和计算预算共同约束的决策问题。本节建立一个可操作的选型框架：先看方程属于守恒型还是非守恒型，再看非线性项出现在通量、源项还是系数中，最后用稳定性与收敛性要求筛掉不合适的候选格式。
四类方法的本质差异
有限差分（FDM）直接在网格节点上逼近导数，实现最简单，对结构化网格上的光滑解效率很高。但它的守恒性依赖差分格式的对称构造，遇到激波或陡梯度时容易产生非物理振荡。有限体积（FVM）以控制体通量平衡为出发点，局部守恒是格式的内禀性质，特别适合对流主导问题和守恒律，非线性通量的处理（如Riemann求解器、通量分裂）是其核心工作。有限元（FEM）基于弱形式，对复杂几何和变系数问题适应性强，高阶单元可达到高精度，但非线性项会使刚度矩阵变为解相关，需要处理非线性迭代与质量矩阵求逆。谱方法用全局基函数展开，对光滑解具有指数收敛速度，但对间断和局部陡变极其敏感，且非线性项的卷积或配点计算会引入混叠误差。
非线性项如何改变格式选择
线性PDE中，离散矩阵是常数矩阵，稳定性分析可以一次性完成。非线性项打破了这个前提，带来三类具体影响。第一，通量非线性（如Burgers方程中的 u²/2）使特征速度依赖解本身，迎风方向随解变化，FVM中必须用Riemann求解器或通量分裂来确定界面通量，FDM中则需采用守恒型差分或加入人工粘性。第二，系数非线性（如扩散系数 D(u)）使刚度矩阵变为解相关，FEM组装后得到非线性代数方程组，必须用Newton法或Picard迭代求解，每次迭代重新组装矩阵。第三，源项非线性可能引入刚性，例如反应项为指数型时，显式格式的时间步会被反应时间尺度严重限制，此时IMEX格式（对流显式、反应隐式）是常见折中。
谱方法还要额外面对混叠误差：非线性乘积在配点上的离散计算会把高频模态折叠到低频，表现为虚假的能量增长。工程上的标准做法是3/2去混叠规则——将模态数扩展1.5倍再截断，代价是每次非线性评估的计算量增加约2.25倍。
稳定性与收敛性的实际权衡
Lax等价定理告诉我们：相容格式收敛当且仅当稳定。对非线性问题，稳定性不再是简单的von Neumann条件，而需要能量方法或熵条件来保证。实践中常见的权衡是：显式格式每步便宜但受CFL限制，非线性使最大特征速度随时间变化，固定步长要么保守浪费、要么越界失稳；隐式格式无条件稳定但每步需要非线性迭代，且大时间步会引入过多数值耗散，抹平陡梯度。一个实用的折中是自适应时间步：每步根据当前解的|u|max重新计算CFL允许的Δt，并设置增长上限（如每步最多放大1.2倍）以避免突变。
实战：一维非线性对流扩散方程的离散策略
考虑方程 ut + u·ux = ν·uxx，x∈[0,1]，周期边界，ν=0.01，初始条件 u(x,0)=sin(2πx)。这是Burgers方程的标准形式，非线性对流项会产生陡梯度，扩散项负责耗散。下面给出一套可执行的选型与实现方案。
选型判断：方程是守恒型，非线性在通量中，解在后期会出现陡梯度但不会形成真正的激波（ν>0使激波厚度有限）。FVM能保证守恒且通量处理成熟，是首选；若追求更高精度，可用FEM的间断Galerkin（DG）格式。这里选择FVM配合Rusanov通量，空间二阶MUSCL重构，时间上用显式SSP-RK2，并配合自适应CFL。
import numpy as np
N = 200; L = 1.0; dx = L/N
nu = 0.01; CFL = 0.4
x = np.linspace(0, L, N, endpoint=False)
u = np.sin(2*np.pi*x)
def flux_rusanov(uL, uR):
# f(u)=u^2/2, 最大特征速度 max(|uL|,|uR|)
a = max(abs(uL), abs(uR))
return 0.25*(uL**2 + uR**2) - 0.5*a*(uR - uL)
def rhs(u):
# MUSCL重构（minmod限制器）
uL = np.roll(u, 1); uR = np.roll(u, -1)
slope = np.where((u-uL)*(uR-u) > 0,
np.sign(u-uL)*np.minimum(np.minimum(abs(u-uL), abs(uR-u)), 0.5*abs(uR-uL)),
0.0)
u_minus = u - 0.5*slope # 界面左侧值
u_plus = u + 0.5*slope # 界面右侧值
F = flux_rusanov(u_plus, np.roll(u_minus, -1))
# 扩散项中心差分
diff = nu*(np.roll(u,1) - 2*u + np.roll(u,-1))/dx**2
return -(F - np.roll(F,1))/dx + diff
t = 0.0; t_end = 0.5
while t < t_end:
umax = np.max(np.abs(u))
dt = CFL*dx/max(umax, 1e-8)
dt = min(dt, t_end - t)
k1 = rhs(u)
k2 = rhs(u + dt*k1)
u = u + 0.5*dt*(k1 + k2)
t += dt
关键参数与判读：CFL=0.4给显式RK2留出安全裕度；MUSCL重构在陡梯度处自动降为迎风以保证无振荡；Rusanov通量虽然耗散偏大，但对Burgers这类标量问题足够稳健。运行后应观察：t=0.5时解在x≈0.5附近形成陡峭但光滑的过渡层，峰值从1.0衰减到约0.6（扩散耗散），且全程无伪振荡。若把CFL提到0.9，可能在陡梯度形成阶段出现振荡甚至发散，这正是非线性使CFL约束"动态化"的直接体现。若把ν改为0，则需改用更精确的Riemann求解器并考虑熵修正，否则激波位置和强度会出现偏差。