The numdiff function is a quasi-duplicate for derivative
Reported by Michael BAUDIN
Originally assigned to Paul BIGNIER
-- Bug description --
The "numdiff" primitive and the "derivative" primitives
both compute the numerical derivatives of a given function.
But their source code is not shared so that difference are to
expect.
The "derivative" macro has a step which is computed with
formula :
case 1 , h = sqrt(%eps)
case 2 , h = %eps^(1/3)
case 4 , h = %eps^(1/4)
which is not good if x is lower than %eps, because in that case x+h may
be so that h dominates x so that sqrt(%eps)*abs(%x) should be
prefered.
The "numdiff" macro has the problem that
the step is computed with the sub-optimal formula
%dx=sqrt(%eps)*(1+1d-3*abs(%x))
That formula is good when x is lower than %eps, because in that
case the scaling will be close to 1. But if x is large, the
value of the step is not optimal (but reasonable),
and sqrt(%eps)*abs(%x) should be prefered.
The conclusion is that
* numdiff and derivative do not produce the same results
* the user has no option to control the strategy to compute the
automatic step
* the different strategies are not consistent
* no strategy is better than the other
* the documentation lacks of details
-- Scilab error message --
-- How to reproduce the bug --
function y = myfunction (x)
y=x^3;
mprintf("myfunction(%e)=%e\n",x,y);
endfunction
format(25);
// numdiff is better for small values of x (5 digits more accurate)
x = 10^-32;
numdiff(myfunction, x) // 0.0000000000000002220446 (exact = 0.0)
derivative(myfunction, x) // 0.0000000000366685286250 (exact = 0.0)
// numdiff is better for large values of x ...
x = 10^32;
numdiff(myfunction, x) // 2.999996991115918983D+64 (exact = 3*10^64)
derivative(myfunction, x) // 0.0 (exact = 3*10^64)
// ... but it could be improved (with 2 more digits)
h=sqrt(%eps)*abs(x);
numdiff(myfunction, x,h) // 3.000000038712711568D+64 (exact = 3*10^64)
// derivative is better for medium values of x (3 more digits)
x = 1.0;
numdiff(myfunction, x) // 3.0000000507324253717911 (exact = 3.0)
derivative(myfunction, x) // 3.0000000000248241427414 (exact = 3.0)