How could I improve this approximation?

Question

In a computer application, I need to solve trillions of times an equation which can be reduced to $$f(x)=\sin(x)-a x=0$$ Newton methods (quadratic and higher orders) are used for the solution. Parameter $a$ is random and the solution $x$ looked for is the one between $0$ and $\pi$ if it exists. The retained solution is immediate ($x=0$) if $a \geq 1$ or $a \leq 0$.

For the other cases ($0 \lt a \lt 1$), since I need to save as many iterations as I can, I focused on how to establish a good and unexpensive approximation of the solution. After some empirical experiments, what I found is that writing $$\sin(x) \simeq \frac{4}{\pi^2} x(\pi-x)$$ is a quite good approximation if $0 \leq a \leq 0.7$ leading to $$x \simeq \pi -\frac{\pi ^2 a}{4}$$

For the remaining interval, using Pade approximation $$\sin(x) \simeq \frac{x-\frac{7 x^3}{60}}{1+\frac{x^2}{20}}$$ which leads to $$x \simeq \frac{2 \sqrt{15} \sqrt{1-a}}{\sqrt{3 a+7}}$$ seems interesting.

I wonder if this could be improved, the goal being a simple explicit expression for the estimate of the solution.

Any idea and/or suggestion would really be welcome.

Added later

In comments and answers, Bhoot suggested me to look at the approximation $$\sin(x) \simeq \frac{16 (\pi -x) x}{5 \pi ^2-4 (\pi -x) x}$$ proposed by Mahabhaskariya of Bhaskara I, a seventh-century Indian mathematician. This spendid approximation leads to $$x \simeq \frac{2 \sqrt{-\pi ^2 a^2+2 \pi a+4}+\pi a-4}{2 a}$$ which is effectively very good except very close to $a=1$. This already makes a significant improvement.

Temporary improvement

Extremely impressed by the quality of the approximation proposed by Mahabhaskariya of Bhaskara I (more than 1400 years ago), I tried to undertand why it was very good except in the vicinity of $a=1$. I suspected that the derivative could be in error at the end points. Effectively, this formula gives a slope equal to $\frac{16}{5\pi} \simeq 1.01859$ instead of $1$. On the other side, the area under the curve is given by $$A=\pi \left(-4+\pi +\tan ^{-1}\left(\frac{3116}{237}\right)\right) \simeq 1.99955$$ So, modestly, I built an approximation which is $$\sin(x) \simeq \frac{\pi x(\pi-x)}{\pi^2+(\pi-4)x(\pi-x)}$$ which allows to match exactly the function and derivative values at $x=0,\frac{\pi}{2},\pi$; with respect to accuracy, it is not as good as the original showing a maximum error of $0.0052$ instead of $0.0016$. The area under the curve is then given by $$A=\frac{\pi \left(-4 \pi +\pi ^2+4 \sqrt{(4-\pi ) \pi } \tan ^{-1}\left(\sqrt{\frac{4}{\pi }-1}\right)\right)}{(\pi -4)^2} \simeq 1.99161$$ quite significantly worse than the original.

The estimate of the solution is given by $$x \simeq \frac{\pi \left((\pi -4) a-\sqrt{(\pi -4) a (\pi a-2)+1}+1\right)}{2 (\pi -4) a}$$ which makes Newton scheme converging in less than two iterations for the whole range.

Added after Christian Blatter's answer

I used what has been kindly proposed by Christian Blatter in his answer and set

$$\tilde f^2(a):={p(a)\over q(a)},\qquad p(a):=c_0+c_1 a+ c_2 a^2,\quad q(a):=1+d_1a +d_2 a^2\ $$ Using nonlinear regression, I adjusted the five involved parameters in order to minimize $$SSQ=\sum_{i=1}^n \Big(\sin(\tilde f(a_i))-a_i \tilde f(a_i)\Big)^2$$ The values of the $a_i$ were generated using $1000$ equally spaced values of the $x_i$ between $0$ and $\pi$. I have not been able to compute formally $$\int_0^1 \Big(\sin(\tilde f(a_i))-a_i \tilde f(a_i)\Big)^2 da$$

Starting with the coefficients given in Christian Blatter's answer, the initial $SSQ=5.968\times 10^{-5}$ which is already very good. I arrived to $SSQ= 3.800\times 10^{-6}$. The corresponding parameters are $$c_0=9.86774920$$ $$c_1=4.91765690$$ $$c_2=-14.77935381$$ $$d_1=2.48744104$$ $$d_2=0.63396306$$ For these values, the largest error is $0.000295$ and the average error is $0.000052$ which is incredibly good. As a result, a single Newton iteration is basically required for the desired accuracy. In the following plot the function $\tilde f$ is denoted $g$:

