数值直线方法

引言边界条件
空间导数逼近空间误差估计

引言

数值直线方法是一种求解偏微分方程的技术，它将除某一维度以外的其他维离散化，把得到的半离散问题作为常微分方程组或微分代数方程组积分求解. 这个方法的显著优点是它可以用上已经针对数值求解常微分方程和微分代数方程组开发的成熟的，一般用途的方法和软件. 对适用直线方法的偏微分方程，通常可以证明该方法是非常高效的.

该方法要求偏微分方程问题必须至少在一个维度上是个适定的初值（柯西）问题，因为所用的常微分方程和微分代数方程的积分程序是针对初值问题的求解程序. 这就排除了纯椭圆方程，如拉普拉斯方程等，但仍然留下了一大类演化方程，可以相当有效的求解.

简单的例子比文字更能展示该方法的基本思想. 考虑下面的问题（一个土壤热量随季节变化的简单模型）：

这可以用直线方法求解因为有初值 .

问题（1）将利用二阶有限差分把 离散化，特别是使用了下面这个近似公式：

尽管有限差分离散化是最常见的，但直线方法要用的离散化并不一定要用有限差分完成；有限体积甚至有限元离散化也都可以使用.

若要使用前述的离散化方法，可选择一个均匀网格 且间隔为 从而满足 . 记 为 的值. 出于演示问题配置的目的，这里选择了一个特别的 .

这里为后续命令定义了一个特别的 值以及对应的 值. 可以改变这些值来进行更精细或更粗略的空间逼近：
定义向量 ：

对 ，可以利用中心差分公式（2）来得到常微分方程组. 但是，在这样做之前，纳入边界条件是非常有用的.

在 处的狄利克雷边界条件可以把 简单定义为 的函数来轻松处理. 另一种方法是对边界条件做关于时间的微分后加入边界条件. 这个例子将使用后一种方法.

在 处的诺伊曼边界条件要更困难些. 关于二阶差分，一种处理办法是通过反射：设想在区间 上用在 和 处相同的边界条件求解问题. 因为初始条件和边界条件关于 对称，所以解应该始终关于 对称，于是对称等价于 处的诺伊曼边界条件. 因此 即 .

通过 ListCorrelate 应用差分公式. 填充用的 {un–1[t]} 实现了诺伊曼边界条件：
设置零初始条件：

现在偏微分方程已部分离散化变成常微分方程初值问题，可以用 NDSolve 中的常微分方程积分程序求解.

求解常微分方程初值问题：
把解 作为 和 的函数绘制展示：

该图形表明了为什么这一技术被称为数值“直线方法”.

直线之间的解可用插值法求得. 当 NDSolve 计算偏微分方程的解时，结果是二维的 InterpolatingFunction.

这里使用 NDSolve 直接计算热传导方程（1）的解：
绘制解的曲面图：

设定 时并没有给出非常精确的解. 当 NDSolve 计算解时，它使用初始条件下的空间误差估计来确定网格间距该是多少. 时间（或至少是类时间）变量的误差则由自适应的常微分方程积分程序处理.

在例（1）中，时间和空间的区别在问题上下文中很明显. 即使在没有明显区别的时候，本教程也将称之为“空间”和“时间”变量. “空间”变量是那些要离散化的变量. “时间”变量是在常微分方程组中剩下的那个有待积分的变量.

选项名称
默认值
TemporalVariableAutomatic在导出的常微分方程组或微分代数方程组中仍然保留导数项的那些变量
MethodAutomatic使用何种方法对常微分方程或微分代数方程积分
SpatialDiscretizationAutomatic使用何种方法进行空间离散化
DifferentiateBoundaryConditionsTrue是否对边界条件做关于时间变量的微分
ExpandEquationsSymbolicallyFalse是否对有效函数进行符号式扩展
DiscretizedMonitorVariablesFalse是否把由 WhenEvent 给出的，由诸如 StepMonitor 这样的监视器给出的，或者由 Projection 等方法的方法选项给出的相关变量解释为空间变量的函数或者表示空间离散值向量

NDSolve`MethodOfLines 的选项.

部分选项的使用需要进一步了解直线方法的工作方式，这将在后续的章节中解释.

目前实现的空间离散化的方法是 "TensorProductGrid" 和 "FiniteElement". 两种方法都有自己的一组选项可以用来控制离散化过程的细节. "FiniteElement" 方法是非常普适的，可用于任意形状的区域；它有自己的教程. 本教程的后续章节将主要使用 "TensorProductGrid" 离散化方法.

空间导数逼近

有限差分

有限差分概念的实质是体现在导数的标准定义中的

不同之处在于导数是 趋向于零时的极限，有限差分使用的是与相邻点 之间的有限间距，于是得到的是导数近似值

差分公式也可以由泰勒公式导出

这种形式更有用，因为它提供了误差估计（假设足够平滑）

此公式的一个要点是 必须位于 和 之间，使得误差局限在包含采样点的区间内部. 对于有限差分公式，相对于模板或样本点集的误差是局部的这点在一般情况下都成立. 对于收敛性或其它分析，通常把误差表示为渐近形式：

此公式通常被称为一阶前向差分. 后向差分使用的是 .

可以用泰勒公式方便的推导高阶近似公式. 例如，从

中减去

并求解 ，可得到一阶导数的二阶中心差分公式，

如果把刚才的泰勒公式再多扩展一阶并相加，结合刚才给出的公式，不难得到二阶导数的中心公式

注意到点之间一致为 的步长使得写出上述公式比较方便，但这并不是必须的. 比如一般情形下的二阶导数近似公式是

其中 对应的是最大的局部网格间距. 请注意请注意三点公式的渐近阶数已降至一阶；在均匀网格上是二阶则缘自对应各项恰好相消.

一般而言，只要样本点的数量足够，渐进误差为任意阶的任意阶导数公式都可由泰勒公式推导得到. 然而若把该方法用于比给出的简单例子更复杂的问题，则会变得繁琐而低效. 另一种方法则是基于多项式插值的：泰勒公式对于足够低阶的多项式来说是精确的（无误差项），有限差分公式也同样如此. 不难证明有限差分公式等价于插值多项式的导数. 例如，一种简单的推导刚才的二阶导数公式的方法就是做一个二次多项式的插值并找到它的二阶导数（其实就是首项系数）.

这里对三个点 、 和 做多项式插值，并通过求导得到了二阶导数的三点有限差分公式：

通过这种形式的公式，容易看出它实际上是前向和后向一阶导数逼近之差. 有时用这种形式的有限差分是有其优点的，特别是对于导数内部含有系数的项，如通常出现在偏微分方程中的 .

在考虑插值公式时另一个性质也很明显：得到导数近似值的点不一定要在网格上. 用到此公式的一个常见情况是交错网格，其中想要的导数可能在网格点的中点处.

在均匀交错网格 上生成其一阶导数的四阶近似。主网格点 在 处， 为奇数：

此式的四阶误差系数为 而下面导出的标准四阶公式的则是 . 误差的减少大部分是因为模板尺寸的减小.

从均匀网格上的一点生成一阶导数的四阶近似：

一般而言，使用 个点的有限差分公式在函数为 次多项式，且渐进阶至少为 时是精确的. 在均匀网格上，可以期望得到更高的渐进阶数，特别是对于中心差分而言.

使用高效的多项式插值技术是生成系数的合理方法，但 B. Fornberg 开发了一种快速算法用于有限差分权重生成 [F92]、[F98]，这种方法要快得多.

在 [F98]，Fornberg 给出了显式有限差分的一行 Wolfram 语言公式.

在均匀网格上生成权重的 Fornberg 的简单公式. 这里，此式被稍作修改，改为函数定义：

这里 是导数的阶， 是模板内的网格区间数，而 是导数近似所在的点与模板最左边的点之间的网格区间数. 不必是整数；非整数值会直接生成交错网格近似公式. 令 为 总会生成中心公式.

使用 Fornberg 公式生成一阶导数的交错四阶逼近. 与之前用 InterpolatingPolynomial 计算得到的结果相同：

下面的表格列出了一些常用的有限差分公式以作参考.

公式
误差项

一阶导数在均匀网格上的有限差分公式.

公式
误差项

二阶导数在均匀网格上的有限差分公式.

这个表格中要注意的一点是，公式与中心偏离越远，误差项系数越大，有时大几百倍. 因此有时在需要单边的导数公式时（如在边界处），会用更高阶的公式来抵消额外的误差.

NDSolve`FiniteDifferenceDerivative

