Решение должно быть точным с точностью до машинного эпсилона

>>> 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
>>> 

Вы можете использовать scipy.optimize для этой задачи:

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

Затем определите свой цель как так:

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

И вы можете решить эту проблему численно:

>>> 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.])

Теперь давайте попробуем с новым значением

>>> 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.])
>>> 

• Во-первых, из границ синуса вы знаете, что любое решение должно быть в интервале [-abs(a),abs(a)]. Если abs(a)\le 1, то единственный корень в [-1,1] - это x=0

• Помимо интервала, содержащего ноль, вы также знаете, что в любом из интервалов есть ровно один корень корни из cos(x)=1/a, которые являются экстремумами из a*sin(x)-x. Установите phi=arccos(1/a) in [0,pi], тогда этими корнями являются -phi+2*k*pi и phi+2*k*pi.

• Интервал для k=0 может содержать 3 корня, если 1<a<0.5*pi. Для положительного корня известно x/a=sin(x)>x-x^3/6, так что x^2>6-6/a.

• И, наконец, проблема симметрична, если x является корнем, то же самое происходит с -x, поэтому все, что вам нужно сделать, это найти положительные корни.

Итак, чтобы вычислить корни,

• Запустите корневой список с корнем 0.
• В случае abs(a)<=1 дальнейших корней нет, возвратите. Можно также использовать -pi/2<=a<=1.
• в случае 1<a<pi/2, применить выбранный метод брекетинга к интервалу [sqrt(6-6/a), pi/2], добавить корень в список и вернуться.

• В остальных случаях, когда abs(a)>=0.5*pi:

  • Вычислить phi=arccos(1/a).
  • Затем для любого натурального числа k примените метод брекетинга к интервалам [2 * (k-1) * pi + phi, 2 * k * pi-phi] и [2 * k * pi- phi, 2 * k * pi-phi, так что (k-0.5)*pi < abs(a) [(k-0.5)*pi, (k+0.5)*pi] до тех пор, пока нижняя граница интервала меньше, чем abs(a) и функция имеет изменение знака в течение интервала.
  • Добавить найденный корень в список. Вернитесь со списком после окончания цикла.

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"/>

Простой способ оценить точность sin() И cos() для данного аргумента x будет выглядеть так:

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

Возможно, вы захотите сбросить последнюю 2, просто чтобы быть на «безопасная сторона» (ну, есть значения x, для которых это приближение не выполняется очень хорошо, в частности для x, близкого к -90 deg). Я бы предложил протестировать около x=pi/4

Объяснение:

Основная идея этого подхода заключается в следующем ... Допустим, наши sin(x) и cos(x) отклоняются от точных значений одно «значение ошибки» eps. То есть exact_sin(x) = sin(x) + eps (то же самое для cos(x)). Кроме того, назовем delta измеренным отклонением от пифагорейской тригонометрической идентичности :

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

Для точных функций delta должно быть равно нулю:

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

или, переходя к неточным функциям:

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

Пренебрежение последним слагаемым 2*eps**2 (допустим небольшие ошибки):

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

Если мы выберем x так, чтобы sin(x)+cos(x) ] колеблется около 1 (или где-то в диапазоне 0.5-2), мы можем приблизительно оценить, что eps = |1 - sin(x)**2 - cos(x)**2|/2.

Точность синусоидальной функции здесь не так важна, вам лучше провести исследование уравнения.

Если вы напишите это в форме sin x / x = sinc x = 1 / a, вы сразу увидите, что число решений - это число пересечений кардинального синуса с горизонталью. Это число зависит от ординат экстремумов последних.

Экстремумы находятся там, где x cos x - sin x = 0 или x = tan x, и соответствующие значения равны cos x. Это снова трансцендентное уравнение, но оно не имеет параметров, и вы можете решить его раз и навсегда. Также отметим, что при увеличении значений x решения становятся все ближе и ближе к (k+1/2)π.

Теперь для заданного значения 1 / a вы можете найти все экстремумы ниже и выше, и это даст вам стартовые интервалы, где искать корни. Секущий метод будет полезен.

Ответ

np.sin в общем случае будет настолько точным, насколько это возможно, учитывая точность переменных double (т.е. 64-битных float), в которых вход, выход и промежуточный значения сохраняются. Вы можете получить разумную меру точности np.sin, сравнив ее с версией произвольной точности sin из 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} 
 Вывод: 
 
     
 
 Таким образом, разумно сказать, что как относительные, так и абсолютные ошибки np.sin имеют верхнюю границу 2e-16. 
 
 Лучший ответ 
 
 Существует отличный шанс, что если вы сделаете increment достаточно маленьким, чтобы ваш подход был точным, ваш алгоритм будет слишком медленным для практического использования. Подходы к решению стандартных уравнений не будут работать для вас, поскольку у вас нет стандартной функции. Вместо этого у вас есть неявная многозначная функция. Вот пример общего подхода для получения всех решений этого вида уравнения: 
 
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 
 Вывод: 
 
    [1119]  
 
 Вам придется настроить initial_guess в зависимости от вашего значения a. initial_guess должно быть не меньше фактического количества решений. 
)
plt.plot(x, 10*np.sin(x), label=r'$y = 10 \cdot sin(x) 
 Вывод: 
 
    [1119]  
 
 Вам придется настроить initial_guess в зависимости от вашего значения a. initial_guess должно быть не меньше фактического количества решений. 
)
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()
 
 Вывод: 
 
    [1119]  
 
 Вам придется настроить initial_guess в зависимости от вашего значения a. initial_guess должно быть не меньше фактического количества решений. 
)
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()

Вывод:

Таким образом, разумно сказать, что как относительные, так и абсолютные ошибки np.sin имеют верхнюю границу 2e-16.

Лучший ответ

Существует отличный шанс, что если вы сделаете increment достаточно маленьким, чтобы ваш подход был точным, ваш алгоритм будет слишком медленным для практического использования. Подходы к решению стандартных уравнений не будут работать для вас, поскольку у вас нет стандартной функции. Вместо этого у вас есть неявная многозначная функция. Вот пример общего подхода для получения всех решений этого вида уравнения:

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 
 Вывод: 
 
    [1119]  
 
 Вам придется настроить initial_guess в зависимости от вашего значения a. initial_guess должно быть не меньше фактического количества решений. 
)
plt.plot(x, 10*np.sin(x), label=r'$y = 10 \cdot sin(x) 
 Вывод: 
 
    [1119]  
 
 Вам придется настроить initial_guess в зависимости от вашего значения a. initial_guess должно быть не меньше фактического количества решений. 
)
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()

Вывод:

[1119]

Вам придется настроить initial_guess в зависимости от вашего значения a. initial_guess должно быть не меньше фактического количества решений.