Calcolo numerico - Appunti con tutte le esercitazioni svolte

A

% La struttura della matrice A che si ottiene alla 100-esima iterazione

% è quasi triangolare. Verifichiamo che gli autovalori che si deducono

% dalla matrice finale A siano effettivamente approssimazioni degli

% autovalori di B.

a = zeros(n,1);

a(1) = A(1,1);

a(2:3) = eig(A(2:3,2:3));

err = abs(sort(a)-sort(b))

% La forma base del metodo QR ha dunque fornito un'approssimazione degli

% autovalori della matrice B.

% Ripetiamo lo stesso esercizio per la seguente matrice

B2 = [0 1 0 1; 2 0 4 -2; -1 0 -1 0; -1 2 1 0];

n = size(B2,1);

b = eig(B2);

A = B2;

for m = 1:m_max

[Q,R] = qr(A);

A = R*Q;

end

A

a = zeros(n,1);

a(1:2) = eig(A(1:2,1:2));

a(3) = A(3,3);

a(4) = A(4,4); 44

err = abs(sort(a)-sort(b))

% Anche in questo caso, la forma base del metodo QR ha

% fornito un'approssimazione degli autovalori della matrice B.

Il numero delle iterazioni impiegate dal metodo delle potenze inverse per raggiungere la tolleranza

desiderata dipende dalla scelta di p. Scegliendo un valore di p molto vicino all'autovalore desiderato, il

metodo delle potenze inverse potrebbe impiegare un numero di iterazioni inferiore rispetto al metodo delle

potenze. Ciononostante il metodo delle potenze inverse potrebbe risultare più costoso del metodo delle

potenze in quanto ad ogni passo richiede la risoluzione di due sistemi triangolari. Per esempio, scegliendo

p=5.4, il tempo impiegato dal metodo delle potenze per eseguire 35 iterazioni (sul mio PC) è inferiore a quello

impiegato dal metodo delle potenze inverse per eseguire 12 iterazioni.

clear all

close all

format long e

for n = 5:5:100

A = -ones(n);

A = triu(A,1)+diag(ones(n,1));

[U,S,V] = svd(A);

determinante = det(A)

rango = rank(A)

s = diag(S);

pause

end

Il determinante della matrice A è 1 (prodotto degli elementi della diagonale) e, pertanto, il rango effettivo

della matrice è n. Il rango numerico (che è calcolato a meno di una tolleranza) è invece n-1, a partire da un

certo ordine n della matrice in poi! Tale risultato va interpretato nel seguente modo: al crescere di n la

matrice si avvicina sempre di più a una matrice singolare. Infatti, l'ultimo valore singolare s(n) decresce al

crescere di n, e ciò conferma che la distanza della matrice assegnata dall'insieme delle matrici di rango n-1

diminuisce all'aumentare di n. Si osservi che tale comportamento non si deduce dal valore del determinante,

che è costantemente uguale a 1, ma dal rango numerico ovvero dai valori singolari della matrice. 45

clear all

close all

format short e

A = [3 -2 1 2; -1 0 2 1; 0 5 -6 -1; 1 1 -1 1; 1 -1 -1 -1; 8 -1 -5 2];

b = [1; -3; 7; 0; -6; 2];

r = rank(A); % la matrice ha rango 3

rangoAb = rank([A b]); % la matrice orlata ha rango 4

% Il sistema non ammette soluzione classica e

% si calcola la soluzione ai minimi quadrati.

% Poiché la matrice A non ha rango massimo 4,

% esistono infinite soluzioni ai minimi quadrati

% ma solo una di esse ha norma 2 minima. Tale soluzione

% si ottiene con la decomposizione ai valori

% singolari di A.

[U,S,V] = svd(A);

valori_singolari = diag(S);

ystar = zeros(4,1);