Fornberg [F92]、[F98] 还给出了一个算法，虽然不是非常优雅简单，但更具有一般性且特别适用于非均匀网格. 在 Wolfram 语言中对它编程并不难，但为了使之尽可能高效，我们还是提供了一个新的内核函数作为更简单的接口（还有一些额外的特性）.

NDSolve`FiniteDifferenceDerivative[Derivative[m],grid,values]
求网格上取值的函数的 m 阶导数的近似
NDSolve`FiniteDifferenceDerivative[Derivative[m1,m2,…,mn],{grid1,grid2,…,gridn},values]
求 n 变量函数的 (m1, m2, …, mn) 阶偏导数的近似，函数在由 (grid1, grid2, …, gridn) 外积定义的张量积网格上取值
NDSolve`FiniteDifferenceDerivative[Derivative[m1,m2,…,mn],{grid1,grid2,…,gridn}]
计算所需的有限差分权重，以求 n 变量函数的 (m1, m2, …, mn) 阶偏导数的近似，函数在由 (grid1, grid2, …, gridn) 外积定义的张量积网格上取值；结果以 NDSolve`FiniteDifferenceDerivativeFunction 形式返回，可重复作用于网格上的值

找寻导数的有限差分近似.

定义了一个均匀网格，各点之间距离为符号 h：
给出了符号函数 f 在网格上的一阶导数公式：

端点处的导数是用单边公式计算的. 前一个例子中的公式精确到四阶，这是默认设置. 一般而言若使用符号网格和/或数据，得到的就是符号公式. 这对于分析相关方法非常有用；然而对于实际的数值网格，NDSolve`FiniteDifferenceDerivative 用数值网格比用符号公式更快也更精确.

定义了一个 到 之间随机间距的网格：
在网格的每个点上求正弦函数导数的近似值：
显示近似值的误差：

在多维情况下，NDSolve`FiniteDifferenceDerivative 适用于张量积网格，您只需指定每个维度上的网格点.

定义在 和 方向上的网格 xgrid 和 ygrid，得到了高斯函数在 xgrid 和 ygrid 张量积上的混合 偏导数的逼近，并绘制误差的曲面图：

请注意需要在矩阵中给出的值应该对应于网格坐标的外积.

NDSolve`FiniteDifferenceDerivative 不计算导数和的权重. 这意味着对于一般的算符如拉普拉斯算符，您需要合并两个近似值.

定义一个函数，在张量积网格上逼近拉普拉斯算符：
使用该函数对在前一例子中所用的相同网格和高斯函数求拉普拉斯算符的近似值：
选项名称
默认值
"DifferenceOrder"4误差的渐进阶数
PeriodicInterpolationFalse是否考虑周期等于网格所包含的区间的周期函数的值

NDSolve`FiniteDifferenceDerivative 的选项.

假设函数周期性重复，在之前定义的随机网格上求正弦函数的导数的近似值：

当使用 PeriodicInterpolation->True 时，可以忽略值中的最后一个点，因为它总是与第一个点相同. 这个特点在求解带周期性边界条件的偏微分方程时很有用.

生成符号函数一阶导数的二阶有限差分式：

四阶差分通常能够均衡机器精度的截断（逼近）误差和舍入误差. 然而在有些应用场景下四阶差分会产生过度振荡（Gibb 现象），因此二阶差分更好. 另外高精度情况下高阶差分可能更合适. "DifferenceOrder" 的偶数值使用中心公式，这时的误差系数比非中心公式时要小，因此若情况合适还是推荐使用偶数值.

NDSolve`FiniteDifferenceDerivativeFunction

在计算偏微分方程的解时，通常需要将同一个有限差分近似值反复应用于相同网格上的不同值. 通过存储必要的权重计算并将其作用于变化的数据，可以显著节省计算量. 如果在 NDSolve`FiniteDifferenceDerivative 中忽略带有函数值的（第三个）参数，结果将是 NDSolve`FiniteDifferenceDerivativeFunction，一个把权重计算高效存储以便将来重复利用的数据对象.

NDSolve`FiniteDifferenceDerivative[{m1,m2,…},{grid1,grid2,…}]
计算 n 变量函数的 (m1, m2, …) 阶偏导数近似所需要的有限差分权重，函数定义在由 (grid1, grid2, …) 外积定义的张量积网格上返回结果是 NDSolve`FiniteDifferenceDerivativeFunction 对象
NDSolve`FiniteDifferenceDerivativeFunction[Derivative[m],data]
含有快速逼近函数 m 阶导数所需的权重及其它数据的数据对象；在标准输出形式中，仅显示它所逼近的 Derivative[m] 运算符
NDSolve`FiniteDifferenceDerivativeFunction[data][values]
计算函数的导数近似值，函数在用于确定数据的网格上取值

计算可重复使用的有限差分权重.

定义了一个在单位区间上带有 25 个点的均匀网格，并计算网格上一个周期的正弦函数值：
定义了一个 NDSolve`FiniteDifferenceDerivativeFunction，它可重复作用于网格上的不同值以逼近二阶导数：

注意标准输出形式被缩写了，仅显示被逼近的导数运算符.

计算正弦函数的二阶导数近似值：

此函数仅适用于在构建函数的特定网格上定义的值. 如果问题要求改变网格，则需要使用 NDSolve`FiniteDifferenceDerivative 来生成对应于每次网格变化的权重. 然而，若能使用 NDSolve`FiniteDifferenceDerivativeFunction 对象，则计算速度将显著加快.

比较用刚才定义的函数和用之前一节的函数计算拉普拉斯运算符所需的时间. 因为速度过快，不能通过 Timing 来显示时间的不同，故而在每个情形下都用了循环来重复计算：

NDSolve`FiniteDifferenceDerivativeFunction 可重复用于多种情形. 作为一个简单的例子，考虑在单位区间上求解边值问题

的搭配法. （这个简单方法未必是求解该特定问题的最好方法，但作为示例非常有用.）

定义一个函数，该函数的所有分量在边界值问题的近似解处将为零. 使用中间向量 并将其端点设置为 0（{1,-1} 部分）是一个快速而简单的强制满足边界条件的技巧. 除了数值 以外的所有求值都被阻止，因为这在其它情况下不可行.（另外，由于 Times 属性为 Listable，Sin[2 Pi grid] u 将逐分量被乘.）
用 FindRoot 找到近似特征函数并绘制其图像，该特征函数的初始值用的是常系数情形：

由于此问题的设置非常简单，很容易比较其它不同方案. 例如，要将上面的使用默认四阶差分的解与通常使用二阶差分的作比较，只需要改变 "DifferenceOrder".

用二阶差分求解该边界值问题，并将它与四阶解的差异绘制成图：

一种确定哪个解更好的方式是研究在网格不断精化下的收敛性. 这在“微分矩阵”一节中有一定程度的讨论.

尽管 NDSolve`FiniteDifferenceDerivativeFunction 对象的最关键信息——导数的阶——已经显示在输出形式中，但把该信息及其它信息从 NDSolve`FiniteDifferenceDerivativeFunction 中提取出来有时也是很有用的，比如用在某程序中. 数据存储方式的结构可能随Wolfram 语言版本的不同而变化，因此不推荐用部分表达式来提取信息. 更好的替代方式是使用任一种专为此目的提供的方法函数.

设 FDDF 表示一个 NDSolve`FiniteDifferenceDerivativeFunction[data] 对象.

FDDF@"DerivativeOrder"得到 FDDF 所逼近的导数阶
FDDF@"DifferenceOrder"得到每一维上用于逼近的带差分阶数的列表
FDDF@"PeriodicInterpolation"得到带 True 或 False 的列表，指明周期性插值是否用于每一维
FDDF@"Coordinates"得到带有每一维中网格坐标的列表
FDDF@"Grid"形成网格点的张量；这是网格坐标的外积
FDDF@"DifferentiationMatrix"计算稀疏微分矩阵 mat 使得 mat.Flatten[values] 等价于 Flatten[FDDF[values]]

从 NDSolve`FiniteDifferenceDerivativeFunction[data] 对象提取信息的方法函数.

任一方法函数，只要它返回的是带有每一维元素的列表，都可以结合一个整数参数 dim, 它将仅返回那个特定维度上的值，满足 FDDF@method[dim]=(FDDF@method)[[dim]].

下面的例子展示了如何使用这些方法.

这是一个由每一维上有 10 到 16 个点的随即网格创建的 NDSolve`FiniteDifferenceDerivativeFunction 对象：
显示外积网格的维数：

