NDSolve 方法的插件框架

引言调用其他方法
经典 Runge–KuttaAdams 方法
Crank–Nicolson

引言

为 NDSolve 所设置的控制机制使用户能够自定义数值积分算法并将其作为规范用于 NDSolve 的选项 Method.

NDSolve 以面向对象的方式对它的数值算法以及所需信息进行访问. 在数值积分的每一步，NDSolve 以一种能够根据需要保存私有数据的形式将方法保存起来.

AlgorithmIdentifier [data]一个算法对象，该对象包含某一常微分方程的数值积分算法可能需要使用的任何 data；数据为算法所专有；AlgorithmIdentifier 应是一个 Wolfram 系统符号，并且算法通过使用选项 Method->AlgorithmIdentifier 从 NDSolve 访问

用于 NDSolve 的方法数据的结构.

NDSolve 不直接对与算法相关的数据进行访问，因此用户可以以任何方便或高效使用的形式保留所需的信息. 可能保存在私有数据中的算法和信息只能通过该算法对象的方法函数才能访问. 方法算法对象和 NDSolve 的重要性是通过 "Step" 函数实现的.

obj["Step"[inputs]]使用数值算法尝试单个时间步骤

NDSolve 要求使用的 "Step" 函数.

下面是方法的重要属性和它应该被定义的 "Step" 函数组成的列表.

obj["StepInputs"]步骤函数应该给出的输入参数列表
obj["StepOutput"]步骤函数返回的输出列表
obj["StepMode"]设置步骤模式
obj["DifferenceOrder"]设置算法的当前渐进差分阶数

使用 NDSolve 中的方法算法对象 obj，应该定义的属性.

下面的表格给出可以用于步骤输入和/或输出的规范.

简称
全称
分量
"F""Function"给出 ODE 系统的右侧和 DAE 系统的残差的函数
"T""Time"当前时间
"H""TimeStep"时间步骤
"X""DependentVariables"除了离散变量的应变量的当前值
"XI""DependentVariablesIncrement"，使得在步骤后 给出应变量的数值；如果 给出 $Failed，则拒绝该步骤并且重新计算（返回的时间步骤应该更小）
"XP""TemporalDerivatives"应变量的时间导数
"D""DiscreteVariables"采用数值连续范围的离散变量
"ID""IndexedDiscreteVariables"采用数值有限集合的离散变量
"P""Parameters"计算灵敏度需要的其他参数
"SP""SensitivityParameters"计算灵敏度参数
"SD""SolutionData"解数据列表
"Obj""MethodObject"用于下一个步骤的 obj 的新设置
"State""StateData"NDSolve`StateData 对象

步骤输入和输出规范.

该解数据是与不同变量类型相关的分量组成的数值列表. 激活方向下的解数据可以使用 state@"SolutionData"["Active"] 从 NDSolve`StateData 对象 state 得到. 上面列出的许多规范直接对应于解数据分量.

对于输入，函数规格应该使用给定的参数. 这些参数应该是解数据分量. All 可用于表明 NDSolve 所使用的所有参数都应该有.

对于输出，如果该方法仅返回某些情况下的某些数值（例如，该方法对象没有再每个步骤改变），那么您可以对任何输出使用 Null，除了增量 ，以表明不应该使用它.

经典 Runge–Kutta

这个例子说明如何设置和访问一个简单的数值算法.

这里定义一个方法函数，使用经典四阶 Runge–Kutta 方法朝向一个常微分方程的积分取一个单一步进. 由于该方法非常简单，不必保存任何私有数据：
下面对步骤函数指定了 NDSolve 的输入和输出. ___ 模板用于使 CRK4 方法对象的所有实例可以使用：
这是使用完整名称的等价规范：
这里定义一个方法函数，使得 NDSolve 可以获得用于该方法的适当差分阶：
这里为步进模式定义一个方法函数，使得 NDSolve 知道如何控制时间步. 这种算法没有任何步进控制，因此步进模式定义为 Fixed：
这里用固定步长对简谐振子方程进行积分：

