形式陈述
考虑带初边值数据的演化型偏微分方程 公理库 偏微分方程 Partial differential equation · PDE 含未知多元函数及其偏导数的方程。
u t = L u + N ( u , t ) in Ω , B u = g ( t ) on ∂ Ω , u ( ⋅ , t 0 ) = u 0 . 线法先在固定空间网格 公理库 网格与网格剖分 Grid and mesh discretization · Structured grid · Unstructured mesh 用节点、单元及其几何映射离散计算区域,并以局部尺度、形状正则性和边界拟合质量描述一族网格。 上离散 L , N ,保留时间连续,得到半离散系统
U h ′ ( t ) = L h U h ( t ) + N h ( t , U h ( t ) ) + b h ( t ) = F h ( t , U h ( t ) ) . 向量 U h ( t ) ∈ R M h 表示节点值、单元平均或其他空间自由度;L h 继承局部空间耦合,通常是大规模稀疏矩阵 公理库 稀疏矩阵表示与运算 Sparse matrix · CSR matrix · CSC matrix 只存非零项及其位置,以稀疏模式组织矩阵运算、图结构、重排序与填充成本。 ,b h ( t ) 则承载消元后的边界贡献和外部源项。有限元等离散还可能产生质量矩阵 M h U h ′ = F h ( t , U h ) ;只有在 M h 可逆并且求解策略明确时,才可将它形式化为普通 ODE,奇异质量矩阵则可能给出微分代数系统。
第二阶段把半离散系统交给ODE 时间步进 公理库 常微分方程时间步进框架 Time-stepping method for ODE · ODE integrator · Numerical time integration 将常微分方程初值问题离散为时间网格上的数值状态与一步或多步更新,并分离流映射、误差层和退出状态。 。显式方法主要需要计算 F h ( t , U ) ,每步成本低,但稳定步长受 L h 谱尺度限制。隐式方法在一步中求解含 L h 与非线性 Jacobian 的代数系统,可以跨越快速衰减尺度,却要承担稀疏线性求解、Newton 收敛和内层容差。IMEX 方法把易引起刚性的线性扩散等项隐式处理,把反应或输运等项显式处理;拆分方式必须与实际稳定性和求解成本相匹配。
semi-discrete 分析研究连续时间系统 U h ( t ) 是否稳定并逼近 PDE;fully discrete 分析研究时间离散后的状态 U h n 。两层误差可由恒等式分开:
u ( t n ) − P h U h n = [ u ( t n ) − P h U h ( t n ) ] + P h [ U h ( t n ) − U h n ] . 第一项是空间半离散误差,通常由局部 h K 、空间逼近阶、边界几何和解的空间正则性控制;第二项是对维数为 M h 的 ODE 所产生的时间误差,依赖 Δ t n 、时间积分阶和稳定域。在额外的一致稳定条件下,总误差可能写成 O ( h p ) + O ( Δ t q ) ,但时间误差常数可随 h 变化,刚性边界条件还会造成 order reduction,不能把这条形式分解无条件当作独立阶数相加。
空间加密会移动半离散谱。对扩散算子,最负特征值的模通常按 h − 2 增长;对一阶输运算子,相关尺度通常按 h − 1 增长。因此同一个显式时间方法在更细网格上必须同步减小 Δ t ,即使时间方向的真解看起来同样平滑。CFL 限制正是空间离散谱与时间积分稳定域的接口,而不是纯粹的空间或时间属性。
随时间变化的边界数据必须进入 ODE 右端。若 Dirichlet 边界自由度被消去,它们通过 b h ( t ) 影响邻近内部节点;若保留为未知量,约束 U ∂ ( t ) = g ( t ) 是代数条件,不能让 ODE 求解器自行推进。使用边界 lifting 时还会出现 lifting 的时间导数,漏掉它会改变所求 PDE,而不是只降低数值阶。
直觉
线法把空间中无穷多个自由度压缩成一条高维状态轨道:每个网格自由度随时间形成一条“线”,所有线通过离散空间算子相互耦合。这样可以复用成熟的 ODE 误差控制和隐式求解框架,同时保留 PDE 网格产生的局部稀疏结构。
空间网格越细,状态向量越长,相邻自由度之间的快速交换尺度也越短。扩散解的宏观轮廓可能变化缓慢,网格上最高频的锯齿模态却能以 O ( h − 2 ) 的速率衰减;这正是半离散 PDE 会变成刚性 ODE 公理库 刚性常微分方程 Stiff ordinary differential equation · Stiff ODE · Stiffness 刻画快速衰减模态迫使显式方法使用远小于精度需求的步长这一方法相关现象,并说明隐式路线的能力与代价。 的来源。
把 ODE 求解器当黑箱会丢掉最有价值的信息。矩阵的非零模式、Jacobian–vector product、边界强迫和物理分量尺度决定每一步能否高效可靠地完成;只提供一个稠密数值 Jacobian 接口,会把原本局部的 PDE 耦合膨胀成不必要的存储与计算。
例子与边界
对一维齐次 Dirichlet 热方程
u t = κ u x x , 0 < x < 1 , u ( 0 , t ) = u ( 1 , t ) = 0 , 取 m 个内部节点和 h = 1 / ( m + 1 ) ,中心差分给出
U ′ ( t ) = L h U ( t ) , L h = κ h 2 [ − 2 1 1 − 2 ⋱ ⋱ ⋱ 1 1 − 2 ] . 其特征值为
λ j = − 4 κ h 2 sin 2 ( j π 2 ( m + 1 ) ) , j = 1 , … , m . 最负模态的大小趋近 4 κ h − 2 。前向 Euler 要求每个 Δ t λ j 落在负实轴稳定区间 [ − 2 , 0 ] ,因而
Δ t ≤ h 2 2 κ cos 2 ( π 2 ( m + 1 ) ) ∼ h 2 2 κ . 空间步长减半时,内部未知量约翻倍,显式稳定步数还约增加四倍。这个成本增长来自半离散谱,不是前向 Euler 的局部时间误差估计能够看见的。
后向 Euler 对同一系统给出
( I − Δ t L h ) U n + 1 = U n . 它移除了负实谱上的显式稳定上限,但每步要解稀疏线性系统。实现应求解该系统,而不是形成 ( I − Δ t L h ) − 1 ;若 Δ t 或 L h 不变,可以复用分解或预条件器。对非线性 N h ,Newton 矩阵还会包含 I − Δ t ( L h + J N h ) ,迭代失败必须触发拒步、缩步或求解器诊断。
若热方程改用时变边界 u ( 0 , t ) = g 0 ( t ) , u ( 1 , t ) = g 1 ( t ) ,只保留内部未知量时,半离散系统的第一和最后一行分别含有
κ h 2 g 0 ( t ) , κ h 2 g 1 ( t ) . 把这两项漏掉,ODE 求解器会稳定而精确地求解错误的齐次边界问题。若 g 在时间上有跳跃,时间积分器还必须在跳点重启或显式对齐步长;用跨越跳点的高阶插值会失去原有光滑性假设。
半离散能量衰减也不是全离散稳定性的结论。上面的 L h 满足 U T L h U ≤ 0 ,故连续时间半离散解不增大 Euclidean 能量;前向 Euler 若违反 Δ t = O ( h 2 ) 仍会放大最高频模态。反过来,一个 A-stable 时间方法也无法修复不一致的空间 stencil、错误边界或非物理解振荡。
自适应 ODE 求解器估计的是当前空间网格上的时间误差。空间分辨率不足时,把相对和绝对容差收得更紧只会更精确地积分错误的半离散模型。可靠实验应分别做空间加密与时间加密,并把迭代残差控制在两者以下;若同时改变 h 与 Δ t ,必须报告二者关系,不能用一次斜率拟合声称分别验证了空间阶和时间阶。
推论与应用
实用工作流先固定空间离散、自由度语义与边界处理,再根据 L h 的谱尺度和非线性结构选择显式、隐式或 IMEX 时间方法。大规模接口应公开稀疏模式、Jacobian 或 Jacobian–vector product、预条件器更新时机、分量尺度和实际边界强迫;这些信息决定 ODE 算法的成本与失败模式。
计算记录应区分空间自由度 M h 、局部尺度 h K 、时间步 Δ t n 、空间残差、时间误差估计和代数求解残差。终止状态也应说明是到达终点、时间步预算耗尽、稳定限制、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.