[MATLAB] quali tra questi algoritmi è giusto??

[gandelf]

Salve a tutti, ho questo problema:

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

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

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

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

[gandelf]

UP!


ps la parte dove si differenziano i 3 algoritmi è questa:

% M è il valore max dei rapporti, p indice relativo a k in cui si trova M
[M,p] = max(abs(A(k:n,k)));

[gandelf]

ç_ç