2

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
Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
Jimmy Sabater
  • 127
  • 3

2 Answers2

4

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
Hunter Jiang
  • 1,300
  • 3
  • 14
  • 23
  • 1
    Nice plots! At 10^-13, the error increases a little already, so the last two dots before the big jump are already increasing. You can see this with `set(gca,'yscale','log')`. Actually I see a minimum at `k=8`. – Cris Luengo Jan 05 '19 at 10:18
4

You are computing x+h numerically. with x=1, the smallest number larger than x that you can make is x+eps(x) = 1+eps(1). eps(1) is 2.2204e-16. So adding h=1e-16 doesn't change x. Now, you are symbolically computing tan(1)-tan(1), which is 0. Hence your finite difference approximation to the derivative is 0.

But even with a larger h you get differences of 0. I believe this is because of rounding errors that occur in computing the tangent. Note that 

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

returns tan(1). That is, the symbolic toolbox thinks of 1+1e-14 within the context of the tan function as sufficiently close to 1 to warrant rounding. Funnily enough,

f(x+h)-f(x)

returns 0, whereas

tan(x+h)-tan(x)

returns 3.4195e-14 (which is close to the expected 3.4255e-14).

Note that, as you can tell by the graph that Hunter plots, the best approximation is for h=10^-8, as you decrease h the rounding errors in x+h start to increase (use set(gca,'yscale','log') to see this).

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120