时间戳

首页 > 解决方案 > 尝试在MATLAB中实现差分公式

matlab - 尝试在MATLAB中实现差分公式

问题描述

我试图实现差异公式

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

使用 MATLABx=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

标签： matlabdifferencederivativeapproximation

解决方案

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

  • a/b当两者a和b都非常小时，通过计算产生的错误。

  • 计算MATLAB 给出的a-b时间a和b非常接近的时间0。

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

在此处输入图像描述

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

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

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

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

推荐阅读