组件和数据结构

引言

NDSolve 分为几个基本步骤. 对于高级的使用,有时候通过访问这些组件来分别执行这些步骤是有利的.

  • 方程处理和方法选择
  • 方法初始化
  • 数值解
    • 解的处理

    NDSolve 是在内部执行这些步骤的,对于临时用户隐藏了详细信息.

    以下是用来划分这些步骤的底层函数.

    把微分系统分为初值问题、边值问题、微分代数问题、偏微分方程问题,等等. 它也选择适当的默认积分方法,并且构建主要的 数据结构.

    推进了数值解. 它的第一次调用(可以有多次)初始化数值积分方法.

    把数值数据转换成一个 InterpolatingFunction 来表示每个解.

    请注意 可能会用掉全部时间的大部分来求解一个微分系统. 在这种情况下,只执行一次该步骤并且使用 对不同的选项或者初始条件来反复求解是有益的.

示例

处理方程并且建立求解微分系统的数据结构.
In[887]:=
Click for copyable input
Out[887]=
初始化方法 并且把系统积分到时间10. 的返回值是 Null 以避免多余的引用导致不必要的复制.
In[888]:=
Click for copyable input
Out[888]=
把每组解数据转换为一个 InterpolatingFunction.
In[889]:=
Click for copyable input
Out[889]=
把解表示为一个 InterpolatingFunction 允许连续的输出,甚至对于不是数值解网格的一部分的点也是这样.
In[890]:=
Click for copyable input
Out[890]=

ProcessEquations

任何使用 NDSolve 的解的第一个阶段是对被指定为一种可以被实际积分算法有效访问的形式的方程进行处理. 该阶段最小包含决定每个变量的微分阶数,作必要的替换以获取一个一阶系统,以函数的形式求解关于函数的时间导数,并且把结果形成为一个 对象. 如果用户想要节省对同一组方程重复这一过程的时间,或者用户想要对数值积分过程更多的控制,处理阶段可以用 分开执行.

NDSolve`ProcessEquations[{eqn1,eqn2,},{u1,u2,},t]
将函数 的微分方程 处理为一个标准形式;返回 对象的一个列表,其中包含解以及与用函数形式表达的函数的时间导数的每个解相关的数据;t 可以像在 NDSolve 中那样在一个数值范围内的一个列表中指定
NDSolve`ProcessEquations[{eqn1,eqn2,},{u1,u2,},{x1,x1min,x1max},{x2,x2min,x2max},]
将函数 的偏微分方程 处理为一个规范形式;返回 对象的一个列表,其中包含解以及与用函数形式表达的函数的时间导数的每个解相关的数据;如果 是时间变量,它不需要使用边界 来指定

NDSolve 的方程处理.

这里创建一个由两个 对象组成的列表,因为对于 有两个以 表达的可能解.
In[891]:=
Click for copyable input
Out[891]=

重新初始化

一个更加复杂问题的解包含反复求解相同的微分方程,这种情况并不罕见,但是它们具有不同的初始条件. 在某些情况下,处理方程可能和对微分方程数值积分一样耗时. 在这种情况下,能够简单地提供新的初始值更加具有优势.

NDSolve`Reinitialize[state,conditions]假设方程和变量与用来创建 对象 state 的那些相同的情况下,形成新的 对象的一个列表,每个对象对应函数的初值方程 conditions 的每个可能解

对处理过的方程的再利用.

这里创建谐振子的一个 对象.
In[892]:=
Click for copyable input
Out[892]=
这里创建三个新的 对象,每个都具有一个不同的初始条件.
In[893]:=
Click for copyable input
Out[893]=

当用户需要对许多不同的初始条件求解相同的微分方程时,使用 可以节省计算时间,就如同在边值问题的打靶法中一样.

NDSolve 选项的一个子集可以被指定为 的选项.

这里创建一个新的 对象,指定了一个初始步长.
In[894]:=
Click for copyable input
Out[894]=

迭代解

对象的一个重要用途是对积分有更多的控制. 对于某些问题,一个适当的做法是先检查解,然后根据一定的条件重新开始或者改变参数.

NDSolve`Iterate[state,t]从当前时间直到时间 t,计算已经被赋值给变量 state 的一个 对象中微分方程的解

对微分方程的解进行迭代.

