Integrazione delle equazioni del moto

In generale, le equazioni del moto della meccanica newtoniana si presentano nella forma di sistemi di equazioni differenziali del secondo ordine in cui le derivate seconde delle variabili dinamiche (le coordinate) sono espresse mediante il rapporto tra la legge di forza e la massa del corpo.

La più generale legge di forza compatibile con un sistema di equazioni in forma normale (in cui cioè le accelerazioni non dipendono da derivate delle coordinate di ordine superiore al primo) esprime le accelerazioni come funzioni (in genere non lineari) delle posizioni e velocità di tutti i gradi di libertà ed eventualmente funzione esplicita del tempo (il caso di forzanti esterne).

Formalmente significa lavorare con un sistema del tipo

    x1(t)=a1(x1(t),…  ,xN(t),ẋ1(t),…  ,ẋN(t),t) 
    x2(t)=a2(x1(t),…  ,xN(t),ẋ1(t),…  ,ẋN(t),t) 
    … 
    xN(t)=aN(x1(t),…  ,xN(t),ẋ1(t),…  ,ẋN(t),t)

Una volta imposte le condizioni iniziali (posizioni e velocità all’ istante iniziale t0):

    x1(t0)=x0(1)  
    ẋ1(t0)=vx0(1)  
    x2(t0)=x0(2)  
    ẋ2(t0)=vx0(2)  
    … 
    xN(t0)=x0(N)  
    ẋN(t0)=vx0(N),

sotto opportune condizioni sulle accelerazioni, sappiamo che esiste una soluzione unica per il problema.

Risolvere numericamente il problema delle equazioni del moto significa trasformare le equazioni differenziali in relazioni che permettano di calcolare l’ evoluzione di posizioni e velocità su piccoli intervalli discreti di tempo con accuratezza controllata. L’ idea di base è che se abbiamo uno schema numerico che permette di passare dalle condizioni iniziali alle posizioni e velocità al tempo t0+Δ t, questo potrà essere reiterato un numero arbitrario di volte. Naturalmente è cruciale che si possa tenere sotto controllo gli errori numerici così introdotti.

Esistono molti algoritmi per la soluzione di generici sistemi di equazioni differenziali. In questo corso ci limiteremo ad un unico algoritmo che ha alcune limitazioni (non consente di trattare direttamente forze dipendenti dalle velocità) ma che ha molti vantaggi uniti ad un’ estrema semplicità di derivazione e di implementazione.

L’ algoritmo è noto in contesti diversi con nomi diversi. Qui ci riferiremo ad esso come algoritmo di Verlet. Una derivazione elementare è la seguente. Per semplicità consideriamo il caso di un solo grado di libertà. L’ estensione a più gradi di libertà è banale.

Indichiamo per brevità con un indice 0 le quantità calcolate al tempo generico t e con indice + quelle calcolate al tempo t+Δ t. Abbreviamo anche Δ t con Δ. Se per piccoli intervalli di tempo Δ approssimiamo la forza con una forza costante, possiamo approssimare l’ evoluzione locale mediante le formule

     
    x+ =
 x0  + Δ v0  + 
Δ2
2
 a0  
 
    v+ = v0 + Δ a0       (1)

Questa discretizzazione delle equazioni del moto, nota come formula di Eulero, corrisponde a prendere i termini fino a quelli contenenti l’ accelerazione in un’ espansione in serie di Taylor della posizione e della velocità al tempo t+Δ attorno al tempo t.

La qualità dell’ integrazione delle equazioni fornita dall’ algoritmo di Eulero non è eccezionale. Uno dei problemi collegati a questo algoritmo è la mancanza di reversibilità per inversione temporale dello stesso nelle situazioni in cui ci si aspetterebbe che la dinamica esatta sia reversibile.

Possiamo misurare la mancata reversibilità immaginando di avere a che fare con interazioni conservative e non dipendenti dal tempo. Invertendo la velocità al tempo t+Δ e seguendo il sistema per un ulteriore intervallo Δ il sistema dovrebbe tornare al punto iniziale con la stessa velocità iniziale ma invertita. Se facciamo questo con le (1) troviamo (indicando con un apice le quantità ottenute quando il sistema “torna indietro”):

     
    x′0 =
 x0 + 
Δ2
2
 ( a+ − a0)  
 
    v′0 = −v0   − Δ ( a0 −a+ )        (2)

Per cercare di risolvere il problema (e sperabilmente migliorare l’ algoritmo) si può procedere nel seguente modo:

l’ equazione meno precisa è quella per le velocità. Infatti il termine successivo dell’ espansione in serie sarebbe proporzionale a Δ2. Aggiungiamo allora all’ equazione per la velocità un termine correttivo cercando di ripristinare anche per l’ algoritmo discretizzato l’ esatta simmetria per inversione temporale. Avremo allora:

     
    x+ =
 x0  + Δ v0  + 
Δ2
2
 a0  
 
    v+ = v0 + Δ a0 + Δ c0       (3)

Se adesso invertiamo le velocità e seguiamo il sistema per un ulteriore Δ troviamo i nuovi valori:

     
    x′0 =
 x+  − Δ v+  + 
Δ2
2
 a+  
 
    v′0 = −v+ + Δ a+ + Δ c+   

ovvero

     
    x′0 =
 x0 + 
Δ2
2
 (a+ − a0 ) −Δ2 c0  
    (4)
    v′0 = −v0 − Δ (a0 − a+) − Δ ( c0 − c+ )       (5)

se adesso richiediamo che x′0 = x0 e v′0 = −v0 abbiamo dalla (4) c0 = a+ − a0 /2 e dalla (5) c+ = a0 − a+ /2. Quindi l’ algoritmo (Verlet) è:

     
    x+ =
 x0  + Δ v0  + 
Δ2
2
 a0  
 
    v+ =
 v0 + 
Δ
2
 (a0 + a+)  
    (6)

che differisce da quello di Eulero solo per aver sostituito all’ accelerazione che appare nel termine di velocità il valor medio tra l’ accelerazione all’ istante iniziale (a0) e quella all’ istante t+Δ (a+). Da notare che l’ algoritmo resta esplicito perché l’ accelerazione a+ può essere calcolata appena è stata aggiornata la posizione al valore x+.

Questo spiega anche la limitazione a forze che non dipendano dalla velocità (anche se una generalizzazione è possibile per forze che dipendano linearmente dalla velocità).

Un’ analisi elementare dell’ errore introdotto dalla discretizzazione può essere svolta a partire dalla differenza tra le quantità ottenuta al tempo Δ mediante le (6) e lo sviluppo in serie di Taylor di posizione e velocità. In tal modo si trova che deve valere:

     
    x(t+Δ) − x+ = A Δ3 ȧ(τ)   
    v(t+Δ) − v+ = B Δ3  a(τ′)       (7)

Dove A e B sono due costanti numeriche e τ e τ′ sono due istanti di tempo appartenenti all’ intervallo [t,t+Δ ]. Se integriamo la traiettoria per un intervallo di tempo finito il numero di step di integrazione sarà inversamente proporzionale a Δ e quindi l’ errore globale su un tempo finito sarà maggiore di un ordine di grandezza rispetto a quello locale.

This document was translated from LATEX by HEVEA.