Search code examples
c++floating-pointfloating-accuracy

How to increase accuracy of floating point second derivative calculation?


I've written a simple program to calculate the first and second derivative of a function, using function pointers. My program computes the correct answers (more or less), but for some functions, the accuracy is less than I would like.

This is the function I am differentiating:

float f1(float x) {
    return (x * x);
}

These are the derivative functions, using the central finite difference method:

// Function for calculating the first derivative.

float first_dx(float (*fx)(float), float x) {
    float h = 0.001;
    float dfdx;

    dfdx = (fx(x + h) - fx(x - h)) / (2 * h);
    return dfdx;
}

// Function for calculating the second derivative.

float second_dx(float (*fx)(float), float x) {
    float h = 0.001;
    float d2fdx2;

    d2fdx2 = (fx(x - h) - 2 * fx(x) + fx(x + h)) / (h * h);
    return d2fdx2;
}

Main function:

int main() {
    pc.baud(9600);
    float x = 2.0;

    pc.printf("**** Function Pointers ****\r\n");
    pc.printf("Value of f(%f): %f\r\n", x, f1(x));
    pc.printf("First derivative: %f\r\n", first_dx(f1, x));
    pc.printf("Second derivative: %f\r\n\r\n", second_dx(f1, x));
}

This is the output from the program:

**** Function Pointers ****
Value of f(2.000000): 4.000000
First derivative: 3.999948
Second derivative: 1.430511

I'm happy with the accuracy of the first derivative, but I believe the second derivative is too far off (it should be equal to ~2.0).

I have a basic understanding of how floating point numbers are represented and why they are sometimes inaccurate, but how can I make this second derivative result more accurate? Could I be using something better than the central finite difference method, or is there a way I can get better results with the current method?


Solution

    • Go analytical. ;-) probably not an option given "with the current method".
    • Use double instead of float.
    • Vary the epsilon (h), and combine the results in some way. For example you could try 0.00001, 0.000001, 0.0000001 and average them. In fact, you'd want the result with the smallest h that doesn't overflow/underflow. But it's not clear how to detect overflow and underflow.