How do I avoid catastrophic cancellation for small numbers in f(x) = (1-cos(x))/x**2 in Python 3.7?

user_15252134 0
Viewed 664

How do I avoid catastrophic cancellation for small numbers in f(x) = (1-cos(x))/x**2 in Python 3.7?

This is what I tried so far (the key, I know, is some trigonometric identity that enables you to avoid the cancellation, and I also know, having used L'Hopital's rule, that the limit→0 for f(x) is 0.5, so the correct program output is something very close to 0.5, which is what you do get if you use x = 1.2e-4 for example, but you get cancellation with smaller numbers like 1.2e-8, and I need to make it so this doesn't happen).

from math import *
def f(x):     #these are all the same function using different identities  
   a = (1-(sin(x)/tan(x)))/(x**2)
   b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
   c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
   d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
   e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
   return a, b, c, d, e

print(k(1.2e-8))
#Output: (0.0, 0.7709882115452477, 0.0, 0.0, 0.0) - whereas need 0.5000...
floating-accuracy trigonometry python-3.x function python
4 Answers

Like this:

sin(x)/x * tan(x/2)/x

It does the job right up to the end, x = 1e-308 is still OK.

Unfortunately I cannot offer much insight into why it works well.

user_5550450

Use numpy.float128 instead. numpy is the standard for data analysis and more complex math. It can be installed with the following command in the terminal.

pip install numpy

from numpy import *
def f(x):     #these are all the same function using different identities  
   a = (1-(sin(x)/tan(x)))/(x**2)
   b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
   c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
   d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
   e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
   return a, b, c, d, e
print(f(float128(1.2e-8)))

This prints

(0.5003141275115400736, 0.49956120933620291774, 0.49993766842387149567, 0.49993766842387149567, 0.49956120933620291774)
user_109578440

You can conditionally return the limit for small x, like this:

import matplotlib.pyplot as plt
from numpy import *
epsilon=1e-8

def f(x):
   if x<epsilon
      return 0.5

   return (1-sin(x)/tan(x))/x**2
   #note: same as (1-cos(x))/x**2

x=arange(0,6,0.01)
y=vectorize(f)
plt.plot(x,y(x))
plt.show()

The curve plotted looks smooth

Note: I prefer numpy to math. Vectorize makes it possible to call the function with an array (not very efficient, but it is easy to use).

user_7501720

THe problem is the limited precision of float and double. You have to use arithmetic with more precision such as mpfr. It can be used in Python through binding such as

https://pypi.org/project/gmpy2/

Here is an example where I'm using it through a more high level environment called Sagemath: I'm using 100 bits of precision:

sage: R = RealField(100)  // 100 bits of precision
sage: def f(x):     #these are all the same function using different identities  
....:    a = (1-(sin(x)/tan(x)))/(x**2)
....:    b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
....:    c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
....:    d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
....:    e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
....:    return a, b, c, d, e
....: 
sage: f(R(1.2e-8))
(0.50000000000000264647624775223,
 0.49999999999999716827551705076,
 0.49999999999999716827551705076,
 0.49999999999999716827551705076,
 0.49999999999999169007478634929)
user_16921070
Related