In his article "Algorithms for High-Precision Finite Differences" Michael Herstine describes a convoluted procedure which fails to achieve its objectives on a simple cubic polynomial function. There is clearly something wrong with the method. Whilst it is true that computing numerical derivatives is rather like making a pact with the devil, if there are no alternatives then it is important for the numerical analyst to do his best with the available data. I would normally use either golden ratio or simplex optimisation for functions where analytic derivatives are not available. These are very robust methods.
However, returning to the problem in hand which is to compute the accurate numerical derivative of a given external function. The method for doing this has been known since before the era of the digital computer (in the days when tabular data was common). We must assume the given function is continuous and analytic over the range under consideration. The classic Taylor expansion of the function is
f(x0+h) = f(x0) + f'(x0).h + f"(x0).h²/2! + f"'(x0).h³/3! + ...
An expression for the derivative f'(x0) may be obtained from this by use of
finite difference algebra beyond the scope of a short note. See for example:
Milne-Thomson, The Calculus of Finite Differences (MacMillan, London 1933) or
Froberg, Introduction to Numerical Analysis (Wiley, London, 1965)
We start by choosing a suitable small step size
h =
|x|.sqrt(eps(M))
where
eps(M) = 10^-7 for single precision
eps(M) = 10^-16 for double precision
And then tabulating the given function at x-2h, x-h, x, x+h, x+2h and
computing the various differences as follows:
| Function | Diff 1 | Diff 2 | Diff 3 | Diff 4 |
| f(x-2h) | ||||
| D1(-3h/2) | ||||
| f(x-h) | D2(-h) | |||
| D1(-h/2) | D3(-h/2) | |||
| f(x) | D2(0) | D4(0) | ||
| D1(h/2) | D3(h/2) | |||
| f(x+h) | D2(h) | |||
| D1(3h/2) | ||||
| f(x+2h) |
Where the forward differences
D1(-h/2) = f(x)-f(x-h),
D2(0) = D1(h/2)-D1(-h/2) etc.
Some fiddly algebra leads to two useful results:
1. Symmetric difference derivative which requires 5 function
evaluations:
f' = 1/(2h).[ f(x+h)-f(x-h) - 1/6.(D2(h)-D2(-h)) +
1/30.(D4(h)-D4(-h)) - ...]
2. Assymetric difference derivative which requires 4 function
evaluations:
f' = 1/h.[ f(x+h)-f(x) -1/6.D2(h) -1/3.D2(0) + 1/30.D4(h) + 1/20.D4(h)
- ...]
The expressions here are shown upto and including fourth
order terms. Usually for single precision data only second order terms are
needed. The second order term in D2 correcting the initial raw estimate. If the
function is double precision it is worth considering the fourth order
corrections when 6 or 7 function evaluations will be required. Each of these
formulae has advantages depending on how difficult your function is to evaluate,
and the location of any nearby singularities. In general the symmetric
difference form should be preferred as it usually results in a much smaller
second order correction. Practical implementation of these expressions to
compute derivatives would require some additional management code, but to
illustrate the performance I will re-examine the "difficult" cubic
polynomial.
| x | Analytic f'(x) | f(x) = (x-100)²+1e-6*(x-300)³ | Trunc(6dp) | delta 1 | delta 2 |
| 98.001 | -3.8755892 | -4.24628458860602 | -4.246284 | ||
| -2.876194 | |||||
| 99.001 | -1.8767982 | -7.12247879760301 | -7.122478 | 1.998793 | |
| -0.877401 | |||||
| 100.001 0 | 0.1219988 | -7.9998790006000 | -7.999879 | 1.998801= | |
| 1.121400 | |||||
| 101.001 | 2.1208018 | -6.87847919759699 | -6.878479 | 1.998806 | |
| 3.120206 | |||||
| 102.001 | 4.1196108 | -3.75827338859398 | -3.758273 |
Using a step length of h = 1 we have after truncating to 6 decimal places.
Symmetric difference version 5 function evaluations
| (f(x+h)-f(x-h))/2h | 0.1219995 |
| Correction -(D2(h)-D2(-h))/12h | -0.0000011 |
| 0.1219984 |
Asymmetric difference version 4 function evaluations
| (f(x+h)-f(x))/h | 1.1214000 |
| Correction -(D2(h)+2.D2(0))/6h | -0.9994013 |
| 0.1219987 |
Comparison of symmetric and asymmetric difference derivatives with step length h
| f'(x | Symmetric Difference | Asymmetric Difference | ||
| h | D/2h | f'(x) | D1/h | f'(x) |
| 0.000001 | 0.5000000 | 0.5833333 | 1.0000000 | 0.833333 |
| 0.00001 | 0.1500000 | 0.1583333 | 0.2000000 | 0.183333 |
| 0.0001 | 0.1250000 | 0.1258333 | 0.1300000 | 0.128333 |
| 0.001 | 0.1220000 | 0.1220000 | 0.1230000 | 0.122000 |
| 0.01 | 0.1220000 | 0.1220000 | 0.1320000 | 0.122000 |
| 0.1 | 0.1219950 | 0.1219942 | 0.2219400 | 0.121996 |
| 1 | 0.1219995 | 0.1219984 | 1.1214000 | 0.121998 |
| 10 | 0.1220988 | 0.1219988 | 10.1160988 | 0.121998 |
| 100 | 0.1319988 | 0.1219988 | 100.0719991 | 0.121998 |
| 1000 | 1.1219988 | 0.1219988 | 1000.5220018 | 0.121998 |
Analytic target value f'(x) = 0.12199880000301
The table shows results of numerical derivatives as a function of step size
Symmetric difference gives decent results (f(x+h)-f(x-h))/2h when
0.001 < x < 5
f' sym corrected when 0.001 < x
Asymmetric
difference gives decent results (f(x+h)-f(x))/h when 0.001 = x (still very
poor)
f' asym corrected when 0.001 < x
Note how for the smooth cubic polynomial function the results are exact to available machine precision after correction with both methods. And comparing the uncorrected forward difference value 1.1214 for h=1 with results in Table 2 of the article the systematic error is revealed. I deduce the asymmetric difference approximation was used uncorrected. The answers computed in the article are approximate values of f'(x+h/2) The sin(x) test also presents no problems for the corrected derivative expressions.
To test the algorithm aggressively consider evaluating the derivative of the pathological function sin(1/x) at x = 0.11 Table data shown for h=0.012 to emphasise the effects of corrections.
| x | Analytic f'(x) | f(x) =sin(1/x) | Trunc(6dp) | delta 1 | delta 2 |
| 0.1076 | 85.6313313 | 0.13072246547674 | 0.130722 | ||
| 0.100760 | |||||
| 0.1088 | 82.1832274 | 0.23148268739985 | 0.231482 | -0.004542 | |
| 0.096218 | |||||
| 0.11 | 78.0811228 | 0.32770070881350 | 0.327700 | -0.005256 | |
| 0.090962 | |||||
| 0.1112 | 73.4419339 | 0.41866265767427 | 0.418662 | -0.005834 | |
| 0.085128 | |||||
| 0.1124 | 68.3744253 | 0.50379013846981 | 0.503790 |
Symmetric difference version 5 function evaluations
| (f(x+h)-f(x-h))/2h | 77.9916667 |
| Correction -(D2(h)-D2(-h))/12h | 0.0897222 |
| 78.0813889 |
Asymmetric difference version 4 function evaluations
| (f(x+h)-f(x))/h | 75.8016667 |
| Correction -(D2(h)+2.D2(0))/6h | 2.2702778 |
| Value | 78.0719444 |
Comparison of symmetric and asymmetric difference derivatives with step length h
| Symmetric Difference | Asymmetric Difference | |||
| h | D/2h | f'(x) | D1/h | f'(x) |
| 0.000001 | 78.000000 | 78.000000 | 78.000000 | 78.000000 |
| 0.00001 | 78.100000 | 78.108333 | 78.100000 | 78.116667 |
| 0.0001 | 78.080000 | 78.080000 | 77.900000 | 78.081667 |
| 0.0002 | 78.080000 | 78.082917 | 77.715000 | 78.083333 |
| 0.0005 | 78.066000 | 78.081667 | 77.152000 | 78.081000 |
| 0.001 | 78.019000 | 78.081167 | 76.193000 | 78.075833 |
| 0.002 | 77.832500 | 78.083542 | 74.196500 | 78.040083 |
| 0.005 | 76.508600 | 78.156217 | 67.703400 | 77.488200 |
| 0.01 | 71.565750 | 78.918567 | 55.959400 | 74.029117 |
| 0.05 | 7.852680 | 11.290922 | -7.217580 | 4.844980 |
The table shows results of numerical derivatives as a function of step size
Analytic target value f'(x) = 78.0811228185973
Symmetric difference gives decent results (f(x+h)-f(x-h))/2h when 0.0001 < x < 0.001
f' sym corrected when 0.00001 < x < 0.002
Asymmetric difference gives decent results (f(x+h)-f(x))/h when 0.00001 = x (still poor)
f' asy corrected when 0.0001 < x < 0.001
This is a very severe test and the algorithm is struggling to get 4 or 5 significant figures accuracy for the numerical derivative. It is clear that by using the fully corrected forms of the numerical derivative truncated to either second or fourth order it is possible to obtain highly accurate results for difficult functions. A practical implementation would use a short step length initially and double it to enable 3 out 5 function evaluations to be reused. The results are not very sensitive to the exact choice of h as the tabulated difference results show.