Métodos Práctica 5

FUNCIONES/PROCEDIMIENTOS ÚTILES:
- spec(A) devuelve los autovalores de A
-   x = poly(0, 'x')
    pol = det(x*eye(A) - A)
    roots(pol)
Hace lo mismo




Ejercicio 1 - Apartado A:
- Ya revisamos que ambas matrices de coeficientes son inversibles.

- Primer sistema: La matriz N no es inversible pues tenemos un 0 en la entrada (1,1) de A. Permutando para que las filas de A estén en el orden (2,3,1), da que el radio espectral de A permutada es menor a 1, lo cual por un Corolario permite afirmar que converge.
- Segundo sistema: El segundo sistema converge pues la matriz de coeficientes es diagonal dominante (Un teorema dice que la sucesión converge si pasa eso).

Ejercicio 1 - Apartado B:
- Primer sistema: Con la permutación de filas (2,3,1) el radio espectral es menor a 1.
- Segundo sistema: La mtriz sin permutar da un radio espectral menor a 1.

Ejercicio 1 - Apartado C:
//Somos gente simple y suponemos que la diagonal de A no tiene elementos nulos.
function x2 = metodoJacobi(A, b, inicio, tolerancia)
    N = size(A)(1)
    x1 = inicio

    for i = 1:N
        sumatoria = 0
        for j = 1:N
            if(j ~= i)
               sumatoria = sumatoria + A(i, j) * x1(j)
            end
        end
        x2(i) = (b(i) - sumatoria) / A(i,i)
    end

    while (norm(x2-x1,2) > tolerancia)
        x1 = x2
        for i = 1:N
            sumatoria = 0
            for j = 1:N
                if(j ~= i)
                   sumatoria = sumatoria + A(i, j) * x1(j)
                end
            end
            x2(i) = (b(i) - sumatoria) / A(i,i)
        end
    end
endfunction


Ap1 = [1 -1 -1; 1 -1 2; 0 2 4]
 Ap1  =
   1.  -1.  -1.
   1.  -1.   2.
   0.   2.   4.
--> b = [0.375; 0; 0]
 b  =
   0.375
   0.
   0.
--> metodoJacobi(Ap1, b, [0; 0; 0], 0.01)
 ans  =
   0.4980469
   0.2460938
  -0.1230469
--> Ap1 * ans
 ans  =
   0.375
   0.0058594
   0.





--> Bp1
 Bp1  =
   1.  -1.   0.
  -1.   2.  -1.
   0.  -1.   1.1
--> b = [0; 1; 0]
 b  =
   0.
   1.
   0.
--> metodoJacobi(Bp1, b, [0; 0; 0], 0.01)
 ans  =
   10.789086
   10.798673
   9.8082602
--> Bp1 * ans
 ans  =
  -0.009587
   1.
  -0.009587





Ejercicio 2:

function [x2, iteraciones] = metodoGaussSeidel(A, b, inicio, tolerancia)
    N = size(A)(1)
    x1 = inicio
    for i = 1:N
        sumatoria = 0
        for j = 1:i-1
            sumatoria = sumatoria + A(i, j) * x2(j)
        end
        for j = i+1:N
            sumatoria = sumatoria + A(i, j) * x1(j)
        end
        x2(i) = (b(i) - sumatoria) / A(i,i)
    end

    iteraciones = 1

    while (norm(x2-x1,2) > tolerancia)
        x1 = x2
        for i = 1:N
            sumatoria = 0
            for j = 1:i-1
                sumatoria = sumatoria + A(i, j) * x2(j)
           end
            for j = i+1:N
                sumatoria = sumatoria + A(i, j) * x1(j)
            end
            x2(i) = (b(i) - sumatoria) / A(i,i)
        end
        iteraciones = iteraciones + 1
    end
endfunction


--> A = [10 1 2 3 4; 1 9 -1 2 -3; 2 -1 7 3 -5; 3 2 3 12 -1; 4 -3 -5 -1 15]
 A  =
   10.   1.   2.   3.    4.
   1.    9.  -1.   2.   -3.
   2.   -1.   7.   3.   -5.
   3.    2.   3.   12.  -1.
   4.   -3.  -5.  -1.    15.
--> b = [12; -27; 14; -17; 12]
 b  =
   12.
  -27.
   14.
  -17.
   12.
--> [res, iter] = metodoJacobi(A, b, [0;0;0;0;0], 10**-6)
 iter  =
   67.
 res  =
   1.0000016
  -2.0000015
   2.9999973
  -1.9999996
   0.9999981
--> [res2, iter2] = metodoGaussSeidel(A, b, [0;0;0;0;0], 10**-6)
 iter2  =
   38.
 res2  =
   1.0000009
  -2.0000007
   2.9999987
  -1.9999999
   0.9999992
--> A * res - b
 ans  =
   0.0000028
  -0.0000023
  -0.0000033
   0.0000009
  -0.0000051
--> A * res2 - b
 ans  =
   0.0000029
  -0.0000018
  -0.000002
  -0.0000004
   0.





Ejercicio 3:

--> N = [2 0 0 0 0; -1 2 0 0 0; 0 -1 2 0 0; 0 0 -1 2 0; 0 0 0 -1 2]
 N  =
   2.   0.   0.   0.   0.
  -1.   2.   0.   0.   0.
   0.  -1.   2.   0.   0.
   0.   0.  -1.   2.   0.
   0.   0.   0.  -1.   2.
