Libro Álgebra lineal numérica: archivos descargables

function [B,Q]=givens(A,m,n,columna,fila1,fila2)
Q=eye(m);
a=A(fila1,columna);
b=A(fila2,columna);
if b==0
    c=1;
    s=0;
else
    if abs(b)>abs(a)
        t=-a/b;
        s=1/sqrt(1+t^2);
        c=s*t;
    else
        t=-b/a;
        c=1/sqrt(1+t^2);
        s=c*t;
    end
end

Q(fila1,fila1)=c;
Q(fila1,fila2)=-s;
Q(fila2,fila2)=c;
Q(fila2,fila1)=s;

B=Q*A;
end
 

function [Q,R]=myQR(A)
    [m,n] = size(A);
        Q = eye(m,m);
        for j=1:n
            for i=m:-1:j+1
                [A,Q1]=givens(A,m,n,j,i-1,i);
                Q=Q1*Q;
            end
        end
        R=A;
        Q=Q';
end
 

clear all;
C=imread('Imagen.jpg');
D=double(C(:,:,1));
A=color2negro(D);
size(D)
hold on
subplot(3,5,1);image(A);title('Original');colormap('gray');
r=rank(A)
[U,S,V] = svd(A); B=255*ones(size(A));
colormap('gray');
p=0;
for i=1:r
    B=U(:,1:i)*S(1:i,1:i)*V(:,1:i)';
    if i==1 || mod(i,5)==0 ||i==r
        B=color2negro(B);
        figure; 
        %B=color2negro(B);
        colormap('gray');
        image(U*S*V')
        image(B)     
        title(['r=',num2str(i)]);
        p=p+1;
    end
end

hold off  

La siguiente rutina cambia una imágen a color a otra en blanco y negro.

Archivo: color2negro.m

function B = color2negro(A)
[m,n]=size(A);
for i=1:m
    for j=1:n
        if abs(A(i,j)-255)<128
            B(i,j)=255;
        else
            B(i,j)=0;
        end
    end
end
end  

function b = epsilon()
  a=1;
  while a+1>1
    b=a;
    a=a/2;
   end
end

function [P,L,U] = pivparcial(A)
m=length(A);
U=A;
P=eye(m);
L=eye(m);
for k=1:m-1
    [mx,fila]=max(abs(U(k:m,k)));
    p=(k-1)+fila;
    U=permutar(U,k,p,k,m);
    P=permutar(P,k,p,1,m);
    L=permutar(L,k,p,1,k-1);
    for j=k+1:m
        L(j,k)=U(j,k)/U(k,k);
        U(j,k:m)=U(j,k:m)-L(j,k)*U(k,k:m);
    end
    
end

end
  

Archivo: permutar.m

function [A] = permutar(A,i,j,col1,col2)
    Vec=A(i,col1:col2);
    A(i,col1:col2)=A(j,col1:col2);
    A(j,col1:col2)=Vec;
end
  

function [L,D] = fac_simetrica(A)
clc
m=length(A)
D=zeros(m)
L=eye(m)
r=zeros(m,1)
for k=1:m
    for p=1:k-1
        r(p)=L(k,p)*D(p,p)
    end
    r(k)=A(k,k)-L(k,1:k-1)*r(1:k-1)
    D(k,k)=r(k)
    A(k+1:m,k)
    L(k+1:m,1:k-1)
    r(1:k-1)
    L(k+1:m,k)=(A(k+1:m,k)-L(k+1:m,1:k-1)*r(1:k-1))./r(k);
end
    

end

A=[34 4;8 20];
x=[1;1];
clc
fprintf('p \t x1 \t x2 \n');
while A*x-for i=1:5
p=(x'*A*x)/(x'*x);
B=A-p*eye(2);
x1=B\x;
x=(1/norm(x1))*x1;
fprintf('%2.0f \t %f \t %f \t %f \n',i,p,x(1),x(2));
end