Método de Newton para hallar raíces

Este es un método iterativo inventado por Isaac Newton alrededor de 1664. Sin embargo, este método a veces también se llama método de Raphson, porque Raphson inventó el mismo algoritmo unos años después de Newton, pero su artículo se publicó mucho antes.

La tarea es la siguiente. Se da la siguiente ecuación:

$f(x) = 0$

Queremos resolver la ecuación. Más precisamente, queremos hallar una de sus raíces (se asume que la raíz existe). Se asume que $f(x)$ es continua y diferenciable en un intervalo $[a, b]$ .

Algoritmo

Los parámetros de entrada del algoritmo consisten no solo en la función $f(x)$ sino también en la aproximación inicial: algún $x_0$ , con el que empieza el algoritmo.

plot_f(x)

Supongamos que ya hemos calculado $x_i$ ; calculamos $x_{i+1}$ como sigue. Trazamos la tangente al gráfico de la función $f(x)$ en el punto $x = x_i$ , y hallamos el punto de intersección de esta tangente con el eje $x$ . $x_{i+1}$ se pone igual a la coordenada $x$ del punto hallado, y repetimos todo el proceso desde el principio.

No es difícil obtener la siguiente fórmula,

$x_{i+1} = x_i - \frac{f(x_i)}{f^\prime(x_i)}$

Primero, calculamos la pendiente $f’(x)$ , derivada de $f(x)$ , y luego determinamos la ecuación de la tangente, que es

$y - f(x_i) = f’(x_i)(x - x_i)$

La tangente se interseca con el eje x en la coordenada $y = 0$ y $x = x_{i+1}$ ,

$- f(x_i) = f’(x_i)(x_{i+1} - x_i)$

Ahora, resolviendo la ecuación obtenemos el valor de $x_{i+1}$ .

Intuitivamente está claro que si la función $f(x)$ es “buena” (suave), y $x_i$ está lo bastante cerca de la raíz, entonces $x_{i+1}$ estará aún más cerca de la raíz deseada.

La velocidad de convergencia es cuadrática, lo que, hablando de forma condicional, significa que el número de dígitos exactos en el valor aproximado $x_i$ se duplica con cada iteración.

Aplicación para calcular la raíz cuadrada

Usemos el cálculo de la raíz cuadrada como ejemplo del método de Newton.

Si sustituimos $f(x) = x^2 - n$ , entonces después de simplificar la expresión, obtenemos:

$x_{i+1} = \frac{x_i + \frac{n}{x_i}}{2}$

La primera variante típica del problema es cuando se da un número racional $n$ , y hay que calcular su raíz con alguna precisión eps:


double sqrt_newton(double n) {
	const double eps = 1E-15;
	double x = 1;
	for (;;) {
		double nx = (x + n / x) / 2;
		if (abs(x - nx) < eps)
			break;
		x = nx;
	}
	return x;
}

Otra variante común del problema es cuando necesitamos calcular la raíz entera (para el $n$ dado hallar el mayor $x$ tal que $x^2 \le n$ ). Aquí es necesario cambiar un poco la condición de terminación del algoritmo, porque puede ocurrir que $x$ empiece a “saltar” cerca de la respuesta. Por lo tanto, añadimos una condición de que si el valor $x$ disminuyó en el paso anterior, e intenta aumentar en el paso actual, entonces el algoritmo debe detenerse.


int isqrt_newton(int n) {
	int x = 1;
	bool decreased = false;
	for (;;) {
		int nx = (x + n / x) >> 1;
		if (x == nx || nx > x && decreased)
			break;
		decreased = nx < x;
		x = nx;
	}
	return x;
}

Por último, se da la tercera variante: para el caso de aritmética de enteros grandes (bignum). Como el número $n$ puede ser lo bastante grande, tiene sentido prestar atención a la aproximación inicial. Obviamente, cuanto más cerca esté de la raíz, más rápido se alcanzará el resultado. Es suficientemente simple y efectivo tomar la aproximación inicial como el número $2^{\textrm{bits}/2}$ , donde $\textrm{bits}$ es el número de bits del número $n$ . Aquí está el código Java que demuestra esta variante:


public static BigInteger isqrtNewton(BigInteger n) {
	BigInteger a = BigInteger.ONE.shiftLeft(n.bitLength() / 2);
	boolean p_dec = false;
	for (;;) {
		BigInteger b = n.divide(a).add(a).shiftRight(1);
		if (a.compareTo(b) == 0 || a.compareTo(b) < 0 && p_dec)
			break;
		p_dec = a.compareTo(b) > 0;
		a = b;
	}
	return a;
}

Por ejemplo, este código se ejecuta en $60$ milisegundos para $n = 10^{1000}$ , y si quitamos la selección mejorada de la aproximación inicial (empezando simplemente con $1$ ), entonces se ejecutará en unos $120$ milisegundos.

Problemas de práctica

UVa 10428 - The Roots