NDSolve 的 IDA 方法

IDA 程序包是在劳伦斯利弗莫尔国家实验室的应用科学计算的中心研制出来的 SUNDIALS（SUite of Nonlinear and DIfferential/ALgebraic equation Solvers）的一部分. 如 IDA 用户指南中所描述的 [HT99]，“IDA 是微分-代数方程组（DAE）的初值问题的通用求解器. 名称 IDA 代表隐式微分-代数求解器（Implicit Differential-Algebraic solver）. IDA 是以 DASPK 为基础的...” DASPK [BHP94]，[1] 是用来求解大规模微分-代数系统的 Fortran 代码.

Wolfram 语言为 IDA 程序包供给了一个界面，因而用户无需以 C 编写函数来计算残差和编译程序，而是由 Wolfram 语言从用户输入给 NDSolve 的方程来自动产生函数.

IDA 求解指标为1的 DAE 系统的形式如下：

DAE 的指标是对DAE进行求导以获得ODE系统所需的求导次数. 有关微分代数方程及其来源的更多信息，请查阅教程微分代数方程的数值解.

IDA 使用一个以可变步长形式实现的，阶数为1到5的向后微分公式（BDF）（Backward Differentiation Formula） 来求解系统 (2) . 在时间 上的阶数为 的 BDF 由下面公式给出

系数 取决于阶数 和过去的步长. 把 BDF 应用到 DAE (3) 给出一个非线性方程组来求解：

方程组的解由牛顿-类型的方法实现，通常使用对雅可比的一个近似

“IDA 最显著的特点是，在每个时间步进的内部非线性系统的解中，它提供了牛顿/直接方法或者一个不精确牛顿/Krylov（迭代）方法的一个选择.” [HT99] 在 Wolfram 语言中，用户可以使用方法选项选用这些求解器，或者 Wolfram 语言默认的LinearSolve，它会自动切换到对于大型问题的直接稀疏求解器.

在解的每个步进上，IDA 计算局部截断误差的一个估计 ，并且选择步长和阶数使得加权范数 Norm[En/wn] 小于1. 的第 部分 由下式给出

值 prec 和 acc 采用 PrecisionGoal->prec 和 AccuracyGoal->acc 的 NDSolve 设置.

因为 IDA 提供了很大的灵活度，尤其是在求解非线性方程方面，有大量的方法选项供用户来控制这一灵活性的实现. 用户可以通过给出 NDSolve 选项 Method->{IDA,ida method options} 的方式对IDA 使用这些方法选项.

