Métodos iterativos para resolver ecuaciones (VI): el método del punto fijo

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, parte IV: método de la secante, parte V: método de regula falsi).

Otro método iterativo para la resolución numérica de ecuaciones es el denominado método del punto fijo. Este método e basa en convertir la ecuación f(x)=0 en otra de la forma g(x)=x, despejando alguna de las ocurrencias de x en la ecuación original.

En el ejemplo que estamos utilizando en esta serie de entradas

 

podemos despejar el término de mayor grado, quedando la ecuación

De forma gráfica podemos ver la ecuación g(x)=x como el punto de corte de la función g(x) con la recta y=x.

Comenzando en un punto inicial x[0], las iteraciones consisten en ir aplicando la función g(x) al resultado de la anterior iteración (x[1]=g(x[0]), x[2]=g(x[1]), ...) hasta que la diferencia entre dos iteraciones sucesivas sea tan pequeña como queramos. El método se llama del punto fijo porque una solución de la ecuación original debe cumplir que sea invariante al aplicarle la función g(x) (es decir, g(x)=x).

Este método algorítmicamente es muy sencillo, en wxMaxima se implementa con las instrucciones:

while abs(x[0]-x[1])>prec do(x[1]:x[0], x[0]:float(g(x[0])));

Donde x[0] es el punto inicial, prec la precisión definida y x[1] una variable auxiliar para controlar la diferencia entre iteraciones consecutivas.

Gráficamente podemos ver el recorrido del algoritmo en las siguientes animaciones de ejemplo:

 

Sin embargo, el método del punto fijo tiene una gran exigencia para asegurar la convergencia: que la derivada de la función g(x) sea, en valor absoluto, menor que 1 en un entorno del punto de cada iteración.

Podemos observar gráficamente un caso de no convergencia en la siguiente animación:

Incluimos esta consideración en el algoritmo (en negrita) y aprovechamos para ejecutarlo completo con el ejemplo de esta serie.

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

q(x):=(-3*x^4+2*x^3+x^2-5*x+3)^(1/5);

define(q1(x),diff(q(x),x,1));

prec:10^(-4);

x[0]:2;

x[1]:x[0]+2*prec;

while abs(x[0]-x[1])>prec do(if abs(q1(x[0]))>=1 then print("El método puede no converger"), x[1]:x[0], x[0]:float(q(x[0])));

    El método puede no converger

    done

float(x[0]);

    −3.350772866871043

Obsérvese que en una de las iteraciones la derivada del punto era mayor o igual que uno por lo que ha salido el mensaje avisando que el método puede no converger. A pesar de esto, el método ha convergido hacia una de las soluciones de la ecuación inicial.

Personalmente, de todos los métodos analizados hasta el momento en esta serie, el método del punto fijo es el que menos me ha gustado, pero siempre está bien conocer un método más para poder valorar su idoneidad en función de cada caso.

Para finalizar esta entrada dejo algunas preguntas abiertas por si alguien quiere investigar un poco más. ¿Variando el punto inicial el método converge hacia otra solución? ¿Qué pasa si despejamos la x de cualquier otro término de la ecuación?

Métodos iterativos para resolver ecuaciones (V): el método de "regula falsi"

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, parte IV: método de la secante).

Existe otro método iterativo para resolver numéricamente ecuaciones denominado "regula falsi". Una vez entendidos el método de bisección y el método de la secante, el método de "regula falsi" es sencillo de comprender porque se basa en ambos. Por una parte, realiza el mismo proceso de calcular el punto de corte de la secante a la función en dos puntos con el eje X. Una vez hallado el punto de corte, utiliza el criterio de los signos opuestos de las imágenes de las abscisas (del método de bisección) para decidir cuál de las dos abscisas anteriores se mantiene en la nueva iteración junto con la nueva abscisa proveniente del punto de corte calculado.

Conceptualmente no tengo mucho más que añadir, sin embargo sí que quiero compartir algunos detalles de la algoritmia del método. En un principio pensé que escribir el bucle del procedimiento sería poco más que un "copiar y pegar" de los dos procedimientos que ya están implementados en entradas anteriores (bisección y secante). Sin embargo solamente añadir el criterio de los signos opuestos de las imágenes de las abscisas al final del bucle del método de la secante no funciona y a continuación explico la razón.

La razón principal es que con este método, a diferencia del método de bisección, el intervalo en el que trabajamos en las sucesivas iteraciones no tiene por qué estrecharse tanto como se quiera para alcanzar una precisión predefinida. En la figura anterior puede observarse un caso en el que el intervalo se estrecha sólo por uno de los dos extremos. ¿Qué consecuencia tiene esto? Pues que el criterio que estábamos utilizando en métodos anteriores sobre la cercanía de los extremos de los intervalos para decidir cuándo se había alcanzado la precisión deseada no sirve.

Así que introduciendo la comprobación inicial de que las abscisas con las que empezamos tengan imágenes de signo opuesto, entramos en un bucle que se repite mientras los extremos de los intervalos disten más que la precisión definida. Dentro del bucle se calcula el punto de corte de la secante con el eje X y se determina el nuevo intervalo en función de dónde quede la alternancia de signo de las imágenes de las abscisas. Aquí se añade una comprobación más que en caso de cumplirse se sale del bucle: si dos iteraciones consecutivas producen abscisas tan cercanas como se haya definido en la precisión (¿podría esto parar el bucle en algunos casos antes de estar cerca de la solución?).

Las instrucciones en wxMaxima quedan de la siguiente manera:
x[0]:0; x[1]:1;
x[3]:x[0];
prec:10^(-4);
if p(x[0])*p(x[1]) < 0 then 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]), if p(x[0])*p(x[2]) < 0 then x[1]:x[2] else x[0]:x[2], if abs(x[3]-x[2])
float(x[2]);

Que en ejemplo que estamos resolviendo en esta serie

devuelve como solución .6629131443657456

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.