16

（注意：这是一个社区 Wiki。）

假设我有一组点xi = { x0 , x1 , x2 ,... xn } 和相应的函数值fi = f(xi) = { f0 , f1 , f2 ,..., fn }，其中f ( x ) 通常是一个未知函数。（在某些情况下，我们可能会提前知道f ( x )，但我们通常希望这样做，因为我们通常不提前知道f ( x )。）什么是近似f (x ) 在每个点xi ? 也就是说，我如何估计每个点xi的dfi == d/d x fi == d f ( xi )/d x的值？

不幸的是，MATLAB 没有很好的通用数值微分例程。造成这种情况的部分原因可能是因为选择一个好的例程可能很困难！

那么有哪些方法呢？存在哪些套路？我们如何为特定问题选择一个好的例程？

在选择如何在 MATLAB 中进行微分时有几个注意事项：

  1. 你有一个符号函数还是一组点？
  2. 你的网格是均匀的还是不均匀的？
  3. 您的域是周期性的吗？你能假设周期性边界条件吗？
  4. 您在寻找什么级别的准确度？您是否需要在给定的公差范围内计算导数？
  5. 在定义函数的相同点上评估导数对您来说是否重要？
  6. 您是否需要计算多阶导数？

最好的方法是什么？

4

2 回答 2

20

这些只是一些简单粗暴的建议。希望有人会发现他们有帮助！

1.你有一个符号函数还是一组点？

  • 如果你有一个符号函数，你可以分析地计算导数。（很可能，如果它那么容易，你就会做到这一点，而且你不会在这里寻找替代品。）
  • 如果您有一个符号函数并且无法解析地计算导数，您始终可以在一组点上评估该函数，并使用此页面上列出的其他方法来评估导数。
  • 在大多数情况下，您有一组点 (xi,fi)，并且必须使用以下方法之一....

2. 你的网格是均匀的还是不均匀的？

  • 如果您的网格是均匀分布的，您可能需要使用有限差分方案（请参阅此处或此处的维基百科文章），除非您使用周期性边界条件（见下文）。这是在求解网格上的常微分方程的背景下对有限差分方法的一个不错的介绍（尤其参见幻灯片 9-14）。这些方法通常计算效率高，易于实现，并且该方法的误差可以简单地估计为用于推导它的泰勒展开式的截断误差。
  • 如果您的网格间隔不均匀，您仍然可以使用有限差分方案，但表达式更加困难，并且精度随网格的均匀程度而变化很大。如果您的网格非常不均匀，您可能需要使用较大的模板尺寸（更多相邻点）来计算给定点的导数。人们经常构造一个插值多项式（通常是拉格朗日多项式）并对该多项式进行微分以计算导数。例如，参见这个StackExchange 问题。使用这些方法通常很难估计误差（尽管有些人试图这样做：这里和这里）。Fornberg 方法在这些情况下通常非常有用....
  • 必须注意域的边界，因为模板通常涉及域外的点。有人引入“鬼点”或将边界条件与不同阶的导数结合起来，以消除这些“鬼点”，简化模板。另一种方法是使用右侧或左侧有限差分方法。
  • 这是一个极好的有限差分方法“备忘单”，包括低阶的居中、右侧和左侧方案。我在我的工作站附近保留了一份打印输出，因为我发现它非常有用。

3. 你的域名是周期性的吗？你能假设周期性边界条件吗？

  • 如果您的域是周期性的，您可以使用傅里叶谱方法计算导数到非常高的阶精度。这种技术在一定程度上牺牲了性能以获得高精度。事实上，如果您使用 N 个点，您对导数的估计大约是 N^th 阶准确的。有关更多信息，请参阅（例如）本 WikiBook。
  • 傅里叶方法通常使用快速傅里叶变换 (FFT) 算法来实现大致 O(N log(N)) 的性能，而不是简单实现的离散傅里叶变换 (DFT) 可能采用的 O(N^2) 算法。
  • 如果您的函数和域不是周期性的，则不应使用傅里叶谱法。如果您尝试将其与非周期性函数一起使用，则会出现大错误和不希望的“振铃”现象。
  • 计算任何阶的导数需要 1) 从网格空间到光谱空间的变换 (O(N log(N)))，2) 傅立叶系数乘以其光谱波数 (O(N))，以及 2)从光谱空间到网格空间的逆变换（同样为 O(N log(N))）。
  • 将傅立叶系数乘以其光谱波数时必须小心。FFT 算法的每个实现似乎都有自己的光谱模式和归一化参数排序。例如，请参阅Math StackExchange 上此问题的答案，了解有关在 MATLAB 中执行此操作的说明。

4. 您在寻找什么级别的准确度？您是否需要在给定的公差范围内计算导数？

  • 对于许多目的，一阶或二阶有限差分方案可能就足够了。为了获得更高的精度，您可以使用更高阶的泰勒展开，去掉高阶项。
  • 如果您需要在给定的容差内计算导数，您可能需要四处寻找具有所需误差的高阶方案。
  • 通常，减少误差的最佳方法是减少有限差分方案中的网格间距，但这并不总是可行的。
  • 请注意，高阶有限差分方案几乎总是需要更大的模板尺寸（更多的相邻点）。这可能会导致边界问题。（参见上面关于鬼点的讨论。）

