Práctica 2
Resolución de sistemas de ecuaciones lineales (SEL)
Ax=b     x=LinearSolve[A,b]

Métodos Directos

Método de Gauss y Gauss-Jordan

Sistema triangular superior, Método de sustituciones regresivas

U = {{1, 2, -1, 3, -4}, {0, 6, 0, 0, 1}, {0, 0, 2, 2, 0}, {0, 
            0, 0, 1, -4}, {0, 0, 0, 0, 5}};
MatrixForm[U]
b = {25, 24, -15, 40, 61};

n = Length[U];
x = b;
For[i = n, i ≥ 1, i--, 
      x[[i]] = (b[[i]] - 
              Sum[U[[i, k]]*x[[k]], {k, i + 1, n}])/
          U[[i, i]]]
x
LinearSolve[U, b]
U . x == b

Método de Gauss

A = {{1, 2, 3, 4}, {2, 20, 18, 16}, {3, 18, 19, 21}, {4, 16, 21, 
            33}};
A = N[A];
b = {4, 6, 8, 2};

n = Length[A];
For[i = 1, i ≤ n, i++, 
      For[j = i + 1, j ≤ n, j++, 
        pivote = A[[j, i]]/A[[i, i]]; 
A[[j]] = 
          A[[j]] - pivote*A[[i]]; 
b[[j]] = 
          b[[j]] - pivote*b[[i]]]]
MatrixForm[A]
MatrixForm[b]

Método de Gauss-Jordan

A = {{1, 2, 3, 4}, {2, 20, 18, 16}, {3, 18, 19, 21}, {4, 16, 21, 
            33}};
A = N[A];
b = {4, 6, 8, 2};
LinearSolve[A, b]

n = Length[A];
For[i = 1, i ≤ n, i++, 
      For[j = i + 1, j ≤ n, j++, 
        pivote = -A[[j, i]]/
            A[[i, i]]; 
A[[j]] = 
          A[[j]] + pivote*A[[i]]; 
b[[j]] = 
          b[[j]] + pivote*b[[i]]]; 
For[j = 1, 
        j ≤ i - 1, j++, 
        pivote = -A[[j, i]]/
            A[[i, i]]; 
A[[j]] = 
          A[[j]] + pivote*A[[i]]; 
b[[j]] = 
          b[[j]] + 
            pivote*b[[i]]]
]
MatrixForm[A]

x = b;
For[i = 1, i ≤ n, i++, 
        x[[i]] = b[[i]]/A[[i, i]]];
x

Condicionamiento

A = {{10, 7, 8, 7}, {7, 5, 6, 5}, {8, 6, 10, 9}, {7, 5, 9, 
            10}};
A = N[A];
radioespectraldeA = Max[Abs[Eigenvalues[A]]]
radioespectraldesuinversa = 
      Max[Abs[Eigenvalues[Inverse[A]]]]
condicionamientooptimo = 
      radioespectraldeA*radioespectraldesuinversa

Métodos de descomposición directa

Descomposición LU

A = {{1, 2, 4}, {2, 9, 7}, {-1, 
            8, -2}};
L = {{1, 0, 0}, {2, 1, 0}, {-1, 2, 
            4}};
U = {{1, 2, 4}, {0, 5, -1}, {0, 0, 
            1}};
L . U == A

b = {10, 9, -2};
y = LinearSolve[L, b];
x = LinearSolve[U, y];
A . x == b

Descomposición LU de Doolittle

n = Dimensions[A][[1]];
U = 0*A;
L = IdentityMatrix[n];
For[k = 1, 
        k <= n, k++, 
U[[k, k]] = 
          A[[k, k]] - Sum[L[[k, t]]*U[[t, k]], {t, 1, k - 1}];  
        For[j = k + 1, 
          j <= n, j++, 
U[[k, j]] = 
            A[[k, j]] - 
              Sum[L[[k, t]]*U[[t, j]], {t, 1, 
                  k - 1}]]; 
     For[i = k + 1, 
          i <= n, i++, 
L[[i, 
                k]] = (A[[i, k]] - 
                  Sum[L[[i, t]]*U[[t, k]], {t, 1, k - 1}])/
              U[[k, k]]]
  ];
MatrixForm[L]
MatrixForm[U]
A == L . U

Sobre el comando LUDecomposition

LUDecomposition[A]

desc = LUDecomposition[A][[1]];
auxL = {{0, 0, 0}, {1, 0, 0}, {1, 1, 
            0}};
auxU = {{1, 1, 1}, {0, 1, 1}, {0, 0, 
            1}};
U = desc*auxU;
L = desc*auxL + IdentityMatrix[n];
A == L . U

Descomposición LL^t de Cholesky

A = {{16, -4, 0, 4, 0, -4}, {-4, 2, 2, 0, 0, 3}, {0, 2, 
            13, 5, 0, 13}, {4, 0, 5, 4, -3, 0}, {0, 0, 0, -3, 45, 
            12}, {-4, 3, 13, 0, 12, 31}};
A = A // N;
A == Transpose[A]     
Min[
        Eigenvalues[
          A]]>0

