NDSolve 的Projection方法

引言

当一个微分系统具有一定的结构,那么如果数值积分方法保留该结构将是有利的. 在某些情况下,求解具有约束解的微分方程是有用的. 投影方法通过采用数值积分方法的一个时间步骤,并且把近似解投影到真正的解演化的流形上来工作的.

NDSolve 包含一个微分代数求解器,它可能是适当的并且在 "微分代数方程的数值解" 中有更详细的描述.

有时候,方程的形式可能不会简化为 DAE 求解器所需要的形式. 此外,所谓的指标简化技术可以破坏某些微分系统可能拥有的结构特性,如辛结构 (symplecticity)(见 [HW96] 和 [HLW02]). 一个说明这一点的例子可以在 DAE 的文档中找到.

在这种情况下,往往可能求解一个微分系统,然后使用一个投影程序来确保约束仍被满足. 这就是方法 "Projection" 背后的思想.

如果微分系统是 -可逆的,那么一个对称投影过程可能更好(见 [H00]). 对称投影一般比投影的代价更加昂贵,还没有在 NDSolve 中得以实现.

不变量

考虑一个微分方程

其中 可以是一个向量或者是一个矩阵.

定义:一个非常量函数 被称为 (1) 的一个不变量,如果对于所有 .

这意味着 (2) 的每个解 都满足 .

不变量通常也被称作初积分,守恒量或者运动常数,意思完全一样.

流形

给定一个 维子流形,满足

给定一个微分方程 (3) 那么 意味着,对于所有 . 这是一个比不变性较弱的假设,因而 被称为弱不变量(见 [HLW02]).

投影算法

表示从一步数值积分器得到的解. 考虑一个约束条件下极小化问题将产生以下方程组(见 [AP91]、[HW96] 和 [HLW02]):

为了节省工作, 被近似为 . 把第一个关系代入 (4) 中的第二个关系将产生以下关于 的简化牛顿方案:

其中 .

第一个增量 具有 的量级,因此 (5) 通常快速收敛.

使用一个高阶积分方法的额外代价可以被投影步骤中的较少的牛顿迭代来抵消

对于在方法 "Projection" 中的终止标准,选项 "IterationSafetyFactor" 是与 NDSolve 所使用的工作精度中的最后一位中的一个单位相结合的.    

示例

加载一些功用程序包

线性不变量

定义一个模拟化学反应的刚性系统:
该系统具有一个线性不变量:
数值积分(见 [S86]),包括默认的 NDSolve 方法,通常可保持线性不变量不变 ,正如在不变量的误差的图形中可以看到的:

因此在这个例子中没有必要使用 "Projection" 方法.

某些数值方法准确保持二次型不变量(例如见 [C87]). 隐式中点规则,或者一阶段的 Gauss 隐式 RungeKutta 方法,就是这样的一种方法.    

谐振子

定义谐振子:
谐振子具有以下不变量:
使用方法 "ExplicitRungeKutta" 求解该系统. 不变量上的误差大致呈线性增长,这是耗散方法应用到哈密顿系统时的典型行为:
这里也使用方法 "ExplicitRungeKutta" 求解该系统,但是它也在每一步的末尾对解进行投影. 不变量误差的图形表明它在舍入误差内守恒:
由于该系统是哈密顿的(不变量是哈密顿),一个辛(symplectic) 积分器对这个问题表现很好,只给出一个小的有界误差:

摄动开普勒问题

这里加载一个被称为摄动开普勒问题的哈密顿系统,设置积分区间和采用的步长,以及定义在哈密顿公式中的位置变量:
该系统具有两个不变量,它们被定义为 HL

现在用一个实验说明在投影过程中使用所有可用的不变量的重要性(见 [HLW02]). 考虑使用以下方法获得的解:

Lotka Volterra

LotkaVolterra 系统的约束投影的一个例子在"求解 LotkaVolterra 方程的数值方法" 中给出.

欧拉方程

欧拉方程的约束投影的一个例子在 "刚体求解器" 中给出.    

选项总结

选项名
默认值
"Invariants"None指定微分系统的不变量
"IterationSafetyFactor"指定在不变量的迭代解中使用的安全因子
MaxIterationsAutomatic指定在不变量的迭代解中所使用的最大迭代次数
Method"StiffnessSwitching"指定对微分系统数值积分所使用的方法

方法 "Projection" 的选项.