Cos'è ElectroYou | Login Iscriviti

ElectroYou - la comunità dei professionisti del mondo elettrico

2
voti

Un metodo per la risoluzione dei sistemi lineari

Indice

Introduzione

Con questo nuovo articolo vorrei introdurre uno dei più semplici metodi numerici per risolvere un sistema di n equazioni lineari.

Sistemi lineari

Si consideri il seguente sistema lineare:

[A][x] = [b]

nel quale:

[A] è la matrice quadrata dei coefficienti;

[x] il vettore colonna delle incognite;

[b] il vettore colonna dei termini noti.

La risoluzione formale del sistema, come è noto, è:

[x] = [A] − 1[b]

Esistono diverse formulazioni per risolvere il problema, quella magica è la regola di Cramer:


x_{i}=\frac{\left | A_{i} \right |}{\left | A \right |}

con i=1,2...n.

La regola di Cramer è numericamente inutilizzabile poiché richiede la definizione di n+1 determinanti di ordine n, e ciascun determinante richiede n! operazioni per essere calcolato. Il costo computazionale del metodo è esorbitante: (n+1)! operazioni, anche se con un piccolo accorgimento è possibile abbassarlo. Bisogna abbandonare questa strada se desideriamo risolvere il sistema lineare di ordine n numericamente. I metodi numerici utilizzati per la risoluzione dei sistemi lineari sono:

  • Metodi diretti: si giunge alla soluzione con un numero finito di passi (abbiamo errori di arrotondamento).
  • Metodi iterativi: si giunge alla soluzione con un numero infinito di passi. E’ chiaro che bisogna accettare un errore di troncamento.

I metodi iterativi, generalmente, sono utilizzati quando la matrice [A] è sparsa. Le matrici sparse si incontrano quando si cerca di discretizzare un sistema di equazioni alle derivate parziali (per esempio un fenomeno retto dalla equazione di Laplace). In questo articolo introdurrò il metodo diretto noto come Metodo di eliminazione di Gauss con pivot parziale. Consideriamo la matrice aumentata

[A | b]

Affinché l'elemento pivot sia un elemento non nullo, si utilizza la tecnica del pivoting parziale o totale. Il codice Fortran Metodo_di_Gauss_con_pivot_parziale utilizza, come appare chiaro dal nome, la tecnica del pivoting parziale, pertanto si effettueranno solo scambi di righe e non di colonna.

Una piccola precisazione

Con la tecnica del pivoting parziale si cerca l'elemento più grande in modulo in tutta la colonna in esame. Se vogliamo parlare in termini più marziani, la tecnica del pivoting parziale calcola la norma infinito del vettore (colonna). Con la tecnica del pivoting totale si cerca l'elemento più grande in tutta la sottomatrice considerata; ciò comporta scambi di righe e colonna. Il pivoting totale è migliore rispetto al pivoting parziale, ma computazionalmente più costoso.

Il codice Fortran 90/95

!***************************************************************
!*													      *
!* Programmatore: EdmondDantes						                      *
!* 													      *
!* Copyright: Nessuno								                      *
!*													                                                                                              *
!* Data: 4 ottobre 2009	                  		              *
!*													      *
!***************************************************************
!													     !*
!***************************************************************
!* Descrizione:										              *
!*                                                                                                           *
!* Questo codice Fortran risolve un sistema lineare                                         *
!                                                                                                             *
!* di n incognite, utilizzando il metodo diretto di                                             *
!                                                	                                                      *
!* Gauss e pivot parziale.							                      *
!* 													      *
!***************************************************************
													     !*
!==================Inizio codice=======================*
													       !*
PROGRAM Metodo_di_Gauss_con_pivot_parziale			                       !*
													          !*
IMPLICIT NONE									                          !*
													          !*
!===================================================*
!													          !*
!-------------------Definisco le variabili---------------------                                                 *
													          !*
INTEGER*4   :: i, j, k, n, riga_pivot				                                          !*
REAL*8      :: a(100,100), pivot, c, v, det			                                 !*
													          !*
!===================================================!*
! Dizionario:										                             !*
!													          !*
! a(i,j) ==> elemento della matrice [A] dei coeff.	                                         !*
! det ==> variab. ove si immagazzina il determinante.                                           !*
! riga_pivot ==> variab. che tiene conto del numero	                                     !*
!                di riga dell'elemento pivot.		                                                         !*
! pivot ==> è l'elemento pivot.						                                !*
! n ==> numero di equazioni.						                               !*
!=====================================================!*
!													                  !*
!====================================================!*
OPEN(200,FILE='Sistema.txt',STATUS='OLD',ACTION='READ')                                  !*
OPEN(300,FILE='Soluzione.txt',STATUS='REPLACE',ACTION='WRITE')                         !*
WRITE(*,*) 'Introduci il numero n di equazione'		                                         !*
READ(*,*) n                                                                                                        !*
WRITE(*,*) 'Lettura matrice da file'                                                                         !*
READ(200,*) ((a(i,j), j=1,n+1), i=1,n)                                                                     !*
WRITE(*,*) 'Leggi il risultato nel file Soluzione'                                                           !*
!===================================================!***
													                   !*
det=1.D00 ! Inizializzo il determinante			                                                      !*
													                  !*
DO 100 k=1,n-1										                         !*
													                   !*
	pivot=DABS(a(k,k))							                         	!*
	riga_pivot=k									                       !*
													              !*
	DO 101 i=k+1,n									                     !*
													                  !*
		v= DABS(a(i,k))								                          !*
		IF(v > pivot) THEN							                          !*
													                !*
			pivot=DABS(a(i,k))						                           !*
			riga_pivot=i							                        !*
													                    !*
		END IF										                 !*
													                       !*