n = Length[A];
L = 0*A;
For[k = 1, k ≤ n, k++, 
      L[[k, k]] = 
        Sqrt[A[[k, k]] - 
            Sum[L[[k, j]]^2, {j, 1, k - 1}]]; 
For[
        i = k + 1, i ≤ n, i++, 
        L[[i, k]] = (A[[i, k]] - 
                Sum[L[[i, j]]*L[[k, j]], {j, 1, k - 1}])/
            L[[k, k]]]]
MatrixForm[L]
A == L . Transpose[L]

Descomposición QR

A = {{1, 2, 4}, {2, 9, 7}, {-1, 
            8, -2}};
{Q, R} = QRDecomposition[A];
Q = Transpose[Q];
A == Q . R

Métodos Iterativos

Descomposición  previa     A=M-L-U

A = {{2, 1, -1}, {-2, 4, 3}, {6, 5, 
            9}};
A = N[A];
b = {14/3, 67, 24};

n = Length[A];
auxL = 
        Table[If[i > j, 1, 0], {i, n}, {j, n}];
auxU = 
        Table[If[i < j, 1, 0], {i, n}, {j, n}];
L = (-A)*auxL;                    
U = (-A)*auxU;
d = A + L + U;
A == d - L - U
Print[MatrixForm[A], " = ", MatrixForm[d], " - ", 
      MatrixForm[L], " - ", MatrixForm[U]]

Jacobi

A = N[A];
n = Length[A];
tolerancia = 10^(-5);
error = 10*tolerancia;
maxiter = 100;
paso = 0;
x = Table[0, {n}];
y = x;
While[
      paso < maxiter  &&  error > tolerancia, 
Do[
        y[[i]] = (b[[i]] - 
                Sum[A[[i, j]]*x[[j]], {j, 1, i - 1}] - 
                Sum[A[[i, j]]*x[[j]], {j, i + 1, n}])/
            A[[i, i]], {i, n}]; 
error = 
        Sqrt[(x - y) . (x - y)]; 
x = 
        y; 
paso = paso + 1
]
Print["Solución aproximada: ", x];
Print["Iteraciones realizadas: ", 
        paso, ",      Tope de iteraciones: ", maxiter];
Print["Error permitido: ",  
        tolerancia, ",        Error cometido: ", error];
Print["Residuo (= |Ax-b|): ", 
        Sqrt[(A . x - b) . (A . x - 
                b)]];

Convergencia, estudio matricial

B = Inverse[d] . (L + U);
Max[Abs[Eigenvalues[B]]]

Gauss-Seidel

A = N[A];
n = Length[A];
iteraciones = 50;
x = Table[0, {n}];
For[paso = 0, 
      paso ≤ iteraciones, paso++, 
Do[
        x[[i]] = (b[[i]] - 
                Sum[A[[i, j]]*x[[j]], {j, 1, i - 1}] - 
                Sum[A[[i, j]]*x[[j]], {j, i + 1, n}])/
            A[[i, i]], {i, 
          n}]
]
Print["Solución aproximada: ", 
        x];
Print["Iteraciones realizadas: ", 
        iteraciones];
Print["Residuo (= |Ax-b|): ", 
        Sqrt[(A . x - b) . (A . x - b)]];

Convergencia, estudio matricial

B = Inverse[d - L] . U;
Max[Abs[Eigenvalues[B]]]

Relajación

A = N[A];
w = 0.8;
n = Length[A];
iteraciones = 50;
x = Table[0, {n}];
For[paso = 0, paso ≤ iteraciones, paso++, 
      Do[x[[i]] = 
          w*(b[[i]] - Sum[A[[i, j]]*x[[j]], {j, 1, i - 1}] - 
                    Sum[A[[i, j]]*x[[j]], {j, i + 1, n}])/
                A[[i, i]] + (1 - w)*x[[i]], {i, 
          n}]]
Print["Solución aproximada: ", 
        x];
Print["Iteraciones realizadas: ", 
        iteraciones];
Print["Residuo (= |Ax-b|): ", 
        Sqrt[(A . x - b) . (A . x - b)]];

Convergencia,estudio matricial

B = 
        Inverse[d/w - L] . ((1 - w)*d/w + 
              U);
Max[Abs[Eigenvalues[B]]]

f[w_] := 
      Max[Abs[Eigenvalues[
            Inverse[d/w - L] . ((1 - w)*d/w + 
                  U)]]]
Plot[f[w], {w, 0, 2}];

Ejercicios

1.- Programa la resolución de un sistema triangular inferior.
2.- Calcula el condicionamiento de una matriz usando la norma 1.
3.- Programa la descomposición LU de una matriz con la elección de Crout.
4.- Reescribe el programa para que detecte si una matriz NO puede ser descompuesta LU.
5.- Reescribe el programa de Cholesky para que detecte si la matriz no admite descomposición.
6.- Modifica el programa de Gauss-Seidel para que se detenga cuando la distancia eucl´idea entre 2 iteraciones sea inferior a 10^(-5).
7.- ¿Cómo detectar´ias numéricamente que un método no converge?
8.- Repite el ejercicio 6 para Relajación.