一般地，与允许步长随积分局部困难程度而变化相比，使用固定步长的效率要降低. 现代显式 Runge–Kutta 方法 （在 NDSolve 中通过 Method->"ExplicitRungeKutta" 访问）具有一种所谓嵌入式误差估计程序，使得非常高效地确定大致步长成为可能. 另一种方法是使用内置的步进控制程序，这种程序使用的是外推法. 方法 "DoubleStep" 使用的是一种基于时间步积分的外推法，其积分的单一步的步长为 ，两步步长为 . 方法 "Extrapolation" 进行的是一种更为复杂的外推，随积分进行自动地修改外推度，但一般与差分阶1和2的基本方法一起使用.

这里利用步进由 "DoubleStep" 方法控制的经典四阶 Runge–Kutta 方法对简谐振子进行积分：
这里绘出一幅图形来比较步进终点计算得到的解的误差. 使用 "DoubleStep" 方法得到的误差用蓝色表示：

固定步长的整体误差较小，主要是因为步长非常小，所需步数超过另一种方法的3倍. 对于一个局部解结构变化较为显著的问题，这种差距可能会更大.

一种刚度检测的工具在 "NDSolve 的 DoubleStep 方法" 一节中介绍.

对于更加复杂的方法，为所用方法设置一些数据可能是必要或更加有效的. 当 NDSolve 首次使用一种特别的数值算法时，将调用一个初始化函数. 您可以定义初始化规则，为您的方法设置适当的数据.

NDSolve`InitializeMethod[Algorithm Identifier,stepmode,sd,f,state,Algorithm Options]
NDSolve 在首次使用一种用于特定积分的算法时，进行初始化所计算的表达式，根据方法是否在一个步进控制程序或另一种方法框架内被调用，stepmode 是 Automatic 或 Fixed；sd 是当前解数据，state 是由 NDSolve 所用的 NDSolveState 对象，Algorithm Options 是一个包含特别指定使用某种特定算法的任何选项的列表，例如 Method->{Algorithm Identifier,opts} 中的 {opts}
Algorithm Identifier/:NDSolve`InitializeMethod[Algorithm Identifier,stepmode_,sd_,f_NumericalFunction,state_NDSolveState,{opts___?OptionQ}]:=initialization
使得规则与算法相联系的初始化定义，initialization 应该返回一个 Algorithm Identifier[data] 形式的算法对象

从 NDSolve 中初始化算法.

