Informática de Sistemas
Cálculo Numérico con Mathematica

Tema 3
Resolución numérica de sistemas de ecuaciones lineales

Nota: Se escribe al principio de la práctica el comando que permite borrar todos los símbolos almacenados (en concreto, valores de variables o funciones) desde que se inició la sesión con Mathematica. Así se facilita su uso siempre que sea necesario.

Clear["Global`*"]

La siguiente definición será necesaria en varias ocasiones a lo largo de esta práctica. La anotamos aquí no para que se ejecute ahora, sino para tenerla a la vista y cargarla siempre que se necesite.

ρ[A_]:=Abs[Eigenvalues[A, 1]][[1]]

Sistemas lineales: órdenes directas para su resolución

Mathematica ofrece varias posibilidades para resolver sistemas de ecuaciones lineales. Mostramos aquí las dos más comunes.

Utilizando el comando LinearSolve[A,b], Mathematica permite obtener una solución de un sistema lineal compatible de la forma A x = b, donde A es una matriz m×n, x es el vector cuyas componentes son las n incógnitas y b es un vector con m componentes. Si el sistema es compatible indeterminado, para obtener todas las soluciones hay que utilizar el comando Reduce[Ecuaciones,Incógnitas] (que también es válido para el caso compatible determinado).
Observe que el primer comando exige que se introduzcan las matrices del sistema; pero el segundo exige las ecuaciones.

Ejemplo 1

Resolvemos el sistema lineal Ax = b, donde A=(({{1, 3, -1}, {0, 3, 1}, {-2, 1, -1}})) y b=(({{1}, {3}, {2.5}})), que es compatible determinado, utilizando ambas instrucciones.

A = {{1, 3, -1}, {0, 3, 1}, {-2, 
          1, -1}}; 
b = {1, 3, 5/2}; 
LinearSolve[A, b]

Reduce[{x + 3  y - z == 1, 
        3  y + z == 3, -2 x + y - z == 5/2}, {x, y, 
        z}]

Ejemplo 2

Resolvemos el sistema lineal Ax = b, donde A=(({{2, 2, 2}, {-2, 5, 2}, {-7, 7, 1}})) y b=(({{0}, {0}, {0}})), que es compatible indeterminado.

LinearSolve[{{2, 2, 2}, {-2, 5, 2}, {-7, 7, 1}}, {0, 0, 
        0}]

El comando LinearSolve tan solo ha dado la solución trivial. El comando Reduce proporciona, sin embargo, las infinitas soluciones (que dependen en este caso de un solo parámetro):

Reduce[{2  x + 2  y + 2  z == 0, -2 x + 5  y + 2  z == 
          0, -7 x + 7  y + z == 0}, {x, y, z}]

Ejemplo 3

Veamos cómo responden los comandos LinearSolve y Reduce ante un sistema incompatible, como es el caso de
				 (({{2, -3}, {2, 1}, {3, 2}}))(({{x}, {y}}))=(({{-2}, {1}, {1}})):

LinearSolve[{{2, -3}, {2, 1}, {3, 2}}, {-2, 1, 1}]

Reduce[{2  x - 3  y == -2, 2  x + y == 1, 
        3  x + 2  y == 1}, {x, y}]

Valores propios, vectores propios y radio espectral de una matriz cuadrada: órdenes directas para su obtención. Matrices diagonalizables.

Mathematica dispone de una orden directa para calcular los valores propios de una matriz cuadrada A: el comando Eigenvalues[A]. La orden para calcular un sistema de generadores del espacio engendrado por los vectores propios de A es el comando Eigenvectors[A]. La orden Eigensystem[A] proporciona tanto los valores como los sistemas de generadores de los vectores propios asociados. Por ejemplo:

Ejemplo 1

Obtenga los valores y vectores propios de la matriz B=(({{3, -1, 1}, {0, 2, 1}, {0, 1, 2}})). Como consecuencia, determine también si la matriz B es o no diagonalizable.

B = {{3, -1, 1}, {0, 2, 1}, {0, 1, 
          2}};  
Eigenvalues[B]
Eigenvectors[B]

En el resultado anterior, los valores propios se corresponden con una base de los respectivos vectores propios. Eso se observa mejor con la orden Eigensystem[B] y expresando el resultado en forma de tabla:

Eigensystem[B] // TableForm

El valor propio doble (y dominante) tiene una base de vectores propios formada por dos vectores, luego la matriz es diagonalizable.

Si la matriz no fuera diagonalizable, entonces algún valor propio múltiple (con multiplicidad algebraica mayor que 1) se le haría corresponder al menos un vector propio nulo. Esto es lo que ocurrirá en el ejemplo 3.

El radio espectral de la matriz, como se observa a simple vista, es 3.

Radio espectral

La instrucción Eigenvalues[A,1] proporciona un valor propio dominante de una matriz numérica A: si hay varios sólo elige uno. La instrucción  Eigenvectors[A,1] proporciona un vector propio asociado al valor propio dominante dado por Eigenvalues[A,1]. Por ejemplo, para la matriz del ejemplo 1 se tiene

Eigenvalues[B, 1]
Eigenvectors[B, 1]

Para que el resultado aparezca sin llaves escribimos:

Eigenvalues[B, 1][[1]]
Eigenvectors[B, 1][[1]]

Luego definimos la función radio espectral, que a cada matriz cuadrada numérica A le asocia su radio espectral ρ(A), de la siguiente forma:

ρ[A_]:=Abs[Eigenvalues[A, 1][[1]]]

Por ejemplo:

ρ[B]

Ejemplo 2

Obtenga los valores y vectores propios de la siguiente matriz y deduzca su radio espectral. Determine también si la matriz es o no diagonalizable.
			 A=(({{-9, 1, 2, 1}, {1, 2, 3, 0}, {1, 1, 0, 1}, {1, 2, 1, -3}})).

A = {{-9, 1, 2, 1}, {1, 2, 3, 0}, {1, 1, 0, 1}, {1, 2, 
          1, -3}}; 
Eigensystem[A] // TableForm

N[Eigensystem[A]] // TableForm

Evidentemente, -9.36163731749222` (redondeado a 6 dígitos) es el valor propio estrictamente dominante; y el radio espectral es 9.36163731749222`. Comprobación de esto último:

ρ[A] // N

Observe también que la matriz A es diagonalizable, ya que sus cuatro valores propios son distintos.

Otro modo algo más directo de obtener los resultados de este ejemplo es el siguiente, en el que algún elemento de la matriz se expresa con punto decimal:

A = {{-9. , 1, 2, 1}, {1, 2, 3, 0}, {1, 1, 0, 1}, {1, 2, 
          1, -3}}; 
Eigensystem[A] // TableForm
ρ[A]