5. 你的导数在定义函数的相同点上进行评估对你来说重要吗？

  • MATLAB 提供了diff计算相邻数组元素之间差异的函数。这可用于通过一阶前向差分（或前向有限差​​分）方案计算近似导数，但估计是低阶估计。如diff( link ) 的 MATLAB 文档中所述，如果您输入一个长度为 N 的数组，它将返回一个长度为 N-1 的数组。当您使用此方法在 N 个点上估计导数时，您将只能估计 N-1 个点的导数。（请注意，如果它们按升序排序，这可以用于不均匀的网格。）
  • 在大多数情况下，我们希望在所有点都评估导数，这意味着我们希望使用diff方法之外的东西。

6. 需要计算多阶导数吗？

  • 可以建立一个方程组，其中网格点函数值和这些点的一阶和二阶导数都相互依赖。这可以通过像往常一样在相邻点处组合泰勒展开式来找到，但保留导数项而不是取消它们，并将它们与相邻点的导数项联系在一起。这些方程可以通过线性代数求解，不仅给出一阶导数，还给出二阶导数（或更高阶，如果设置正确）。我相信这些被称为组合有限差分方案，它们经常与紧致有限差分方案结合使用，这将在下面讨论。
  • 紧凑的有限差分方案（链接）。在这些方案中，建立一个设计矩阵并通过矩阵求解同时计算所有点的导数。它们被称为“紧凑型”，因为它们通常被设计为比具有可比精度的普通有限差分方案需要更少的模板点。因为它们涉及一个将所有点连接在一起的矩阵方程，所以某些紧凑的有限差分方案据说具有“类似光谱的分辨率”（例如Lele 1992 年的论文——优秀！），这意味着它们通过依赖于所有节点值来模拟光谱方案，因此，它们在所有长度尺度上都保持精度。相比之下，典型的有限差分方法仅是局部准确的（例如，点 #13 的导数通常不依赖于点 #200 的函数值）。
  • 当前的一个研究领域是如何最好地解决紧凑型模板中的多个衍生问题。尽管许多研究人员倾向于根据特定需求（性能、准确性、稳定性或特定研究领域，如流体动力学）来调整它们，但此类研究的结果（组合的、紧凑的有限差分方法）是强大且广泛适用的。

准备就绪的例程

  • 如上所述，可以使用该diff函数（文档链接）来计算相邻数组元素之间的粗略导数。

  • MATLAB 的gradient例程（文档链接）对于许多用途来说都是一个不错的选择。它实现了二阶中心差分方案。它具有计算多维导数和支持任意网格间距的优点。（感谢@thewaywewalk 指出这个明显的遗漏！）

  • 我使用 Fornberg 的方法（见上文）开发了一个小程序（nderiv_fornberg）来计算任意网格间距的一维有限差。我觉得它很容易使用。它在边界处使用 6 点的双面模板，在内部使用居中的 5 点模板。可在此处的 MATLAB 文件交换处获得它。

结论

数值微分领域非常多样化。对于上面列出的每种方法，都有许多变体，它们各有优缺点。这篇文章几乎不是对数值微分的完整处理。

每个应用程序都是不同的。希望这篇文章为感兴趣的读者提供一个有组织的考虑因素和资源列表，以选择适合他们自己需求的方法。

可以使用 MATLAB 特有的代码片段和示例来改进这个社区 wiki。

于 2015-04-06T20:09:10.230 回答
2

我相信这些特定问题还有更多内容。所以我进一步阐述了这个主题如下：

(4) 问：您希望达到什么级别的准确度？您是否需要在给定的公差范围内计算导数？

答：数值微分的准确性取决于感兴趣的应用。通常它的工作方式是，如果您在前向问题中使用 ND 来近似导数以估计感兴趣信号的特征，那么您应该注意噪声扰动。通常这样的伪影包含高频分量，并且根据微分器的定义，噪声效应将以$i\omega^n$ 的数量级被放大。因此，提高微分器的精度（提高多项式精度）根本没有帮助。在这种情况下，您应该能够消除噪声的影响以进行区分。这可以按级联顺序完成：首先平滑信号，然后微分。但更好的方法是使用“低通微分器”。在这里。

但是，如果不是这种情况，并且您在逆问题中使用 ND，例如求解 PDE，那么微分器的全局精度非常重要。根据适合您的问题的边界条件 (BC) 类型，将相应地调整设计。重击规则是增加已知的数值精度，即全频带微分器。您需要设计一个衍生矩阵来处理合适的 BC。您可以使用上述链接找到此类设计的综合解决方案。

(5) 在定义函数的同一点上评估导数对您来说是否重要？ 答：绝对是的。对同一网格点的 ND 评估称为“集中式”和离点“交错”方案。请注意，使用奇数阶导数，集中式 ND 会偏离微分器的频率响应精度。因此，如果您在逆问题中使用这种设计，这将扰乱您的近似。同样，相反的情况适用于交错方案使用的均匀微分顺序的情况。您可以使用上面的链接找到有关此主题的全面说明。

(6) 需要计算多阶导数吗？ 这完全取决于您手头的应用程序。您可以参考我提供的相同链接并处理多个衍生设计。

于 2017-09-05T17:27:16.967 回答