ystar(1:r) = (U(:,1:r)'*b)./valori_singolari(1:r);

xstar = V*ystar;

% La soluzione calcolata si può ottenere anche

% semplicemente con la seguente istruzione.

x = pinv(A)*b;

err = norm(xstar-x)

%

% Si osservi che tutte le soluzioni ai minimi quadrati del sistema

% sono date da xstar+z con z appartenente al Ker(A). Il

% Ker(A) è generato dal quarto vettore colonna di V e, pertanto,

% le soluzioni sono del tipo xstar+c*V(:,4), con c costante. 46

Il condizionamento in norma due si ottiene dal rapporto tra valore singolare massimo e minimo, quindi 20

A=[1 5 0 3 9; 7 8 4 0 1; 2 5 3 9 0; 1 -1 2 1 1; 7 3 -2 0 1];

[Q,R]=qr(A);

for i=1:100

A=R*Q;

[Q,R]=qr(A);

end

A

% Si osserva che la matrice ha un blocco 2x2 diagonale in posizione 3:4

clear all

close all

format long e

n = 100;

% per generare A1 si procede nel seguente modo:

% 1) si genera il vettore p i cui elementi sono n,n-1,,...,2,1

% 2) con il comando p(ones(n,1), :) si ripete n volte il vettore p

% e si genera una matrice le cui n righe sono tutte uguali al 47

% vettore p;

% 3) con il comando diag( ones(n-1,1), -1) si genera una matrice i cui

% elementi sono tutti uguali a zero tranne quelli della codiagonale

% inferiore che sono uguali a 1;

% 4) infine, con triu(A,-1) si estraggono gli elementi di A a partire

% dalla codiagonale inferiore e si pongono uguali a 0 quelli

% sotto la codiagonale inferiore

p = n:-1:1;

A1 = triu( p( ones(n,1), :) - diag( ones(n-1,1), -1), -1 );

% si calcolano gli autovalori di A1

lambda1 = eig(A1);

% e si rappresentano graficamente

plot(real(lambda1),imag(lambda1),'r*','markersize',6)

hold on

% si genera la matrice perturbata A1p, perturbando

% gli elementi dell'ultima riga di A1

A1p = A1;

A1p(n,:) = A1p(n,:)+1.0e-10;

% si calcolano gli autovalori di A1p

lambda1p = eig(A1p);

% e si rappresentano graficamente

plot(real(lambda1p),imag(lambda1p),'ko','markersize',6)

pause

% dal grafico si deduce che il problema è mal condizionato in quanto

% a una piccola perturbazione negli elementi della matrice non è corrisposta

% una perturbazione negli autovalori dello stesso ordine di grandezza

% il cattivo condizionamento è segnalato dal comando condeig:

% le componenti del vettore output, che rappresentano il

% condizionamento di ciascun autovalore, risultano molto maggiore di 1

condizionamento_autovalori_A1 = condeig(A1)

hold off

close

pause

% si ripetono le stesse istruzioni precedenti per la matrice A2

A2 = triu(A1)+triu(A1,1)';

lambda2 = eig(A2);

plot(real(lambda2),imag(lambda2),'r*','markersize',6)

hold on

A2p = A2;

A2p(n,:) = A2p(n,:)+1.0e-10;

lambda2p = eig(A2p);

plot(real(lambda2p),imag(lambda2p),'ko','markersize',6)

pause

condizionamento_autovalori_A2 = condeig(A2)

err_mat_A2 = norm(A2-A2p)

err_eig_A2 = abs(sort(lambda2)-sort(lambda2p))

% In questo caso il problema è ben condizionato (la matrice è simmetrica)

% e gli autovalori della matrice perturbata sono affetti da un errore che è

% dello stesso ordine di grandezza (anzi più piccolo) dell'errore associato 48

% alla matrice perturbata

clear all

close all

z = (1:3)';

A = [0.1 3.8 0; 1 0 0; 0 1 0];

for m_max = 500:300:1100

[lambda,w,m] = potenze_no_norma(A,z,m_max);

w

pause

end

% Se non si normalizzano i vettori iterate si hanno problemi di overflow!

APPUNTI MATLAB

VARIABILI E SCALARI

Nella command window posso scrivere direttamente delle assegnazioni e la classe

di appartenenza viene assegnata automaticamente da matlab

>> a = 6x3

a=18

>> 6x3

ans=18

ans è una variabile temporanea che contiene il risultato più recente

>>whos

mostra variabili, classe appartenenza, memoria, attributo