101 CONTINUE										                           !*
													                    !*
	IF(riga_pivot > k) THEN							                    !*
													               !*
		DO 102 j=1,n+1								                       !*
													                !*
			c=a(k,j)							                         	!*
			a(k,j)=a(riga_pivot,j)					                               !*
			a(riga_pivot,j)=c						                         !*
													!*
102		CONTINUE									!*
													!*
		det=-det									!*
													!*
	END IF											!*
													!*
	det=det*a(k,k)									!*
													!*
	DO 100 i=k+1,n									!*
													!*
		a(i,k)=a(i,k)/a(k,k)						!*
													!*
		DO 100 j=k+1,n+1							!*
													!*
			a(i,j)=a(i,j)-a(i,k)*a(k,j)				!*
													!*
100 CONTINUE						                !*
													!*
	det=det*a(n,n)									!*
	a(n,n+1)=a(n,n+1)/a(n,n)						!*
													!*
	DO 103 i=n-1,1,-1								!*
		DO 104 j=i+1,n								!*
													!*
			a(i,n+1)=a(i,n+1)-a(i,j)*a(j,n+1)		!*
													!*
104     CONTINUE									!*
													!*
		a(i,n+1)=a(i,n+1)/a(i,i)					!*
													!*
103 CONTINUE 										!*
                                                    !*
DO i=1,n                                            !*
  													!*
	WRITE(300,*) i, a(i,n+1)                        !*
                                                    !*
END DO                                              !*
						                            !*
CLOSE(200)                                          !*
CLOSE(300)                                          !*
  										            !*
STOP										           !*
END PROGRAM Metodo_di_Gauss_con_pivot_parziale		   !*
											   !*
!*****************************************************

Circuito in corrente continua

Consideriamo il circuito in corrente continua (a regime) proposto da admin come esercizio numero 9 nell'articolo Reti in c.c. - esercizi da svolgere.

circuito.jpg

circuito.jpg

I dati del problema sono:

  • R_1= 5 \, \Omega
  • R_2= 20 \, \Omega
  • R_3= 40 \, \Omega
  • R_4= 40 \, \Omega
  • E_1= 100 \, V
  • E_2= 200 \, V

L'esercizio deve essere risolto utilizzando la sovrapposizione degli effetti, ma admin ci perdonerà sicuramente se scriveremo il sistema di equazioni:

\begin{array}{*{20}c}
   {I_1 } & { + I_2 } & { + I_3 } &  =  & 0 & {}  \\
   {R_1 \cdot I_1 } & {} & { - (R_3+R_4) \cdot I_3 } &  =  & {E_1} & {}  \\
   {} & { - R_2 \cdot I_2 } & { + (R_3+R_4) \cdot I_3 } &  =  & {E_2} & {}  \\
\end{array}

Risolvendo manualmente il sistema di equazioni si ottiene:

  • I_1= 12,4 \, A
  • I_2= -11,9 \, A
  • I_3= -0,48 \, A

Soluzione numerica con Fortran

Seguendo le istruzioni riportate nell'articolo Using Fortran 90/95 (è possibile utilizzare qualsiasi compilatore; quello riportato è solo un esempio), l'uso del codice è molto semplice. Salviamo un file Sistema.txt nella stessa cartella che contiene il nostro codice fortran.f90 (possiamo nominarlo gauss.f90). Sostituiamo i valori numerici al sistema di equazioni precedente:


\begin{array}{*{20}c}
   {I_1 } & { + I_2 } & { + I_3 } &  =  & 0 & {}  \\
   {5I_1 } & {} & { - 80I_3 } &  =  & {100} & {}  \\
   {} & { - 20I_2 } & { + 80I_3 } &  =  & {200} & {}  \\
\end{array}


A questo punto, nel file Sistema.txt riportiamo i seguenti dati:


\begin{array}{*{20}c}
   1 & 1 & 1 & 0  \\
   {5} & 0 & { - 80} & {100}  \\
   0 & { - 20} & {80} & {200}  \\
\end{array}


ovvero la matrice aumentata già menzionata nell'articolo (matrice dei coefficienti più la colonna dei termini noti). Apriamo e compiliamo il file Fortran.

Ci verrà richiesto di inserire il numero di equazioni n (in questo caso 3 ==> Invio).

Dopo l'esecuzione del codice, apriamo il file Soluzione.txt (creato dal codice)

            1          12.3809523810
            2         -11.9047619048
            3        -0.476190476190

che conterrà la soluzione del problema.

Ovviamente i segni meno indicano che il verso delle correnti da noi scelto è errato.

Conclusioni

Il mio codice non può sostituire lo script di RenzoDF dell'articolo Simulazione reti in Scilab, però cerca di evidenziare un fatto:

La semplice divisione a sinistra oppure un comando di inversione del tipo inv(b) nasconde una serie di "conti" illustrati nel codice precedente, anzi ancora più precisi e quindi più laboriosi.

Il codice è in grado di lavorare solo con numeri reali, quindi possiamo inserire sistemi che derivano da un circuito in corrente continua o sistemi generici con coefficienti reali e soluzioni nel campo dei numeri reali.

1

Commenti e note

Inserisci un commento

di viola,

lo trovo davvero molto interessante e sicuramente utile. Complimenti!!!

Rispondi

Inserisci un commento

Per inserire commenti è necessario iscriversi ad ElectroYou. Se sei già iscritto, effettua il login.