形式陈述
画出一条近似轨道 z ( t ) 后,能否证明真解在整个时间窗口内都离它不远?答案需要同时管住方程缺陷、误差传播和解是否存在,不能先假设真解已经留在待证明的区域。
考虑实初值问题 理路 常微分方程 Ordinary differential equation · ODE 未知函数及其单一自变量导数组成的方程。
x ′ = f ( t , x ) , x ( a ) = ξ , a < T . 设 f 定义在开集 D ⊂ R × R d 上,f 与状态Jacobian D x f 联合连续。固定一个向量范数及其对数范数 理路 对数范数与瞬时增长界 Logarithmic norm · Matrix measure 从单位映射的单侧增率定义带符号的矩阵增长量,证明常用公式与最小指数传播界,并追踪换尺度的代价。 μ 。给定候选 z ∈ C 1 ( [ a , T ] , R d ) 、连续正函数 ρ ,要求闭管
(1) K = { ( t , y ) : a ≤ t ≤ T , ‖ y − z ( t ) ‖ ≤ ρ ( t ) } 是 D 的紧子集。初值集合 X 0 满足 ‖ ξ − z ( a ) ‖ ≤ r 0 ,其中 0 ≤ r 0 < ρ ( a ) 。
设已证明连续函数 m : [ a , T ] → R 、δ : [ a , T ] → [ 0 , ∞ ) 满足
(2) μ ( D x f ( t , y ) ) ≤ m ( t ) ( ( t , y ) ∈ K ) , ‖ z ′ ( t ) − f ( t , z ( t ) ) ‖ ≤ δ ( t ) ( a ≤ t ≤ T ) . m 可以为负。令
(3) M ( t ) = ∫ a t m ( s ) d s , R ( t ) = e M ( t ) r 0 + ∫ a t e M ( t ) − M ( s ) δ ( s ) d s . 若整个区间上都有可检查的严格余量
(4) R ( t ) < ρ ( t ) ( a ≤ t ≤ T ) , 则对每个 ξ ∈ X 0 ,初值问题存在唯一解到 T ,且
(5) ‖ x ( t ; ξ ) − z ( t ) ‖ ≤ R ( t ) , a ≤ t ≤ T . 结论同时覆盖全部允许初值,不只某个采样初值。式(2)第一行必须在管内所有状态成立,第二行则沿候选中心曲线计算。
常数输入时直接复算
若 m , δ 为常数,写 τ = t − a ,则
(6) R ( t ) = { e m τ r 0 + δ e m τ − 1 m , m ≠ 0 , r 0 + δ τ , m = 0. m < 0 时,持续缺陷的贡献趋向 δ / | m | ,而不会指数增长。式(6)只是误差半径的公式;还需验证式(4),以及输入界确实覆盖式(1)的整条管。
直觉
可以先在候选轨道周围放一根足够宽的“工作管”。在这根管里计算最坏增长量,然后用残差算出一根更细的“误差管”。若细管始终严格包含在工作管内部,真解就无法成为第一个走出工作管的轨道。
工作管承担的是空间条件,误差管承担的是结论。把两根管混成一根,容易出现循环论证:因为真解离得近,所以Jacobian有界;又因为Jacobian有界,所以真解离得近。严格余量将这两个方向接起来。
从向量误差到标量半径
局部存在唯一性 理路 Picard–Lindelöf 存在唯一性定理 Picard-Lindelof theorem · Cauchy-Lipschitz theorem 连续且对状态变量局部一致 Lipschitz 的向量场给出常微分方程初值问题的唯一局部解。 先给出每个初值附近的解。只在它还留在工作管时,令
e = x − z , r = z ′ − f ( t , z ) , B ( t ) = ∫ 0 1 D x f ( t , z + θ e ) d θ . 球是凸集,故整条线段 z + θ e 仍在同一截面的管内。沿线段积分导数给
(7) e ′ = B ( t ) e − r ( t ) , μ ( B ( t ) ) ≤ ∫ 0 1 μ ( D x f ( t , z + θ e ) ) d θ ≤ m ( t ) . 中间一步使用对数范数的凸性,而不是只在未知真解端点代入一个Jacobian。
对正 h ,由可微性
e ( t + h ) = ( I + h B ( t ) ) e ( t ) − h r ( t ) + o ( h ) . 取范数、减去 ‖ e ( t ) ‖ 并除以 h ,得到上右导数估计
(8) D + ‖ e ( t ) ‖ ≤ m ( t ) ‖ e ( t ) ‖ + δ ( t ) . 这里 D + u ( t ) = lim sup h ↓ 0 ( u ( t + h ) − u ( t ) ) / h ,只用于说明范数在零点或折角处也可处理。
在解存在的任何紧子区间上,e 为 C 1 ,其范数是Lipschitz函数,因而绝对连续 理路 绝对连续函数 Absolutely continuous function · Absolute continuity on an interval 把有限组总长度很小的区间送到总振幅很小的函数类,并满足 Lebesgue 版微积分基本定理。 。式(8)在普通导数存在处成立。用带符号微分不等式的积分因子 理路 Grönwall 不等式 Gronwall inequality · Grönwall lemma 将受自身积分控制的非负函数封闭为显式指数上界。 ,对 E = ‖ e ‖ 积分
( e − M E ) ′ ≤ e − M δ 便得 E ( t ) ≤ R ( t ) 。无需假设范数处处可微,也无需把负的 m 换成零。
首次出管与有限端点
若解在某个最早时刻 t ∗ 碰到工作管边界,则在之前的区间上已有 E ≤ R 。连续性给
ρ ( t ∗ ) = E ( t ∗ ) ≤ R ( t ∗ ) < ρ ( t ∗ ) , 矛盾。因此只要解存在,它就不会出管。
还须排除“解没出管,却在 T 前停止存在”。若最大存在区间的右端 b ≤ T 有限,轨道在紧集 K 内;f 在这里有统一上界,积分式使 x ( t ) 当 t ↑ b 时具有有限极限。极限点仍在 D 内,局部存在定理可从该点续接,且唯一性保证与原轨道相接。这与最大性矛盾。于是解确实延续到 T ,而不是仅在一个事先未知的短窗口满足式(5)。
例子与边界
一条非线性轨道,无需知道精确解
在 0 ≤ t ≤ 2 上考虑
(9) x ′ = − x − x 3 + 1 + t 10 , x ( 0 ) ∈ [ − 1 / 100 , 1 / 100 ] . 取候选 z ( t ) = t / 10 。代回方程得
r ( t ) = z ′ − f ( t , z ) = t 3 1000 , | r ( t ) | ≤ 1 125 . 标量Jacobian为 − 1 − 3 x 2 ≤ − 1 ,所以可取 m = − 1 。选择常数工作半径 ρ = 1 / 50 、初值误差 r 0 = 1 / 100 ,则
R ( t ) = 1 125 + 1 500 e − t ≤ 1 100 < 1 50 . 式(4)有至少 1 / 100 的余量。于是全部允许初值的解存在到二,并落在
t 10 − R ( t ) ≤ x ( t ) ≤ t 10 + R ( t ) . 特别是 x ( 2 ) 落在以 1 / 5 为中心、半径 1 / 125 + e − 2 / 500 的区间内。这份证书只用了多项式代入、全区间单调界和标量指数,没有把另一个数值解当作真解。
只检查中心Jacobian会漏掉什么
对 x ′ = x 2 ,取 z = 0 、初值 x ( 0 ) = 1 / 10 。中心的Jacobian为零,但真实解是
x ( t ) = 1 / 10 1 − t / 10 . 若误用 m = 0 、δ = 0 ,就会给出恒定误差 1 / 10 ,与任何正时间的真解矛盾。
正确地在工作管 | y | ≤ 1 / 4 上取 m = sup 2 y = 1 / 2 ,得到 R ( t ) = e t / 2 / 10 。对 T = 1 ,e 1 / 2 < 2 ,故 R < 1 / 5 < 1 / 4 ,认证成功。对 T = 2 ,R ( 2 ) = e / 10 > 1 / 4 ,同一证书失败;但真实终值仅为 1 / 8 ,仍在工作管内。失败只表示这组界没有闭合,不表示解不存在或必已逃出。
小残差不约束漏掉的初始条件
即使 r ≡ 0 ,两条不同初值的精确轨道也可能不同。式(3)的初值项不能省略。同样,只在若干时间点看到残差小,不能推出连续区间上的第二行式(2);采样间可能有窄峰,必须有覆盖整段的导数、区间或多项式界。
推论与应用
将有限数据交给复核器
实际证书可以把时间分成有限多个闭区间。每段保存候选多项式、工作半径、初值半径,以及可靠的残差和Jacobian界。用区间包含运算 理路 验证数值计算与区间算术 Verified numerics · Validated numerics · Interval arithmetic 用向外舍入的区间运算构造包含真实结果的可检查外包络,并分辨可靠包含与界的紧致程度。 在时间子区间与状态管的外包盒上评估,才能支持式(2)的全称量词;扩大到盒可以变松,但不能漏掉工作管。
例如采用固定加权无穷范数时,‖ y − z ( t ) ‖ w ≤ ρ 等价于逐分量 | y i − z i ( t ) | ≤ w i ρ 。对区间Jacobian,保留对角上端点、对非对角项取最大绝对值,再按 w j / w i 加权,可以直接构造 m 。指数和式(6)的机器值也要用可靠外包络,不能把最近舍入的小数当作严格余量。
候选若连续且逐段 C 1 ,在每段应用定理,将前段的末端半径作为下一段的初值界即可。若候选在连接点换中心产生跳跃 Δ z = z + − z − ,真实解仍连续,新的起始半径至少取
r new ≥ R − + ‖ Δ z ‖ . 不记录这次中心移动,就会丢失一部分误差。不同范数或权重的段间切换也要显式换算半径。
验收状态有明确含义
只有全部时间段连续覆盖 [ a , T ] ,每段式(1)–(4)均通过,才能返回“已认证到 T ”。某段管宽不足、区间界太松、指数上界超预算或时间覆盖有缺口时,应返回未决区间和最后认证时间。缩步、换中心、改善尺度或提高精度可能修复失败,但一般不保证搜索总会成功。
这与局部误差阶 理路 ODE 的局部误差与全局误差 Local and global error for ODE · ODE truncation error 区分从精确状态迈出一步的局部缺陷与累计数值轨道的全局误差,并用离散 Grönwall 估计误差传播。 及嵌入式误差估计 理路 自适应步长与嵌入式 Runge–Kutta Adaptive Runge-Kutta step size · Embedded Runge-Kutta method · Adaptive Runge-Kutta 用共享 stages 的嵌入式 Runge–Kutta 对估计单步误差,并通过尺度化接受、拒绝与受限控制器分配时间步。 配合使用:已有求解器负责提出便宜的候选轨道,本页负责检查一份全时段缺陷证书。候选怎样算出,不影响已经完整通过的包含证明;候选算得很精细,也不能代替该证明。
若模型读入过去状态,残差仍可作证据,但误差传播还须保留整个历史。延迟耗散比较 理路 Halanay 不等式与延迟误差底线 Halanay inequality · Forced Halanay bound · 时滞微分不等式 比较当前阻尼与过去窗口的最坏反馈,用完整上解证明指数衰减及持续扰动留下的可计算误差底线。 在全局适用的标量线性模型中给出持续缺陷底线,连续延迟残差检查 理路 连续 Euler 延迟方法与残差证书 Continuous Euler method for delay equations · Continuous delayed Euler residual certificate 以同一条连续折线读取离网格历史,并合并平移切点检查全时段残差,交付有限窗及耗散延迟模型的误差保证。 还须在平移后的旧节点处分段。这里的全局存在条件不能被只有当前中心曲线的局部Jacobian采样替代。
参考资料
Ordinary Differential Equations ,Uppsala大学课程存档的数值方法第13章章节稿 ,§13.1.3 Lemma 13.1.7、§13.1.5 Theorems 13.1.22–13.1.23,章内pp.13–14、27–29:Jacobian平均及缺陷控制。本页另明确紧管延拓与整窗口严格余量。
N. S. Nedialkov、K. R. Jackson、G. F. Corliss,“Validated Solutions of Initial Value Problems for Ordinary Differential Equations”,Applied Mathematics and Computation 105 (1999), pp.21–68。作者1997技术报告入口 ,§§3–5:先验证存在和全步包围,再细化端点的接口。目录中该报告文件名为 vsode.97.ps.gz。