注意网格点张量的秩比张量积的维数多一. 这是由于定义每个点的三个坐标自身在一个列表中. 如果有一个依赖网格变量的函数，则可以使用 Apply[f,fddf["Grid"],{n}]，其中 n=Length[fddf["DerivativeOrder"]] 是所要逼近导数的空间维数.

定义了一个三变量高斯函数并作用于定义了 NDSolve`FiniteDifferenceDerivativeFunction 的网格上：
这里显示了网格点的三维图像，根据导数的数量值着色：

对于像这个例子一样中等大小的张量积网格，用 Apply 还算快. 然而当网格尺寸更大时，这种方法可能不是最快的，因为 Apply 和 Wolfram 语言编译器只有寥寥几种结合使用的方式，因此于压缩数组配合使用的方式也很有限. 若能自己定义函数则可用 Map 代替 Apply，这样有可能可以用 CompiledFunction 加速，因为在 Wolfram 语言编译器中 Map 比 Apply 有更广的应用性.

用 Map 定义了一个得到网格上值的 CompiledFunction.（如果第一个网格维度大于系统选项 "MapCompileLength"，则不必构造 CompiledFunction 因为网格是压缩数组时编译将自动进行.）

如果可能的话，还有一种更好的方法，就是当函数是由具备 Listable 属性的运算和函数组成时，可利用可列表性. 技巧是把张量积网格上各点处的 、 和 值分离. 做到这点最快的方式是使用 Transpose[fddf["Grid"]],RotateLeft[Range[n+1]]]，其中 n=Length[fddf["DerivativeOrder"]] 是要逼近导数所在空间的维数. 这将返回一个长度为 n 的列表，该列表对于每一个分量维度都具有网格上的值. 由于具有 Listable 属性，作用到该列表上的函数将遍及整个网格.

定义了一个函数，利用 Exp 具备 Listable 属性来找到网格上的值：
比较三种方法所用的时间. 命令被多次重复以得到更准确的计时：

这些例子表明，用 CompiledFunction 通常比用 Apply 要快得多而利用可列表性也可以加快一些速度.

拟谱导数

差分阶数可取的最大值由网格点数决定. 如果超过极限，将给出警告信息并自动降低阶数.

使用最大阶数来逼近随机网格上正弦函数的一阶导数：

使用一个极限阶数通常被称作拟谱导数. 一个常见问题是人为振荡（Runge 现象）会非常大. 然而有两种情形除外：具有周期重复性的均匀网格，和格点在切比雪夫多项式 零点处或称为 Chebyshev–Gauss–Lobatto 点处的网格 [F96a]、[QV94]. 这两种情况均可以使用快速傅立叶变换进行计算，这种变换非常有效且能够最小化舍入误差.

"DifferenceOrder"->n使用 n 阶有限差分来逼近导数
"DifferenceOrder"->Length[grid]使用最高阶数有限差分在 grid 上逼近导数（一般不推荐）
"DifferenceOrder"->"Pseudospectral"使用拟谱导数逼近；仅适用于网格点间距与 Chebyshev–Gauss–Lobatto 点对应，或网格均匀且 PeriodicInterpolation->True 的情况
"DifferenceOrder"->{n1,n2,…}分别在维度 1, 2, … 上使用差分阶 n1, n2, …

选项 "DifferenceOrder" 的设置.

给出均匀网格上正弦函数二阶导数的拟谱逼近：
计算每一个点处的误差. 逼近精确到舍入误差，因为周期函数在均匀网格上的拟谱导数的有效基是三角函数：

Chebyshev–Gauss–Lobatto 点是 的零点. 利用性质 ，可以看出它们在 的这些点处.

定义了一个简单的函数，生成由 n 个点组成的网格，最左边的点在 x0 处而区间长度 L 具有 Chebyshev–Gauss–Lobatto 点的间距：
计算高斯函数的拟谱导数：
显示近似值与准确值的图像：
用同样点数的均匀网格在最大差分阶数下计算的导数图像：

尽管在中心处的逼近效果较好（因为均匀网格上中心处点较密集），该图还是很清楚的显示出阶数过高时有限差分逼近所常见的灾难性振荡. 使用 Chebyshev–Gauss–Lobatto 间距可以使振荡最小化.

此图显示的是利用周期性重复的均匀网格计算得到的拟谱导数逼近：

周期性假设使逼近有了明显的改进. 周期性拟谱逼近的精确度极高，于是有充分的理由在某些情况下用较大的计算域来模拟周期，比如像此例一样的脉冲. 虽然这些逼近有极高的精确度，但并非没有缺陷：最糟糕的一种是重叠误差，即一个频率过高的振荡函数分量可能被错误逼近或完全消失.

有限差分逼近的精确性和收敛性

在使用有限差分时，很重要的一点是要明白截断误差，即由截断泰勒级数逼近导致的渐近逼近误差，并不是误差的唯一来源. 在应用有限差分公式时误差还有另外两个来源：条件误差和舍入误差 [GMW81]. 舍入误差来自算术运算要求的舍入. 条件误差则来自任意函数值误差的放大，通常来自与步长的幂相除，因此随步长减小而增大. 这意味着在实践中，即使截断误差随 减小而趋向于零，实际误差将在超过某点之后开始增长. 下面的图形表示的是一个平滑函数在 变小时的典型行为.

这是在机器精度下，把高斯函数 在覆盖了区间 的网格上各点处的一阶导数近似的最大误差作为网格点数 的函数而形成的函数的对数坐标图. 在均匀网格上 2、4、6 及 8 阶的有限差分分别用红、绿、蓝及洋红色表示. 均匀（周期性）间距及切比雪夫间距的拟谱导数分别用黑色和灰色表示.

这是高斯函数 在覆盖了区间 的网格上各点处的一阶导数近似的截断误差（虚线）以及条件和舍入误差（实线）作为网格点数 的函数的对数坐标图. 在均匀网格上2、4、6 及 8 阶的有限差分分别用红、绿、蓝及洋红色表示. 均匀（周期性）间距及切比雪夫间距的拟谱导数分别用黑色和灰色表示. 截断误差通过计算非常高精度的近似值得到. 舍入和条件误差的估计值通过高精度近似值减去机器精度近似值得到. 舍入和条件误差偏向于线性增长（因为一阶导数有限差分公式常见的 因子），对于较高阶导数往往也有所增大. 拟谱导数表现更多样化这是因为 FFT 计算的误差随长度变化而变化. 请注意均匀（周期性）拟谱导数的截断误差不会降低到 以下. 这是因为数学上，高斯函数并不是一个周期函数；这个误差在本质上给出了与周期性的偏差.

这是高斯函数 在覆盖了区间 的 45 点网格上各点处的一阶导数的近似误差作为 的函数的半对数函数坐标图. 在均匀网格上 2、4、6 及 8 阶的有限差分分别用红、绿、蓝及洋红色表示. 均匀（周期性）间距及切比雪夫间距的拟谱导数则分别用黑色和灰色表示. 除了带有切比雪夫间距的拟谱导数外，其它全部用均匀间距 计算. 显然拟谱导数的误差并不是局部性的；这不奇怪因为在任何点上的近似都基于整个网格上的值. 有限差分逼近的误差是局部性的且误差量级随高斯函数的大小而变化（在半对数图形中是抛物线型）.

从第二幅图可以看出，显然存在一个能得到最近导数近似的大小；对更大的 ，截断误差起主导作用，对更小的 ，条件和舍入误差起主导作用. 最优的 对高阶差分往往也能给出较好的近似. 对于偏微分方程的空间离散化，这通常不是问题，因为达到那个精确度水平的计算代价实在太高. 然而，当使用低阶差分来逼近，比如，雅可比矩阵时，这种误差的平衡就变得至关重要了. 为避免额外的函数计算，通常使用一阶前向差分，而误差均衡与单位舍入的平方根成正比，故而挑选一个良好的 值就很重要了 [GMW81].

这些图形说明的是没有真正的边界效应的平滑函数的典型情况. 如果改变高斯函数的参数，使函数变得更平，边界效应将开始出现.

这是高斯函数 在覆盖了区间 的 45 点网格上各点处的一阶导数的近似误差作为 的函数的半对数函数坐标图. 在均匀网格上 2、4、6 及 8 阶的有限差分分别用红、绿、蓝及洋红色表示. 均匀（周期性）间距及切比雪夫间距的拟谱导数则分别用黑色和灰色表示. 除了带有切比雪夫间距的拟谱导数外，其它全部用均匀间距 计算. 有限差分逼近的误差是局部性的且误差量级随高斯函数的大小而变化. 均匀间距拟谱（45 阶多项式）逼近在边界附近的误差变得很大；而随着 的减小，此误差是无界的. 另一方面，切比雪夫间距拟谱的误差更均匀，整体相当小.

到目前为止展示的例子似乎显得逼近时的阶数越高越好. 然而还有两个额外的问题需要考虑. 逼近的阶数越高，函数求值的代价也越高，如果需要隐式迭代（比如对于刚性问题），不仅计算雅可比行列式的代价更高，矩阵的特征值也趋向于更大，从而导致更大的刚性以及迭代求解时更大的困难. 这是拟谱方法的极端情形，其中雅可比行列式基本没有非零项 [F96a]. 当然，这些问题对于较小的系统（及矩阵）是个折衷.

另一个问题与不连续性有关. 通常情况下多项式逼近的阶数越高，逼近的效果越差. 更糟糕的是对于真正的不连续点，误差会随网格间距的减小而放大.

这是不连续的单位阶越函数 f(x)=UnitStep(x - 1/2) 在覆盖了区间 的 128 点网格上各点处的一阶导数近似作为 的函数的图像. 在均匀网格上 2、4、6 及 8 阶的有限差分分别用红、绿、蓝及洋红色表示. 均匀（周期性）间距及切比雪夫间距的拟谱导数则分别用黑色和灰色表示. 除了带有切比雪夫间距的拟谱导数外，其它全部用均匀间距 计算. 它们全部表现出振荡的行为，但显然切比雪夫拟谱导数在这方面做得更好.

有许多替代方法可用于不连续点附近，比如前端跟踪. 一阶前向差分可以最小化振荡，但却引入了人为的粘性项. 比较好的替代方法是所谓的实质不振荡（ENO）方案，它全阶偏离不连续点，但在不连续点附近引入界限，限制了逼近阶数及振荡行为. 目前 ENO 方案并未在 NDSolve 中实现.

总之，选择合适的差分阶数在很大程度上取决于问题的结构. 默认选择 4 阶一般情况下对很大一部分偏微分方程都是合适的，但针对特定问题，您可能想尝试其它设置，以获得更好的结果.

微分矩阵

由于微分和有限差分逼近一样，非常自然的是线性运算，所以一种表现 FiniteDifferenceDerivativeFunction 行为的替代方式是使用矩阵. 表示微分算符近似的矩阵被称为微分矩阵 [F96a]. 虽然微分矩阵并不总是应用有限差分逼近的最佳方式（尤其是当可以用 FFT 来减小复杂度和误差时），但它们对于辅助分析以及有时作为求解偏微分方程时常常需要的线性求解程序而言，是非常宝贵的.

设 FDDF 表示 NDSolve`FiniteDifferenceDerivativeFunction[data] 对象.

