NDSolve 方法的插件框架

引言

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

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

AlgorithmIdentifier [data] 一个算法对象,该对象包含某一常微分方程的数值积分算法可能需要使用的任何 data;数据为算法所专有;AlgorithmIdentifier 应是一个 Mathematica 符号,并且算法通过使用选项 TemplateBox[{Method, paclet:ref/Method}, RefLink, BaseStyle -> InlineFormula]->AlgorithmIdentifier 从 TemplateBox[{NDSolve, paclet:ref/NDSolve}, RefLink, BaseStyle -> InlineFormula] 访问;

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

NDSolve 不直接对与算法相关的数据进行访问,因此用户可以以任何方便或高效使用的形式保留所需的信息. 可能保存在私有数据中的算法和信息只能通过该算法对象的方法函数才能访问.The key interaction between a method algorithm object and NDSolve is through a function.

obj["Step"[inputs]]attempt to take a single time step using the numerical algorithm

Required function used by NDSolve.

The following is a list of important properties of the method and its function that should be defined.

obj["StepInputs"]the list of input arguments that the step function should be given
obj["StepOutput"]the list of outputs that the step function will return
obj["StepMode"]set to the step mode
obj["DifferenceOrder"]set to the current asymptotic difference order of the algorithm

Properties that should be defined for using a method algorithm object obj in NDSolve.

The table below gives specifications that can be used for step inputs and/or outputs.

short name
full name    
component
"F""Function"the function given the right-hand side of the ODE system or the residual for a DAE system
"T""Time"current time
"H""TimeStep"time step
"X""DependentVariables"current values of dependent variables other than discrete
"XI""DependentVariablesIncrement" such that gives the values of the dependent variables after the step; if is given as $Failed, then the step is rejected and will be recomputed (the time step returned should be smaller)
"XP""TemporalDerivatives"the time derivatives of the dependent variables
"D""DiscreteVariables"the discrete variables that take on a continuous range of values
"ID""IndexedDiscreteVariables"the discrete variables that take on a finite set of values
"P""Parameters"the parameters other than ones for which sensitivity is being computed
"SP""SensitivityParameters"the parameters for which sensitivity is being computed
"SD""SolutionData"the solution data list
"Obj""MethodObject"a new setting for obj that will be used for the next step
"State""StateData"the object

Step input and output specifications.

The solution data is the list of values with components relating to different variable types. The solution data in the active direction can be obtained from an object state using . Many of the specifications listed above correspond directly to solution data components.

For input, the specification for the function should be given with the arguments it will be given. These arguments should be solution data components. All can be used to indicate that all arguments used by NDSolve should be expected.

For output, if the method only returns some of the values in some cases (e.g. the method object does not change every step) then you can use Null for any output except the increment to indicate that it should not be used.

经典 Runge-Kutta

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

这里定义一个方法函数,使用经典四阶 Runge-Kutta 方法朝向一个常微分方程的积分取一个单一步进. 由于该方法非常简单,不必保存任何私有数据.
In[1]:=
Click for copyable input
This specifies for NDSolve the inputs and outputs for the step function. The template is used to make this work for all instances of CRK4 method objects.
In[2]:=
Click for copyable input
This is an equivalent specification using the full names.
In[4]:=
Click for copyable input
这里定义一个方法函数,使得 NDSolve 可以获得用于该方法的适当差分阶.
In[6]:=
Click for copyable input
这里为步进模式定义一个方法函数,使得 NDSolve 知道如何控制时间步. 这种算法没有任何步进控制,因此步进模式定义为 .
In[7]:=
Click for copyable input
这里用固定步长对简谐振子方程进行积分.
In[8]:=
Click for copyable input
Out[8]=

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

这里利用步进由 方法控制的经典四阶 Runge-Kutta 方法对简谐振子进行积分.
In[9]:=
Click for copyable input
Out[9]=
这里绘出一幅图形来比较步进终点计算得到的解的误差. 使用 方法得到的误差用蓝色表示.
In[10]:=
Click for copyable input
Out[11]=

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

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

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

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

NDSolve 中初始化算法.

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

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

这里定义一个方法函数,使用经典四阶 Runge-Kutta 方法向着常微分方程积分的方向取一单一步进,使用保存的数值作为所需的系数.
In[12]:=
Click for copyable input
这里定义一个规则,使用后来将要用到的数据初始化算法对象.
In[13]:=
Click for copyable input

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

Crank-Nicolson

This is an example of how to set up an implicit method. To solve the system of ODEs , the scheme for a time step of size is , where and . In general, for nonlinear , the equations need to be solved with Newton iteration.

Instead of using a true Newton iteration where the Jacobian is updated at each step, the Jacobian will be computed once and its LU decomposition will be used repeatedly to solve the linear equations. As long as time steps are not too large, this will generally converge as well as the full Newton method and each iteration beyond the first is substantially cheaper.

To avoid possible problems with the Newton iteration not converging, the method will be set up to take method options that control the maximum number of iterations and the convergence tolerance.

