0

I want write a script in MATLAB that computes the quotient of the derivative (f(x+h)-f(x))/h for the function x^2 at x=2 and h starting at 1 and decreasing by a factor of 10, and the book said that the effect of the rounding error becomes apparent when h gets to small (less than 10^-12).

Then I wrote the following code:

for i=0:-2:-16
    h=10^i;
    ((2+h)^2-(2)^2)./h
end

Then my question is, How can I improve my code? because It gives me indeed an error saying that the last approximation to the derivative is zero.

Thanks a lot in advance :)

user162343
  • 209
  • 1
  • 3
  • 11
  • 3
    You can start by **not** using i as a variable: http://stackoverflow.com/questions/14790740/using-i-and-j-as-variables-in-matlab – Adriaan Sep 13 '15 at 20:51
  • I am sorry I have moved my question to code review because The thing is that I want a verification :) Thanks – user162343 Sep 13 '15 at 21:04
  • 3
    In code review told me that I have to verify my code here :) what can be done then ? – user162343 Sep 13 '15 at 21:25

1 Answers1

1

Due to the limit of floating point arithmetic you are limited in how small you can choose h. An reasonable option for choosing a "safe" small h is using h=x*sqrt(eps) where eps is the distance from 1.0 to the next largest double-precision number, that is, eps = 2^-52

This Wikipedia page contains some more information

If you REALLY wanted higher precision arithmetic you could consider using an implementation of multi-precision arithmetic. For example this is a MATLAB wrapper around the well-established GMP (GNU Multiple Precision Arithmetic Library) and MPFR (multiple-precision floating-point computations with exact rounding)

gregswiss
  • 1,456
  • 9
  • 20
  • Thanks a lot, the thing is that I don't want to fix the error I want a verification of my code, If my code is right or someone that knows what to expect that tell that my code is doing the things right :) – user162343 Sep 14 '15 at 00:53
  • 1
    looks right. In fact if you develop your expression y = ((x+h)^2-(x)^2)./h = (x^2 + 2*h*x + h^2 - x^2) / h = 2*x + h -> When h is small y=2x which is the derivative of x^2, and is 4 at x=2 – gregswiss Sep 14 '15 at 01:14
  • Thnks :), but what do you mean by "If you develop your expression"? – user162343 Sep 14 '15 at 01:15
  • Right :) then my code is indeed having a rounding error when I less than 10^-12 right? – user162343 Sep 14 '15 at 01:22
  • Ok thanks a lot, That is what I was searching for, a verification, that someone could tell me if it was right or I was missing something and how to fix it :) – user162343 Sep 14 '15 at 01:30