En este caso no ha normalizado la última componente de los vectores propios, que antes era igual a 1.

Ejemplo 3

Obtenga los valores y vectores propios y el radio espectral de la matriz B=(({{3, 0, 0}, {0, 2, 0}, {0, 1, 2}})). Determine también si la matriz B es o no diagonalizable.

B = {{3, 0, 0}, {0, 2, 0}, {0, 1, 2}};  
Eigensystem[B] // TableForm

El valor propio doble tiene una base de vectores propios formada por un solo vector, luego la matriz no es diagonalizable.

El radio espectral de la matriz, como se observa a simple vista, es 3, que coincide con el valor propio estrictamente dominante.

Ejemplo 4

Obtenga los valores y vectores propios y el radio espectral de la matriz (Subscript[B, 2])=(({{3, 0, -2}, {0, 3,    4}, {0, 0, -3}})). Determine también si la matriz (Subscript[B, 2]) es o no diagonalizable.

B2 =  {{3, 0, -2},  {0, 3, 4},  {0, 
          0, -3}}; 
Eigensystem[B2] // TableForm

El valor propio doble tiene una base de vectores propios formada por dos vectores, luego la matriz es diagonalizable.

El radio espectral de la matriz, como se observa a simple vista, es 3.
No hay un valor propio estrictamente dominante: es más, los tres son dominantes.

Matrices simétricas definidas positivas

Vamos a mostrar los menores asociados a una matriz (simétrica o no), como por ejemplo la siguiente (que en este caso es simétrica):

A = {{1, 2, 3}, {2, 5, 6}, {3, 6, 7}}; MatrixForm[A]

Los menores de órdenes 1, 2 y 3 se obtienen con las siguientes instrucciones:

Minors[A, 1]
Minors[A, 2]
Minors[A, 3]

Luego los menores principales de la matriz son:.

Table[Minors[A, i][[1, 1]], {i, 3}]

Una matriz simétrica de orden n es definida positiva si y sólo si sus menores principales (Subscript[m, k])(A), k = 1, 2,  ... , n, son todos positivos. Por tanto, la matriz anterior no es definida positiva.

Ejemplo 1

A = {{10, 1, 2, 3}, {1, 9, -1, 2}, {2, -1, 7, 3}, {3, 2, 3, 
          12}}; MatrixForm[A]

La matriz A es definida positiva ya que es simétrica y los menores principales son todos positivos:

n = Length[A]; 
    Table[Minors[A, i][[1, 1]], {i, n}]

Ejemplo 2

A = {{9, 5, 2, 2}, {5, 3, 1, 1}, {2, 1, 2, 2}, {2, 1, 2, 2}}; 
    MatrixForm[A]
n = Length[A]; 
    Table[Minors[A, i][[1, 1]], {i, n}]

La matriz A no es definida positiva ya que no todos los menores principales son positivos (uno de ellos es igual a cero).

Normas vectoriales y matriciales: órdenes directas para su cálculo

