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

enter image description here

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.

  • 4
    $\begingroup$ Maybe of interest: en.wikipedia.org/wiki/… $\endgroup$ – user105475 Jul 19 '14 at 7:03
  • $\begingroup$ @Bhoot. This is a superb approximation (I did not know it). Please, put it as an answer to my post. I shall edit and refer to you. $\endgroup$ – Claude Leibovici Jul 19 '14 at 7:20
  • $\begingroup$ @ClaudeLeibovici: You are right! It is a nice approximation formula. $\endgroup$ – Mhenni Benghorbal Jul 19 '14 at 9:22
  • 1
    $\begingroup$ @MhenniBenghorbal. When you think that is was done almost $1400$ years ago ! Cheers :) $\endgroup$ – Claude Leibovici Jul 19 '14 at 9:24
  • $\begingroup$ @ClaudeLeibovici: He is a seventh-century Indian mathematician. $\endgroup$ – Mhenni Benghorbal Jul 19 '14 at 9:25

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

  • $\begingroup$ This looks very nice ! Where did you find this expansion of $\sin(x)$ ? I never saw it before. Thanks. $\endgroup$ – Claude Leibovici May 31 '15 at 3:32

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.

  • $\begingroup$ Table lookup and interpolation are quite expensive. Again, thank you very much for this superb reference. Cheers :) $\endgroup$ – Claude Leibovici Jul 19 '14 at 7:33
  • $\begingroup$ You're welcome. I think the lookup can be made to be quite cheap, since you know the structure of the table. But yes, for trillions of calls, even a small cost per call can mean a lot. An interesting problem! $\endgroup$ – user105475 Jul 19 '14 at 7:45
  • $\begingroup$ If you are interested, have a look to what I added. Thanks to you, I have been able to make just incredible improvements in few hours. Thanks again. Cheers :) $\endgroup$ – Claude Leibovici Jul 19 '14 at 15:33

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)$:

enter image description here

  • $\begingroup$ I thank you very much for this very interesting answer. I shall start working on this basis right now. Thanks again. Cheers :-) $\endgroup$ – Claude Leibovici Sep 27 '14 at 5:34
  • $\begingroup$ I worked based on your ideas and I made a quite good approximation. I added that to my question and let it for your review (if you wish). I would have liked to add a plot similar to your; the problem is that I am blind and producing graphics is very very hard to me. May be, if you wish, could you add one ? Thanks again. $\endgroup$ – Claude Leibovici Sep 29 '14 at 4:59
  • $\begingroup$ My wife just tells me that the plot is very nice and illustrative. Thanks again. Cheers :) $\endgroup$ – Claude Leibovici Sep 29 '14 at 10:22
  • $\begingroup$ May be you should be amazed by the following :increasing the degree of the polynomials makes the solution better and better, not to say almost exact. The most funny is that this can be done without requiring initial estimates ! Thanks again. $\endgroup$ – Claude Leibovici Oct 1 '14 at 9:13

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}$:

enter image description here

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}$:

enter image description here

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:

enter image description here

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

  • $\begingroup$ This is an very intersting answer and I thank you very much. I shall start working on this basis right now. Moreover, from a purely esthetical point of view, it is beautiful ! Thanks again. Cheers :-) $\endgroup$ – Claude Leibovici Sep 27 '14 at 5:30

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)}$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.