Diferenciación numérica Matlab

Question

Estoy tratando de calcular la segunda derivada de $\sin(x)$ $0.4$ $h=10^{-k}, \ \ k=1, 2, ..., 20$ el uso de:

$$\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)=-\pecado(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)$ e de 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

Pero tengo todos los ceros para el error. Y la estimación no mira a la derecha después de las 7 de la computación:

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

Así que, ¿qué hay de malo en mi código? Y ¿por qué los errores igual a cero? Pensé que el error vs pasos el tamaño de la parcela debe ser algo como Esto.

Cualquier ayuda es muy apreciada.

Answer 1

Para $f(x)=\sin(x)$ que se obtenga una bonita fórmula $$ f^{(k)}(x)=\sin\Bigl(x+k\frac\pi2\Bigr) $$ Por lo tanto $f^{(4)}(x)=\sin(x)$.

En el error de la fórmula, el argumento de la 4ª derivada es que algunos no se conoce el punto dentro del intervalo de $(x-h,x+h)$. Como primera aproximación se puede tomar el valor en $x$ si $h$ es pequeña.

Aparte del hecho de que no tiene lugar en la fórmula, (x-2/5) $x=0.4$ evalúa a cero. La correcta estimación del error es, siguiendo la fórmula

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

(Aparte, si no uso vectorizados de operaciones, no es necesario el punto de modificator.)

---

La evaluación de las $f(x\pm h)$ se desvía de su valor exacto por el acumulado de los errores de las operaciones de punto flotante. Primero hay un error de magnitud $|x|\mu$ en la adición $x\pm h$ $|h|\ll|x|$ lo que se traduce en un error de la contribución de $|f'(x)|·|x|\mu$ en el valor. A continuación, la evaluación de las $f$ sí va a contribuir $|f(x)|·m_f\mu$ algunos $m_f$, lo que indica la evaluación de la complejidad de $f$, y suponiendo que no hay grandes cancelaciones en esta evaluación algoritmo de $f$ ($m_{\sin}=1$ como es una FPU operación). El mensaje de error completo obligado estimación combina los errores de punto flotante y la discretización de los errores a algo como $$ \frac{2(m_f|f(x)|+|f'(x)x|)\mu}{h^2}+\frac{|f^{(4)}(x)|}{12}h^2 $$ El mínimo de que límite se alcanza cuando ambos términos son iguales que ocurre alrededor de $h\sim \sqrt[4]\mu\simeq 10^{-4}$. Por lo tanto la mejor de derivados valor es el 4º con alrededor de 7-9 correcto de dígitos, de forma que la inspección visual confirma.

Por otro lado, como ejemplo, el 8 de valor, tiene un error de $\sim 10^{-16}·(10^8)^2=1$ desde el primer plazo, los errores de punto flotante, por lo tanto el error es mayor que el valor esperado, dando un totalmente inútil resultado.

---

Los valores exactos de obedecer $f(x+h)=f(x)+f'(x)h+0.5f''(x)h^2+…$. Para $h^2<\mu$, es decir, para $|h|<10^{-8}$ double, el segundo grado, término cae por debajo del umbral de errores de redondeo y por lo tanto no juega ningún papel en la diferencia de la fórmula. Y desde la segunda diferencia de una función lineal es cero, esto explica el cero resultados en la lista de los resultados numéricos.

Answer 2

Sugerencia Conseguir efectos de cancelación. Una diferencia que es bastante cercana a cero perderá precisión casi todos que los términos podrían tener. Un ejemplo: decir que tiene precisión de 16 bits y resta dos números que son de igual magnitud e igualan en los primeros 12 bits. Esto deja sólo 4 bits posible de izquierda de precisión (los bits menos significativos). El problema rápido obtiene catastrófico como usted comienza a acercarse a todos los bits de la sustracción se cancelan hacia fuera.