Numerical Differentiation

Question

How can I calculate the numerical second derivative of a function involving an exponential and a singularity at infinity. Unfortunately, the numerical derivative by Ridder's methods provided in "Numerical Recipes in C" can only calculate the first derivative (It requires analytical expression of the function beforehand.) Furthermore I have tried Chebyshev approximation and differentiating the function afterwards but the values given were way off the actual values. I have also tried some finite difference algorithms provided in a mathematical paper yet they were error prone too. The function is e^(x/2) / x^2. I would appreciate any help on the matter.

Thanks in advance

Latest Edit: The issue was solved the FADBAD libraries available in C++ did an extremely good job. They are available via http://www.fadbad.com/fadbad.html

EDIT:

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

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

#define LIM1 20.0
#define a -5.0
#define b 5.0
#define pre 100.0 // This defines the pre
/* 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);
}

int main(void)
{
    FILE *fp = fopen("Q3data.dat", "w+"), *fp2 = fopen("Q3results.dat", "w+");
    int i; // Declaring our loop variable
    float x, y, min, max, err, nd1, nd2;
    // Define the initial value of the func to be the minimum
    min = func(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;
    for(i = 0, x = a; x < b; i++)
    {   
        x = a + i / pre;
        y = func(x);
        nd1 = dfridr(func, x, 0.1, &err); 
        //nd2 = dfridr((*func), x, 0.1, &err);
        fprintf(fp, "%f \t %f \t %f \t %f \n", x, y, nd1);
        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;
}

EDIT: Chebyshev

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


#include <stdio.h>
#include <math.h>
#include "nr.h"
#define NVAL 150 // Degree of Chebyshev polynomial
#define LIM1 20.0
#define a -5.0
#define b 5.0
#define pre 100.0 // This defines the pre
/* 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);
}

int main(void)
{
    FILE *fp = fopen("Q3data.dat", "w+"), *fp2 = fopen("Q3results.dat", "w+");
    int i; // Declaring our loop variable
    float x, y, min, max;
    float nd1, nd2, c[NVAL], cder[NVAL], cder2[NVAL];
    // Define the initial value of the func to be the minimum
    min = func(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;
    // We make a Chebyshev approximation to our function our interval of interest 
    // The purpose is to calculate the derivatives easily
    chebft(a,b,c,NVAL,func);
    //Evaluate the derivatives
    chder(a,b,c,cder,NVAL); // First order derivative
    chder(a,b,cder,cder2,NVAL); // Second order derivative
    for(i = 0, x = a; x < b; i++)
    {   
        x = a + i / pre;
        y = func(x);
        nd1 = chebev(a,b,cder,NVAL,x);
        nd2 = chebev(a,b,cder2,NVAL,x);

        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

That function is differentiable so using a numeric method is likely not the best. The second derivative is:

6*exp(x/2)/(x^4)-2*exp(x/2)/x^3 + exp(x/2)/(4*x^2)

The above can be simplified of course to speed up computation. Edit: had original formula wrong the first time.