Implementé un código C++ para resolver numéricamente la n-ésima derivada de una función en un punto x_0:
double n_derivative( double ( *f )( double ), double x_0, int n ) { if( n == 0 ) return f( x_0 ); else { const double h = pow( __DBL_EPSILON__, 1/3 ); double x_1 = x_0 - h; double x_2 = x_0 + h; double first_term = n_derivative( f, x_2, n - 1 ); double second_term = n_derivative( f, x_1, n - 1); return ( first_term - second_term ) / ( 2*h ); } }Me preguntaba si esto es para usted una buena implementación o si puede haber una forma de escribirlo mejor en C++. El problema es que noté que la n-ésima derivada diverge para valores de n mayores a 3. ¿Sabes cómo resolver esto?
No es una buena implementación.
Al menos estos problemas.
matemáticas enteras
Use matemáticas FP ya que 1/3 es cero.
1/3 --> 1.0/3
Usando la raíz cúbica óptima para n==1
Pero no ciertamente otro n . @eugene
épsilon incorrecto
El siguiente código solo es útil para |x_0| alrededor de 1.0. Cuando x_0 es grande, x_0 - h puede ser igual a x_0 . Cuando x_0 es pequeño, x_0 - h puede ser igual a -h .
OP's +/- some epsilon es bueno para el punto fijo , pero double es un punto flotante .
// Bad const double h = pow( __DBL_EPSILON__, 1.0/3 ); double x_1 = x_0 - h;Se necesita una escala relativa .
#define EPS cbrt(DBL_EPSILON) // TBD code to well select this if (fabs(x_0) >= DBL_MIN && isfinite(x_0)) { double x_1 = x_0*(1.0 - EP3); double x_2 = x_0*(1.0 + EPS); double h2 = x_2 - x_1; ... } else { TBD_Code for special cases }Codigo invalido
f es double ( *f )( int, double ) , pero call es f( x_0 )
Menor: nombres confusos
¿Por qué first_term con x_2 y second_term con x_1 ?