Métodos iterativos para resolver ecuaciones (II): el método de Newton

Esta entrada pertenece a la serie "Métodos iterativos para resolver ecuaciones" (parte I).

El método de Newton comienza en una abscisa inicial (x[0]) y consiste en ir haciendo aproximaciones a una solución de la ecuación f(x)=0 mediante el cálculo del punto de corte del eje X con la recta tangente a la función f(x) en dicha abscisa. Por tanto, necesitamos que la función sea derivable para poder aplicar este método.

Vayamos paso a paso. Partimos de la ecuación punto-pendiente de la recta tangente a f(x) en la abscisa inicial x[0]. (Nota: recuerda que la pendiente de la recta tangente a f(x) en la abscisa x[0] es la derivada de la función en x[0])

(y - f(x[0])) = f'(x[0]) · (x - x[0])

Queremos encontrar su punto de corte con el eje X (y=0) por lo que resolvemos el sistema. Nos queda lo siguiente al despejar la x:

-f(x[0]) = f'(x[0]) · (x - x[0])

-f(x[0])/f'(x[0]) = x - x[0]

x[0] - f(x[0])/f'(x[0]) = x

Esa x, el punto de corte de la recta tangente con el eje X, será la nueva abscisa para la iteración. Es decir, cada iteración del método de Newton se calcula de la siguiente manera:

x[n+1] = x[n] - f(x[n]) / f'(x[n])

A modo de ejemplo, vamos a implementar en wxMaxima dicho método con la ecuación que aparece en la primera entrada de esta serie

Definimos la función y calculamos su primera derivada:

p(x):=x^5+3*x^4−2*x^3−x^2+5*x−3;

define(p1(x),diff(p(x),x,1));

Escogemos una abscisa inicial, por ejemplo x0=1 y aplicamos la primera iteración del método de Newton:
x[0]:1;
x[1]:x[0]-p(x[0])/p1(x[0]);
11/14

Gráficamente, lo que realizamos en la primera iteración es lo siguiente:

Así pues, la nueva abscisa devuelta por la primera iteración es 11/14, cuya aproximación decimal es:
float(%);
.7857142857142857

Si la utilizamos para realizar otra iteración obtenemos:
x[2]:x[1]-p(x[1])/p1(x[1]);
1363457/2003603
float(%);
.6805025746118368

Y realizando unas pocas iteraciones más:
x[3]:x[2]-p(x[2])/p1(x[2]);
float(%);
.6632754087092675

x[4]:x[3]-p(x[3])/p1(x[3]);
float(%);
.6629399905083025

x[5]:x[4]-p(x[4])/p1(x[4]);
float(%);
.6629398707510754

Podemos observar que las abscisas se estabilizan (las cifras .662939 son comunes en las últimas dos iteraciones).

¿Cuándo paramos la iteración? Cuando la diferencia entre dos pasos de la iteración sea tan pequeña como nosotros toleremos. O visto de otra manera, cuando dos abscisas consecutivas de la iteración coincidan en tantos decimales como queramos tener de precisión en la aproximación decimal de la solución.

Podemos automatizar todo el proceso anterior en wxMaxima de la siguiente manera. Primero definimos la precisión del metodo, por ejemplo queremos una precisión hasta diezmilésima:
prec:10^(-4);

Ahora definimos la abscisa inicial, en nuestro ejemplo es x[0]=1
x[0]:1;

Y también definimos una variable auxiliar en la que se irán guardando los valores devueltos por las iteraciones. Nos aseguramos de que el valor inicial de dicha variable diste de la abscisa inicial más que la precisión definida anteriormente. Por ejemplo:
x[1]:x[0]+2*prec;

Definimos el bucle que debe ejecutar mientras la diferencia entre dos iteraciones consecutivas diste más que la precisión definida:
while abs(x[1]-x[0])>prec do (x[0]:x[1], x[1]:x[0]-p(x[0])/p1(x[0]));

Una vez ejecutado el bucle pedimos a wxMaxima que nos muestre el resultado
float(x[1]);
.6629398707510757

Finalmente comprobamos que la imagen de la solución obtenida está cercana a 0:
float(p(x[1]));
8.530092782798483*10^(−14)

Faltaría en el bucle while definido anteriormente contemplar la posibilidad de que el algoritmo detecte la entrada en ciclos que no convergen hacia ninguna solución. No voy a entrar en esta ocasión en este aspecto dado que todo lo anterior tiene fines didácticos, con el objetivo de comprender qué es lo que hace en cada iteración el método de Newton.

En realidad wxMaxima ya lleva este procedimiento implementado en sus librerías y no es necesario definir las instrucciones anteriores. Para ello primero debemos cargar el paquete "mnewton" (sólo es necesario cargarlo una vez):

load(mnewton);

Y para usarlo utilizamos la instrucción

mnewton(función, variable, abscisa inicial)

En nuestro ejemplo haríamos:
load(mnewton);
mnewton(p(x),x,1);
[[x=.6629398707510602]]

La última pregunta que cabe hacerse es: ya tengo una aproximación decimal de una solución de la ecuación, ¿hay más soluciones?

Pues a la vista del gráfico sí que tiene más soluciones. ¿Quieres encontrar también sus aproximaciones decimales? ¿Qué tal si pruebas a poner otras abscisas iniciales en el método de Newton?