在MATLAB中计算数值导数的最佳方法是什么? [英] What's the best way to calculate a numerical derivative in MATLAB?

问题描述

(注意:这旨在成为社区Wiki.)

(Note: This is intended to be a community Wiki.)

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

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. (In some situations, we might know f(x) ahead of time, but we want to do this generally, since we often don't know f(x) in advance.) What's a good way to approximate the derivative of f(x) at each point xi? That is, how can I estimate values of dfi == d/dx fi == df(xi)/dx at each of the points xi?

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

Unfortunately, MATLAB doesn't have a very good general-purpose, numerical differentiation routine. Part of the reason for this is probably because choosing a good routine can be difficult!

那么有哪些方法?存在哪些例程?我们如何为特定问题选择良好的例行程序?

So what kinds of methods are there? What routines exist? How can we choose a good routine for a particular problem?

选择如何在MATLAB中进行区分时需要考虑以下几点:

1. 您是否具有符号功能或一组要点?
2. 您的网格是均匀分布还是不均匀分布?
3. 您的域是否定期?您可以假设周期性边界条件吗?
4. 您正在寻找什么水平的准确性?您是否需要在给定的公差范围内计算导数?
5. 对您的导数在定义函数的相同点上求值对您来说是否重要?
6. 您需要计算多个衍生产品订单吗?

最好的进行方式是什么?

What's the best way to proceed?

推荐答案

这些只是一些简单的建议.希望有人能对他们有所帮助！

These are just some quick-and-dirty suggestions. Hopefully somebody will find them helpful!

1.您具有符号功能或一组点吗?

• 如果您具有符号函数，则 可能能够解析地计算导数. (有可能，如果那么简单，您可能会这样做，并且您不会在这里寻找替代方案.)
• 如果您具有符号函数并且无法解析计算导数，则可以始终对一组点求值函数，并使用此页面上列出的其他一些方法来求导数.
• 在大多数情况下，您具有一组点(xi，fi)，并且必须使用以下方法之一....

• If you have a symbolic function, you may be able to calculate the derivative analytically. (Chances are, you would have done this if it were that easy, and you would not be here looking for alternatives.)
• If you have a symbolic function and cannot calculate the derivative analytically, you can always evaluate the function on a set of points, and use some other method listed on this page to evaluate the derivative.
• In most cases, you have a set of points (xi,fi), and will have to use one of the following methods....

2.您的网格是均匀分布还是不均匀分布?

• 如果网格间隔 均匀，则可能需要使用有限差分方案(请参阅Wikipedia文章此处)，除非您是使用周期性边界条件(请参见下文). 此处是对有限项的不错介绍在求解网格上的常微分方程的上下文中的差分方法(尤其参见幻灯片9-14).这些方法通常计算效率高，易于实现，并且该方法的误差可以简单地估计为用于推导该方法的泰勒展开式的截断误差.
• 如果网格间隔 不均匀，您仍然可以使用有限差分方案，但是表达式更加困难，并且精度随网格的均匀性变化很大.如果网格非常不均匀，则可能需要使用较大的模具尺寸(更多相邻点)来计算给定点的导数.人们通常构造一个插值多项式(通常是 Lagrange多项式)，然后对该多项式求微分以计算导数.例如，参见这个 StackExchange问​​题.使用这些方法通常很难估计错误(尽管有些方法已尝试这样做:此处和此处). Fornberg方法在这些情况下通常非常有用....
• 必须在域的边界处进行护理，因为模具通常涉及域外的点.有人引入重影点"或将边界条件与不同阶数的导数组合以消除这些重影点"并简化模板.另一种方法是使用右侧或左侧的有限差分方法.
• 这是一个优秀作弊者"有限差分法的表格"，包括低阶的居中，右侧和左侧方案.我将其打印输出保存在我的工作站附近，因为我发现它非常有用.

• If your grid is evenly spaced, you probably will want to use a finite difference scheme (see either of the Wikipedia articles here or here), unless you are using periodic boundary conditions (see below). Here is a decent introduction to finite difference methods in the context of solving ordinary differential equations on a grid (see especially slides 9-14). These methods are generally computationally efficient, simple to implement, and the error of the method can be simply estimated as the truncation error of the Taylor expansions used to derive it.
• If your grid is unevenly spaced, you can still use a finite difference scheme, but the expressions are more difficult and the accuracy varies very strongly with how uniform your grid is. If your grid is very non-uniform, you will probably need to use large stencil sizes (more neighboring points) to calculate the derivative at a given point. People often construct an interpolating polynomial (often the Lagrange polynomial) and differentiate that polynomial to compute the derivative. See for instance, this StackExchange question. It is often difficult to estimate the error using these methods (although some have attempted to do so: here and here). Fornberg's method is often very useful in these cases....
• Care must be taken at the boundaries of your domain because the stencil often involves points that are outside the domain. Some people introduce "ghost points" or combine boundary conditions with derivatives of different orders to eliminate these "ghost points" and simplify the stencil. Another approach is to use right- or left-sided finite difference methods.
• Here's an excellent "cheat sheet" of finite difference methods, including centered, right- and left-sided schemes of low orders. I keep a printout of this near my workstation because I find it so useful.

3.您的域是定期的吗?您可以假设周期性边界条件吗?

• 如果您的域是周期性的，则可以使用傅立叶频谱方法以很高的精度计算导数.该技术为了获得高精度而在某种程度上牺牲了性能.实际上，如果您使用的是N个点，则您对导数的估计大约是N ^阶准确.有关详细信息，请参阅(例如)此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中执行此操作的说明.

• If your domain is periodic, you can compute derivatives to a very high order accuracy using Fourier spectral methods. This technique sacrifices performance somewhat to gain high accuracy. In fact, if you are using N points, your estimate of the derivative is approximately N^th order accurate. For more information, see (for example) this WikiBook.
• Fourier methods often use the Fast Fourier Transform (FFT) algorithm to achieve roughly O(N log(N)) performance, rather than the O(N^2) algorithm that a naively-implemented discrete Fourier transform (DFT) might employ.
• If your function and domain are not periodic, you should not use the Fourier spectral method. If you attempt to use it with a function that is not periodic, you will get large errors and undesirable "ringing" phenomena.
• Computing derivatives of any order requires 1) a transform from grid-space to spectral space (O(N log(N))), 2) multiplication of the Fourier coefficients by their spectral wavenumbers (O(N)), and 2) an inverse transform from spectral space to grid space (again O(N log(N))).
• Care must be taken when multiplying the Fourier coefficients by their spectral wavenumbers. Every implementation of the FFT algorithm seems to have its own ordering of the spectral modes and normalization parameters. See, for instance, the answer to this question on the Math StackExchange, for notes about doing this in MATLAB.