FDDF@"DifferentiationMatrix"把 FDDF 作为线性运算符重新进行线性运算

形成微分矩阵.

创建 FiniteDifferenceDerivativeFunction 对象：
生成代表底层线性算子的矩阵：

矩阵以稀疏形式给出，因为一般而言，微分矩阵的非零项相对较少.

转换成常规稠密矩阵并用 MatrixForm 显示：
这三种表示方法就其对于数据的作用而言，基本上是等价的：

正如之前所提，矩阵形式对分析很有用. 例如它可以用于直接求解程序，或者用于找寻特征值，从而用在线性稳定性分析中.

计算微分矩阵的特征值：

对于可通过快速傅立叶变换计算得到的拟谱导数而言，若规模较小，使用微分矩阵可能较快，但最终，在更大的网格上，FFT 由于其较好的复杂度及数值性质表现要好得多.

对于多维导数，使用矩阵可以使其在扁平化数据，即一维导数矩阵的 KroneckerProduct 上运算. 通过一个例子来理解这一点是最容易的.

计算网格上的高斯函数，该网格是在 和 方向上网格的外积：
定义 NDSolve`FiniteDifferenceDerivativeFunction 用于通过四阶差分计算函数的混合 - 偏导数：
计算伴随的微分矩阵：

注意该微分矩阵是 1353×1353 的矩阵. 数字 1353 是张量积网格上的总点数；当然也就是 和 网格上的点数之积. 微分矩阵作用于数据向量上，该向量来自于张量积网格上的扁平化数据. 此矩阵也十分稀疏；只有大概百分之一点五的项是非零的. 这点很容易就可以从非零值的位置图上看出来.

显示微分矩阵中非零值的位置图：
比较计算混合 - 偏导数的两种方法：

该矩阵是 KroneckerProduct，即一维矩阵的直积.

得到一维微分矩阵并形成它们的矩阵直积：

使用微分矩阵会导致机器数值略有不同，因为不同的运算顺序导致了不同的舍入误差.

当想要得到导数的线性组合时，微分矩阵是非常有优势的. 例如，拉普拉斯算子的计算可以被放入单个矩阵中.

定义一个函数在张量积网格上逼近拉普拉斯算子：
计算与在 和 方向上的导数相关的微分矩阵：
将两个稀疏矩阵相加，得到拉普拉斯算子的单一矩阵：
显示微分矩阵非零值的位置图：
比较逼近拉普拉斯算子的两种不同方法得到的值和计算时间：

离散化因变量的解释

当一个因变量在 WhenEvent 中给出，在监视器（例如 StepMonitor）选项中给出，或在一个需要解释因变量的方法（例如 Projection）中给出时，对于常微分方程，一般有清晰的解释：在时间（或自变量）的某一特定值处，使用因变量所对应的那一部分解的值.

对于偏微分方程，解释并不这样明确. 从数学角度而言，在一特定时刻的因变量是空间的函数. 这导致了默认的解释，它使用 InterpolatingFunction 作为空间区域上的逼近函数来代表因变量.

对于偏微分方程的另一种可能的解释是把在一个特定时刻的因变量看成是在那个时刻空间上离散化的值的表示——也就是说既在时间上也在空间上离散化. 可以通过设置 MethodOfLines 选项的 DiscretizedMonitorVariables->True 来要求监视器和方法使用这种完全离散化的解释.

看出这两种解释之间区别的最好方法还是通过例子.

求解伯格斯方程. 设置 StepMonitor 使之能每 10 个时间步绘制一副解的图，产生渐近色彩的曲线序列. 可以把 Show 替换为 ListAnimate 来产生动画效果；注意动画中的波动并不反映真实的波动速度，因为它实际上包括了在 NDSolve 中使用的步长：

在执行上述命令时，StepMonitor 中的 u[t,x] 实际上是关于 x 的函数，所以可以用 Plot 绘制. 也可以对它进行其它运算，例如数值积分.

求解伯格斯方程. 设置 StepMonitor 使之能每 10 个时间步绘制一系列空间离散化解的图. 可以把 Show 替换为 ListAnimate 来产生动画效果：

在这一例中，每一步 u[t,x] 都作为空间网格上解的离散值向量给出. 在这个例子里显示离散化的点成了一个提供信息更丰富的监视器，因为可以看清楚形成时前部的分裂情形.

值向量不含有网格自身的信息；在这个例子中，图像时相对于指标值绘制，显示的是均匀网格的正确间距. 注意当 u 被解释为函数时，网格将被包含在用于表示空间解的 InterpolatingFunction 中，所以若需要网格，得到它最容易的方法是从表示 u[t,x] 的 InterpolatingFunction 中提取.

最后要注意的是离散化表示要快得多. 若在 WhenEvent 中或在如 Projection 这样的解方法中使用该表示时，这可能很要紧. 在一个使用事件检测防止解集超越计算区域的例子中，离散化的解释使计算大大加快.

边界条件

通常在偏微分方程领域，对于特定的方程和边界条件，是有可能找到对于该边界条件适用的数值方法的. 之前“数值直线方法”引言中给出的例子就是如此. 然而，找到一般性的算法要困难得多，而且边界条件对刚性和整体稳定性的影响也让问题更加复杂.

周期性边界条件特别易于处理：周期性插值被用于有限差分. 由于均匀网格的拟谱逼近是精确的，通常可以很有效找到解.

NDSolve[{eqn1,eqn2,…,u1[t,xmin]==u1[t,xmax],u2[t,xmin]==u2[t,xmax],…},{u1[t,x],u2[t,x],…},{t,tmin,tmax},{x,xmin,xmax}]
求解函数 u1, u2, … 的偏微分方程组，其中空间变量 x 带有周期性边界条件（假设 t 是时间变量）
NDSolve[{eqn1,eqn2,…,u1[t,x1min,x2,…]==u1[t,x1max x2,…],u2[t,x1min,x2,…]==u2[t,x1max x2,…],…},{u1[t,x1,x2,…],u2[t,x1,x2,…],…}, {t,tmin,tmax},{x,xmin,xmax}]
求解函数 u1, u2, … 的偏微分方程组，其中空间变量 x1 带有周期性边界条件（假设 t 是时间变量）

给定偏微分方程的边界条件.

如果您要求解若干个函数 u1、u2、……，则当任一函数具有周期性边界条件时，它们全体都必须有（只需要为一个函数指定条件）. 如果空间维度大于一，则可以只在一些自变量上有周期性边界条件而其它变量没有.

用拟谱方法求解带周期性边界条件的二维空间上的广义正弦-Gordon 方程. 如果不用周期性所允许的拟谱法，求解可能需要更长时间：

在作为解返回的 InterpolatingFunction 对象中，记法 {…,xmin,xmax,…} 中的省略号表示表示这一维度的周期性重复.

绘制从 处的周期性延续中推导得到的部分解的曲面图：

NDSolve 使用了两种方法来处理非周期性边界条件. 两者都有其优点和缺点. 第一种方法是相对于时间变量对边界条件微分，并在边界处求解由此产生的微分方程或方程组. 第二种方法是离散化每一个边界条件. 这通常会生成一个边界解的分量的代数方程，因此必须用微分代数方程求解程序求解. 这是用 DifferentiateBoundaryConditions 的选项 MethodOfLines 控制的.

要了解微分法的工作原理，请再次考虑直线法引言部分中的简单例子. 在第一个表达式中， 处的狄利克雷边界条件是通过对 微分来处理的. 诺伊曼边界条件则用反射的想法来处理，该想法对于二阶有限差分逼近来说是可行的，但并不那么容易推广到更高阶（虽然对该问题可以很容易的通过计算整个反射来完成）. 然而微分方法可用于在 处诺伊曼边界条件的任意阶微分. 作为例子该问题的解可以用四阶微分求得.

设置空间点数量和之间间距. 故意设置得很小是为了让您看到结果方程式. 您可以以后修改它来提到逼近的精确性. sf 是边界条件的比例因子：
定义向量 ：
在 处离散化边界条件很简单就是 . 微分该方程并加入到条件中产生下式：
微分此方程并加上比例因子 sf 乘以离散化边界方程，可以得到新方程：

该方程的优点是包含显式的导数 u0'[t]，可以把离散化系统作为一组常微分方程来处理. 通过加上离散化边界条件，该方程的解本身就渐进收敛于正确的边界条件，即便有暂时的偏离或初始条件不一致. 这个可以从它的精确解中看到.

找到导出边界方程在 处的精确通解：

请注意比例因子的值将决定实际边界条件的解的收敛速率. 该因子可以用 NDSolve 中的选项 "DifferentiateBoundaryConditions"->{True,"ScaleFactor"->sf} 设置. 在 "ScaleFactor" 的默认值为 Automatic 的情况下，狄利克雷边界条件所用的比例因子是 1 而其它情形则为 0.

在空间方向上离散化 处的诺伊曼边界条件：
求离散化边界条件关于 的微分并用来替代离散化边界条件：

从技术上讲，离散化边界条件没有必要和其余微分方程有相同的差分阶数；事实上，由于单边导数的误差项要大得多，在边界处附近增加阶数往往是有必要的. NDSolve 并没有这样做因为微分阶数和 InterpolatingFunction 的插值阶数在空间方向上一致是有可能的.

另一种用 NDSolve`FiniteDifferenceDerivative 生成方程的办法. 第一个和最后一个方程将被替换为源于边界条件的适当方程：
现在，可以用边界条件替换第一个和最后一个方程：
NDSolve 可以像求解适当的导数一样求解该方程组，因此现在可以求解这些常微分方程：
这幅图形显示的是边界条件在 处的满足程度：

