Introduzione

Qualsiasi fenomeno elettromagnetico è governato dalle equazioni di Maxwell, alle quali bisogna aggiungere le condizioni al contorno per individuarlo univocamente. La soluzione analitica del problema, descritto dal modello matematico precedente, può essere ottenuta solo in casi semplici caratterizzati da elevate simmetrie.

E' ovvio che l'approccio ai problemi di interesse pratico non può che essere di tipo numerico, il quale non presenta alcuna limitazione di carattere teorico sulla forma dei contorni della regione da studiare. Utilizzando un metodo numerico, il sistema di equazioni differenziali è sostituito da un sistema di equazioni algebriche risolvibile con un metodo diretto o iterativo. Qualora si adottasse un metodo di risoluzione iterativo, occorre porre particolare attenzione alla stabilità e convergenza dello schema numerico utilizzato, nonché al tempo di calcolo. Nel seguito si farà riferimento al metodo delle differenze finite.

Il metodo delle differenze finite: caso bidimensionale

In questo articolo applicheremo il metodo delle differenze finite all'equazione di Laplace:

Se la regione da studiare è illimitata, è necessario fissare un contorno (con condizioni assegnate) a grande distanza rispetto alla zona di reale interesse. Alla regione si associa un reticolo, preferibilmente regolare, individuando così dei nodi in cui si valuterà la funzione incognita (la funzione potenziale). Consideriamo la seguente stella di reticolo a bracci equali con passo k:

Stella.jpg

Sviluppiamo la funzione potenziale in serie di Taylor nell'intorno del nodo 0 secondo l'asse x:

Sommando membro a membro le due espressioni precedenti, trascurando i termini con potenze di k superiori alla seconda, si ottiene:

dalla quale si ricava:

Sviluppando secondo l'asse y con procedimento identico al precedente, si ricava:

Sommando le due espressioni precedenti (tenendo conto dell'equazione di Laplace), si ottiene:

Questa equazione conduce alla formulazione del seguente sistema matriciale:

che presenta un numero di equazioni algebriche pari al numero di nodi interni.

Il problema viene affrontanto con la tecnica dei residui, assegnando arbitrariamente i valori di potenziali nei punti interni del reticolo. Ovviamente non si arriverà mai alla soluzione esatta del problema, ma, dopo un certo numero di iterazioni, si giungerà ad un insieme di valori di potenziali che presentano un errore, rispetto al valore vero, inferiore ad una soglia prestabilita (x% del valore medio dei valori di potenziali nei nodi di reticolo). E' chiaro che lo schema deve essere stabile e convergente affinchè si possa giungere alla soluzione. Il procedimento iterativo adottato è lo schema di Gauss-Seidel che presenta elevata rapidità di convergenza e risparmio di memoria.

Codice Fortran 90/95

!***************************************************************
!*													      *
!* Programmatore: EdmondDantes       							      *
!* 										        		      *
!* Copyright: Nessuno									              *
!*									                		      *
!* Data: 1 ottobre 2009									      *
!*										                	      *
!***************************************************************
!
!***************************************************************
!* Descrizione:											      *		             
!*                                                                  	                                      *                        
!* Questo codice risolve l'equazione di Laplace nel dominio bidimensionale  	      *
!*                                                      	                                              *
!* risolvendo il sistema matriciale ottenuto col metodo delle differenze   	      *
!* 								                                              *
!* finite utilizzando il metodo iterativo di Gauss-Seidel.     				      *
!*													      *
!***************************************************************

PROGRAM Laplace

IMPLICIT NONE

REAL*8    :: alfa, u(100,100), R(100,100), media(100), &
                  max, Rmax 	
INTEGER*4 :: M, N, i, j, numero_iterazioni

OPEN(1,FILE='Potenziale.txt',STATUS='REPLACE',ACTION='WRITE')

!=================================================

! Dimensioni del reticolo associato

WRITE(*,*) 'Inserisci il numero di elementi'
WRITE(*,*) 'secondo y'
READ(*,*) N
WRITE(*,*) 'secondo x'
READ(*,*) M

WRITE(*,1000) M*N
1000 FORMAT(1X,'Il reticolo e'' costituito da',2X,I5,2X,'elementi')


!==================================================

! Coefficiente di rilassamento
! 0<alfa<1 sottorilassamento
! 1<alfa<2 sovrarilassamento

alfa=1.8

!==================================================

! Creo la matrice dei potenziali

DO 1 i=1,M
	DO 2 j=1,N

		u(i,j)=0

	2 CONTINUE
1 CONTINUE

!==================================================

! Condizioni al contorno

DO i=1,N

	u(1,i)=1
	u(M,i)=0

END DO

DO i=1,M

	u(i,1)=0
	u(i,M)=0

END DO

!=================================================


! Creo la matrice dei residui

DO 3 i=1,M
	DO 4 j=1,N

		R(i,j)=0

	4 CONTINUE
3 CONTINUE

!=================================================

numero_iterazioni=0 ! Costante che mi indicherà 
                    ! il numero di iterazioni eseguite
                    ! dal codice

!=================================================


DO 

	DO 5  j=1,N
	
		media(j)=0	

		DO 6  I=1,M
		
			media(j)=media(j)+u(i,j)
					
		6 CONTINUE
	5 CONTINUE


	DO 7 i=1,N

		media(i)=media(i)/M

	7 CONTINUE

 ! Con media(i) indico la media dei potenziali

	max=media(1)

	DO i=2,N

		IF(media(i)>max) max=media(i)
	
	END DO

	Rmax=0	! Definisco il residuo massimo
	DO 8 i=2,M-1
		DO 9 j=2,N-1

			R(i,j)=u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)-4*u(i,j)

			IF(R(i,j) > Rmax) Rmax=R(i,j)

			u(i,j)=u(i,j)+alfa*(R(i,j)/4)

		9 CONTINUE
	8 CONTINUE

	IF(Rmax < (0.001*max)) EXIT

numero_iterazioni=numero_iterazioni+1

END DO  


CALL stampa(M,N,u)
	 
WRITE(*,*) 'Numero di iterazioni',numero_iterazioni

CLOSE(1)

STOP
END PROGRAM Laplace

!**************** FINE PROGRAMMA PRINCIPALE ***********************!


!*************** Sezione Subroutine *********************************!

SUBROUTINE stampa(M,N,MAT)

REAL*8 MAT(100,100)
INTENT(IN) M,N
INTENT(OUT) MAT

DO 100 i=1,M
WRITE(1,1100) (MAT(i,j),j=1,N)
100 CONTINUE
1100 FORMAT(2X,100E20.10)

RETURN
END SUBROUTINE

Istruzioni

Il codice Fortran 90/95 deve essere salvato in una cartella al cui interno si creerà automaticamente il file Potenziale.txt, contenente la matrice MxN con i valori dei potenziali di tutti i nodi del reticolo associato alla regione piana studiata. I dati del file possono essere inviati in un software qualsiasi in grado di creare una superficie che evidenzia la distribuzione del potenziale.

Esempio

Fornendo al codice i seguenti dati:

• M=30 (a video)
• N=30 (a video)
• alfa=1.8 (modificando il valore presente nel codice)

si ottiene un file Potenziale.txt con 30x30 elementi.

La distribuzione del potenziale è la seguente:

Distribuzione del potenziale

Conclusioni

Spero di contribuire, con questo articolo, alla diffusione dell'uso "intelligente" del computer, auspicato da g.schgor, nell'articolo Il cane e la lepre...in moto qualsiasi