% rutero para pregunta 2, laboratorio 5, 521230, S2-2011

% construir matriz con m=10, n=5
m = 10; n= 5;
[A, t] = matriz_ej2_lab5(m,n);
% calculando descomposicion QR de A
[Q, R] = qr(A);
% Q debe ser ortogonal, Q' es inversa de Q
Q*Q'- eye(m+1)
% R contiene a una matriz triangular superior de
% en sus primeras n+1 filas
spy(sparse(R))
% rango de R debe ser n+1
rank(R)
% R1 es la matriz triangular superior de tamanno n+1 en primeras n+1 filas
% de R
R1 = R(1:n+1,:);
Q1 = Q(:,1:n+1);
Q2 = Q(:,n+2:end);
% A debe ser igual a Q1*R1
A - Q1*R1
% matriz de sistema de ecuaciones normales
B = A'*A;
% B es peor condicionada que R1
cond(B)
cond(R1)^2
H = chol(B)
H
R1

% ajustando exp(t)
% parte derecha de sistema sobredeterminado
y = exp(t);
% coeficientes de polinomio se obtienen de resolver R1x = Q1'*y
coefpol1 = R1\(Q1'*y);
% mostrando puntos
plot(t,y,'*')
% mostrando polinomio
hold on
plot(t,polyval(coefpol1,t));
% comprobando que norma de este vector es norma de y-Ac
norm(Q2'*y)
norm(y-A*coefpol1)

% ajustando sin(pi t) con errores aleatorios de tamanno 0.05
y1 = sin(pi*t);
y2 = sin(pi*t) + (2*rand(m+1,1)-1)*0.05; % valores con errores
% coeficientes de polinomio se obtienen de resolver R1x = Q1'*y
coefpol2 = R1\(Q1'*y2);
figure(2)
% mostrando puntos errores, polinomio (que ajusta valores con errores)
% y funcion
plot(t,y2,'*',t,polyval(coefpol2,t),t,y1,'r');
% comprobando que norma de este vector es norma de y-Ac
norm(Q2'*y2)
norm(y2-A*coefpol2)






