Skip to content

线法:先空间离散再时间积分

Method of lines · MOL · Semi-discretization in space

先把演化型 PDE 的空间算子离散成大规模 ODE,再依据谱尺度、稀疏 Jacobian 和误差目标选择时间积分器。

形式陈述

考虑带初边值数据的演化型偏微分方程

ut=Lu+N(u,t)in Ω,Bu=g(t)on Ω,u(,t0)=u0.

线法先在固定空间网格上离散 L,N,保留时间连续,得到半离散系统

Uh(t)=LhUh(t)+Nh(t,Uh(t))+bh(t)=Fh(t,Uh(t)).

向量 Uh(t)RMh 表示节点值、单元平均或其他空间自由度;Lh 继承局部空间耦合,通常是大规模稀疏矩阵bh(t) 则承载消元后的边界贡献和外部源项。有限元等离散还可能产生质量矩阵 MhUh=Fh(t,Uh);只有在 Mh 可逆并且求解策略明确时,才可将它形式化为普通 ODE,奇异质量矩阵则可能给出微分代数系统。

第二阶段把半离散系统交给ODE 时间步进。显式方法主要需要计算 Fh(t,U),每步成本低,但稳定步长受 Lh 谱尺度限制。隐式方法在一步中求解含 Lh 与非线性 Jacobian 的代数系统,可以跨越快速衰减尺度,却要承担稀疏线性求解、Newton 收敛和内层容差。IMEX 方法把易引起刚性的线性扩散等项隐式处理,把反应或输运等项显式处理;拆分方式必须与实际稳定性和求解成本相匹配。

semi-discrete 分析研究连续时间系统 Uh(t) 是否稳定并逼近 PDE;fully discrete 分析研究时间离散后的状态 Uhn。两层误差可由恒等式分开:

u(tn)PhUhn=[u(tn)PhUh(tn)]+Ph[Uh(tn)Uhn].

第一项是空间半离散误差,通常由局部 hK、空间逼近阶、边界几何和解的空间正则性控制;第二项是对维数为 Mh 的 ODE 所产生的时间误差,依赖 Δtn、时间积分阶和稳定域。在额外的一致稳定条件下,总误差可能写成 O(hp)+O(Δtq),但时间误差常数可随 h 变化,刚性边界条件还会造成 order reduction,不能把这条形式分解无条件当作独立阶数相加。

空间加密会移动半离散谱。对扩散算子,最负特征值的模通常按 h2 增长;对一阶输运算子,相关尺度通常按 h1 增长。因此同一个显式时间方法在更细网格上必须同步减小 Δt,即使时间方向的真解看起来同样平滑。CFL 限制正是空间离散谱与时间积分稳定域的接口,而不是纯粹的空间或时间属性。

随时间变化的边界数据必须进入 ODE 右端。若 Dirichlet 边界自由度被消去,它们通过 bh(t) 影响邻近内部节点;若保留为未知量,约束 U(t)=g(t) 是代数条件,不能让 ODE 求解器自行推进。使用边界 lifting 时还会出现 lifting 的时间导数,漏掉它会改变所求 PDE,而不是只降低数值阶。

直觉

线法把空间中无穷多个自由度压缩成一条高维状态轨道:每个网格自由度随时间形成一条“线”,所有线通过离散空间算子相互耦合。这样可以复用成熟的 ODE 误差控制和隐式求解框架,同时保留 PDE 网格产生的局部稀疏结构。

空间网格越细,状态向量越长,相邻自由度之间的快速交换尺度也越短。扩散解的宏观轮廓可能变化缓慢,网格上最高频的锯齿模态却能以 O(h2) 的速率衰减;这正是半离散 PDE 会变成刚性 ODE的来源。

把 ODE 求解器当黑箱会丢掉最有价值的信息。矩阵的非零模式、Jacobian–vector product、边界强迫和物理分量尺度决定每一步能否高效可靠地完成;只提供一个稠密数值 Jacobian 接口,会把原本局部的 PDE 耦合膨胀成不必要的存储与计算。

例子与边界

对一维齐次 Dirichlet 热方程

ut=κuxx,0<x<1,u(0,t)=u(1,t)=0,

m 个内部节点和 h=1/(m+1),中心差分给出

U(t)=LhU(t),Lh=κh2[2112112].

其特征值为

λj=4κh2sin2(jπ2(m+1)),j=1,,m.

最负模态的大小趋近 4κh2。前向 Euler 要求每个 Δtλj 落在负实轴稳定区间 [2,0],因而

Δth22κcos2(π2(m+1))h22κ.

空间步长减半时,内部未知量约翻倍,显式稳定步数还约增加四倍。这个成本增长来自半离散谱,不是前向 Euler 的局部时间误差估计能够看见的。

后向 Euler 对同一系统给出

(IΔtLh)Un+1=Un.

它移除了负实谱上的显式稳定上限,但每步要解稀疏线性系统。实现应求解该系统,而不是形成 (IΔtLh)1;若 ΔtLh 不变,可以复用分解或预条件器。对非线性 Nh,Newton 矩阵还会包含 IΔt(Lh+JNh),迭代失败必须触发拒步、缩步或求解器诊断。

若热方程改用时变边界 u(0,t)=g0(t),u(1,t)=g1(t),只保留内部未知量时,半离散系统的第一和最后一行分别含有

κh2g0(t),κh2g1(t).

把这两项漏掉,ODE 求解器会稳定而精确地求解错误的齐次边界问题。若 g 在时间上有跳跃,时间积分器还必须在跳点重启或显式对齐步长;用跨越跳点的高阶插值会失去原有光滑性假设。

半离散能量衰减也不是全离散稳定性的结论。上面的 Lh 满足 UTLhU0,故连续时间半离散解不增大 Euclidean 能量;前向 Euler 若违反 Δt=O(h2) 仍会放大最高频模态。反过来,一个 A-stable 时间方法也无法修复不一致的空间 stencil、错误边界或非物理解振荡。

自适应 ODE 求解器估计的是当前空间网格上的时间误差。空间分辨率不足时,把相对和绝对容差收得更紧只会更精确地积分错误的半离散模型。可靠实验应分别做空间加密与时间加密,并把迭代残差控制在两者以下;若同时改变 hΔt,必须报告二者关系,不能用一次斜率拟合声称分别验证了空间阶和时间阶。

推论与应用

实用工作流先固定空间离散、自由度语义与边界处理,再根据 Lh 的谱尺度和非线性结构选择显式、隐式或 IMEX 时间方法。大规模接口应公开稀疏模式、Jacobian 或 Jacobian–vector product、预条件器更新时机、分量尺度和实际边界强迫;这些信息决定 ODE 算法的成本与失败模式。

计算记录应区分空间自由度 Mh、局部尺度 hK、时间步 Δtn、空间残差、时间误差估计和代数求解残差。终止状态也应说明是到达终点、时间步预算耗尽、稳定限制、Newton 失败、线性求解失败还是非有限状态。只有把这些接口分开,才能判断下一步应加密空间、缩小时间步、改用隐式方法,还是修正边界与模型。

参考资料
  • Willem Hundsdorfer and Jan Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer, 2003, method of lines and IMEX integration.
  • Randall J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007, Chs. 9–10.
  • Ernst Hairer and Gerhard Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd rev. ed., Springer, 1996, stiff systems and implicit solvers.
  • MIT OpenCourseWare, 18.336, Numerical Methods for Partial Differential Equations, time-dependent PDE notes.