Search code examples
matlabdifferencederivativeapproximation

Trying to implement the difference formula in MATLAB

Im trying to implement the difference formula 

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

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

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
Solution

  • 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:

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

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

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

    enter image description here

    Thanks @CrisLuengo 's insightful comment!

    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