El post que buscas se encuentra eliminado, pero este también te puede interesar

Varios Ejemplos Matlab Codigos

Anuncios

%-----------------------------------------------------------------------------
% MIS CODIGOS PARA PRACTICAR IF FOR ELSE WHILE
%-----------------------------------------------------------------------------
(Al final puede verse el formato que se dio a esa práctica)


%------------------------------------------------
function [L,U,p]=lupa(A,piv)
%------------------------------------------------
% function [L,U,p] = lupa(A,piv)
% Halla la factorización LU de A (que puede ser rectangular),
% salvo que encuentre un pivote nulo.
% Lo hace SIN pivotar si el argumento 'piv' es 0 (o no se da)
% y PIVOTANDO si se da 'piv' distinto de 0 . En ese caso,
% halla L,U tales que L*U = P*A , pero devuelve simplemente
% un vector 'p' que da el orden de las filas en P*A.
% Si size(A)=[m,n] , entonces size(L)=[m,m] , U es como A .

if nargin<2, piv=0; end % si se ha llamado 'lupa(A)'
[m,n] = size(A);
q = 1:m; % filas del bloque restante
p = q; % orden inicial de las filas

%-------------------------------- loop de ELIMINACION:
for k = 1:min(m-1,n)
if piv % escoge la fila ip del PIVOTE:
[a,ip] = max(abs(A(q,k))); ip = q(ip);
p([k,ip]) = p([ip,k]);
A([k,ip], = A([ip,k],;
end
if A(k,k)==0, error('pivote nulo'), end % !!

q = k+1:m; col = k+1:n; % el bloque restante es A(q,col)
A(q,k) = A(q,k)/A(k,k); % multiplicadores
A(q,col) = A(q,col) - A(q,k)*A(k,col);
end
%--------------------------------
U = triu(A); L = A-U;
L(m,m) = 0; % hay que agrandar L, si n < m
L = L(:,1:m)+eye(m);



%------------------------------------------------
function X = luaxb(A,B,piv)
%------------------------------------------------
% function X = luaxb(A,B,piv)
% Usa la factorización LU de A (llamando a 'lupa(A,piv)'),
% para resolver el S.E.L. A*X = B .
% B puede tener varias columnas, pero A debe ser cuadrada.
% Lo hace SIN pivotar si el argumento 'piv' es 0 (o no se da)
% y PIVOTANDO si se da 'piv' distinto de 0 .

[m,n]=size(A);
if m ~= n, error('A no es cuadrada'),
elseif nargin<3, piv=0;
end
% ahora llama a la factorización LU
[L,U,p]=lupa(A,piv);
if prod(diag(U))==0, error('U no es regular'), end

X=B(p,; % dará error si B tiene <m filas;
% si tiene >m , las ignora

for k=1:m-1 % ETAPA HACIA DELANTE ...
q=k+1:m;
q, = q, - L(q,k)*k,;
end

for k=m1 % ... Y ETAPA HACIA ATRAS
k, = k,/U(k,k);
q=1:k-1; q, = q, -U(q,k)*k,;
end


%------------------------------------------------
function [Q,R] = qra(A,met)
%------------------------------------------------
% function [Q,R]=qra(A,met)
%
% Produce los factores QR de A usando Gram-Schmidt clásico si met=0,
% y en caso contrario, Gram-Schmidt modificado.
% Admite A de cualquier forma y rango, y devuelve factores "reducidos":
% con r columnas en Q y filas en R, r=rango(A).

if nargin<2, met=1; end
[m,n] = size(A);
Q=zeros(m,0); % para empezar en k=1 sin que " Q'*v " dé error
R=0*A;
K=[]; % K = índices de filas R y cols Q que se conservarán

if met==0
%------------------------------------ Gram-Schmidt clásico
for k=1:n
v = A(:,k); j=1:k-1;
R(j,k) = Q'*v; v = v -Q*R(j,k); d = norm(v);
if d > 0
Q(:,k) = v/d; R(k,k) = d; K=[K,k];
end
end

else
%------------------------------------ Gram-Schmidt modificado
for k=1:n
q = A(:,k); d = norm(q);
if d > 0
q = q/d; Q(:,k) = q;
j= (k+1):n; p = q'*A(:,j); R(k,k:n) = [d,p];
A(:,j) = A(:,j) -q*p; K=[K,k];
end
end
%------------------------------------
end
% retira lineas inútiles:
Q = Q(:,K); R = R(K,;


%------------------------------------------------
function [W,R,Q]=house(A)
%------------------------------------------------
%
% function [W,R,Q]=house(A)
% Halla la factorización QR de A mediante reflexiones de Householder,
% y devuelve el factor R y los vectores unitarios W_j que realizan
% las reflexiones H_j , como columnas de la matriz W .
% Si se la llama con el tercer argumento de salida, devuelve en él
% el factor ortogonal Q , producido a partir de la W .

[m,n] = size(A);
%------------------------------------ Householder
K=[]; % K = índices de las W_k no nulas
for k= 1:min(m-1,n) % es el número de columnas a "limpiar"
i=k:m; % sólo hay que operar con estas filas
v = A(i,k); d = norm(v);
if d > 0 % si fuese v=0, nada que hacer esta vez
K=[K,k];
sign = 1-2*(v(1)>0); % signo = -1 si v(1)>0
A(k,k) = sign*d;
v(1) = v(1) - A(k,k); v = v/norm(v);
W(i,k) = v;
j= (k+1):n; % columnas en las que opera H_k :
A(i,j) = A(i,j) - 2*v*(v'*A(i,j));
end
end

%------------------------------------
W=W(:,K); % limpia sobrantes en W y R
R=triu(A);
if nargout>2 % si se ha pedido, calcula Q ...
Q = eye(m); % ... aplicando a I las H_k ...
for j = size(W,2)1 % ... en orden inverso
v = W(:,j);
Q = Q - 2*v*(v'*Q);
end
end


%------------------------------------------------
function X = qraxb(A,B,met)
%------------------------------------------------
% function X = qraxb(A,B,met)
%
% Usa la factorización QR de A para resolver el S.E.L. A*X = B ,
% llamando
% a 'qra(A,met)' si se da met = 0 ó 1
% a 'house(A)' si no se da ese argumento de entrada.
%
% Si A es rectangular, pero con rango = no.de cols , calcula por
% ese método la solución LS .

if nargin<3
%------------------------------------
[W,R]=house(A);
for j = 1:size(W,2) %
v = W(:,j); B = B - 2*v*(v'*B);
end
X = RB;

elseif met==0
%------------------------------------ Gram-Schmidt clásico
[Q,R]=qra(A,met);
B = Q'*B; X = RB;

else
%------------------------------------ Gram-Schmidt modificado
[Q,R]=qra(A);
B = Q'*B; X = RB;

%------------------------------------
end


%------------------------------------------------
function [X,res,normrk] = iteraxb(A,B,w,tol)
%------------------------------------------------
% function [X,res] = iteraxb(A,B,w)
%
% Usa un método iterativo para el S.E.L. A*X = B :
% Jacobi, si se da w = 0,
% SOR, si se da un valor 0 < w < 2 ; en particular,
% Gauss-Seidel si se da w = 1 .
% Devuelve la aproximación X y su correspondiente res = B - A*X .
%
% Condiciones de parada:
% que el 'res' crezca a un tamaño muy superior al de B ;
% que descienda hasta una fracción 'tol' de ese tamaño, que
% puede darse como otro argumento de entrada (por defecto, 1/10^6).
% Si se la llama con un 3er argumento de salida, devuelve el
% historial de la convergencia: las sucesivas norm(r_k)/norm(B).

if nargin < 4, tol=10^-6; end
m = size(B,1);
if any( size(A) ~= m ), error('datos inadecuados'), end
tol = tol* norm(B); % tamaño de parada para el residual
it = 1; M = 100; % contador y máximo de iteraciones
if w==0
S = diag(diag(A)); %----------------------- Jacobi
else
S = tril(A); %------------------------ SOR
end

X = 0*B; % valor inicial x = 0
res = B; % residual inicial
normrk = norm(res);
while normrk(it) > tol & it < M
if w==0, X = X + (Sres);
else, X = X + w* (Sres);
end
res = B - A*X;
it = it+1; normrk(it) = norm(res);
end
if nargout > 2, normrk = normrk/norm(B); end




%------------------------------------------------ SESION 3 DE LAB (J 25- L 29/10)

%--------------------------------------------------------------------
% EXPERIMENTOS CON 'lu'

% Producir una A a partir de sus factores LU,
% de tal modo que NO HAYA QUE PIVOTAR ??

m=5; % m = tamaño escogido
format bank
d=5; % Cortaremos'LU' para producir L,U ,
LU = ceil(d*rand(m))/d % ... y no queremos un |l_ik|>1 ,
% ... ni ceros en la diagonal de U.

A = (tril(LU,-1)+eye(m))*(triu(LU))

[L,U] = lu(A) % y Matlab recupera los factores L,U
format short


%--------------------------------------------------------------------
% RESOLVIENDO UN S.E.L. Ax=b con 'Ab'...

% Producimos b=A_1+A_2 :
z=zeros(m,1); z(1:2)=1;
b = A*z;
x = Ab % y Matlab recupera x ...
err = x-z % ... con errores ...
res = b-A*x % consistentes con los k(L), k(U), k(A):
cnLUA = [cond(L),cond(U),cond(A)]

% (!!) Si se toma d=4 al producir A (bloque anterior) los errores ... ??


%--------------------------------------------------------------------
% ... Y CON LOS METODOS ITERATIVOS CLASICOS

m=5;
A=rand(m)/m+eye(m) % Ahora necesitamos A para las que
cn = cond(A) % ... el método converja
z=zeros(m,1); z(1:2)=1;
b = A*z;

% Jacobi:
S = A.*eye(m); T = S-A; normaDG = norm(ST)
x = 0*b;
for it=1:10
xj = S(T*x+b); err(it) = norm(xj-x), res(it) = norm(b-A*xj)
x=xj, pause
end

% Gauss-Seidel:
S = triu(A); T = S-A; normaDG = norm(ST)
x = 0*b;
for it=1:10
xj = S(T*x+b);
errG(it) = norm(xj-x); resG(it) = norm(b-A*xj); x=xj;
end
error = norm(x-z), residual = norm(b-A*x)

it=1:10;
semilogy(it,err,'*k',it,res,'.k',it,errG,'*r',it,resG,'.r')


% Comparemos el rendimiento de ambos,
m=400; % ... con m "grande"
A=rand(m)/m+eye(m);
cn = cond(A)
z=zeros(m,1); z(1:2)=1;
b = A*z;

% El método DIRECTO: (Matlab usa LU)
tic, x=Ab; toc

% ... frente a Gauss-Seidel:
S = triu(A); T = S-A;
normaDG = norm(ST)
tic
for it=1:10, x = S(T*x+b); end
toc
error = norm(x-z), residual = norm(b-A*x)


%--------------------------------------------------------------------
% PERDIDA INESPERADA DE RANGO: las L de LU son MUY especiales!!

for m=40:20:120
for it = 1:20
i = it+m-40;
A = rand(m); % (!!) probar con A = randn(m);

[L,U,P] = lu(A); kn(i, = [cond(L),cond(U),cond(A)];
ra = [rank(L),rank(U),rank(A)]; if min(ra)<m, ra, end
LA = tril(A); K(i) = cond(LA); r(i) = rank(LA);
end
end
i = 1:i; m = 40+fix((i-1)/20)*20;

% figure
subplot(2,1,1), semilogy(i,kn,i,m,'k')
ylabel('rango y cond de las L, U, A')

subplot(2,1,2), semilogy(r,K,'r.'), grid
xlabel('rango de las "L" recortadas: ligeramente deficitario')
ylabel('sus "cond"')



%------------------------------------------------ SESION 5 DE LAB (L 12, J 15/11).

%--------------------------------------------------------------------
% A) fase de prueba del código con mis propias versiones de las 6 funciones:
% pruebas con matrices A pequeñas (que recomiendo repetir con las de cada uno),
% para ver simplemente si funcionan como es debido.
%--------------------------------------------------------------------

m=4; A=rand(m); clc

%----------------------- Factorización LU :
[L0,U0,P] = lu(A) % con la de Matlab
[L,U,p] = lupa(A,1) % con la nuestra
bkerr = A(p, - L*U
help lupa, pause, disp('... y sin pivotar:')
[L,U] = lupa(A);
bkerr = A - L*U , pause, clc

%----------------------- Factorización QR :
[Q0,R0] = qr(A) % la de Matlab
[Q,R] = qra(A) % la que probamos: Gr-Sch modificado, ...
disp(' OJO A LOS SIGNOS !!'),
bkerr = A - Q*R
pause, disp('... ahora, Gr-Sch clásico:')
[Q,R] = qra(A,0)
bkerr = A - Q*R , help qra
pause, disp('... y ahora, Householder:')
[W,R,Q] = house(A)
bkerr = A - Q*R , help house, pause, clc

%----------------------- Resolver Ax=b :
x0=zeros(m,1); x0(1:2)=1, b=A*x0;
xM = Ab; % con Matlab
errMatlab = xM - x0
x = luaxb(A,b,1); % con LU
errpiv = x - x0, help luaxb, pause
x = luaxb(A,b,0);
errNOpiv = x - x0, pause
x = qraxb(A,b,1); % llamando a 'qra'
errGSch = x - x0, help qraxb, pause
x = qraxb(A,b); % llamando a 'house'
errHouse = x - x0, pause, clc

%----------------------- M»todos iterativos
A=A+sqrt(m)*eye(m); b=A*x0;
format short e, clc

[X,res,hist] = iteraxb(A,b,0), help iteraxb, pause
tol = 10^-10
w=1, [X,res,hist] = iteraxb(A,b,w,tol); hist

S = tril(A); DG = S(S-A);
ep=10^-6; plot(eig(DG)+ep*i,'*r'), axis equal, hold on
format short, rho = max(abs(eig(DG)))
t = 0:0.01:1; z = exp(j*2*pi*t); plot([z,rho*z]), hold off


%--------------------------------------------------------------------
% B) fase de exploración (comienzo):
%--------------------------------------------------------------------

%----------------------- Métodos iterativos
m=40; d=5;
U=triu(ones(m),1); U=U^d; U=U+U'; P=1-(U>0);
A=randn(m)+sqrt(m)*eye(m); A=A.*P; spy(A), pause
x0=zeros(m,1); x0(1:2)=1; b=A*x0;
tol = 10^-10
w=1, [X,res,hist] = iteraxb(A,b,w,tol); hist, % rhow = max(abs(1+w*(eig(DG)-1)))

S = tril(A); DG = S(S-A);
format short, rho = max(abs(eig(DG))),
ep=10^-6; plot(eig(DG)+ep*i,'*r'), axis equal, hold on
t = 0:0.01:1; z = exp(j*2*pi*t); plot([z,rho*z]), hold off


%----------------------- Factorización QR :
m=100; A=rand(m); [Q0,R0]=qr(A);
% R0=R0*diag(sign(diag(R0))); A=Q0*R0;
% R0=triu(A)*diag(sign(diag(A))); A=Q0*R0;
kA = cond(A)
disp('ERRORES:')
disp('norm(Q*R-A)/norm(A) , cond(Q)-1 , norm(Q-Q0) , norm(R-R0)/norm(R0)')
format short e
disp(' -------------------- Gr-Sch CLASICO:')
[Q,R] = qra(A,0); errores0 = [norm(Q*R-A)/norm(A), cond(Q)-1, norm(Q-Q0), norm(R-R0)/norm(R0)];
disp(errores0), pause
disp(' -------------------- modificado:')
[Q,R] = qra(A); errores1 = [norm(Q*R-A)/norm(A), cond(Q)-1, norm(Q-Q0), norm(R-R0)/norm(R0)];
disp(errores1), pause
disp(' -------------------- Householder:')
[W,R,Q] = house(A); erroresH = [norm(Q*R-A)/norm(A), cond(Q)-1, norm(Q-Q0), norm(R-R0)/norm(R0)];
disp(erroresH), pause
disp(' -------------------- qr de Matlab:')
[Q,R] = qr(A); erroresM = [norm(Q*R-A)/norm(A), cond(Q)-1, norm(Q-Q0), norm(R-R0)/norm(R0)];
disp(erroresM), format short


%----------------------- Resolver Ax=b :
x0=zeros(m,1); x0(1:2)=1; B=A*x0;
disp(' -------------------- Householder:')
X = qraxb(A,B); errH = norm(X-x0), pause
disp(' -------------------- Gr-Sch CLASICO:')
X = qraxb(A,B,0); err0 = norm(X-x0), pause
disp(' -------------------- Modificado:')
X = qraxb(A,B,1); err1 = norm(X-x0), pause
disp(' -------------------- LU:')
X = luaxb(A,B); err0 = norm(X-x0), pause
disp(' -------------------- LU con pivotaje:')
X = luaxb(A,B,1); err1 = norm(X-x0), pause
disp(' -------------------- Ab de Matlab:')
errMatlab = norm((AB)-x0)


%--------------------------------------------------------------------
% PLAZO Y FORMA PARA LA PRACTICA 2:
%--------------------------------------------------------------------

Deberá estar entregada al mediodía (12:00) del V 23/11, dejando cada
equipo EN LA CARPETA 'practica2**', QUE SE LES VA A CREAR PARA ELLO,
los documentos siguientes:

* 6 M-files con las funciones cuyos nombres y uso se especifica abajo;

* un documento de texto que explique (con toda la concisión que lo
complejo del asunto permite) las pruebas que han hecho con el código,
comentando los resultados que consideren más relevantes o hayan
llamado su atención, en relación con los siguientes aspectos:

- costes: tiempos de ejecución, comparados entre sí y/o con los de
rutinas equivalentes de Matlab ('lu(A), qr(A), Ab');

- precisión y estabilidad: errores "back-/forward" y su relación con
los 'cond' de las diversas matrices;

- y la evolución de ambas cosas al crecer el tamaño de las matrices.

En la sesión de esta semana próxima ofreceré ejemplos de ese tipo de
indagación (usando mis M-files); naturalmente, se pueden escribir
varios volúmenes sobre algo así (más aún con todos los métodos que
abarcamos en la Práctica), de modo que una selección prudente será
parte importante del mérito del trabajo.

Las 6 funciones deben tener los nombres, forma y comportamiento que
se explica abajo. Las especificaciones precedidas por
  • son
  • "OPTATIVAS": se debe conseguir como mínimo que los códigos funcionen
    para una matriz cuadrada y regular, y añadir esas otras opciones sólo
    si sobra tiempo y ganas.

    Para verificar que funcionan como es debido, se recomienda usar
    datos de tamaño pequeño (n=4, por ejemplo) producidos con 'rand' y/o
    'randn' para poder visualizar los resultados, además de compararlos
    con los que calcula Matlab.
    Esa 'fase A' de las pruebas NO DEBE formar parte del "informe" que
    acompañe a los códigos, en el que debe comentarse sólo lo que pasa
    al probar con datos más "exigentes".

    ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
    NO OLVIDAR incluir el NUMERO DE EQUIPO Y NOMBRES DE SUS COMPONENTES
    en las lineas iniciales de comentario de CADA FUNCION, que deben
    seguir inmediatamente a la linea de declaración de la misma, para
    que Matlab las reconozaca y use como su 'help'.
    ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
    Y naturalmente, redactar esos comentarios pensando en que sirvan para
    eso: copiar los que pongo a continuación y añadir, además de los DATOS
    citados, lo que proceda explicar si se han programado más opciones de
    entrada o salida.

    %--------------------------------------------------------------------
    function [L,U,p]=lupa(A,piv)
    %
    % Halla la factorización LU de A (cuadrada),
    % pero para y devuelve 'error' si encuentra un pivote nulo.
    % Lo hace SIN pivotar si el argumento 'piv' es 0 (o no se da)
    % y PIVOTANDO si se da 'piv' distinto de 0 . En ese caso,
    % halla L,U tales que L*U = P*A , pero devuelve simplemente
    % un vector 'p' que da el orden de las filas en P*A.

    %--------------------------------------------------------------------
    function X = luaxb(A,B,piv)
    %
    % Usa la factorización LU de A (llamando a 'lupa(A,piv)'),
    % para resolver el S.E.L. A*X = B .
    % Por lo tanto, lo hace SIN pivotar si el argumento 'piv' es 0
    % (o no se da) y PIVOTANDO si se da 'piv' distinto de 0 ,
    % y devuelve 'error' si 'lupa' encuentra algún pivote nulo.
    %
    % B puede tener varias columnas, pero A debe ser cuadrada.

    %--------------------------------------------------------------------
    function [Q,R]=qra(A,met)
    %
    % Produce los factores QR de A usando Gram-Schmidt clásico si met=0,
    % y en caso contrario, Gram-Schmidt modificado.

    %
  • Admite A de cualquier forma y rango, y devuelve factores "reducidos":
  • % con r columnas en Q y filas en R, r=rango(A).

    %--------------------------------------------------------------------
    function [W,R]=house(A)
    %
    % Halla la factorización QR de A mediante reflexiones de Householder,
    % y devuelve el factor R y los vectores unitarios W_j que realizan
    % las reflexiones, como columnas de la matriz W .

    %
  • Si se la llama con un tercer argumento de salida, devuelve en él
  • % el factor ortogonal Q , producido a partir de la W .

    %--------------------------------------------------------------------
    function X = qraxb(A,B,met)
    %
    % Usa la factorización QR de A para resolver el S.E.L. A*X = B ,
    % llamando
    % a 'qra(A,met)' si se da met = 0 ó 1
    % a 'house(A)' si no se da ese argumento de entrada.

    %
  • Si A es rectangular, pero con rango = no.de cols , calcula por
  • % ese método la solución LS .

    %--------------------------------------------------------------------
    function [X,res] = iteraxb(A,B,w)
    %
    % Usa un método iterativo para el S.E.L. A*X = B :
    % Jacobi, si se da w = 0,
    % SOR, si se da un valor 0 < w < 2 ; en particular,
    % Gauss-Seidel si se da w = 1 .
    % Devuelve la aproximación X y su correspondiente res = B - A*X .
    %
    % Condiciones de parada:
    % que el 'res' crezca a un tamaño muy superior al de B ;
    % que descienda hasta una fracción fijada de ese tamaño.

    %
  • Esa fracción puede darse como un 4o. argumento de entrada 'tol'.
  • % Si se la llama con un 3er argumento de salida, devuelve el historial
    % de la convergencia: la norma de los sucesivos residuales /norm(B) .

    Anuncios

    4 comentarios - Varios Ejemplos Matlab Codigos

    @gisss
    te dejé 10 y no te doy un abrazo porque no sé dónde vivís... muchas gracias!!!
    @len_520317
    xvre broooo.... me sirvio thanks!!!