Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab Numerical Differentiation

  1. Mar 19, 2016 #1
    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? :confused:

    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:

    2016_03_19_19_15_33.jpg

    So, what do I need to do here?

    Any help would be greatly appreciated.
     
  2. jcsd
  3. Mar 19, 2016 #2

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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.
     
  4. Mar 20, 2016 #3
    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
  5. Mar 20, 2016 #4

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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
  6. Mar 20, 2016 #5

    phyzguy

    User Avatar
    Science Advisor

    The reason the error is always equal to zero is that you are calculating it incorrectly. You have a function[itex]f_{num}(x,h)[/itex] that is calculating the derivative at x, given a value of h. To calculate the error, you should take:
    [tex] abs(f_{num}(x,h) - f_{exact}(x))[/tex]
    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 [itex] f^{(4)}[/itex] refers to the 4th derivative of f, which in this case is sin(x), so the predicted error should be:
    [tex] \frac{h^2 sin(0.4)}{12}[/tex]
    It would be interesting to compare this predicted error with your calculated error as calculated above.
     
  7. Mar 20, 2016 #6

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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}##

    upload_2016-3-20_23-14-53.png
    Not unreasonable: we are subtracting ##\approx \cos(x+h/2)## from ## \approx \cos(x - h/2)##
     
  8. Mar 20, 2016 #7
    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##:

    error.jpg

    er2.jpg

    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?
     
  9. Mar 21, 2016 #8

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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
     
     
  10. Mar 21, 2016 #9
    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];
     
  11. Mar 21, 2016 #10

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Looks like you plotted ##\
    \frac{h^2 sin(0.4)}{12}\ ## ? That is not the error ! That's the predicted error !
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Numerical Differentiation
  1. Numerically Solve (Replies: 2)

  2. Numerical Integration (Replies: 3)

  3. Numerical Jacobian (Replies: 1)

Loading...