Corso di Teoria ed applicazione delle macchine calcolatrici

Elaborato induviduale :

Fattorizzazione LU con pivoting parziale (FORTRAN)

Simulazione ecosistema a quattro componenti (MATLAB)

INDICE

Fattorizzazione LU con pivoting parziale (FORTRAN) 545b16f   3

MANUALE D'USO 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f   4

DESCRIZIONE DEL PROBLEMA 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f     4

MAIN PROGRAM LINSIST 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f     6

SUBROUTINE FATLU ( a , lda , n , punt ) 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f    8

SUBROUTINE TRDWN ( a , lda , n , punt , b , y ) 545b16f 545b16f 545b16f 545b16f 545b16f   10

SUBROUTINE TRUP ( a , lda , n , y, x ) 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f   12

FUNCTION PIV ( x , n ) 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f   14

TEST (contollo eseguito con MATLAB) 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f    15

1° test - Matrice 5x5 ben condizionata - 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f   15

2° test - Matrice 8x8 ben condizionata - 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f   16

3° Test - Matrice 8x8 mal condizionata - 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f     17

SIMULAZIONE ECOSISTEMA A QUATTRO COMPONENTI (MATLAB) 545b16f 18

TESTO 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f     19

M-File Ecosist.m 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f    20

Function Sisteco.m 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f 545b16f    22

Fattorizzazione LU con pivoting parziale

(FORTRAN)

MANUALE D'USO

DESCRIZIONE DEL PROBLEMA

Il software di seguito sviluppato si propone di risolvere sistemi di equazioni lineari del tipo :

Ax = b

In realtà il sistema di equazioni precedente viene semplificato e diviso in due sottosistemi. Ciò avviene fattorizzando , cioè triangolarizzando , la matrice dei coefficienti A utilizzando l'algoritmo di Gauss , quindi :

A = LU

dove L è la matrice dei matrice dei moltiplicatori ed U la matrice finale dell'algoritmo di Gauss. Inoltre poiché l'algoritmo di Gauss è sensibile alla precisione finita si utilizza il cosiddetto pivoting parziale . Esso consiste , al generico passo k , nella ricerca nella colonna del pivot del massimo , al fine di ottimizzarne la stabilità. Poiché dopo la ricerca dela massimo si effettua uno scambio di righe la fattorizzazione precedente si trasforma in :

PA = LU

dove PA è la matrice dei coefficienti permutata.

Nota la fattorizzazione precedente il sistema Ax = b si trasforma in :

LUx = b

I due sistemi equivalenti saranno quindi:

Ly = b

Ux = y

che risolti l'uno dopo l'altro , rispettivamente con una forward ed una back sobstituton , forniscono la soluzione x.

Il software è stato sviluppato in Fortran 77 ed è costituito da 5 File sorgente :

Linsist . FOR  Main program Linsist

Fatlu . FOR  Subroutine Fatlu per la fattorizzazione LU

Trdwn . FOR  Subroutine per la risoluzione di Ly=Pb

Trup . FOR  Subroutine per la risoluzione di Ux=y

Piv . FOR  Function per la ricerca del massimo pivot

i rispettivi . OBJ di compilazione , ed uno eseguibile :

Linsist . EXE

MAIN PROGRAM LINSIST

Scopo : risolve il sistema lineare Ax=b (con pivoting parziale, in place) ;

Descrizione dei parametri :

A : - input - array bidimensionale che contiene la matrice dei coefficienti del sistema ;

n : - input - ordine del sistema ;

b : - input - array monodimensionale che contiene il vettore dei termini noti del sistema ;

sol : - output - array monodimensionale che contiene la soluzione calcolata dalla procedura .

Subroutine chiamate : Fatlu ,Trdwn , Trup.

Program linsist

c Main program chiamante le subroutine fatlu,trdwn,trup

integer lda,punt(100),n

c dichiarazione dell'array dei puntatori e della dimensione n di A

real a(100,100),b(100),x(100),y(100),sol(100)

c dichiarazione degli array della matrice A ,del vettore dei termini noti b

c ( permutato y ), delle incognite x e delle soluzioni sol

lda=100

c leading dimension

print *,'Inserire la dimensione N della matrice quadrata A'

