Métodos iterativos de resolución de sistemas lineales

En esta práctica nos ocuparemos de tres métodos iterativos de resolución de sistemas lineales: Jacobi, Gauss-Seidel y relajación.

Jacobi

Vamos a considerar un sistema lineal, especificado a partir de su matriz de coeficientes y de su vector de términos independientes, y construiremos aproximaciones de la solución mediante el método de Jacobi.

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}; \)

Determinamos el orden de la matriz del sistema.

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

Especificamos el vector inicial.

x = Table[0., {i, n}];

Establecemos el número de iteraciones.

iteraciones = 50;

Producimos las aproximaciones deseadas.

\(y = Table[0., {i, n}]; 

    For[m = 1, m ≤ iteraciones, \(m++\), 
      
\ \ \ \ \ \ \ 
      For[i = 1, i ≤ n, \(i++\), 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        y[[i]] = 
          \(b[[i]] - 
              ∑\+\(j = 1\)\%\(i - 1\)A[[i, j
                    ]]\ 
                  x[[j]] - 
              ∑\+\(j = i + 1\)\%n
                   A[[i, j]]\ 
                  x[[j]]\)\/A
              [[i, i]]]; 
      
\ \ \ \ \ \ \ x = y]\)

La aproximación de la única solución del sistema considerado, tras 50 iteraciones es

x

La solución exacta es

LinearSolve[A, b]

Como a priori no podemos conocer el número de iteraciones que deben realizarse, podemos proceder con otro criterio para determinar el momento en que hay que detener el cálculo. Establecemos como criterio de parada, por ejemplo, que la norma euclídea de la diferencia entre los resultados de dos iteraciones consecutivas sea inferior a un valor dado, tolerancia, prefijando también un número máximo de iteraciones, maximoiteraciones.

\(maximoiteraciones = 50; 
tolerancia = 1\/10.\^5; \)

Para generar las iteraciones debemos volver a definir el vector inicial, pues fue modificado en el procedimiento de cálculo anterior.

\(x = Table[0., {i, n}]; 
y = Table[0., {i, n}]; 
m = 0; 

    dif = 2\ tolerancia; 

    While[m ≤ maximoiteraciones\  && \ 
        dif ≥ tolerancia, 
      
For[i = 1, i ≤ n, \(i++\), 
        
\ \ \ \ \ \ \ 
        y[[i]] = 
          \(b[[i]] - 
              ∑\+\(j = 1\)\%\(i - 1\)A[[i, j
                    ]]\ 
                  x[[j]] - 
              ∑\+\(j = i + 1\)\%n
                   A[[i, j]]\ 
                  x[[j]]\)\/A
              [[i, i]]]; 
      
dif = 
        \@\(∑\+\(i = 1\)\%n
              \((y[[i]] - 
                  x[[i]])\)\^2\); 
      
x = y; 
m = m + 1]\)\(Print["METODO DE JACOBI"]; 
Print["================"]; 

    Print["APROXIMACION DE LA SOLUCION:"]; 
Print[x]; 

    Print["NUMERO DE ITERACIONES REALIZADAS:    ", m - 1]; 

    Print["----------------"]; 

    Print["NUMERO MAXIMO DE ITERACIONES:        ", maximoiteraciones]; 

    Print["TOLERANCIA PREFIJADA:                ", tolerancia]; 

    Print["----------------"]\)\(If[dif < tolerancia, Print["Objetivo de error:      LOGRADO"], 
      
\ \ 
      Print["Objetivo de error:      NO LOGRADO"]]; 

    Print["Diferencia entre las dos aproximaciones finales: ", \ dif]; 

    \)

Gauss-Seidel

De nuevo operaremos como en el caso anterior. Todas las órdenes se incluyen en una única celdilla.

