| |
- Numerical differentiation using finite differences means the
approximation of a derivative of a function at a given abscissa x using
a linear expression in function values "near" to this abscissa.
- Here we restrict ourselves to the first derivative and to three choices of
formulas:
- ( f(x+h) - f(x) )/h , the forward difference quotient
with error term
(1/2)*h*f''(ξ)
- ( f(x+h) - f(x-h)) /2h , the symmetric difference quotient
with error term
(1/6)*h2f(3)(ξ)
- Richardson extrapolation using the symmetric difference quotient.
This is obtained from the series development
( f(x+h) - f(x-h)) /2h = f'(x)+ ∑i=1 to m h2*i/((2i+1)!)f(2*i+1)(x)
+ O(h2*m+1)
for sufficiently smooth f by combining two, three or more terms of this type linearly in order
to eliminate the powers 2,4,... in h
- Theoretically all these formulas give approximations which converge for grid size h
going to zero to the true derivative value and the actual error is of the form
Ci(h)*horder where Ci(h) is bounded
and order is 1, 2, resp. 2,4,6,8. This requires that the derivatives
of order 2,3, resp. 3,5,7,9 exist and are bounded near x.
In the double logarithmic plot we are using here the order can be seen as the slope
of the error curve.
- Because of the roundoff effects the approximation error of the formula is
superimposed by an irregular behaved error term which grows like Const/h linearly as the
stepsize decreases and which finally dominates the computation. Indeed, if
h has become so small that in the computers arithmetic x+h=x
(for example in IEEE754 double with x=1 and h=2-53)
then any function looks like a constant, since the computed derivative will be zero.
This has the effect that for any such formula there exists an optimal stepsize hopt
which depends on the formulas coefficients and the magnitude of the derivative which appears
in the remainder term. In practice this latter one cannot be estimated without additional effort
and there remains the rule of thumb that
hopt = machine precision 1/(order+1)
with an optimal precision of machine precision order/(order+1).
In the computation shown here h varies in the range 2-i, i=0,...,40.
Both axis are in the logarithm base 10.
For the case of the Richardonextrapolation 4 error curves are displayed, which belong to the orders 2,4,6,8.
These require 2,4,6,8 function evaluations respectively.
For the evaluation of a numerical gradient of a function of n variables one must take
this number of evaluations times n, that means that this process is quite costly.
- There are cases where the above statements seem not to apply, especially in combination with
special inputs x, mostly 0 or 1, where the roundoff behaviour deviates from the usual one,
depending on the internal implementation of a function. And sometimes the results simply appear to
be chaotic. We have used the error value -100 (in the log10 scale) as indicator for "exact".
The "exact" value of the derivative (rounded in machine arithmetic, of course) is computed
"analytically" (i.e. evaluation of a formula for it in machine arithmetic) using a tool
called "automatic differentiation", in this case the code "JAKEF" from netlib.
Hence your input is restricted by this. For example you have no
intrinsic functions other than dexp, dlog, dlog10,
dsin, dcos, datan, dsqrt at your disposal.
|