Recursive calculation of second order derivative

Question

I am writing code to do some numerical task using the routines of the book Numerical Recipes. One of my objectives is to calculate the second derivative of a function and have a routine that calculates the first derivative of a function in a nice manner within a specified accuracy. However, I would like to generalize this method and use this function recursively to find the second derivatives or higher order derivatives if needed. I have made some changes in my code and defined a reduced function which takes one argument and achieved my aim more or less. The function that calculates the first derivatives are declared as follows:

float dfridr(float (*func)(float), float x, float h, float *err);

My aim is to set the h value by calling another function for both of the nested dfridr functions at once. For clarity I have enclosed my code. I would appreciate it if you could comment on my method and provide suggestions for improvement and general feedback. As far as I've checked, it works better than the finite difference algorithms, but I think it may be refined.

// The compilation command used is given below
// gcc Q3.c nrutil.c DFRIDR.c -lm -o Q3

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "nr.h"

#define LIM1 20.0
#define a -5.0
#define b 5.0
#define pre 100.0 // This defines the precision
/* This file calculates the func at given points, makes a 
 * plot. It also calculates the maximum and minimum of the func
 * at given points and its first and second numerical derivative.
*/
float func(float x)
{
    return exp(x / 2) / pow(x, 2);
}

// We define the following functions to aid in the calculation of the
// second derivative
float reddfridr(float x)
{
    float err;
    return dfridr(func, x, 0.1, &err);
}
float dfridr2(float x, float h)
{
    float err;
    return dfridr(reddfridr, x, h, &err);
}
int main(void)
{
    FILE *fp = fopen("Q3data.dat", "w+"), *fp2 = fopen("Q3results.dat", "w+");
    int i; // Declaring our loop variable
    float min, max, err, nd1, nd2, x, y;
    // Define the initial value of the func to be the minimum
    min = func(0); 
    // Initialize x
    x = 0;
    for(i = 0; x < LIM1 ; i++)
    {   
        x = i / pre; // There is a singularity at x = 0 
        y = func(x);
        if(y < min)
            min = y;
        fprintf(fp, "%f \t %f \n", x, y);
    }
    fprintf(fp, "\n\n");
    max = 0;
    // Never forget to initialize x again
    x = a;
    // Since i is incremented at the end of the loop
    for(i = 0; x < b  ; i++) 
    {   
        x = a + i / pre;
        y = func(x);
        nd1 = dfridr(func, x, 0.1, &err); 
        nd2 = dfridr2(x, 0.1);
        fprintf(fp, "%f \t %f \t %f \t %f\n", x, y, nd1, nd2);
        if(y > max)
            max = y;
    }
    fprintf(fp2, "The minimum value of f(x) is %f when x is between 0 and 20. \n", min);
    fprintf(fp2, "The maximum value of f(x) is %f when x is between -5 and 5. \n", max);
    fclose(fp);
    fclose(fp2);
    return 0;
}

Answer 1

If you need compile-time controlled order of derivation you probably should play with C++ templates.
If you need run-time (i.e. write own calculator) you'll need to extend that dfridr function (use same approximation of derivation).

Naive implementation and modification can look like:

float dfridr(float (*f)(float), float x, float h)
{
    return f(x+h) - f(x);
}

/* changed */
float dfridr(float (*f)(float), float x, float h, size_t order)
{
  if (order == 0) return f(x);
  else if (order == 1) return f(x+h) - f(x);
  /* extend same algorithm for second derivative */
  else return dfridr(f, x+h, h, order-1) - dfridr(f, x, h, order-1);
}

Answer 2

Recommendations:

1. Your output files are delimited by both spaces and tabs, which is weird. You probably want to use tabs only.

2. There are a couple of points related to floating-point handling that I think deserve commenting:

   • Evaluating func(0.0) has a side effect of setting some floating-point error flags.
   • You actually intend to step x by Δx=0.01 between some lower bound and upper bound. However, repeated addition of 0.01 would result in accumulation of round-off errors, so a workaround is used.
3. pow(x, 2) could be simplified to x * x.

4. Naming could be better. I suggest:

   • pre → inv_delta_x (to emphasize its relationship to x)
   • fp → table, and fp2 → extrema (to explain the contents of the files)
5. Since you are already using C99-style comments, you should also take advantage of C99-style variable declarations, which allow you to declare variables closer to the point of use to enhance readability.

I've defined a FOR_RANGE() macro to make the loops in main() more readable, but that may be a controversial improvement. I've also initialized min and max to INFINITY and -INFINITY, respectively, as a personal preference.

I agree that using recursion will probably yield unsatisfactory results, though I am unable to advise you further. One of the comments above links to a related math question. You will probably get a better answer about numerical methods on Computational Science SE.

#include <math.h>
#include <stdio.h>

// Singularity at x=0.0 raises exception FE_DIVBYZERO and/or sets errno to EDOM
float func(float x)
{
    return exp(x / 2) / (x * x);
}

#define FOR_RANGE(type, var, lower_bound, upper_bound, inv_delta) \
        for (type iteration_dummy = (lower_bound) * (inv_delta), var; \
             (var = iteration_dummy / (inv_delta)) <= (upper_bound); \
             iteration_dummy++)

int main(void)
{
    // We want to increment x in steps of 0.01, but repeatedly adding 0.01
    // would cause round-off errors to accumulate.  Therefore, we will divide
    // successive floating-point integers by 100 instead.
    const float inv_delta_x = 100.0;

    FILE *table   = fopen("Q3data.dat", "w+"),
         *extrema = fopen("Q3results.dat", "w+");

    float min = INFINITY;
    FOR_RANGE(float, x, 0.0, 20.0, inv_delta_x)
    {
        float y = func(x);
        if (y < min) min = y;
        fprintf(table, "%f\t%f\n", x, y);
    }

    fprintf(table, "\n\n");

    float max = -INFINITY;
    FOR_RANGE(float, x, -5.0, 5.0, inv_delta_x)
    {
        float y = func(x);
        float err;
        float nd1 = dfridr(func, x, 0.1, &err); 
        float nd2 = dfridr2(x, 0.1);
        if (y > max) max = y;
        fprintf(table, "%f\t%f\t%f\t%f\n", x, y, nd1, nd2);
    }

    fprintf(extrema, "The minimum value of f(x) is %f when x is between 0 and 20.\n", min);
    fprintf(extrema, "The maximum value of f(x) is %f when x is between -5 and 5.\n", max);
    fclose(table);
    fclose(extrema);
    return 0;
}