# Recursive calculation of second order derivative

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

• Where is dfridr function defined?
– haccks
Jan 12 '14 at 21:33
• @haccks It is defined in another file called DFRIDR.c and the declaration follows from the header file "nr.h". Jan 12 '14 at 21:35
• What's the question? I think this should be in code review
– A Person
Jan 12 '14 at 21:57
• @APerson The question is if refinement of the above code is possible so that the second numerical derivative is evaluated more accurately. Jan 12 '14 at 22:07
• math.stackexchange.com/questions/130192/… : recursion may not be the nor the fastest nor the most accurate method !
– francis
Jan 12 '14 at 22:39

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

• creating the derivatives of user-entered expressions is trivial and thus there is also no need for numerical derivatives.
– user101516
May 12 '16 at 20:28

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:

• preinv_delta_x (to emphasize its relationship to x)
• fptable, and fp2extrema (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;
}

• I would recommend parenthesizing the parameters in the body of the macro, just in case the caller passes x + 1 for lower_bound or whatever. Feb 4 '14 at 0:57
• @GarethRees Yes, of course! Incorporated in Rev 2. Feb 4 '14 at 1:00