Numerical Differentiation

Radovan Omorjan pointed me to this page of the „Kitchin Group“. It is all about numerical stuff done with Python. One trick is very interesting. For numerical differentiation, we can use the simple differential quotient as an estimate. Using the Taylor formula, we can get the error as follows.

\(\dfrac{f(x+h)-f(x)}{h}  = f'(x) + \dfrac{h}{2} f'(\xi)\)

So it is in the order of |h|. The problem, however, is the evaluation of f(x+h)-f(x). If f(x) is not exactly zero, the two numbers have many equal digits and the subtraction cancels these digits leaving much less correct digits. Here is an example computed in EMT.

>format(20,16); ...
>h=1e-8; exp(1+h), exp(1), exp(1+h)-exp(1)
>(exp(1+h)-exp(1))/h, E,

In the result, only 8 digits are correct. If we take a smaller h, then the result is much worse.

>h=1e-14; exp(1+h), exp(1), exp(1+h)-exp(1)
>(exp(1+h)-exp(1))/h, E,

The  optimal h must keep a subtle balance between the accuracy of the formula and the accuracy of the subtraction.

The difference formula can be improved by taking two other points, or even taken more points. The Taylor expansion shows

\(\dfrac{f(x+h)-f(x-h)}{2h} = f'(x) + \dfrac{h^2}{6} f“(\xi)\)

The optimal h for this is the third root of the machine epsilon.

>h=1e-5; (exp(1+h)-exp(1-h))/(2*h), E,

The function diff() of EMT uses four points, which makes things a bit better, but not much.

>diff("exp(x)",1), E,

For analytic functions, there is as simple trick to avoid the cancellation. We take a complex value for h! This is just as good an approximation as the simple differential with real h. Indeed a purely imaginary step yields

\(\text{re}\, \dfrac{f(x+ih)-f(x)}{ih} = \text{im}\, \dfrac{f(x+ih)-f(x)}{h} = \dfrac{\text{im}\,f(x+ih)}{h}\)

since the imaginary part of f(x) is zero. Now, there is no subtraction to be done, and we can use very small h.

>h=1e-16; im(exp(1+I*h))/h, E

The trick is only good for analytic functions, and we assume that we can compute complex values easily. If the function involves complicated real computations we cannot use the trick at all.

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht.

Diese Website verwendet Akismet, um Spam zu reduzieren. Erfahre mehr darüber, wie deine Kommentardaten verarbeitet werden.