Las órdenes directas para el cálculo de las normas vectoriales y matriciales \(TraditionalForm\`\(\(||\)\(\(·\)\( || \_1\)\)\)\), \(\(\(||\)\(\(·\)\( || \_2\)\)\)\)  y  \(\(\(||\)\(\(·\)\( || \_∞\)\)\)\)  son:

v = {7. , 3, 1, -2}; 
{Norm[v, 1], Norm[v], Norm[v, 2], Norm[v, ∞]}

Observe que, en las instrucciones anteriores se tiene que Norm[v] = Norm[v,2]. Lo mismo ocurre con las normas matriciales:

A =  {{1. , 0, 1, 3},  {-1, 0, 1, 1},  {0, 3, 1, 0},  {7, 3, 
          1, -2}}; 
{Norm[A, 1], Norm[A], Norm[A, 2], Norm[A, ∞]}

El radio espectral de una matriz cuadrada A es el ínfimo de todas sus normas matriciales; de ahí que ρ(A) sea menor o igual que todas esas normas:

ρ[A]

Condicionamiento de una matriz

Las siguientes fórmulas definen el condicionamiento de una matriz inversible A según las normas matriciales \(TraditionalForm\`\(\(||\)\(\(·\)\( || \_1\)\)\)\), \(TraditionalForm\`\(\(||\)\(\(·\)\( || \_2\)\)\)\) y \(TraditionalForm\`\(\(||\)\(\(·\)\( || \_∞\)\)\)\).

condic1[A_] := 
      Norm[A, 1] Norm[Inverse[A], 1]
condic2[A_] := Norm[A] Norm[Inverse[A]]
condic∞[A_] := 
      Norm[A, ∞] Norm[Inverse[A], ∞]

Por ejemplo:

A =  {{1, 0, 1, 3},  {-1, 0, 1, 1},  {0, 3, 1, 0},  {7, 3, 
          1, -2}}; 
N[{condic1[A], condic2[A], condic∞[A]}]

El condicionamiento no depende del tamaño de los números de la matriz. Por ejemplo, si todos los elementos de la matriz se multiplican por el mismo número, el condicionamiento no varía.

N[{condic1[100  A], condic2[100  A], 
        condic∞[100  A]}]

Sin embargo, en el siguiente ejemplo los coeficientes de la matriz son bastante pequeños en relación a su condicionamiento:

A = {{2.00021, 2.9, 1}, {1, 0.94, 1}, {-0.01, 
          1, -1}}; 
{condic1[A], condic2[A], condic∞[A]}

Supongamos que se desea resolver un sistema n×n compatible determinado A x = b, en el que se ha producido un error en la obtención de los términos independientes, de  modo que el sistema que se resuelve realmente es \(TraditionalForm\`A\ x = \(b\^*\)\). Si el condicionamiento de la matriz A es muy grande (como ha ocurrido en el ejemplo anterior), esto podría influir negativamente en el error que se comete en la resolución del sistema, de manera que la solución x del sistema A x = b, diste mucho de la solución \(TraditionalForm\`\(x\^*\)\) del sistema \(TraditionalForm\`A\ x = \(b\^*\)\). Veamos algunos ejemplos:

Ejemplo

En este ejemplo el vector de los términos independientes es b = (2, -301, -501), pero por error se ha tomado \(TraditionalForm\`\(b\^*\) = \((2, \(-300\), \(-500\))\)\), en lugar de b. La matriz de coeficientes A es la siguiente:

A = {{-1, 2, 
          1}, {150, -300, -151}, {250, -500.1, -250.9}};

Las solución x del sistema A x = b y la solución \(TraditionalForm\`\(x\^*\)\) del sistema aproximado \(TraditionalForm\`A\ x = \(b\^*\)\) son en este caso las siguientes:

b = {2, -301, -501}; x = LinearSolve[A, b]

\(\(b\^*\) = {2, \(-300\), \(-500\)}; \(x\^*\) = 
      LinearSolve[A, \(b\^*\)]\)

Las soluciones de ambos sistemas son, por tanto, muy diferentes. Esto está relacionado con el  condicionamiento de la matriz A, que obtenemos para las tres normas consideradas:

{condic1[A], condic2[A], condic∞[A]}

El condicionamiento es bastante grande, lo que justifica que el error en sus términos independientes —que es más bien pequeño— haya afectado mucho a la solución del sistema. Este error relativo \(TraditionalForm\`r\_b = \(\(||\)\(\(b\^*\) - b\)\(||\)\)\/\(\(||\)\(b\)\(||\)\)\) lo obtenemos a continuación para las tres normas consideradas:

\({Norm[\(b\^*\) - b, 1]\/Norm[b, 1], Norm[\(b\^*\) - b]\/Norm[b], 
      Norm[\(b\^*\) - b, ∞]\/Norm[b, ∞]}\)
N[%]

Luego el error relativo en los términos independientes está en torno al 0.2  %. Sin embargo, el error relativo \(TraditionalForm\`r\_x = \(\(||\)\(\(x\^*\) - x\)\(||\)\)\/\(\(||\)\(x\)\(||\)\)\) cometido en el vector solución es bastante grande a simple vista; en concreto, según las normas consideradas, es igual a:

\({Norm[\(x\^*\) - x, 1]\/Norm[x, 1], Norm[\(x\^*\) - x]\/Norm[x], 
      Norm[\(x\^*\) - x, ∞]\/Norm[x, ∞]}\)

Así el error relativo (Subscript[r, x]) en la solución multiplica entre 670 y 1503 veces al error relativo (Subscript[r, b]) en el vector de términos independientes, según la norma que se use:

\({\(Norm[\(x\^*\) - x, 1]\/Norm[x, 1]\)÷\(Norm[\(b\^*\) - b, 
              1]\/Norm[b, 1]\), \(Norm[\(x\^*\) - x]\/Norm[
              x]\)÷\(Norm[\(b\^*\) - b]\/Norm[
              b]\), \(Norm[\(x\^*\) - x, ∞]\/Norm[
              x, ∞]\)÷\(Norm[\(b\^*\) - 
                b, ∞]\/Norm[b, ∞]\)}\)

El error relativo (Subscript[r, x]) sabíamos que estaba acotado como sigue:  ((Subscript[r, b])/(cond(A))) ≤ 
        (Subscript[r, x]) ≤ (Subscript[r, b]) cond(A). Para las normas consideradas, las cotas inferior y superior del error relativo, y el propio error (Subscript[r, x]), son:

\(TableForm[{{"Cota inferior", "Error relativo", "Cota superior"}, {Norm[\(b\^*\) - b, 1]\/\(condic1[A] Norm[b, 1]\), 
          Norm[\(x\^*\) - x, 1]\/Norm[x, 1], \(condic1[A] Norm[\(b\^*\) - b, 1]\)\/Norm[b, 1]}, {Norm[\(b\^*\) - b]\/\(condic2[A] Norm[b]\), 
          Norm[\(x\^*\) - x]\/Norm[x], \(condic2[A] Norm[\(b\^*\) - b]\)\/Norm[b]}, {Norm[\(b\^*\) - b, ∞]\/\(condic∞[A] Norm[b, ∞]\), 
          Norm[\(x\^*\) - x, ∞]\/Norm[x, ∞], \(condic∞[A] Norm[\(b\^*\) - b, ∞]\)\/Norm[b, ∞]}}]\)

Ejemplo (continuación del Ejemplo anterior)

Supongamos ahora que en el vector b hay un error distinto del anterior, y que el sistema se resuelve para \(TraditionalForm\`\(b\^*\) = \((2.01, \(-301.01\), \(-500.99\))\)\):

\(\(b\^*\) = {2.01, \(-301.01\), \(-500.99\)}; \)
\(\(x\^*\) = LinearSolve[A, \(b\^*\)]\)

La solución del sistema aparece tremendamente alterada (recuerde que la verdadera es x = (1, 1, 1)), a pesar de que el error de partida es ahora más pequeño: ahora se ha notado más que el condicionamiento de la matriz es superior a (10^6). El error relativo \(TraditionalForm\`r\_b = \(\(||\)\(\(b\^*\) - b\)\(||\)\)\/\(\(||\)\(b\)\(||\)\)\) lo obtenemos de nuevo para las tres normas consideradas:

\({Norm[\(b\^*\) - b, 1]\/Norm[b, 1], Norm[\(b\^*\) - b]\/Norm[b], 
      Norm[\(b\^*\) - b, ∞]\/Norm[b, ∞]}\)

Luego el error relativo en el vector de términos independientes está entre el 0.0019  % y  el 0.0038  %, según las normas consideradas. Sin embargo, el error relativo \(TraditionalForm\`r\_x = \(\(||\)\(\(x\^*\) - x\)\(||\)\)\/\(\(||\)\(x\)\(||\)\)\) cometido en el vector solución es enorme para las tres normas:

\({Norm[\(x\^*\) - x, 1]\/Norm[x, 1], Norm[\(x\^*\) - x]\/Norm[x], 
      Norm[\(x\^*\) - x, ∞]\/Norm[x, ∞]}\)

Así el error relativo (Subscript[r, x]) en la solución multiplica entre 340.000 y 1.250.000 veces al error relativo (Subscript[r, b]) en el vector de términos independientes, según la norma que se use:

\({\(Norm[\(x\^*\) - x, 1]\/Norm[x, 1]\)÷\(Norm[\(b\^*\) - b, 
              1]\/Norm[b, 1]\), \(Norm[\(x\^*\) - x]\/Norm[
              x]\)÷\(Norm[\(b\^*\) - b]\/Norm[
              b]\), \(Norm[\(x\^*\) - x, ∞]\/Norm[
              x, ∞]\)÷\(Norm[\(b\^*\) - 
                b, ∞]\/Norm[b, ∞]\)}\)

El error relativo (Subscript[r, x]) sabíamos que estaba acotado como sigue:  ((Subscript[r, b])/(cond(A))) ≤ 
        (Subscript[r, x]) ≤ (Subscript[r, b]) cond(A). Para las normas consideradas, las cotas inferior y superior del error relativo, y el propio error (Subscript[r, x]), son:

\(TableForm[{{"Cota inferior", "Error relativo", "Cota superior"}, {Norm[\(b\^*\) - b, 1]\/\(condic1[A] Norm[b, 1]\), 
          Norm[\(x\^*\) - x, 1]\/Norm[x, 1], \(condic1[A] Norm[\(b\^*\) - b, 1]\)\/Norm[b, 1]}, {Norm[\(b\^*\) - b]\/\(condic2[A] Norm[b]\), 
          Norm[\(x\^*\) - x]\/Norm[x], \(condic2[A] Norm[\(b\^*\) - b]\)\/Norm[b]}, {Norm[\(b\^*\) - b, ∞]\/\(condic∞[A] Norm[b, ∞]\), 
          Norm[\(x\^*\) - x, ∞]\/Norm[x, ∞], \(condic∞[A] Norm[\(b\^*\) - b, ∞]\)\/Norm[b, ∞]}}]\)

Resuelva los ejercicios 5 y 6 de la sección 3.7 de los apuntes de Gestión

Métodos de Jacobi y de Gauss-Seidel

El método de Jacobi es atribuido a Carl Jacobi (1804-1851) y la iteración de Gauss-Seidel a Johan Carl Friedrich Gauss (1777-1855) y Philipp von Seidel (1821-1896).

Recuerde que para un sistema de ecuaciones lineales n×n compatible determinado

(∑_(i=1)^n)(Subscript[a, ij])(Subscript[x, j])=(Subscript[b, i]),   i=1,..., n,

escrito en forma matricial como A x = b, el método de Jacobi consiste en efectuar la iteración

(x^((k + 1))) = (D^(-1)) (L + U) 
            (x^((k))) + (D^(-1)) b

o, equivalentemente,

(x^((k + 1)))=((Subscript[I, n])-(D^(-1))A)(x^((k)))+(D^(-1)) b

para k = 0, 1,  ... , a partir de un vector inicial cualquiera (x^((0))), donde A=D-L-U (D es la matriz diagonal extraída de la matriz A; -L coincide con la parte situada bajo la diagonal principal de A y -U es la parte situada sobre la diagonal principal de A; los demás elementos de las matrices L y U son ceros). El desarrollo de la expresión anterior da lugar a las siguientes igualdades:

((Subscript[x, i])^((k + 1))) = (1/((Subscript[a, ii]))) ((Subscript[b, i])-(∑_(j=1)^(i-1))(Subscript[a, ij]) ((Subscript[x, j])^((k)))- (∑_(j=i+1)^n)(Subscript[a, ij]) 
                      ((Subscript[x, j])^((k)))),       i=1,..., n.

Por otra parte, el método de Gauss-Seidel consiste en efectuar las iteraciones

(x^((k + 1)))=(((D - L))^(-1))U (x^((k)))+(((D - L))^(-1))b

o, equivalentemente,

(x^((k + 1)))=[(Subscript[I, n])-(((D - L))^(-1))A](x^((k)))+(((D - L))^(-1))b

para k = 0, 1,  ... , a partir de un vector inicial cualquiera (x^((0))). El desarrollo de la expresión anterior da lugar a las siguientes igualdades:

((Subscript[x, i])^((k + 1))) = (1/((Subscript[a, ii]))) ((Subscript[b, i])-(∑_(j=1)^(i-1))(Subscript[a, ij]) ((Subscript[x, j])^((k + 1)))- (∑_(j=i+1)^n)(Subscript[a, ij]) 
                      ((Subscript[x, j])^((k)))),       i=1,..., n.

Observe la pequeña diferencia de esta expresión con la del método de Jacobi.

En esta práctica, las matrices utilizadas en las iteraciones de los métodos de Jacobi y Gauss-Seidel las denotaremos así:

J=(D^(-1)) (L + U)=(Subscript[I, n])-(D^(-1))A        		J2=(D^(-1))b
GS=(((D - L))^(-1))U=(Subscript[I, n])-(((D-L))^(-1))A        	GS2=(((D - L))^(-1))b

A continuación introducimos un procedimiento para obtener las matrices D, L y U a partir de la matriz A. Luego lo aplicaremos para implementar los métodos de Jacobi y Gauss-Seidel.

Procedimiento para obtener las matrices D, L y U a partir de la matriz A

Partimos de una matriz cuadrada cualquiera. Por ejemplo:

A = {{5, 6, 7, 3, 4}, {1, 2, 3, 4, 5}, {2, 3, -4, 5, 0}, {3, 
          4, -5, 6, 1}, {4, 5, 1, 7, 5}}; MatrixForm[A]

Entonces la matriz D se obtiene así:

n = Length[A]; 
d = DiagonalMatrix[Table[A[[i, i]], {i, n}]]; 
    MatrixForm[d]

Para obtener las matrices L y U el procedimiento comienza por construir matrices triangulares genéricas:

Clear[a, i, j]; 
L = SparseArray[{i_, j_} /; i > j -> a[i, j], {n, n}]; 
    MatrixForm[L]
U = SparseArray[{i_, j_} /; i < j -> a[i, j], {n, n}]; 
    MatrixForm[U]

Y luego les asignamos los valores que les corresponden:

For[i = 1, i ≤ n, i++, 
      For[j = 1, j ≤ n, j++, 
        a[i, j] = -A[[i, j]]]]
MatrixForm[L]
MatrixForm[U]

En resumen:

Clear[a, i, j]; 
A = {{5, 6, 7, 3, 4}, {1, 2, 3, 4, 5}, {2, 3, -4, 5, 0}, {3, 
          4, -5, 6, 1}, {4, 5, 1, 7, 5}}; 
    MatrixForm[A]
n = Length[A]; 
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, n}]]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
U = 
      SparseArray[{i_, j_} /; i < j -> a[i, j], {n, 
          n}]; 
For[i = 1, i ≤ n, i++, 
      For[j = 1, j ≤ n, j++, 
        a[i, j] = -A[[i, j]]]]
{MatrixForm[d], MatrixForm[L], MatrixForm[U]}

Ejemplo: método de Jacobi

Consideramos el sistema cuya matriz de coeficientes y su vector de términos independientes son:

A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; 
    b = {16, 44, 59}; 
{MatrixForm[A], MatrixForm[b]}

Seguidamente se obtiene el orden n del sistema, y se establece un número de iteraciones máximo para el método (podrá variarse en cada ejercicio según convenga).

n = Length[A]; iter = 50;

Se introduce el vector inicial x de la iteración (por ejemplo, el vector nulo). Se asigna un valor inicial cualquiera al vector y, para que luego se puedan asignar nuevos valores a sus componentes.

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

Para mostrar las respectivas aproximaciones de la solución que proporcionan las iteraciones del método de Jacobi utilizamos la siguiente instrucción:

For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] x[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])]; x = y; Print[k, "  ", InputForm[x]]]

Alternativamente, podemos escribir la aplicación del método de Jacobi mediante la fórmula  matricial, obteniendo resultados muy parecidos pero no idénticos, por cuestiones de redondeo:

x = Table[0. , {i, n}];  
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, n}]]; 
J = IdentityMatrix[n] - Inverse[d] . A; 
J2 = Inverse[d] . b; 
For[k = 1, k ≤ iter, k++,  x = J . x + J2; 
      Print[k, "  ", InputForm[x]]]

Sabemos obtener una cota del error de aproximación siempre que encontremos una norma matricial tal que la norma de la matriz del método sea menor que 1: \(TraditionalForm\`‖J‖ < 1\). Para una matriz de gran tamaño usaríamos las normas \(TraditionalForm\`‖·‖\_1\) y \(TraditionalForm\`‖·‖\_∞\). Para nuestros ejercicios podremos usar también la norma espectral, pues son matrices de tamaño pequeño. Veamos los valores de las tres normas para el ejercicio que estamos  haciendo:

{Norm[J, 1], Norm[J], Norm[J, ∞]} // N

Luego las tres normas son utilizables en nuestro ejemplo. Optaremos en tal caso por la norma espectral.

Así, todas las instrucciones necesarias para aplicar el método de Jacobi, usando la primera fórmula introducida y acotando el error de cada aproximación serían las siguientes:

A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; 
    b = {16, 44, 59}; 
n = Length[A]; iter = 50; 
    x = Table[0. , {i, n}]; 
    y = x; 
d = 
      DiagonalMatrix[Table[A[[i, i]], {i, n}]]; 
    J = IdentityMatrix[n] - 
        Inverse[d] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[J, 
          1]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[J]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[J, ∞]]]\)
Print["n\.ba \
iterac", "               aproximación Jacobi",  "            \
         cota error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] x[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[J]/(1.  - Norm[J]) 
              Norm[y - x]; Print[k,  "   ", InputForm[y],  "   ", cotaerr];x = y ]

Con la cota del error podríamos establecer un criterio de parada del algoritmo, con lo que damos por concluido el programa. Por ejemplo, si se desea que el error sea menos que (10^(-6)):

A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; 
    b = {16, 44, 59}; 
n = Length[A]; iter = 50; 
    cota = (10^(-6)); x = Table[0. , {i, n}]; 
    y = x; 
d = 
      DiagonalMatrix[Table[A[[i, i]], {i, n}]]; 
    J = IdentityMatrix[n] - 
        Inverse[d] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[J, 
          1]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[J]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[J, ∞]]]\)
Print["n\.ba \
iterac", "               aproximación Jacobi",  "            \
         cota error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] x[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[J]/(1.  - Norm[J]) 
              Norm[y - x]; Print[k,  "   ", InputForm[y],  "   ", cotaerr];If[cotaerr < cota, Break[]];x = y ]

Comentario 1

Como el orden de la matriz de coeficientes es pequeño, se puede obtener sin problema la solución exacta del sistema, que es:

z = LinearSolve[A, b]

La usamos para ilustrar los errores cometidos por el método de Jacobi en este ejemplo. El error de la última aproximación es:

y - z

cuya norma euclídea es, efectivamente, inferior a la cota obtenida (6.997041677531698`*^-7):

