Codici matlab per l'esame di Calcolo numerico

QUADRATI ECC.)

% CONDIZIONAMENTO DELLA MATRICE: ESEMPIO 1

% In un problema ben condizionato (cond(A) piccolo):

% piccole variazioni sui dati provocano piccole variazioni sui risultati.

% In un problema mal condizionato (cond(A) grande):

% piccole variazioni sui dati possono provocare grandi variazioni sui risultati.

% ESEMPIO 1: PROBLEMA BEN CONDIZIONATO

% Ax=b

A=[2 0.3 0 0.1 0

0 2 0 0 0

0 0 2 0 0.1

0.2 0 0 1 0

0 0 -0.1 0 1];

b1=[1 2 3 4 5]';

x=ones(5,1);

c=cond(A); % =2.2: relativamente basso, il problema è ben condizionato

% Soluzione del problema con Jacobi (A è diag dom)

[xsol_j,ifail_j,iter_j,res_j] = jacobi_matrix(A,b1,x,100,1e-3,1e-5);

% xsol=[0.152, 1; 1.244, 3.970, 5.124]

% Soluzione del problema con Gauss-Seidel

[xsol_g,ifail_g,iter_g,res_g] = gs_matrix(A,b1,x,100,1e-3,1e-5);

% xsol=[0.152, 1; 1.244, 3.970, 5.124]

% Perturbiamo i dati e vediamo come varia il risultato

db=[0.1 0.3 -0.5 0.8 -0.3 ]'; % errore sul dato

b2= b1 + db;

[xsol_j2,ifail_j2,iter_j2,res_j2] = jacobi_matrix(A,b2,x,100,1e-3,1e-5);

% xsol=[0.139, 1.15, 1.01, 4.77, 4.8]

[xsol_g2,ifail_g2,iter_g2,res_g2] = gs_matrix(A,b2,x,100,1e-3,1e-5);

% xsol=[0.139, 1.15, 1.01, 4.77, 4.8]

% Il problema è ben condizionato: l'errore relativo sul risultato è rimasto

% piccolo: Ex < 2.2 * Eb

dx= xsol_g2 - xsol_g; % errore sul risultato

Ex= (norm(dx))/(norm(xsol_g)); % errore relativo sul risultato: 0.1361

Eb= (norm(db))/(norm(b1)); % errore relativo sul dato: 0.1401

% L'errore sul risultato è più piccolo dell'errore sul dato

% CONDIZIONAMENTO DELLA MATRICE: ESEMPIO 2

% In un problema ben condizionato (cond(A) piccolo):

% piccole variazioni sui dati provocano piccole variazioni sui risultati.

% In un problema mal condizionato (cond(A) grande):

% piccole variazioni sui dati possono provocare grandi variazioni sui risultati.

% ESEMPIO 1: PROBLEMA MAL CONDIZIONATO

% Ax=b

% Adesso costruisco una matrice simmetrica e definita positiva con

% autovalori prefissati, in modo che:

% il metodo di Gauss-Seidel converga;

% il condizionamento sia dato dal rapporto tra il massimo e il minimo

% autovalore (e voglio una matrice mal condizionata)

% Utilizzo la funzione fatta da me:

d=[1 31 53 24 571]; % autovalori di A

[A]=mat_eig_sp(d); % costruzione di A

b1=[1 2 3 4 5]';

x=ones(5,1);

c=cond(A); %=571: alto, il problema è mal condizionato

% Soluzione del problema con Gauss-Seidel

[xsol_g,ifail_g,iter_g,res_g] = gs_matrix(A,b1,x,1000,1e-3,1e-8);

% xsol=[-0.573;1.213;0.081;0.086;0.075]

% Perturbiamo i dati e vediamo come varia il risultato

db=[-0.5 2.3 -0.7 0.9 0.8 ]'; % errore sul dato

b2= b1 + db;

[xsol_g2,ifail_g2,iter_g2,res_g2] = gs_matrix(A,b2,x,1000,1e-3,1e-8);

% xsol=[-1.58;3.26;0.056;0.113;0.083]

% Il problema è mal condizionato: l'errore relativo sul risultato è

% diventato grande, aumentando di 5 volte rispetto all'errore sul dato

dx= xsol_g2 - xsol_g; % errore sul risultato

Ex= (norm(dx))/(norm(xsol_g)); % errore relativo sul risultato: 1.69

Eb= (norm(db))/(norm(b1)); % errore relativo sul dato: 0.37

% L'elevato condizionamento della matrice avrebbe potuto portare ad un

