T.1 Oscillatore
In questo tutorial scriveremo un codice per risolvere l'equazione differenziale dell'oscillatore armonico generica, che implementa i tre metodi presentati nella dispensa (pagine rif): Eulero, Eulero-Cromer, Verlet.
dove con indichiamo il passo temporale e con indichiamo l'accelerazione, o la forza diviso la massa, a cui è sottoposto il nostro sistema: nel nostro caso avremo:
La traduzione in codice di questi schemi è immediata. Vediamola di seguito e poi ci costruiamo il programma intorno:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | #Accelerazione
F0=-g/l*sin(x0)
#Eulero
v1=v0+tau*F0
x1=x0+tau*v0
#Eulero-Cromer
v1=v0+tau*F0
x1=x0+tau*v1
#velocity-Verlet
x1=x0+v0*tau+0.5*F0*tau**2
F1=-g/l*sin(x1)
v1=v0+0.5*tau*(F0+F1)
|
Come notiamo facciamo uso di due variabili una denominata con 0 l'altra con 1 che si riferiscono al tempo e rispettivamente. Se ora si definiscono le costanti necessarie (, , ) e le condizioni iniziali (,), la serie di righe scritte sopra ci permettono di calcolare la posizione e la velocità dopo un tempo , con il metodo che desideriamo. Vediamo di seguito tali definizioni:
1 2 3 4 5 | g=1
l=1
x0=deg2rad(5)
v0=0.0
tau=0.01
|
Facciamo notare quà l'uso della funzione deg2rad(alpha) che permette di convertire l'angolo alpha espresso in gradi sessagesimali in radianti. Questo è necessario perchè la funzione sin(alpha) che troviamo per calcolare la forza prende in input un angolo alpha in radianti. Esiste anche la funzione inversa rad2deg().
Ovviamente vorremo prolungare la nostra simulazione, calcolando la posizione e la velocità del nostro oscillatore, per un periodo di tempo piu lungo di in solo passo , per esempio per , con N a piacere. E facile capire che dovremo ripetere lo schema di integrazione scelto per un numero di volte, mediante un ciclo for, come vediamo di seguito:
1 2 3 4 5 6 7 8 | N=3000
for t_step in range(N):
#Eulero-Cromer
F0=-g/l*sin(x0)
v1=v0+tau*F0
x1=x0+tau*v1
x0=x1
v0=v1
|
Nelle ultime due righe ci ricordiamo di aggiornare la variabile della posizione e della velocità denominata con 0 con quelle appena calcolate, denominate con 1.
A questo punto il nostro programma calcola la posizione e la velocità dell'oscillatore per gli N istanti di durata , ma è ancora poco utile, perchè dopo l'uscità dal ciclo for avremo a disposizione solo l'ultimi valori di posizione e velocità. Sarebbe utile quindi registrare l'intera dinamica per poi studiarne l'andamento mediante un grafico, per esempio. Facciamo in questo modo:
- Nelle righe iniziali del programma, quindi prima del ciclo for, creiamo due arrays:
- registriamo nel primo elemento degli array le condizioni iniziali:
- all'interno del ciclo for registriamo ogni nuovo valore di x e v all'interno dei rispettivi arrays:
Ora una volta eseguite le righe di codice precedenti potremo fare il grafico della posizione (velocità) in funzione del tempo, creando prima un array con tutti gli istanti temporali usati nella integrazione, mediante la funzione linspace():
1 2 | time_arr=linspace(0,N*tau,N+1)
plot(time_arr,x_arr)
|
Se ora vogliamo riprodurre la fig 13 della dispensa, dovremo lasciare aperta la finestra con il grafico appena creato, cambiare l'angolo da 5 a 90 gradi ed eseguire nuovamente il programma. Cambiando ancora da 90 a 179, otteniamo:
Provate a fare il grafico mostrato in figura 14, dove si mostra la velocita' in funzione dell'angolo ( la nostra posizione ).
Passiamo ora al calcolo dell'energia cinetica e potenziale dell'oscillatore. L'implementazione e' del tutto analoga a quanto abbiamo fatto per salvare la dinamica:
- aggiungiamo all'inizio del codice due array:
- dentro il ciclo for calcoliamo i valori delle due energie usando la posizione e la velocita' ad ogni iterazione:
Ora quindi, ogni volta che eseguiamo il nostro programma, avremo a disposizione l'energia potenziale e cinetica per tutti gli istanti considerati e potremo facilmente farne il grafico:
1 2 3 | plot(time_arr,Ek_arr)
plot(time_arr,Ep_arr)
plot(time_arr,Ep_arr + Ek_arr)
|
In questo modo si ottiene uno dei due grafici di figura 16, a seconda dell'angolo scelto. Notiamo come nell'ultima riga possiamo fare il grafico dell'energia totale sommando i vettori direttamente dentro la funzione plot():
Il programma sin qui costruito permette di fare tutti i confronti mostrati nella dispensa, tra metodi differenti considerando lo stesso angolo di partenza (quindi oscillatore armonico o anarmomico) oppure tra oscillatori diversi usando lo stesso metodo.
Per evitare di dover cambiare il codice ogni volta che vogliamo provare un metodo o delle condizioni iniziali differenti, mostriamo di seguito come inserire dei parametri dall'esterno, subito dopo che il programma viene eseguito:
1 | x0 = input('Inserisci angolo di partenza in gradi')
|
La funzione input('testo') mostra il testo scelto e aspetta che inseriate un valore e premiate invio; il valore inserito verra' memorizzato nella variabile x0. Lo stesso metodo lo si puo' usare per inserire tutti i parametri che si vogliono.
Di seguito, invece, vediamo come usare il costrutto condizionale if per scegliere quale dei tre schemi di integrazione usare:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | method = input('0 - Eulero; 1 - Cromer; 2 - Verlet')
for t_step in range(N):
if method == 1:
#codice Eulero
else if method == 2:
#codice Eulero-Cromer
else if method == 3:
#codice Verlet
x0=x1
v0=v1
...
...
|
- nella prima riga supponiamo di chiedere all'utente un numero per identificare il metodo da usare;
- dentro il ciclo for usiamo il costrutto if method == valore:, vogliamo che il programma confronti la variabile method con un valore e se questa e' uguale allora esegue le righe di codice indentate che si trovano sotto;
- poiche' method ha un solo valore, una sola delle tre condizioni sara' verificata, quindi uno solo dei tre metodi verra' di fatto eseguito;
- al posto del commento andranno inserite le righe di codice relative al metodo viste all'inizio del tutorial;
- subito dopo il costrutto if, possiamo lasciare tutte le righe gia usate in precedenza, poiche' sono comuni ai tre metodi.