Struttura del programma

Come indicato sopra, il potenziale $V(x)$ è definito dalla subroutine potential separata, e quindi è in un certo senso esterno al programma principale. Per un dato $x$, potential ritorna il potenziale $V(x)$ e la corrispondente forza $F(x)=-dV/dx$. Solo quest'ultima è, strettamente parlando, necessaria per portare avanti l'integrazione; ma conoscere anche il potenziale permette di verificare la legge di conservazione dell'energia passo per passo.

Il programma chiede i seguenti dati di input: time step $\Delta t$, posizione iniziale $x_1$, energia totale $E_1$, numero di time steps (verrà seguita l'evoluzione temporale per un tempo totale $N\Delta t$), e nome del file di uscita (è anche possibile convogliare l'output sullo standard output).

L'energia totale serve a calcolare la velocità iniziale $v_1$, secondo

$$E_1 = \frac{1}{2} m v_1^2 + V(x_1)$$ (1.7)

Naturalmente è necessario che sia $E_1 \ge V(x_1)$ (se così non è il programma si arresta). È stato scelto di usare l'energia, piuttosto che direttamente la velocità iniziale, perchè così è immediatamente chiaro quale regione del potenziale verrà esplorata nel corso della simulazione: quella in cui $V(x)\le E_1$. I punti in cui $V = E_1$ sono i punti di inversione, in cui la velocità è nulla e la particella inverte la sua direzione di moto. Inoltre, come ben noto l'energia totale

$$E = \frac{1}{2} m v^2 + V(x)$$ (1.8)

è supposta essere una costante del moto. La conservazione dell'energia costituisce per questi problemi un eccellente strumento di verifica che non vi siano errori grossolani o siano stati dati input errati, come un $\Delta t$ troppo grande.

Assegnare la velocità iniziale significa, nell'ottica dell'algoritmo di Størmer-Verlet, assegnare la posizione al "tempo precedente". Tale algoritmo però opera su posizioni ed accelerazioni. La velocità alla posizione $n$ può essere in generale ottenuta come

$$v_n = \frac{x_{n+1} - x_{n-1}}{2\Delta t}$$ (1.9)

con un errore dell'ordine di $(\Delta t)^2$ (come subito si vede dalle (1.2), sottraendo un'equazione dall'altra). Tuttavia, questa formula coinvolge tre punti e come condizioni iniziali noi abbiamo a disposizione la sola posizione e velocità iniziali al medesimo punto $x_1$. Si è allora scelto di utilizzare l'espressione

$$v_1 \simeq \frac{x_1 - x_0}{\Delta t}$$ (1.10)

per definire la "posizione precedente" $x_0$ a partire da $x_1$ e $v_1$. Se $x_0$ e $x_1$ fossero dati, ad una tale formula sarebbe associato un errore dell'ordine di $\Delta t$ nella velocità, quindi non si tratta di una buona stima. Tuttavia, quello che al più succede è che una volta che la simulazione è partita, l'energia totale si assesterà ad un valore lievemente diverso da quello da noi richiesto. In uno schema più raffinato, si potrebbe pensare di effettuare una correzione al secondo passo di integrazione per ottenere una migliore corrispondenza fra energia richiesta e energia effettiva.

La velocità iniziale è sempre definita positiva dal programma, ossia il punto inizia a muoversi verso gli $x$ positivi. Data la simmetria del potenziale, questa non è una limitazione (invertire la velocità è equivalente a lasciarla invariata e cambiare il segno della coordinata di posizione iniziale).

Il loop sul numero di time steps è la parte centrale del programma. Il punto su cui va posta particolare attenzione è il trasferimento dell'informazione della "posizione precedente" da una iterazione alla successiva, effettuato usando tre variabili x, xprev, xprevsav.

A fini di controllo sulla conservazione dell'energia, viene calcolata anche la velocità e l'energia cinetica. Per ottenere una precisione accettabile, le velocità vengono definite utilizzando la (1.9). Questo però comporta che per definire $v_n$ è necessario conoscere $x_{n+1}$, e quindi che alla conclusione di ogni iterazione le velocità calcolate sono riferite all'iterazione precedente, ossia alla posizione che aveva la particella all'inizio dell'iterazione, prima di essere spostata.

I dati emessi sul file di output sono, nell'ordine, il numero dell'iterazione, la posizione, la velocità, la forza, l'energia cinetica, l'energia potenziale e l'energia totale: tutte riferite alla posizione che aveva la particella all'inizio dell'iterazione.