如何表示二阶数值微分的截断误差是一阶数值微分的平方

Question

我试图证明二阶数值微分的截断误差确实比一阶数值微分的精度高出一倍(考虑了机器误差/舍入误差eps())

下面是我在Julia中的代码：

function first_order_numerical_D(f)
  function df(x)
    h = sqrt(eps(x))        
    (f(x+h) - f(x))/h
  end
  df
end

function second_order_numerical_D(f)
  function df(x)
    h = sqrt(eps(x))          
    (f(x+h) - f(x-h))/(2.0*h)
  end
  df
end

function analytical_diff_exp(x)
    return exp(x)
end
function analytical_diff_sin(x)
    return cos(x)
end
function analytical_diff_cos(x)
    return -sin(x)
end
function analytical_diff_sqrt(x)
    return 1/(2.0*sqrt(x))
end

function first_order_error_exp(x)
    return first_order_numerical_D(exp)(x) - analytical_diff_exp(x) 
end
function first_order_error_sin(x)
    return first_order_numerical_D(sin)(x) - analytical_diff_sin(x) 
end
function first_order_error_cos(x)
    return first_order_numerical_D(cos)(x) - analytical_diff_cos(x) 
end
function first_order_error_sqrt(x)
    return first_order_numerical_D(sqrt)(x) - analytical_diff_sqrt(x) 
end

function second_order_error_exp(x)
    return second_order_numerical_D(exp)(x) - analytical_diff_exp(x)
end
function second_order_error_sin(x)
    return second_order_numerical_D(sin)(x) - analytical_diff_sin(x)
end
function second_order_error_cos(x)
    return second_order_numerical_D(cos)(x) - analytical_diff_cos(x) 
end
function second_order_error_sqrt(x)
    return second_order_numerical_D(sqrt)(x) - analytical_diff_sqrt(x) 
end

function round_off_err_exp(x)
    return 2.0*sqrt(eps(x))*exp(x)
end
function round_off_err_sin(x)
    return 2.0*sqrt(eps(x))*sin(x)
end
function round_off_err_cos(x)
    return 2.0*sqrt(eps(x))*cos(x)
end
function round_off_err_sqrt(x)
    return 2.0*sqrt(eps(x))*sqrt(x)
end

function first_order_truncation_err_exp(x)
    return abs(first_order_error_exp(x)+round_off_err_exp(x))
end
function first_order_truncation_err_sin(x)
    return abs(first_order_error_sin(x)+round_off_err_sin(x))
end
function first_order_truncation_err_cos(x)
    return abs(first_order_error_cos(x)+round_off_err_cos(x))
end
function first_order_truncation_err_sqrt(x)
    return abs(first_order_error_sqrt(x)+round_off_err_sqrt(x))
end

function second_order_truncation_err_exp(x)
    return abs(second_order_error_exp(x)+0.5*round_off_err_exp(x))
end
function second_order_truncation_err_sin(x)
    return abs(second_order_error_sin(x)+0.5*round_off_err_sin(x))
end
function second_order_truncation_err_cos(x)
    return abs(second_order_error_cos(x)+0.5*round_off_err_cos(x))
end
function second_order_truncation_err_sqrt(x)
    return abs(second_order_error_sqrt(x)+0.5*round_off_err_sqrt(x))
end

如果我减去round_off_err_f项(这里我使用add，因为实际的泰勒展开显示舍入误差和截断误差前面都有一个负号)，这应该会给出正确的截断误差。

有关分析推导/证明的信息，请参阅：https://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h10/kompendiet/kap11.pdf http://www2.math.umd.edu/~dlevy/classes/amsc466/lecture-notes/differentiation-chap.pdf

但研究结果表明：

first_order_truncation_err_exp(0.5), first_order_truncation_err_sin(0.5), first_order_truncation_err_cos(0.5), first_order_truncation_err_sqrt(0.5)

(4.6783240139052204e-8, 1.2990419187857229e-8, 2.8342226290287478e-9, 4.364449135429996e-9)

second_order_truncation_err_exp(0.5), second_order_truncation_err_sin(0.5), second_order_truncation_err_cos(0.5), second_order_truncation_err_sqrt(0.5)

(1.8874426561390482e-8, 7.938850300905947e-9, 4.1240999200086055e-9, 7.45058059692383e-9)

其中：

eps(0.5)=1.1102230246251565e-16

second_order_truncation_err_f()应该是1e-16的顺序，而不是1e-8的顺序，我不知道为什么这不起作用。

Answer 1

这是因为你的计算是由舍入误差决定的。即：

julia> round_off_err_sqrt(0.5)
1.490116119384766e-8

为了看到你的“二阶”导数之间的差异(我通常看到术语中心差异)，你需要选择一个更大的步长。在文献中通常可以看到h = cbrt(eps())。

function second_order_numerical_D(f)
    function df(x)
        h = cbrt(eps(x))          
        (f(x+h/2) - f(x-h/2))/h
    end
    df
end
julia> second_order_error_exp(0.5)
1.308131380994837e-11