这里创建一个 对象,它包含对一个具有变化系数的振荡器使用显式RungeKutta方法求解方程所需的信息.
In[895]:=
Click for copyable input
Out[895]=

请注意,当用户使用 时,不必明确提供 变量的范围,因为对于以一种可求解的方式建立方程来说,该信息不是必要的. (然而对于偏微分方程,用户必须给出所有空间变量的范围,因为这些信息对于确定一个适当的离散化是必要的.)

这里计算到达时间 的解.
In[896]:=
Click for copyable input
Out[896]=

不返回值,因为它修改赋给变量state的 对象. 因此,该命令以一种与设置列表的部分相似的方式影响变量值,如 "通过索引操作列表" 中描述的. 用户可以看到 state 的值已经改变,因为它现在显示的是积分所到达的当前时间.

state 的输出形式显示解的积分经过的时间范围.
In[897]:=
Click for copyable input
Out[897]=

如果用户想要进一步积分,用户可以再次调用 ,但是使用时间的一个较大值.

这里计算解到时间 .
In[898]:=
Click for copyable input
Out[898]=

用户可以指定一个早于第一当前时间的时间,在这种情况下,积分对时间反向进行.

这里从初始条件向后计算解,至 .
In[899]:=
Click for copyable input
Out[899]=

允许用户指定停止积分的中间时间. 这会很有用,例如,可以避免不连续性. 通常情况下,该策略对于所谓的 one-step 方法更加有效,例如在该例子中使用的显式 RungeKutta方法.但是,它通常也适用于默认的NDSolve 方法.

这里计算解至 ,确保解对于在 , , 上的系数的不连续点没有问题.
In[900]:=
Click for copyable input
Out[900]=

获取解函数

一旦把系统积分到一定时间,通常用户希望能够看到当前的解值,并且产生一个表示至今计算的解的近似函数. 命令 允许用户做到这两点.

NDSolve`ProcessSolutions[state]给出已经在 state 中计算的解,作为具有 InterpolatingFunction 对象的一个规则列表

获取作为 InterpolatingFunction 对象的解.

这里作为一个 InterpolatingFunction 对象提取 前一节 中计算所得的解.
In[901]:=
Click for copyable input
Out[901]=
这里画出解的图形.
In[902]:=
Click for copyable input
Out[902]=

正如直接使用 NDSolve 一样,在 的第二个变量中指定的每个函数都将会有一个规则. 只有解的指定分量被以一个 InterpolatingFunction 对象可以创建的方式保存.

NDSolve`ProcessSolutions[state,dir]以具有函数及其导数值的规则的一个列表的形式,给出在 中在方向 上最近计算的解

获取当前解的值.

这里给出在向前方向上的当前解的值和导数.
In[903]:=
Click for copyable input
Out[903]=

对于方向 dir 用户可以给出的选择是 ,它们指的是从初始条件向前和向后积分.

"Forward"在时间变量的值增加的方向上积分
"Backward"在时间变量的值减少的方向上积分
"Active"在当前积分的方向上积分;通常,该值应该只能从在一个active的积分中使用的方法初始化调用

积分方向的规格.

给出的输出总是以应变量形式给出,要么对应于自变量的某个特定值,要么对应于所有保存值上的插值. 这意味着,当一个偏微分方程被积分时,用户将会获得代表在空间变量上的应变量的结果.

这里计算从时间 的热方程的解.
In[904]:=
Click for copyable input
Out[905]=
这里给出 处的解.
In[906]:=
Click for copyable input
Out[906]=

解作为一个在空间变量 上插值的 InterpolatingFunction 对象给出.

这里给出 处的解.
In[907]:=
Click for copyable input
Out[907]=

当用户处理偏微分方程的当前解时,空间误差估计将被核查. (它一般不被核查,除非当解产生时,因为这样做会非常耗时.)因为量过大,所以发出 NDSolve::eerr 消息. 将词向后特别与热方程相联系以显示不稳定性,这告诉我们在这个例子中错误之处的线索.

下面是 上的解的图形.
In[908]:=
Click for copyable input
Out[908]=

解的图形表明不稳定性的确是个问题.

虽然热方程的例子足够简单到可以知道在时间上反向的解是有问题的,使用 来监测一个偏微分方程的解可以用来挽救一个不如预期准确的解的计算. 另一种简单的监测形式如下.    

