# Matlab Numerical Differentiation

1. Mar 19, 2016

### roam

I am trying to estimate the second derivative of $sin(x)$ at $0.4$ for $h=10^{-k}, \ \ k=1, 2, ..., 20$ using:

$$\frac{f(x+h)-2f(x)+f(x-h)}{h^2}, \ \ (1)$$

I know the exact value has to be $f''(0.4)=-sin(0.4)= -0.389418342308651.$

I also want to plot the error as a function $h$ in Matlab.

According to my textbook equation (1) has error give by $h^2f^{(4)}/12$. So I am guessing $f^{(4)}$ refers to the fifth term in the Taylor expansion which we can find using Matlab syntax taylor(f,0.4):

sin(2/5) - (sin(2/5)*(x - 2/5)^2)/2 + (sin(2/5)*(x - 2/5)^4)/24 + cos(2/5)*(x - 2/5) - (cos(2/5)*(x - 2/5)^3)/6 + (cos(2/5)*(x - 2/5)^5)/120

So substituting this in, here is my code for $f''(0.4)$ and error so far:

Code (Text):
format long
x=0.4;

for k = 1:1:20
h=10.^(-k);
ddf=(sin(x+h)-2*sin(x)+sin(x-h))./(h.^2)
e=((h.^2)*((cos(2/5)*(x - 2/5).^3)/6))./12
plot(h,e)
end
But this does not converge and I get some strange values for $f''$ which are very different from the exact value I found above:

Code (Text):
ddf = -0.389093935175844
ddf = -0.389415097166723
ddf = -0.389418309876266
ddf = -0.389418342017223
ddf = -0.389418497448446
ddf = -0.389466237038505
ddf = -0.388578058618805
ddf = -0.555111512312578
ddf = 0.00
ddf = 0.00
ddf = -5.551115123125782e+05
ddf = 0.00
ddf = 0.00
ddf = 0.00
ddf = -5.551115123125782e+13
ddf = 0.00
ddf = 0.00
ddf = 0.00
ddf = 0.00
ddf = 0.00
What is wrong with my code?

And for the error I get all zeros, so it is not possible to plot error vs h. But I believe the graph should look something like this image:

So, what do I need to do here?

Any help would be greatly appreciated.

2. Mar 19, 2016

### BvU

Hello roam,

what I see looks a lot like crossing the limits of what you can do with double precision: by the time $h^2 < 10^{-15}$ things go awry. Not in the denominator, but in the numerator.

3. Mar 20, 2016

### roam

Hi BvU,

Yes, that's true. Do you know by any chance why the error is always equal zero? Am I using the wrong formula or the wrong code?

Last edited: Mar 20, 2016
4. Mar 20, 2016

### BvU

Formula is fine, code is fine too, but not usable for this purpose: the code is hampered by the selfsame issue you try to demonstrate. As far as I can find, Matlab always uses double precision (about 15 significant digits), so once the relative difference becomes too small, you get noise and then zero. Print sin(x-h), sin(x) and sin(x+h) to see it happen.

Study the way double precision numbers are stored to understand why this happens.

You could try to improve precision, but it will only help you a few steps further -- you'll need more than 80 digits of precision for your whole graph...

Last edited: Mar 20, 2016
5. Mar 20, 2016

### phyzguy

The reason the error is always equal to zero is that you are calculating it incorrectly. You have a function$f_{num}(x,h)$ that is calculating the derivative at x, given a value of h. To calculate the error, you should take:
$$abs(f_{num}(x,h) - f_{exact}(x))$$
Your error term has a factor of (x-2/5) in it, and since x = 2/5, of course it gives zero. To calculate the predicted error correctly, the term $f^{(4)}$ refers to the 4th derivative of f, which in this case is sin(x), so the predicted error should be:
$$\frac{h^2 sin(0.4)}{12}$$
It would be interesting to compare this predicted error with your calculated error as calculated above.

6. Mar 20, 2016

### BvU

Phyz is correct, of course. (I wasn't looking any further than ddf -- where the results were already not as expected).

In fact the exercise can be done with excel comfortably as well and it gives the (Phyz) expected result -- up to k = 4 , so when $h^4 = 10^{-16}$

Not unreasonable: we are subtracting $\approx \cos(x+h/2)$ from $\approx \cos(x - h/2)$

7. Mar 20, 2016

### roam

Thank you very much Phyz, I have now managed to calculate $e$ correctly. The values seem reasonable when I compared them with the calculated error (as given by your equation).

Thank you BvU, I think your explanation is why the result seems very erratic, when you have a difference which is close enough to zero.

Anyway, I have now plotted the error function against $h$:

Is it about right? Would you say it looks anything like the generic shape of the error function as illustrated by the figure in my first post?

8. Mar 21, 2016

### BvU

I miss the upswing when h is close to zero. Don't you plot h < 0.001 ?

Code (C):

k   h         | sin(0.4) - ddif |
1 0.1          0.000324407
2 0.01         3.24514E-06
3 0.001        3.24324E-08
4 0.0001       2.91428E-10
5 0.00001      7.10251E-07
6 0.000001     4.78947E-05
7 0.0000001    0.000840284
8 0.00000001   0.16569317

9. Mar 21, 2016

### roam

Why is that? I plotted all the values, here's are the values used:

Code (Text):
e=[3.245152852572088e-04
3.245152852572088e-06
3.245152852572088e-08
3.245152852572088e-10
3.245152852572087e-12
3.245152852572088e-14
3.245152852572088e-16
3.245152852572088e-18
3.245152852572088e-20
3.245152852572088e-22
3.245152852572088e-24
3.245152852572088e-26
3.245152852572088e-28
3.245152852572088e-30
3.245152852572088e-32
3.245152852572087e-34
3.245152852572088e-36
3.245152852572088e-38
3.245152852572087e-40
3.245152852572088e-42];

h=[0.100000000000000
0.010000000000000
1.000000000000000e-03
1.000000000000000e-04
9.999999999999999e-06
1.000000000000000e-06
1.000000000000000e-07
1.000000000000000e-08
1.000000000000000e-09
1.000000000000000e-10
1.000000000000000e-11
1.000000000000000e-12
1.000000000000000e-13
1.000000000000000e-14
1.000000000000000e-15
1.000000000000000e-16
9.999999999999999e-18
1.000000000000000e-18
1.000000000000000e-19
1.000000000000000e-20];

10. Mar 21, 2016

### BvU

Looks like you plotted $\ \frac{h^2 sin(0.4)}{12}\$ ? That is not the error ! That's the predicted error !