% errore molto più grande di quello ottenuto in quest'esempio, fino a 3

% ordini di grandezza in più (c ha 3 ordini di grand in più rispetto ad Eb).

% Nel nostro caso si è comunque registrato un errore sul risulato che è di

% 1 ordine di grandezza maggiore rispetto all'ordine sul dato.

% COSTRUZIONE DI UNA MATRICE SIMMETRICA E DEFINITA POSITIVA

A=input("Enter a matrix: ");

AS=A*A';

% COSTRUZIONE DI UNA MATRICE CON AUTOVALORI PREFISSATI

d=input("Enter a row vector that contains the eigevalues: ");

[A]=mat_eig(d);

% Posso usare questa funzione per creare una matrice con tutti

% autovalori positivi, e quindi definita positiva

% COSTRUZIONE DI UNA MATRICE CON AUTOVALORI PREFISSATI, SIMMETRICA E DEFINITA POSITIVA

d=input("Enter a row vector that contains the eigevalues: ");

% La matrice sarà definita positiva se tutti gli autovalori sono positivi

[AM]=mat_eig_sp(d);

% PROBLEMA DI MINIMI QUADRATI: ESEMPIO 1 (Sistema Sovradeterminato)

% In un sistema sovradeterminato (numero di equazioni superiore al numero

% di incognite) la soluzione si può ottenere risolvendo un'equazione normale.

% In quest'esempio si considera una serie di dati x, da cui dipende una

% seconda serie di dati y. Si vuole trovare la funzione lineare che meglio

% approssima la dipendenza di queste serie di dati.

% Generazione di dati esempio

x = linspace(0, 10, 10)'; % 10 punti dati

y = 2 * x + 3 + rand(size(x)); % y = 2x + 3 + un certo disturbo casuale

% Costruzione della matrice del sistema

X = [x, ones(size(x))];

% In questo modo abbiamo ottenuto un sistema del tipo: Xb=y, dove X (10x2)

% e y (10x1) sono noti, e b (2x1) rappresenta i coefficienti incogniti del

% polinomio di primo grado che vogliamo trovare