输入以下命令生成了一系列图示,显示了正在计算的 sine-Gordon 方程的推广的解.
In[909]:=
Click for copyable input
Out[911]=

当用户以这种方式对解进行监控时,通常可以中断计算,如果用户看到找到的解已足够的话. 用户仍然可以使用 对象以获取已计算的解.

NDSolve`StateData 属性

一个 对象包含了很多信息,但是它的组成方式是使得解的迭代更加容易,而不是使得对有关信息保存的位置的理解更加容易. 但是,有时用户会希望从 state 数据对象获取信息:为此,我们定义了一些属性来使得对信息的获取更加容易.

来自 对象的信息的一个最重要的用途是设置积分方法. 范例在 "NDSolve 方法的插件框架" 中给出.

state@"Variables"给出对应于 解数据 的所有变量的列表组成的列表
state@"VariableDimensions"给出每个变量的维度组成的列表
state@"VariablePositions"给出每个变量的 解数据 for 中的位置组成的列表
state@"NumericalFunction"给出计算对应于时间变量 t 的解向量的导数所用的 NumericalFunction 对象
state@"ProcessExpression"[args,expr,dims]
使用与 生成 state 相同的变量转换来处理表达式 expr,以给出 NumericalFunction 对象 用于数值计算 expr; args 是数值函数的变量,它要么是 All 要么是系统的应变量的一个变量列表; dims 应该是 Automatic 或者是一个给出数值函数结果的期望维度的显式列表
state@"SystemSize"给出正在求解的一阶常微分方程的有效数目
state@"MaxSteps"给出迭代微分方程所允许的最大步骤数目
state@"WorkingPrecision"给出用来求解方程所用的工作精度
state@"Norm"用于仪表误差的规范范式scaled norm)

对象 state 的属性.

下面是钟摆方程的 对象.
In[912]:=
Click for copyable input
Out[912]=
获取系统规模.
In[913]:=
Click for copyable input
Out[913]=

因为存在两个标量常微分方程,系统规模是2.

现有的信息很大程度上取决于当前的解值. 从初始条件时间开始,NDSolve 既对向前方向的解,也对向后方向的解保留解的信息. 这些都可以通过下列方向属性访问.

state@"SolutionData"[dir]给出积分方向dir上的解数据的当前值
state@"TimeStep"[dir]给出在积分方向dir上下一步的时间步长
state@"TimeStepsUsed"[dir]给出到达积分方向dir上的当前时间所需的时间步骤的数目
state@"MethodData"[dir]给出在积分方向 dir 上所用的方法数据对象

一个 对象 state 的定向方法函数.

如果方向变量被省略,该函数将返回一个具有两个方向上的数据的列表(该列表在初始条件下具有单个元素). 否则,方向可以是 或者 ,如 上一个子章节 中指明的.

时间步骤直至积分开始才初始化,除非已经使用 StartingStepSize 选项明确指定.
In[914]:=
Click for copyable input
Out[914]=
积分至 并且在向前方向获取时间步骤.
In[915]:=
Click for copyable input
Out[916]=
用于向后和向前方向的时间步骤的数目.
In[917]:=
Click for copyable input
Out[917]=
获取钟摆范例的解数据.
In[918]:=
Click for copyable input
Out[918]=

解数据

解数据包括每个类型的解分量的数值列表组成的列表. 解数据分量是不同类型的变量,例如,时间、应变量和离散变量. 不用于特定问题的任意分量以空列表 () 或者 None 的形式给出.

NDSolve`SolutionDataComponent[sd,c]从解数据 sd 获取解分量 c
NDSolve`SetSolutionDataComponent[sd,c,v]在赋给符号 sd 的解数据中,把解分量 c 设为 v
nf@"SolutionDataComponents"NDSolve 中使用的NumericalFunction nf 的参数给出解数据分量

获取并且设置解数据分量.

解数据的分量在下列表格中给出.

简短名称
完全名称    
分量
"T""Time"当前时间
"S""Space"空间离散化
"X""DependentVariables"非离散应变量的当前值
"X'""TemporalDerivatives"应变量的时间导数
"D""DiscreteVariables"使用连续范围的数值的离散变量
"ID""IndexedDiscreteVariables"使用有限集合上的离散变量
"P""Parameters"除了计算敏感度以外的参数
"SP""SensitivityParameters"计算敏感度的参数

解数据分量.

下面是 对象,用于具有计算到 的离散变量的钟摆方程的解.
In[919]:=
Click for copyable input
Out[921]=
这显示问题中使用的变量.
In[922]:=
Click for copyable input
Out[922]=
显示应变量.
In[923]:=
Click for copyable input
Out[923]=
下面从向后和向前两个方向获取当前解数据.
In[924]:=
Click for copyable input
Out[924]=
下面获取向前方向应变量分量的数值.
In[925]:=
Click for copyable input
Out[925]=
下面获取向后方向的离散变量.
In[926]:=
Click for copyable input
Out[926]=
下面获取向后方向的索引离散变量. 结果是一个空列表,因为在该计算中没有使用索引离散变量.
In[927]:=
Click for copyable input
Out[927]=

注意,解数据是 NDSolve 用于计算解的原始数据. 在上面的例子中,应变量分量是长度为2的列表,因为二阶方程自动化简为一阶系统. 如果您想要当前解的部分的清晰识别,您可以使用 .

通过规则获取与变量相关联的解数据.
In[928]:=
Click for copyable input
Out[928]=

NumericalFunction

当进行时间积分时,方法 NDSolve 使用常微分方程右边函数的操作,以满足 ,或者 NumericalFunction 的残差函数,以满足 ,其中 是实数,而 和其他解数据分类所使用的是向量. 由向量 表示的实际变量可能具有更复杂的结构,因此 NDSolve 使用具有向量参数的中间NumericalFunction 对象,并且返回基于任意内在函数的向量.

获取具有向量和标量变量的系统的 对象.
In[929]:=
Click for copyable input
Out[929]=
获取计算导数的 NumericalFunction.
In[930]:=
Click for copyable input
Out[930]=

NDSolve 对计算构建的 NumericalFunction 只作为 solution data 分量实际在问题中的参数使用. 您可以找到通过使用 属性它实际使用哪些分量.

求计算需要哪些解数据分量.
In[931]:=
Click for copyable input
Out[931]=
获取初始解数据和计算参数.
In[932]:=
Click for copyable input
Out[932]=
In[933]:=
Click for copyable input
Out[933]=

注意,索引离散变量通过 Sign 函数的自动不连续处理以 变量引入时. (查看下面给出的解).

由于系统是一个常微分方程,计算函数给出时间导数.
In[934]:=
Click for copyable input
Out[934]=
In[935]:=
Click for copyable input
Out[935]=

若要使计算给定问题的解数据变得更轻松,两个方便函数已经被包含,以使用解数据计算.

NDSolve`EvaluateWithSolutionData[nf,sd]使用解数据列表 sd 的参数计算 NumericalFunction nf
NDSolve`EvaluateJacobianWithSolutionData[nf,sd,c]关于解数据分量 c 中的所有变量,计算雅克比导数

使用解数据计算.

使用解数据列表计算函数.
In[936]:=
Click for copyable input
Out[936]=

NumericalFunction 所做的一件事是自动处理导数. 默认情况下,如果可以找到,则使用符号导数,否则使用有限差值.

关于应变量计算雅克比.
In[937]:=
Click for copyable input
Out[937]=

在这个例子中,使用有效差值,因为向量参数和 Dot 乘积.

积分并且绘制 的解.
In[938]:=
Click for copyable input
Out[939]=

有时候设置 NumericalFunction 来计算 NDSolve 变量的某些函数是有用的,这可以使用 完成.

在方法初始化过程中,许多内置方法完成这些,因此函数可以有效地重复计算. 例如,投射方法设置函数以计算不变量从初始值的偏离,并且在每个步骤后面使用导数完成投射.

设置振荡器的 对象.
生成能量积分的 NumericalFunction.
In[941]:=
Click for copyable input
Out[941]=
在初始条件处计算.
In[942]:=
Click for copyable input
Out[942]=
生成在积分过程中在规则区间上监控的能量图线.
In[943]:=
Click for copyable input
Out[943]=
考虑到解在 处快速振荡,能量偏离实际上很小.
In[944]:=
Click for copyable input
Out[944]=