Norm[y - z]

También se puede mostrar el error en cada paso del método de Jacobi:

A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; b = {16, 44, 59}; 
    z = LinearSolve[A, b]; 
n = Length[A]; 
    iter = 50; cota = (10^(-6)); x = Table[0. , {i, n}]; 
    y = x; 
d = 
      DiagonalMatrix[Table[A[[i, i]], {i, n}]]; 
    J = IdentityMatrix[n] - 
        Inverse[d] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[J, 
          1]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[J]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[J, ∞]]]\)
Print["n\.ba \
iterac", "               aproximación Jacobi",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] x[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[J]/(1.  - Norm[J]) 
              Norm[y - x]; r = Norm[y - z];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];If[cotaerr < cota, Break[]];x = y ]

Comentario 2

Otro modo de medir el error de aproximación es comprobando hasta qué punto Ay - b se aproxima al vector 0:

A . y - b

Y la norma euclídea de este error es:

Norm[A . y - b]

Comentario 3

Sabemos que el método de Jacobi converge si y sólo si los elementos diagonales de la matriz A son nulos y el radio espectral de la matriz del método J = (Subscript[I, 3]) - (D^(-1)) A es menor que 1 (cosa que sabemos desde que calculamos su norma matricial). Observaremos que cuanto menor sea el radio espectral, el método convergerá con más rapidez. Mostramos siempre también la matriz J del método para que se pueda comprobar si se ha obtenido bien en los ejercicios hechos sin el ordenador.

