Mejor algoritmo para encontrar el valor numérico de uno sobre la raíz cuadrada

Question

Descripción: Para hallar el valor de $1/\sqrt{a}$ básicamente encontramos la raíz $x^{\star}$ de la función $f(x)$

$$ f(x) = x^2 - \dfrac{1}{a} $$

Una de las formas de hallarlo es utilizando el método de Newton e iterando a partir de un valor inicial $x_0$ .

$$ x_{k+1} = x_{k} - \dfrac{f(x_{k})}{f'(x_k)} = \dfrac{1}{2}\left(x_k + \dfrac{1}{ax_k}\right) $$ $$ \lim_{k\to \infty}x_{k} = \dfrac{1}{\sqrt{a}} $$

Esto se conoce como método de Newton que tiene convergencia cuadrática, para alguna constante positiva $M$ :

$$ \left|x_{k+1} - \dfrac{1}{\sqrt{a}}\right| \le M \cdot \left|x_{k} - \dfrac{1}{\sqrt{a}}\right|^2 $$

Pregunta: ¿Existe un método numérico mejor, cuya convergencia sea superior a la del método de Newton, para hallar $x^{\star} = \dfrac{1}{\sqrt{a}}$ ?

EDITAR: He pensado en utilizar la expansión de Taylor hasta segundo orden, una vez que el método de Newton se realiza utilizando la expansión de primer orden:

$$ f(x) = f(x_k) + (x-x_k) \cdot f'(x_k) + \dfrac{(x-x_k)^2}{2} \cdot f''(x_k) + \dfrac{(x-x_k)^3}{6} \cdot f'''(\varepsilon(x)) $$

Y luego encontrar el valor de $x_{k+1}$ que

$$0 = f(x_k) + (x_{k+1} - x_k) \cdot f'(x_k)+ \dfrac{(x_{k+1}-x_k)^2}{2} \cdot f''(x_k)$$

Mediante la función $f(x) = x - \dfrac{1}{ax}$ reescribimos

$$x_{k+1}^2 - \left(3+ax_k^2\right) \cdot x_k \cdot x_{k+1} + 3x_k^2 = 0$$

¿Qué solución es

$$x_{k+1} = \dfrac{x_k}{2}\left(3ax_k^2 - \sqrt{\left(3+ax_{k}^2\right)^2 - 12}\right)$$

Como hay una raíz cuadrada, puede ser difícil de calcular y entonces hacemos una aproximación de Taylor de segundo orden alrededor de $\frac{1}{\sqrt{a}}$ de esta raíz cuadrada:

$$\sqrt{\left(3+ax^2\right)^2 - 12}\approx 6x\sqrt{a} - ax^2 -3$$

Una vez $\sqrt{a} \approx \dfrac{1}{2}\left(y+\dfrac{a}{y}\right) = \dfrac{1}{2}\left(\dfrac{1}{x} + ax\right)$ entonces

$$\sqrt{\left(3+ax^2\right)^2 - 12}\approx 6x \cdot \dfrac{1}{2}\left(\dfrac{1}{x} + ax\right) -ax^2 - 3 = 2ax^2$$ $$\boxed{x_{k+1} = \dfrac{x_k(3-ax_k^2)}{2}}$$

Que es una función simple de polinomios.

Pero no tengo ni idea de la convergencia de utilizar esta función de iteración. Su orden sería $3$ porque utilizamos la expansión de taylor de segundo orden de $f$ pero utilizando aproximaciones lineales de la raíz cuadrada alrededor de $1/\sqrt{a}$ parece que es abajo $3$ .

1. Entonces, ¿qué hay de otras funciones como $f=x^4-a^{-2}$ y utilizando los mismos trucos?

Motivación: Existe un algoritmo conocido como Raíz cuadrada inversa rápida (más detalles aquí ) que parte de una buena estimación (utilizando el logaritmo y la manipulación de bits), pero utiliza la iteración de Newton para refinar el valor en el paso final.

Answer 1

Método de orden de Householder $d$ resolver $f(x) = 0$ da

$$x \to H_d(x) = x + d \frac {\left(1/f\right)^{(d-1)}(x)}{\left(1/f\right)^{(d)}(x)}$$

con tasa de convergencia $d + 1$ .

Utilizando el sistema de álgebra computacional (wx)Maxima, para

$$f(x) = x^{-2} - a$$

Recibo

$$\begin{aligned} H_1(x) &= \frac{(3 - a x^2)x}{2} \\ H_2(x) &= \frac{(3 + a x^2)x}{3 a x^2 + 1} \\ H_3(x) &= \frac{a^2x^4 + 6ax^2 + 1}{4ax(ax^2+1)} \\ H_4(x) &= \frac{x(a^2x^4 + 10ax^2 + 5)}{5a^2x^4 + 10ax^2 + 1} \\ \end{aligned}$$

Supongamos que la precisión objetivo es $N$ dígitos exactos. La última iteración que da $N$ necesidades de dígitos exactos $N/(d+1)$ dígitos precisos como entrada, por lo que hay alrededor de $\log_{d+1} N$ pasos de iteración necesarios.

Supongamos que el coste de $D$ -La multiplicación de dígitos (y por tanto la división) viene dada por alguna función $M(D)$ (por ejemplo $M(D) = D^{\log_2 3}$ cuando se utiliza el algoritmo de Karatsuba), y suponer que la multiplicación y la división por enteros pequeños, y la suma y la resta de $D$ -son insignificantes en comparación. Entonces el coste depende de si el número de dígitos de $a$ es pequeño o grande. Muchos cálculos pueden reutilizarse/compartirse para reducir el número de multiplicaciones.

Para los pequeños $a$ :

$$\begin{aligned} \operatorname{cost}(H_1) = 2 (M(N/1) + M(N/2) + M(N/4) + \cdots + M(1)) &< \frac{2}{\log 2} M(N) \log(N) \\ \operatorname{cost}(H_2) = 3 (M(N/1) + M(N/3) + M(N/9) + \cdots + M(1)) &< \frac{3}{\log 3} M(N)\log(N) \\ \operatorname{cost}(H_3) = 4 (M(N/1) + M(N/4) + M(N/16) + \cdots + M(1)) &< \frac{4}{\log 4} M(N) \log(N) \\ \operatorname{cost}(H_4) = 4 (M(N/1) + M(N/5) + M(N/25) + \cdots + M(1)) &< \frac{4}{\log 5} M(N) \log(N) \\ \end{aligned}$$

que los pedidos más elevados mejoran el coste global, al menos cuando $N$ es lo suficientemente grande como para que las adiciones, etc., sean realmente insignificantes.

Para $a$ con $N$ dígitos:

$$\begin{aligned} \operatorname{cost}(H_1) = 3 (M(N/1) + M(N/2) + M(N/4) + \cdots + M(1)) &< \frac{3}{\log 2} M(N) \log(N) \\ \operatorname{cost}(H_2) = 4 (M(N/1) + M(N/3) + M(N/9) + \cdots + M(1)) &< \frac{4}{\log 3} M(N) \log(N) \\ \operatorname{cost}(H_3) = 5 (M(N/1) + M(N/4) + M(N/16) + \cdots + M(1)) &< \frac{5}{\log 4} M(N) \log(N) \\ \operatorname{cost}(H_4) = 5 (M(N/1) + M(N/5) + M(N/25) + \cdots + M(1)) &< \frac{5}{\log 5} M(N) \log(N) \\ \end{aligned}$$

Referencia: https://en.wikipedia.org/wiki/Householder%27s_method

Código fuente de (wx)Maxima:

f(x) := x^(-2) - a;
H(d, x) := factor(ratsimp(x + d * (diff(1/f(x), x, d - 1)) / (diff(1/f(x), x, d))));
H(1, x);
H(2, x);
H(3, x);