把边界条件当作代数条件可以在使用微分代数方程求解程序时节省好几步处理过程.

用边界条件对应的代数条件替换（之前的）第一个和最后一个方程：
用 NDSolve 求解微分代数方程组：
显示边界条件在 处的满足程度：

对于这个例子，边界条件在两种情况下的容差范围内均得到很好的满足，微分法做得稍好一点. 但并非总是如此；当边界条件被微分时，解可能会局部的偏离实际的边界条件.

绘制图形比较两种方法中狄利克雷边界条件在 处的满足程度. 由微分边界条件得到的解用黑色表示：

使用 NDSolve 时，很容易通过 DifferentiateBoundaryConditions 选项在两种方法之间切换. 应该记住的是当使用 DifferentiateBoundaryConditions->False 时并不能自由选择积分方法；所需方法必须是一个微分代数方程求解程序.

对于偏微分方程组或边界条件更加复杂的具有高阶导数的方程组，两种方法一般都可以工作良好.当在一端有多个边界条件时，可能需要在内点处附加一些条件. 这里有个例子，偏微分方程在空间区间的每一端都具有两个边界条件.

求解在空间区间的每一端都具有两个边界条件的微分方程. 使用 StiffnessSwitching 积分方法是为了避免四阶导数稳定性的潜在问题：

如何理解空间误差的提示信息将在下一节中介绍. 现在请忽略此信息而只考虑边界条件.

形成 InterpolatingFunction 列表，微分阶数与每一个边界条件相同：
绘制四个边界条件对 NDSolve 计算得到的解的满足程度作为 的函数的对数图像：

很明显边界条件在 AccuracyGoal 和 PrecisionGoal 选项允许的公差范围内得到了很好的满足. 通常，具有高阶导数的条件不会像具有低阶导数的条件那样好满足.

对于带有空间导数的边界条件，乘以原始边界条件的比例因子为零的原因有两个. 第一，强加条件到离散方程只是一个空间近似，所以任何时候都让它尽可能精确并不是总有意义. 第二，当包括条件时，特别对于高阶空间导数，单边有限差分的大系数可能是潜在的不稳定源. 上面的例子中，在 的边界条件上就发生了这样的情况.

求解所有边界条件比例因子为 1 的微分方程：
绘制四个边界条件对 NDSolve 计算得到的解的满足程度作为 的函数的对数坐标图. 请注意最后一个的指数增长是不稳定的信号：

不一致的边界条件

有一点很重要，您指定的边界条件要和初始条件及偏微分方程都一致. 若不是这种情况，NDSolve 会发出关于不一致性的警告信息. 这种情况下解可能不满足边界条件，甚至更糟的是可能会出现不稳定性.

在这个热传导方程的例子中，在 处的边界条件明显与初始条件不一致：
这里绘制了 处的解作为 的函数时的图像. 边界条件 并未被满足：

边界条件没有满足的原因是因为常微分方程的积分起始于给定的初始条件. 在边界的微分方程一旦被微分，就变成了 ，于是解接近正确的边界条件但在积分区间中并不是非常接近.

如果边界条件不被微分，则微分代数方程求解程序实际上会修改初始条件来满足边界条件：

微分代数方程求解程序并不总能找到产生像这样高效正确解的一致初始条件. 一种改进用常微分方程方法求解的方法是在 "DifferentiateBoundaryConditions"->{True,"ScaleFactor"->sf} 中使用比默认值 1 更高的边界条件比例因子值.

在这个热传导方程的例子中，在 处的边界条件显然与初始条件不一致：
这里绘制了在 处的解作为 的函数与边界条件 的值的差异. 显然随着比例因子的增长，边界条件越来越快的在公差范围内被满足：

当比例因子很大时，求解程序更难对方程积分因为它引入一个可能的刚度源. 可以看到在比例因子为 100 的解中，就有一个变化很快的初始瞬间.

这里显示了比例因子为 100 的解的曲面图：

更保险的处理这个问题方法是给出一个与边界条件一致的初始条件，即使它不连续也行. 在这种情况下，单位阶跃函数足矣.