read *,n

print *,'Inserire per colonna gli elementi della matrice A'

read *,((a(i,j),i=1,n),j=1,n)

print *, 'Inserire gli N elementi del vettore dei termini noti b'

read *,(b(i),i=1,n)

c inserimento della matrice A e del vettore dei termini noti b

call fatlu(a,lda,n,punt)

c chiamata alla subroutine fatlu che da come output la matrica A fattorizzata e

c l'array dei puntatori in cui sono memorizzati gli scambi di riga

call trdwn(a,lda,n,punt,b,y)

c chiamata alla subroutine trdwn che risolve la parte del sistema Ly=Pb

call trup(a,lda,n,y,sol)

c chiamata alla subroutine trup che risolve Ux=y restituendo l'array delle soluzioni

print *,'Soluzione del sistema lineare impostato'

print *,(sol(i),i=1,n)

c stampa a video delle soluzioni del sistema di equazioni lineari

stop

end

SUBROUTINE FATLU ( a , lda , n , punt )

Scopo :calcolo di PA=LU (pivoting parziale , in place) ;

Descrizione dei parametri :

A : - input - array bidimensionale contenente la matrice dei coefficienti del sistema ;

n : - input - ordine del sistema ;

punt : - array monodimensionale contenente il vettore dei puntatori .

Function chiamate: Piv.

Complessità di tempo T(n) : O

subroutine fatlu(a,lda,n,punt)

integer n,punt(*),piv,jt,ir

c dichiarazione delle variabili locali e dell'array del puntatori

real a(lda,*)

do 10 j=1,n

punt(j)=j

10 continue

c inizializzazione del vettore puntatori

do 20 k=1,n-1

ir=piv(a(k,k),n-k+1)+k+1

c ricerca del max nella colonna con richiamo alla function PIV

if ( a(ir,k).eq.0) then

print *, 'La matrice immessa e SINGOLARE'

go to 100

endif

c struttura di controllo sulla singolarità della matrice A

jt=punt(k)

punt(k)=punt(ir)

punt(ir)=jt

c scambio degli indici nel vettore puntatori

do 30 j=1,n

t=a(k,j)

a(k,j)=a(ir,j)

a(ir,j)=t

30 continue

c scambio tra le righe per la massimizzazione del coefficiente pivot

do 40 i=k+1,n

a(i,k)=a(i,k)/a(k,k)

c calcolo del moltiplicatore

do 50 j=k+1,n

a(i,j)=a(i,j)-a(i,k)*a(k,j)

c aggiornamento coefficienti matrice attiva

50 continue

40 continue

20 continue

100 return

end

SUBROUTINE TRDWN ( a , lda , n , punt , b , y )

Scopo : risoluzione del sottosistema equivalente Ly = Pb ( Triangolare Inferiore )

Viene utilizzato l'algoritmo della forward substitution.

Descrizione dei parametri :

A : - input - array bidimensionale contenente la matrici dei coefficiente del sistema ;

n : - input - ordine del sistema ;

punt : - input - array monodimensionale contenente il vettore puntatori ;

b : - input - array monodimensionale che contiene il vettore di termini noti ;

y : - output - array monodimensionale contenente il vettore termini noti dopo la permutazione.

Subroutine chiamate : nessuna .

Complessità di tempo T(n) : O

subroutine trdwn(a,lda,n,punt,b,y)

integer n,lda,punt(*)

real a(lda,*) , b(*) , y(*)

c dichiarazione variabili I/O

y(1)=b(punt(1))/a(1,1)

do 10 i=2,n

y(i)=b(punt(i))

do 20 k=1,i-1

y(i)=y(i)-a(i,k)*y(k)

c scambi tra gli elementi del vettore x memorizzati ne vettore puntatori

c risoluzione del sistema Ly=Pb

20 continue

y(i)=y(i)/a(i,i)

10 continue

return

end

SUBROUTINE TRUP ( a , lda , n , y, x )

Scopo : risoluzione del sottosistema equivalente Ux = y  (Triangolare Superiore).

Viene utilizzato l'algoritmo della back substitution.

Descrizione dei parametri :

A : - input - array bidimensionale contenente la matrice dei coefficienti del sistema ;