NDSolve`InitializeMethod 是受保护的，因此如果想在其上附加规则，您需要首先将其取消保护. 最好是使规则与您的方法相关联. 要做到这一点，一个简洁的方法是使用上面提到的 TagSet 来进行初始化定义.

举例来说，假设您想要重新定义先前提到的 Runge–Kutta 方法，使用精确度适当的数字值而不是精确系数2、1/2 和1/6，从而使得算法稍微有所加快.

这里定义一个方法函数，使用经典四阶 Runge–Kutta 方法向着常微分方程积分的方向取一单一步进，使用保存的数值作为所需的系数：
这里定义一个规则，使用后来将要用到的数据初始化算法对象：

保存数字的数值对于使用 "DoubleStep" 的较长积分来说，速度提高了5%到10%.

Crank–Nicolson

下面是如何设置隐式方法的实例. 若要求解常微分方法组 ，尺寸为 的时间步骤的方案是 ，其中 和 . 一般来说，对于非线性 , 方程需要使用牛顿迭代求解.

不适用每步更新雅可比的真正的牛顿方法，只计算一次雅可比，然后重复使用它的 LU 分解来求解线性方程. 只要时间步骤没有太大，就通常会收敛并且完全牛顿方法和除了第一次迭代以外的每次迭代都会显著便宜.

若要避免牛顿迭代不收敛的可能问题，将使用控制最大迭代次数和收敛容差的方法选项.

定义 Crank–Nicolson 方法的选项和方法初始化：
定义步骤函数：
对于 Crank–Nicolson 方法，下面指定了NDSolve 的输入、输出、差分阶数和步骤模式：
下面显示方法可以使用选项设置如何调用：

对于线性问题，解由第一次迭代给出. 但是，对于非线性问题，收敛可能需要更长时间.

这没有成功，因为非线性问题不能在如此大的时间步长下，在一次迭代内收敛：
使用更多次迭代，解是成功的：

由于 Crank–Nicolson 方法通常用来求解偏微分方程，下面是一个求解波动方程的例子. 由于这是一个线性方程，收敛出现在1次迭代中，因此该方法速度相对快.

这里使用 Crank–Nicolson方法来求解具有周期性边界条件的波动方程：
解的密度图：

调用其他方法

NDSolve 框架的一个强项来自以各种方式将方法嵌套的功能.

NDSolve`InitializeSubmethod[caller,submspec,stepmode,sd,f,state]在方法 caller 的初始化中，初始化由 Method->submspec 指定的子方法；stepmode、sd、f 和 state 是 NDSolve`InitializeMethod 中的
NDSolve`InvokeMethod[submobj,f,h,sd,state]使用方法对象 submobj 调用子方法，其中函数为 f、时间步长为 h，解数据为 sd，以及 NDSolve`StateData 对象 state，返回一个列表 {hnew,sdnew,submobjnew}，其中 hnew 是新的时间步骤，sdnew 是设置为 Null 的未改变的组件的步骤结束时的解数据，而 submobjnew 是新子方法对象

使用子方法.

在把方法名称从给出的任何子选项分隔出来后，NDSolve`InitializeSubmethod 实际上调用 NDSolve`InitializeMethod. 如果成功，它返回您后来在 NDSolve`InvokeMethod 中使用的方法对象.

NDSolve`InvokeMethod 实际上对子方法构建 "Step" 函数. 由于函数的参数由子方法的 "StepInput" 规范指定，对 NDSolve`InvokeMethod 给出的函数 f 应该具有由 NDSolve 使用的全部参数，并且可以在 "Function"[All] 的调用器 "StepInput" 规范中指定. 注意，对于内置方法，可能没有显式的 "Step" 函数和 NDSolve`InvokeMethod 使用有效的内部代码处理这些. 由 NDSolve`InvokeMethod 返回的方法总数从子方法实际输出中构建，以匹配 {"TimeStep","SolutionData","MethodObject"} 的一个 "StepOutput".

下面的例子简单地调用一个子方法，以提供步骤和差分阶数的简单监控.

定义方法选项和方法的初始化：
该方法的步骤函数简单地调用子方法，并且将监控函数应用于结果和差分阶数：
设置输入和输出.
从子方法推导所有其他属性：
使用该方法：

对于更复杂的积分，该方法可以与 Reap 一起使用，用于收集数据.

收集步长的对数和差分阶数，作为时间的函数，用于对 Van der Pol 振荡器求积分，并且绘制该解：
显示步长的对数和突发改变附近步骤的差分阶数：

Adams 方法

就 NDSolve 框架而言，写一个自动控制步进的算法并非更难，但误差估计和确定适当步长的要求通常使得这一任务的难度无论从数学角度还是从编程角度来说大为增加. 下面的例子是对 Shampine 和 Watts 的 Fortran DEABM 代码的局部改编，以放入 NDSolve 框架. 根据 [SG75] 中所介绍的判据，该算法适应性地选择步长和阶数.

第一步是定义系数. 积分方法使用步长可变的系数. 给定一定次序排列的步长值 （其中 是当前所取的步进），使得对于阶数为 的 Adams–Bashforth 预测程序，以及阶数为 的 Adams–Moulton 校正程序 的方法系数，