% Risoluzione con i minimi quadrati: (X'X)b_min=X'y equazione normale

b = (X' * X) \ (X' * y);

% Calcolo dei valori della funzione lineare

y2 = X * b; % combinazione lineare delle colonne di X

% Visualizzazione dei risultati

figure;

plot(x, y, 'ro', 'DisplayName', 'Dati Originali'); % Dati originali

hold on;

plot(x, y2, 'b-', 'LineWidth', 1.5, 'DisplayName', 'Retta dei Minimi Quadrati'); % Polinomio

approssimante

legend;

title('Approssimazione con il Metodo dei Minimi Quadrati');

xlabel('x');

ylabel('y');

grid on;

hold off;

% PROBLEMA DI MINIMI QUADRATI: ESEMPIO 2 (Sistema Sovradeterminato con Fattorizzazione di

Cholesky)

% In un sistema sovradeterminato (numero di equazioni superiore al numero

% di incognite) la soluzione si può ottenere risolvendo un'equazione normale.

% In quest'esempio si considera una serie di dati x, da cui dipende una

% seconda serie di dati y. Si vuole trovare la funzione lineare che meglio

% approssima la dipendenza di queste serie di dati.

% L'equazione normale è risolta sfruttando la FATTORIZZAZIONE DI CHOLESKY.

% Generazione di dati esempio

x = linspace(0, 10, 10)'; % 10 punti dati

y = 2 * x + 3 + rand(size(x)); % y = 2x + 3 + un certo disturbo casuale

% Costruzione della matrice del sistema

X = [ones(size(x)), x];

% In questo modo abbiamo ottenuto un sistema del tipo: Xb=y, dove X (10x2)

% e y (10x1) sono noti, e b (2x1) rappresenta i coefficienti incogniti del

% polinomio di primo grado che vogliamo trovare

% Risoluzione con i minimi quadrati usando la fattorizzazione di Cholesky

S = X' * X;

c = X' * y;

% Fattorizzazione di Cholesky

R = chol(S); % A=R'R con R triang sup

% Risoluzione del sistema lineare usando la fattorizzazione di Cholesky

% Primo passo: R' * ystar = c

ystar = R' \ c;

% Secondo passo: R * b = ystar

b = R \ ystar;

% Calcolo dei valori della funzione lineare

y2 = X * b; % combinazione lineare delle colonne di X

% Visualizzazione dei risultati

figure;

plot(x, y, 'ro', 'DisplayName', 'Dati Originali'); % Dati originali

hold on;

plot(x, y2, 'b-', 'LineWidth', 1.5, 'DisplayName', 'Retta dei Minimi Quadrati'); % Polinomio

approssimante

legend;

title('Approssimazione con il Metodo dei Minimi Quadrati');

xlabel('x');

ylabel('y');

grid on;

hold off;

% RAGGIO SPETTRALE (Metodo di Gauss-Seidel): ESEMPIO 1

% Condizione necessaria (e sufficiente) affinché un metodo iterativo del tipo

% x(k+1)=Tx(k)+c converga è che il raggio spettrale della matrice di

% iterazione T sia < 1.

% In questo esempio la matrice dei coefficienti è presa in modo che il raggio

% spettrale della matrice di iterazione per il metodo di G-S [T=(D-L)^(-1)*U] sia > 1.

% Ci aspettiamo quindi che il metodo non converga.

A = [

3 -1 -0.5 0.2 0.8;

0 4 -1 0.3 0.5;

5 2 3 -0.2 0.1;

-4 3.5 0 2 -0.7;

1 0.7 0 -8 1

];

D=diag(diag(A)); % Matrice diagonale

L=-(tril(A)-D); % Triangolare inferiore (tutta nulla) cambiata di segno

U=-(triu(A)-D); % Triangolare superiore cambiata di segno

T= ((D-L)^(-1)*U); % Matrice di iterazione (G-S)

% Calcolo il raggio spettrale di T

rho=max(abs(eig(T))); % rho=1.08

% Il raggio spettrale è > 1 quindi non avremo convergenza

% Definisco termine noto e punto iniziale

b=[1 2 3 4 5]';

x=ones(5,1);

[xsol,ifail,iter,res] = gs_scalar(A,b,x,500,1e-3,1e-5);

% ifail=1, iter=500, è stato raggiunto il massimo numero di iterazioni,

% senza però arrivare alla condizione di arresto sul residuo: il metodo non

% è andato a convergenza.

%--------------------------------------------------------------------------

% In questo esempio la matrice dei coefficienti è presa in modo che il raggio

% spettrale della matrice di iterazione [T=((D-L)^(-1)*U)] sia < 1.

% Ci aspettiamo quindi che il metodo converga

A2 = [

3 -1 -0.5 0.2 0.8;

0 4 -1 0.3 0.5;

5 2 3 -0.2 0.1;

0 3 0 2 -0.7;

1 0.7 0 0 1

];

D2=diag(diag(A2)); % Matrice diagonale

L2=-(tril(A2)-D2); % Triangolare inferiore (tutta nulla) cambiata di segno

U2=-(triu(A2)-D2); % Triangolare superiore cambiata di segno

T2= ((D2-L2)^(-1)*U2); % Matrice di iterazione (G-S)

% Calcolo il raggio spettrale di T2

rho2=max(abs(eig(T2))); % rho2=0.3264

% Il raggio spettrale è < 1 quindi avremo convergenza

% Definisco termine noto e punto iniziale

b2=[1 2 3 4 5]';

x2=ones(5,1);

[xsol2,ifail2,iter2,res2] = gs_scalar(A2,b2,x2,500,1e-3,1e-5);

% ifail2=0, iter=9 : il metodo converge, come previsto

% è anche possibile prevedere il numero di iterazioni: iter>=k, dove k:

k=(log10(1e-3))/(log(rho2))

% RAGGIO SPETTRALE (Metodo di Gauss-Seidel): ESEMPIO 2

% Confronto velocità di convergenza del metodo di G-S al variare del raggio spettrale.

% La matrice dei coefficienti è presa in modo che il raggio spettrale della

% matrice di iterazione [T=((D-L)^(-1))*U] sia < 1.

A = [

3 -1 -0.5 0.2 0.8;

0 4 -1 0.3 0.5;

2 0 3 -0.2 0.1;

0 3 0 2 -0.7;

1 0.7 0 0 1

];

D=diag(diag(A)); % Matrice diagonale

L=-(tril(A)-D); % Triangolare inferiore (tutta nulla) cambiata di segno

U=-(triu(A)-D); % Triangolare superiore cambiata di segno

T=((D-L)^(-1))*U; % Matrice di iterazione (G-S)

% Calcolo il raggio spettrale di T

rho=max(abs(eig(T))); % rho=0.2462

% Il raggio spettrale è < 1 quindi avremo convergenza

% Definisco termine noto e punto iniziale

b=[1 2 3 4 5]';

x=ones(5,1);

[xsol,ifail,iter,res] = gs_scalar(A,b,x,100,1