计算匿名函数的二阶导数答案

【问题标题】：Calculate second derivative of anonymous function计算匿名函数的二阶导数
【发布时间】：2019-04-10 15:01:01
【问题描述】：

我想在 Matlab 中计算一个匿名函数的二阶导数。我已经知道一些公式（数值微分），但它们似乎不起作用。

我可以计算一阶导数：

f = @(x) (x^3);
h = 1e-10;

df = @(x) (f(x+h) - f(x))/h;

但是当我尝试用以下公式计算二阶导数时，我没有得到预期的结果：

f = @(x) (x^3);
h = 1e-10;

d2f = @(x) (f(x+h) - 2*f(x) + f(x-h))/(h^2);

对于 d2f，我应该得到一个类似于 d2f = 6x 的函数，但是如果绘制 d2f，我会得到： plot d2f

我做错了什么？

【问题讨论】：

如果您提供有关哪些功能不起作用的信息会有所帮助。你在尝试什么，你期望什么结果，你会得到什么？
另见：stackoverflow.com/questions/22160712/…
这里的问题在于h 和catastrophic cancellation 的大小，这是由于浮点的数值精度。试试h=1e-7，或更大。但是，花在评估函数的位置上，您可能仍然会看到很多噪音。
您也可以使用 vpa:h = vpa(1e-10) 增加h 的浮点精度。您现在有 32 位精度，因此 h 应该能够与 1e-16 一样小。
@horchler 成功了！谢谢！！

标签： matlab numerical-methods differentiation

【解决方案1】：

除法差公式的理论误差为O(h^2)。函数的浮点计算将各自产生一个机器精度μ左右的相对误差。然后除以 h^2。两个误差的最佳总和是在它们达到平衡的地方，即 h^4=mu 或 h=1e-4 的地方。

如果误差项的系数（即 f 的 4 次导数）为零，这当然是无效的，因为 f(x)=x^3 会发生这种情况。那么唯一的误差贡献是浮点误差，对于较大的 h，浮点误差最小，即使 h=1 也会产生最小的误差。

对于像 f(x)=sin(x) 这样不那么简单的函数，不同 h 的误差如下图所示（其中标记为 x 的变量是步长 h）

【讨论】：

【解决方案2】：

我不确定你做错了什么，但下面的代码有效

f=@(x) x.^3;

x = (0:1E-12:1E-6)' ;
d2y = secondDerivative(f,x(1),x(end),x(2)-x(1))';

fit(x,d2y,'poly1')

ans = 

     Linear model Poly1:
     ans(x) = p1*x + p2
     Coefficients (with 95% confidence bounds):
       p1 =           6  (6, 6)
       p2 =   1.352e-15  (-4.734e-13, 4.761e-13)

函数定义

function d2y = secondDerivative(f, x1, x2, dx)

y = f(x1:dx:x2);

d2y = nan(size(y));
d2y(2:end-1) = y(1:end-2) - 2*y(2:end-1) + y(3:end);

if length(d2y) == 3
    d2y(1) = y(1) - 2*y(2) + y(3);
    d2y(2) = y(end-2) - 2*y(end-1) + y(end);
elseif length(d2y) > 4
    d2y(1) = 2*y(1) - 5*y(2) + 4*y(3) - y(4);
    d2y(end) = -y(end-3) + 4*y(end-2) - 5*y(end-1) + 2*y(end);
end

d2y = d2y / dx^2 ;
end

【讨论】：