--> N \ eye(N)
 ans  =
   0.5       0.       0.      0.     0.
   0.25      0.5      0.      0.     0.
   0.125     0.25     0.5     0.     0.
   0.0625    0.125    0.25    0.5    0.
   0.03125   0.0625   0.125   0.25   0.5
--> A = [2 -1 0 0 0; -1 2 -1 0 0; 0 -1 2 -1 0; 0 0 -1 2 -1; 0 0 0 -1 2]
 A  =
   2.  -1.   0.   0.   0.
  -1.   2.  -1.   0.   0.
   0.  -1.   2.  -1.   0.
   0.   0.  -1.   2.  -1.
   0.   0.   0.  -1.   2.
--> N \ A
 ans  =
   1.  -0.5       0.       0.      0.
   0.   0.75     -0.5      0.      0.
   0.  -0.125     0.75    -0.5     0.
   0.  -0.0625   -0.125    0.75   -0.5
   0.  -0.03125  -0.0625  -0.125   0.75
--> B = N/eye(N)
 B  =
   2.   0.   0.   0.   0.
  -1.   2.   0.   0.   0.
   0.  -1.   2.   0.   0.
   0.   0.  -1.   2.   0.
   0.   0.   0.  -1.   2.
--> B = N \ eye(N)
 B  =
   0.5       0.       0.      0.     0.
   0.25      0.5      0.      0.     0.
   0.125     0.25     0.5     0.     0.
   0.0625    0.125    0.25    0.5    0.
   0.03125   0.0625   0.125   0.25   0.5
--> C=B*A
 C  =
   1.  -0.5       0.       0.      0.
   0.   0.75     -0.5      0.      0.
   0.  -0.125     0.75    -0.5     0.
   0.  -0.0625   -0.125    0.75   -0.5
   0.  -0.03125  -0.0625  -0.125   0.75
--> eye(C) - C
 ans  =
   0.   0.5       0.       0.      0.
   0.   0.25      0.5      0.      0.
   0.   0.125     0.25     0.5     0.
   0.   0.0625    0.125    0.25    0.5
   0.   0.03125   0.0625   0.125   0.25





Ejercicio 4:
Para contar el tiempo, es algo como:
    tic()
    //Inserte función para timear aquí.
    Tiempo = toc()

--> N=500;
--> A = 8*eye(N,N) + 2*diag(ones(N-1,1),1) + 2*diag(ones(N-1,1),-1) + diag(ones(N-3,1),3) + diag(ones(N-3,1),-3);
--> b = ones(N, 1);
--> [sol, timer] = metodoGauss(A,b);
--> timer
 timer  =
   3.5102211

--> N=1000;
--> A2 = 8*eye(N,N) + 2*diag(ones(N-1,1),1) + 2*diag(ones(N-1,1),-1) + diag(ones(N-3,1),3) + diag(ones(N-3,1),-3);
--> b2 = ones(N, 1);
--> [sol2, timer2] = metodoGauss(A2,b2);
--> timer2
 timer2  =
   22.402927

--> N = 500;
--> help zeros
--> ceros = zeros (N, 1);
--> [sol3, timer3] = gaussSeidelTimeado(A, b, ceros, 10**-6);
--> timer3
 timer3  =
   11.317839
--> [sol4, timer4] = gaussSeidelTimeado(A, b, ceros, 10**-12);
--> timer4
 timer4  =
   28.087232

--> N = 1000;
--> ceros2 = zeros(N, 1);
--> [sol5, timer5] = gaussSeidelTimeado(A2, b2, ceros2, 10**-6);
--> timer5
 timer5  =
   45.916292
--> [sol6, timer6] = gaussSeidelTimeado(A2, b2, ceros2, 10**-12);
--> timer6
 timer6  =
   113.07043




Ejercicio 5:

function x2 = metodoSOR(A, b, factorCarga,  inicio, tolerancia)
    N = size(A)(1)
    x1 = inicio

    for i = 1:N
        sumatoria = 0

        for j = 1:i-1
            sumatoria = sumatoria + A(i, j) * x2(j)
        end

        for j = i+1:N
            sumatoria = sumatoria + A(i, j) * x1(j)
        end

        x2(i) = (1-factorCarga) * x1(i) + factorCarga * (b(i) - sumatoria) / A(i,i)
    end

    while (norm(x2-x1,2) > tolerancia)
        x1 = x2
        for i = 1:N
            sumatoria = 0

            for j = 1:i-1
                sumatoria = sumatoria + A(i, j) * x2(j)
            end

            for j = i+1:N
                sumatoria = sumatoria + A(i, j) * x1(j)
            end

            x2(i) = (1-factorCarga) * x1(i) + factorCarga * (b(i) - sumatoria) / A(i,i)
        end
    end
endfunction

// Suponemos ya se sabe que A es definida positiva y tridiagonal.
function solucion = metodoSORDisney(A, b, inicio, tolerancia)
    tamano = size(A)(1)
    N = zeros(A)
    for i=1:tamano
        N(i,i) = A(i,i)
    end

    Pj = eye(A) - N \ A
    radioE = max(abs(spec(Pj)))
    w = 2/(1+sqrt(1-radioE**2))
    solucion = metodoSOR(A, b, w, inicio, tolerancia)
endfunction



--> [sol1, iter1] = metodoGaussSeidelContador(A, b, [0; 0; 0], 10**-7)
 iter1  =
   36.
 sol1  =
   3.0000001
   3.9999999
  -5.
--> [sol2, iter2] = metodoSORDisney(A, b, [0; 0; 0], 10**-7)
 iter2  =
   46.
 sol2  =
   3.
   4.
  -5.