n : - input - ordine del sistema ;

y : - input - array monodimensionale contenente il vettore dei termini noti permutato;

x : - output - array monodimensionale contente il vettore delle soluzioni del sistema ;

Subroutine chiamate : nessuna .

Complessità di tempo T(n) : O

subroutine trup(a,lda,n,y,x)

integer lda,n

real a(lda,*) , y(*) , x(*)

c dichiarazione variabili I/O

x(n)=y(n)/a(n,n)

do 10 i=n-1,1,-1

x(i)=y(i)

do 20 k=i+1,n

x(i)=x(i)-a(i,j)*x(k)

20 continue

x(i)=x(i)/a(i,i)

c procede alla back substitution calcolando le soluzioni del sistema Ux=y

c e restituendo nel programma principale l'array delle soluzioni

10 continue

return

end

FUNCTION PIV ( x , n )

Scopo : ricerca del massimo nella colonna dove si elimina durante la fattorizzazione .

La function utilizza un array monodimensionale ( x ) per memorizzare la

colonna dove ricercare il max , memorizzando in piv l'indice di riga che

verrà poi utilizzato nella subroutine chiamante per lo scambio di righe.

Descrizione dei parametri :

x : - input - array monodimensionale contenente la colonna dove ricercare il massimo ;

n : ordine del sistema ;

Subroutine chiamate : nessuna .

Complessità di tempo T(n) : O(n)

Integer function piv(x,n)

integer n,k,g

real x(*),max

g=1

c inizializzazione di g

max=x(1)

do 10 k=2,n

if (abs(x(k)).gt.abs(max)) then

c se al posto k il valore di x(k) è maggiore del max

c pone g=k

g=k

max=x(g)

c e pone il valore trovato uguale al massimo

endif

10 continue

piv=g

c memorizza il numero di riga del pivot massimo

return

end

TEST (controllo eseguito con MATLAB)

1° test - Matrice 5x5 ben condizionata -

Matrice dei coefficienti A

cond(a) = 17.1443

Det(a) = 0.1622

Vettore dei termini noti b

( 0.8278 0.1254 0.0159 0.6885 0.8682)^T

Vettore soluzioni x=A\b

Risultato MATLAB 545b16f 545b16f   Risultato LINSIST

-0.6531 545b16f 545b16f 545b16f -0.6516

0.3435 545b16f 545b16f 545b16f    0.3559

-0.0180 545b16f 545b16f 545b16f -0.1530

1.6046 545b16f 545b16f 545b16f    1.5751

-0.3632 545b16f 545b16f 545b16f -0.3593

OK !

2° test - Matrice 8x8 ben condizionata -

Matrice dei coefficienti A

cond(a) = 31.7362

Det(a) = 0.0290

Vettore dei termini noti b

( 0.0737 0.5007 0.3841 0.2771 0.9138 0.5297 0.4644 0.9410)^T

Vettore delle soluzioni x=a\b

Risultato MATLAB 545b16f 545b16f   Risultato LINSIST

1.2061 545b16f 545b16f 545b16f    1.2059

-0.7482 545b16f 545b16f 545b16f -0.7479

-1.3227 545b16f 545b16f 545b16f -1.3229

0.1052 545b16f 545b16f 545b16f    0.1049

0.3386 545b16f 545b16f 545b16f    0.3388

0.8584 545b16f 545b16f 545b16f -0.8586

-0.9028 545b16f 545b16f 545b16f -0.9030

1.1888 545b16f 545b16f 545b16f 1.1892

OK !
3° Test - Matrice 8x8 mal condizionata -

Matrice dei coefficienti A

cond(a) = 2.6106e+003 (mal condizionata cond > 10³/10⁴)

Det(a) = -8.0436e-005

Vettore dei termini noti b

( 0.4414 0.4314 0.0590 0.3832 0.6443 0.5180 0.2089 0.7378)^T

Vettore delle soluzioni x=a\b

Risultato MATLAB 545b16f 545b16f   Risultato LINSIST

35.4356 545b16f 545b16f 545b16f   36.7736

11.9939 545b16f 545b16f 545b16f   12.2816

191.2157 545b16f 545b16f     199.5888

