How precise is numpy's sin(x) ? How do I find out? [need it to numerically solve x=a*sin(x)]

Viewed 3147

I'm trying to numerically solve the equation x=a*sin(x), where a is some constant, in python. I already tried first solving the equation symbolically, but it seems this particular shape of expression isn't implemented in sympy. I also tried using sympy.nsolve(), but it only gives me the first solution it encounters.

My plan looks something like this:

x=0
a=1
rje=[]
while(x<a):
    if (x-numpy.sin(x))<=error_sin:
        rje.append(x)
    x+=increment

print(rje)

I don't want to waste time or risk missing solutions, so I want to know how to find out how precise numpy's sinus is on my device (that would become error_sin).

edit: I tried making both error_sin and increment equal to the machine epsilon of my device but it a) takes to much time, and b) sin(x) is less precise that x and so I get a lot of non-solutions (or rather repeated solutions because sin(x) grows much slower than x). Hence the question.

edit2: Could you please just help me answer the question about precision of numpy.sin(x)? I provided information about the purpose purely for context.

5 Answers

The answer

np.sin will in general be as precise as possible, given the precision of the double (ie 64-bit float) variables in which the input, output, and intermediate values are stored. You can get a reasonable measure of the precision of np.sin by comparing it to the arbitrary precision version of sin from mpmath:

import matplotlib.pyplot as plt
import mpmath
from mpmath import mp

# set mpmath to an extremely high precision
mp.dps = 100
x = np.linspace(-np.pi, np.pi, num=int(1e3))

# numpy sine values
y = np.sin(x)

# extremely high precision sine values
realy = np.array([mpmath.sin(a) for a in x])

# the end results are arrays of arbitrary precision mpf values (ie abserr.dtype=='O')
diff = realy - y
abserr = np.abs(diff)
relerr = np.abs(diff/realy)

plt.plot(x, abserr, lw=.5, label='Absolute error')
plt.plot(x, relerr, lw=.5, label='Relative error')
plt.axhline(2e-16, c='k', ls='--', lw=.5, label=r'$2 \cdot 10^{-16}$')
plt.yscale('log')
plt.xlim(-np.pi, np.pi)
plt.ylim(1e-20, 1e-15)
plt.xlabel('x')
plt.ylabel('Error in np.sin(x)')
plt.legend()

Output:

enter image description here

Thus, it is reasonable to say that both the relative and absolute errors of np.sin have an upper bound of 2e-16.

A better answer

There's an excellent chance that if you make increment small enough for your approach to be accurate, your algorithm will be too slow for practical use. The standard equation solving approaches won't work for you, since you don't have a standard function. Instead, you have an implicit, multi-valued function. Here's a stab at a general purpose approach for getting all solutions to this kind of equation:

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as spo

eps = 1e-4

def func(x, a):
    return a*np.sin(x) - x

def uniqueflt(arr):
    b = arr.copy()
    b.sort()
    d = np.append(True, np.diff(b))
    return b[d>eps]

initial_guess = np.arange(-9, 9) + eps
# uniqueflt removes any repeated roots
roots = uniqueflt(spo.fsolve(func, initial_guess, args=(10,)))
# roots is an array with the 7 unique roots of 10*np.sin(x) - x == 0:
#     array([-8.42320394e+00, -7.06817437e+00, -2.85234190e+00, -8.13413225e-09,
#            2.85234189e+00,  7.06817436e+00,  8.42320394e+00])

x = np.linspace(-20, 20, num=int(1e3))
plt.plot(x, x, label=r'$y = x$')
plt.plot(x, 10*np.sin(x), label=r'$y = 10 \cdot sin(x)$')
plt.plot(roots, 10*np.sin(roots), '.', c='k', ms=7, label='Solutions')

plt.ylim(-10.5, 20)
plt.gca().set_aspect('equal', adjustable='box')
plt.legend()

Output:

enter image description here

