1 Potenziali
In questo tutorial scriveremo un codice per la risoluzione di equazioni differenziali ellitiche come quella nota in fisica come eq. di Laplace-Poisson per un potenziale generato da una carica :
I metodi iterativi costituiscono un modo efficiente per calcolare la soluzione di queste eqs. Come visto nelle dispense, si tratta di applicare il metodo FTCS alla eq. di Poisson eguagliata ad un termine temporale:
come notiamo ci permette di calcolare il potenziale all'iterazione successiva in un certo punto, facendo la media del potenziale nei punti attorno all'iterazione precedente. Questo ci dice anche che per iniziare l'iterazione abbiamo bisogno di sapere le condizioni al contorno della griglia e il potenziale iniziale su tutta la griglia.
1 2 3 4 | for i in range(1,N-1):
for j in range(1,N-1):
#Jacobi
U1[i,j] = 0.25*(U[i-1,j] + U[i+1,j] + U[i,j+1] + U[i,j-1])
|
Questa rappresenta la prima iterazione. Poichè noi cerchiamo una soluzione al limite, per la quale la variazione del potenziale sia nulla, dovremo implementare un modo per definire la convergenza. Un metodo rapido è quello di calcolare la media del poteziale su tutta la griglia ad ogni iterazione e confrontarla con il valore medio trovato alla iterazione precedente. Quando questa differenza è minore di un certo valore potremo considerare la soluzione trovata quella finale. Il codice seguente esegue ciò che vogliamo:
1 2 3 4 5 6 7 8 9 10 11 | conv=1
mean_U0=mean(U)
niter=0
while conv > 1E-04:
#ciclo for i
#ciclo for j
conv = abs(mean(U1)-mean_U0)
mean_U0=mean(U1)
U=U1.copy()
niter+=1
|
- le iterazioni vengono fatte mediante l'uso di un ciclo while, il quale esegue le righe al suo interno (quelle indentate) sintanto che la condizione conv>1E-04 è verificata;
- è necessario quindi aggiornare il valore di conv, quindi la differenza tra le medie a iterazioni successive;
- la funzione mean(arr) calcola la media su tutti i valori dell'array in input;
- conserviamo il valore attuale della media in mean_U0;
- conserviamo il potenziale attuale in U copiando U1 con la funzione copy();
- poichè non sappiamo a priori quante iterazioni saranno eseguite le contiamo incrementando la variabile niter;
- nelle prime righe impostiamo i valori iniziali per la media, il parametro di convergenza e il numero di iterazioni, mean_U0, conv, niter rispetivamente.
Non ci crederete ma questo è quanto basta. Anzi, possiamo anche semplificare ulteriormente, se utilizziamo il metodo di Gauss-Seidel:
che sfrutta il fatto che nei punti i-1 e j-1 il potenziale è gia stato calcolato. Questo oltre a velocizzare la convergenza, semplifica anche il codice perche possiamo usare la stessa matrice leggendo le nove componenti e scrivendo le nuove. Vediamo come il codice precedente diventa:
1 2 3 4 5 6 7 8 9 10 11 12 | conv=1
mean_U0=mean(U)
niter=0
while conv > 1E-04:
for i in range(1,N-1):
for j in range(1,N-1):
#Seidel
U[i,j] = 0.25*(U[i-1,j] + U[i+1,j] + U[i,j+1] + U[i,j-1])
conv = abs(mean(U)-mean_U0)
mean_U0=mean(U)
niter+=1
|
Resta quasi tutto invariato, tranne l'uso della sola matrice U nei cicli i,j, e non abbiamo quindi bisogno di salvare il potenziale ad ogni iterazione.
A questo punto implementiamo anche il metodo SOR, semplicemente traducendo in codice la seguente iterazione:
che andrà a sostituire la riga n del codice precedente. Per questo metodo dovremo definire il peso W nel range , solitamente vicino a 2.
1 2 | #SOR
U[i,j] = W*0.25*(U[i-1,j] + U[i+1,j] + U[i,j+1] + U[i,j-1])+ (1-W)*U[i,j]
|
Ora che abbiamo implentato i tre metodi, fissiamo le dimensioni del sistema, le condizioni al contorno e un valore iniziale del potenziale per tutta la griglia e confrontiamo il numero di iterazioni necessarie per la convergenza dei tre metodi:
1 2 3 4 5 6 7 8 | N=100
U=zeros((N,N))
U[:,0]=10.
#Jacobi
#U1=U.copy()
#SOR
#W.19
|
Usiamo una griglia 100x100 e supponiamo di avere un potenziale fisso pari a 10 lungo un lato e zero lungo gli altri tre lati. Stiamo anche partendo da un potenziale nullo su tutto il resto della griglia. Eseguendo i tre metodi separatamente, con queste condizioni iniziali i tre metodi eseguono il seguente numero di iterazioni:
1 2 3 4 5 | #Jacobi
#Seidel
#SOR
|
che possono essere facilmente confrontati con il grafico in figura 27. Se volessimo replicare interamente quel grafico, basterebbe racchiudere il nostro codice dentro un ciclo for come questo:
1 2 3 4 5 6 | niter_list=[]
for N in range(20,200,20):
#codice precedente per uno dei 3 metodi
niter_list.append(niter)
plot(range(20,200,20),niter_list)
|
Salviamo il numero di iterazioni per ogni valore di N mediante la lista niter_list, su cui aggiungiamo i nuovi valori mediante la funzione append(); il plot è eseguito al solito modo.
Sfruttando un ciclo simile ma fatto sul peso W e mantenendo la griglia fissa su 100, possiamo riprodurre il secondo grafico in fig 27, dove si mostra l'andamento del numero di terazioni con il peso:
1 2 3 4 5 6 7 | niter_list=[]
for W in arange(1,2,0.02):
N=100
#codice precedente per SOR
niter_list.append(niter)
plot(arange(1,2,0.02),niter_list)
|