IDA 方法的选项与 NDSolve` 语境中的符号 IDA 相关联：
IDA 方法选项名
默认值
"ImplicitSolver""Newton"如何求解隐式方程
"MaxDifferenceOrder"5采用的最大阶数 BDF

IDA 方法选项.

当使用从 NDSolve 返回的 InterpolatingFunction 对象计算得到的中间值的严格准确性很重要的时候，用户将要使用 NDSolve 方法选项设置 InterpolationOrder->All，该设置基于方法的阶数使用内插法，有时候称为稠密输出，来表示在时间步进之间的解. 默认情况下，NDSolve 存储一个足够用图形来表示解的最小量的数据. 保持小的数据量为更复杂的解节省了内存和时间.

作为一个强调最小输出和完全方法内插阶数的不同的例子，考虑跟踪一个从一个简单线性方程的解推导出来的量 ，对这个简单线性方程，精确解可以用 DSolve 计算.

这里定义函数 f，并且把该量作为基于解 x[t] 和 y[t] 的时间的函数：
这里定义具有初始条件的线性方程：
f 作为一个时间函数的精确值可以使用 DSolve 符号计算：

精确解将与由稠密输出和非稠密输出计算所得的解比较.

跟踪该数量的一个简单方式是创建一个函数，可以从微分方程的数值解中把它推导出来：
它也可以使用稠密输出求解：
该图显示了两个计算解的相对误差. 在时间步长上的计算解以黑点表示. 默认输出错误以灰色显示，密集输出错误以黑色显示:

从该图形中，很明显看出当时间步长变大时，默认解输出在时间步进之间具有一个很大的误差. 稠密输出解表示由求解器计算得到的解，甚至在时间步进之间. 不过，记住稠密输出解将使用更多空间.

当用户想要从解推导出来的量复杂的时候，用户可以通过把它作为一个代数量给出，确保它保持在局部容差范围内，使得求解器可以把误差保持在可控的范围内.

这里用一个代数方程增加了一个应变量，该代数方程把应变量设置为等于感兴趣的量：
这里做出图形，比较 DAE 解的相对误差. IDA 的时间步长用蓝点显示：

DAE 解需要更多步骤来控制感兴趣的量的误差，但与使用密集输出相比，所用的内存仍然更少.

这将比较默认、密集输出和 DAE 解的大小：

本文档的剩余部分将重点介绍 "ImplicitSolver" 选项的两个可能设置的子选项，"Newton" 和 "GMRES". 使用 "Newton"，可以计算雅可比或近似值，并求解线性系统以找到牛顿阶跃. 而 "GMRES" 是一种迭代非线性求解器方法，它并不计算整个雅可比行列式，而是为每个迭代步骤计算方向导数.

"Newton" 方法具有一个方法选项，"LinearSolveMethod"，用户可以使用它来告诉 Wolfram 语言如何求解要求的线性方程以找到牛顿步. 该选项有几个可能的值.

Automatic这是默认的：基于雅可比的带宽，在使用 Wolfram 语言 LinearSolve 和 Band 方法间自动切换；对于具有较大带宽的系统，这将对具有稀疏雅可比的大型系统自动切换到一个直接稀疏求解器
"Band"使用 IDA 带方法（参阅 IDA 用户手册以获取更多信息）
"Dense"使用 IDA 稠密方法（参阅 IDA 用户手册以获取更多信息）

"LinearSolveMethod" 选项的可能设置.

"GMRES" 方法可能实质上更快，但是通常使用起来棘手，因为要真正有效，它通常需要一个预设条件，它可以通过一个方法选项提供. 我们还有一些其它的方法选项来控制 Krylov 子空间的过程. 要使用这些，请参阅 IDA 用户指南 [HT99].

GMRES 方法选项名
默认值
"Preconditioner"Automatic返回作为预设条件的另一个函数的一个 Wolfram 语言函数
"OrthogonalizationType""ModifiedGramSchmidt"这也可以是 "ClassicalGramSchmidt"（参阅 IDA 用户指南中的变量 gstype）
"MaxKrylovSubspaceDimension"Automatic最大子空间维度（参阅 IDA 用户指南中的变量 maxl）
"MaxKrylovRestarts"Automatic重启的最大数目（参阅 IDA 用户指南中的变量 maxrs）

"GMRES" 方法选项.

作为一个实例问题，考虑一个二维 Burgers 方程

通常可以使用常微分方程求解器求解此问题，但是在本示例中，使用 DAE 求解器可以实现两个目的. 首先，边界条件被强制为代数条件. 其次，NDSolve 被迫使用代数项来进行保守差分. 为了比较，已知的精确解将用于初始条件和边界条件. 下面的计算将显示两种方法之间的差异.

这里定义了一个满足 Burger 方程的函数：
这里对单元正方形定义与精确解一致的初始和边界条件：
这里定义了该微分方程：
这里设置扩散常数 为某一个值，以使用户可以在一个合理的时间内找到解，并且显示在 的解的图形：

从图形中可以看出，使用 ，这里有一个相当陡峭的前面（front）. 它以恒定速度移动.

这里使用 NDSolve 和 IDA 方法的默认设置求解该问题，"MethodOfLines" 的 "DifferentiateBoundaryConditions" 选项除外，因为这导致 NDSolve 将边界条件作为代数条件处理：

由于有一个可比较的精确解，总体解的准确性也可以比较.

这里定义一个函数，以寻找在所有时间步进上精确解和计算所得解在网格点的最大偏差：
这里对计算所得的解计算最大误差：

接下来，将对 IDA 方法的选项的不同设置进行比较. 为了强调选项设置，将需要定义一个函数来对解的计算计时，并且给出最大误差.

这里定义一个用来比较不同 IDA 选项设置的函数：
没有选项使用前面的默认值：
这里使用 "Band" 方法：

"Band" 方法不是非常有效，因为雅可比的带宽相对较大，部分原因是由于采用了四阶导数，部分原因是因为在靠近边界处采用了单边模板（stencil）增加了顶部和底部的宽度. 用户可以明确指定带宽.

这里使用 "Band" 方法，带宽的设置只包括在一个方向的差异的模板（stencil）：

虽然求解时间更短并且误差不大，但总的时间步长更大. 如果问题更加严峻，则迭代可能不会收敛，因为它缺少来自另一个方向的信息. 理想情况下，带宽不应将整个维度的信息消除.

这里计算在 和 方向使用的网格 ，并且显示采用的点的数目：
这里使用 "Band" 方法时，宽度的设置至少包括在两个方向上的模板（stencil）的一部分：

使用更合适的带宽设置，求解时间仍比选择默认带宽稍快，但比使用默认线性求解器慢. "Band" 方法有时会对二维问题有效，但通常对一维问题最有效.

在没有预设条件的情况下，这里使用 "GMRES" 隐式求解器计算解：

与默认方法相比，使用不带前置条件的 "GMRES" 方法会导致计算时间更长，步骤更多. 因此，不建议在没有前置条件的情况下使用此方法. 然而，找到一个好的预处理器通常并不容易. 对于此例，将使用对角预处理器.

"Preconditioner" 选项的设置应该是函数 f，它接受由 NDSolve 给出的四个自变量，使得 f[t,x,x′,c] 返回另一个函数将预处理器应用于残差矢量. （参阅 IDA 用户指南 [HT99] 了解如何使用该预处理器的更多信息.) 自变量 t、x、x′、c 是当前时间、解向量、解导数向量和上面公式 (4) 中的常数 c. 例如，如果可以确定一个过程来生成适当的前置条件矩阵 作为这些自变量的函数，则可以使用

"Preconditioner"->Function[{t,x,xp,c},LinearSolve[P[t,x,xp,c]]]

来产生一个 LinearSolveFunction 对象，它实际上将对预设条件矩阵 求逆. 通常，每次当预设条件函数建立的时候，它就被几次应用到残差向量中，所以使用某种因式分解，比如在 LinearSolveFunction 中包含的，是一个好办法.

对于对角情况，逆可以简单通过使用倒数实现. 建立一个对角预设条件最难的部分是要时时记住边界上的值不应该被它影响.

这里寻找微分矩阵的对角元素，以便计算预条件：
这里获取位置信息，其中在边界处满足一个简单的代数条件的元素位于压平（flattened）的解向量中：
这里定义了函数，来设置调用来获取预设条件的有效逆的函数. 对于对角线的情况下，逆的完成只需要求倒数：
这里使用具有预设条件的 "GMRES" 方法来计算解：

因此，即使使用一个粗略的预设条件，"GMRES" 方法计算解的速度也比使用直接稀疏求解器快.

对于高阶时间导数或者偏微分方程组的偏微分方程离散，用户可能需要查看相应的 NDSolve`StateData 对象来确定变量是如何进行排序的，以便用户可以正确地获取预条件的结构形式.