Calcolo numerico - Funzione sislin

ACCURATEZZA:

L'accuratezza dell'algoritmo è massima, nel caso in cui la matrice A dei coefficienti sia una matrice ben condizionata. Al contrario, se la matrice è malcondizionata, la soluzione potrebbe essere poco accurata. Si invita a controllare il condizionamento della matrice A con la routine matlab cond(A).

ESEMPIO D'USO:

A=rand(4);
b=rand(4,6);
x=sislin(A,b)
x =
-0.2703 0.0620 0.0987 0.7654 -0.2380 0.8251
0.4077 0.4003 0.1056 0.4601 0.5383 0.4422
0.1431 0.1797 0.0561 0.7122 -0.5151 0.7930
0.2749 0.4722 0.2308 -0.6636 1.1931 -0.5502

sislin.m

%sislin
% x=sislin(A,b)
% La funzione sislin calcola, a seconda del numero di
% parametri in input, l'inversa della matrice A (un solo
% parametro di input), oppure la risoluzione di m sistemi
% lineari, con A matrice dei coefficienti e b (nxm)
% matrice dei termini noti (in caso di due parametri in
% input)

ESEMPIO D'USO:

A=rand(4)
A =
0.2077 0.8443 0.2277 0.4302
0.3012 0.1948 0.4357 0.1848
0.4709 0.2259 0.3111 0.9049
0.2305

0.1707 0.9234 0.9797
b=rand(4,6)
b =
0.4389 0.5949 0.2217 0.4242 0.8010 0.4886
0.1111 0.2622 0.1174 0.5079 0.0292 0.5785
0.2581 0.6028 0.2967 0.0855 0.9289 0.2373
0.4087 0.7112 0.3188 0.2625 0.7303 0.4588

x=sislin(A,b)
x =
-0.2703 0.0620 0.0987 0.7654 -0.2380 0.8251
0.4077 0.4003 0.1056 0.4601 0.5383 0.4422
0.1431 0.1797 0.0561 0.7122 -0.5151 0.7930
0.2749 0.4722 0.2308 -0.6636 1.1931 -0.5502

function [x]=sislin(A,b)
if nargin==0
error('Numero di dati in input insufficiente')
end
[n m]=size(A);
%controllo per il calcolo dell'inversa o delle soluzioni di m sistemi
%lineari
if nargin==1
b=eye(n);
elseif nargin==2
if size(b)~=n
error('b deve avere lo stesso numero di righe di A')
end
end
%controllo che la matrice A non sia di caratteri
if (ischar(A)==1)
error('A deve essere una matrice di reali')
end
%Controllo errore per dimensioni nulle
if isempty(A)
error ('La matrice è vuota,nulla da calcolare')
end
%Controllo per matrici non

quadrate
if(n~=m)
    error('Attenzione matrice non quadrata');
endif

if isreal(A)==0
    error('A deve essere una matrice di reali')
end

%controllo che la matrice b non sia di caratteri
if (ischar(b)==1|| isreal(b)==0)
    error('b deve essere una matrice di reali')
end

zero=eps*norm(A,inf); %calcola lo zero per la verifica
piv=(1:n);%ciclo di fattorizzazione LU con pivoting parziale e scambio virtuale delle righe

for k=1:(n-1)
    r=min(find(abs(A(piv(1:n),k))==max(abs(A(piv(k:n),k))))); %cerca il pivot maggiore in modulo
    
    if abs(A(piv(r),k))>zero %controllo che il pivot considerato sia diverso da zero
        piv([r k])=piv([k r]); %scambio virtuale delle righe
        A(piv((k+1):n),k)=A(piv((k+1):n),k)/A(piv(k),k); %genera moltiplicatori, elisalva in place
        A(piv((k+1):n),(k+1):n)=A(piv((k+1):n),(k+1):n)-(A(piv((k+1):n),k)*A(piv(k),(k+1):n)); %sottrae riga a riga moltiplicata
    else %altrimenti avvisa l'utente che la matrice è singolare
        error ('Attenzione il sistema è singolare.Il')
    endif
endfor

programma verrà arrestato!')endendif abs(A(piv(n),n))<=zero %verifica singolarità matrice dei coefficienti per l'ultimo valoreerror('Attenzione sistema singolare su ultimo passo, Il programma verrà arrestato!')endt=size(b);%ciclo per calcolo della soluzione per ogni vettore colonna della matrice dei termini noti%forward substitutiony(1,:)=b(piv(1),:);for i=2:ny(i,:)=b(piv(i),:)-(A(piv(i),1:(i-1))*(y(1:(i-1),:)));end%back substitutionx(n,:)=y(n,:)/A(piv(n),n); %salvataggio dell'ultimo valore del vettore sol.for i=(n-1):-1:1x(i,:)=(y(i,:)-(A(piv(i),(i+1):n)*(x((i+1):n,:))))/A(piv(i),i);%calcolo delle soluzioniend