4.您正在寻找什么水平的精度?您是否需要在给定的公差范围内计算导数?

• 出于许多目的，一阶或二阶有限差分方案可能就足够了.为了获得更高的精度，您可以使用更高阶的泰勒展开式，从而舍弃更高阶的项.
• 如果需要在给定的公差范围内计算导数，则可能需要查看具有所需误差的高阶方案.
• 减少误差的最佳方法通常是在有限差分方案中减小网格间距，但这并不总是可能的.
• 请注意，高阶有限差分方案几乎总是需要较大的模具尺寸(更多的相邻点).这可能会在边界处引起问题. (请参见上面有关幻影点的讨论.)

5.对您的导数在定义函数的相同点上进行评估对您来说是否重要?

• MATLAB提供了diff函数来计算相邻数组元素之间的差异.这可以用于通过一阶前向微分(或前向有限差​​分)方案来计算近似导数，但估计值是低阶估计值.如diff的MATLAB文档中所述(链接 )，如果您输入长度为N的数组，它将返回长度为N-1的数组.当使用此方法在N点上估计导数时，您将仅在N-1点上得到导数的估计. (请注意，如果它们以升序排序，则可以在不均匀的网格上使用.)
• 在大多数情况下，我们希望在所有点上对导数进行求值，这意味着我们希望使用diff方法之外的方法.

• MATLAB provides the diff function to compute differences between adjacent array elements. This can be used to calculate approximate derivatives via a first-order forward-differencing (or forward finite difference) scheme, but the estimates are low-order estimates. As described in MATLAB's documentation of diff (link), if you input an array of length N, it will return an array of length N-1. When you estimate derivatives using this method on N points, you will only have estimates of the derivative at N-1 points. (Note that this can be used on uneven grids, if they are sorted in ascending order.)
• In most cases, we want the derivative evaluated at all points, which means we want to use something besides the diff method.

