尝试在MATLAB中实现差异公式 [英] Trying to implement the difference formula in MATLAB

问题描述

我正在尝试实现差异公式

Im trying to implement the difference formula

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

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

using MATLAB for x=1 and h=10^-k, where k=0,...,16. Furthermore, I want to plot the error.

下面是我的代码.我看到错误大约是3，我认为它太大了.它应该接近于0.

Below is my code. I see that the error is around 3, which I believe it too big. It should be close to 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

推荐答案

您的代码中没有任何错误，有限差分公式运行良好，并且您得到的错误在于按数字计算项目:

There is nothing wrong in your code, the finite difference formula works just well, and the error you get lies in compute items following numerically:

当a和b都非常小时，通过计算a/b生成的

• 错误.

• errs generated by calculating a/b when both a and b are very very small.

计算a-b，MATLAB将给出0.

calculating a-b when a and b are very very close that MATLAB will give 0.

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

Here is the result when k changes from 1 to 15:

感谢@CrisLuengo的有见地的评论！

Thanks @CrisLuengo 's insightful comment!

表示err立即下降到几乎为零，并在h变为1e-9时再次上升(此后发生情况1).最后，当h变为1e-14时，df变为0(情况2在这里发生).

which shows that err drops to nearly zero instantly, and rise again when h becomes 1e-9(situation 1 happens after this). Finally df becomes 0 when h becomes 1e-14(situation 2 happens here).

我在您的代码中添加了几行代码来显示这一点:

I added few lines of code to yours to show this:

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

这篇关于尝试在MATLAB中实现差异公式的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！