Práctica 3: Resolución numérica de SEL

Resolución mediante comandos directos de Mathematica

(*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 (*Idem,  pero dando el tiempo empleado*)

At = Transpose[A];                       
Ab = Transpose[AppendTo[At, b]];                  
Dimensions[Ab]                         
Ab = RowReduce[Ab];                   
x2 = Transpose[Ab][[n + 1]] ;      
Max[Abs[x1 - x2]](*comparamos resultados*)

Clear[A, b, n]
n = 12;                                       
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 de descomposición directa

n = 5;                                        
A = Table[E^(-1. /i*j), {i, n}, {j, n}];
{q, r} = QRDecomposition[A];
  
MatrixForm[
        q] // N
MatrixForm[
      r // N]
Transpose[q] . r - A // N    (*debería ser la matriz cero ...  ¿lo es?*)

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

Resolución por métodos iterativos

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]]
      ]  (*hemos programado una iteración;  
          ahora hacemos varias usando otro For*)

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]];    (*calculamos la distancia*)
                   Print["iteración n. = ", 
              i,  "    distancia = ", distancia];
                  If[distancia < tol, Break[]]   (*hacemos el test*)      ]

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

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

Ejercicio 6.
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