<- previous index next ->
Numerical differentiation is typically considered poor at best.
In general, high order techniques are needed to get reasonable
values.
Similar to numerical integration, the more values of the function
that are used, then the more accuracy can be expected.
Similar to numerical integration, given a function that can be
fit accurately by an n th order polynomial, an equation for
the derivative exists to give very good accuracy.
Consider a function f(x) that you can compute but do not know
a symbolic representation. To find the derivatives at a few points,
say 4 for example, and computing only 4 evaluations of f(x):
f(x2)
* f(x3)
| *
f(x1) | |
* | |
| | |
f(x0) | | |
* | | |
| | | |
| | | |
x0 x1 x2 x3 h = x3-x2 = x2-x1 = x1-x0
We use f'(x0) to represent the derivative of f(x) at x0.
f'(x0) = 1/6h (-11 f(x0) +18 f(x1) -9 f(x2) 2 f(x3) )
f'(x1) = 1/6h ( -2 f(x0) -3 f(x1) 6 f(x2) -1 f(x3) )
f'(x2) = 1/6h ( 1 f(x0) -6 f(x1) 3 f(x2) 2 f(x3) )
f'(x3) = 1/6h ( -2 f(x0) +9 f(x1) -18 f(x2) 11 f(x3) )
The error estimate is 1/h^4 f''''(z) for worse z in x0..x3
Thus, if the four points are fit by a third degree polynomial,
f''''(z), the fourth derivative is always zero and the numerical
derivatives are exact within roundoff error.
To check that the coefficients such as -11 18 -9 2 are correct,
we symbolically fit the 4 points with a third order polynomial,
then symbolically take the derivatives. Too much work to do by
hand, thus I used Maple to check and found the coefficients correct.
Using 5 evaluations of f(x) provides the five derivatives:
f'(x0) = 1/12h (-25 f(x0) 48 f(x1) -36 f(x2) 16 f(x3) -3 f(x4) )
f'(x1) = 1/12h ( -3 f(x0) -10 f(x1) 18 f(x2) -6 f(x3) 1 f(x4) )
f'(x2) = 1/12h ( 1 f(x0) -8 f(x1) 0 f(x2) 8 f(x3) -1 f(x4) )
f'(x3) = 1/12h ( -1 f(x0) 6 f(x1) -18 f(x2) 10 f(x3) 3 f(x4) )
f'(x4) = 1/12h ( 3 f(x0) -16 f(x1) 36 f(x2) -48 f(x3) 25 f(x4) )
The above were typed by hand and could have an error.
Some quick checking: If all function evaluations are the same,
then the function is flat and the derivative must be zero at
all points. Thus, the sum of the coefficients must be zero for
every derivative. Note that the first and last coefficients have
the same values, in reverse order with reverse sign.
Always try to cut-and-past from the computer printout:
nderiv.out tabulates a few cases.
These were generated by the computer program:
nderiv.c computes these cases.
nderiv.f90 in Fortran 95
nderiv.adb in Ada 95
Note that the coefficients can be generated in you program by
copying the function deriv and calling it for a specific
derivative order, number of evaluations and point where derivative
is to be computed. For the second derivative using 4 evaluations
and computing the derivative at point zero,
deriv( 2, 4, 0, a, &bb); returns a={2, -5, 4, -1) and bb=1
thus f''(x0) = 1/bb h^2 ( a[0] f(x0) + a[1] f(x1) + a[2] f(x2) + a[3] f(x3) )
or f''(x0) = 1/h^2 (2 f(x0) -5 f(x0+h) +4 f(x0+2h) -1 f(x0+3h) )
Note that integers are returned and should be cast to appropriate
floating point type.
This computation of derivatives will be used extensively in the
lectures that follow on ordinary and partial differential equations.
<- previous index next ->