成立，其中 为均差.

定义一个函数，计算系数 和 ，以及用于误差估计的 . 公式来自[HNW93]，并使用相同的符号表示：

hlist 是过去的步进的步长列表 . 常系数 Adams 系数可以大为简便地一次计算. 由于步长恒定的Adams–Moulton 系数用于改变方法阶数的误差预计，采用保存值的规则一次进行定义是合理的.

定义一个函数，计算并保存固定步长 Adams–Moulton 系数的值：

下一阶段是设置一个数据结构，保存步进之间的必要信息，并对如何初始化数据进行定义. 需要保存的关键信息是过去的步长列表 hlist 以及均差 . 由于该方法进行误差估计，需要从 NDSolve`StateData 对象中得到所用的正确范数是. 为了优化和方便，一些其它数据如精确度等被保存.

定义一个规则，对 NDSolve 中的 AdamsBM 方法进行初始化：

由于在初始化时刻不存在步进，hlist 的初识形式为 {}. 的初识值是初始条件下解向量的导数，因为0阶均差正好是函数值. 注意 是一个矩阵. 初始化为一个零向量的第三个元素，用于确定下一步的最佳阶数. 它事实上是一个额外的均差. 数据其它部分的用法在步进函数的定义中被阐明.

初始化的设置还可以用来得到一个选项的值，该值可用于限制所用方法的最大阶数. 在代码 DEABM 中，它限定在12以内，因为这是机器精度计算的一个实用限度. 但在 Wolfram 系统中可以进行更高精度的计算，使用阶数更高的方法可能更为有利，因此设置这样一个绝对限度是没有必要的. 因此，您设置的选项默认值为 .

这里设置 AdamsBM 方法的 MaxDifferenceOrder 选项的默认值：

在进入更加复杂的 "Step" 方法函数前，有必要定义简单的 "StepMode" 和 "DifferenceOrder" 方法函数.

定义 AdamsBM 方法的步进模式恒为 Automatic. 这意味着它不能从步进控制程序调用，因为这要要求固定步长来改变尺寸：
定义 AdamsBM 方法的差分阶. 这个值随着值保存的过去值的变化而变化：

最后，这里是 "Step" 函数的定义. 采取一个步进的实际过程仅有几行. 剩下的代码应对的是自动步长和阶数选择，与 Shampine 和 Watts 的 DEABM 代码非常相似.

定义 AdamsBM 的 "Step" 方法函数，根据先前介绍的模板返回步进数据：

这里的代码与 DEABM 有一些差异. 最显著的一点是系数在每一步被重新计算，而 DEABM 仅计算需要更新的系数. 做这一改变的目的是为了简化代码，但却是导致了一种性能的损失，对于小型或中型系统来说尤为如此. 另外一个重要的改变是用于限制被拒绝步进的大部分代码留给了 NDSolve，导致这部分代码不对步长是否过小或公差是否过大进行检查. 钢度检测试探法也被忽略了. 确定阶数和步长的代码已经被模块化到不同的函数中.

定义一个函数，构建 、 及 时的误差估计 ，并确定是否需要降低阶数：
定义一个函数，确定一个成功步进后所用的最佳阶数：
定义一个函数，确定步长为 的一个成功步进后所用的最佳步长：

一旦输入了这些定义，您只需使用 Method->AdamsBM 就可以访问 NDSolve 中的方法.

通过先前定义的 Adams 方法求解谐振子方程：
这里显示出计算得到的解的误差. 很明显可以看出误差被限制在一个合理的界限内. 注意经过开始的几个点后，步长增大：

当进行严格公差的高精度计算时，这种方法有可能在性能上超越一些内置方法. 原因是内置方法来自阶数限定在12以内的代码.

系数计算需要大量时间.

这里所用的阶数并不像所预计的那样高.

在任何情况下，大概一般的时间用于系数的生成，因此为了改善计算，您需要找到系数的更新.