[Calcolo Numerico]Risoluzione sistemi lineari

[Otherguy2k]

Ragazzi ci è stato assegnato un elaborato per la progettazione di una function, in matlab , per la risoluzione di sistemi lineari, mediante:
1)L'algoritmo di forward substitution se il sistema è triangolare inferiore;
2)L'algoritmo di back substitution se il sistema è triangolare superiore;
3)L'algoritmo di Gauss con pivoting parziale virtuale, in seguito al quale si applica back substitution;
Ora tutta la parte sui controlli sull'input l'ho realizzata, cosi come l'algoritmo di forward e back substitution se la matrice data in input è triag inf(sup).
Ho dei problemi per la realizzazione di gauss con pivoting parziale virtuale, quello che in sostanza si deve fare è l'implementazione di un vettore che contenga gli indici delle righe della matrice ,ad esempio detto piv,sul quale effettuo gli scambi richiesti dal pivoting parziale, e quindi al posto di lavorare sugli elmenti A(i,j), lavoro sugli elementi A8piv(i),j), questo per evitare lo spostamento fisico delle righe della matrice.
Ora l'algoritmo io avrei anche realizzata ma il mat quando inserisco la matrice piena mi dice che la mia funzione non torno alcuni argoemnti di output , quindi devo aver sbagliato qualcosa.
Ora posto il codice relativo all'algoritmo di gauss + back substitution da me realizzato,cosi magari potete dirmi se ho fatto qualche errore proprio nella functionc.

    n=length(A);
    piv=1:n;
    for k=1:n-1
        [m,p]=max(abs(A(piv(k:n),k)));
        if abs(A(piv(p),k))<eps*norm(A)
           error('Sistema singolare')
        end
        piv(k)=piv(p);
        piv(p)=piv(k);
        for i=k+1:n
            A(piv(i),k)=A(piv(i),k)/A(piv(k),k);
            for j=k+1:n
                A(piv(i),j)=A(piv(i),j)-A(piv(i),k)*A(piv(k),j);
            end
            B(piv(i))=B(piv(i))-A(piv(i),k)*B(piv(k));
        end   
    end
    if A(piv(n),n)<eps*norm(A)
       error('Errore:Sistema singolare')
    end
    X(piv(n))=B(piv(n))/A(piv(n),n);
    for i=n-1:-1:1
        somma=0;
        for f=i+1:n
            somma=somma+A(piv(i),k)*X(piv(f));
        end
        X(piv(i))=(B(piv(i))-somma)/A(piv(i),i);
    end          
return 

Se c'e qualche parte non chiara ditemelo che la commento , spero che mi possiate aiutare.
Grazie 1000 a ti quelli che eventualmente mi risponderanno
PS:Il codice è venuto senza indentatura perche non so come si fa a postare codice nel formato diciamo originale:(

[Otherguy2k]

Ho rieditato il post sfruttando ,[code] ora è leggibile
Grazie sergio per l'info ^_-

[Otherguy2k]

Nessuno può aiutarmi?

[ilyily87]

potresti postare anche il codice sulla back e sulla forward per piacere?

[Kroldar]

Ho fatto un elaborato simile per l'esame di Calcolo Numerico. Se vuoi darci un'occhiata te lo mando.

[Otherguy2k]

"Kroldar":Ho fatto un elaborato simile per l'esame di Calcolo Numerico. Se vuoi darci un'occhiata te lo mando.

Si grazie magari^^

"ilyily87":potresti postare anche il codice sulla back e sulla forward per piacere?

Certo ora posto proprio tutto il codice allora , ho risolto il problema ma in maniera poco elegante, in pratica ora il mio algoritmo funziona , solo che non mi da le soluzioni nell 'ordine giusto,dunque ho riodidanto il vettore delle soluzioni alla fine , però non mi piace come soluzione è inelegante e brutale , magari qualcuno mi potrebbe consigliare una soluzione migliore?^^
Intanto posto il codice:

   n=length(A);
    piv=1:n;
    for k=1:n-1
        [m,p]=max(abs(A(piv(k:n),k)));
        p=p+(k-1);
        if abs(A(piv(p),k))>eps*norm(A)
            piv(:,[k p])=piv(:,[p k]);
           for i=k+1:n
               A(piv(i),k)=A(piv(i),k)/A(piv(k),k);
               for j=k+1:n
                   A(piv(i),j)=A(piv(i),j)-A(piv(i),k)*A(piv(k),j);
               end
               B(piv(i))=B(piv(i))-A(piv(i),k)*B(piv(k));
           end   
        else
           error('Sistema singolare')
        end
    end
    if abs(A(piv(n),n))<eps*norm(A)
       error('Errore:Sistema singolare')
    end
    X(piv(n))=B(piv(n))/A(piv(n),n);
    for i=n-1:-1:1
        somma=0;
        for f=i+1:n
            somma=somma+A(piv(i),f)*X(piv(f));
        end
        X(piv(i))=(B(piv(i))-somma)/A(piv(i),i);
    end
    X(1:n)=X(piv(1:n));
    end
    X=X';
end

Questo è gauss con pivoting parziale e back substitution "virtuale".

   for j=1:length(A)
       if A(j,j)==0
          error('Errore:matrice incompelta singolare')
       end
   end
   X(1)=B(1)/A(1,1);
   for i=2:length(A)
       somma=0;
       for k=1:i-1
           somma=somma+A(i,k)*X(k);
       end
       X(i)=(B(i)-somma)/A(i,i);
   end
   return

Questo è forward classico

   for j=1:length(A)
       if A(j,j)==0
          error('Errore:matrice incompleta singolare')
       end
   end
   n=length(A);
   X(n)=B(n)/A(n,n);
   for i=n-1:-1:1
       somma=0;
       for k=i+1:n
           somma=somma+A(i,k)*X(k);
       end
       X(i)=(B(i)-somma)/A(i,i);
   end
   return

Infine questo è back classico
Magari poi provo anche a ridurre l'uso dei cicli for però intanto vorrei una soluzione piu elegante a gauss con pivoting per il problema dell'ordine delle soluzioni^^

[Otherguy2k]

Ragazzi ora ho risolto i problemi relativi agli algoritmi,ora funge tutto benissimo ^_^

[lostinyou]

qualcuno potrebbe postarmi il codice dell'algoritmo di gauss con pivoting virtuale parziale che funge io utilizzo come back e forward i seguenti codici

function [x]=forward(A,b)
x=[];
n=length(b);
if prod(diag(A))==0
    error('ERROR: il sistema è singolare');
    
end
x(1)=b(1)/A(1,1);
r=b;
for i=2:n
    r(i:n)=r(i:n)-A(i:n,i-1)*x(i-1);
    x(i)=r(i)/A(i,i)
end

Questa è la forward

function [x]=backsub(A,b)
x=[];
n=length(b);
if prod(diag(A))==0
    error('Il sistema è singolare');
end
x(n)=b(n)/A(n,n);
r=b;
for i=(n-1:-1:1)
    r(1:i)=r(1:i)-A(1:i,i+1)*x(i+1);
    x(i)=r(i)/A(i,i);
end

Questa è la Back substitution

[/quote]

[Axel861]

Anke a me è stato assegnato un elaborato simile ma ci sono alcuni problemi che non riesco a risolvere proprio..non è che qualcuno per farvore può postare l'intero codice della function???