I would like to thank all people who contribute to this work. You have been extremely helpful.

Added later

Continuing working the problem, I set $$\tilde f^2(a):={p(a)\over q(a)},\qquad p(a):=\sum_{i=0}^n c_i a^i,\quad q(a):=1+\sum_{i=1}^n d_i a^i\ $$ and played with $n$. The first result is that moving to cubic polynomials changes a lot the result : $SSQ=4.9609379\times 10^{-10}$, maximum error $= 0.000004$, average error $< 0.000001$. This means that one single Newton iteration is required for a high accuracy. Mowing to fourth oder gives the solution without any Newton iteration.

Answer 1

Hi intersting question usimg hermite polynomial $$\text{Solve}\left[\sum _{n=0}^3 \frac{\left(i (-1)^{2 n} i^n 2^{-n-1} \left((-1)^n-1\right)\right) H_n(x)}{\sqrt[4]{e} n!}=a x,x\right]$$ you could get $$\left\{\{x\to 0\},\left\{x\to -\sqrt{\frac{3}{2}} \sqrt{5-4 \sqrt[4]{e} a}\right\},\left\{x\to \sqrt{\frac{3}{2}} \sqrt{5-4 \sqrt[4]{e} a}\right\}\right\}$$ give arounf three digits of precistion

Answer 2

Bhaskara's sine approximation looks a bit like a Padé approximant and is eerily accurate.

Another option---since you need to do this trillions of times---is to save the solutions for many values of $a$, and then look up/interpolate solutions as you need them.

Answer 3

We have to solve the equation ${\rm sinc}(x)=a$ for $x\in[0,\pi]$ in terms of the parameter $a\in[0,1]$. As Kirill has remarked, for $a\to1\!-\ $ we obtain $x\doteq\sqrt{6(1-a)}$. Therefore, if we want a single simple expression giving accurate values over the whole $a$-interval $[0,1]$ we have to take a square root at the end.

For this reason we put $x:=\sqrt{u}$ and solve the equation $$\bigl(s(u):=\bigr)\quad {\rm sinc}\bigl(\sqrt{u}\bigr)=a\tag{1}$$ for $u$ in terms of $a$. Let $$a\mapsto u:=f(a)\qquad(0\leq a\leq1)$$ be the solution of $(1)$, i.e. the inverse function of $s$. From known values of $s$ and $s'$ we can, e.g., deduce the values $$f(0)=\pi^2,\quad f\left({2\over\pi}\right)={\pi^2\over4},\quad f\left({\sqrt{27}\over 4\pi}\right)={4\pi^2\over9},\quad f(1)=0\ ,\tag{2}$$ and $$f'(0)=-2\pi^2\tag{3}$$ (one could take more such values into account and determine more coefficients $c_k$, $d_k$ in the following).

We now make the "Ansatz" $$\tilde f(a):={p(a)\over q(a)},\qquad p(a):=\pi^2+c_1 a+ c_2 a^2,\quad q(a):=1+d_1a +d_2 a^2\ ,$$ and determine the coefficients $c_1$, $c_2$, $d_1$, $d_2$ such that $(2)$ and $(3)$ are satisfied. This leads to the approximation $$x=g(a):=\sqrt{\tilde f(a)}=\sqrt{{\pi^2+5.95839 a - 15.828 a^2\over 1 + 2.60371 a + 0.690687 a^2}}\ .$$ The following figure shows a plot of $a\mapsto \sin\bigl(g(a)\bigr)-a\> g(a)$:

Answer 4

Here is a way to solve this equation quickly that I think is more efficient that the methods proposed so far. The reason I think this is more efficient is that allows you to get about 15 digits with some multiplications (which are cheap), one square root, and at most one sin-cos evaluation (which are expensive).

Your equation reminds me of the Kepler equation, that I once tried to analyse, and used a similar approach.

Instead of approximating $\sin x$ by some algebraic function, and then solving the approximate equation, let's approximate the true solution instead. Denote the solution $y(a)$, $0<y(a)<\pi$, $0<a<1$, of the equation $$ \sin y(a)-a y(a) = 0, $$ then $y(a)$ has a root, $y(1)=0$, and it behaves at $a\approx1$ as $$ y(a) \sim \sqrt{6}(1-a)^{1/2}. $$ Here is the graph of $y(a)/\sqrt{1-a}$:

Now the red function $y(a)/\sqrt{1-a}$ is smooth, non-singular, has no roots, so we can approximate it with a Chebyshev series $\sum_{k\geq0}c_k T_k(x)-\frac12 c_0$, with these coefficients:

