% SETTIMANA 11
% ESERCIZIO 1.1

X=sym([3;-1;2]) 
X1=sym([1;2;-1])
X2=sym([2;0;1])

% in questo caso siamo in R^3 e quindi la proiezione ortogonale di X su U si potrebbe trovare con il prodotto vettoriale. Utilizziamo invece l'algoritmo di Gram-Schmidt per costruire una base ortonormale di U a partire da X1 ed X2:

F1=X1
a12=sym((X2'*F1)/(F1'*F1))
F2=X2-a12*F1

% {F1, F2} è una base ortogonale di U. La proiezione di X su U è data da

a=sym(X'*F1)/sym(F1'*F1)
b=sym(X'*F2)/sym(F2'*F2)
ProjX=a*F1+b*F2
Y=X-ProjX

% Quindi Y appartiene al complemento ortogonale di U e

X=ProjX+Y

%
% ESERCIZIO 1.5

syms a b c d
X=[a;b;c;d]
X1=sym([1;-1;2;0])
X2=sym([-1;1;1;1])

% Utilizziamo l'algoritmo di Gram-Schmidt per costruire una base ortonormale di U a partire da X1, X2 ed X3:

F1=X1
a12=sym((X2'*F1)/(F1'*F1))
F2=X2-a12*F1

% {F1, F2}={X1,X2} è una base ortonormale di U. La proiezione di X su U è data da

s=sym(F1'*X)/sym(F1'*F1)
t=sym(F2'*X)/sym(F2'*F2)
ProjX=s*F1+t*F2
Y=X-ProjX

% Quindi Y appartiene al complemento ortogonale di U e

X=ProjX+Y

%
% ESERCIZIO 2.1

syms x y
eqn1= x-y==1
eqn2= 2*x+y==-1
eqn3= x+5*y==-4
[A,b]=equationsToMatrix(eqn1,eqn2,eqn3)

% Vediamo se il sistema è compatibile:

linsolve(A,b)

% Non lo è! Cerchiamo le soluzioni approssimate. Abbiamo visto due metodi: 1) Risolvere il sistema A'AX=A'b. 2) Trovare la decomposizione QR di A=QR, e risolvere RX=Q'b.
% Intanto vediamo se c'è un'unica soluzione approssimata, studiando il rango di A

rank(A)

% Poichè A ha due colonne, ed il suo rango è due, ne segue che A'A è invertibile e quindi c'è un'unica soluzione approssimata.
% Metodo 1
% Poniamo

C=A'*A
d=A'*b

% e risolviamo:

Z=linsolve(C,d);

% Per cui l'unica soluzione approssimata è

Z

% La distanza tra b e la sua proiezione su Im(A) è

norm(b-A*Z)

% Metodo 2
% Troviamo la decomposizione QR di A

A=sym(A)

% Dobbiamo ortogonalizzare le colonne C1 e C2 di A:

C1=A(:,1)
C2=A(:,2)
F1=C1
a12=sym(C2'*F1)/sym(F1'*F1)
F2=C2-a12*F1

% Adesso normalizziamo F1 ed F2 e troviamo una base ortonormale {E1,E2} di Im(A)

N1=sqrt(F1'*F1)
N2=sqrt(F2'*F2)
E1=1/N1*F1
E2=1/N2*F2

% La matrice Q ha per colonne E1 ed E2

Q=[E1,E2]

% La matrice R è ottenuta dalla matrice [1,a12; 0, 1] moltiplicando la prima riga per N1 e la seconda riga per N2:

R=[N1,N1*a12;0,N2]

% Verifichiamo:

Q*R==A

% Torna! Adesso dobbiamo risolvere il sistema RX=Q'b

d=Q'*b
Z=linsolve(R,d)

% ... e ritroviamo la soluzione ottenuta con il metodo 1.
%
% ESERCIZIO 2.3

syms x y z
eqn1= x+y+z==1
eqn2= 2*x-y-z==0
eqn3= -x+2*y+3*z==2
eqn4= 2*x+2*y+3*z==1
[A,b]=equationsToMatrix(eqn1,eqn2,eqn3,eqn4)

% Vediamo se il sistema è compatibile:

linsolve(A,b)

% Non lo è! Cerchiamo le soluzioni approssimate. Abbiamo visto due metodi:
% 1) Risolvere il sistema A'AX=A'b.
% 2) Trovare la decomposizione QR di A=QR, e risolvere RX=Q'b.
% Intanto vediamo se c'è un'unica soluzione, studiando il rango di A

rank(A)

% Poichè A ha tre colonne, ed il suo rango è tre, ne segue che A'A è invertibile e quindi c'è un'unica soluzione approssimata.
% Metodo 1
% Poniamo

C=A'*A
d=A'*b

% e risolviamo:

Z=linsolve(C,d);

% Per cui l'unica soluzione approssimata è

Z

% La distanza tra b e la sua proiezione su Im(A) è

norm(b-A*Z)

% Metodo 2
% Troviamo la decomposizione QR di A

A=sym(A)

