Métodos Práctica 4

// CONSULTAR EJERCICIO 2: SI ESTÁ BIEN Y SI TENGO QUE PIVOTEAR. Fuck me
// CONSULTAR EJERCICIO 3: ¿CÓMO DEVOLVER LAS MATRICES? Asignarle un par de valor de retorno. CAMBIAR TODAS LAS FUNCIONES. Fuck me ^ 2


// Ejercicio 1
function x = sustitucionRegresiva(A)
    dimension = size(A)(1)
    sol(dimension) = A(dimension, dimension+1) / A(dimension,dimension)
    for i = dimension-1:-1:1
        sol(i) = (A(i,dimension+1) - A(i, i+1:dimension) * sol(i+1:dimension)) / A(i,i)
    end
    x = sol
endfunction

function matrizExpandida = eliminacionGaussSinPivoteo(A, b)
    dimension = size(b)(1)
    for i = 1 : dimension-1
        for j = i+1 : dimension
            coeficiente = A(j,i)/A(i,i)
            A(j,i:dimension) = A(j,i:dimension) - A(i,i:dimension) * coeficiente
            b(j) = b(j) - b(i) * coeficiente
        end
    end
    matrizExpandida = [A,b]
endfunction

function matrizExpandida = eliminacionGaussConPivoteo(A, b)
    dimension = size(b)(1)
    for i = 1 : dimension-1
        if(A(i,i) == 0)
            for j = i+1:dimension
                if(A(j,i) ~= 0)
                    VectorAux = A(j, :)
                    A(j, :) = A(i, :)
                    A(i, :) = VectorAux

                    Auxiliar = b(j)
                    b(j) = b(i)
                    b(i) = Auxiliar
                    break
                end;
            end;
        end;
        for j = i+1:dimension
            coeficiente = A(j,i)/A(i,i)
            A(j,i:dimension) = A(j,i:dimension) - A(i,i:dimension) * coeficiente
            b(j) = b(j) - b(i) * coeficiente
        end;
    end;
    matrizExpandida = [A,b]
endfunction

function solucion = metodoGauss(A,b)
    solucion = sustitucionRegresiva(eliminacionGaussConPivoteo(A,b))
endfunction


--> A = [1 1 0 3; 2 1 -1 1; 3 -1 -1 2; -1 2 3 -1];
--> b = [4; 1; -3; 4];
--> metodoGauss(A,b)
 ans  =
  -1.
   2.
   0.
   1.


--> A = [1 -1 2 -1; 2 -2 3 -3; 1 1 1 0; 1 -1 4 3];
--> b = [-8; -20; -2; 4];
--> metodoGauss(A,b)
 ans  =
  -7.
   3.
   2.
   2.


--> A = [1 1 0 4; 2 1 -1 1; 4 -1 -2 2; 3 -1 -1 2];
--> b = [2; 1; 0; -3];
--> metodoGauss(A,b)
 ans  =
  -4.
   0.6666667
  -7.
   1.3333333



// Ejercicio 2
function x = sustitucionEspecial(A)
    dimension = size(A)(1)
    sol(dimension) = A(dimension, dimension+1) / A(dimension,dimension)
    for i = dimension-1:-1:1
        sol(i) = (A(i,dimension+1) - A(i, i+1) * sol(i+1)) / A(i,i)
    end
    x = sol
endfunction

//Supongo que no tengo que pivotear.
function matrizExpandida = eliminacionGaussEspecial(A, b)
    dimension = size(b)(1)
    for i = 1 : dimension-1
        coeficiente = A(i+1,i)/A(i,i)
        A(j,i:i+1) = A(j,i:i+1) - A(i,i:i+1) * coeficiente
        b(j) = b(j) - b(i) * coeficiente
    end;
    matrizExpandida = [A,b]
endfunction

function solucion = GaussEspecial(A,b)
    solucion = sustitucionEspecial(eliminacionGaussEspecial(A,b))
endfunction




// Ejercicio 3
//Supongo que no tengo que pivotear. A es rectangular m x n.
function matrices = factorizacionLU(A)
    m = size(A)(1)
    n = size(A)(2)
    L = eye(m,m)
        for i = 1 : m
            for j = i+1 : m
                coeficiente = A(j,i)/A(i,i)
                L(j,i) = coeficiente
                A(j,i:n) = A(j,i:n) - A(i,i:n) * coeficiente
            end
        end
    U = A
    matrices = [L, U]
endfunction


--> factorizacionLU(A1)
 ans  =
   1.   0.   0.   0.   1.   1.   0.
   2.   1.   0.   0.   0.  -1.  -1.
   4.   5.   1.   0.   0.   0.   3.
   3.   4.   1.   1.   0.   0.   0.