这里用了不连续的初始条件来匹配边界条件，给出了相对于空间离散化的分辨率而言正确的解：

一般而言，对不连续的初始条件，空间误差估计是不能被满足的. 由于它们是根据平滑度估计的，因此一般情况下，最好是对您所希望模拟的不连续性效果进行选择，要么给出逼近不连续性的平滑函数，要么显式指定空间离散化所用的点的数目. 关于空间误差估计和离散化的更多细节可参见“空间误差估计”.

一种更微妙的不一致出现在当时间变量有更高阶的导数而边界条件可能被多次微分时.

考虑波动方程

初始条件 满足边界条件，所以您可能对 NDSolve 会报告 NDSolve::ibcinc 信息而感到惊讶.

在这个例子中，边界和初始条件初看上去似乎是一致的，但实际上在微分时会表现得不一致：

不一致性出现在对第二个初始条件做关于 的微分时，得到 ，而对第二个边界条件做关于 的微分时，得到 . 在 时这两者不一致.

偶尔，NDSolve 会在边界条件实际上不一致时给出 NDSolve::ibcinc 这样的警告信息. 这归因于在逼近诺伊曼边界条件或任何涉及到空间导数的边界条件时的离散化误差. 原因是用于确定有多少点需要离散化的空间误差估计（参见“空间误差估计”）是基于偏微分方程和初始条件，而不是边界条件. 用于逼近边界条件的单边有限差分公式比相同阶数的中心公式有更大的误差，从而导致边界处额外的离散化误差. 通常来说这不是问题，但构造例子使其发生是可能的.

在这个例子中，因为离散化误差，NDSolve 错误的发出不一致边界条件的警告信息：
边界条件的图像表明误差虽然不大，却在默认的公差范围之外：

当边界条件一致时，修正此误差的方法是指定 NDSolve 使用更精细的空间离散化.

用更精细的空间离散化，没有警告信息而且边界条件被更好的满足：

空间误差估计

概述

当 NDSolve 求解偏微分方程时，除非已经以显式给出或给出 MinPoints 和 MaxPoints 选项相等值的方式指定了它的空间网格，否则 NDSolve 需要做空间误差估计.

理想情况下，空间误差估计会随着时间推移而被监测，而空间网格则会随着解的演变而更新. 网格自适应性问题对于指定类型的偏微分方程是极为困难的，当然也没有一般性的求解方法. 而且，诸如局部细化等技术用于直线方法的时候可能是有问题的，因为改变网格点的数量需要完全重新启动常微分方程方法. 移动网格技术对这种方法而言似乎是有前途的，但目前，NDSolve 使用的是静态网格. 所用的网格由基于初始条件的先验误差估计确定. 在时间区间的末端会进行一个后验检查来考察一致性是否合理，如果检查没有通过则发出警告信息. 这一点当然可以被愚弄，但实际上它提供了合理的折衷. 失败最常见的原因是是在初始条件没有多大变化时，所以估计本质上是没有意义的. 在这种情况下，可能需要自行选择合适的网格设置.

加载用于从 InterpolatingFunction 对象中提取数据的软件包：

先验误差估计

当 NDSolve 利用直线方法求解偏微分方程时，需要在合适的空间网格上作出决定. NDSolve 用基于初始条件（也就是先验条件）的误差估计来做到这一点.

最容易的是通过例子来说明其工作原理. 出于演示目的，考虑带周期性边界条件的一维正弦-Gordon 方程.

求解带高斯初始化条件的正弦-Gordon 方程：
分别给出所用的空间和时间点的数目：

时间点是由常微分方程方法基于局部误差控制而自适应的选择的. NDSolve 用了 97（如果包括周期性的端点就是 98）个空间点. 这个选择将在后续步骤中演示.

在 NDSolve 的方程处理阶段，首先发生的事情之一是有二阶（或高阶）时间导数的方程被只有一阶时间导数的方程组所替代.

这是与先前正弦-Gordon 方程等价的一阶方程组：

下一阶段是求解时间导数.

这是时间导数的解，方程右手侧是正常（常微分方程）形式：

现在问题是选择一个均匀网格来逼近导数，局部误差的容差由 AccuracyGoal 和 PrecisionGoal 指定. 对于这个例子，使用默认的 "DifferenceOrder"（4）和默认的 AccuracyGoal 和 PrecisionGoal（对于偏微分方程都是 4）. 用于积分从离散化得到的常微分方程组的方法，其自身的误差估计是建立在有足够准确的函数值的假设上的. 这一估计的目标在于寻找一个空间网格，使得空间误差与局部时间误差（至少在初始条件下）能够保持平衡.

这里将变量设置得能够看出 "DifferenceOrder"、AccuracyGoal 和 PrecisionGoal 的默认设置：

误差估计基于理查德森外推法. 如果知道误差是 且在 的两个不同值 和 处有两个近似值 和 ，则事实上可以外推至极限 而得到一个误差估计

因此 的误差估计为

这里 和 是长度不同的向量，而 是个函数，故而需要选择一个适当的范数. 如果选择了 ，那么可以简单的在两个向量公有的分量上使用一个缩放后的范数，即所有的 和 的所有间隔点. 这是一个好选择，因为它不需要任何网格间的插值.

若要在一个给定区间上设置均匀网格，则可以定义一个函数 ，其中 是使网格为 的区间长度，而 .

为正弦-Gordon 方程定义返回步长 和网格点列表的 的函数：

