以下是快速且简单的建议,希望能对某些人有所帮助!
1. 您是否拥有一个符号函数或一组点?
- 如果您拥有一个符号函数,则可能可以通过解析计算导数。(很有可能,如果这么容易,您早就已经这么做了,也不会在这里寻找替代方法。)
- 如果您拥有符号函数但无法解析计算导数,则始终可以在一组点上评估函数,并使用此页面上列出的其他方法之一来计算导数。
- 在大多数情况下,您拥有一组点(xi,fi),并且必须使用以下方法之一....
2. 您的网格间距是否均匀?
如果你的网格是均匀分布的,你可能会想要使用有限差分方法(参见维基百科上的
这里或
这里),除非你在使用周期性边界条件(见下文)。
这里是解决普通微分方程时在网格上使用有限差分方法的很好的介绍(特别是第9-14张幻灯片)。这些方法通常计算效率高,实现简单,并且该方法的误差可以简单地估计为用于导出方法的Taylor级数截断误差。
如果你的网格不均匀,你仍然可以使用有限差分方法,但表达式会更加困难,而且精度会随着网格的均匀性变化而变化。如果你的网格非常不均匀,你可能需要使用大的模板大小(更多邻近点)来计算给定点处的导数。人们经常构造一个插值多项式(通常是
Lagrange多项式),并对该多项式求导以计算导数。例如,请参见
此StackExchange问题。使用这些方法通常很难估计误差(虽然一些人已经尝试过:
这里和
这里)。在这些情况下,
Fornberg的方法通常非常有用。
必须注意领域边界,因为模板通常涉及到在领域外的点。一些人引入“幽灵点”或将边界条件与不同阶数的导数相结合,以消除这些“幽灵点”并简化模板。另一种方法是使用右侧或左侧的有限差分方法。
这里是一个极好的有限差分方法“备忘单”,包括低阶的中心、右侧和左侧方案。我将它的打印件放在我的工作站附近,因为我发现它非常有用。
3. 您的域是否是周期性的?您可以假定周期性边界条件吗?
- 如果您的域是周期性的,您可以使用傅里叶谱方法计算高阶精度导数。这种技术为获得高精度而牺牲了一定的性能。事实上,如果您使用N个点,则导数的估计值约为N的阶数精确。有关更多信息,请参见(例如)this WikiBook。
- 傅里叶方法通常使用快速傅里叶变换(FFT)算法,以实现大约O(N log(N))的性能,而不是naively-implemented离散傅里叶变换(DFT)可能采用的O(N ^ 2)算法。
- 如果您的函数和域不周期性,则不应使用傅里叶谱方法。如果您尝试在非周期函数中使用它,则会产生大误差和不良的“环形”现象。
- 计算任意阶导数需要1)从网格空间到谱空间的变换(O(N log(N))),2)将傅里叶系数乘以它们的谱波数(O(N)),以及2)从谱空间到网格空间的逆变换(再次是O(N log(N)))。
- 在将傅里叶系数乘以它们的谱波数时,必须小心。FFT算法的每个实现似乎都有其自己的谱模式排序和归一化参数。例如,请参见Math StackExchange上this question的答案,了解在MATLAB中执行此操作的注意事项。
4. 您需要什么精度水平?您需要在给定的容差范围内计算导数吗?
- 对于许多目的,第一或第二阶有限差分方案可能足够。为了更高的精度,您可以使用更高阶的泰勒展开,省略高阶项。
- 如果您需要在给定的容差范围内计算导数,则可能需要寻找具有所需误差的高阶方案。
- 通常,减小有限差分方案中的网格间距是减小误差的最佳方法,但这并非总是可行。
- 请注意,更高阶的有限差分方案几乎总是需要更大的模板大小(更多相邻点)。这可能会在边界处引起问题。(请参见上面关于幽灵点的讨论。)
5. 您是否关心您的导数在与函数定义相同的点上评估?
MATLAB提供了
diff
函数来计算相邻数组元素之间的差异。这可以用于通过一阶前向差分(或前向有限差分)方案计算近似导数,但是估计是低阶估计。如MATLAB对
diff
的文档所述(
link),如果您输入长度为N的数组,则它将返回长度为N-1的数组。当您在N个点上使用此方法估计导数时,您只会在N-1个点处得到导数的估计值。(请注意,如果它们按升序排序,则可以在不均匀的网格上使用此方法。)
在大多数情况下,我们希望在所有点处评估导数,这意味着我们想使用除
diff
方法之外的东西。
6. 您需要计算多个导数阶数吗?
- 可以建立一个方程系统,其中网格点函数值和这些点的一阶和二阶导数都彼此依赖。这可以通过将相邻点的泰勒展开组合在一起来找到,但保留导数项而不是抵消它们,并将它们与相邻点的导数链接在一起。这些方程可以通过线性代数求解,不仅可以得到一阶导数,还可以得到二阶导数(如果正确设置了高阶导数也可以得到)。我认为这些被称为组合有限差分方案,通常与紧致有限差分方案一起使用,下面将讨论它们。
- 紧致有限差分方案(link)。在这些方案中,人们设置一个设计矩阵,并通过矩阵求解同时计算所有点的导数。它们之所以被称为"紧致",是因为它们通常被设计为需要比具有相同精度的普通有限差分方案更少的模板点。由于它们涉及将所有点链接在一起的矩阵方程,某些紧致有限差分方案被称为具有"类似谱解析"的分辨率(例如Lele's 1992 paper--非常好的!),这意味着它们通过依赖所有节点值来模仿谱方法,因此在所有长度尺度上保持准确性。相比之下,典型的有限差分方法仅在局部精确(例如,点#13处的导数通常不取决于点#200处的函数值)。
- 当前的研究领域是如何最好地解决紧致模板中的多个导数。这种研究的结果,组合,紧致有限差分方法,是功能强大且广泛适用的,尽管许多研究人员倾向于为特定需求进行调整(性能、精度、稳定性或特定研究领域,如流体动力学)。
准备就绪的例程
- 如上所述,可以使用
diff
函数(链接到文档)计算相邻数组元素之间的粗略导数。
MATLAB的gradient
例程(链接到文档)是许多目的的绝佳选择。它实现了一个二阶中心差分方案。它具有在多个维度中计算导数和支持任意网格间距的优点。(感谢@thewaywewalk指出这个明显的遗漏!)
我使用Fornberg的方法(见上文)开发了一个小例程(nderiv_fornberg
),用于计算任意网格间距下一维有限差分。我发现它很容易使用。它在边界处使用6个点的侧向模板,在内部使用中心化的5点模板。它可在MATLAB文件交换中心 这里获得。
结论
数值微分领域十分广泛。对于上述每种方法,都有许多变体具有自己的优缺点。这篇文章并不是关于数值微分的完整介绍。
每个应用都是不同的。希望这篇文章为感兴趣的读者提供了一个组织良好的考虑因素和资源列表,以选择适合自己需求的方法。
这个社区维基可以通过特定于MATLAB的代码片段和示例进行改进。