math - 如何评估matlab中函数的导数？

Question

这应该很简单。我有一个函数 f(x)，我想在 MATLAB 中为给定的 x 计算 f'(x)。

我所有的搜索都得出了符号数学，这不是我需要的，我需要数值微分。

例如如果我定义:fx = inline('x.^2')

我想找到 say f'(3)，即 6，我不想找到 2x

Answer 1

如果已知您的函数是二次可微的，请使用

f'(x) = (f(x + h) - f(x - h)) / 2h

这是 h 中的二阶精度。如果只有一次可微，使用

f'(x) = (f(x + h) - f(x)) / h     (*)

这是 h 中的一阶。

这是理论。实际上，事情非常棘手。我将采用第二个公式(一阶)，因为分析更简单。做第二个订单作为练习。

首先要注意的是，您必须确保 (x + h) - x = h，否则会出现巨大错误。事实上，f(x + h) 和 f(x) 彼此接近(比如 2.0456 和 2.0467)，当你减去它们时，你会失去很多有效数字(这里是 0.0011，少了 3 个有效数字比 x)。因此，h 上的任何错误都可能对结果产生巨大影响。

因此，第一步，确定一个候选 h(稍后我将向您展示如何选择它)，然后将 h' = (x + h) - x 作为 h 用于计算。如果您使用的是 C 之类的语言，则必须注意将 h 或 x 定义为可变的，以免该计算被优化掉。

接下来是h的选择。 (*) 中的误差有两部分:截断误差和舍入误差。截断错误是因为公式不准确:

(f(x + h) - f(x)) / h = f'(x) + e1(h)

其中 e1(h) = h/2 * sup_{x in [0,h]} |f''(x)|。

舍入误差来自 f(x + h) 和 f(x) 彼此接近的事实。可以粗略估计为

e2(h) ~ epsilon_f |f(x) / h|

其中 epsilon_f 是计算 f(x)(或接近的 f(x + h))的相对精度。这必须根据您的问题进行评估。对于简单的函数，epsilon_f 可以作为机器 epsilon。对于更复杂的，它可能比它差几个数量级。

所以你需要 h 来最小化 e1(h) + e2(h)。将所有内容整合在一起并优化 h 产量

h ~ sqrt(2 * epsilon_f * f / f'')

必须根据您的函数进行估算。你可以粗略估计一下。如有疑问，请取 h ~ sqrt(epsilon)，其中 epsilon = 机器精度。对于 h 的最佳选择，已知导数的相对精度为 sqrt(epsilon_f)，即。一半的有效数字是正确的。

简而言之:a h 太小 => 舍入误差，a h 太大 => 截断误差。

对于二阶公式，同样的计算结果

h ~ (6 * epsilon_f / f''')^(1/3)

导数的分数精度为 (epsilon_f)^(2/3)(假设 double ，通常比一阶公式好一个或两个有效数字)。

如果这太不精确，请随时询问更多方法，有很多技巧可以获得更好的准确性。理查森外推法是平滑函数的良好开端。但是这些方法通常会计算 f 很多次，如果您的函数很复杂，这可能是您想要的，也可能不是您想要的。

如果您打算在不同点多次使用数值导数，那么构建切比雪夫近似值会变得很有趣。