Dear all, consider the routine fdjac as listed in the NR book

#include <math.h>
#include "nrutil.h"
#define EPS 1.0e-4

void fdjac(int n, float x[], float fvec[], float **df,
void (*vecfunc)(int, float [], float []))
{
int i,j;
float h,temp,*f;
f=vector(1,n);
for (j=1;j<=n;j++) {
temp=x[j];
h=EPS*fabs(temp);
if (h == 0.0) h=EPS;
x[j]=temp+h; Trick to reduce finite precision error.
h=x[j]-temp;
(*vecfunc)(n,x,f);
x[j]=temp;
for (i=1;i<=n;i++) df[i][j]=(f[i]-fvec[i])/h; Forward difference formula.
}
free_vector(f,1,n);
}


It is a routine copied from Netlib. It is used in the MINPACK package for example (the "Trick to reduce finite precision error." is adde by NR authors). Well the problem is in the line

if (h == 0.0) h=EPS;

(considering double prec from now)

If the values of the abscissas x[i] are close to zero, say 1.0e-15,
the product EPS by fabs(temp) is less that machine precision.
For the example h will be 1.0 e-19, that is NOT zero. If the function evaluated by
vecfun has values close to 1.0 the new value f[i] can be identical to old fvec[i] and the jacobian element will be zero.
Just an example. If the function is

f[1] = -2/5 (x1+x2) + 1 + x1
f[2] = -2/5 (x1+x2) + 1 + x2

the real jacobian is

(-2/5 +1 -2/5)
(-2/5 -2/5+1)

if you try to compute the jacobian with the routine above in
x1=1e-20 and x2=1e-20 you will observe the reported issue. The computed jacobian will be wrong, zero values will appear.

The right statement should be

if (h <= EPS) h=EPS;

or something similar. Never check for identity (==) with float or double variable.

Please note that the problem is not in NR only but in MINPACK as well.

Let me know any further opinion/idea.

Regards

Massimo