2009-10-13 31 views
15

¿Cómo se calcula programáticamente la derivada de f(x) para garantizar la máxima precisión?Implementación de la derivada en C/C++

Estoy implementando el método Newton-Raphson, y requiere tomar la derivada de una función.

+0

Muéstranos lo que has hecho hasta ahora. – Lazarus

+4

¿quieres que me despidan :)? – vehomzzz

+5

La precisión del método de Newton no depende (solo) de la precisión de la derivada, una aproximación aproximada también lo hará (en el caso extremo se obtiene el método secante), pero será un poco más lento obtener la misma precisión (más iteraciones necesarias). – fortran

Respuesta

51

Estoy de acuerdo con erikkallen que (f (x + h) - f (x - h))/2h es el método habitual para aproximar numéricamente derivadas. Sin embargo, obtener el tamaño de paso correcto h es un poco sutil.

El error de aproximación en (f (x + h) - f (x - h))/2h disminuye a medida que h se hace más pequeño, lo que indica que debe tomar h lo más pequeño posible. Pero a medida que h se hace más pequeño, el error de la resta de coma flotante aumenta ya que el numerador requiere restar números casi iguales. Si h es demasiado pequeño, puede perder mucha precisión en la resta. Por lo tanto, en la práctica, debe elegir un valor no demasiado pequeño de h que minimice la combinación de error de aproximación y error numérico.

Como regla general, puede probar h = SQRT (DBL_EPSILON) donde DBL_EPSILON es el número de precisión doble más pequeño e de manera que 1 + e! = 1 en la precisión de la máquina. DBL_EPSILON es aproximadamente de 10^-15, por lo que podría usar h = 10^-7 o 10^-8.

Para obtener más detalles, consulte estos notes al elegir el tamaño de paso para ecuaciones diferenciales.

+0

+1 Muy interesante, investigando eso. – vehomzzz

+4

+1 para mencionar y explicar el dilema del tamaño de paso – sellibitze

+3

Creo que su regla empírica supone que utiliza una regla de primer orden para aproximar la derivada. Sin embargo, la regla de diferencia central que mencionas es de segundo orden, y la regla de oro correspondiente es h = EPSILON^(1/3) que es aproximadamente 10^(- 5) cuando se usa precisión doble. –

5
fprime(x) = (f(x+dx) - f(x-dx))/(2*dx) 

para algunos dx pequeños.

+2

¿Qué tan pequeño? gracias – vehomzzz

+1

Recetas numéricas tiene algunos comentarios sobre eso http://books.google.co.uk/books?id=1aAOdzK3FegC&printsec=frontcover&dq=related:ISBN0521437202#v=onepage&q=newton-Raphson&f=false – Mark

10

Newton_Raphson supone que puede tener dos funciones f (x) y su derivada f '(x). Si no tiene la derivada disponible como una función y tiene que estimar la derivada de la función original, entonces debe usar otro algoritmo de búsqueda de raíz.

Wikipedia root finding ofrece varias sugerencias como cualquier texto de análisis numérico.

5

¿Qué sabes acerca de f (x)? Si solo tiene f como caja negra, lo único que puede hacer es aproximar numéricamente la derivada. Pero la precisión generalmente no es tan buena.

Puede hacer mucho mejor si puede tocar el código que calcula f. Pruebe "automatic differentiation". Hay algunas buenas bibliotecas para eso disponible. Con un poco de magia de biblioteca puede convertir fácilmente su función a algo que calcule la derivada automáticamente. Para un ejemplo simple de C++, consulte la discusión alemana source code in this.

2

Además de la respuesta anterior de John D. Cooks, es importante no solo tener en cuenta la precisión del punto flotante, sino también la robustez de la función f (x). Por ejemplo, en finanzas, es un caso común que f (x) es en realidad una simulación de Monte Carlo y el valor de f (x) tiene algo de ruido. Usar un tamaño de paso muy pequeño puede en estos casos degradar severamente la precisión de la derivada.

8

alt text

alt text

1) Primer caso:

alt text

alt text - error de redondeo relativa, aproximadamente 2^{- 16} para el doble y 2^{- 7} para flotar.

Podemos calcular el error total:

alt text

Suponga que está utilizando doble operación flotante. Por lo tanto, el valor óptimo de h es 2sqrt (DBL_EPSILON/f '' (x)). Usted no sabe f '' (x). Pero debes estimar este valor. Por ejemplo, si f '' (x) es aproximadamente 1, entonces el valor óptimo de h es 2^{- 7} pero si f '' (x) es aproximadamente 10^6, entonces el valor óptimo de h es 2^{- 10}!

2) Segundo caso:

alt text

Tenga en cuenta que el segundo error de aproximación tiende a 0 más rápido que el primero. Pero si f '' '(x) es muy lagre entonces primera opción es más preferible:

alt text

Tenga en cuenta que en el primer caso h es proporcional a e pero en el segundo caso h es proporcional a e^{1/3}. Para operaciones flotantes dobles e^{1/3} es 2^{- 5} o 2^{- 6}. (Supongo que f '' '(x) es aproximadamente 1).


¿Qué camino es mejor? No se conoce si no conoce f '' (x) y f '' '(x) o no puede estimar estos valores. Se cree que la segunda opción es preferible. Pero si sabes que f '' '(x) es muy grande usa primero.

¿Cuál es el valor óptimo de h? Supongamos que f '' (x) y f '' '(x) son aproximadamente 1. Supongamos también que utilizamos operaciones flotantes dobles. Luego, en el primer caso h es aproximadamente 2^{- 8}, en el primer caso h es aproximadamente 2^{- 5}. Corrija estos valores si conoce f '' (x) o f '' '(x).

+0

épsilon debe ser 2^-53 para el doble, y 2^-24 para el flotante (que es aproximadamente 10^-16 y 10^-7, respectivamente). –

+0

epsilon es ** relativo ** error de redondeo (no absoluto). Siempre es aproximadamente 10^{- 16} para el doble y 10^-7 para el flotador –

+0

Sí, lo sé. En su respuesta, dice "épsilon - error de redondeo relativo, alrededor de 2^{- 16} para el doble y 2^{- 7} para el flotante," que es claramente incorrecto. El error de redondeo relativo (hacia adelante) también es ** no ** siempre en esa escala, sino más bien el error hacia atrás. El error de reenvío puede ser mucho, mucho mayor cuando ocurre la cancelación, como es probable que ocurra aquí. –

4

Definitivamente, debe tener en cuenta la sugerencia de John Cook para elegir h, pero normalmente no desea usar una diferencia centrada para aproximar la derivada. La razón principal es que cuesta una evaluación de la función adicional, si se utiliza una diferencia hacia adelante, es decir,

f'(x) = (f(x+h) - f(x))/h 

entonces obtendrá el valor de f (x) de forma gratuita porque es necesario para calcular ya de Newton método. Esto no es tan importante cuando tienes una ecuación escalar, pero si x es un vector, entonces f '(x) es una matriz (el jacobiano), y tendrás que hacer n evaluaciones de funciones adicionales para aproximarlo usando el enfoque de diferencia centrada.

1

Normalmente, el ruido de la señal afecta la calidad de las derivaciones más que cualquier otra cosa. Si tiene ruido en su f (x), Savtizky-Golay es un excelente algoritmo de suavizado que a menudo se utiliza para calcular buenos derivados. En pocas palabras, SG ajusta un polinomio localmente a sus datos, entonces este polinomio se puede usar para calcular la derivada.

Paul