% Loesung Aufgabe 23: 

close all 
clear all

% Vektor, der die Messpunkte enthaelt: 
t = zeros(6,1);
t(1)=0.21;
t(2)=0.62;
t(3)=1.19;
t(4)=2.01;
t(5)=2.42;
t(6)=4.18;

% Vektor, der die Messwerte zu den Messpunkten enthaelt: 
b = zeros(6,1);
b(1)=2.23;
b(2)=2.49; 
b(3)=4.22; 
b(4)=6.13; 
b(5)=6.95;
b(6)=1.01;

% Systemmatrix
A = [ones(6,1) t t.^2 t.^3];
disp('                         ')
disp('Die Systemmatrix ist A = ')
disp(A)
disp('Die rechte Seite ist b = ')
disp(b)
disp('                         ')

% es gibt mehrere Moeglichkeiten zur Loesung von Ax=b.

% 1: Der Matlab Operator \ berechnet die Loesung von Ax=b via x=A\v. Falls
% A ueberbestimmt ist, liefert A\b automatisch die kleinste Quadrate
% Loesung. 

% 2: Bestimme QR-Zerlegung von A via [Q,R]=qr(A) und loese Rx = Q*b 

% 3: Ermittle mit pinv(A) die Pseudoinverse von A und berechne x =
% pinv(A)*b. 

% 4: Loese System ueber die Normalengleichungen, dh loese A.'A x = A.'b

% 1:
x1=A\b;
alpha1=x1(1);
beta1=x1(2);
gamma1=x1(3);
delta1=x1(4);

% 2: 
[Q,R]=qr(A); 
b1 = Q'*b;
x2 = R\b1;
alpha2=x2(1);
beta2=x2(2); 
gamma2=x2(3);
delta2=x2(4);

% 3: 
x3 = pinv(A)*b;
alpha3=x3(1);
beta3=x3(2);
gamma3=x3(3);
delta3=x3(4);

% 4: 
x4 = A.'*A\A.'*b;
alpha4=x4(1);
beta4=x4(2);
gamma4=x4(3);
delta4=x4(4);




disp('                                    ')
disp('Das System reduziert sich zu dem LGS')
disp('                                    ')
disp('           Rx=(Q^*)*b               ')
disp('                                    ')
disp('mit R =                             ')
disp(R)
disp('                                    ')
disp('und (Q^*)*b =                       ')
disp(b1)
disp('                                    ')


Fehler = norm(A*x2-b,2);


disp(['Der Fehler betraegt ||Ax-b||_2  = ',num2str(Fehler)]) 
disp(' ')
disp(['Die Ausgleichsfunktion ist y(t) = (',num2str(alpha2),') + (',num2str(beta2),')*t + (',num2str(gamma2),')*t^2 + (',num2str(delta2),')*t^3' ])
N = 40;
L1=-1;
L2=5;
h = (L2-L1)/(N-1);
t_plot=[L1:h:L2]';

y = alpha2*ones(N,1)+beta2*t_plot+gamma2*t_plot.^2+delta2*t_plot.^3;


figure(1)
plot(t_plot,y,'k')
hold on
plot(t(1),b(1),'x')
plot(t(2),b(2),'x')
plot(t(3),b(3),'x')
plot(t(4),b(4),'x')
plot(t(5),b(5),'x')
plot(t(6),b(6),'x')
hold off
title(['Ausgleichsfunktion y(t) = (',num2str(alpha2),') + (',num2str(beta2),')*t + (',num2str(gamma2),')*t^2 + (',num2str(delta2),')*t^3' ])
xlabel('t')
ylabel('y')
legend('Ausgleichsfunktion','Messpunkte',4)