L’idea è di scrivere un codice per risolvere l’equazione differenziale dell’oscillatore armonico generica, implementando i metodi di : Eulero, Eulero-Cromer e Verlet.
Gli schemi di integrazione sono:
- Eulero: \(v_{n+1}=v_{n}+\tau F_{n} ; x_{n+1}=x_{n}+\tau v_{n}\)
- Eulero-Cromer: \(v_{n+1}=v_{n}+\tau F_{n} ; x_{n+1}=x_{n}+\tau v_{n+1}\)
- Verlet: \(x_{n+1}=x_{n}+\tau v_{n}+F_{n} \frac{\tau^{2}}{2} ; v_{n+1}=v_{n}+\frac{\tau^{2}}{2}\left(F_{n}+F_{n+1}\right) \)
dove con \(\tau\) indichiamo il passo temporale e con \(F_n\) indichiamo l’accelerazione, o la forza diviso la massa, a cui è sottoposto il nostro sistema. Nel nostro caso avremo:
\(\frac{d^{2} \theta}{d t^{2}}=-\frac{g}{l} \sin \theta=-\omega^{2} \sin \theta\)
La traduzione in codice di questi schemi è immediata. Vediamola di seguito e poi ci costruiamo il programma intorno:
#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)
Nel codice, facciamo uso di due variabili una denominata con 0 l’altra con 1 che si riferiscono al tempo \(n\) e \(n+1\) rispettivamente. Se ora si definiscono le costanti necessarie (g, l, τ) e le condizioni iniziali (x0,v0), 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:
import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi
g = 1
l = 1
x0 = np.deg2rad(5)
v0 = 0.0
tau = 0.01
Facciamo notare 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 più lungo di in solo passo τ, per esempio per T = N · τ, con N a piacere. E facile capire che dovremo ripetere lo schema di integrazione scelto per un numero N di volte, mediante un ciclo for,
come vediamo di seguito:
N = 3000
for t_step in range(N):
#Accelerazione
F0 = -g/l*sin(x0)
#Eulero-Cromer
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’uscita dal ciclo for avremo a disposizione solo l’ultimo 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 array:
N = 3000
x_vet = np.zeros(N+1)
v_vet = np.zeros(N+1)
• registriamo nel primo elemento degli array le condizioni iniziali:
x_vet[0] = x0
v_vet[0] = v0
• all’interno del ciclo for registriamo ogni nuovo valore di x e v all’interno dei rispettivi array:
for timestep in range(N):
#Accelerazione
F0 = -g/l*sin(x0)
x1 = x0+v0*tau+0.5*F0*tau**2
F1 = -g/l*sin(x1)
v1 = v0+0.5*tau*(F0+F1)
x0 = x1
v0 = v1
x_vet[timestep+1] = x0
v_vet[timestep+1] = v0
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():
time_vet = np.linspace (0,N*tau, N+1)
plt.xlabel('Tempo (sec)')
plt.ylabel('Ampiezza (rad)')
plt.plot(time_vet, x_vet)
plt.show()
Ecco il codice completo:
import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi
g = 1
l = 1
x0 = np.deg2rad(5)
v0 = 0.0
tau = 0.01
N = 3000
x_vet = np.zeros(N+1)
v_vet = np.zeros(N+1)
x_vet[0] = x0
v_vet[0] = v0
for timestep in range(N):
#Accelerazione
F0 = -g/l*sin(x0)
x1 = x0+v0*tau+0.5*F0*tau**2
F1 = -g/l*sin(x1)
v1 = v0+0.5*tau*(F0+F1)
x0 = x1
v0 = v1
x_vet[timestep+1] = x0
v_vet[timestep+1] = v0
time_vet = np.linspace (0,N*tau, N+1)
plt.xlabel('Tempo (sec)')
plt.ylabel('Ampiezza (rad)')
plt.plot(time_vet, x_vet)
plt.show()
e questo è il diagramma ottenuto:
checksum: 38957cccdef0229b17dfce1dfb41c983f71f59d88a0e338d8718d21943258cae965609e9d16032c9fa87a0271e120d2c465784b7048fa0c9cd7dfecbe94e7680