¿Cómo podría mejorar esta aproximación?

Question

En una aplicación informática, necesito solucionar miles de millones de veces una ecuación que puede ser reducido a $$f(x)=\sin(x)-a x=0$$ Newton methods (quadratic and higher orders) are used for the solution. Parameter $un$ 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 $\geq 1$ or $\leq 0$.

Para el resto de casos ($0 \lt a \lt 1$), ya que la necesito para guardar tantas iteraciones como puedo, me centré en cómo establecer una buena y unexpensive aproximación de la solución. Después de algunos experimentos empíricos, lo que he encontrado es que escribir $$\sin(x) \simeq \frac{4}{\pi^2} x(\pi-x)$$ is a quite good approximation if $0 \leq \leq 0.7$ leading to $$x \simeq \pi -\frac{\pi ^2 a}{4}$$

Para el resto de intervalo, usando la aproximación de Pade $$\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}}$$ parece muy interesante.

Me pregunto si esto podría ser mejorado, el objetivo de ser un simple expresión explícita para la estimación de la solución.

Cualquier idea y/o sugerencia, sería muy bienvenida.

Añadió más tarde

En los comentarios y respuestas, Bhoot me sugirió buscar en la aproximación $$\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$. Esto ya se hace una mejora significativa.

Mejora temporal

Muy impresionado por la calidad de la aproximación propuesta por Mahabhaskariya de Bhaskara I (más de 1400 años atrás), traté de entender por qué era muy buena, excepto en las proximidades de $a=1$. Yo sospechaba que la derivada podría estar en un error en los puntos finales. Efectivamente, esta fórmula se obtiene una pendiente igual a $\frac{16}{5\pi} \simeq 1.01859$ en lugar de $1$. Por otro lado, el área bajo la curva está dada por $$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 } \bronceado ^{-1}\left(\sqrt{\frac{4}{\pi }-1}\right)\right)}{(\pi -4)^2} \simeq 1.99161$$ muy significativamente peor que el original.

La estimación de la solución está dada por $$x \simeq \frac{\pi \left((\pi -4) a-\sqrt{(\pi -4) a (\pi a-2)+1}+1\right)}{2 (\pi -4) a}$$ lo que hace Newton esquema de convergencia en menos de dos iteraciones para toda la gama.

Añadido después de que Christian la respuesta de Blatter

Yo lo que ha sido amablemente propuesto por Christian Blatter en su respuesta y establecer

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

Comenzando con los coeficientes dados en Cristiano la respuesta de Blatter, la inicial $SSQ=5.968\times 10^{-5}$ que ya es muy bueno. Llegué a $SSQ= 3.800\times 10^{-6}$. Los parámetros correspondientes son $$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$:

Me gustaría agradecer a todas las personas que contribuyeron a este trabajo. Ustedes han sido de gran ayuda.

Añadió más tarde

Continuar trabajando el problema, me puse $$\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$. Esto significa que una sola iteración de Newton se requiere para una alta precisión. La siega de la cuarta oder da la solución sin ningún tipo de iteración de Newton.

Answer 1

Tenemos que resolver la ecuación de ${\rm sinc}(x)=a$ $x\in[0,\pi]$ en términos del parámetro de $a\in[0,1]$. Como Kirill ha comentado, para $a\to1\!-\ $ obtenemos $x\doteq\sqrt{6(1-a)}$. Por lo tanto, si queremos que una simple expresión dando valores precisos sobre la totalidad de la $a$-intervalo de $[0,1]$ tenemos que tomar una raíz cuadrada en el extremo.

Por esta razón ponemos a $x:=\sqrt{u}$ y resolver la ecuación de $$\bigl(s(u):=\bigr)\quad {\rm sinc}\bigl(\sqrt{u}\bigr)=a\tag{1}$$ para $u$ en términos de $a$. Deje $$a\mapsto u:=f(a)\qquad(0\leq a\leq1)$$ ser la solución de $(1)$, es decir, la función inversa de la $s$. A partir de los valores conocidos de $s$ $s'$ podemos, por ejemplo, deducir los valores $$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}$$ y $$f'(0)=-2\pi^2\tag{3}$$ (uno podría tomar más en cuenta los valores y de determinar los coeficientes de $c_k$, $d_k$ en el siguiente).

Ahora hacemos el "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\ ,$$ y determinar los coeficientes $c_1$, $c_2$, $d_1$, $d_2$ tal que $(2)$ $(3)$ está satisfecho. Esto conduce a la aproximación $$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}}\ .$$ La siguiente figura muestra un gráfico de $a\mapsto \sin\bigl(g(a)\bigr)-a\> g(a)$:

Answer 2

Hola interesante pregunta usimg hermite % polinomio $$\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]$$ $$\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\}$ $ que usted podría conseguir dar surge tres cifras de precistion

Answer 3

Aquí es una manera de resolver esta ecuación rápidamente que creo que es más eficiente que los métodos propuestos hasta el momento. La razón por la que creo esto es más eficiente es la que permite obtener alrededor de 15 dígitos con algunos multiplicaciones (que son baratos), una raíz cuadrada, y en la mayoría de uno el pecado-cos de evaluación (que son costosos).

Su ecuación me recuerda a la ecuación de Kepler, que una vez traté de analizar, y se utiliza un enfoque similar.

En lugar de aproximar $\sin x$ por algunos algebraica de la función y, a continuación, solución aproximada de la ecuación, vamos a aproximar la solución verdadera en su lugar. Denotar la solución $y(a)$, $0<y(a)<\pi$, $0<a<1$, de la ecuación $$ \sin y(a)-a y(a) = 0, $$ a continuación, $y(a)$ tiene una raíz, $y(1)=0$, y se comporta en $a\approx1$ $$ y(a) \sim \sqrt{6}(1-a)^{1/2}. $$ Aquí está el gráfico de $y(a)/\sqrt{1-a}$:

Ahora el rojo de la función $y(a)/\sqrt{1-a}$ es suave, no-singular, no tiene las raíces, por lo que podemos aproximar con una serie de Chebyshev $\sum_{k\geq0}c_k T_k(x)-\frac12 c_0$, con estos coeficientes:

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

Dependiendo de cuánta precisión que necesita, podría ser suficiente para sumar el 34 primeros términos de uso de la Clenshaw recurrencia (esto lleva 33 multiplicaciones, y los usos no caro operaciones, como divisiones). Aquí es la trama de error relativo para las sucesivas Chebyshev aproximaciones a $y(a)/\sqrt{1-a}$:

Dependiendo del costo relativo de las operaciones trigonométricas en su de la máquina (esto debe ser establecido experimentalmente), puede ser más rápido tomar sólo lo suficiente de Chebyshev términos que una iteración de Newton es suficiente para obtener precisión precisión de la máquina. Aquí es de 8 términos seguido por uno de Newton iteración (azul) o un Halley iteración (rojo). La única caro la operación aquí sería la evaluación simultánea del pecado y de la cos para la la función y su derivada:

El costo de este es de 7 multiplicaciones por el de Chebyshev de la serie, y 4 multiplicaciones, 2 divisiones y un pecado-cos para el Halley iteración, y una raíz cuadrada y una multiplicación para volver $y(a)$.