对给定的网格，方程可以用有限差分离散化. 这点可用 NDSolve`FiniteDifferenceDerivative 轻松做到.

这里定义了正弦-Gordon 方程右手侧的符号离散化作为网格的函数. 它返回一个关于 和 的函数，给出 和 近似值的列表.（请注意原则上这对于任何网格均适用，不论是否均匀，尽管下面仅使用了均匀网格.）

对给定的步长和网格，也可以离散化 和 的初始条件.

这里定义了一个函数用于离散化 和 的初始条件. 最后一个网格点被丢弃，这是因为周期性延续下它被认为与第一个网格点相同：

我们感兴趣的量是在此初始条件下， 取特定值时右手侧的近似.

这里定义了一个函数，在给定的步长和网格的初始条件下，返回一个由方程右手侧近似值组成的向量. 向量是扁平化的，以简化后续分析：

从 的特定值开始，可以通过生成 和 个点的右手侧来获得误差估计.

这里给了一个 10 点网格的右手侧近似向量的例子：
这里给了一个 20 点网格的右手侧近似向量的例子：

如之前提到的，在 个点的网格上每隔一个点都在 个点的网格上. 因此为简单起见，可以使用一个范数仅比较两个网格的公共点. 因为目标最终是满足绝对和相对的公差标准，所以使用缩放的范数是合理的. 除了要考虑到缩放时右手侧的大小，包括网格上 和 相应分量的大小也很重要，因为右手侧的误差最终会包括在 和 中.

这里对右手侧逼近的差分定义了一个范数函数：
把范数函数应用于得到的两个近似值：

根据理查德森外推法公式（3），为从距离得到误差估计，只需简单将此值除以 .

计算 时的误差估计. 因为它基于缩放后的范数，所以结果小于 1 即满足公差标准：
创建了一个结合了之前函数以给出误差估计的 的函数：

我们的目的是为了找到 的最小值，使得误差估计小于或等于 1（因为是基于缩放后的范数）. 原则上，在这上面可能使用求根算法，但因为 只能是整数，求根算法可能过于复杂且必须调整终止条件. 更简单的解决方案是用理查德森外推法公式预测什么样的 值是恰当的并重复预测过程直至找到合适的 .

需要满足的条件是

而已经估计到

因此可以预测

或者以 的形式，它与 的倒数成比例

基于之前 时计算的误差估计来计算预测最佳 值：
计算 的新值的误差估计：

通常情况下基于非常粗糙的网格的预测是很不充分的. 粗糙的网格可能会完全丢失某些解特征，从而在更精细的网格上导致误差. 此外，由于误差估计是基于渐近公式的，因此对于粗糙间距，估计本身可能不是很好，即使在所有解特征都在一定程度上得到分辨的情况下也是如此.

实践中，计算这些误差估计可能相当耗费资源. 此外也没有必要非找到最佳的 ，能满足误差估计即可. 记住，一切都可能随偏微分方程的演变而改变，因此根本不值得只为初始时间寻找最佳间隔而付出太多额外努力. 一个简单的解决方案是在新 的预测公式中包含一个额外的大于 1 的因子. 即使有额外的因子，可能仍然需要进行一些迭代来得到一个可接受的值. 然而它确实在通常情况下能加快进程.

定义了一个函数，给出满足误差估计的，网格点数目的预测值：
对预测进行迭代直至找到满足误差估计的值：

重要的是要注意，此过程必须有一个限制值，因为它可能无法满足误差的公差，比如在有个不连续的初始函数时. 在 NDSolve 中，MaxSteps 选项提供了这个限制值；对空间离散化，默认在所有空间维度上的总值为 10000.

拟谱导数不能使用此项误差估计，因为它们是指数收敛而不是多项式收敛. 基于之前用的公式在 p->Infinity 时可以对误差做一个估计. 这么做事实上认为结果在更精细的网格上是精确的且误差估计是基于差分的，因为 趋向于 . 更好的方法是利用这么一个事实：在给定的 个点的网格上拟谱方法的复杂度是 . 在比较两个网格时，对 使用较小的 是合适的. 就确定网格大小这一目的而言，这种方法提供了一个不完美但却充分的估计.

修改误差估计函数使之能用于拟谱导数：

类似的，预测公式也可以改用 代替 .

修改函数预测适用于拟谱导数的适当的 的值. 这一公式化并不试图挑选一个高效的 FFT 长度：

在最终确定拟谱方法所用的 时，还要考虑的是选择一个不仅仅满足公差条件，更是能计算 FFT 高效率长度的值. 在Wolfram 语言中，高效率的 FFT 不要求长度为二的幂次，因为 Fourier 有一个内置的素因子算法.

通常差分阶数对满足误差估计所需要的点数有着深刻的影响.

创建一个表格，表示的是满足先验误差估计所需要的点数作为差分阶数的函数：

这一表示需要的点数作为差分阶的函数的表格对解释为什么直线方法默认设置 "DifferenceOrder"->4 大有帮助. 从 2 到 4 的改进通常是最显著的，且在默认的公差范围内四阶差分不大可能产生在更高阶时可能出现的大的舍入误差. 拟谱差分常常也是个好选择，特别是带有周期性边界条件时，它们不是好的默认设置因为它们会导致满雅可比矩阵，这对于可能的刚性方程的生成和求解都非常耗资源.

对非周期性网格，仅用内点即可完成误差估计. 原因是边界附近的导数的误差系数是不同的. 这可能会丢失一些边界附近的特征，但主要的想法是使估计较为简单易行，因为无论如何偏微分方程的发展可能会改变一切.

对多个空间维度，一次只对一个维度做出决定. 因为一个维度上更好的决策可能会改变对另一个维度的要求，整个过程逆序重复以改进选择.

后验误差估计

当用 NDSolve 计算偏微分方程的解时，最后一步是对演化后的解做空间误差估计，若误差太大则发出警告消息.

这些误差估计的方式与之前描述的先验估计类似. 唯一真正的区别是估计误差不用 和 个点的网格而是 和 个点的网格. 这是因为在 个点的网格上生成值必须用插值而插值自身会引入误差，在 个点的网格上的值是现成的，间隔取值即可. 令 ，很容易根据理查德森外推公式得到

定义了一个函数（基于前一节定义的函数），为 和 的解计算关于正弦-Gordon 方程的解的误差估计，以向量的形式. 它被定义为网格的函数因为它作用在一个已经构建好的网格上.（注意正如这里定义的那样，这仅适用于偶数长度的网格. 处理奇数长度并不困难，但多多少少会使函数更加复杂.）
求解带高斯初始条件的正弦-Gordon 方程：
这是用于空间方向的网格，是用在 InterpolatingFunction 中的第一组坐标. 由于周期性延续，用于得到值的网格丢弃了最后一个点：
创建一个函数，给出 取某一特定数值时的后验误差估计：
绘制后验误差估计作为 的函数的图像：

在这个函数中看到的大量的局部变化是很典型的.（在这一例中，这个大的峰值正好出现在原来的单峰分裂成独立的波时.） 由于这个原因，NDSolve 并没有警告误差过高，除非误差超过 10（而不是 1，那是用于选择基于初始条件的网格的）. 后验误差估计不如先验误差估计准确这一事实进一步证明了额外因子 10 的合理性. 因此，当 NDSolve 基于后验误差估计发出警告信息时，通常是因为解的新特征已经出现，或因为求解过程中的不稳定性.

这个例子用的是与之前例子相同的初始条件，但 NDSolve 基于后验误差估计给出警告信息：
这里显示的是 时解的图像. 显然警告信息是合理的，因为峰值附近的振荡是非物理的：

如果出现 NDSolve::eerr 信息，则可能有必要用选项来控制网格选择过程，因为很可能默认设置并没有找到准确的解.

控制空间网格选择

直线方法的 NDSolve 实现有若干种方式来控制空间网格的选择.

选项名称
默认值
AccuracyGoalAutomatic确定网格间距的绝对公差位数
PrecisionGoalAutomatic确定网格间距的相对公差位数
"DifferenceOrder"Automatic用于空间离散化的有限差分逼近阶数
CoordinatesAutomatic自变量 在各个空间维度上的列表 这将覆盖此列表中随后所有选项的设置
MinPointsAutomatic网格中每个维度上的最少点数；设为 Automatic 时，值将通过给定差分阶数下计算误差估计所需的最少点数确定
MaxPointsAutomatic用于网格中的最大点数
StartingPointsAutomatic用先验误差估计开始网格细化的点数
MinStepSizeAutomatic所用最小网格间距
MaxStepSizeAutomatic所用最大网格间距
StartingStepSizeAutomatic用先验误差估计开始网格细化的网格间距

直线方法的张量积网格选项.

所有用于张量积网格离散化的选项可以用一个列表给出，列表长度等于空间维数，在这种情况下，各个空间维度的参数由列表中相应的分量确定.

除了非周期性问题的拟谱方法是个例外，离散化都是用均匀网格实现的，所以在区间长度 上求解问题时，在 Points 和 StepSize 选项之间有直接的对应关系.

StepSize 选项被有效转化为等价的 Points 值. 这样只是为了提供方便，因为有时问题的设定与步长能更自然的联系在一起，而不是与离散点的个数. 当为 Points 和对应的 StepSize 选项指定除了 Automatic 之外的值时，一般来说限制条件更严格.

在之前一节所示的例子中解的分辨率并不充分，这是因为解的图像逐渐变陡. 下面的例子将展示一些修改网格参数的不同方法，以使接近激波的解能得到更好的分辨.

一种避免这种解随图线变陡而出现振荡的方法是确保在图线最陡的地方使用足够多的点来分辨解. 在伯格斯方程

的单峰解中，可以证明 [W76] 在 时激波宽度与 成正比. 超过 95% 的变化包含在一个宽度为 的层内. 因此如果选择一半激波宽度的最大步长，总是会有一点在激波的陡峭部分，因而有希望在没有显著振荡的条件下对其进行分辨.

计算伯格斯方程的解，使得有足够多的点来分辨激波：

注意单单分辨波形绝不足以满足 NDSolve 的默认公差，它的精确度需要达到 . 然而，一旦有足够多的点来分辨基本图线，根据 NDSolve::eerr 信息中的后验误差估计进行预测并无不妥（可以再增加 10%，因为说到底它只是预测）.

计算伯格斯方程的解，最大步长选择得足够小，以满足基于先前计算的误差预测的默认公差：

要比较像这样的解，观察空间格点处解的图像就非常有用. 因为网格点已被作为 InterpolatingFunction 的一部分存储，定义这样一个函数相当简单.

定义了一个函数，仅绘制时刻 时在空间网格点上的解：
绘制图像比较在 时的三个解：

在这个例子里，区域的左手侧真的不需要这么多点. 点需要在图线陡峭的地方聚集，因此考虑显式指定一个在陡峭处有更多点的网格是合理的.

在指定的，大部分点在 右侧的网格上求解伯格斯方程：
在这个指定的空间网格上绘制解的值：

许多同样的原则适用于多个空间维度. 在两个维度上各向异性的伯格斯方程提供了一个好例子.

求解在 和 方向上速度不同的二维伯格斯方程的变种：
这是解在 时的前缘曲面图：

与一维情形相似，前缘变得陡峭. 由于速度项（）较大，陡峭显得并不怎么极端，这个默认的解实际上对前端分辨得相当好. 因此应该可以从误差估计来预测以满足默认公差. 一个简单的缩放参数表明，图线宽度在 方向上会比在 方向上窄 倍. 因此 方向的步长比 方向的步长大这样一个因数是合理的，或者最小点数可以相应的小上 .

求解伯格斯方程的二维变种，它在 和 方向上具有根据之前计算的后验误差估计预测的适当步长限制，在 方向上是 69 个点：

这个解用了大量时间计算，这不奇怪因为这个解涉及到求解 18000 个常微分方程组成的方程组. 在许多情况下，特别是在一个以上的空间维度上，满足默认公差可能是不切实际的，所以您可能需要适当使用 AccuracyGoal 和/或 PrecisionGoal 来降低公差要求. 有时，特别是当网格因不太严格的公差而较粗糙时，该方程组不是刚性的，可能使用显式方法避免可能有问题的数值线性代数，尤其在处理高维问题时. 对这个例子，用 Method->ExplicitRungeKutta 求解只需要大约一半的时间.

网格的其它任何选项都可以用列出每个维度上值的列表指定. 当只给出一个值时，它将用于所有空间维度. 有两个例外，一个是 MaxPoints，单个值被当作外积网格的总点数，另一个是 Coordinates，必须在每个维度上指定网格.

这里选择了之前解的部分网格，在前端较陡峭的部分网格更密集：

需要记住的很重要的一点是，空间的后验误差估计只是计算空间导数时对局部误差的估计，可能无法反映给定解实际的空间累积误差. 一种得到实际空间误差估计的方法是计算在非常严格的时间容差下不同空间网格上的解. 为了说明这一工作原理，再次考虑更简单的一维伯格斯方程.

用 {33,65,…,4097} 的空间网格点计算伯格斯方程在差分阶为 2、4、6 和拟谱时的解的列表. 时间准确度和精度公差都设得非常高以使几乎全部误差都来自于空间离散化. 请注意，通过在 NDSolve 中指定 {t,4,4}，只有在 时的解被保存. 若没有这种预防措施，在更精细的网格上（需要更多的时间步）的某些解就可能耗尽可用内存. 即便如此，解的列表还是需要的大量时间来计算：

给了两个解，就需要在两者之间作比较. 为了排除解自身以外的所有误差来源，最好使用用于创建 InterpolatingFunction 的插值数据. 这可以通过使用两个解的公共点来完成.

定义了一个函数通过在公共点上比较两个不同的解来估计误差. 参数 coarse 和 fine 分别对应粗糙和精细网格上的解. 这适用于先前网格间距以 2 的幂变化时所生成的解：

为了得到误差的一般趋势的指标（在不稳定的情况下，解不收敛，所以不考虑这种情况）， 可以比较连续的解对的差异.

定义了一个函数，绘制对给定差分阶数时得到的连续的解的误差估计序列的图像，并用它作出估计误差作为网格点数量的函数的对数坐标图：

这是伯格斯方程在 时的近似解的最大空间误差作为网格点数的函数的对数坐标图. 在均匀的网格上阶数为 2、4 和 6 的有限差分分别显示为红色、绿色和蓝色. 均匀（周期性）间距的拟谱倒数显示为黑色.

图像左上角是图线没有充分分辨的网格，所以差异的数量级只是 1（如果有不稳定性情况会更糟）. 然而，一旦有足够数量的点来分辨没有振荡的波形，收敛将非常迅速. 毫不奇怪，对数直线的斜率为 ，这对应于 NDSolve 默认使用的差分阶数. 如果网格足够精细从而进入渐近收敛的部分，可能像在前面两节中那样通过使用理查德森外推法作简单的误差估计，但这是在整体解上，而不是局部误差上. 另一方面，计算更多的值并观察图像能够指示是否在渐近部分.

从图像中显而易见，最好的解是 2049 个点的拟谱解（具有更多点数的解没有计算，因为其空间准确度远远超过了设定的时间公差）. 这个解实际上几乎可视为精确解， 至少在 左右的公差范围内是这样.

要判断如何最佳的求解这个问题，进行下列操作是非常有用的：对每一个求得的，至少存在一个合理逼近的解，用和解可能的空间准确性差不多的时间准确性公差重新计算，并将得到的准确度作为求解时间的函数绘制图像. 下面（有点复杂）的命令可以完成这些.

这里确定了将要使用的“最佳”解，实际也是后面计算的精确解. 和它比较的列表中去掉了它因为与自身的比较没有意义：
定义了一个函数，给定差分阶数 do 和利用该差分阶数计算得到的解 sol 后重新计算该解，其中局部时间公差比所得到的实际空间准确度稍微严格一些，如果该准确度已足够的话. 函数输出的是列表 {网格点数, 差分阶数, 计算所用秒数, 重新计算的解的实际误差}：
将函数应用于之前计算的每一个解.（差分阶数要合适！）
去除没有重新计算的例子，并绘制准确度作为计算时间的函数的对数坐标图：

表示伯格斯方程在 时的近似解的误差作为计算时间的函数的对数坐标图. 显示的每个点标明了计算解所用的空间网格点数. 在均匀网格上阶数为 2、4 和 6 的有限差分分别显示为红色，蓝色和绿色. 具有均匀（周期性）间距的拟谱导数用黑色显示. 注意拟谱方法的成本从 513 大幅跳跃到了 1025. 这是因为该方法已切换到了刚性求解程序，它因离散化产生的密集雅可比矩阵而变得耗费巨大.

所得的图像相当有力的表明，在适用的，比如这种情况下周期性拟谱逼近是惊人的高效. 不适用的情况下，在某种程度上，差分阶数越高，逼近效果一般越好. 这些都是平滑问题的特征，伯格斯方程就是一个特定的实例. 然而在极限 时，更高阶的解一般都非常差.

最后要注意的一点是，上图在时间方向上的计算用到了 Automatic 方法. 这用到了 LSODA，它根据解的演变方式在刚性和非刚性方法之间切换. 对于粗糙网格，严格显式的方法通常稍快，并且除了拟谱情况，隐式 BDF 的方法对于较精细的网格更快. NDSolve 中有多种常微分方程的替代方法可用.

边界处的误差

先验误差估计是在计算区域内部计算的，因为那里所用的差分都具有一致的误差项，可以用来有效地估计要用的点数. 将边界包括在误差估计内将使这一过程变得更为复杂，超出了这种先验估计所能保证的合理性. 通常情况下，这种方法能成功的把误差控制在合理范围内. 然而，有一些情况下可能会导致困难.

偶尔可能会因为边界处所用单边导数的误差项较大，导致 NDSolve 检测到边界条件和初始条件之间的不一致性，这是离散误差的假缺陷.

求解左端恒温，右端向自由空间辐射的一维热传导方程：

在这种情况下，NDSolve::ibcinc 发出提示信息完全是因为在右边界有较大的离散误差. 对于这个特定的例子，额外的误差并不是一个问题，因为误差会因为方程的性质而逐渐减弱并消失. 当然，多用几个空间点也可以使这则消息消失.

计算与上面方程相同的解，但在 方向至少用 50 个点：

另一种边界处的误差问题可能意外影响到离散化的情形是当周期性边界条件由并不是真正的周期函数给出时，这时一个无意的不连续就被引入到了计算中.

开始计算带高斯初始条件和周期性边界条件的正弦-Gordon 方程的解. NDSolve 命令被包在 TimeConstrained 里，因为求解给定的这个问题可能需要很长的时间和大量的系统内存：

这里的问题是，当考虑周期性连续时，初始条件实际上是不连续的.

这里显示的是三个完整周期上初始条件的图像：

由于在尖点处总是有一个大的导数误差，NDSolve 被迫使用最大点数试图满足先验误差界限. 更糟糕的是，极端的变化使求解由此产生的常微分方程更加困难，导致极长的求解时间和大量的内存消耗.

如果确实需要不连续性，那么通常要指定空间网格中若干个点或间距，使之足以处理兴趣范围内不连续性的各个方面. 模拟高精度的不连续性通常需要专门的方法，这超出了 NDSolve 所提供的一般方法的范围.

另一方面，如果确实不需要不连续性，如这个例子中不连续只是由于选择的计算域过小而无意造成，这时问题通常可以通过延长区域或增加项使周期之间平滑过渡而轻松解决.

求解正弦-Gordon 问题，计算域足够大所以初始条件中的不连续性与默认公差所允许的误差相比可以忽略不计：