Define options and method initialization for the Crank-Nicolson method.
In[14]:=
Click for copyable input
In[15]:=
Click for copyable input
Define the step function.
In[16]:=
Click for copyable input
This specifies for NDSolve the inputs, outputs, difference order, and step mode for the Crank-Nicolson method.
In[17]:=
Click for copyable input
Here is how the method can be called with option settings.
In[21]:=
Click for copyable input
Out[21]=

For a linear problem, the solution is given by the first iteration. However, for a nonlinear problem, convergence may take longer.

This fails because the nonlinear problem cannot be expected to converge in 1 iteration with such large time steps.
In[22]:=
Click for copyable input
Out[22]=
With more iterations the solution is successful.
In[23]:=
Click for copyable input
Out[23]=

Since the Crank-Nicolson method is usually considered regarding solving PDEs, here is an example of the method solving the wave equation. Since this is a linear equation, convergence occurs in 1 iteration so the method is quite fast.

This uses the Crank-Nicolson method to solve the wave equation with periodic boundary conditions.
In[24]:=
Click for copyable input
Out[24]=
A density plot of the solution.
In[25]:=
Click for copyable input
Out[25]=

Invoking Another Method

One of the strengths of the NDSolve framework comes from the capability to nest methods in various ways.

NDSolve`InitializeSubmethod[caller,submspec,stepmode,sd,f,state]initializes the submethod specified by Method->submspec in the initialization of the method caller; stepmode, sd, f, and state are as in InitializeMethod
NDSolve`InvokeMethod[submobj,f,h,sd,state]invoke the submethod with method object submobj with function f, time step h, solution data sd, and object state, returning a list where hnew is the new time step, sdnew is the solution data at the step end with unchanged components set to Null, and submobjnew is the new submethod object

Using submethods.

InitializeSubmethod effectively calls after separating out the method name from any suboptions that are given. If successful, it returns the method object that you can later use in InvokeMethod.

InvokeMethod effectively constructs a function for the submethod. Since the arguments for the function are specified by the specification for the submethod, the function f given to InvokeMethod should have all the arguments used by NDSolve and can be specified in the caller specification by "Function"[All]. Note that for built-in methods, there may not be an explicit function and InvokeMethod handles these using efficient internal code. The result returned by InvokeMethod is always constructed from the submethod's actual output to match a of .

Here is an example that simply calls a submethod to provide simple monitoring of the steps and difference order.

Define the method options and the initialization for the method.
In[1]:=
Click for copyable input
In[2]:=
Click for copyable input
The step function for the method simply invokes the submethod and applies the monitor function to the result and the difference order.
In[3]:=
Click for copyable input
Set the inputs and outputs.
In[4]:=
Click for copyable input
Derive all other properties from the submethod.
In[6]:=
Click for copyable input
Use the method.
In[7]:=
Click for copyable input
Out[7]=

For a more complicated integration the method can be used with Reap to collect data.

Collect the logarithm of the step size and difference order as a function of time for integrating the Van der Pol oscillator and plot the solution.
In[8]:=
Click for copyable input
Out[9]=
Show the logarithm of the step size and the difference order for steps near the sudden change.
In[10]:=
Click for copyable input
Out[10]=

Adams 方法

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

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

成立,其中 为均差.

定义一个函数,计算系数 ,以及用于误差估计的 . 公式来自[HNW93],并使用相同的符号表示.
In[29]:=
Click for copyable input

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

定义一个函数,计算并保存固定步长 Adams-Moulton 系数的值.
In[30]:=
Click for copyable input

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

定义一个规则,对 NDSolve 中的 方法进行初始化.
In[32]:=
Click for copyable input

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

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

这里设置 方法的 选项的默认值.
In[33]:=
Click for copyable input

在进入更加复杂的 方法函数前,有必要定义简单的 方法函数.

定义 方法的步进模式恒为 Automatic. 这意味着它不能从步进控制程序调用,因为这要要求固定步长来改变尺寸.
In[34]:=
Click for copyable input
定义 方法的差分阶. 这个值随着值保存的过去值的变化而变化.
In[35]:=
Click for copyable input

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

定义 方法函数,根据先前介绍的模板返回步进数据.
In[36]:=
Click for copyable input

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

定义一个函数,构建 时的误差估计 ,并确定是否需要降低阶数.
In[37]:=
Click for copyable input
定义一个函数,确定一个成功步进后所用的最佳阶数.
In[38]:=
Click for copyable input
定义一个函数,确定步长为 的一个成功步进后所用的最佳步长.
In[40]:=
Click for copyable input

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

通过先前定义的 Adams 方法求解谐振子方程.
In[41]:=
Click for copyable input
Out[41]=
这里显示出计算得到的解的误差. 很明显可以看出误差被限制在一个合理的界限内. 注意经过开始的几个点后,步长增大.
In[42]:=
Click for copyable input
Out[42]=

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

In[43]:=
Click for copyable input

系数计算需要大量时间.

In[45]:=
Click for copyable input
Out[45]=

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

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

In[46]:=
Click for copyable input
Out[46]=
New to Mathematica? Find your learning path »
Have a question? Ask support »