patterncMinor

Recursive calculation of second order derivative

Submitted by: @import:stackexchange-codereview·Mar 10, 2026·

Viewed 0 times

derivativecalculationorderrecursivesecond

Problem

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
#include
#include
#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,

Solution

Recommendations:

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

-
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.

pow(x, 2) could be simplified to x * x.

-
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)

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

// 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))  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;
}

Code Snippets

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

Context

StackExchange Code Review Q#39106, answer score: 4

Revisions (0)

No revisions yet.