l5.tar

disco.m0000664000175000017500000000047313373642005011624 0ustar  cccsarcccsar#Dibujo de disco de Gershgorin
function [] =disco(A) 
  circulo=0:0.01:2*pi;
  n=length(A);
  x=zeros(length(circulo),n);
  r=sum(abs(A),2)-abs(diag(A)); 
  for i=1:n
    x(:,i)=r(i)*cos(circulo)+real(A(i,i));
    y(:,i)=r(i)*sin(circulo)+imag(A(i,i));
  end
  e=eig(A); 
  plot(x,y,real(e),imag(e),'x');
end


    gauss_seidel.m0000664000175000017500000000073313374045644013200 0ustar  cccsarcccsar  

function [xn,k]=gauss_seidel(A,b,xv,maxit,tol)
  n=length(A);
  xn=NaN*ones(n,1);
  for k=0:1:maxit
  
    for i=1:1:n
      suma1=0;
      
      for j=1:1:i-1
        suma1+=A(i,j)*xn(j); # Unica diferencia
      endfor
      
      suma2=0;
      for j=i+1:1:n
        suma2+=A(i,j)*xv(j);
      endfor
    
      xn(i)=(b(i)-suma1-suma2)/A(i,i);
    
    endfor  
    
    if (norm(xv-xn,inf))<tol
      return
    endif
    
    xv=xn;
    
  endfor
endfunction
    jacobi.m0000664000175000017500000000070013374045645011753 0ustar  cccsarcccsar

function [xn,k]=jacobi(A,b,xv,maxit,tol)
  n=length(A);
  xn=NaN*ones(n,1);
  for k=0:1:maxit
  
    for i=1:1:n
      suma1=0;
      
      for j=1:1:i-1
        suma1+=A(i,j)*xv(j);
      endfor
      
      suma2=0;
      for j=i+1:1:n
        suma2+=A(i,j)*xv(j);
      endfor
    
      xn(i)=(b(i)-suma1-suma2)/A(i,i);
    
    endfor  
    
    if (norm(xv-xn,inf))<tol
      return
    endif
    
    xv=xn;
    
  endfor
endfunction
    lab5.m0000664000175000017500000000240113374053605011342 0ustar  cccsarcccsar#Ejercicio1
%{
#Calculos A1: 

X=-1+(2)*rand(10,10);

A1=X+10*eye(10,10);

b1=A1*ones(10,1);

x0=0*ones(10,1);
maxit=1000;
tol=10^-13;

#Calculos A2: 

Y=-1+(2)*rand(20,20);

A2=Y.'*Y+(10^-3)*eye(20,20);

b2=A2*ones(20,1)
x1=0*ones(20,1);

#Soluciones:

[xn1 k1]=jacobi(A1,b1,x0,maxit,tol)
Mi1=matriziteracion(A1,0);
re1=max(abs(eig(Mi1)))

[xv1 k11]=gauss_seidel(A1,b1,x0,maxit,tol)
Mi2=matriziteracion(A1,1);
re2=max(abs(eig(Mi2)))


[xn2 k2]=jacobi(A2,b2,x1,maxit,tol)
Mi3=matriziteracion(A2,0);
re3=max(abs(eig(Mi3)))

[xv2 k22]=gauss_seidel(A2,b2,x1,maxit,tol)
Mi4=matriziteracion(A2,1);
re4=max(abs(eig(Mi4)))



#Ejercicio2
%}

#Matrices dadas:
M1=[7 3 -1 2 ; 3 8 1 -4 ; -1 1 4 -1 ; 2 -4 -1 6];
j=[-1;0;-3;1];

M2=[4 3 0;3 4 -1; 0 -1 4];
k=ones(3,1);

tol=10^-13;

xaux=0*ones(4,1);
xaux1=0*ones(3,1);

#Graficas:

%{
  Se utiliza el siguiente for para crear dos lista de los resultados 
  de la funcion SOR sobre cada valor dado de W. Cada lista corresponde 
  a una de las matrices dadas
%}
imag=[]
for w=100:1:19900
  imag(w)=sor(M1,j,xaux,w/(10^4),tol);
  imag2(w)=sor(M2,k,xaux1,w/(10^4),tol);
  #plot(w,sor(M2,k,xaux1,w,tol))
endfor

w=100:1:19900

#Grafica 1:

tic();
plot(w,imag(w))
tiempog1=toc()

hold on

#Grafica 2:

