Práctica 2: Resolución numérica de SEL (28-nov-06)
1º. B Curso 2006/07

Ax=b     x=LinearSolve[A,b]    x=(A^(-1))b

EJERCICIOS PARA ENTREGAR:
Para la evaluación continua hay que realizar al menos los ejercicios 2 y 3 de esta práctica. El resto de ejercicios son voluntarios y se considerarán evaluación adicional. 

Forma y plazo de entrega.
En formáto electrónico, diskete, CD, DVD o memoria USB (que se recuperarán después de corregidos) o por email al profesor (jjmnieto ARROBA ugr.es) con la palabra TUTORIA en el asunto para evitar pérdidas. La corrección será en presencia del alumno, en horario de tutorías antes del día 22 de diciembre. Los ficheros enviados por email deben llamarse con el/los apellidos del alumno para evitar confusiones.

Los alumnos pueden decidir realizar estos ejercicios por parejas; en este caso ambos deben estar presentes para su evaluación.

Resolución mediante comandos directos de Mathematica

El comando LinearSolve

(*comenzamos con un ejemplo de entradas aleatorias*) 
n = 50;                                          (*número de ecuaciones e incógnitas*) 
A = Table[2*Random[], {i, n}, {j, n}];     
b = Table[Random[], {i, n}];                      
LinearSolve[A, b] ;   
x1 = %;                                    
LinearSolve[A, b] // Timing

Clear[A, b, n]
n = 6;                                       
A = Table[E^(-i*j), {i, n}, {j, n}];  
b = Table[1, {i, n}];                   
LinearSolve[A, b]; // Timing  
A = 
        N[A];
LinearSolve[A, 
          b];//Timing

Ejercicio 1.
Repetir el proceso anterior varias veces aumentando el n para comparar los tiempos
(Elimina algunas salidas que seguramente se salgan de la pantalla)
(En un PIV 1500Mhz, 256Mb para n=2000  el tiempo es 42 segundos)

Ejercicio 2.
Resuelve numéricamente Ax = b  para la matriz

 A=(({{ 2, -1,  0, ..., ...,  0}, {-1,   2, -1,  0, ...,  0}, { 0, -1,  2, -1, ...,  :}, { :,  0, -1, ..., -1,  0}, { :,  :,  :, -1,  2, -1}, { 0,  0, ...,  0,  -1,  2}})) de dimensión (50 x 50) y b=(({{1}, {2}, {3}, {:}, {:}, {50}})).

Métodos directos clásicos

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 (sin normalizar)

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

Sistema Triangular Superior (*Sustituciones regresivas*)

U = 
        Table[If[i ≤ j, Random[], 0], {i, 7}, {j, 
            7}];
b = Table[Random[], {i, 7}];

n = Length[U];
x = Table[0, {i, n}];
For[i = n, i ≥ 1, i = i - 1, 
      x[[i]] = (b[[i]] - 
              Sum[U[[i, j]]*x[[j]], {j, i + 1, n}])/
          U[[i, i]]]   
Print["errores = ", U . x - b]

Ejercicio 3.
Reescribe el programa para un sistema triangular inferior usando sustituciones progresivas.

Métodos de factorización directa

Ejemplo de resolución por factorización

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 QR

n = 5;
A = 
        Table[E^(-1. /i*j), {i, n}, {j, 
            n}];
{q, r} = 
        QRDecomposition[A];
  
MatrixForm[
      q]
MatrixForm[r]
Transpose[q] . r-A

Ejercicio 4.
Verifica que Q es ortogonal en el ejemplo propuesto.

Descomposición LU (de Doolittle)   (*sin usar LUDecomposition*)

A = Table[Random[], {i, 7}, {j, 7}];

n = Length[A];
L = IdentityMatrix[n];     
U = 
        L;
For[i = 1,i ≤ n,i++,
      U[[i, i]] = (A[[i, i]] - 
                Sum[L[[i, k]]*U[[k, i]], {k, 1, i - 1}]);  
      If[U[[i, i]] == 0, 
            Print["No admite descomp. LU"];  Break[]];  
       For[j = i + 1, 
            j ≤ 
              n, j++, 
                  L[[j, i]] = (A[[j, i]] - 
                    Sum[L[[j, k]]*U[[k, i]], {k, 1, i - 1}])/
                U[[i, 
                    i]]; 
                  U[[i, j]] = (A[[i, j]] - 
                  Sum[L[[i, k]]*U[[k, j]], {k, 1, 
                      i - 1}]);
                 ]
      ]
Print["L = ", 
      MatrixForm[L]]
Print["U = ", 
      MatrixForm[L]]
Print["error = ", 
      MatrixForm[A - L . U]]

Ejercicio 5.
Reescribe el programa para el caso en que A sea simétrica la descomposición obtenida sea la de Cholesky.

b = Table[Random[], {i, n}];   
y = LinearSolve[L, b];
x = LinearSolve[U, y];
Print["residuo = ", 
          A . x - b]

Resolución por métodos iterativos

Descomposición  previa     A=M-L-U, estudio matricial

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]]

Convergencia, Jacobi

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

Convergencia, Gauss-seidel

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

Convergencia, relajación w=0.8

w = 0.8;
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}];

Gauss-Seidel

Clear[A, b, n, k, i, 
          j];
n = 50;
A = Table[Random[], {i, n}, {j, n}];   
Do[A[[i, i]] = n, {i, n}]   
b = Table[Random[], {i, n}];

x = Table[0, {i, n}];    
 (*recordamos ((Subscript[x, k])^((n + 1)))= ((Subscript[b, k])-(∑_(j = 1)^(k - 1))(Subscript[a, kj])((Subscript[x, j])^((n + 1)))-(∑_(j = k + 1)^n)(Subscript[a, kj])((Subscript[x, j])^((n))))/(Subscript[a, kk  ])para G-S*)
For[k = 1, 
          k <= n, k++, 
        x[[k]] = (b[[k]] - Sum[A[[k, j]]*x[[j]], {j, 1, k - 1}] - 
                  Sum[A[[k, j]]*x[[j]], {j, k + 1, n}])/
              A[[k, k]]
      ]

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

    For[i = 1, 
      i <= 10, i++, 
                  For[k = 1, 
        k ≤ 
          n, k++, 
               x[[k]] = (b[[k]] - Sum[A[[k, j]]*x[[j]], {j, 1, k - 1}] - 
                Sum[A[[k, j]]*x[[j]], {j, k + 1, n}])/
            A[[k, k]]
                ]       
      ]

A . x - b

x = Table[0, {i, n}];
tol = 10^(-6);
For[i = 1,i <= 100,i++,
                  y = x;                 
                  For[k = 1, 
                k ≤ 
                  n, k++, 
               x[[k]] = (b[[k]] - Sum[A[[k, j]]*x[[j]], {j, 1, k - 1}] - 
                        Sum[A[[k, j]]*x[[j]], {j, k + 1, n}])/
                    A[[k, 
                        k]]
                      ];
                  distancia = Max[Abs[x - y]];      
                   Print["iteración n. = ", 
                i,  "    distancia = ", distancia];
                  If[distancia < tol, Break[]]           ]

Ejercicio 6.
Modifica el programa para que en cada iteración muestre la solución aproximada y aplícalo con una matriz 5x5.

Ejercicio 7.
Usa una matriz aleatoria cualquiera y verifica que G-S no tiene por qué converger

Ejercicio 8.
Programa tú el método de Jacobi y el de Relajación para un peso genérico ω.
Para la matriz del ejercicio 2, compara la velocidad