Esercitazioni matlab

Laboratorio 3

3n=20; A=zeros(n,n); b=zeros(n,1); % è vettore colonna for i=1:n for j=1:n if i==j A(i,j)=-1; else if i==j+1 A(i,j)=1; else A(i,j)=0; end end end end for i=1 for j=1:n A(1,j)=1; end end A b(1)=2; for i=2:n b(i) = 0; end

Calcolo della fattorizzazione LU

Si calcoli la fattorizzazione LU della matrice A, mediante la funzione Matlab lu. Verificare che la tecnica del pivoting non è stata usata in questo caso.

[L,U,P] = lu(A) Non c'è pivoting in quanto P mi viene la matrice identità, e se ci fosse pivoting non dovrebbe venire matrice identità.

Algoritmo di sostituzione in avanti

function [x] = fwsub(L,f)

% algoritmo sostituzioni in avanti [x] = function fwsub(L,f)

n=length(f); x= zeros(n,1); % così ci viene vettore colonna

% la matrice deve essere triangolare superiore

if (L==tril(L))

% La matrice deve essere non singolare. Essendo triangolare, i suoi autovalori si trovano sulla diagonale principale

if (prod(diag(L)) == 0) % almeno un elemento diagonale nullo

error('Matrice singolare')

end

x(1)= f(1)./L(1,1);

for i=2:n

sum = 0; % ogni volta che cerca soluzione la somma deve resettarsi

for j=1:i-1

sum= sum + L(i,j).*x(j);

x(i)= (f(i) - sum )./ L(i,i);

end

end

else error('Matrice non triangolare superiore')

end

Algoritmo di sostituzione all'indietro

function [x] = bksub (U,f)

% algoritmo sostituzioni all'indietro function [x] = bksub (U,f)

n=length(f); x= zeros(n,1);

% la matrice deve essere triangolare inferiore

if (U==triu(U))

% La matrice deve essere non singolare. Essendo triangolare, i suoi autovalori si trovano sulla diagonale principale

if (prod(diag(U)) == 0) % almeno un elemento diagonale nullo

error('Matrice singolare')

end

x(n)= f(n)./ U(n,n);

for i= n-1:-1:1

sum=0;

for j= i+1:n

sum=sum + U(i,j).*x(j);

x(i)= (f(i) - sum)./ U(i,i);

end

end

else error('Matrice non triangolare superiore')

end

Risoluzione numerica

Risolvere numericamente, utilizzando le funzioni fwsub.m e bksub.m implementate al punto precedente, i due sistemi triangolari necessari per ottenere la soluzione del sistema di partenza Ax =b.

y= fwsub (L,P*b)

x = bksub (U,y)

questo perché Ly=Pb e Ux=y

k=cond(A) numero condizionamento Ax_ex= 2/n*ones(n,1);

norma_err_rel=norm((x_ex-x),2)./norm(x_ex,2)

norma_residuo= (norm(b-A*x,2)./norm(b,2))

Valutazione del numero di condizionamento

N=[10 20 40 80 160];

k=[];

norma_err_rel=[];

norma_residuo=[];

for n=N

A=zeros(n,n); b=zeros(n,1); % è vettore colonna

for i=1:n

for j=1:n

if i==j

A(i,j)=-1;

else if i==j+1

A(i,j)=1;

else

A(i,j)=0;

end

end

end

end

for i=1

for j=1:n

A(1,j)=1;

end

end

b(1)=2;

for i=2:n

b(i) = 0;

end

[L,U,P] = lu(A);

y = fwsub (L,P*b);

x = bksub (U,y);

x_ex= 2/n*ones(n,1);

norma_err_rel=[norma_err_rel ; norm((x_ex-x),2)./norm(x_ex,2)];

norma_residuo= [norma_residuo; (norm(b-A*x,2)./norm(b,2))];

k= [k; cond(A)];

end