tic();
plot(w,imag2(w))
tiempog2=toc()






LU.m0000664000175000017500000000035713373642005011044 0ustar  cccsarcccsarfunction A = LU(A)

  [m,n] = size(A);
  for k=2:1:n
    for i=k:1:n
      alpha = A(i,k-1)/A(k-1,k-1);
      A(i,k-1) = alpha;
      
      for j=k:1:n
        A(i,j) = A(i,j) - (A(i,k-1)*A(k-1,j));
      end
    end
  end
  
endfunction
matriziteracion.m0000664000175000017500000000112113374046171013721 0ustar  cccsarcccsar%{
  Funcion que calcula el radio espectral de la matriz de iteracion para los 
  metodos de Jacobi o Gauss Seidel. Se basa en los teoremas: 
  
    GaussSeidel:  respectral(D^-1(L+U)<1) == La matriz de Seidel converge
    
    Jacobi: respectral(N^-1*P<1) == La matriz de Jacobi converge
    
    En donde: 
    
    D=diagonal principal
    L=diagonal inferior 
    U=diagonal superior 
    N=(D+U)^-1
    P=-L
   
%}

function r=matriziteracion(A,m)
    D=triu(A)-triu(A,1);
    L=triu(A);
    U=tril(A);
  if m==0
    r=inv(D)*(L+U);
  elseif m==1
    r=inv(D+U)*(-L);
  endif
endfunction
octave-workspace0000664000175000017500000000000013373643322013532 0ustar  cccsarcccsarpotenciainversadez.m0000664000175000017500000000050613373642005014415 0ustar  cccsarcccsar
function [x,k,r]=potenciainversadez(A,error,u,itmax) 
  n=length(A); 
  x=ones(n,1);
  dezpla=u*eye(n,n);
  A=A-dezpla;
  A=LU(A); 
  for k=1:1:itmax 
    y=sistemaLU(A,x);
    r=y(1)/x(1);
    if norm(x-(y/norm(y,inf)),inf)<error
      r=1/r+u;
      return 
    endif
    x=y/norm(y,inf);
  endfor
  r=1/r+u;
endfunction
  potenciainversa.m0000664000175000017500000000035513373642005013714 0ustar  cccsarcccsarfunction [x,k,r]=potenciainversa(A,error,x,itmax) 
  A=LU(A); 
  for k=1:1:itmax 
    y=sistemaLU(A,x) ;
    r=y(1)/x(1);
    if norm(x-(y/norm(y,inf)),inf)<error
      return 
    endif
    x=y/norm(y,inf);
  endfor
  r=1/r;
endfunctionpotencia.m0000664000175000017500000000042413373642005012321 0ustar  cccsarcccsar#Busca el maximo autovalor 

function [x,k,r]=potencia(A,error,x,itmax)
  for k=1:1:itmax
    y=A*x; # x^(k+1)=A*x^(k)
    r=y(1)/x(1); # phi(x^(k+1))
    if norm(x-(y/norm(y,inf)),inf)<error
      return 
    endif
    x=y/norm(y,inf); #rescatar a u^(k-1)
 endfor
endfunctionsistemaLU.m0000664000175000017500000000021713373642005012425 0ustar  cccsarcccsarfunction x=sistemaLU(A,b)
  n=length(A);
  Lu=LU(A);
  L=tril(Lu,-1)+eye(n);
  U=triu(Lu);
  y=sust_adel(L,b);
  x=sust_atras(U,y);
endfunctionsor.m0000664000175000017500000000112713374045611011324 0ustar  cccsarcccsar#Implementacion del metodo de sort

%{
  La funcion toma al error relativo como referencia de parada,
  haciendo que la funcion se detenga cuando el error relativo es
  menor a la tolerancia 
 %}
 
function [k]=sor(A,b,x0,w,tol)
  #Separacion de A en A=D+L+U
  n=length(A);
  D=triu(A)-triu(A,1);
  U=triu(A,1);
  L=tril(A,-1);
  cuenta=1;
  xn=(inv(D-w*L))*((w*U+(1-w)*D)*x0+w*b);
  erel=norm(xn-x0,inf)/norm(xn,inf);
  
  while erel>tol 
    xn=(inv(D-w*L))*((w*U+(1-w)*D)*x0+w*b);
    erel=norm(xn-x0,inf)/norm(xn,inf)  ;
    x0=xn;
    cuenta+=1;
  endwhile
  
  x=xn;
  k=cuenta;
endfunction