--> ans(:, 1:4) * ans(:, 5:7)
 ans  =
   1.   1.   0.
   2.   1.  -1.
   4.  -1.  -2.
   3.  -1.  -1.
--> A1
 A1  =
   1.   1.   0.
   2.   1.  -1.
   4.  -1.  -2.
   3.  -1.  -1.


--> factorizacionLU(A2)
 ans  =
   1.     0.   4.  -1.    -2.    2.
   0.75   1.   0.  -0.25   0.5   0.5
--> ans(:, 1:2) * ans(:, 3:6)
 ans  =
   4.  -1.  -2.   2.
   3.  -1.  -1.   2.
--> A2
 A2  =
   4.  -1.  -2.   2.
   3.  -1.  -1.   2.


--> ApartadoA = [1.012 -2.132 3.104; -2.132 4.096 -7.013; 3.104 -7.013 0.014]
 ApartadoA  =
   1.012  -2.132   3.104
  -2.132   4.096  -7.013
   3.104  -7.013   0.014
--> factorizacionLU(ApartadoA)
 ans  =
   1.          0.          0.   1.012  -2.132       3.104
  -2.1067194   1.          0.   0.     -0.3955257  -0.4737431
   3.0671937   1.1977555   1.   0.      0.         -8.9391408
--> ans(:, 1:3) * ans(:, 4:6)
 ans  =
   1.012  -2.132   3.104
  -2.132   4.096  -7.013
   3.104  -7.013   0.014


--> ApartadoB = [2.1756 4.0231 -2.1732 5.1967; -4.0231 6 0 1.1973; -1 5.2107 1.1111 0; 6.0235 7 0 4.1561]
 ApartadoB =
   2.1756   4.0231  -2.1732   5.1967
  -4.0231   6.       0.       1.1973
  -1.       5.2107   1.1111   0.
   6.0235   7.       0.       4.1561
--> factorizacionLU(ApartadoB)
 ans  =
         column 1 to 6
   1.          0.          0.          0.   2.1756   4.0231
  -1.849191    1.          0.          0.   0.       13.43948
  -0.4596433   0.5253098   1.          0.   0.       8.882D-16
   2.7686615  -0.3079436   2.1497102   1.   0.       0.

         column 7 to 8
  -2.1732      5.1967
  -4.0186619   10.806991
   2.2232457  -3.2883901
   0.          0.1652261
--> ans(:, 1:4) * ans(:, 5:8)
 ans  =

   2.1756   4.0231  -2.1732   5.1967
  -4.0231   6.       0.       1.1973
  -1.       5.2107   1.1111   0.
   6.0235   7.       0.       4.1561



// Ejercicio 4
// Para la primer factorización, tenemos: P34  E42(-8)  E32(5)  E41(-6)  E31(-5)  E21(-4)  A  = U
// Es decir, antes de hacer la factorización, debería permutar las filas 3 y 4 de A.

//Siendo U cuadrada.
function x = sustitucionHaciaArriba(U, b)
    dimension = size(U)(1)
    sol(dimension) = b(dimension) / U(dimension,dimension)
    for i = dimension-1:-1:1
        sol(i) = (b(i) - U(i, i+1:dimension) * sol(i+1:dimension)) / U(i,i)
    end
    x = sol
endfunction


function x = sustitucionHaciaAbajo(L, b)
    dimension = size(L)(1)
    sol(1) = b(1) / L(1,1)
    for i = 2:dimension
        sol(i) = (b(i) - L(i, 1:i-1) * sol(1:i-1)) / L(i,i)
    end
    x = sol
endfunction


function matrices = factorizacionLU(A)
    m = size(A)(1)
    n = size(A)(2)
    L = eye(m,m)
        for i = 1 : m
            for j = i+1 : m
                coeficiente = A(j,i)/A(i,i)
                L(j,i) = coeficiente
                A(j,i:n) = A(j,i:n) - A(i,i:n) * coeficiente
            end
        end
    U = A
    matrices = [L, U]
endfunction


//Suponiendo no tengo que permutar filas de A, y que A es cuadrada.
//LUx = b ==> Ly = b , Ux = y
function x = ecuacionFactorizacionLU(A, b)
    dimension = size(A)(1)
    LU = factorizacionLU(A)
    y = sustitucionHaciaAbajo(LU(:, 1:dimension), b)
    x = sustitucionHaciaArriba(LU(:, dimension+1:dimension*2), y)
endfunction




--> A_Real = [1 2 -2 1; 4 5 -7 6; 5 25 -15 -3; 6 -12 -6 22;]
 A_Real  =
   1.   2.   -2.    1.
   4.   5.   -7.    6.
   5.   25.  -15.  -3.
   6.  -12.  -6.    22.
--> b_Real = [2; 2; 1; 0]
 b_Real  =
   2.
   2.
   1.
   0.


