# Sesión 4: Resolución de sistemas de 
ecuaciones lineales: métodos directos.

## Manejando listas y matrices

Definiendo matrices:

kill(all) $

A : matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]);

A[1,2];

B : genmatrix(b,3,3);

b[i,j]:=i-j;

B;

B : genmatrix(b,3,3);

Operaciones básicas con matrices:

B.A;

B*A;

AA : matrix([1, 2, 1], [2, 5, 0], [-3, 2, 6]);
cc1 : matrix([4, -2, 1]);
cc2 : matrix([4],[-2],[1]);

AA.cc1;

AA.cc2;

cc1.AA;

cc2.AA;

A^2;

A^^2;

AA^^-1;

transpose(A);

determinant(A);

Manipulación de filas, columnas y submatrices:

submatrix(1,A,1);

submatrix(1,A,2,3);

row(A,2);

col(A,3);

addcol(A,[1,2,3]);

A;

addrow(A,[1,2,3]);

length(cc1);

length(cc2);

## Método de Gauss

cc : matrix([4, -2, 1]) $

### Gauss evitando ceros

La expresión ``A[i,k] # 0´´ sólo tiene sentido si operamos con valores exactos. Por tanto, usaremos 
``abs(A[i,k]) > tol´´ siendo tol una tolerancia predefinida.

kill(A,b,c);

A : copymatrix(AA) $
c : copymatrix(cc) $
b : A.c $ 
n : length(A) $
tol : 10^(-15) $

print("Sistema original:") $
print(A,".x=",b) $

Algoritmo básico de triangulación sin control de los pivotes:

for k:1 while k<n do(
    for i:k+1 while i<=n do (
        m : -A[i,k]/A[k,k],
        A[i] : A[i] + m*A[k],
        b[i] : b[i] + m*b[k]
    ),
    print(A,".x=",b)
) $

EJERCICIO: buscar una situación en la que este método falle (aparezcan pivotes nulos)

Resolución del sitema triangular resultante:

x: genmatrix(lambda([i,j], 0), n, 1) $
/* Otra opción: x: zeromatrix (n, 1) */
for i:n thru 1 step -1 do 
x[i] : (b[i] - sum(A[i,j]*x[j],j,i+1,n))/A[i,i] $
print("x=",x) $

Comprobemos que la solución obtenida es la correcta.

A.x - b;

Algoritmo de triangulación con control de los pivotes:

A : copymatrix(AA) $
c : copymatrix(cc) $
b : A.c $ 
n : length(A) $
tol : 10^(-15) $

print("Método de Gauss:") $
for k:1 while k < n do (
    filabuena : 0,
    for i:k while i <= n do (
        if abs(A[i,k]) > tol then (
            filabuena : i,
            i : n+1
        )
    ),
    if filabuena = 0 then return(print("ERROR: Matriz singular")),
    if filabuena > k then (
        print("cero en fila ", k, "; cambio con fila ", filabuena),
        tmp : A[k],
        A[k] : A[filabuena],
        A[filabuena] : tmp,
        tmp : b[k],
        b[k] : b[filabuena],
        b[filabuena] : tmp
    ),
    for i:k+1 while i<=n do (
        m : -A[i,k]/A[k,k],
        A[i] : A[i] + m*A[k],
        b[i] : b[i] + m*b[k]
    ),
    print(A,".x=",b)
) $

EJERCICIO: Resuelve, usando el Método de Gauss, el sistema de ecuaciones lineales A.X=b, donde A y b se describen a continuación.
A es una matriz 8x8, cuyas entradas son:
- en la diagonal principal  a[i,i] = 1/i,
- por encima de la diagonal a[i,j] = j/i,
- por debajo de la diagonal a[i,j] = i/j.
El vector de términos independientes tiene coordenadas b[i] = cos(i).

## Descomposición LU

### Algoritmo de factorización de Doolittle

A : copymatrix(AA) $
c : copymatrix(cc) $
b : A.c $
[n,n2] : matrix_size (A) $ 
L : ident(n) $
U : zeromatrix(n,n) $

for k:1 thru n do (
    U[k,k] : A[k,k]-sum(L[k,r]*U[r,k],r,1,k-1),
    for j:k+1 thru n do (
        U[k,j] : A[k,j]-sum(L[k,r]*U[r,j],r,1,k-1),
        L[j,k] : (A[j,k]-sum(L[j,r]*U[r,k],r,1,k-1))/U[k,k]
    )
) $
LD : copymatrix(L) $
UD : copymatrix(U) $
print(A, "=", LD, UD, "=", LD.UD) $