Resolución de Sistemas Lineales
mediante métodos iterativos.

En esta práctica nos ocuparemos principalmentalmente de tres de los  métodos iterativos más empleados para la resolución de sistemas lineales: Jacobi, Gauss-Seidel y relajación. Pero antes definiremos lo que se entiende por norma matricial (en concreto las inducidas por una norma vectorial) e  introduciremos el concepto fundamental de número de condición de una matriz, así como el radio espectral de la misma (que no es más que el máximo módulo o valor absoluto de sus valores propios).

Normas matriciales

Nos centraremos ahora sólo en los conceptos útiles desde el punto de vista de los métodos numéricos. Entre éstos resultarán fundamentales las denominadas normas matriciales, que nos permitirán medir distancias entre matrices y comprobar la convergencia de diferentes métodos para el cálculo aproximado, que veremos a lo largo de este libro.

Por ejemplo la p-norma de un vector {(Subscript[v, 1]), (Subscript[v, 2]), ...} viene dada por la fórmula    ((((∑_i)|(Subscript[v, i])(|^p)))^((1/p))),  si  1 ≤ p < ∞, y por  \(TraditionalForm
      \`max\+i\ \( | \(v\_i | \)\)\) si p = ∞.  Por otro lado, para una matriz matA, la  p-norma se define como  el  máximo, sobre todos los vectores vec de p-norma  unitarios, de la p-norma del vector producto "\!\(\*StyleBox[\"\\\"mat\\\"\", \"TI\"]\)"A   "\!\(\*StyleBox[\"\\\"vec\\\"\", \"TI\"]\)".  Para la norma de la inversa de una matriz, se puede calcular explícitamente la norma de dicha inversa, pero resulta mucho más rápido trabajar con la correspondiente descomposición de tipo  LU  de dicha matriz. Así pues, tanto matA o alguna descomposición  LU de la misma pueden emplearse como argumento de la sentencia InverseMatrixNorm.  Otras sentencias de Mathematica  relacionadas con el cálculo de la norma de un vector, una matriz o de su inversa se muestran en la siguiente tabla.  Para emplear cualquiera de ellas hay que cargar previamente el paquete  LinearAlgebra`MatrixManipulation`.

({{"\!\(\*StyleBox[\"\\\"VectorNorm[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"vec\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", "norma del m\[AAcute]ximo de un vector de componentes\n\
num\[EAcute]ricos aproximados,"}, {"\!\(\*StyleBox[\"\\\"VectorNorm[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"vec\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\",\\\"\", \"MR\"]\) \!\(\*StyleBox[\"\\\"p\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", p"\[Hyphen]norma  (1\[LessEqual]p<\[Infinity]) de un \
vector, "" "}, {"\!\(\*StyleBox[\"\\\"MatrixNorm[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"mat\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", "norma del m\[AAcute]ximo de una matriz de \
n\[UAcute]meros\naproximados,"}, {"\!\(\*StyleBox[\"\\\"MatrixNorm[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"mat\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\",\\\"\", \"MR\"]\) \!\(\*StyleBox[\"\\\"p\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", p"\[Hyphen]norma (p=1,2 o \[Infinity]) de una matriz \
de n\[UAcute]meros\naproximados,  "" "}, {"\!\(\*StyleBox[\"\\\"InverseMatrixNorm[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"mat\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", "norma del m\[AAcute]ximo de la inversa de una matriz \
de\nn\[UAcute]meros\:f39faproximados, "}, {"\!\(\*StyleBox[\"\\\"InverseMatrixNorm[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"lu\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", norma del máximo  de la inversa de una 
                        matriz cuya 
descomposición LU es  lu,}, {"\!\(\*StyleBox[\"\\\"InverseMatrixNorm[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"mat\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\",\\\"\", \"MR\"]\) \!\(\*StyleBox[\"\\\"p\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\), \!\(\*StyleBox[\"\\\"InverseMatrixNorm[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"lu\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\",\\\"\", \"MR\"]\) \!\(\*StyleBox[\"\\\"p\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", p"\[Hyphen]norma de la inversa de una matriz" numérica}, {"\!\(\*StyleBox[\"\\\"MatrixConditionNumber[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"mat\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", "n\[UAcute]mero de condici\[OAcute]n (respecto a la \
norma \[Infinity]) de\nuna matriznum\[EAcute]rica"}, {"\!\(\*StyleBox[\"\\\"MatrixConditionNumber[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"mat\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\",\\\"\", \"MR\"]\) \!\(\*StyleBox[\"\\\"p\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", "n\[UAcute]mero de condici\[OAcute]n (para la p-norma, \
p=1,2 o \[Infinity] )\nde una matriz num\[EAcute]rica " " " " "}})

<< LinearAlgebra`MatrixManipulation`

Veamos algunos ejemplos de normas vectoriales y matriciales:

v = {1, 2, 3, 4, -10};

VectorNorm[v]

Nótese que esta sentencia no devuelve resultado alguno en el caso de introducir un vector cuyas componentes no estén dadas en coma flotante (con decimales). Lo mismo pasará en el caso de las matrices.

También dará error si introducimos el vector en forma de columna en vez de como un vector fila:

v=(({{1.}, {2}, {3}, {4}, {-10}}));

VectorNorm[v]

No obstante, la orden Flatten elimina las posibles sublistas. Aplicada a v, evitará el problema anterior.

VectorNorm[Flatten[v]]

Análogamente, si el argumento de VectorNorm se introduce como vector fila con  alguna componente con el punto decimal, no se producirá ningún problema, obteniendo el valor absoluto máximo de las componentes.

v = {1., 2, 3, 4, -10};

VectorNorm[v]

No obstante, no es imprescindible usar la orden VectorNorm para hallar una norma, pues se puede utilizar directamente la definición:

p = 2;  n = Length[v]; 
(((∑_(k = 1)^n) Abs[v[[k]]]^p)^(1/p))

Todo esto lo podemos incluir ahora en nuestro propio módulo, mediante la sentencia  Module que permite definir variables locales y realizar un bloque de programa de una vez, pudiéndolo llamar todas las veces que necesitemos.

normap[v_List, p_] := 
      Module[{w, n}, 
w = Flatten[v]; 
        
n = Length[w]; 
        
If[p == ∞, Max[Abs[w]], 
          (((∑_(k = 1)^n) Abs[w[[k]]]^p)^(1/p))]]

normap[v, 2]

normap[v, ∞]

Veamos ahora algunos ejemplos con matrices:

B=(({{2, 2, 3}, {4, 6, 6}, {7, 8, 10}}));

MatrixNorm[B]

Nótese que esta sentencia no devuelve resultado alguno porque ninguna componente a una representación decimal (defina como tal una cualquiera de ellas y calcule las normas correspondientes a p = 1, 2 e ∞).

En esta ocasión también podemos optar por programar una subrutina propia para calcular la correspondiente norma matricial inducida por las normas vectoriales con  p = 1, 2  o  ∞.

normap[matriz_List, p_] := 
      Module[{m, n, valor}, 
        
m = Dimensions[matriz][[1]]; 
        
n = Dimensions[matriz][[2]]; 
        
If[m == n, 
     
          If[p == 1, 
            valor = Max[
                Table[(∑_(i = 1)^n) Abs[ matriz[[i, j]] ], {j, 
                    n}]]]; 
   
          If[p == 2, 
            valor = √Max[Abs[Eigenvalues[Transpose[matriz].matriz]]]]; 
          
 
          If[p == ∞, 
            valor = Max[
                Table[(∑_(j = 1)^n) Abs[matriz[[i, j]] ], {i, 
                    n}]]]
]; 
valor
        
]

Veamos algunos ejemplos de aplicación.

normap[B, 2]

normap[B // N, 1]

normap[B // N, ∞]

O bien, dar la matriz siempre con componentes en coma flotante y emplear las órdenes directas incluidas en el paquete   LinearAlgebra`MatrixManipulation`.

MatrixNorm[B // N]

MatrixNorm[B // N, ∞]

MatrixNorm[B // N, 1]

MatrixNorm[B // N, 2]

Inverse[B // N] // MatrixNorm

InverseMatrixNorm[B // N]

MatrixConditionNumber[B // N]

MatrixConditionNumber[B // N, ∞]

MatrixConditionNumber[B // N, 1]

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.