Funzione 'sislin'

Funzione 'sislin'

SCOPO:

Calcolare la soluzione di m sistemi lineari di n equazioni in n incognite, nel caso di due parametri in input, o dell'inversa di una matrice, nel caso di un unico parametro di input.

SPECIFICHE:

x=sislin(A,b)

DESCRIZIONE ALGORITMO:

Per il calcolo della soluzione, la funzione fa uso di tre algoritmi.

·L' algoritmo di fattorizzazione LU con pivoting parziale e scambio virtuale delle righe, per calcolare le matrici triangolari inferiori e superiori L ed U, e la matrice P delle permutazioni (virtuale).

·L'algoritmo di forward subtitution per calcolare y tale che Ly=Pb

·L'algoritmo di back subtitution per calcolare la soluzione x tale che Ux=y.

LISTA PARAMETRI:

Input:

.A matrice quadrata dei coefficienti reali delle n equazioni, o matrice quadrata di cui calcolare l'inversa.

.b facoltativo matrice nxm dei termini noti, per la risoluzione di m sistemi lineari di n equazioni in n incognite. Se b non viene dato in input, la funzione calcola l'inversa della matrice A.

Output:

.x matrice inversa della matrice A, nel caso di un solo dato in input

.x matrice nxm, le cui colonne sono le m soluzioni degli m sistemi lineari Ax=b, nel caso di due dati di input

INDICATORI DI ERRORI:

'Numero di dati in input insufficiente'

Non ci sono parametri di input

.'b deve avere lo stesso numero di righe di A'-

Il vettore/matrice b deve avere tanti termini noti quante sono le equazioni del sistema lineare.

.'A deve essere una matrice di reali'

La matrice A inserita è una stringa

.'La matrice è vuota,nulla da calcolare'

A è una matrice vuota

.'Attenzione matrice non quadrata'

A non è una matrice quadrata nxn

'A deve essere una matrice di reali'

A deve essere una matrice di reali

.'b deve essere una matrice di reali'

La matrice b inserita è una stringa o una matrice di complessi

.'Attenzione il sistema è singolare. Il programma verrà arrestato!'

Il sistema da risolvere è un sistema singolare, quindi la soluzione non è unica

.'Attenzione sistema singolare su ultimo passo.  Il programma verrà arrestato!'

Il sistema da risolvere è un sistema singolare, quindi la soluzione non è unica

COMPLESSITA':

La complessità spaziale di tale algoritmo è S(n)=O(n²), in            quanto l'algoritmo è di tipo in place, cioè l'unico spazio di memoria che occupa è quello della matrice nxn, usando per le variabili interne lo spazio di memoria già allocato e non più utile. La complessità temporale invece è: T(n)=O(n³/3+mn²) con m numero dei sistemi lineari da risolvere.

ACCURATEZZA:

L'accuratezza dell'algoritmo è massima, nel caso in cui la            matrice A dei coefficienti sia una matrice bencondizionata. 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'

end

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 pivotmaggiore 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, e lisalva 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 programma verrà arrestato!'

end

end

if abs(A(piv(n),n))<=zero %verifica singolarità matrice dei coefficienti perl'ultimo valore

error('Attenzione sistema singolare su ultimo passo, Il programma verrà arrestato!'

end

t=size(b);

%ciclo per calcolo della soluzione per ogni vettore colonna

%della matrice dei termini noti

%forward substitution

y(1,:)=b(piv(1),:);

for i=2:n

y(i,:)=b(piv(i),:)-(A(piv(i),1:(i-1))*(y(1:(i-1),:)));

end

%back substitution

x(n,:)=y(n,:)/A(piv(n),n); %salvataggio dell'ultimo valore del vettore sol.

for i=(n-1):-1:1

x(i,:)=(y(i,:)-(A(piv(i),(i+1):n)*(x((i+1):n,:))))/A(piv(i),i);%calcolo dellesoluzioni

end

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.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

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 32

Numero di dati in input insufficiente

>> A=rand(4);

>> b=rand(7,8);

>> sislin(A,b)

??? Error using ==> sislin at 42

b deve avere lo stesso numero di righe di A

>> sislin('ciao')

??? Error using ==> sislin at 48

A deve essere una matrice di reali

>> sislin([])

??? Error using ==> sislin at 52

La matrice è vuota,nulla da calcolare

>> sislin(rand(4,7))

??? Error using ==> sislin at 56

Attenzione matrice non quadrata

>> sislin(rand(4),'ciao')

??? Error using ==> sislin at 60

b 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 74

Attenzione il sistema è singolare. Il programma verrà arrestato!

>>A =

0.5181 0.2407 0.6951 0.6678

0.9436 0.6761 0.0680 0.8444

0.6377 0.2891 0.2548 0.3445

0 0 0 0 0

>> sislin(A)

??? Error using ==> sislin at 78

Attenzione sistema singolare su ultimo passo. Il programma verrà arrestato!

>> a=rand(10);

>> a=a+i;

>> sislin(a)

??? Error using ==> sislin at 61

A deve essere una matrice di reali

>> b=rand(10);

>> b=b+i;

>> sislin(rand(10),b);

??? Error using ==> sislin at 65

b deve essere una matrice di reali

Test 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-

>> err=norm(x-xc)/norm(x)

err =

1.392538378092079e-

Con 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'errore relativo 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 bencondizionati2:

>> 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 malcondizionati1:

>> 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.

Risoluzione m sistemi malcondizionati2:

>> a=hilb(10);

>> cond(a);

>> ans

ans =

1.602467527403655e+013

>> xc=sislin(a);

>> ic=xc*a;

>> i=eye(10);

>> err=norm(i-ic)/norm(i)

err =

3.693850179255893e-003

Nel caso di una matrice malcondizionata, il calcolo dell'inversa subisce una perdita di 13 cifre significative.