Métodos iterativos para resolver sistemas de ecuaciones (lineales y no lineales)

Métodos iterativos para resolver sistemas de ecuaciones (lineales y no lineales)

Los métodos iterativos se caracterizan porque buscan acercarse al conjunto solución por medio de aproximaciones que van llevando sistemáticamente a la convergencia. A diferencia de los métodos directos (Eliminación de Gauss, Eliminación de Gauss-Jordan, Inversa-Multiplicación, Descomposición LU, Cramer) donde se aplica un algoritmo de pasos finitos para llegar a la solución, los métodos iterativos inician con un valor supuesto de la solución y por medio de aproximaciones y la evaluación de un criterio de convergencia se determina que se ha llegado a la solución más cercana, es decir que se puede tener un rango de error permitido de la solución.
¿cuándo aplicar métodos iterativos en lugar de métodos directos? si tomamos en cuenta la cantidad de operaciones que se deben realizar en un método directo para llegar a la solución el valor es de n^3 operaciones donde n es el número de lineas en la matriz o el número de ecuaciones en el sistema. Si hablamos de un sistema de 100 o 200 ecuaciones entonces se requiere de 1,000,000 u 8,000,000 operaciones para resolver esos sistemas respectivamente. Es de considerar el trabajo computacional que deben hacer las máquinas o el ser humano, si lo resolvemos a papel y lápiz, para llegar a la solución, debemos entonces buscar alternativas para llegar a la solución de una manera más eficiente además de eficaz.
La estrategia que siguen los métodos iterativos es encontrar una ecuación que nos vaya llevando sistemáticamente a la solución en cada iteración dado un valor inicial, la ecuación iterativa se parece un poco al método de sustitución o punto fijo que se aplica para ecuaciones no lineales, donde se tiene una ecuación con la forma x=g(x), se da un valor inicial de x que se sustituye en g(x), el resultado se vuelve a dar como parámetro a g(x) hasta que el valor sustituido y el valor obtenido tienen una mínima diferencia |x-g(x)| es menor que una tolerancia.

Método de Jacobi

En el método de Jacobi se descompone la matriz de coeficientes A en dos matrices, una se compone de los elementos de la diagonal D(i,i)=A(i,i). La matriz D no debe tener elementos cero, si fuera así se deben intercambiar los renglones para evitar ese caso. La otra matriz es el resto de los elementos, es decir los elementos que no son la diagonal R(i,j)=A(i,j) de tal manera que A = D + R.
Por lo tanto la ecuación que describe un sistema de ecuaciones lineales Ax = b ahora se expresa como [D+R]x = b.
La ecuación iterativa de Jacobi se inicia con valores supuestos x(k) y la ecuación iterativa obtiene los siguientes valores x(k+1). El proceso se repite hasta que la diferencia absoluta entre x(k+1) y x(k) sea mínima.

El siguiente código es el de la función de Jacobi, se debe crear un archivo .m en la opción de New Script para que en la ventana de edición se copie el código de abajo y se gurade con el nombre jacobi.m y debe quedar en la carpeta de trabajo o current folder

function x = jacobi(A,b,x,imax,tol)
if nargin==3
   imax=100;
   tol=1e-8;
end
D=diag(A);
R=A-diag(D);
cumple=0;
k=0;
while ( ~cumple && k<imax ) 
   xk1 =(b-R*x)./D; 
   norma = norm(x-xk1 ); 
   fprintf('iteracion :%3.0f norma: %e\n',k, norma );
   cumple = norma < tol; 
   x= xk1; 
   k =k+1; 
end 
if k >= imax
   error('El sistema no converge ');
end
end

ejemplo del uso del método de Jacobi (ventana command window de matlab)

A = [-512,12,0,0,0;
     500,-512,12,0,0;
     0,500,-512,12,0;
     0,0,500,-512,12;
     0,0,0,-500,512];
b = [0,0,0,0,9]';
x = [0,0,0,0,0]'; %valores iniciales
format longE
x = jacobi(A,b,x)

Un método iterativo puede converger a la solución pero también puede diverger, lo que significa que la diferencia absoluta entre x(k+1) y x(k) no se minimice sino que sea cada vez más grande por lo tanto no llegará a la solución. Para que un sistema llegue a la convergencia debe ser diagonalmente dominante, esto quiere decir que los elementos de su diagonal (la matriz D) debe ser mayor en valor absoluto a la suma de los absolutos de su mismo renglón o columna de la matriz R.