You'll have to tweak the initial_guess depending on your value of a. initial_guess has to be at least as large as the actual number of solutions.

The accuracy of the sine function is not so relevant here, you'd better perform the study of the equation.

If you write it in the form sin x / x = sinc x = 1 / a, you immediately see that the number of solutions is the number of intersections of the cardinal sine with an horizontal. This number depends on the ordinates of the extrema of the latter.

The extrema are found where x cos x - sin x = 0 or x = tan x, and the corresponding values are cos x. This is again a transcendental equation, but it is parameterless and you can solve it once for all. Also note that for increasing values of x, the solutions get closer and closer to (k+1/2)π.

Now for a given value of 1 / a, you can find all the extrema below and above and this will give you starting intervals where to look for the roots. The secant method will be handy.

enter image description here

A simple way to estimate the accuracy of sin() AND cos() for a given argument x would be:

eps_trig = np.abs(1 - (np.sin(x)**2 + np.cos(x)**2)) / 2

You may want to drop last 2 just to be on the "safe side" (well, there are values of x for which this approximation does not hold very well, in particular for x close to -90 deg). I would suggest to test at around x=pi/4

Explanation:

The basic idea behind this approach is as follows... Let's say our sin(x) and cos(x) deviates from exact values by a single "error value" eps. That is, exact_sin(x) = sin(x) + eps (same for cos(x)). Also, let's call delta to be the measured deviation from the Pythagorean trigonometric identity:

delta = 1 - sin(x)**2 - cos(x)**2

For exact functions, delta should be zero:

1 - exact_sin(x)**2 - exact_cos(x)**2 == 0

or, going to inexact functions:

1 - (sin(x) + eps)**2 - (cos(x) + eps)**2 == 0 =>
1 - sin(x)**2 - cos(x)**2 = delta = 2*eps*(sin(x) + cos(x)) + 2*eps**2

Neglecting last term 2*eps**2 (assume small errors):

2*eps*(sin(x)+cos(x)) = 1 - sin(x)**2 - cos(x)**2

If we choose x such that sin(x)+cos(x) hovers around 1 (or, somewhere in the range 0.5-2), we can roughly estimate that eps = |1 - sin(x)**2 - cos(x)**2|/2.

To the precision you already got good answers. To the task itself, you can be faster by investing some calculus.

  • First, from the bounds of the sine you know that any solution must be in the interval [-abs(a),abs(a)]. If abs(a)\le 1 then the only root in [-1,1] is x=0

  • Apart from the interval containing zero, you also know that there is exactly one root in any of the intervals between the roots of cos(x)=1/a which are the extrema of a*sin(x)-x. Set phi=arccos(1/a) in [0,pi], then these roots are -phi+2*k*pi and phi+2*k*pi.

  • The interval for k=0 might contain 3 roots if 1<a<0.5*pi. For the positive root one knows x/a=sin(x)>x-x^3/6 so that x^2>6-6/a.

  • And lastly, the problem is symmetric, if x is a root, so is -x so all you have to do is find the positive roots.

