[Math] Numerical Differentiation Matlab

derivativesMATLABnumerical methods

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}\tag{1}$$

And then plot the error as a function of $h$ in Matlab.

Attempt:

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

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:

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 I get all zeros for the error. And the estimation does not look right after the 7th computation:

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

So, what is wrong with my code? And why does the errors equal zero? I thought the error vs. steps size plot should look something like This.

Any help is greatly appreciated.

Best Answer

For $f(x)=\sin(x)$ you get the nice formula $$ f^{(k)}(x)=\sin\Bigl(x+k\frac\pi2\Bigr) $$ Thus $f^{(4)}(x)=\sin(x)$.

In the error formula, the argument of the 4th derivative is some not known point inside the interval $(x-h,x+h)$. As a first estimate one can take the value at $x$ if $h$ is small.

Apart from the fact that it has no place in the formula, (x-2/5) for $x=0.4$ evaluates to zero. The correct error estimate is, following the formula

 e = sin(x)*h^2/12

(Aside, if you do not use vectorized operations, you do not need the point modificator.)

The evaluation of $f(x\pm h)$ deviates from its exact value by the accumulated errors of the floating point operations. First there is an error of magnitude $|x|\mu$ in the addition $x\pm h$ for $|h|\ll|x|$ which translates into an error contribution of $|f'(x)|·|x|\mu$ in the value. Then the evaluation of $f$ itself will contribute $|f(x)|·m_f\mu$ for some $m_f$ indicating the evaluation complexity of $f$ and assuming that there are no great cancellations in this evaluation algorithm of $f$ ($m_{\sin}=1$ as it is one FPU operation). The complete error bound estimate combines the floating point errors and the discretization errors to something like $$ \frac{2(m_f|f(x)|+|f'(x)x|)\mu}{h^2}+\frac{|f^{(4)}(x)|}{12}h^2 $$ The minimum of that bound is reached when both terms are equal which happens at about $h\sim \sqrt[4]\mu\simeq 10^{-4}$. Thus the best derivative value is the 4th with about 7-9 correct digits, which visual inspection confirms.

On the other hand, as example, the 8th value has an error of $\sim 10^{-16}·(10^8)^2=1$ from the first term, the floating point errors, thus the error is larger than the expected value, giving a totally useless result.

The exact values obey $f(x+h)=f(x)+f'(x)h+0.5f''(x)h^2+…$. For $h^2<\mu$, that is, for $|h|<10^{-8}$ for double, the second degree term falls below the threshold of rounding errors and thus plays no role in the difference formula. And since the second difference of a linear function is zero, this explains the zero results in the list of numerical results.

Related Question