D(i,i) debe ser mayor que R(j,i) j desde 1 hasta n
D(i,i) debe ser mayor que R(i,j) j desde 1 hasta n

para toda i desde 1 hasta n

Tomemos por ejemplo la matriz A que tiene los valores

La matriz D son los elementos de la diagonal, es decir lo elementos de color naranja

La matriz R son el resto de los elementos que no son la diagonal

De tal manera que A=D+R. Cuando se dice que el sistema debe ser diagonalmente dominante se refiere a que los elementos de la diagonal deben ser mayores, en valor absoluto, a los elementos del mismo renglón y misma columna. Por ejemplo si revisamos es renglón 1. El elemento de la diagonal es 512 y la suma de los elementos del mismo renglón son 12+0+0+0=12, por lo tanto 512>12 y es diagonalmente dominante.

Ahora analicemos la columna 1, aquí vemos que el elemento de la diagonal es el mismo 512 y la suma de la misma columna es 500+0+0+0+0, entonces 512>500 y es diagonalmente dominante

El mismo análisis se hace para el renglón 2 y columna 2 y todos los demás

En este caso el elemento de la diagonal no es mayor sino igual, eso quiere decir que no es diagonalmente dominante y que mi sistema no va a converger? la respuesta que es si va a converger pero se va a llevar más iteraciones para llegar a la convergencia. Los sistemas que se generan en ingeniería son principalmente diagonalmente dominantes con coeficientes reales, así que no hay mucho tema aquí.

Método de Gauss-Seidel

Revisemos ahora el método de Gauss-Seidel o el método de Liebmann es un método iterativo que separa la matriz de coeficientes A en dos matrices, una matriz triangular inferior L que contiene los elementos de A que están en la diagonal inferior (la diagonal no debe contener ceros) y la matriz U que contiene los elementos de A que están sobre la diagonal.
De tal manera que A = L+U.La ecuación que describe un sistema de ecuaciones lineales Ax=b ahora se expresa como [L+U]x = b.

La ecuación iterativa de Gauss-Seidel se inicia con valores supuestos x(k) y la ecuación iterativa obtiene los siguientes valores x(k+1). El proceso se repite hasta que la diferencia absoluta entre x(k+1) y x(k) sea mínima.

El siguiente código es el de la función de Gauss-Seidel, se debe crear un archivo .m en la opción de New Script para que en la ventana de edición se copie el código de abajo y se gurade con el nombre gauss_seidel.m y debe quedar en la carpeta de trabajo o current folder

function x = gauss_seidel(A,b,x,imax,tol)
if nargin==3
   imax=100;
   tol=1e-8;
end
L= tril(A);
U=A-L;
Linv = inv(L);
cumple = 0;
k=0;
while ( ~cumple && k<imax ) 
   xk1 = Linv*(b-U*x); 
   norma = norm(x-xk1); 
   fprintf('iteracion: %3.0f norma: %e\n',k,norma);
   cumple = norma < tol; 
   x=xk1; 
   k=k+1; 
end 
if k >= imax
   error('El sistema no converge ');
end
end

ejemplo del uso del método de Gauss-Seidel en la ventana command window

A = [-512,12,0,0,0;
     500,-512,12,0,0;
     0,500,-512,12,0;
     0,0,500,-512,12;
     0,0,0,-500,512];
b = [0,0,0,0,9]';
x = [0,0,0,0,0]'; %valores iniciales
format longE
x = gauss_seidel(A,b,x)

Recordemos que el sistema debe ser diagonalmente dominante.

Método SOR

El método de Gauss-Seidel se puede mejorar usando la modificación sugerida por Southwell por la introducción de un factor (llamado factor de relajación), lo que da como resultado el método de sobre-relajación sucesiva (SOR).

Donde w es un factor con los siguientes rangos de valores:
Sub-relajación 0 < w < 1 se logra la convergencia en sistemas divergentes.
Sobre-relajación 1 < w < 2 acelera la convergencia de sistemas convergentes.