Test dell'algoritmo:

Caso funzionante1:

A=rand(4);b=rand(4,6);x=sislin(A,b)x =-0.2703 0.0620 0.0987 0.7654 -0.2380 0.82510.4077 0.4003 0.1056 0.4601 0.5383 0.44220.1431 0.1797 0.0561 0.7122 -0.5151 0.79300.2749 0.4722 0.2308 -0.6636 1.1931 -0.5502

Caso funzionante2:

>> A=rand(4);>> x=sislin(A)x =-0.6950 2.6552 1.5099

-1.5902

1.3591 -0.0776 -0.4835 -0.1355

-0.1598 1.2550 -1.1877 0.9303

0.0773 -1.7939 0.8484 0.5416

Casi d'errore:

>> sislin??? Error using ==> sislin at 32Numero di dati in input insufficiente

>> A=rand(4);>> b=rand(7,8);>> sislin(A,b)??? Error using ==> sislin at 42b deve avere lo stesso numero di righe di A

>> sislin('ciao')??? Error using ==> sislin at 48A deve essere una matrice di reali

>> sislin([])??? Error using ==> sislin at 52La matrice è vuota,nulla da calcolare

>> sislin(rand(4,7))??? Error using ==> sislin at 56Attenzione matrice non quadrata

>> sislin(rand(4),'ciao')??? Error using ==> sislin at 60b deve essere una matrice di reali

>>A =0.7184 0 0.0908 0.4401

0.9686 0 0.2665 0.5271

0.5313 0 0.1537 0.4574

0.3251 0 0.2810 0.8754

>> sislin(A)??? Error using ==> sislin at 74Attenzione il sistema è singolare. Il programma verrà arrestato!

>>A =0.5181 0.2407

0.6951 0.66780.9436 0.6761 0.0680 0.84440.6377 0.2891 0.2548 0.34450 0 0 0 0>> sislin(A)??? Error using ==> sislin at 78Attenzione sistema singolare su ultimo passo. Il programma verrà arrestato!>> a=rand(10);>> a=a+i;>> sislin(a)??? Error using ==> sislin at 61A deve essere una matrice di reali>> b=rand(10);>> b=b+i;>> sislin(rand(10),b);??? Error using ==> sislin at 65b deve essere una matrice di realiTest per l'accuratezza:Risoluzione m sistemi bencondizionati1:>> a=rand(7);x=rand(7,9);b=a*x;>> cond(a)ans =2.112908744886480e+001>> xc=sislin(a,b);>> norm(b-a*xc)/norm(a*xc)ans =8.913894804626601e-017>> err=norm(x-xc)/norm(x)err =1.392538378092079e-015Con una matrice bencondizionata, il residuo, che è nell'ordine dell'errore di round-off,dimostra che l'algoritmo trova una soluzione con la massima accuratezza. L'errorerelativo anch'esso dell'ordine

dell'errore di round off, invece, dimostra che la soluzione trovata è molto vicina a quella reale, con la perdita al massimo di un'unica cifra significativa.

Risoluzione m sistemi ben condizionati 2:

>> a=rand(10);
>> cond(a)
ans = 6.750871549975865e+001
>> xc=sislin(a);
>> ic=xc*a;
>> i=eye(10);
>> err=norm(i-ic)/norm(i)
err = 9.728661661915149e-015

L'errore relativo dimostra che anche per il calcolo dell'inversa l'algoritmo fornisce la soluzione con la massima accuratezza, con un errore al più sull'ultima cifra significativa.

Risoluzione m sistemi mal condizionati 1:

>> a=hilb(10);
x=rand(10,2);
b=a*x;
>> cond(a)
ans = 1.602467527403655e+013
>> xc=sislin(a,b);
>> norm(b-a*xc)/norm(a*xc)
ans = 1.357295860779872e-016
>> err=norm(x-xc)/norm(x)
err = 2.668192003144126e-004

Anche questa volta, il residuo dell'ordine dell'errore di round off dimostra che l'algoritmo risolve il problema in modo accurato.

Però l'errore relativo mostra che c'è stata una perdita di circa 12 cifre significative, questo perché la matrice è malcondizionata, quindi la soluzione trovata è lontana da quella reale.