6.您需要计算多个衍生产品订单吗?

可以建立一个方程组，在该方程组中，网格点函数值以及这些点的一阶和二阶导数都相互依赖.这可以通过像往常一样在相邻点处合并泰勒展开式，但保留导数项而不是将其抵消掉并将它们与相邻点的项链接在一起来找到.这些方程式可以通过线性代数求解，不仅给出一阶导数，而且给出二阶导数(如果设置正确，也可以给出更高阶的阶数).我相信这些被称为 combined 有限差分方案，并且它们通常与 compact 有限差分方案结合使用，这将在下面进行讨论.
• 紧凑的有限差分方案(链接).在这些方案中，建立一个设计矩阵，并通过矩阵求解同时计算所有点的导数.它们之所以被称为紧凑型"，是因为它们通常被设计为比可比较精度的普通有限差分方案需要更少的模板点.因为它们包含将所有点链接在一起的矩阵方程，所以某些紧凑的有限差分方案被称为具有类似光谱的分辨率"(例如

• One can set up a system of equations in which the grid point function values and the 1st and 2nd order derivatives at these points all depend on each other. This can be found by combining Taylor expansions at neighboring points as usual, but keeping the derivative terms rather than cancelling them out, and linking them together with those of neighboring points. These equations can be solved via linear algebra to give not just the first derivative, but the second as well (or higher orders, if set up properly). I believe these are called combined finite difference schemes, and they are often used in conjunction with compact finite difference schemes, which will be discussed next.
• Compact finite difference schemes (link). In these schemes, one sets up a design matrix and calculates the derivatives at all points simultaneously via a matrix solve. They are called "compact" because they are usually designed to require fewer stencil points than ordinary finite difference schemes of comparable accuracy. Because they involve a matrix equation that links all points together, certain compact finite difference schemes are said to have "spectral-like resolution" (e.g. Lele's 1992 paper--excellent!), meaning that they mimic spectral schemes by depending on all nodal values and, because of this, they maintain accuracy at all length scales. In contrast, typical finite difference methods are only locally accurate (the derivative at point #13, for example, ordinarily doesn't depend on the function value at point #200).
• A current area of research is how best to solve for multiple derivatives in a compact stencil. The results of such research, combined, compact finite difference methods, are powerful and widely applicable, though many researchers tend to tune them for particular needs (performance, accuracy, stability, or a particular field of research such as fluid dynamics).

即用型例程

• 如上所述，可以使用diff函数(链接到文档)以计算相邻数组元素之间的粗导数.

• MATLAB的gradient例程(链接到文档)是多种用途的绝佳选择.它实现了二阶中央差分方案.它具有在多个维度上计算导数并支持任意网格间距的优点. (感谢@thewaywewalk指出了这一明显的遗漏！)

• As described above, one can use the diff function (link to documentation) to compute rough derivatives between adjacent array elements.

• MATLAB's gradient routine (link to documentation) is a great option for many purposes. It implements a second-order, central difference scheme. It has the advantages of computing derivatives in multiple dimensions and supporting arbitrary grid spacing. (Thanks to @thewaywewalk for pointing out this glaring omission!)

我使用Fornberg方法(参见上文)开发了一个小例程(nderiv_fornberg)，以计算任意网格间距的一维有限差分.我发现它易于使用.它在边界处使用6个点的侧模板，在内部使用一个居中的5点模板.可在MATLAB文件交换中找到此处.

I used Fornberg's method (see above) to develop a small routine (nderiv_fornberg) to calculate finite differences in one dimension for arbitrary grid spacings. I find it easy to use. It uses sided stencils of 6 points at the boundaries and a centered, 5-point stencil in the interior. It is available at the MATLAB File Exchange here.

结论

数值微分的领域非常多样化.对于上面列出的每种方法，都有许多变体，各有其优缺点.这篇文章很难完全解决数值微分问题.

The field of numerical differentiation is very diverse. For each method listed above, there are many variants with their own set of advantages and disadvantages. This post is hardly a complete treatment of numerical differentiation.

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

Every application is different. Hopefully this post gives the interested reader an organized list of considerations and resources for choosing a method that suits their own needs.

可以通过MATLAB特有的代码片段和示例来改进此社区Wiki.

This community wiki could be improved with code snippets and examples particular to MATLAB.

这篇关于在MATLAB中计算数值导数的最佳方法是什么?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！