Métodos iterativos de resolución 
de sistemas lineales

En esta práctica nos ocuparemos de los métodos iterativos para la  resolución de sistemas de ecuaciones lineales, en particular de los métodos de Jacobi, Gauss-Seidel .

Optamos por utilizar las versiones matriciales de los métodos indicados, a fin de utilizar las posibilidades que ofrece Mathematica.

Vamos a considerar un sistema lineal de ecuaciones 

                                         2 x - y       = 2
                                         x+ 6 y -2 z = -4
                                         4 x+3 y-8 z = 5
                                         
para poder clasificarlo, esto es, saber si tiene solución o no y si es única,  introducimos la matriz de coeficientes y el vector de téminos independientes.

A={{2,-1,0},{1,6,-2},{4,3,-8}};
b={2,-4,5};

Det[A]

La matriz de coeficientes tiene rango tres, por lo tanto,  el sistema tiene solución única. Esta solución la podemos calcular directamente con la
 
sentencia    		    LinearSolve[A,b] .

LinearSolve[A, b]

Ya conocemos la solución, lo que nos permitirá comprobar si las aproximaciones que nos dan los métodos son  buenas aproximaciones o no.

Método de Jacobi

Nuestro propósito es aproximar la solución de un sistema de ecuaciones lineales mediante el método de Jacobi. 

Tomamos el sistema anterior

                                         2 x - y       = 2
                                         x+ 6 y -2 z = -4
                                         4 x+3 y-8 z = 5

Definimos la matriz de coeficientes, que no debe tener ningún elemento diagonal nulo, y el vector de términos independientes.

A={{2,-1,0},{1,6,-2},{4,3,-8}};
b={2,-4,5};

n = Dimensions[A][[1]];

El siguiente bucle nos avisa si la matriz no es diagonalmente estrictamente dominante, es decir no hay garantía de convergencia.

For[i = 1, 
      i ≤ n, i++, 
   If[
        Abs[A[[i, i]]] <= 
          Sum[Abs[A[[i, j]]], {j, i - 1}] + 
            Sum[Abs[A[[i, j]]], {j, i + 1, 
                n}], 
Print["A no es diagonalmente \
estrictamente dominante"]; i = n, True]]

\(Vamos\ a\ constuir\ el\ método\ iterativo\ de\ Jacobi\ de\ la\ forma\ \ x\^\((k + 1)\) = Bx\^\((k)\) + c\)

Hay que formar la matriz B y el vector c del método iterativo de Jacobi. Necesitamos formar las matrices D, L y U. El ordenador lo hará automáticamente a partir de la matriz A ya dada.

Ya estamos en condiciones de formar D ,que llamaremos Di , pues recuerde que la letra D está reservada en Mathematica.

Comenzamos con la matriz diagonal. Para formar una matriz diagonal se puede hacer directamente con la sentencia
 
		DiagonalMatrix[lista] 

sólo hay que darle en lista el vector de elementos que van en la diagonal.

Di=DiagonalMatrix[Table[A[[i,i]],{i,n}]];

Continuamos formando las matrices triangulares inferior L y superior  U, respectivamente, de manera que se cumpla la igualdad

                                              A=Di-L-U
                                              
La U la generamos automáticamente a partir de la matriz A y la L despejando de la igualdad anterior

U=Table[0,{i,n},{j,n}];
For[i=1,i<=n-1,i++,
    For[j=i+1,j<=n,j++,
    U[[i,j]]=-A[[i,j]]
       ]
   ]
L=Di-U-A;

Mostramos las matrices.

MatrixForm[A]
MatrixForm[Di]
MatrixForm[L]
MatrixForm[U]

Construimos la matriz del método de Jacobi y el vector correspondiente, y los mostramos.

B=Inverse[Di].(L+U);
c=Inverse[Di].b;
MatrixForm[B]
MatrixForm[c]

Ya estamos en condiciones de producir las iteraciones. Debemos introducir la aproximación inicial, generalmente será (x^((0)))=(0.,0.,0.) y el número de iteraciones que se desea, aquí suponemos que son cinco.

{x[0], y[0], z[0]} = {0. , 0. , 0. };

iteraciones = 5;

For[i = 1, 
      i ≤ 
        iteraciones, i++, 
    {x[i], y[i], z[i]} = 
        B . {x[i - 1], y[i - 1], z[i - 1]} + 
          c; 
    Print[
        i, "   ", {x[i], y[i], 
          z[i]}]
     ]

Cada i-ésima aproximación de la solución va siendo calculada y mostrada en pantalla antes de utilizarla para la siguiente iteración.

Se puede comparar con la solución verdadera.

LinearSolve[A,b]//N

Método de Gauss-Seidel

Vamos a constuir el método iterativo de Gauss - 
        Seidel  de la forma  (x^((k + 1))) = 
      (Bx^((k))) + c

De nuevo operaremos como en el caso anterior. Las matrices Di, L y U siguen siendo válidas, pues la matriz A del sistema es la misma. Recordemos que la matriz del método (Subscript[B, GS])=(((D - L))^(-1))U y el vector (Subscript[c, GS])=(((D - L))^(-1))b

Las siguientes órdenes proporcionan las sucesivas aproximaciones de la única solución del sistema inicialmente considerado tras definir la matriz y el vector de términos independientes del método:

Formación de la matriz y vector del método de Gauss-Seidel y pedimos que los muestre en forma matricial

mGS=Inverse[Di-L].U;
cGS=Inverse[Di-L].b;
MatrixForm[mGS]
MatrixForm[cGS]

Solución inicial de la que se parte y número de iteraciones que se quieren realizar

{x[0], y[0], z[0]} = {0. , 0. , 0. };

iteraciones = 5;

Ya estamos en condiciones de  pedir que se forme una lista de aproximaciones a la solución a partir de la inicial dada

For[i = 1, 
      i ≤ 
        iteraciones, i++, 
{x[i], y[i], z[i]} = 
        mGS . {x[i - 1], y[i - 1], z[i - 1]} + cGS; 
Print[
        i, "   ", {x[i], y[i], z[i]}]]

Ejercicios

1. Dado el sistema
			    x+(1/2)y+(1/2)z=2
(1/2) x + y + (1/2) z  = 2
(1/2) x+(1/2)y+z=2
			
			
(a) ¿Es la matriz de coeficientes del sistema diagonalmente estrictamente dominante?

(b) Partiendo de la aproximación inicial (({{(x^((0)))}, {(y^((0)))}, {(z^((0)))}}))= (({{0.8}, {0.8}, {0.8}})) compruebe que las iteraciones del método de Jacobi oscilan entre los valores 0.8 y 1.2; es decir, el método no convergerá.

(c) Muestre que las iteraciones obtenidas con el método de Gauss-Seidel convergen a la solución exacta x=y=z=1, calculando las ocho primeras iteraciones con 4 decimales.

2. Se considera el sistema de ecuaciones lineales
					
				 4  x+3y=24
3  x + 4  y - z  = 30

        -y + 4  z = -24	
			
tiene como solución exacta x=3 , y=4, z=-5.

(a) Escriba las ecuaciones de los métodos de Jacobi y Gauss-Seidel para tal sistema.

(b) Partiendo de la aproximación inicial  (((1, 1, 1))^t)  confecccione las tablas numéricas con siete iteraciones para cada uno de los métodos utili   zando 6 cifras decimales en los cálculos.

(c) Según se ve en los resultados del apartado anterior, ¿cuál de los métodos citados proporcionará antes la solución del sistema con las cifras usadas?.