\(x = Table[0., {i, n}]; 
y = Table[0., {i, n}]; 

    maximoiteraciones = 10; 
tolerancia = 1\/10.\^4; 
m = 0; 

    dif = 2\ tolerancia; 

    While[m ≤ maximoiteraciones\  && \ 
        dif ≥ tolerancia, 
      
\ \ \ \ \ \ \ \ \ \ \ \ \ 
      For[i = 1, i ≤ n, \(i++\), 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        y[[i]] = 
          \(b[[i]] - 
              ∑\+\(j = 1\)\%\(i - 1\)A[[i, j
                    ]]\ 
                  y[[j]] - 
              ∑\+\(j = i + 1\)\%n
                   A[[i, j]]\ 
                  x[[j]]\)\/A
              [[i, i]]]; 
      
\ \ \ \ \ \ \ \ \ \ \ \ \ 
      dif = \@\(∑\+\(i = 1\)\%n
              \((y[[i]] - 
                  x[[i]])\)\^2\); 
      
\ \ \ \ \ \ \ \ \ \ \ \ \ x = y; 
      
\ \ \ \ \ \ \ \ \ \ \ \ \ m = m + 1]\)\(Print["METODO DE GAUSS-SEIDEL"]; 

    Print["======================"]; 

    Print["APROXIMACION DE LA SOLUCION:"]; 
Print[x]; 

    Print["NUMERO DE ITERACIONES REALIZADAS:        ", m - 1]; 

    Print["----------------"]; 

    Print["NUMERO MAXIMO DE ITERACIONES:            ", 
      maximoiteraciones]; 

    Print["TOLERANCIA PREFIJADA:                    ", tolerancia]; 

    Print["----------------"]\)\(If[dif < tolerancia, Print["Objetivo de error:     LOGRADO"], 
      
\ \ \ 
      Print["Objetivo de error:     NO LOGRADO"]]; 

    Print["Diferencia entre las dos aproximaciones finales: ", \ dif]; 
    \)

Relajación

Con el método de relajación procedemos de la misma forma, especificando el valor del parámetro w deseado.

Un caso concreto

\(w = 1.1; 
x = Table[0., {i, n}]; 
y = Table[0., {i, n}]; 

    maximoiteraciones = 25; 
tolerancia = 1\/10.\^3; 
m = 0; 

    dif = 2\ tolerancia; 

    While[m ≤ maximoiteraciones\  && \ 
        dif ≥ tolerancia, 
      
\ \ \ \ \ \ \ \ \ \ \ \ 
      For[i = 1, i ≤ n, \(i++\), 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        y[[i]] = 
          \(w\ \((b[[i]] - 
                    ∑\+\(j = 1\)\%\(i - 1\)A[[i, j
                          ]]\ 
                        y[[j]] - 
                    ∑\+\(j = i + 1\)\%n
                         A[[i, j]]\ 
                        x[[j]])\)\)\/A
                [[i, i]] + 
            
\((1 - w)\)\ 
              x[[i]]]; 
      
\ \ \ \ \ \ \ \ \ \ \ \ 
      dif = \@\(∑\+\(i = 1\)\%n
              \((y[[i]] - 
                  x[[i]])\)\^2\); 
      
\ \ \ \ \ \ \ \ \ \ \ \ x = y; 
      
\ \ \ \ \ \ \ \ \ \ \ \ m = m + 1]\)\(

Print["METODO DE RELAJACION"]; 

    Print["======================"]; 

    Print["APROXIMACION DE LA SOLUCION:"]; 
Print[x]; 

    Print["NUMERO DE ITERACIONES REALIZADAS:    ", m - 1]; 

    Print["----------------"]; 

    Print["NUMERO MAXIMO DE ITERACIONES:        ", maximoiteraciones]; 

    Print["TOLERANCIA PREFIJADA:                ", tolerancia]; 

    Print["----------------"]\)\(If[dif < tolerancia, Print["Objetivo de error:      LOGRADO"], 
      
\ \ 
      Print["Objetivo de error:      NO LOGRADO"]]; 
\ \ \ \ 
    Print["Diferencia entre las dos aproximaciones finales: ", \ dif]; 
    \)

Ejercicios

1.-Obtenga aproximaciones  tanto por Jacobi como por Gauss-Seidel, efectuando 20 iteraciones, un sistema de ecuaciones de orden 30 para el cual todos los elementos de la diagonal valen 9, y los elementos de las dos diagonales siguientes ( a ambos lados de la diagonal principal) son todos iguales a 2. Los restantes coeficientes son nulos. Las componentes del vector de términos independientes son (Subscript[b, i]) = i - 15, 1 ≤ i ≤ 30. Debe dar algunos valores con punto decimal al objeto de trabajar con representaciones decimales y que el problema se resuelva más rápidamente.

2.- Obtenga aproximaciones por Jacobi y por Gauss-Seidel el sistema cuya matriz de coeficientes es 

		A=(({{2, 2, 3, 1, 3, (39/8)}, {2, 1, 2, 2, 5, (39/8)}, {1, 1, 1, 1, 4, (39/8)}, {0, 0, 1, 1, 5, (39/8)}, {-1, 0, -1, 1, -1, 3}, {1, -1, 3, -4, 5, 39}}))

y cuyos términos independientes son

		\(TraditionalForm
      \`b = \((5.1, \(-5.1\), 4.1, \(-4.1\), 111.1, 3)\)\^T\),

partiendo de la aproximación inicial  x = (((6, 6, 6, 6, 6, 6))^T) y efectuando 100 iteraciones.

3.- Obtenga las soluciones de los ejercicios 1 y 2 usando la sentencia LinearSolve y también invirtiendo la matriz de coeficientes.