4

I implemented a C++ code to numerically solve the n-th derivative of a function in a point 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 );
   }
 }

I was wondering if this is for you a good implementation or if there can be a way to better write it in C++. The problem is that I noticed that the n-th derivative diverges for values of n higher than 3. Do you know how to solve this?

Gianluca Bianco
  • 656
  • 2
  • 11
  • 2
    Perhaps you could implement one of these? https://math.stackexchange.com/q/130192 – konsolas Jan 10 '22 at 22:37
  • Have a look at https://en.wikipedia.org/wiki/Numerical_differentiation for techniques – Ahmed Masud Jan 10 '22 at 22:59
  • 7
    Note `1/3 == 0` – Ranoiaetep Jan 10 '22 at 23:11
  • Obligatory question: do you actually *need* numerical differentiation, or are you only using it because you haven't been given the function `f` in a convenient form for some kind of exact differentiation? – o11c Jan 11 '22 at 00:32
  • If you are looking for feedback (and if your code functions correctly), [Code Review](https://codereview.stackexchange.com/help/on-topic) might be a better choice of venue. On the other hand, if the code does **not** function correctly (as suggested by "The problem is [...]"), then you should focus your question on the malfunction, giving a concrete case along with expected and actual output. – JaMiT Jan 11 '22 at 00:35
  • 1
    `pow( __DBL_EPSILON__, 1/3 );` won't work as expected, and would be very costly to calculate in every function call like that instead of storing in a static constant – phuclv Jan 11 '22 at 00:57
  • @o11c I know the function shape and I just want to obtain its derivative in a point x_0 – Gianluca Bianco Jan 11 '22 at 08:55
  • @phuclv what do you mean with "instead of storing in a static constant"? – Gianluca Bianco Jan 11 '22 at 08:56
  • 1
    @GianlucaBianco I mean as that. Store a constant in a static variable or even better a compile-time constant that's evaluated only once instead of every function call – phuclv Jan 11 '22 at 09:15

1 Answers1

7

It is not a good implementation

At least these problems.

Integer math

Use FP math as 1/3 is zero.

1/3 --> 1.0/3

Using the cube root optimal for n==1

But not certainly other n. @Eugene

Wrong epsilon

Below code is only useful for |x_0| about 1.0. When x_0 is large, x_0 - h may equal x_0. When x_0 is small, x_0 - h may equal -h.

OP's +/- some epsilon is good for fixed point, but double is a floating point.

// Bad
const double h = pow( __DBL_EPSILON__, 1.0/3 );
double x_1 = x_0 - h;

A relative scaling is needed.

#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
}

Invalid code

f is double ( *f )( int, double ), but call is f( x_0 )

Minor: confusing names

Why first_term with x_2 and second_term with x_1?

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • 2
    the `1.0/3` exponent is only optimal for n==1. – Eugene Jan 11 '22 at 00:23
  • 1
    The step `DBL_EPSILON * SCALE` is way too small. A correct choice of the step is a trade-off between roundoff errors and truncation errors. W/o more info on the function, `pow(DBL_EPSILON, 1.0/3)` is recommended instead of `1024*DBL_EPSILON` (for n==1; even higher for n>1). – Eugene Jan 11 '22 at 00:31
  • 2
    Nitpick: Use `cbrt(x)` instead of `pow(x, 1.0/3.0)` – njuffa Jan 11 '22 at 00:32
  • @njuffa `cbrt()` it is. – chux - Reinstate Monica Jan 11 '22 at 00:53
  • @chux-ReinstateMonica thanks, I understand some stuff, but not all. Can you please post a code in which you implement the derivative with your corrections? thanks. – Gianluca Bianco Jan 11 '22 at 09:24
  • @GianlucaBianco The corrections are already posted in the answer. So we can narrow the issue, what parts did you understand? Which part you did not? – chux - Reinstate Monica Jan 11 '22 at 15:36
  • What does it mean ```if (fabs(x_0) >= DBL_MIN && isfinite(x_0))``` and also ```TBD_Code for special cases```? Thanks. – Gianluca Bianco Jan 11 '22 at 15:40
  • 1
    @GianlucaBianco `fabs(x_0)` takes the absolute value of `x_0`. Compare that to `DBL_MIN`, which is the minimum _normal_ `double`. Between `DBL_MIN` and `DBL_TRUE_MIN` are the sub-normals and 0.0. (e.g. 1e-309 to 1e-324). With those wee numbers, code should use a coarser epsilon. `isfinite(x_0)` test if `x_0` is not infinite and not a "not-a-number". For your learner code, you can skip these aside from `x_0 == 0`. For 0.0, try `cbrt(DBL_MIN)`. 0 is special and really has no easy general +/- h that is _best_. – chux - Reinstate Monica Jan 11 '22 at 15:54