[5.468893253088218, -0.3315308418632746, 0.05711908619285262, -0.0133386391330606, 0.003604590132654065, -0.001061080932244192, 0.0003303619584939474, -0.0001070040549514469, 3.56907424708586e-5, -1.217702854008323e-5, 4.22996639149332e-6, -1.491025149248465e-6, 5.319859653014493e-7, -1.917577069128536e-7, 6.972596591016018e-8, -2.554518680989035e-8, 9.420605435988877e-9, -3.494307402618534e-9, 1.302781037141137e-9, -4.879450499229215e-10, 1.835088955571608e-10, -6.927187870873964e-11, 2.62374324329517e-11, -9.968294055829801e-12, 3.797891128370018e-12, -1.450731732722623e-12, 5.554781896411238e-13, -2.131593799567149e-13, 8.196531167958859e-14, -3.15778010597506e-14, 1.218720260971691e-14, -4.711359031839824e-15, 1.824159314975387e-15, -7.073133344264407e-16, 2.746327946190114e-16, -1.067647609986269e-16, 4.153961510615797e-17, -1.613912600904692e-17, 6.171147771603768e-18, -2.090308692683401e-18]

Depending on how much accuracy you need, it might be enough to sum the first 34 terms using the Clenshaw recurrence (this takes 33 multiplications, and uses no expensive operations, like divisions). Here is the plot of relative error for successive Chebyshev approximations to $y(a)/\sqrt{1-a}$:

Depending on the relative cost of trigonometric operations on your machine (this should be established experimentally), it may be faster to take just enough Chebyshev terms that one Newton iteration suffices to get accuracy to machine precision. Here is 8 terms followed by one Newton iteration (blue) or one Halley iteration (red). The only expensive operation here would the simultaneous evaluation of sin and cos for the function and its derivative:

The cost of this is 7 multiplications for the Chebyshev series, and 4 multiplications, 2 divisions and one sin-cos for the Halley iteration, and one square root and a multiplication to get back $y(a)$.

Answer 5

Back to the problem almost five years later.

Still working it but now requiring less accuracy (looking basically for a short cut method), I considered the Taylor expansion of $\text{sinc}(x)$ around $x=0$ and used sries reversion.

This leads to $$x=t+\frac{t^3}{40}+\frac{107 t^5}{67200}+\frac{3197 t^7}{24192000}+\frac{49513 t^9}{3973939200}+\frac{196309411 t^{11}}{154983628800000}++O\left(t^{13}\right)\tag 1$$ where $t=\sqrt{6(1-a)}$ (already mentioned in Kirill's answer).

The expansion given in $(1)$ can be shortened using Padé approximants and the simplest would be $$x=t \,\frac{1680-65 t^2}{1680-107 t^2}\tag 2$$ $$x=t\,\frac{138524601600-14754907440 t^2+193184137 t^4}{138524601600-18218022480 t^2+428067253 t^4}\tag 3$$ which seem to be quite accurate for my current needs. $$\left( \begin{array}{cccc} a & \text{using }(2) & \text{using }(3) & \text{exact}\\ 0.00 & 3.04416 & 3.12587 & 3.14159 \\ 0.05 & 2.92159 & 2.98242 & 2.99146 \\ 0.10 & 2.80196 & 2.84708 & 2.85234 \\ 0.15 & 2.68478 & 2.71807 & 2.72114 \\ 0.20 & 2.56956 & 2.59395 & 2.59574 \\ 0.25 & 2.45585 & 2.47354 & 2.47458 \\ 0.30 & 2.34316 & 2.35585 & 2.35644 \\ 0.35 & 2.23102 & 2.23999 & 2.24033 \\ 0.40 & 2.11893 & 2.12516 & 2.12535 \\ 0.45 & 2.00634 & 2.01058 & 2.01067 \\ 0.50 & 1.89264 & 1.89544 & 1.89549 \\ 0.55 & 1.77712 & 1.77891 & 1.77894 \\ 0.60 & 1.65892 & 1.66002 & 1.66003 \\ 0.65 & 1.53696 & 1.53761 & 1.53761 \\ 0.70 & 1.40983 & 1.41018 & 1.41019 \\ 0.75 & 1.27552 & 1.27570 & 1.27570 \\ 0.80 & 1.13103 & 1.13110 & 1.13110 \\ 0.85 & 0.97133 & 0.97135 & 0.97135 \\ 0.90 & 0.78668 & 0.78668 & 0.78668 \\ 0.95 & 0.55191 & 0.55191 & 0.55191 \\ 1.00 & 0.00000 & 0.00000 & 0.00000 \end{array} \right)$$

Inspired by $(2)$, quick and dirty regression gives $$x=\frac{87 \sqrt{1-a} (1553 a+728 \pi )}{364 (199 a+174)}$$