计算匿名函数的二阶导数

Question

我想在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

我做错了什么？

Answer 1

分割差分公式的理论误差为O（h ^ 2）。函数的浮点评估将产生关于机器精度mu的相对误差。然后将其除以h ^ 2。两个误差的最佳总和是在它们关于平衡的情况下实现的，即，其中h ^ 4 = mu或h = 1e-4。

如果误差项的系数（f是f的四阶导数）为零，这当然是无效的，因为它发生在f（x）= x ^ 3处。那么唯一的误差贡献是浮点误差，对于较大的h，它是最小的，即使h = 1也会给出最小的误差。

对于像f（x）= sin（x）这样不太重要的函数，不同h的误差表现如下图（其中标记x的变量是步长h）

Answer 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