Métodos iterativos para resolver ecuaciones (IV): el método de la secante

Esta entrada pertenece a la serie "Métodos iterativos para resolver ecuaciones" (parte I: introducción, parte II: método de Newton, parte III: método de bisección).

En el método de Newton que vimos en una entrada anterior calculábamos la recta tangente a la función en los sucesivos puntos de las iteraciones. Para ello necesitábamos que la función fuera derivable.

Existe un método parecido que en lugar de utilizar la recta tangente a la función en un punto calcula la recta secante a la función en dos puntos. Una ventaja es que no necesitamos que la función sea derivable.

Dicho método se conoce con el nombre de método de la secante. Entonces partiendo de dos abscisas, a y b, calculamos la recta secante a la función en dichos puntos. Es decir, buscamos la recta que pasa por los puntos (a,f(a)) y (b,f(b)).

Hay varias formas de calcular dicha recta secante. Por ejemplo, podemos hallar la pendiente de la recta a partir de los puntos por los que pasa: (f(b)-f(a))/(b-a) y utilizar la ecuación punto-pendiente:
y-f(a) = [(f(b)-f(a))/(b-a)]·(x-a)
[Nota: también podíamos haber escogido el punto (b,f(b)) y la ecuación hubiera quedado y-f(b) = [(f(b)-f(a))/(b-a)]·(x-b), pero la recta sigue siendo la misma]

Si ahora buscamos el punto de corte de la recta secante con el eje X (y=0) y despejamos la x tenemos:
-f(a) = [(f(b)-f(a))/(b-a)]·(x-a)
-f(a)/[(f(b)-f(a))/(b-a)] = (x-a)
a - f(a)/[(f(b)-f(a))/(b-a)] = x
a - f(a)·(b-a)/(f(b)-f(a)) = x

Esta x es la nueva abscisa, la llamamos c, que substituye en la próxima iteración a la abscisa a. Es decir, ahora se repetiría el proceso con los puntos (b,f(b)) y (c,f(c)).

Escrito en forma de sucesión, la iteración general queda definida, empezando con x[0]=a y x[1]=b, de la siguiente manera:
x[n]= x[n-2] - f(x[n-2])·(x[n-1]-x[n-2])/(f(x[n-1])-f(x[n-2]))
[Nota: si cogemos la otra ecuación de la recta secante la expresión general de la iteración queda x[n]= x[n-1] - f(x[n-1])·(x[n-1]-x[n-2])/(f(x[n-1])-f(x[n-2])), ambas expresiones son igual de válidas]

Vamos a utilizar wxMaxima para aplicar el método con la ecuación de la primera entrada de la serie:

Definimos el polinomio y los puntos iniciales (por ejemplo 1 y 2):
p(x):=x^5+3*x^4−2*x^3−x^2+5*x−3;
x[0]:1; x[1]:2;

Realizamos la primera iteración y mostramos la aproximación decimal del resultado:
x[2]: x[0] - p(x[0])*(x[1]-x[0])/(p(x[1])-p(x[0]));
float(%);
0.953125

Gráficamente lo que hemos hecho en la iteración anterior es:

Realizamos unas cuantas iteraciones más:
x[3]: x[1] - p(x[1])*(x[2]-x[1])/(p(x[2])-p(x[1]));
float(%);
.9144356310802643

x[4]: x[2] - p(x[2])*(x[3]-x[2])/(p(x[3])-p(x[2]));
float(%);
.7451112080395356

x[5]: x[3] - p(x[3])*(x[4]-x[3])/(p(x[4])-p(x[3]));
float(%);
.6868819480267324

x[6]: x[4] - p(x[4])*(x[5]-x[4])/(p(x[5])-p(x[4]));
float(%);
.6651550727881196

x[7]: x[5] - p(x[5])*(x[6]-x[5])/(p(x[6])-p(x[5]));
float(%);
.6629972438289449

x[8]: x[6] - p(x[6])*(x[7]-x[6])/(p(x[7])-p(x[6]));
float(%);
.6629400060916464

Podemos observar que las 3 últimas iteraciones mantienen dos cifras decimales (0.66) por lo que damos por válida dicha aproximación (por truncamiento) de la solución de la ecuación.

El proceso anterior puede automatizarse en wxMaxima de la siguiente manera:

- Definimos las abscisas iniciales (en nuestro ejemplo a=1 y b=2)
x[0]:1; x[1]:2;

- Definimos la precisión deseada en la aproximación decimal de la solución:
prec:10^(-4);

- Mientras las abscisas disten más que la precisión definida se van haciendo iteraciones, salvo que en alguna iteración se encuentre la solución exacta:
while abs(x[1]-x[0])>prec do (x[2]: float(x[0] - p(x[0])*(x[1]-x[0])/(p(x[1])-p(x[0]))), if p(x[2])=0 then return(x[2]), x[0]:x[1], x[1]:x[2]);

- Si el bucle no para con la solución de la ecuación, pedimos a wxMaxima que nos muestre la aproximación decimal de la última iteración:
float(x[2]);
.6629400060916464

Al igual que en los otros métodos, te propongo que pruebes con otros valores iniciales para encontrar las otras soluciones de la ecuación.