Suppose you’d like to evaluate the function

for small values of *z*, say *z* = 10^{−8}. This example comes from [1].

The Python code

from numpy import exp def f(z): return (exp(z) - 1 - z)/z**2 print(f(1e-8))

prints -0.607747099184471.

Now suppose you suspect numerical difficulties and compute your result to 50 decimal places using `bc -l`

.

scale = 50 z = 10^-8 (e(z) - 1 - z)/z^2

Now you get *u*(*z*) = .50000000166666667….

This suggests original calculation was completely wrong. What’s going on?

For small *z*,

*e*^{z} ≈ 1 + *z*

and so we lose precision when directly evaluating the numerator in the definition of *u*. In our example, we lost *all* precision.

The mean value theorem from complex analysis says that the value of an analytic function at a point equals the average of the values over a circle centered at that point. If we take the average of 16 values in a circle of radius 1 around our point, we get full accuracy. The Python code

def g(z): N = 16 ws = z + exp(2j*pi*arange(N)/N) return sum(f(ws))/N

returns

0.5000000016666668 + 8.673617379884035e-19j

which departs from the result calculated with `bc`

in the 16th decimal place.

At a high level, we’re avoiding numerical difficulties by averaging over points far from the difficult region.

[1] Lloyd N. Trefethen and J. A. C. Weideman. The Exponentially Convergent Trapezoid Rule. SIAM Review. Vol. 56, No. 3. pp. 385–458.

The post Numerical application of mean value theorem first appeared on John D. Cook.