The functions f(x)=(1−cos2x)/x2 and g(x)=sin2x/x2 are mathematically indistinguishable, but plotted in the region −0.001≤x≤0.001 show a significant difference. Explain the origin of this difference.
The expression 1 - np.cos(x)**2
suffers from catastrophic cancellation close to x=0
resulting in a dramatic loss of precision and wild oscillations in the plot of f(x) (below). Consider, for example, x = 1.e-9
: in this case, the difference 1 - np.cos(x)**2
is indistinguishable from zero (at double precision) so f(x)
returns 0. Conversely, np.sin(x)**2
is indistinguishable from x**2
and g(x)
returns 1.0
correctly.
import numpy as np
import pylab
f = lambda x: (1 - np.cos(x)**2)/x**2
g = lambda x: (np.sin(x)/x)**2
x = np.linspace(-0.0001, 0.0001, 10000)
pylab.plot(x, f(x), label=r'$(1-\cos^2x)/x^2$')
pylab.plot(x, g(x), label=r'$\sin^2x/x^2$')
pylab.ylim(0.99995, 1.00005)
pylab.legend()
pylab.show()