L’algoritmo di Gauss naive per la risoluzione di sistemi lineari Ax=b trasforma il sistema in uno triangolare superiore. Scrivi un algoritmo in Matlab che risolva il sistema lineare Ax=b operando come segue
a) generi a caso un problema del tipo Ax=b di dimensione n la cui soluzione esatta xx è nota;
b) applichi l’algoritmo 1 per triangolarizzare il sistema;
c) trovi la soluzione approssimata x risolvendo il sistema triangolare superiore;
c) stampi la soluzione trovata e la norma euclidea dell’errore xx-x.
Sperimenta l’algoritmo con n=10, 20,40, 60, 80,100.
Suggerimento: Il problema a caso può essere determinato, generando una matrice A casuale di dimensione n i cui elementi sono tra 0 e 1. Si sceglie, ugualmente a caso, il vettore colonna xx di dimensione n e quindi si calcola b=A*xx. n, A e b sono i dati del problema; xx è la soluzione esatta del problema.
Modifica il metodo di Gauss naive in modo che sia applicabile a problemi Ax=b in cui det(A) ≠0 e qualche minore principale sia uguale a zero.
Suggerimento: Il problema a caso può essere determinato come nel caso precedente ed imponendo successivamente che due righe di qualche matrice principale siano uguali.
io ho provato in 3 modi ma nn so quale sia quello giusto
modo1
- Codice:
n=1;
clc % ripulisce lo schermo
while(n)
clear % ripulisce lo Workspace
n=input('Inserisci la dimensione della matrice dei coefficienti, premi ''invio'' per uscire:');
if(n)
clc
A=rand(n); % genera la matrice A
% Annullo il determinante di un minore estratto
% imponendo a 0 la riga 2 fino a l'indice n-1.
for i = 1:n-1
A(2,i) = 0;
end
C=A; % memorizzo in C la matrice originale
xx=rand(n,1); % soluzione esatta
b=A*xx; % genera il vettore dei termini noti
A=[A,b]; % matrice n*(n+1) contenente gli elementi di A e i termini noti
% Algoritmo per la triangolarizzazione superiore
for k = 1:(n-1)
% M è il valore max dei rapporti, p indice relativo a k in cui si trova M
[M,p] = max(abs(A(k:n,k)));
% Ristabilisce l'indice p rispetto alla matrice completa
p = p + k - 1;
if (k ~= p)
% scambio la riga p con la riga k
B = A(k,k:n+1);
A(k,k:n+1) = A(p,k:n+1);
A(p,k:n+1) = B;
end
for i = (k+1):n
m = A(i,k) / A(k,k); % moltiplicatore
for j = (k+1):(n+1)
A(i,j) = A(i,j) - m * A(k,j);
end
end
end
% Algoritmo per la soluzione approssimata x del sistema triangolare superiore
x(n)=A(n,n+1)/A(n,n);
for i=n-1 : -1 : 1
x(i)=A(i,n+1);
for j=i+1 : n
x(i)=x(i)-A(i,j)*x(j);
end
x(i)=x(i)/A(i,i);
end
% Stampa la soluzione
disp(sprintf('La soluzione approssimata del sistema di dimensione %g è x=',n))
disp(x')
disp(sprintf('la norma euclidea dell''errore è = %g',norm(xx-x')))
end
end
modo 2
- Codice:
function A = gaussNaivePivot(A)
% L'algoritmo di Gauss Naive con Pivotting Parziale (scambio di righe)
% consente di risolvere un sistema lineare del tipo Ax=b trasformandolo
% in un sistema triangolare superiore.
% Lo scambio delle righe (pivotting) rende sempre applicabile
% l'argoritmo Gauss Naive se det(A)!=0.
[n, col] = size(A); % numero di righe e colonne
for k = 1:(n-1)
% M è il valore max dei rapporti, p indice relativo a k in cui si trova M
[M,p] = max(abs(A(k:n,k)));
% Ristabilisce l'indice p rispetto alla matrice completa
p = p + k - 1;
if (k ~= p)
% scambio la riga p con la riga k
B = A(k,k:n+1);
A(k,k:n+1) = A(p,k:n+1);
A(p,k:n+1) = B;
end
for i = (k+1):n
m = A(i,k) / A(k,k); % moltiplicatore
for j = (k+1):(n+1)
A(i,j) = A(i,j) - m * A(k,j);
end
end
end
% === EOF =================================================================
modo 3
- Codice:
disp('Risoluzione Ax=b con Gauss Naive Pivotting');
n = input('Inserisci la dimensione del sistema: ');
if isnumeric(n) & (n > 0)
% calcolo una nuova matrice finchè il det(A) = 0 e finchè l'utente
% vuole continuare
reply = 's'; % risposte utente
while strcmpi(reply, 's')
reply = 'n';
disp('genero una matrice casuale');
A = rand(n);
% Annullo il determinante di un minore estratto
% imponendo a 0 la riga 2 fino a l'indice n-1.
for i = 1:n-1
A(2,i) = 0;
end
disp(A);
if (cond(A) >= 1/eps)
% chiede all'utente se generare una nuova matrice
disp('la matrice è mal condizionata');
reply = input('genero una nuova matrice? S/N [N]: ','s');
else
disp('genero la soluzione esatta del problema');
xx = rand(n,1)
disp('calcolo il vettore dei termini noti');
b = A * xx
disp('triangolarizzo il sistema');
A = gaussNaivePivot([A,b])
disp('calcolo la soluzione approssimata');
x = triangSup([A,b])
% Stampa l'errore di approssimazione con la norma euclidea
disp(sprintf('L''errore commesso è (norma euclidea): %d', norm(xx-x)));
end
end
else
disp('ERRORE NEI DATI DI INPUT!');
end
% === EOF =================================================================
Ho l'esame domani.. se qualcuno mi può dare una mano m farebbe davvero piacere
ciaoo