clear all rimuove variabili

clc rimuove log graphic window

clear a rimuove varaibile a

; sopprime l'output: senza le variabili si stampano a video

… per continuare a scrivere un’istruzione nella riga successiva

disp(variabile) stampa solo il contenuto della variabile

help <nome comando> da informazione comando

lookfor <simil nome comando> cerca nel database

ctrl+c per fermare ricerca di lookfor

Gli scalari vengono salvati come matrici 1x1

VARIABILI PREDEFINITE 49

i, j unità immaginaria

pi pigreco

realmax massimo numero di macchina positivo

realmin minimo numero di macchina positivo

eps epsilon di macchina

Inf ossia un numero maggiore di realmax

∞,

Nan not a number

FORMATI DI OUTPUT

format short fixed point con 4 cifre decimali

format long fixed point con 14 cifre decimali

format short e floating point con 4 cifre decimali

format long e floating point con 15 cifre decimali

format rat frazione irriducibile

FUNZIONI

function [outputArg1,outputArg2] = nome_funzione(inputArg1,inputArg2)

%DESCRIZIONE

CORPO;

end

VETTORI E MATRICI

Dichiarare matrici/vettori:

vet = [1,2,3]

vet=[1:3] vuol dire vettore di interi da 1 a 3

vet=[] genera il vettore vuoto

mat = [1,2,3;4,5,6]

mat[1:3;4:6] matrice di riga uno da 1 a 3 e riga due da 4 a 6

vet=[1 2 3] oppure vet=[1, 2, 3] vettore riga

vet=[1; 2; 3] vettore colonna

vet= [1:10] vettore riga di interi da 1 a 10

vet=[1:3:10] 3 è la lunghezza del passo, se non è multiplo preciso l'ultimo si

tronca. La sintassi generale è v=[valore_iniz:passo:valore_fin]

linspace(x1, x2) crea un vettore di 100 valori compresi tra x1 e x2 equidistanti

linspace(x1,x2, n) n numero di punti equispaziati nell'intervallo x1-x2

zeros(n,1) produce un vettore colonna di lunghezza n con elementi tutti nulli

zeros(1,n) produce un vettore riga di lunghezza n con elementi tutti nulli

ones(n,1) (ones(1,n)) genera un vettore colonna(riga) con tutte le componenti pari

a 1

z=v(3) per accedere alla componente di un vettore, ad es la terza, e assegnare

alla variabile z tale valore 50

Posso dichiarare una mat come:

mat = [1:3;4:6]

oppure equivalentemente:

mat = [1:3;linspace(4,6,3)]

mat = zeroes(3) crea matrice 3x3 di zero

mat = ones(3) crea matrice 3x3 di 1

2*ones(3) matrice 3x3 di 2

mat ones(3,2) crea matrice 3x2 di 1

gli indici di matrici iniziano da 1 e non da 0

vet(1) è il primo

vet(end) ultimo del vettore

vet(end-1) penultima posizione

mat(1,3) elemento nelle cordinate 1,3

mat(1,end) riga 1, ultima colonna

size(v) per controllare la dimensione di una variabile, per una matrice genera un

vettore riga contenente il numero di righe e il numero di colonne della matrice v

length(v) restituisce la lunghezza di un vettore, per matrici equivale a calcolare

max(size(A))

OPERAZIONI SU VETTORI E MATRICI

Mat1+mat2 somma tra matrici

mat1*mat2 prodotto tra matrici

mat' matrice trasposta/ vettore trasposto

mat*vet* aggiungere asterisco per dire che è un vettore colonna

norm(vet) restituisce la norma 2 di vettore

v+w somma algebrica di vettori

v*w oppure dot(v,w) prodotto scalre tra vettori

v.*w prodotto componente per componente

v.^3 elevamento a potenza componente per componente

se ho un vettore lungo 3 e assegno alla posizione 4 un valore, il vettore si

allunga e diventa vettore lungo 4

vet(end+1)=-1 allunga di 1 il vettore e gli assegna il valore -1 all'ultima

posizione

mat(1, end+1)= -1 viene creato un elemento in piu nella riga, ma anche altri sulla

stessa colonna (in cui asse