首页 文章

在MATLAB中计算数值导数的最佳方法是什么？

提问于
浏览 次
16

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

Suppose I have a set of points xi = {x0,x1,x2,...xn} and corresponding function values fi = f(xi) = {f0,f1,f2,...,fn}, where f(x) is, in general, an unknown function. （在某些情况下，我们可能提前知道f（x），但我们想要这样做，因为我们通常不提前知道f（x） . ） What's a good way to approximate the derivative of f(x) at each point xi? 也就是说，我怎样才能估算出dfi == d / dx fi == df（xi）/ dx在每个点xi？

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

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

There are several considerations when choosing how to differentiate in MATLAB:

  • 你有一个象征性的功能还是一组点？

  • 你的网格是均匀的还是不均匀的？

  • 您的域名是定期的吗？你能假设周期性边界条件吗？

  • 您在寻找什么样的准确度？您是否需要在给定的容差范围内计算导数？

  • 您的衍 生产环境 品是否与您定义的函数在相同的点上进行评估，这对您有用吗？

  • 您需要计算多个衍生品订单吗？

什么是最好的方法？

matlab numerical-methods derivative differentiation

2 回答

  • 18

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

    (4) Q: What level of accuracy are you looking for? Do you need to compute the derivatives within a given tolerance?

    答：数值微分的准确性对于利益的应用是主观的 . 通常它的工作方式是，如果你在前向问题中使用ND近似导数来估计感兴趣的信号的特征，那么你应该知道噪声扰动 . 通常这种伪像包含高频分量，并且通过微分器的定义，噪声效应将以$ i \ omega ^ n $的大小顺序放大 . 因此，提高微分器的精度（提高多项式精度）根本没有用 . 在这种情况下，您应该能够取消噪声对差异化的影响 . 这可以在casecade顺序中完成：首先平滑信号，然后进行区分 . 但更好的方法是使用"Lowpass Differentiator" . 可以找到MATLAB库的一个很好的例子here .

    但是，如果情况并非如此，并且您在反问题中使用ND，例如solvign PDE，那么微分器的全局精度非常重要 . 根据适合您的问题的适合条件（BC），设计将相应调整 . 砰砰声的规则是增加数值精度，已知是全频带微分器 . 您需要设计一个处理合适BC的衍生矩阵 . 您可以使用上述链接找到此类设计的综合解决方案 .

    (5) Does it matter to you that your derivative is evaluated on the same points as your function is defined? 答：是的 . 在相同网格点上对ND的评估称为"centralized"，并且在点"staggered"方案之外 . 注意，使用奇数阶导数，集中式ND将偏离微分器的频率响应精度 . 因此，如果您在逆问题中使用此类设计，这将扰乱您的近似值 . 而且，相反的情况适用于交错方案所使用的均匀分化顺序的情况 . 您可以使用上面的链接找到关于此主题的全面解释 .

    (6) Do you need to calculate multiple orders of derivatives? 这完全取决于您的应用程序 . 您可以参考我提供的相同链接并处理多种衍生设计 .

    回复于
  • 1

    这些只是一些快速而肮脏的建议 . 希望有人会发现它们有用！

    1. Do you have a symbolic function or a set of points?

    • 如果您有符号功能，您可以分析计算衍生物 . （很可能，如果那么容易的话，你会做到这一点，你不会在这里寻找替代品 . ）

    • 如果您具有符号函数且无法通过分析计算导数，则可以始终在一组点上评估函数，并使用此页面上列出的其他方法来评估导数 .

    • 在大多数情况下，你有一组点（xi，fi），并且必须使用以下方法之一....

    2. Is your grid evenly or unevenly spaced?

    • 如果网格间隔均匀，您可能需要使用有限差分格式（请参阅维基百科文章here或here），除非您使用的是周期性边界条件（见下文） . Here是在解决网格上的常微分方程的背景下对有限差分方法的一个不错的介绍（特别参见幻灯片9-14） . 这些方法通常在计算上有效，易于实现，并且该方法的误差可以简单地估计为用于导出它的泰勒展开的截断误差 .

    • 如果您的网格间距不均匀，您仍然可以使用有限差分格式，但表达式更加困难，并且精度随着网格的统一程度而变化很大 . 如果您的网格非常不均匀，您可能需要使用较大的模板尺寸（更多相邻点）来计算给定的导数点 . 人们经常构造插值多项式（通常是Lagrange polynomial）并区分该多项式来计算导数 . 例如，参见this StackExchange问题 . 使用这些方法通常很难估计错误（尽管有些人试图这样做：here和here） . Fornberg's method在这些情况下通常非常有用....

    • 必须注意您域名的边界，因为模板通常涉及域外的点 . 有些人引入"ghost points"或将边界条件与不同阶数的导数组合以消除这些并简化模板 . 另一种方法是使用右侧或左侧有限差分方法 .

    • Here's一种优秀的有限差分方法，包括低阶的中心，右侧和左侧方案 . 我在工作站附近打印出这个，因为我发现它非常有用 .

    3. Is your domain periodic? Can you assume periodic boundary conditions?

    • 如果您的域是周期性的，您可以使用傅里叶谱方法计算导数到非常高的阶数 . 这种技术在某种程度上牺牲了性能以获得高精度 . 实际上，如果您使用的是N点，那么您对导数的估计大约是第N阶的准确度 . 有关更多信息，请参阅（例如）this WikiBook .
      傅立叶方法通常使用快速傅立叶变换（FFT）算法来实现大致O（N log（N））性能，而不是天然实现的离散傅立叶变换（DFT）可能采用的O（N ^ 2）算法 .

    • 如果您的函数和域不是周期性的，则不应使用傅里叶谱方法 . 如果您尝试将其与非周期性函数一起使用，则会出现大的错误和不良的"ringing"现象 .

    • 计算任何阶数的导数需要1）从网格空间到光谱空间的变换（O（N log（N））），2）傅里叶系数乘以其光谱波数（O（N）），和2）从光谱空间到网格空间的逆变换（再次为O（N log（N））） .

    • 在将傅里叶系数乘以其光谱波数时必须小心 . FFT算法的每个实现似乎都有自己的频谱模式和归一化参数的排序 . 例如，请参阅Math StackExchange上this question的答案，以获取有关在MATLAB中执行此操作的说明 .

    4. What level of accuracy are you looking for? Do you need to compute the derivatives within a given tolerance?

    • 出于许多目的，一阶或二阶有限差分方案可能就足够了 . 为了获得更高的精度，您可以使用更高阶的泰勒展开式，从而降低高阶项 .

    • 如果需要计算给定容差范围内的导数，可能需要查看具有所需错误的高阶方案 .

    • 通常，减少误差的最佳方法是在有限差分方案中减少网格间距，但这并非总是可行的 .

    • 请注意，高阶有限差分格式几乎总是需要更大的模板尺寸（更多的相邻点） . 这可能会导致边界问题 . （参见上面关于鬼点的讨论 . ）

    5. Does it matter to you that your derivative is evaluated on the same points as your function is defined?

    • MATLAB提供 diff 函数来计算相邻数组元素之间的差异 . 这可以用于通过一阶前向差分（或前向有限差分）方案来计算近似导数，但估计是低阶估计 . 如MATLAB的 diff （link）文档中所述，如果输入长度为N的数组，它将返回一个长度为N-1的数组 . 当您在N点上使用此方法估计导数时，您只能估算出N-1点的导数 . （请注意，如果它们按升序排序，则可以在不均匀网格上使用 . ）

    • 在大多数情况下，我们希望在所有点评估衍生物，这意味着我们想要使用除 diff 方法之外的其他东西 .

    6. Do you need to calculate multiple orders of derivatives?

    • 可以 Build 一个方程组，其中网格点函数值和这些点上的一阶和二阶导数都相互依赖 . 这可以通过像往常一样在相邻点处组合泰勒展开来找到，但是保持导数项而不是将它们抵消，并将它们与相邻点的那些连接在一起 . 这些方程式可以通过线性代数求解，不仅可以得到一阶导数，也可以得到第二个导数（或者更高阶，如果设置得当） . 我相信这些被称为组合有限差分方案，它们通常与紧致有限差分方案结合使用，这将在下面讨论 .

    • 紧凑有限差分格式（link） . 在这些方案中， Build 一个设计矩阵并通过矩阵求解同时计算所有点的导数 . 它们被称为"compact"因为它们通常被设计为需要比具有可比精度的普通有限差分方案更少的模板点 . 因为它们涉及将所有点连接在一起的矩阵方程，确定据说紧凑有限差分方案具有"spectral-like resolution"（例如Lele's 1992 paper - 优秀！），这意味着它们通过依赖于所有节点值来模拟频谱方案，并且因此，它们在所有长度尺度上保持准确性 . 相反，典型的有限差分方法仅局部精确（例如，点＃13处的导数通常不依赖于点＃200处的函数值） .

    • 目前的研究领域是如何在紧凑的模板中最好地解决多种衍生物 . 这种研究，组合，紧凑的有限差分方法的结果是强大的和广泛适用的，尽管许多研究人员倾向于针对特定需求（性能，准确性，稳定性或诸如流体动力学的特定研究领域）调整它们 .

    Ready-to-Go Routines

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

    • MATLAB的 gradient 例程（link到文档）是一个很好的选择，用于许多目的 . 它实现了二阶中心差分方案 . 它具有计算多维导数并支持任意网格间距的优点 . （感谢@thewaywewalk指出这个明显的遗漏！）

    • 我用Fornberg的方法（见上文）开发了一个小程序（ nderiv_fornberg ）来计算任意网格间距的一维的有限差分 . 我觉得它很容易使用 . 它在边界处使用6个点的侧面模板，在内部使用中心的5点模板 . 它可以在MATLAB File Exchange here中找到 .

    Conclusion

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

    每个应用程序都不同 . 希望这篇文章能够为感兴趣的读者提供有组织的考虑因素和资源清单，以便选择适合自己需要的方法 .

    可以使用特定于MATLAB的代码片段和示例来改进此社区wiki .

    回复于

相关问题