-52.9895 545b16f 545b16f     -55.1659

-54.0899 545b16f 545b16f     -56.3444

196.0114 545b16f 545b16f     204.5518

-283.1839 545b16f 545b16f    -295.6083

-116.8741 545b16f 545b16f    -121.9154

OK !

TESTO

Si consideri un sistema ecologico costituito da quattro sottosistemi (comparti) interagenti : 1-produzione primaria ; 2-erbivori ; 3-carnivori ; 4-decompositori.

Sia x_i(t) il contenuto energetico del sottosistema  i al tempo t.

I flussi di energia sono modellizzati dal sistema

Risolvere il sistema nell'intervallo [0,800] per i seguenti valori di saturazione della predazione :

k₁=2000 ; k₂=60 e con le condizioni iniziali :

x₁(0)=2800, x₂(0)=190, x₃(0)=35, x₄(0)=66.

Riportare in un grafico le quattro funzioni x_i(t) al variare di t.

Riportare le traiettorie , nello spazio degli stati , del contenuto energetico degli erbivori e della produzione primaria , dei carnivori e della produzione primaria , dei decompositori e della produzione primaria

M-File Ecosist.m

clc

clear

echo on

% Corso di teoria ed applicazione delle macchine calcolatrici

Prof. G.Giunta - A.A. 1995/96

Elaborato individuale

% 545b16f a cura di

% Carlo Pinto (SB/4) e Carlo De Luca (SB/36)

% Sistema ecologico costituito da 4 sottosistemi (comparti) interagenti:

% 1-Produzione primaria 2-Erbivori 3-Carnivori 4-Decompositori

% Premere un tasto per il grafico dei flussi di energia

echo off

pause

global k1 k2

k1=2000;

k2=60;

% Valori di saturazione della predazione

[t,u]=ode45('sisteco',0,800,[2800 190 35 66]);

%Funzione di risoluzione delle 4 eq. differenziali,

%contenute nella function 'sisteco',

%restituisce un vettore t del tempo ed una matrice u delle soluzioni

% N.B.: la prod.primaria è in scala 1:100,gli erbivori 1:10

le condizioni iniziali sono rispettivamente 2800 e 190

figure(1)

plot(t,u(:,1)/100,'+y',t,u(:,2)/10,'-.',t,u(:,3),'o',t,u(:,4),'b')

title('Sistema ecologico a 4 componenti:prod.primaria,erbivori,carnivori,decompositori')

xlabel('Tempo')

ylabel('Contenuto energetico dei sottosistemi in [kcal/(m²* a)]')

text(200,70,'Produzione primaria in scala 1:100')

text(200,12,'Erbivori in scala 1:10')

text(200,25,'Carnivori')

text(200,0,'Decompositori')

echo on

% Premere un tasto per avere i grafici nello spazio degli stati

echo off

pause

figure(2)

plot(u(:,1),u(:,2),'+',u(:,1),u(:,3),'-.',u(:,1),u(:,4),'o')

title('Spazio degli stati')

xlabel('Produzione primaria')

ylabel('Erbivori(1)-Carnivori(2)-Decompositori(3)')

text(1500,500,'(1)')

text(5300,50,'(2)')

text(3500,0,'(3)')

Function Sisteco.m

function uprimo=sisteco(t,u)

global k1 k2

% richiamo alle costanti fissate nell'm-file chiamante (ecosist.m)

uprimo(1)=400+175*sin(2*pi*(t-11)/52)-1.5 * u(1)*u(2)/(k1+u(1))-0.1*u(1);

% d(x₁)/dt = flusso di energia modellizzato per la Produzione primaria

uprimo(2)=0.6* u(1)*u(2)/(k1+u(1))-5.7 *u(3)*u(2)/(k2+u(2))-0.21*u(2);

% d(x₂)/dt = flusso di energia per gli Erbivori

uprimo(3)=0.17 * (u(3)*u(2)/(k2+u(2))-0.12*u(3);

% d(x₃)/dt = flusso di energia per i Carnivori

uprimo(4)=0.03*u(1)+0.06*u(2)+0.02*u(3)-3*u(4);

% d(x₄)/dt = flusso di energia per i Decompositori