% Dobbiamo ortogonalizzare le colonne C1, C2 e C3 di A:

C1=A(:,1)
C2=A(:,2)
C3=A(:,3)
F1=C1
a12=sym(C2'*F1)/sym(F1'*F1)
F2=C2-a12*F1
a13=sym(C3'*F1)/sym(F1'*F1)
a23=sym(C3'*F2)/sym(F2'*F2)
F3=C3-a13*F1-a23*F2

% Adesso normalizziamo F1, F2 ed F3 e troviamo una base ortonormale {E1,E2, E3} di Im(A)

N1=sqrt(F1'*F1)
N2=sqrt(F2'*F2)
N3=sqrt(F3'*F3)
E1=1/N1*F1
E2=1/N2*F2
E3=1/N3*F3

% La matrice Q ha per colonne E1, E2 ed E3

Q=[E1,E2,E3]

% La matrice R è ottenuta dalla matrice [1,a12, a13; 0, 1, a23; 0,0,1] moltiplicando la prima riga per N1, la seconda riga per N2 e la terza riga per N3:

R=[N1,N1*a12, N1*a13;0,N2, N2*a23; 0, 0, N3]

% Verifichiamo:

Q*R==A

% Torna! Adesso dobbiamo risolvere il sistema RX=Q'b

d=Q'*b
Z=linsolve(R,d)

% ... e ritroviamo la soluzione ottenuta con il metodo 1.
%
% ESERCIZIO 3
% Sia X un vettore di R^n e sia P=Proj(X) la proiezione ortogonale di X su U. Allora X=(X-P)+P e (X-P) è ortogonale ad U.
% Se X appartiene ad U, allora (X-P)*(X-P)=(X-P)*X- P*(X-P); poichè sia X che P appartengono ad U e (X-P) è ortogonale ad U, ne segue che (X-P)*(X-P)=0 e quindi (X-P)=0 e quindi X=P.
% Se X=P, allora X appartiene ad U.
%
% ESERCIZIO 4.1
% Importiamo i quattro dati in due vettori X ed Y, in cui X contiene le ascisse dei dati ed Y le ordinate

X=[-1,0,2,3];
Y=[1,0,3,4];

% Il polinomio di grado 1 che meglio approssima i dati è

syms x
p=polyfit(X,Y,1);
p1=[x,1]*p'

% Il polinomio di grado 2 che meglio approssima i dati è

q=polyfit(X,Y,2);
p2=[x^2, x,1]*q'

% Rappresentiamo graficamente i dati (con dei cerchi ROSSI), la loro approssimazione lineare (in BLU) e quadratica (in VERDE)

x=-1:.1:5;
plot(x,polyval(p,x),'b')
hold on
grid on
plot(X,Y,'or')
plot(x,polyval(q,x),'g')
hold off

%
% ESERCIZIO 4.2
% Importiamo i quattro dati in due vettori X ed Y, in cui X contiene le ascisse dei dati ed Y le loro ordinate:

X=[-1,0,1,3];
Y=[0,2,1,-1];

% Il polinomio di grado 1 che meglio approssima i dati è

syms x
p=polyfit(X,Y,1);
p1=[x,1]*p'

% Il polinomio di grado 2 che meglio approssima i dati è

q=polyfit(X,Y,2);
p2=[x^2, x,1]*q'

% Rappresentiamo graficamente i dati (con dei cerchi ROSSI), la loro approssimazione lineare (in BLU) e quadratica (in VERDE)

x=-2:.1:4;
plot(x,polyval(p,x),'b')
hold on
grid on
plot(X,Y,'or')
plot(x,polyval(q,x),'g')
hold off

%
% ESERCIZIO 5.1
% Importiamo i quattro dati in due vettori X ed Y, in cui X contiene le ascisse dei dati ed Y le loro ordinate:

X=[-1,0,1,3];
Y=[1,-1,3,-4];

% Il polinomio di grado 1 che meglio approssima i dati è

syms x
p=polyfit(X,Y,1);
p1=[x,1]*p'

% Il polinomio di grado 2 che meglio approssima i dati è

q=polyfit(X,Y,2);
p2=[x^2, x,1]*q'

% Il polinomio di grado 3 che meglio approssima i dati è

r=polyfit(X,Y,3);
p3=[x^3,x^2, x,1]*r'

% Si noti che dato che i quattro dati hanno ascissa distinta, il polinomio di grado tre (=4-1) è proprio il polinomio interpolatore dei dati.
% Rappresentiamo graficamente i dati (con dei cerchi ROSSI), la loro approssimazione lineare (in BLU), quadratica (in VERDE) e cubica (in NERO)

x=-2:.1:5;
plot(x,polyval(p,x),'b')
hold on
grid on
plot(X,Y,'or')
plot(x,polyval(q,x),'g')
plot(x,polyval(r,x),'k')
hold off

%
% ESERCIZIO 5.2
% Importiamo i cinque dati in due vettori X ed Y, in cui X contiene le ascisse dei dati ed Y le loro ordinate:

X=[-2,0,1,3,4];
Y=[0,3,2,-1,1];

% Il polinomio di grado 1 che meglio approssima i dati è

syms x
p=polyfit(X,Y,1);
p1=[x,1]*p'

% Il polinomio di grado 2 che meglio approssima i dati è

q=polyfit(X,Y,2);
p2=[x^2, x,1]*q'

% Il polinomio di grado 3 che meglio approssima i dati è

r=polyfit(X,Y,3);
p3=[x^3,x^2, x,1].*r

% Rappresentiamo graficamente i dati (con dei cerchi ROSSI), la loro approssimazione lineare (in BLU), quadratica (in VERDE) e cubica (in NERO)

x=-3:.1:5;
plot(x,polyval(p,x),'b')
hold on
grid on
plot(X,Y,'or')
plot(x,polyval(q,x),'g')
plot(x,polyval(r,x),'k')
hold off

%
% Esercizio 6
% Importiamo i tre dati in due vettori X ed Y, in cui X contiene le settimane (1,2 e 3) ed Y il valore dell'azione nella corrispondente settimana:

X=[1,2,3];
Y=[1.22,1.49,2.45];

% Poichè i tre dati hanno ascissa diversa, esiste un unico polinomio di grado due che li interpola. Esso è dato da:

syms x
p=polyfit(X,Y,2);
p2=[x^2,x,1]*p'

% Se utilizziamo questo polinomio come stima dell'andamento del prezzo dell'azione, otteniamo che alla quarta settimana l'azione vale circa (in EURO)

polyval(p,4)

% ed alla decima settimana (in EURO)

polyval(p,10)

% Supponendo di sapere che alla quarta settimana l'azione valga 3.75 Eur (e non 4.1 Eur, come previsto dalla stima quadratica), allora possiamo trovare un unico polinomio di grado 3 che interpola i quattro dati:

X=[X,4]
Y=[Y,3.75]
q=polyfit(X,Y,3);
p3=[x^3,x^2,x,1]*q'

% Rappresentiamo graficamente i quattro dati (con dei cerchi ROSSI), l'interpolazione dei primi tre dati  (in BLU), e l'interpolazione di tutti e quattro i dati (in VERDE)

x=0:.1:10;
plot(x,polyval(p,x),'b')
hold on
grid on
plot(X,Y,'or')
plot(x,polyval(q,x),'g')

% Vediamo che alla decima settimana il valore dell'azione nel modello cubico (in VERDE) sta scendendo drammaticamente, mentre nel modello quadratico (in BLU) cresce rapidamente.
% In effetti il valore dell'azione alla decima settimana nel modello cubico è

polyval(q,10)

% mentre nel modello quadratico è

polyval(p,10)

% La loro differenza (in valore assoluto) è

abs(polyval(q,10)-polyval(p,10))

% Alla quinta settimana, invece la differenza non è poi cosi considerevole:

abs(polyval(q,5)-polyval(p,5))

%
% ESERCIZIO 9

A=sym([2 1 1; 1 2 1; 1 1 2])

% Calcoliamo il polinomio caratteristico di A, determinando la tr(A), tr(A^2) e det(A)

trace(A)
trace(A^2)
det(A)

% Il polinomio caratteristico di A è quindi

syms x
c(x)=x^3-trace(A)*x^2+1/2*(trace(A)^2-trace(A^2))*x-det(A)

% Si poteva anche usare direttamente il comando (valido per ogni matrice quadrata) charpoly

charpoly(A,x)

% Cerchiamo le radici di c tra i divisori di -4:

c(1)

% Per cui 1 è un autovalore per A. Dividiamo c(x) per (x-1)

simplify(c(x)/(x-1))

% Le radici sono 4 e 1. Per cui c(x)=(x-1)^2(x-4)

simplify(c(x)/(x-1)^2)

% e lo spettro di A è {1,4}. Troviamo una base ortonormale di autovettori di autovalore 1 (dato che A è simmetrica, è diagonalizzabile per il teorema spettrale e quindi 1 ha molteplicità geometrica uguale a due):

B=eye(3)-A
N=null(B)

% per cui una base di V_1(A) è composta dalle colonne di questa matrice N. Ortogonalizziamole:

X1=N(:,1)
X2=N(:,2)
F1=X1
a12=sym(X2'*F1)/sym(F1'*F1);
F2=X2-a12*F1

% Adesso normalizziamo i vettori F1 ed F2 trovati

N1=sqrt(sym(F1'*F1));
N2=sqrt(sym(F2'*F2));
E1=1/N1*F1
E2=1/N2*F2

% {E1,E2} è una base ortonormale dell'autospazio relativo all'autovalore 1. Troviamo una base ortonormale dell'autospazio relativo all'autovalore 4:

D=4*eye(3)-A
M=null(D)
N3=sqrt(sym(M'*M));
E3=1/N3*M

% {E1,E2,E3} è una base ortonormale di R^3 composta di autovettori per A. Consideriamo la matrice ortogonale P che li ha compe colonne:

P=[E1,E2,E3]

% otteniamo quindi una diagonalizzazione ortogonale di A:

P'*A*P