function x=sor(A,b,x,w,imax,tol)
if nargin==3
   w=1;
   imax=100;
   tol=1e-8;
end
if nargin==4
   imax=100;
   tol=1e-8;
end
[n,c]=size(A);
cumple=0;
k=0;
while ( ~cumple && k<imax ) 
   xk1=zeros(n,1); 
   for i=1:n 
      s1=A(i,1:i)*xk1(1:i); 
      s2=A(i,i+1:end)*x(i+1:end); 
      xk1(i)=(b(i)-s1-s2)/A(i,i)*w+(1 -w)*x(i); 
   end 
   norma = norm(x-xk1); 
   fprintf('iteracion :%3.0f ->%e\n',k,norma );
   cumple = norma < tol; 
   x= xk1; 
   k =k+1; 
end 
if k > imax
   error ('El sistema no converge ');
end
end

ejemplo del uso del método de Gauss-Seidel

A = [-512,12,0,0,0;
      500,-512,12,0,0
      0,500,-512,12,0;
      0,0,500,-512,12;
      0,0,0,-500,512];
b = [0,0,0,0,9]';
x = [0,0,0,0,0]'; %valores iniciales
w = 1.1;
format longE
x = sor(A,b,x,w)

Sistemas de ecuaciones no lineales

Como vimos, los métodos iterativos resuelven sistemas de ecuaciones lineales, pero se pueden aplicar a los sistemas no lineales, la respuesta es sí. Ya hemos visto que la estrategia del punto fijo o sustitución es eficaz para ecuaciones no lineales porque buscan la intersección de x y g(x) de manera sistemática por lo tanto se pueden aplicar a sistemas de ecuaciones no lineales bajo el mismo concepto.

En el método de Gauss-Seidel multivariable se sigue el mismo procedimiento que en el método de Punto Fijo, los pasos a seguir son:
Despejar xi de la i-ésima ecuación para generar un sistema iterativo.

El vector x(k) es un vector inicial que se sustituye en el sistema para obtener x1(k+1), en el cálculo de x2(k+1) se utiliza el valor recién obtenido de x1(k+1) haciendo la sustitución sucesiva. El proceso se repite hasta cumplir con el criterio de convergencia.

El siguiente código es el de la función de punto fijo multivariable, se debe crear un archivo .m en la opción de New Script para que en la ventana de edición se copie el código de abajo y se gurade con el nombre pfmulti.m y debe quedar en la carpeta de trabajo o current folder

function [x] = pfmulti(g,x,imax,tol)
   if nargin==2
      imax=100;
      tol=1e-8;
   end
   cumple = 0;
   k=0;
   while ( ~cumple && k<imax )
      x0=x;
      x=g(x0);
      norma = norm(x-x0 );
      fprintf('iteracion :%3.0f norma: %e\n',k,norma);
      cumple=norma < tol;
      k =k+1; 
   end 
   if k >= imax
      error ('El sistema no converge ');
   end
end

Para resolver el siguiente sistema

se debe despejar x1 de la primera ecuación y x2 de la segunda para conseguir el sistema

Se debe construir un archivo .m que contenga la definición del sistema. Se crea un archivo g.m con la opción de New Script para contener ambas ecuaciones, el resultado es un vector donde se evalúa cada ecuación del sistema.

function [xk1] = g(x)
  xk1 = x;
  xk1(1) = sqrt(x(2)^2/5) ;
  xk1(2) = 0.25*(sin(x(1))+cos(x(2)) );
end

finalmente se ejecutan las siguientes líneas en la Command window para conseguir el conjunto solución del sistema. El parámetro @g es la función que contiene el sistema a resolver.

x = [0.5, -0.5]'; %Valores iniciales
format longE
x=pfmulti(@g,x)

La salida es el conjunto solución del sistema, que gráficamente es la intersección de ambas ecuaciones, hagamos la gráfica de éste sistema.

syms x1 x2; %variables involucradas en el sistema
sistema=[5*x1^2-x2^2;x2-(sin(x1)+cos(x2))/4]; %creacion del sistema
fimplicit(sistema) %grafica el sistema
hold on
scatter(x(1),x(2),'red','filled') %grafica la solucion del sistema
hold off

Para encontrar la otra solución se deben cambiar los valores iniciales y correr nuevamente el método