semilogy ( N , norma_err_rel , 'm' , N, norma_residuo , 'k' , N ,k,'g')

legend( 'errore rel .' , 'residuo norm .' , 'n. di condizionamento ')

xlabel ( 'dimensione n' )

ylabel ( 'err , r, K' )

grid on

Calcolo con la matrice di Hilbert

n=5; H= hilb (n) x_ex=ones(n,1) b=H*x_ex R=chol(H)

y = fwsub ( R' , b)

x= bksub (R, y)

% b l'abbiamo già trovato prima

[L,U,P]= lu(H) P = 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0

y= fwsub ( L ,P*b) x= bksub (U,y)

% permutazione tra seconda e terza riga K= cond(H) err_rel= norm(x_ex-x)./norm(x_ex,2) residuo = norm (b-H*x,2) ./ norm (b,2)

N=[3:15]; K=[]; err_rel=[]; residuo=[];

for n=N

H= hilb (n); x_ex=ones(n,1); b=H*x_ex; [L,U,P]= lu(H); y= fwsub ( L ,P*b); x= bksub (U,y); % permutazione tra seconda e terza riga K= [K ; cond(H)]; err_rel= [err_rel ; norm(x_ex-x , 2)./norm(x_ex , 2)]; residuo= [residuo; (norm(b-H*x,2)./norm(b,2))]; end

semilogy (N,err_rel, 'm' , N , residuo , 'k' , N , K , 'g') legend( 'errore rel .' , 'residuo' , 'n. di condizionamento ' ) xlabel ('dimensione n') ylabel ('errore relativo, residuo , K') grid on

Laboratorio 4

A=zeros(7,7); b=[7;4;5;5;5;4;7];

A= A + diag(9*ones(7,1))+ diag(-3*ones(6,1),+1)+diag(-3*ones(6,1),-1)+ diag (ones(5,1),+2)+ diag (ones(5,1),-2)

nnz(A) % risultato è 29

flag=0; for i=1:7 if (sum(A(i,:)) - A(i,i)) < A(i,i) flag= flag +1 ; end end if flag==0 display ('no dominanza stretta') else display (' dominanza stretta') end

% A'=A quindi simmetrica L=ones(7,1); if L'*A*L > 0 display ('matrice def positiva') else display ('matr non def positiva') end % esce matrice definita positiva alfa= 2 / max(eig(A)) alfa_opt = 2 / (max(eig(A)) + min(eig(A))) alfa = 1.246860968868881e-01 alfa_opt = 9.549043002954583e-02

Metodo di Richardson

function [x,k] = richardson(A,b,x0,alpha,toll,nmax)

% [x,k] = richardson(A,b,x0,alpha,toll,nmax) La funzione restituisce la soluzione numerica x e il numero di iterazioni effettuate k.

x = x0 ; k=0;

while k < nmax && norm ( b - A*x) / norm (b) > toll

k= k+1;

x = x + alpha * ( b - A * x);

end

end

[x,k] = richardson(A,b,zeros(7,1), alfa_opt,1e-6,1000) x =1.00001.00001.00001.00001.00001.00001.0000 k =22

Metodo di Jacobi e Gauss-Seidel

D=zeros(7,7); D = D + diag( diag(A) ) E = -tril(A , -1); F = -triu(A , 1); B_J = inv(D) * (D-A) B_GS = inv ( D - E) * F % devo calcolare i raggi spettrali raggio_spettrale_BJ = max (abs((eig(B_J)))) raggio_spettrale_BGS = max(abs(eig(B_GS))) →Entrambi <1 condizione okay

function [x, k] = jacobi(A, b, x0, toll, nmax)

% [x, k] = jacobi(A, b, x0, toll, nmax) D = zeros ( size (A)); D = D + diag ( diag(A)); x = x0; k = 0;

while k < nmax && norm(b - A*x)/norm (b)> toll

k = k+1;

x = ( eye(size(A)) - inv(D) * A ) * x + inv(D) * b;

end

end

D = diag ( diag(A)); E = -tril(A,-1) ; F = -triu(A,1) ; x = x0; k=0; while k < nmax && norm(b - A*x) / norm(b) > toll k = k+1; x = inv ( D - E ) * F + inv(D - E) * b; end

end

[x, k] = jacobi(A, b, zeros(7,1) , 1e-6 , 1000) [x, k] = gs(A, b, zeros(7,1), 1e-6 , 1000)

Esercizio 2

A = zeros(100); A = A + diag ( -2 * ones(100,1)) + diag (ones(99,1) , -1); A(1,:) = ones( 1 ,100 ) nnz (A) ans = 298

spdiags ( [1,-2].*ones(100,1) , [-1, 0], 100, 100);

>> full(ans)

D = diag (diag(A)); E = tril(A) - diag(A); F = triu(A) - diag(A); B_J = inv(D) * (D - A); B_gs = inv(D - E) * F; ragg_spet_B_J = max(abs(eig(B_J))) ragg_spett_B_gs= max(abs(eig(B_gs))) [x, k] = jacobi(A, b, zeros(100,1) , 1e-6 , 1000) [x, k] = gs(A, b, zeros(100,1) , 1e-6 , 1000)

Laboratorio 5

A = zeros(50,50) + diag( 4* ones(50,1)) + diag( -1* ones(49,1),+1) + diag( -1* ones(49,1), -1) + diag( -1* ones(48,1), +2)+diag( -1*ones(48,1), -2); b = 0.2*ones(50,1); x0 = zeros(50,1);

if A==A' display('matrice simmetrica') else display ('matrice non simmetrica') end L = ones(50,1); if L'*A*L > 0 display ('matrice definita positiva') else display('matrice non positiva') end K_A = cond(A); ( 336.2412 )

function [x, iter,err] = graddyn(A, b, x0, nmax, tol)

% [x, iter, err] = graddyn(A, b, x0, nmax, tol) x = x0; iter = 0; r= b - A*x; err=[];

while iter < nmax && norm(b - A*x)/ norm(b) > tol

iter = iter + 1;

alpha = ((r)' * r) / ((r)' * A * r);

x = x + alpha * r;

r = r - alpha*A*r;

err = [err ; norm(r)];

end

end

function [x, iter, err] = gradprec(A, b, P, x0, nmax, tol)

% [x, iter, r] = gradprec(A, b, P, x0, nmax, tol) x = x0; r = b - A * x; iter = 0; err=[];

while iter < nmax && norm (b - A*x) / norm (b) > tol

iter = iter + 1;

z = P \ r;

alpha = (z' * r) / (z' * A * z );

x = x + alpha * z;

r = r - alpha * A * z;

err = [err ; norm(r)];

end

end

P = zeros ( 50 , 50 ) + diag ( 2 * ones (50,1)) + diag ( -1 * ones(49 , 1), +1) + diag ( -1 * ones ( 49 ,1) , -1); figure semilogy ( [1:iter_prec], err_prec , 'b-' , [1:iter_nonprec] ,err_non , 'g-', 'linewidth', 2) legend ( 'precondizionato' , 'non precondizionato') xlabel (' n iter') ylabel (' residuo normalizzato') T = zeros (47 , 47) + diag ( 2 * ones(47,1)) + diag ( -1 *ones(46,1), +1) + diag ( -1 * ones(46,1), -1); F = zeros (47 , 47) + diag ( 6 * ones(47,1)) + diag ( -4 *ones(46,1), +1) + diag ( -4 * ones(46,1), -1) + diag ( ones(45,1), +2) + diag ( ones(45,1) , -2); h = 1/50; A = T ./ (h^2) + F ./ (h^4); L = ones(47,1); if A==A' display('matrice simmetrica') else display ('matrice non simmetrica') end if L'*A*L > 0 display ('matrice definita positiva') else display('matrice non positiva') end k_A = cond(A);