--> A = [1 2 -2 1; 4 5 -7 6; 6 -12 -6 22; 5 25 -15 -3]
 A  =
   1.   2.   -2.    1.
   4.   5.   -7.    6.
   6.  -12.  -6.    22.
   5.   25.  -15.  -3.
--> b = [2; 2; 0; 1]
 b  =
   2.
   2.
   0.
   1.
--> ecuacionFactorizacionLU(A, b)
 ans  =
   19.5
  -17.
  -18.
  -19.5


--> A_Real * ans
 ans  =
   2.
   2.
   1.
   0.





// Ejercicio 5
//Siendo U cuadrada.
function x = sustitucionHaciaArriba(U, b)
    dimension = size(U)(1)
    sol(dimension) = b(dimension) / U(dimension,dimension)
    for i = dimension-1:-1:1
        sol(i) = (b(i) - U(i, i+1:dimension) * sol(i+1:dimension)) / U(i,i)
    end
    x = sol
endfunction


function x = sustitucionHaciaAbajo(L, b)
    dimension = size(L)(1)
    sol(1) = b(1) / L(1,1)
    for i = 2:dimension
        sol(i) = (b(i) - L(i, 1:i-1) * sol(1:i-1)) / L(i,i)
    end
    x = sol
endfunction


//Supongo que A es cuadrada.
function LU = factorizacionDoolittle(A)
    dimension = size(A)(1)
    L = eye(A)
    U = zeros(A)

    //Despejo la fila 1 de U.
    for j = 1:dimension
        U(1,j) = A(1,j)
    end
    //Despejo la columna 1 de L.
    for k = 2:dimension
       L(k,1) = A(k,1) / U(1,1)
    end

    for i = 2:dimension
        //Despejo valores de la fila i de U.
        for j = i:dimension
           U(i,j) = A(i,j) - L(i, 1:i-1) * U(1:i-1, j)
        end
        //Despejo valores de la columna i de L.
        for k = i+1:dimension
           L(k,i) = (A(k,i) - L(k, 1:i-1) * U(1:i-1, i)) / U(i,i)
        end
    LU = [L, U]
    end
endfunction


//Suponiendo no tengo que permutar filas de A, y que A es cuadrada.
//LUx = b ==> Ly = b , Ux = y
function x = ecuacionFactorizacionLU(A, b)
    dimension = size(A)(1)
    LU = factorizacionDoolittle(A)
    y = sustitucionHaciaAbajo(LU(:, 1:dimension), b)
    x = sustitucionHaciaArriba(LU(:, dimension+1:dimension*2), y)
endfunction


--> A = [1 2 3 4; 1 4 9 16; 1 8 27 64; 1 16 81 256]
 A  =
   1.   2.    3.    4.
   1.   4.    9.    16.
   1.   8.    27.   64.
   1.   16.   81.   256.
--> b = [2; 10; 44; 190]
 b  =
   2.
   10.
   44.
   190.
--> ecuacionFactorizacionLU(A,b)
 ans  =
  -1.
   1.
  -1.
   1.
--> A * ans
 ans  =
   2.
   10.
   44.
   190.




// Ejercicio 6
function matrices = cholesky(A)
    dimension = size(A)(1)
    R = zeros(A)
    Rt = zeros(A)

    Rt(1,1) = sqrt(A(1,1))
    R(1,1) = Rt(1,1)
    //Despejo la primer columna de Rt.
    for i = 2:dimension
        Rt(i,1) = A(i,1) / R(1,1)
        R(1, i) = Rt(i,1)
    end;
    //Despejo el resto de la matriz.
    for i = 2 : dimension
        for j = 2 : i-1
            Rt(i,j) = (A(i,j) - Rt(i,1:j-1) * R(1:j-1, j)) / R(j,j)
            R(j,i) = Rt(i,j)
        end;
        Rt(i,i) = sqrt(A(i,i) - Rt(i,1:i-1) * R(1:i-1,i))
        R(i,i) = Rt(i,i)
    end;
    matrices = [Rt, R]
endfunction

--> A = [16 -12 8 -16; -12 18 -6 9; 8 -6 5 -10; -16 9 -10 46]
 A  =
   16.  -12.   8.   -16.
  -12.   18.  -6.    9.
   8.   -6.    5.   -10.
  -16.   9.   -10.   46.
--> chol(A)
 ans  =
   4.  -3.   2.  -4.
   0.   3.   0.  -1.
   0.   0.   1.  -2.
   0.   0.   0.   5.
--> cholesky(A)
 ans  =
   4.   0.   0.   0.   4.  -3.   2.  -4.
  -3.   3.   0.   0.   0.   3.   0.  -1.
   2.   0.   1.   0.   0.   0.   1.  -2.
  -4.  -1.  -2.   5.   0.   0.   0.   5.