MatrixForm[J]
ρ[J] // N

Comentario 4

Si se utiliza otra norma distinta de la espectral los resultados son similares.
Para la norma  \(TraditionalForm\`‖·‖\_1\):

A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; 
    b = {16, 44, 59}; 
n = Length[A]; iter = 50; 
    cota = (10^(-6)); x = Table[0. , {i, n}]; 
    y = x; 
d = 
      DiagonalMatrix[Table[A[[i, i]], {i, n}]]; 
    J = IdentityMatrix[n] - 
        Inverse[d] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[J, 
          1]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[J]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[J, ∞]]]\)
Print["n\.ba \
iterac", "               aproximación Jacobi",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] x[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[J, 1]/(1.  - Norm[J, 1]) 
              Norm[y - x, 1]; r = Norm[y - z, 1];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];If[cotaerr < cota, Break[]];x = y ]

Para la norma  \(TraditionalForm\`‖·‖\_∞\).

A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; 
    b = {16, 44, 59}; 
n = Length[A]; iter = 50; 
    cota = (10^(-6)); x = Table[0. , {i, n}]; 
    y = x; 
d = 
      DiagonalMatrix[Table[A[[i, i]], {i, n}]]; 
    J = IdentityMatrix[n] - 
        Inverse[d] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[J, 
          1]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[J]], \*""   \[LeftDoubleBracketingBar]J(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[J, ∞]]]\)
Print["n\.ba \
iterac", "               aproximación Jacobi",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] x[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[J, ∞]/(1.  - 
                    Norm[J, ∞]) Norm[y - x, ∞]; r = Norm[y - z, ∞];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];If[cotaerr < cota, Break[]];x = y ]

Ejemplo: método de Gauss-Seidel

Todos los pasos son similares a los utilizados para el método de Jacobi.

En la siguiente instrucción: (1) se introducen la matriz de coeficientes y los téminos independientes del sistema; (2) se obtiene el orden del sistema; (3) se obtiene la matriz GS del método; (4) \(TraditionalForm\`obtenemos\ las\ normas\ \ ‖·‖\_1\), \(TraditionalForm\`‖·‖\_2\) y \(TraditionalForm\`‖·‖\_∞\) de GS.

Clear[a, i, j]; 
A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; 
    b = {16, 44, 59}; 
n = Length[A]; 
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, n}]]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
For[i = 1, i ≤ n, i++, 
      For[j = 1, j < i, j++, 
        a[i, j] = -A[[i, j]]]]
GS = 
      IdentityMatrix[n] - Inverse[d - L] . A; 
\(Print[\*""\[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[GS, 
          1]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[GS]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[GS, ∞]]]\)

Luego podemos asegurar que el método de Gauss-Seidel converge. A continuación escribimos el resto de instrucciones necesarias: (5) se establece el número de iteraciones máximo para el método (modificable); (6) se introduce el error máximo deseado (en norma); (7) se introduce el vector inicial x de la iteración; (8) se asigna un valor inicial cualquiera al vector y, para que luego se puedan asignar nuevos valores a sus componentes; (9) se muestran las iteraciones y la cota del error hasta que se consigue el error deseado.

iter = 50; cota = (10^(-6)); x = Table[0. , {i, n}]; 
    y = x; 
Print["n\.ba iterac", "         \
aproximación Gauss-Seidel",  "                     cota \
error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] y[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[GS]/(1.  - Norm[GS]) 
              Norm[y - x]; Print[k,  "   ", InputForm[y],  "   ", cotaerr];If[cotaerr < cota, Break[]];x = y ]

Alternativamente, podemos escribir la aplicación del método de Gauss-Seidel mediante la fórmula  matricial, obteniendo resultados muy parecidos pero no idénticos, por cuestiones de redondeo:

x = Table[0. , {i, n}];  
GS2 = Inverse[d - L] . b; 
Print["n\.ba iterac", "         aproximación Gauss-Seidel\
",  "                     cota error"]; 
For[k = 1, k ≤ iter, k++,  y = GS . x + GS2; 
      cotaerr = Norm[GS]/(1.  - Norm[GS]) Norm[y - x];  
      Print[k,  "   ", InputForm[y],  "   ", cotaerr]; 
      If[cotaerr < cota, Break[]]; x = y ]

Comentario 1

Como se conoce la solución exacta del sistema, podemos obtener el error de la última aproximación:

z = LinearSolve[A, b]; 
y - z

cuya norma euclídea es, efectivamente, inferior a la cota obtenida (7.05445×(10^(-7))):

Norm[y - z]

También se puede mostrar el error en cada paso del método de Gauss-Seidel:

Clear[a, i, 
      j]; 
A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 
          11}}; b = {16, 44, 59}; 
    z = LinearSolve[A, b]; 
n = Length[A]; 
    iter = 50; cota = (10^(-6)); x = Table[0. , {k, n}]; 
    y = x; 
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, 
            n}]]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
For[i = 1, 
      i ≤ n, i++, 
      For[j = 1, j < i, j++, 
        a[i, j] = -A[[i, 
              j]]]]
GS = 
      IdentityMatrix[n] - 
        Inverse[d - L] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[GS, 
          1]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[GS]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[GS, ∞]]]\)
Print["n\.ba \
iterac", "         aproximación Gauss-Seidel",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] y[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[GS]/(1.  - Norm[GS]) 
              Norm[y - x]; r = Norm[y - z];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];If[cotaerr < cota, Break[]];x = y ]

Comentario 2

Otro modo de medir el error de aproximación es comprobando hasta qué punto Ay - b se aproxima al vector 0:

A . y - b

Y la norma euclídea de este error es:

Norm[A . y - b]

Comentario 3

Sabemos que el método de Gauss-Seidel converge si y sólo si los elementos diagonales de la matriz A son nulos y el radio espectral de la matriz del método GS=(Subscript[I, 3])-(((D - L))^(-1))A es menor que 1 (cosa que sabemos desde que calculamos su norma matricial). Observaremos que cuanto menor sea el radio espectral, el método convergerá con más rapidez. Mostramos siempre también la matriz GS del método para que se pueda comprobar si se ha obtenido bien en los ejercicios hechos sin el ordenador.

MatrixForm[GS]
ρ[GS] // N

Comentario 4

Si se utiliza otra norma distinta de la espectral los resultados son similares.
Para la norma  \(TraditionalForm\`‖·‖\_1\):

Clear[a, i, j]; A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; 
    b = {16, 44, 59}; 
z = 
      LinearSolve[A, b]; 
n = Length[A]; 
    iter = 50; cota = (10^(-6)); x = Table[0. , {i, n}]; 
    y = x; 
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, 
            n}]]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
For[i = 1, 
      i ≤ n, i++, 
      For[j = 1, j < i, j++, 
        a[i, j] = -A[[i, 
              j]]]]
GS = 
      IdentityMatrix[n] - 
        Inverse[d - L] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[GS, 
          1]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[GS]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[GS, ∞]]]\)
Print["n\.ba \
iterac", "         aproximación Gauss-Seidel",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] y[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[GS, 1]/(1.  - Norm[GS, 1]) 
              Norm[y - x, 1]; r = Norm[y - z, 1];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];If[cotaerr < cota, Break[]];x = y ]

Para la norma  \(TraditionalForm\`‖·‖\_∞\).

Clear[a, i, j]; A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 11}}; 
    b = {16, 44, 59}; 
z = 
      LinearSolve[A, b]; 
n = Length[A]; 
    iter = 50; cota = (10^(-6)); x = Table[0. , {i, n}]; 
    y = x; 
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, 
            n}]]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
For[i = 1, 
      i ≤ n, i++, 
      For[j = 1, j < i, j++, 
        a[i, j] = -A[[i, 
              j]]]]
GS = 
      IdentityMatrix[n] - 
        Inverse[d - L] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[GS, 
          1]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[GS]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[GS, ∞]]]\)
Print["n\.ba \
iterac", "         aproximación Gauss-Seidel",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] y[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[GS, ∞]/(1.  - 
                    Norm[GS, ∞]) Norm[y - x, ∞]; r = Norm[y - z, ∞];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];If[cotaerr < cota, Break[]];x = y ]

Ejemplo

Efectúe 40 iteraciones de los métodos de Jacobi y de Gauss-Seidel para obtener, si es posible, una aproximación de la solución del sistema lineal Ax=b, donde A = (({{15, 17, 19}, {0.25, 0.33, 0.42}, {1, 1.2, 1.6}})) y b = (({{2120}, {43.4}, {164}})). Repita el ejercicio realizando 200 iteraciones. Muestre una cota del error de aproximación así como el verdadero error.

Método de Jacobi

En la siguiente instrucción se elimina la parte referente al error deseado, que no se especifica: lo que se determina es el número de iteraciones a realizar.

A = {{15, 17, 19}, {0.25, 0.33, 0.42}, {1, 1.2, 1.6}}; 
    b = {2120, 43.4, 164}; 
    z = LinearSolve[A, b]; 
n = Length[A]; 
    iter = 40; x = Table[0. , {i, n}]; 
    y = x; 
d = 
      DiagonalMatrix[Table[A[[i, i]], {i, n}]]; 
    J = IdentityMatrix[n] - 
        Inverse[d] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]·(Subscript[\[RightDoubleBracketingBar], 1])="", N[
        Norm[J, 1]], \*""   \[LeftDoubleBracketingBar]·(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[J]], \*""   \[LeftDoubleBracketingBar]·(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[J, ∞]]]\)
Print["n\.ba \
iterac", "               aproximación Jacobi",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] x[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[J]/(1.  - Norm[J]) 
              Norm[y - x]; r = Norm[y - z];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];x = y ]

Como se observa, el método de Jacobi no converge. No tiene sentido que hagamos los casos de 100 y 200 iteraciones. El motivo es que la matriz del método J = (Subscript[I, 3]) - (D^(-1)) A tiene radio espectral mayor que 1, como comprobamos a continuación. Las cotas del error salen negativas porque la norma espectral de J es mayor que 1. Tampoco nos vale elegir alguna de las otras dos normas, que también son mayores que 1.

MatrixForm[J]
ρ[J]

Método de Gauss-Seidel

En la siguiente instrucción se elimina la parte referente al error deseado, que no se especifica: lo que se determina es el número de iteraciones a realizar.

Clear[a, i, 
      j]; 
A = {{9, 1, 1}, {2, 10, 3}, {3, 4, 
          11}}; b = {16, 44, 59}; 
    z = LinearSolve[A, b]; 
n = Length[A]; 
    iter = 50; cota = (10^(-6)); x = Table[0. , {k, n}]; 
    y = x; 
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, 
            n}]]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
For[i = 1, 
      i ≤ n, i++, 
      For[j = 1, j < i, j++, 
        a[i, j] = -A[[i, 
              j]]]]
GS = 
      IdentityMatrix[n] - 
        Inverse[d - L] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[GS, 
          1]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[GS]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[GS, ∞]]]\)
Print["n\.ba \
iterac", "         aproximación Gauss-Seidel",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] y[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];cotaerr = Norm[GS]/(1.  - Norm[GS]) 
              Norm[y - x]; r = Norm[y - z];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];If[cotaerr < cota, Break[]];x = y ]

Clear[a, i, 
      j]; 
A = {{15, 17, 19}, {0.25, 0.33, 
          0.42}, {1, 1.2, 1.6}}; b = {2120, 43.4, 164}; 
    z = LinearSolve[A, b]; 
n = Length[A]; 
    iter = 40; x = Table[0. , {k, n}]; 
    y = x; 
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, 
            n}]]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
For[i = 1, 
      i ≤ n, i++, 
      For[j = 1, j < i, j++, 
        a[i, j] = -A[[i, 
              j]]]]
GS = 
      IdentityMatrix[n] - 
        Inverse[d - L] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[GS, 
          1]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[GS]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[GS, ∞]]]\)
Print["n\.ba \
iterac", "         aproximación Gauss-Seidel",  "            \
         cota error",  "    error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] y[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])]; cotaerr = Norm[GS]/(1.  - Norm[GS]) Norm[y - x];r = Norm[y - z];Print[k,  "   ", InputForm[y],  "   ", 
            cotaerr,  "   ", r];x = y ]

No se observa aún con claridad que el método de Gauss-Seidel converja. Estudiamos el radio espectral de la matriz del método G  S = (Subscript[I, 3]) - (((D - L))^(-1)) A:

MatrixForm[GS]
ρ[GS]

Luego al ser menor que 1 (pero no mucho menor) converge (pero lentamente).

Por otra parte, las cotas del error salieron negativas porque la norma espectral de GS es mayor que 1. Tampoco nos vale elegir alguna de las otras dos normas, que son aún mayores. Por tanto quitamos del programa la parte que se refiere a las cotas del error. Y realizamos las 200 iteraciones que pide el ejercicio:

Clear[a, i, 
      j]; 
A = {{15, 17, 19}, {0.25, 0.33, 
          0.42}, {1, 1.2, 1.6}}; b = {2120, 43.4, 164}; 
    z = LinearSolve[A, b]; 
n = Length[A]; 
    iter = 200; x = Table[0. , {k, n}]; 
    y = x; 
d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, 
            n}]]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
For[i = 1, 
      i ≤ n, i++, 
      For[j = 1, j < i, j++, 
        a[i, j] = -A[[i, 
              j]]]]
GS = 
      IdentityMatrix[n] - 
        Inverse[d - L] . 
          A; 
\(Print[\*""\[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 1])="", 
      N[Norm[GS, 
          1]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], 2])="", 
      N[Norm[GS]], \*""   \[LeftDoubleBracketingBar]GS(Subscript[\[RightDoubleBracketingBar], ∞])="", 
      N[Norm[GS, ∞]]]\)
Print["n\.ba \
iterac", "         aproximación Gauss-Seidel",   "          \
              error"]; 
For[k = 1,k ≤ iter,k++,For[i = 1,i ≤ n,i++,y[[i]]=(1/A)[[i, i]] (b[[i]]-(∑_(j = 1)^(i - 1))A[[i, 
                            j]] y[[j]]- (∑_(j = i + 1)^n) A[[i, 
                            j]]  x[[j]])];r = Norm[y - z]; Print[k,  "   ", InputForm[y], "   ", r];x = y ]

Resuelva el ejercicio 15 de la sección 3.7 de los apuntes de Gestión:

Ejercicio 15 a) de la sección 3.7 de los apuntes de Gestión

Obtenga una aproximación de la solución del sistema lineal Ax=b, donde

A = (({{3, 1, -1}, {0, 3, 1}, {-2, 1, -4}}))    y    b = (({{1}, {3}, {2.5}})),

efectuando el número de iteraciones de los métodos de Jacobi y de Gauss-Seidel que estime suficiente (por ejemplo, que la cota del error sea inferior a (10^(-6))).

Método de Jacobi

A = {{3, 1, -1}, {0, 3, 1}, {-2, 1, -4}}; 
    b = {1, 3, 2.5};

Método de Gauss-Seidel

Ejercicios 15 b) y c) de la sección 3.7 de los apuntes de Gestión

Obtenga una aproximación de la solución del sistema lineal Ax=b, donde

A = (({{3, 1, -1}, {0, 3, 1}, {-2, 1, -1}}))    y    b = (({{1}, {3}, {2.5}})),

efectuando el número de iteraciones de los métodos de Jacobi y de Gauss-Seidel que estime suficiente (por ejemplo, que la cota del error sea inferior a (10^(-6))). Compare los resultados con los del ejercicio 15 a) e intente justificar las diferencias en cuanto a la convergencia de los métodos.

Método de Jacobi

A = {{3, 1, -1}, {0, 3, 1}, {-2, 1, -1}}; 
    b = {1, 3, 2.5};

Método de Gauss-Seidel

Método de las Potencias

Recordemos que un valor propio dominante de una matriz cuadrada es aquél que tiene mayor módulo. Si es único y simple se dice que es estrictamente dominante, y lo representaremos por (Subscript[λ, 1]). 
Para una matriz cuadrada real, el valor propio estrictamente dominante (si existe) es siempre un número real. 
Si además la matriz es diagonalizable, el método de las potencias permite aproximar ese valor propio a partir de un vector inicial no nulo (x^((0))) (aunque no de cualquiera)

El método de las potencias consiste en efectuar las iteraciones

(x^((k + 1))) = A·(x^((k))),    k = 0,1,2,...

Si ((Subscript[x, i])^((k))) denota la i-ésima componente del vector (x^((k))),  entonces

((((Subscript[x, i])^((k + 1))))/(((Subscript[x, i])^((k)))))->(Subscript[λ, 1])

para algunos i = 1, 2,  ...   (al menos para uno).

Ejemplo

Partiendo de un vector inicial x, implemente el método de las potencias para obtener una aproximación del valor propio estrictamente dominante de la matriz

(({{3.9375`, -3.1875`, -5.8125`, 0.25`}, {1.875`, -3.375`, -5.625`, 2.5`}, {-2.3125`, 3.0625`, 8.9375`, -0.75`}, {-5.625`, 5.125`, 16.875`, 1.5`}})).

A = {{3.9375, -3.1875, -5.8125, 
          0.25}, {1.875, -3.375, -5.625, 2.5}, {-2.3125, 3.0625, 
          8.9375, -0.75}, {-5.625, 5.125, 16.875, 
          1.5}}; 
n = Length[A];

x = {1, 1, 1, 1}; 
For[i = 1, i ≤ 20, i++, y = A . x; 
      Print[Table[
          InputForm[
            y[[m]]/x[[m]]], {m, n}]]; x = y]

Como se observa, falla la convergencia del método de las potencias en la segunda coordenada (también si se parte de otros puntos iniciales). El motivo es que la segunda coordenada del vector propio asociado al valor propio dominante es nula.

TableForm[Eigensystem[A]]

TableForm[Chop[Eigensystem[A]]]

Ejemplo

Nota.- Para éste y los siguientes ejercicios, recuerde que el método de Jacobi converge si y sólo si el módulo de cualquier valor propio dominante de la matriz J=I-((D )^(-1))A es menor que 1; y que el método de Gauss-Seidel converge si y sólo si el módulo de cualquier valor propio dominante de la matriz GS=I-(((D - L) )^(-1))A es menor que 1.

Aplique el método de las potencias para comprobar si los métodos de Jacobi y de Gauss-Seidel, aplicados a cualquier sistema lineal con matriz de coeficientes  A = (({{5, 0, -1, 3}, {-1, 5, 2, -1}, {0, 2, 4, -1}, {3, -2, -1, -7}})),  son convergentes o no.

A = {{5, 0, -1, 3}, {-1, 5, 2, -1}, {0, 2, 
          4, -1}, {3, -2, -1, -7}}; 
n = Length[A];

Método de Jacobi

d = 
      DiagonalMatrix[
        Table[A[[i, i]], {i, n}]]; 
J = IdentityMatrix[n] - Inverse[d] . A; MatrixForm[J]

x = {1. , 2, 1, 3}; 
For[i = 1, i ≤ 40, i++, 
      y = J . x; 
Print[
        Table[InputForm[
            y[[m]]/x[[m]]], {m, n}]]; x = y]

Tras probar con varios vectores iniciales (además del anterior) no se consigue que haya convergencia en ninguna coordenada. El motivo en este caso es que la matriz J del método de Jacobi no tiene valor propio estrictamente dominante:

Eigenvalues[J] // N
Abs[%]

Pero lo anterior no significa que el método de Jacobi no converja: de hecho, al obtener sus valores propios, hemos visto que su radio espectral es menor que 1 (0.543032), luego debe converger para cualquier sistema Ax = b: compruébelo para algún ejemplo.

Método de Gauss-Seidel

Clear[a, i, j]; 
L = 
      SparseArray[{i_, j_} /; i > j -> a[i, j], {n, 
          n}]; 
For[i = 1, i ≤ n, i++, 
      For[j = 1, j ≤ n, j++, 
        a[i, j] = -A[[i, j]]]]; 
GS = IdentityMatrix[n] - Inverse[d - L] . A; MatrixForm[GS]

x = {1. , 1, 1, 1}; 
    For[i = 1, i ≤ 50, i++, 
      y = GS . x; 
Print[
        Table[InputForm[
            y[[m]]/x[[m]]], {m, n}]]; x = y]

Por tanto, el método de Gauss-Seidel converge, ya que el módulo del valor propio estrictamente dominante es menor que 1. Comprobación:

ρ[GS] // N
Eigenvalues[GS] // N

Resuelva los ejercicios 19, 21 y 24 de la sección 3.7 de los apuntes de Gestión:

Ejercicio 19 de la sección 3.7 de los apuntes de Gestión

Aplique el método de las potencias para obtener una aproximación del valor propio estrictamente dominante de la matriz A con al menos 6 dígitos de precisión.

A =(({{3, -3.5`, -5.8`, 0, 2}, {1, 0, -3, 2, 2}, {-2, 3, 5, -0.5`, 0}, {-5, 0, 5, 1, 1}, {2, 0, 5, 5, 0}})).

A = {{3, -3.5, -5.8, 0, 2}, {1, 0, -3, 2, 2}, {-2, 3, 
          5, -0.5, 0}, {-5, 0, 5, 1, 1}, {2, 0, 5, 5, 
          0}};

Ejercicio 21 de la sección 3.7 de los apuntes de Gestión

El método de Jacobi converge si y sólo si el módulo de cualquier valor propio dominante de la matriz J = I - (D^(-1)) A es menor que 1. Aplique el método de las potencias para comprobar si el método de Jacobi, aplicado a cualquier sistema lineal con matriz de coeficientes A = (({{5, 2, -3}, {2, 7, 2}, {0, 2, 1}})), es convergente o no. Repita el ejercicio para el método de Gauss-Seidel.

A = {{5, 2, -3}, {2, 7, 2}, {0, 2, 1}};

Método de Jacobi

Método de Gauss-Seidel