So to compute the roots,

  • Start the root list with the root 0.
  • in the case abs(a)<=1, there are no further roots, return. One could also use -pi/2<=a<=1.
  • in the case 1<a<pi/2, apply the chosen bracketing method to the interval [sqrt(6-6/a), pi/2], add the root to the list, and return.
  • In the remaining cases where abs(a)>=0.5*pi:

    • Compute phi=arccos(1/a).
    • Then for any positive integer k apply the bracketing method to the intervals [2*(k-1)*pi+phi,2*k*pi-phi] and [2*k*pi-phi,2*k*pi-phi so that (k-0.5)*pi < abs(a) [(k-0.5)*pi, (k+0.5)*pi] as long as the lower interval boundary is smaller than abs(a) and the function has a sign change over the interval.
    • Add the root found to the list. Return with the list after the loop ends.

let a=10;

function f(x) { return x - a * Math.sin(x); }

findRoots();

//-------------------------------------------------

function findRoots() {
    log.innerHTML = `<p>roots for parameter a=${a}`;
    rootList.innerHTML = "<tr><th>root <i>x</i></th><th><i>x-a*sin(x)</i></th><th>numSteps</th></tr>";
    rootList.innerHTML += "<tr><td>0.0<td>0.0<td>0</tr>";

    if( Math.abs(a)<=1) return;

    if( (1.0<a) && (a < 0.5*Math.PI) ) {
        illinois(Math.sqrt(6-6/a), 0.5*Math.PI);
        return;
    }

    const phi = Math.acos(1.0/a);
    log.innerHTML += `phi=${phi}<br>`;
    let right = 2*Math.PI-phi;
    for (let k=1; right<Math.abs(a); k++) { 
        let left = right;
        right = (k+2)*Math.PI + ((0==k%2)?(-phi):(phi-Math.PI));
        illinois(left, right);
    }
}

function illinois(a, b) {
  log.innerHTML += `<p>regula falsi variant illinois called for interval [a,b]=[${a}, ${b}]`;
  let fa = f(a);
  let fb = f(b);
  let numSteps=2;
  log.innerHTML += ` values f(a)=${fa}, f(b)=${fb}</p>`;

  if (fa*fb > 0) return;

  if (Math.abs(fa) < Math.abs(fb)) { var h=a; a=b; b=h; h=fa; fa=fb; fb=h;}
  
  while(Math.abs(b-a) > 1e-15*Math.abs(b)) {
    let c = b - fb*(b-a)/(fb-fa);
    let fc = f(c); numSteps++;
    log.innerHTML += `step ${numSteps}: midpoint c=${c}, f(c)=${fc}<br>`;

    if ( fa*fc < 0 ) {
      fa *= 0.5;
    } else {
      a = b; fa = fb;
    }
    b = c; fb = fc;
  }
  rootList.innerHTML += `<tr><td>${b}<td>${fb}<td>${numSteps}</tr>`;
}

aInput.addEventListener('change', () => {
  let a_new = Number.parseFloat(aInput.value);
  if( isNaN(a_new) ) {
      alert('Not a number '+aInput.value);
  } else if(a!=a_new) {
     a = a_new;
     findRoots();
  }
});
<p>Factor <i>a</i>: <input id="aInput" value="10" /></p>
<h3>Root list</h3>
<table id="rootList" border = 1>
</table>
<h3>Computation log</h3>
<div id="log"/>

The solution should be precise up to machine epsilon

>>> from numpy import sin as sin_np
>>> from math import sin as sin_math
>>> x = 0.0
>>> sin_np(x) - x
0.0
>>> sin_math(x) - x
0.0
>>> 

You could consider using scipy.optimize for this problem:

>>> from scipy.optimize import minimize
>>> from math import sin
>>> a = 1.0

Then define your objective as so:

>>> def obj(x):
...     return abs(x - a*sin(x))
... 

And you can go ahead and solve this problem numerically by:

>>> sol = minimize(obj, 0.0)
>>> sol
      fun: array([ 0.])
 hess_inv: array([[1]])
      jac: array([ 0.])
  message: 'Optimization terminated successfully.'
     nfev: 3
      nit: 0
     njev: 1
   status: 0
  success: True
        x: array([ 0.])

Now lets try with a new value of a

>>> a = .5
>>> sol = minimize(obj, 0.0)
>>> sol
      fun: array([ 0.])
 hess_inv: array([[1]])
      jac: array([ 0.5])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 315
      nit: 0
     njev: 101
   status: 2
  success: False
        x: array([ 0.])
>>> 

In case you want to find a non-trivial solution to this problem, you need to change x0 iteratively to values greater than zero and also lesser than. Also, manage the bounds of x in minimise by setting bounds in scipy.optimize.minimize, you would be able to walk from -infty to +infty ( or very large numbers ).

Related