Implementando el derivado en C/C++

algorithm math

Question

¿Cómo se calcula normalmente la derivada de una f(x) para asegurar la máxima precisión?

Estoy implementando el Newton-Raphson , y requiere que se tome la derivada de una función.

Answer 1

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

Puede hacerlo mucho mejor si puede tocar el código que calcula f. Prueba la "diferenciación automática" . Hay algunas bibliotecas agradables para eso disponibles. Con un poco de magia de biblioteca puede convertir su función fácilmente a algo que calcula la derivada automáticamente. Para un ejemplo simple de C ++, vea el código fuente en esta discusión alemana.

Answer 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. El uso de un tamaño de paso muy pequeño en estos casos puede degradar gravemente la precisión del derivado.

Answer 3

Definitivamente desea tener en cuenta la sugerencia de John Cook para seleccionar h, pero normalmente no quiere usar una diferencia centrada para aproximar la derivada. La razón principal es que cuesta una evaluación de función adicional, si usa una diferencia de avance, es decir,

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

Luego obtendrás el valor de f (x) de forma gratuita porque ya debes computarlo para el método de Newton. 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 necesitarás hacer n evaluaciones de funciones adicionales para aproximarlo utilizando el enfoque de la diferencia centrada.

Answer 4

Estoy de acuerdo con @erikkallen en que (f(x + h) - f(x - h)) / 2 * h es el enfoque habitual para las aproximaciones numéricas. 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)) / 2 * h disminuye a medida que h se reduce, lo que indica que debe tomar h más pequeño posible. Pero a medida que h se reduce, el error de la resta de punto 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. Entonces, en la práctica, debe elegir un valor no muy 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 doble precisión más pequeño e tal que 1 + e != 1 en la precisión de la máquina. DBL_EPSILON es de aproximadamente 10^-15 así que puedes usar h = 10^-7 o 10^-8 .

Para obtener más detalles, consulte estas notes sobre cómo elegir el tamaño de paso para las ecuaciones diferenciales.

Answer 5

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

El hallazgo de la raíz de Wikipedia ofrece varias sugerencias, como cualquier texto de análisis numérico.

Answer 6

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

Pablo

Answer 7

1) Primer caso:

- error relativo de redondeo, aproximadamente 2 ^ {- 16} para el doble y 2 ^ {- 7} para el flotador.

Podemos calcular el error total:

Supongamos que está utilizando la operación de doble flotación. Por lo tanto, el valor óptimo de h es 2sqrt (DBL_EPSILON / f '''' (x) ). Usted no sabe f '''' (x) . Pero tienes que 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:

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

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 de doble flotación e ^ {1/3} es 2 ^ {- 5} o 2 ^ {- 6}. (Supongo que f '''' ''(x) es aproximadamente 1).

¿Qué camino es mejor? Se desconoce 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, utiliza primero.

¿Cuál es el valor óptimo de h? Supongamos que f '''' (x) y f '''' ''(x) son aproximadamente 1. Supongamos también que usamos operaciones de doble flotación. Entonces, 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).

Answer 8

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

Para algunos pequeños dx.