Esercizio in MATLAB

[deioo]

Calcolo Numerico
C.d.L in Ingegneria Informatica
a. a. 2008/2009
Elaborato 2
1. Scrivere una function Matlab che calcola, a seconda del numero di parametri di
input, l’inversa di una matrice A o la soluzione di un sistema lineare (o m sistemi
lineari con stessa matrice dei coefficienti) Ax=b
utilizzando :
• L’algoritmo di fattorizzazione LU con memorizzazione in place e il pivoting
parziale con scambio virtuale delle righe, e gli algoritmi di forward substitution
e back substitution.
Per lo sviluppo dell’elaborato si seguano le seguenti specifiche:
• parametri di INPUT:
A, matrice quadrata di reali di dimensione n×n
b, facoltativo, vettore di dimensione n(matrice di dimensione n×m), termine
noto del sistema Ax=b,( termini noti di m sistemi con A matrice dei coefficienti).
• parametri di OUTPUT:
Nel caso di un parametro di input:
x, matrice quadrata di reali di dimensione n×n, inversa di A.
Nel caso di due parametri di input:
x, vettore di dimensione n, soluzione del sistema (matrice di dimensione n×m,
le cui colonne sono le soluzioni degli m sistemi)
La function deve arrestarsi con un messaggio di errore per l’utente nel caso di
singolarità del sistema.

qualcuno potrebbe aiutarmi a fare questo esercizio? sn veramente negato in matlab grazie

[Luc@s]

piccolo hint

%fattorizzazione lu senza pivot

function [L U b]= miolu(A,b)

if size(A,1)== size(b,1)

    A=[A b];

else 

    error('Il numero di righe di b non e pari a quello di A')

end

[m n]=size(A);

for k=2:m

    if A(k-1,k-1)==0

        error(['Il minore principale di ordine ' num2str(k-1) 'e''nullo']) % num2str significa "number to strig"

    end

    A(k:end, k-1)=A(k:end,k-1)./A(k-1,k-1);

    for j=k:n

        A(k:end,j)=A(k:end,j)-A(k:end,k-1).*A(k-1,j);

    end

end

if A(m,m)==0

    error(['Il minore principale di ordine ' num2str(m) 'e''nullo']) 

end

L=eye(m);

L=L+tril(A(1:m,1:m),-1);

U=triu(A(1:m,1:m));

b=A(:,m+1:n);

[deioo]

un mio amico l ha fatto cosi voi ke ne pensate?

function x=gausspiv(A,b)
 
[n,m]=size(A);
 
if n~=m
    error('La matrice A deve essere quadrata [n x n]');
end
 
if nargin == 1
    b=eye(n);
end
 
[l,p]=size(b); %#ok<NASGU>
 
 
if n~=l
    error('Dimensione di b non valida');
end
 
 
piv=1:n;
zero=eps(norm(A,inf)); % calcolo dello zero macchina per A
 
for k = 1 : n-1
    
    r = min(find(abs(A(piv(k:n),k))==max(abs(A(piv(k:n),k))))+ (k-1));
    
    if abs(A(piv(r),k)) > zero    
        if r~=k                 
            piv([k r])=piv([r k]);  
        end 
            for z = k+1 : n
                 A(piv(z),k) = A(piv(z),k) / A(piv(k),k);            
                 for s=k+1 : n
                    A(piv(z),s)=A(piv(z),s) - A(piv(z),k) * A(piv(k),s);
                 end
            end       
    else
        error('Il sistema è singolare');
    end 
end
 
if abs(A(piv(n),n)) < zero 
        error('Il sistema è singolare');                                                      
end                                                                                       
%Algoritmo del forward substitution  Ly=Pb
 
y(1,:)=b(piv(1),:);
for t=2:n      
    som=0;     
    for k=1:t-1
        som=som+A(piv(t),k) .* y(k,:);
    end
    y(t,:) = b(piv(t),:) - som;
end
 
%Algoritmo del back sostitution Ux=y
 
x(n,:)=y(piv(n),:) / A(piv(n),n); 
 
for t=n-1:-1:1 
    som=0;
    for k = t+1:n
        som=som+A(piv(t),k).*x(k,:);
    end
    x(t,:)= (y(t,:) - som)/A(piv(t),t);
end