尝试在MATLAB中实现差分公式答案

【问题标题】：Trying to implement the difference formula in MATLAB尝试在MATLAB中实现差分公式
【发布时间】：2019-01-05 18:33:30
【问题描述】：

我正在尝试实现差异公式

f'(x) ≈ [ f(x+h) - f(x) ] / h

将 MATLAB 用于x=1 和h=10^-k，其中k=0,...,16。此外，我想绘制错误。

下面是我的代码。我看到误差在 3 左右，我认为它太大了。它应该接近于 0。

syms f(x)
f(x) = tan(x);
df = diff(f,x);
x = 1;
for k = 0:16
    h = 10^-k;
    finitediff = double((f(x+h)-f(x))/h);
    err = double(abs(finitediff-df(x)));
end

【问题讨论】：

标签： matlab difference derivative approximation

【解决方案1】：

您的代码没有任何问题，有限差分公式运行良好，而您得到的错误在于计算以下数字项：

当a 和b 都非常非常小时，计算a/b 产生的错误。
当a 和b 非常接近时计算a-b，MATLAB 将给出0。

这是 k 从 1 变为 15 时的结果：

感谢@CrisLuengo 富有洞察力的评论！

这表明err 立即下降到几乎为零，当h 变为1e-9 时再次上升（情况1 在此之后发生）。最后，df 变为 0，h 变为 1e-14（情况 2 发生在这里）。

我在你的代码中添加了几行代码来展示这一点：

clc;
clear;
format long
syms f(x)
f(x) = tan(x);
h=1;
df = diff(f,x);
double(df(1));
x=1;


range=1:15;
[finitediff,err]=deal(zeros(size(range)));
for k=range
    h=10^-k;
    finitediff(k)=double((f(x+h)-f(x))/h);
    err(k)=double(abs(finitediff(k)-df(1)));
end

figure(1)

subplot(1,2,1)
hold on
plot(err)
plot(err,'bx')
set(gca,'yscale','log')
title('err')

subplot(1,2,2)
hold on
ezplot(df)
axis([0.5 1.5 0 5])
plot(ones(size(range)),finitediff,'rx','MarkerSize',7)
for ii=range
  text(1,finitediff(ii),['h=1e-' num2str(ii)])
end

【讨论】：

不错的情节！在 10^-13 处，误差已经增加了一点，所以大跳跃前的最后两个点已经在增加。您可以通过set(gca,'yscale','log') 看到这一点。实际上，我在k=8 看到了最低限度。

【解决方案2】：

您正在以数字方式计算x+h。使用x=1，您可以制作的大于x 的最小数字是x+eps(x) = 1+eps(1)。 eps(1) 是 2.2204e-16。所以添加h=1e-16 不会改变x。现在，您正在象征性地计算 tan(1)-tan(1)，即 0。因此，您对导数的有限差分近似为 0。

但即使h 较大，您也会得到 0 的差异。我相信这是因为计算正切时出现舍入误差。请注意，

x = 1;
h = 1e-14;
f(x+h)

返回tan(1)。也就是说，符号工具箱认为 tan 函数上下文中的 1+1e-14 足够接近 1 以保证舍入。有趣的是，

f(x+h)-f(x)

返回 0，而

tan(x+h)-tan(x)

返回3.4195e-14（接近预期的3.4255e-14）。

请注意，正如您可以通过 graph that Hunter plots 判断的那样，最佳近似值是 h=10^-8，随着您减少 h，x+h 中的舍入误差开始增加（使用 set(gca,'yscale','log') 来查看此内容） .

【讨论】：