1Introduzione¶
Per affrontare la complessità dei problemi che caratterizzano la fisica moderna esistono diverse strategie, distribuite lungo uno spettro continuo che unisce la formalizzazione puramente analitica alla risoluzione numerica di forza bruta. Molto spesso la ricerca si colloca in una posizione intermedia, adottando approcci teorico-computazionali ibridi: il modello fisico viene inizialmente semplificato attraverso opportune ipotesi teoriche, per poi validare la soluzione analitica tramite il calcolo numerico. Ma come si traduce, concretamente, un problema fisico in termini computazionali? Sebbene la risposta dipenda strettamente dalla natura del sistema in esame, la fisica computazionale ha sviluppato metodologie di carattere generale applicabili a vastissime classi di fenomeni, dalla meccanica quantistica all’astrofisica. In questa sezione, in particolare, analizzeremo i metodi per l’integrazione numerica di equazioni differenziali, uno strumento pilastro per l’indagine scientifica in ogni ambito della fisica contemporanea[1].
2Equazioni differenziali ordinarie¶
Moltissimi problemi di fisica si possono formalizzare in termini di equazioni differenziali, cioè relazioni che connettono una funzione incognita alle sue derivate. La maggior parte delle equazioni differenziali di interesse non possono essere risolte analiticamente, e richiedono quindi di essere affrontate con metodi numerici. In questo corso ci occuperemo principalmente delle cosiddette equazioni differenziali ordinarie (spesso chiamate ODE, per ordinary differential equations), in cui la funzione incognita è di una sola variabile.
Un’equazione differenziale si dice dell’ordine -esimo se contiene al suo interno derivate -esime della funzione incognita. La soluzione generale di una ODE di ordine -esimo contiene costanti di integrazione indipendenti il cui valore viene fissato specificando condizioni (dette solitamente condizioni iniziali o al bordo) per ottenere una soluzione particolare.
Nel seguito ci occuperemo principalmente di metodi per la risoluzione numerica di ODE del primo ordine, cioè del tipo
dove è la funzione incognita, è la sua derivata prima e è la variable indipendente. Se siete studenti di fisica a Sapienza e state affrontando per la prima volta questo corso, avrete probabilmente già seguito (e magari passato con profitto) il corso di Meccanica. Durante il corso avrete senz’altro incontrato equazioni differenziali del secondo ordine del tipo
dove è la posizione di un punto materiale, la sua accelerazione, la forza a cui è sottoposto e la sua massa[2]. Ricordando il legame tra velocità e posizione, un’equazione del tipo (2) può essere trasformata in nel seguente sistema di equazioni del primo ordine:
cui di solito si affiancano le condizioni iniziali e .
2.1Qualche esempio di sistemi non risolvibili analiticamente¶
Quando si inizia lo studio della fisica teorica, si ha spesso l’illusione che ogni sistema fisico descrivibile tramite equazioni di Newton o di Lagrange possa essere risolto “con carta e penna”, trovando una formula esatta per la traiettoria nel tempo. La realtà, purtroppo, è ben diversa: i sistemi integrabili analiticamente rappresentano una piccolissima eccezione in un oceano di problemi matematicamente intrattabili.
Per capire quanto sia facile imbattersi in equazioni prive di soluzioni analitiche, proviamo a scendere dal complesso al semplice, partendo da un sistema apparentemente elementare: il doppio pendolo.
2.1.1Il doppio pendolo: il regno del caos¶
Immaginiamo di appendere un pendolo rigido (la cui lunghezza, cioè, rimane costante) all’estremità di un altro. Questo sistema, composto da due aste rigide di lunghezza e due masse vincolate a muoversi su un piano verticale, ha solo due gradi di libertà, rappresentati dagli angoli e che le aste formano con la verticale.
Nonostante l’apparente semplicità costruttiva, la dinamica del sistema è determinata da un sistema di due equazioni differenziali del secondo ordine fortemente accoppiate e non lineari:
dove e . Queste equazioni sono impossibili da risolvere in forma chiusa. Non solo: il doppio pendolo è uno dei più celebri esempi di sistema caotico, una proprietà che discuteremo meglio più avanti. Qui basti sapere che con “sistema caotico” si intende un sistema per cui una piccolissima variazione nelle condizioni iniziali o (anche solo dovuta alla precisione finita con cui un computer immagazzina i numeri decimali) produce traiettorie che divergono completamente[3]. Per studiarne la dinamica, l’integrazione numerica al computer non è un’opzione comoda, è l’unica via percorribile. Un esempio di simulazione è mostrato in Figura 1.
Figure 1:Simulazione di un pendolo doppio di parametri m, Kg e Kg e condizioni iniziali , , .
2.1.2Il pendolo semplice¶
Si potrebbe pensare che la difficoltà del doppio pendolo derivi esclusivamente dall’accoppiamento tra i suoi due gradi di libertà. Eliminiamo allora il secondo pendolo, ponendo formalmente . Rimane il pendolo semplice: una massa vincolata a muoversi lungo un arco di circonferenza di raggio . L’equazione del moto è
Il sistema possiede un solo grado di libertà e non presenta alcun accoppiamento. Tuttavia, l’equazione è ancora non lineare a causa del termine . Questa non linearità non rende il problema caotico ma impedisce in generale di esprimere il moto mediante le sole funzioni elementari.
Per comprenderne l’origine, moltiplichiamo l’equazione per :
Integrando rispetto al tempo si ottiene la conservazione dell’energia:
Sia l’ampiezza massima dell’oscillazione. Nel punto di inversione il pendolo è istantaneamente fermo, quindi
e pertanto
Segue che
Questa relazione consente di determinare il tempo mediante un’integrazione. In particolare, il tempo necessario affinché il pendolo vada dall’ampiezza massima alla posizione di equilibrio è pari a un quarto del periodo:
L’integrale a destra è un integrale ellittico di prima specie. Non esiste alcuna manipolazione algebrica o sostituzione trigonometrica in grado di risolverlo usando le funzioni standard (come logaritmi, esponenziali, seni o coseni). Di fatto, le cosiddette funzioni ellittiche usate in matematica avanzata sono definite proprio a partire da questo tipo di integrali, il che equivale a dire che dobbiamo “inventarci” delle nuove funzioni per descrivere la soluzione.
Derivazione dell’integrale ellittico
2.1.3Le piccole oscillazioni¶
Com’è possibile, allora, che in tutti i corsi di fisica scolastici e universitari di base si impari a risolvere il pendolo con una semplice funzione trigonometrica?
Ciò è possibile solo introducendo un’approssimazione fisica cruciale: l’ipotesi di piccole oscillazioni. Se limitiamo lo studio a angoli molto piccoli ( radiante, indicativamente sotto i ), possiamo sviluppare in serie di Taylor la funzione seno attorno a zero, arrestandoci al primo ordine:
Sotto questa assunzione, l’equazione del moto perde la sua natura non lineare e si trasforma in un’equazione differenziale lineare a coefficienti costanti:
Questa equazione è finalmente risolvibile con carta e penna, e la sua soluzione generale è una semplice oscillazione armonica di frequenza :
In questa approssimazione, l’integrale ellittico nell’equazione (11) si riduce a una costante (), e il periodo diventa
Questo modello lineare, noto come oscillatore armonico, è uno dei pilastri della fisica proprio perché rappresenta il “porto sicuro” in cui i fisici rifugiano ogni volta che un sistema non lineare diventa matematicamente inaffrontabile. Ed è proprio dall’oscillatore armonico che partiremo per testare e confrontare i nostri algoritmi di integrazione numerica.
2.2Il sistema modello per definizione: l’oscillatore armonico¶
Abbiamo appena visto come il pendolo in regime di piccole oscillazioni possa essere approssimato con un oscillatore armonico unidimensionale. Nel seguito, come esempio di sistema dinamico utilizzeremo proprio questo modello che, oltre a descrivere direttamente numerosi fenomeni fisici, presenta il vantaggio di possedere una soluzione analitica semplice, che potrà essere utilizzata per valutare l’accuratezza dei diversi algoritmi di integrazione numerica.
Consideriamo una particella di massa soggetta a una forza elastica proporzionale allo spostamento dalla posizione di equilibrio,
dove è la costante elastica della molla. Applicando la seconda legge di Newton, e usando la notazione e per indicare le derivate prime e seconde, rispettivamente, otteniamo
ovvero
dove abbiamo introdotto la pulsazione naturale del sistema
che implica come il periodo del moto sia
La soluzione generale di questa equazione è
oppure, in forma equivalente,
dove le costanti , (oppure e ) sono determinate dalle condizioni iniziali, e . Poiché si tratta di un sistema senza attrito sottoposto a una forza che non dipende esplicitamente dal tempo, l’energia totale, somma di energia potenziale ed energia cinetica , si conserva:
Come accennato precedentemente, per risolvere numericamente l’equazione differenziale del secondo ordine (23) conviene trasformarla nel seguente sistema di due equazioni del primo ordine:
Nel resto del capitolo considereremo anche una versione più generale del problema, che include sia l’attrito viscoso sia una forzante, ovvero una forza esterna dipendente dal tempo:
A seconda della scelta dei parametri si ottengono diversi casi di interesse fisico:
e : oscillatore armonico semplice;
e : oscillatore armonico smorzato;
e : oscillatore armonico forzato.
Questo sistema costituirà il principale banco di prova per gli algoritmi di integrazione numerica discussi nelle sezioni successive.
2.3Dal continuo al discreto¶
I computer sono macchine discrete e, come tali, non possono rappresentare esattamente quantità continue, ma soltanto approssimarle mediante un numero finito di valori. Quando vogliamo studiare numericamente un sistema descritto da equazioni differenziali, dobbiamo quindi trovare un modo per tradurre un problema continuo in una forma compatibile con l’architettura discreta del calcolatore.
Esistono diverse strategie per affrontare questo problema; in queste note ci concentreremo sui cosiddetti metodi delle differenze finite (finite difference methods)[5]. Per semplicità cominciamo la trattazione considerando casi unidimensionali, in cui la funzione incognita è . L’idea fondamentale consiste nel sostituire il dominio continuo della variabile indipendente (ad esempio il tempo ) con una successione discreta di punti separati da un intervallo . In altre parole, invece di descrivere l’evoluzione del sistema in ogni istante, ne consideriamo soltanto una sequenza di “fotogrammi” successivi. Senza perdità di generalità, consideriamo un intervallo temporale e suddividiamolo in intervalli uguali. Definiamo
e i punti della griglia
Nel seguito per alleggerire la trattazione utilizzeremo spesso la notazione .
Una volta discretizzato il dominio, le derivate che compaiono nelle equazioni differenziali possono essere approssimate mediante opportune differenze tra i valori della funzione nei punti della griglia. Ad esempio, la derivata prima di nel punto può essere approssimata come
In questo modo un’equazione differenziale viene trasformata in una relazione algebrica tra valori della funzione calcolati in istanti successivi. Per comprendere l’idea alla base tutti i metodi che introdurremo nelle prossime sezioni, integriamo entrambi i membri dell’equazione (1) tra due istanti consecutivi della griglia temporale, e , ottenendo
ovvero
Introducendo il passo temporale , possiamo riscrivere questa espressione come
dove
rappresenta il valore medio di nell’intervallo .
La relazione (36) è esatta e costituisce il punto di partenza di tutti i metodi di integrazione numerica che vedremo. La difficoltà risiede nel fatto che, in generale, il valore medio non è noto, poiché dipende dall’andamento della soluzione all’interno dell’intervallo stesso. I diversi algoritmi che presenteremo possono essere interpretati come diversi modi di approssimare questa quantità.
3Eulero e Eulero-Cromer¶
3.1Formulazione¶
Il metodo di Eulero è il più semplice algoritmo di integrazione numerica per equazioni differenziali ordinarie. L’idea consiste nell’approssimare il valore medio della funzione nell’intervallo con il suo valore all’inizio dell’intervallo:
Sostituendo questa approssimazione nell’equazione (36) si ottiene
L’algoritmo può quindi essere interpretato come un’estrapolazione lineare della soluzione a partire dalla sua derivata nel punto iniziale dell’intervallo.
Nel caso di ODE del secondo ordine, come l’oscillatore armonico, le equazioni di aggiornamento diventano
dove . Abbiamo quindi approssimato sia l’accelerazione media sia la velocità media nell’intervallo utilizzando i rispettivi valori all’inizio dell’intervallo, cioè
Il metodo di Eulero è semplice da implementare, ma la sua accuratezza è limitata e può produrre risultati qualitativamente scorretti quando viene applicato a sistemi oscillanti per tempi lunghi. Un miglioramente a volte sostanziale si può ottenere utilizzando il metodo di Eulero-Cromer, che è una semplice modifica del metodo di Eulero particolarmente adatta allo studio di sistemi meccanici.
Dal punto di vista concettuale, il metodo di Eulero-Cromer utilizza il valore iniziale dell’accelerazione per stimare l’accelerazione media (, come Eulero), ma il valore finale della velocità per stimare la velocità media (). L’algoritmo completo assume pertanto la forma
Questa semplice modifica produce risultati significativamente migliori in molti problemi meccanici, in particolare nei sistemi oscillanti.
3.2C: Raggruppare i dati con struct e definire nuovi tipi con typedef¶
Per simulare l’oscillatore armonico dobbiamo conservare diverse quantità. Alcune descrivono lo stato del sistema in un particolare istante, come posizione e velocità; altre rimangono costanti durante la simulazione, come la massa, la costante elastica e il passo temporale.
Finora potremmo rappresentare queste quantità mediante variabili indipendenti:
double x;
double v;
double k;
double m;
double dt;Questo approccio funziona, ma non rende esplicito il fatto che x e v descrivono insieme lo stato di un unico punto materiale, mentre k, m e dt caratterizzano il sistema che stiamo simulando.
Il C permette di raggruppare variabili logicamente collegate mediante una struttura, definita dalla parola chiave struct:
struct Punto {
double x;
double v;
};Abbiamo così definito un nuovo tipo di struttura chiamato struct Punto.
Possiamo usarlo per dichiarare una variabile che contiene contemporaneamente
posizione e velocità:
struct Punto p;
p.x = 2.0;
p.v = 1.0;Le variabili contenute in una struttura prendono il nome di membri o campi. Per accedere a un campo si usa l’operatore punto (.): p.x indica la posizione e p.v la velocità.
3.2.1Dare un nome a un tipo: typedef¶
In C, typedef permette di assegnare un nuovo nome a un tipo già esistente.
Per evitare di dover scrivere ogni volta struct Punto, possiamo definire:
typedef struct Punto Punto;Da questo momento Punto è un sinonimo di struct Punto, e possiamo scrivere
più semplicemente:
Punto punto;È molto comune unire la definizione della struttura e il typedef in un’unica
dichiarazione:
typedef struct {
double x;
double v;
} Punto;La parte racchiusa tra parentesi graffe descrive la struttura del dato, mentre
il nome dopo la parentesi graffa, Punto, è il nuovo tipo che stiamo
definendo.
Possiamo anche inizializzare tutti i campi al momento della dichiarazione:
Punto punto = {
.x = 2.0,
.v = 1.0
};Le espressioni .x = 2.0 e .v = 1.0 sono dette inizializzatori
designati. Rendono esplicito quale valore viene assegnato a ciascun campo.
3.2.2Separare lo stato dal sistema¶
Definiamo ora una seconda struttura che contenga i parametri necessari per descrivere l’oscillatore armonico e la sua integrazione numerica:
typedef struct {
// a scelta, potremmo scrivere il problema in funzione della pulsazione
// e tenere in memoria omega_0 invece di k
double k;
double m;
double dt;
} Sistema;Possiamo quindi inizializzare l’intera simulazione nel modo seguente:
Punto punto = {
.x = 2.0,
.v = 1.0
};
Sistema sistema = {
.k = 1.0,
.m = 1.0,
.dt = 0.01
};Le due strutture hanno ruoli diversi:
Puntocontiene lo stato dinamico, che cambia a ogni passo temporale;Sistemacontiene i parametri fisici e numerici, che rimangono costanti durante la simulazione.
3.3Passare una struttura a una funzione¶
Le strutture possono essere passate alle funzioni come le variabili di tipo fondamentale. Per esempio, possiamo scrivere una funzione che calcoli l’accelerazione:
double accelerazione(Sistema sistema, Punto punto) {
return -(sistema.k / sistema.m) * punto.x;
}Questa funzione riceve però delle copie delle due strutture. Per evitare la copia possiamo passarne gli indirizzi:
double accelerazione(Sistema *sistema, Punto *punto) {
return -(sistema->k / sistema->m) * punto->x;
}Quando possediamo un puntatore a una struttura, accediamo ai suoi campi tramite l’operatore freccia (->). Per esempio,
punto->xè una forma più leggibile ma del tutto equivalente a
(*punto).xPossiamo ora riscrivere un passo del metodo di Eulero facendo sì che la funzione modifichi direttamente lo stato del punto:
void eulero(Sistema *sistema, Punto *punto) {
double a = accelerazione(punto, sistema);
punto->x += punto->v * sistema->dt;
punto->v += a * sistema->dt;
}Il metodo di Eulero-Cromer ha la stessa interfaccia e differisce soltanto nell’ordine degli aggiornamenti:
void eulero_cromer(Punto *punto, const Sistema *sistema) {
double a = accelerazione(punto, sistema);
punto->v += a * sistema->dt;
punto->x += punto->v * sistema->dt;
}Il ciclo temporale diventa così particolarmente compatto:
double t = 0.0;
Punto punto;
Sistema sistema;
/* inizializzazione */
while(t < t_max) {
printf("%g %g %g\n", t, punto.x, punto.v);
eulero_cromer(&sistema, &punto);
t += sistema.dt;
}Raggruppare i dati in strutture non cambia l’algoritmo, ma rende più chiaro quali quantità appartengano allo stato del sistema e quali ne definiscano le proprietà. Inoltre, per aggiungere nuove coordinate o nuovi parametri sarà sufficiente estendere la struttura corrispondente, senza modificare la firma di tutte le funzioni che la utilizzano.
3.4Esempio di integrazione numerica¶
Figure 2:Il risultato dell’integrazione del sistema (29) con il metodo di Eulero. Dall’alto verso il basso, i tre pannelli mostrano la posizione , la velocità e l’energia meccanica in funzione del tempo per tre diversi valori del passo temporale (10-1 in blu, 10-2 in arancione e 10-3 in verde), oltre al risultato esatto (riga tratteggiata). Considerando, per comodità, grandezze adimensionali, il sistema simulato ha (e quindi ) e, come condizioni iniziali, e .
Applichiamo i due metodi appena introdotti al sistema di equazioni differenziali (29), cercando di valutare la qualità della soluzione numerica discretizzata al variare della grandezza del passo temporale .
Cominciamo ad analizzare i risultati ottenuti con il metodo di Eulero, mostrati in Figura 2. Notiamo prima di tutto che solo le curva verdi (relative a ) sembrano ricalcare fedelmente, almeno alla scala della figura, la soluzione teorica. Per valori maggiori di tutte le quantità mostrate si discostano anche sensibilmente dalla teoria. È preoccupante non tanto il fatto che ci sia una differenza tra i valori numerici e quelli teorici, quanto che questa differenza aumenti nel tempo. Infatti, una delle principali proprietà dell’oscillatore armonico è la sua periodicità: il moto si ripete esattamente ogni periodo . Come si può vedere dalla figura, questa proprietà non è affatto rispettata dalla soluzione ottenuta con il metodo di Eulero: le oscillazioni di posizione e velocità aumentano di ampiezza col tempo. Questo aumento si riflette nell’energia totale, che a sua volta aumenta monotonicamente: l’errore dovuto alla discretizzazione ha l’effetto netto di immettere energia nel sistema.
Passiamo ad analizzare i risultati ottenuti con Eulero-Cromer e mostrati in Figura 3. Nonostante l’apparente similitudine dei due metodi, il comportamento che si osserva è molto diverso. In questo caso posizione e velocità sembrano venir riprodotte quasi perfettamente per tutti i valori di , almeno alla scala della figura[6]. Per quanto riguarda l’energia, questa sembra comportarsi in una maniera più strana: in tutti i casi (anche se, per , non si vede bene) non è costante nel tempo ma oscilla con periodo uguale a quello di e e ampiezza che decresce al diminuire di . Quindi, se da un lato è vero che l’energia non si conserva, il suo valore medio rimane costante nel tempo: non c’è immissione o dissipazione netta di energia. Questa proprietà di “conservazione media” dell’energia è il massimo che possiamo chiedere a un algoritmo di integrazione numerico.
Figure 4:Simulazione di un oscillatore armonico integrato con Eulero (pallina rossa) ed Eulero-Cromer (pallina blu). I parametri della simulazione sono s, m, m/s e s.
La Figura 4 contiene una simulazione interattiva che mostra come Eulero, a differenza di Eulero-Cromer, non riesca a riprodurre la periodicità dell’oscillatore armonico: si vede chiaramente come l’energia del sistema aumenti via via che il tempo passa, mostrando un comportamento evidentemente non fisico.
Il confronto fatto tra i risultati ottenuti con Eulero ed Eulero-Cromer ci permette di introdurre due proprietà fondamentali degli algoritmi per l’integrazione numerica: stabilità e accuratezza. Questi due concetti non sono necessariamente legati: un algoritmo può essere poco stabile ma molto accurato, un altro molto stabile ma poco accurato.
3.5Stabilità¶
Un algoritmo di integrazione numerica si dice stabile se piccoli errori introdotti durante l’evoluzione (dovuti, ad esempio, all’approssimazione del metodo o all’arrotondamento numerico, che sono fonti di errore sempre presenti su un calcolatore) non vengono amplificati in modo incontrollato al procedere dei passi temporali. Un metodo si dice incondizionatamente stabile se rimane stabile per qualunque valore del passo temporale . Si dice invece condizionatamente stabile se la stabilità è garantita solo quando soddisfa una certa condizione, ad esempio . Sia la proprietà di essere condizionatamente/incodizionatamente stabile che l’eventuale valore di dipendono non solo dall’algoritmo, ma anche dal problema che vogliamo risolvere. Vediamo come studiare la stabilità nel caso dell’oscillatore armonico, un sistema lineare che rende questo tipo di analisi più trasparente.
Le equazioni di aggiornamento del metodo di Eulero, Eq. (40), possono essere riscritte per l’oscillatore armonico come
Introduciamo ora il concetto di spazio delle fasi: questo è l’insieme di tutte le possibili configurazioni (o stati) del sistema. Nel caso dell’oscillatore armonico unidimensionale[7], per identificare una configurazione è sufficiente specificare posizione e velocità , e quindi lo spazio delle fasi comprende l’intero piano . Un punto su questo piano, cioè una configurazione del sistema, si può identificare tramite un vettore . Discretizzando la notazione, possiamo definire lo stato del sistema al generico tempo , , così da poter riscrivere il passo di integrazione temporale (43) in forma compatta:
dove
Utilizzando questo formalismo possiamo scrivere direttamente l’evoluzione del sistema dalle condizioni iniziali ad un generico tempo come
Invece di calcolare la potenza -esima di componente per componente, possiamo utilizzare la decomposizione spettrale in autovalori e autovettori per ottenere direttamente l’operatore che determina l’evoluzione del sistema al tempo voluto. Poiché è una matrice , essa ammette due autovalori e , ai quali corrispondono due autovettori linearmente indipendenti e , tali per cui:
Poiché i due autovettori formano una base dello spazio delle fasi, possiamo esprimere qualsiasi condizione iniziale come una loro combinazione lineare:
dove e sono coefficienti (in generale complessi) che dipendono dallo stato iniziale scelto. Sfruttando la linearità della matrice , l’applicazione ripetuta dell’operatore di evoluzione per passi si riduce a
Se autovalori e autovettori sono complessi, e sono reali?
Un dubbio legittimo sorge spontaneo: se gli autovalori e gli autovettori sono numeri complessi, come fa lo stato fisico del sistema a rimanere composto da coordinate puramente reali (posizione e velocità) a ogni passo?
La risposta risiede in una proprietà fondamentale delle matrici reali. Dimostriamolo in tre passi:
Poiché la matrice di evoluzione ha elementi puramente reali, se i suoi autovalori non sono reali allora sono complessi coniugati. Infatti, se ammette un autovalore complesso , anche il suo complesso coniugato deve essere un autovalore. Se è l’autovettore associato a (ovvero ), coniugando entrambi i membri otteniamo:
Questo mostra che l’autovettore associato a è esattamente il coniugato del primo, cioè .
Esprimiamo la condizione iniziale reale nella base degli autovettori:
Poiché è reale, deve valere . Coniugando l’espressione sopra si ottiene . Uguagliando le due relazioni e sfruttando l’indipendenza lineare di e , deduciamo che i coefficienti devono essere l’uno il coniugato dell’altro, cioè . Possiamo quindi definire e
Sostituiamo ora queste relazioni nella formula generale dell’evoluzione al passo , Eq. (46), ottenendo
Poiché il prodotto di coniugati è il coniugato del prodotto, il secondo termine non è altro che il complesso coniugato del primo: . L’equazione diventa quindi:
Ricordando l’identità algebrica per cui la somma di un numero complesso e del suo coniugato è pari a due volte la sua parte reale (), arriviamo al risultato finale:
Poiché la parte reale di qualunque quantità è, per definizione, un numero reale, lo stato del sistema è garantito essere reale ad ogni istante di tempo.
Poiché per definizione di autovettore si ha , otteniamo l’espressione formale per lo stato del sistema al passo :
Per comprendere a fondo il comportamento di questa equazione senza dover calcolare immediatamente e , analizziamo il sistema da una prospettiva geometrica e strutturale, partendo dal determinante della matrice di evoluzione, che ha un significato geometrico profondo: rappresenta il fattore di scala con cui vengono modificate le aree (o i volumi) nello spazio delle fasi.
Consideriamo prima di tutto l’effetto che l’evoluzione temporale discreta ha sulla propagazione degli errori. Immaginiamo che a un certo passo l’elaboratore introduca un piccolissimo errore di arrotondamento sullo stato del sistema (ad esempio, a causa della rappresentazione a precisione finita dei numeri in virgola mobile, che in doppia precisione hanno errori tipici dell’ordine di ).
Lo stato numerico reale al generico tempo diventa . Se decomponiamo questa perturbazione microscopica nella base degli autovettori possiamo scrivere
Dopo passi di calcolo, l’errore iniziale si sarà evoluto in:
Ipotizziamo che , e consideriamo il caso . In queste condizioni, anche se l’errore iniziale è microscopicamente irrilevante (per esempio ), il fattore , che cresce esponenzialmente con , può portare il termine di errore a diventare dello stesso ordine di grandezza del segnale fisico (si veda il box qui sotto per una dimostrazione rigorosa). In questo regime, detto di instabilità, i risultati dell’integrazione numerica sono del tutto privi di senso.
Passiamo ora ad analizzare come la trasformazione determinata da agisce nello spazio delle fasi. Se consideriamo una regione di condizioni iniziali che racchiude un’area (ad esempio, un quadratino di stati possibili), dopo un passo di integrazione questa regione si deformerà in un parallelogramma la cui area sarà pari a:
Nei sistemi fisici reali conservativi, l’evoluzione temporale non espande né contrae lo spazio delle fasi. Questa proprietà geometrica fondamentale è nota in meccanica classica come teorema di Liouville. Affinché un algoritmo numerico sia un buon modello della fisica reale, deve rispettare questa struttura.
Nel caso di uno spazio delle fasi bidimensionale (come nel nostro caso), la conservazione dell’area corrisponde a richiedere che[9]
Vediamo ora come si collega la conservazione dell’area con il comportamento dei singoli stati descritto dall’equazione (55). Dall’algebra lineare sappiamo che il determinante di una matrice è pari al prodotto dei suoi autovalori:
Discutiamo prima il caso in cui l’equazione (63) non è rispettata. Se , aree (o volumi) dello spazio delle fasi si contraggono man mano che si evolvono nel tempo. In questo caso l’equazione (68) implica che almeno uno degli autovalori è minore di uno. Se l’altro ha modulo maggiore di uno si ricade nell’amplificazione dell’errore discussa prima. Se invece entrambi gli autovalori hanno modulo minore di uno, i termini tenderanno a zero per . L’evoluzione numerica smorzerà artificialmente le oscillazioni, comportandosi come se nel sistema fosse presente un attrito fittizio non fisico.
Di converso, se , almeno uno degli autovalori deve avere modulo maggiore di 1 per via dell’equazione (68), e darà quindi luogo ad un’espansione verso l’infinito di aree (o volumi) dello spazio delle fasi, oltre che ad un’amplificazione incontrollata degli errori. In questo caso, dell’energia viene immessa artificialmente nel sistema.
D’altro canto, se l’equazione (63) è rispettata, allora l’area occupata da un insieme di stati nello spazio delle fasi rimane rigorosamente costante nel tempo. La condizione è quindi una condizione necessaria per garantire la stabilità a lungo termine e la quasi-conservazione[10] dell’energia numerica, la cui violazione porta ad un’alterazione artificiale della fisica del sistema ad ogni passo temporale, con conseguenze più o meno gravi a seconda del sistema studiato. Questa proprietà geometrica è nota come simpletticità (e l’algoritmo di integrazione che ne è provvisto si dice simplettico).
Per un integratore simplettico, gli autovalori sono rigidamente vincolati dalla relazione . Questo vincolo fa sì che esistano diversi scenari da analizzare.
Consideriamo il caso di due autovalori reali e diversi da 1. A causa del vincolo , è impossibile che entrambi abbiano modulo unitario. Uno dei due autovalori (supponiamo ) dovrà essere maggiore di 1 in modulo, mentre l’altro () dovrà essere minore di 1. L’effetto geometrico combinato sulla dinamica del sistema prende il nome di strain (o deformazione a forbice):
Lungo la direzione dell’autovettore , lo stato viene allungato esponenzialmente dal fattore .
Lungo la direzione dell’autovettore , lo stato viene compresso esponenzialmente a ogni passo dal fattore .
L’area totale del parallelogramma nello spazio delle fasi si conserva (poiché la compressione bilancia esattamente l’allungamento), ma la forma si allunga indefinitamente come una striscia infinitamente sottile e lunga. Fisicamente, il sistema diverge ed “esplode”. In questo regime, l’algoritmo è numericamente instabile, in maniera del tutto simile al caso .
Questa divergenza catastrofica viene evitata quando gli autovalori non sono reali ma complessi e coniugati: . In questo caso, il vincolo del determinante si può scrivere come
cioè il modulo di entrambi gli autovalori deve essere esattamente pari a 1. Possiamo quindi scrivere gli autovalori in forma polare come , che mostra esplicitamente come l’evoluzione temporale rappresenti una pura rotazione periodica nel piano complesso. Le traiettorie rimangono limitate e la simulazione è numericamente stabile: gli errori di arrotondamento non vengono amplificati, ma si limitano a oscillare insieme al sistema.
Per visualizzare concretamente il legame profondo tra la conservazione dell’area e la stabilità numerica, analizziamo ora l’animazione mostrata in Figura Figura 5, che confronta l’evoluzione di una regione dello spazio delle fasi secondo i metodi di Eulero ed Eulero-Cromer.
Figure 5:L’evoluzione di un volume di spazio delle fasi (che per l’oscillatore armonico è un piano) delimitato da un rettangolo ottenuto con i metodi di Eulero (in rosso) ed Eulero-Cromer (in blu). I parametri della simulazione sono , (quindi ) e .
L’animazione mostra come nelle condizioni di simulazione (cioè per i valori di e utilizzati), l’algoritmo di Eulero mostra un’espansione dell’area dello spazio delle fasi, che invece non si verifica con Eulero-Cromer. Verifichiamo questi comportamenti calcolando esplicitamente determinanti ed autovalori associati all’oscillatore armonico integrato con i due metodi.
3.5.1Eulero¶
La matrice di propagazione per il metodo di Eulero è
da cui possiamo immediatamente ottenere il determinante:
Poiché e , si ha che per qualunque valore di . Essendo il modulo strettamente maggiore di 1, l’errore globale cresce esponenzialmente a ogni passo temporale. Il metodo è quindi incondizionatamente instabile per l’oscillatore armonico; nello spazio delle fasi, la soluzione numerica descrive una spirale che diverge verso l’infinito, accumulando energia artificiale.
Poiché è un numero strettamente positivo, il determinante della matrice è sempre maggiore di 1. Di conseguenza, il metodo di Eulero è incondizionatamente instabile per l’oscillatore armonico: l’ampiezza delle oscillazioni numeriche crescerà artificialmente all’infinito per qualunque scelta di . Nello spazio delle fasi, questo comportamento si manifesta come mostrato in figura Figura 5: la soluzione numerica descrive una spirale che diverge verso l’infinito, accumulando energia artificiale.
Calcoliamo ora gli autovalori di . Risolvendo il polinomio caratteristico si trova , da cui si ottengono i due autovalori complessi coniugati
Poiché sono complessi coniugati, i due autovalori hanno lo stesso modulo, che vale[11]
cioè un numero maggiore di 1, indipendentemente dal passo di integrazione. Come abbiamo dimostrato precedentemente, se il modulo degli autovalori è strettamente maggiore di 1, l’errore cresce esponenzialmente, dimostrando ancora una volta l’instabilità del metodo di Eulero.
3.5.2Eulero-Cromer¶
Nel caso di Eulero-Cromer, la matrice di propagazione del metodo nello spazio delle fasi è
che ha determinante
Poiché , il metodo conserva l’area nello spazio delle fasi, che per sistemi unidimensionali come l’oscillatore armonico implica simpletticità. Questo garantisce l’assenza di derive energetiche artificiali a lungo termine. In questo caso, il polinomio caratteristico è
da cui si ottengono gli autovalori
Il comportamento del sistema dipende dal segno del radicando ():
. Il radicando è negativo, producendo autovalori complessi coniugati. Poiché il determinante è unitario, e in forza all’equazione (68), i due autovalori devono avere anche modulo 1, e quindi trovarsi sulla circonferenza unitaria. In questo regime il metodo è stabile e genera orbite ellittiche chiuse nello spazio delle fasi.
. Il radicando è positivo, quindi i due autovalori sono reali e distinti. Poiché il loro prodotto deve rimanere pari a 1, uno dei due autovalori sarà necessariamente maggiore di 1 in modulo: il sistema diventa instabile e l’errore diverge esponenzialmente. Questa dipendenza della stabilità dai parametri del sistema (e dell’integrazione numerica) fa sì che il metodo di Eulero-Cromer sia condizionatamente stabile. Nel caso in esame, la condizione di stabilità matematica, che richiede che gli autovalori abbiano modulo 1, impone infatti un limite superiore rigoroso al passo temporale:
3.6Accuratezza¶
Per valutare la bontà (e quindi l’accuratezza) di un metodo di integrazione numerica è fondamentale distinguere tra due definizioni di errore:
Errore di Troncamento Locale (LTE): rappresenta l’errore introdotto dal metodo in un singolo passo temporale , assumendo che tutti i dati al passo precedente siano esatti. Si esprime matematicamente come la differenza tra la soluzione esatta del sistema continuo e quella fornita dallo schema numerico dopo un passo.
Errore Globale: rappresenta l’errore totale accumulato dall’inizio della simulazione fino al tempo finale . Se l’errore locale è dell’ordine di , su un intervallo di tempo limitato (che richiede un numero di passi pari a ) l’errore globale scala come . L’esponente definisce l’ordine di accuratezza del metodo.
3.6.1Eulero¶
La derivazione dell’accuratezza per il metodo di Eulero discende direttamente dallo sviluppo in serie di Taylor di posizione e velocità attorno all’istante :
Confrontando queste espressioni con le equazioni di aggiornamento dello schema di Eulero, eq. (40), si nota immediatamente che lo schema numerico recide i termini di Taylor a partire dal secondo ordine. Definendo l’errore come la differenza tra il dato locale esatto e quello ottenuto numericamente, l’errore locale di troncamento risulta:
Poiché l’errore locale è , l’accumulo globale su passi produce un errore complessivo di ordine . Eulero Esplicito è pertanto un metodo del primo ordine.
3.6.2Eulero-Cromer¶
Nel caso di Eulero-Cromer, lo schema definito in eq. (42) fa uso della velocità aggiornata al tempo successivo per calcolare la nuova posizione. Mentre la relazione per l’aggiornamento di è identica a quella di Eulero Esplicito, e di conseguenza preserva un errore locale pari a , l’analisi della posizione richiede cautela. Sostituendo nella definizione di , possiamo scrivere l’espressione per la variabile in funzione delle sole quantità al tempo :
Confrontiamo ora questa equazione dello schema con lo sviluppo esatto di Taylor di ricavato in precedenza. Calcolando la differenza, si ottiene l’errore di troncamento locale sulla posizione:
Imponendo l’esattezza dei dati al passo , i termini di ordine zero e primo si cancellano, lasciando la discrepanza unicamente sul coefficiente del secondo ordine:
Poiché l’errore locale di troncamento è pari a sia per la velocità che per la posizione, l’integrazione accumula un errore globale proporzionale a , esattamente come per il metodo di Eulero. Quindi, nonostante l’utilizzo di informazioni temporalmente più avanzate per la coordinata spaziale, il metodo di Eulero-Cromer rimane un metodo del primo ordine. L’errore locale sulla posizione ha lo stesso modulo di quello di Eulero Esplicito, ma segno opposto.
Figure 6:La differenza tra la posizione finale teorica e quella ottenuta tramite i due algoritmi di Eulero ed Eulero-Cromer. Il tempo totale di simulazione è , mentre i parametri utilizzati sono , (e quindi ), , .
Se la soluzione teorica è nota (come in questo caso), l’errore si può anche calcolare direttamente dai risultati numerici. Possiamo infatti definire l’errore globale come
dove è il tempo finale, mentre e sono le posizioni ottenute numericamente e teoricamente. La Figura 6 mostra per sistemi con , e , simulati per (in unità adimensionali) con Eulero ed Eulero-Cromer. La figura mostra come l’errore che commette Eulero-Cromer sia, in questo caso, di più di un ordine di grandezza minore rispetto a quello che si ottiene con Eulero, ma la dipendenza da è la stessa. Per grandi valori di l’errore di Eulero sembra crescere più che linearmente, per via dell’instabilità del metodo. Questo comportamento non si verifica con Eulero-Cromer dato che tutti i valori di considerati rispecchiano la relazione di stabilità condizionata, eq. (78).
4Velocity Verlet¶
Introduciamo ora uno dei metodi simplettici, cioè che conserva l’energia e il volume nello spazio delle fasi, migliori e più utilizzati. Sviluppiamo la posizione in serie di Taylor attorno a :
Sommando i due sviluppi notiamo che, per simmetria, i termini con potenze dispari di si elidono e si ottiene
Se trascuriamo i termini di ordine superiore e discretizziamo il tempo, , l’aggiornamento della posizione diventa
Questo è l’algoritmo di Verlet, che permette di calcolare la posizione al passo temporale successivo utilizzando la posizione corrente , la posizione precedente e l’accelerazione corrente . Sebbene scritto in questo modo il metodo di Verlet non coinvolge esplicitamente la velocità, se consideriamo lo sviluppo fino al secondo ordine possiamo scrivere esplicitamente
dove è importante la differenza di accuratezza ( vs ) rispetto a . Inoltre, per sufficientemente piccolo, e saranno molto simili, quindi la loro differenza può dare problemi quando i numeri vengono rappresentanti sul calcolatore. Per questo motivo il metodo è raramente usato in questa forma, ma si modifica per includere un aggiornamento esplicito per la velocità, ottenendo il ben più comune metodo “Velocity Verlet”. Invece di basarsi sulle posizioni dei passi temporali precedente e corrente, l’algoritmo Velocity Verlet aggiorna la posizione e la velocità in un processo a due fasi.
In primo luogo, utilizziamo la velocità e l’accelerazione correnti per aggiornare la posizione al tempo . Questo viene fatto in modo simile al metodo Verlet di base, ma con il termine della velocità esplicitamente incluso:
Questa equazione utilizza la posizione corrente , la velocità corrente e l’accelerazione corrente per calcolare la nuova posizione . Successivamente, dopo aver aggiornato la posizione, dobbiamo calcolare la nuova accelerazione al tempo perché la forza (e quindi l’accelerazione) è cambiata a causa della posizione aggiornata. La nuova accelerazione è data da:
Avendo a disposizione questa nuova accelerazione, possiamo aggiornare la velocità. Invece di usare solo l’accelerazione corrente, il metodo Velocity Verlet usa la media delle accelerazioni corrente e nuova per aggiornare la velocità:
Questa equazione di aggiornamento della velocità tiene conto della variazione dell’accelerazione nell’intervallo di tempo, fornendo un aggiornamento della velocità più accurato rispetto al semplice utilizzo dell’accelerazione corrente. Il metodo Velocity Verlet è lo standard de facto per i codici di Dinamica Molecolare (MD), una tecnica utilizzata per studiare la dinamica e la termodinamica di atomi, molecole, colloidi, ecc. Il modo comune per implementarlo consiste nel suddividere la fase di integrazione della velocità in due, in modo che un passo di integrazione completo diventi:
Aggiornamento della velocità, prima fase: .
Aggiornamento della posizione: (cioè l’eq. (89)).
Calcolo della forza (e quindi dell’accelerazione) utilizzando la nuova posizione: .
Aggiornamento della velocità, seconda fase: (cioè l’eq. (91)).
4.1Stabilità ed accuratezza¶
Consideriamo le equazioni di aggiornamento del metodo Velocity Verlet:
Nel caso dell’oscillatore armonico, e , quindi
Sostituendo la prima equazione nella seconda otteniamo
che, raccogliendo i termini associati a e , diventa:
La matrice di propagazione è quindi
Calcoliamo il determinante della matrice:
Quindi, il determinante è esattamente pari a 1, indipendentemente dal valore del passo temporale : Velocity Verlet è un algoritmo simplettico (conserva l’area nello spazio delle fasi) per qualsiasi parametro di discretizzazione scelto. Non introduce alcuna dissipazione o amplificazione artificiale dell’area di stati iniziali.
Studiato ora gli autovalori. Il polinomio caratteristico di Velocity Verlet è:
Questa equazione vi ricorda qualcosa? Questa equazione è assolutamente identica a quella ottenuta con il metodo di Eulero-Cromer, eq. (76)!
Sebbene le due matrici di evoluzione siano diverse (Velocity Verlet ha coefficienti simmetrici sulla diagonale ed è un metodo del secondo ordine, mentre Eulero-Cromer è asimmetrico ed è del primo ordine), esse condividono lo stesso polinomio caratteristico e di conseguenza hanno gli stessi identici autovalori e quindi la stessa condizione di stabilità, eq. (78). La differenza maggiore tra i due metodi è nella loro accuratezza. Utilizzando la stessa logica applicata a Eulero ed Eulero-Cromer, definiamo l’errore di troncamento locale come la differenza tra la soluzione numerica e quella teorica dopo un passo di integrazione. Sviluppando posizione e velocità fino al quarto ordine, scriviamo le soluzioni esatte come
Se ora sottriamo queste quantità da quelle numeriche, eq. (92), e sostituiamo a il suo sviluppo di Taylor, [12], otteniamo
e quindi sia per la posizione che per la velocità, l’errore locale dell’algoritmo di Verlet (e quindi, equivalentemente, quello di Velocity Verlet) è di ordine , che implica come l’errore globale scali come .
Perché possiamo sviluppare l’accelerazione numerica?
Potrebbe sorgere un legittimo dubbio teorico: l’accelerazione futura è calcolata dall’algoritmo e quindi, per definizione, valutata sulla posizione numerica approssimata () e non sulla posizione reale lungo la traiettoria fisica (). Com’è possibile allora sviluppare in serie di Taylor nel tempo come se ci trovassimo sulla traiettoria esatta?
La giustificazione formale risiede nell’ordine dell’errore locale spaziale. Abbiamo appena dimostrato che l’errore sulla posizione al passo è di terzo ordine:
Se effettuiamo uno sviluppo spaziale in serie di Taylor della funzione continua attorno al punto esatto , otteniamo:
Poiché la differenza tra la coordinata numerica e quella reale è già di ordine , l’errore derivante dal non valutare l’accelerazione sulla traiettoria esatta è di ordine superiore e finisce interamente nel termine di errore generico:
Di conseguenza, per ricavare i termini d’errore fino all’ordine necessari alla nostra dimostrazione, è perfettamente lecito sostituire ad lo sviluppo temporale esatto dell’accelerazione fisica:
L’approssimazione numerica spaziale non altera i coefficienti dei termini di ordine inferiore dello sviluppo.
La Figura 7 illustra vividamente l’enorme impatto del passaggio da un errore globale di ordine a uno di ordine . Per apprezzare concretamente questa differenza, si consideri un passo temporale tipico delle simulazioni reali, ad esempio : in questo scenario, l’accuratezza di Velocity Verlet supera quella di Eulero-Cromer di ben tre ordini di grandezza, riducendo drasticamente l’errore sistematico accumulato sulla traiettoria.
4.2C: Puntatori a funzione, ovvero come scegliere l’algoritmo durante l’esecuzione¶
Finora abbiamo confrontato diversi algoritmi di integrazione, come Eulero, Eulero-Cromer e Velocity Verlet. Un modo poco pratico di farlo consiste nello scrivere un programma diverso per ogni algoritmo, oppure nel disseminare il ciclo di simulazione di istruzioni if. Questo duplica il codice che inizializza il sistema, esegue l’evoluzione e salva i risultati.
Sarebbe più comodo scrivere un solo ciclo di simulazione e scegliere l’algoritmo da usare all’inizio dell’esecuzione. In C possiamo farlo mediante un puntatore a funzione.
Consideriamo tre funzioni che eseguono un singolo passo temporale:
void eulero(Sistema *sist, Punto *p);
void eulero_cromer(Sistema *sist, Punto *p);
void velocity_verlet(Sistema *sist, Punto *p);Le tre funzioni hanno la stessa firma: ricevono argomenti dello stesso tipo e restituiscono tutte void. Possiamo quindi definire un tipo che rappresenta un puntatore a una qualunque funzione con questa firma:
typedef void (*integratore)(Sistema *sist, Punto *p);La dichiarazione può essere letta dall’interno verso l’esterno: integratore è un puntatore (*) a una funzione che riceve un puntatore a un Sistema e uno a un Punto e non restituisce alcun valore.
Possiamo ora dichiarare una variabile di questo tipo e farla puntare all’algoritmo scelto, per esempio utilizzando un argomento da riga di comando
integratore passo;
// vogliamo che il primo argomento determini l'algoritmo di integrazione
int scelta = atoi(argv[1]);
if(scelta == 1) {
passo = eulero;
}
else if(scelta == 2) {
passo = eulero_cromer;
}
else if(scelta == 3) {
passo = velocity_verlet;
}
else {
fprintf(stderr, "Scelta non valida, i valori supportati sono\n");
fprintf(stderr, "1: Eulero\n2: Eulero-Cromer\n3: Velocity Verlet\n");
exit(1);
}Il nome di una funzione, quando viene utilizzato senza parentesi, rappresenta il suo indirizzo. Scrivere
passo = eulero;fa quindi puntare passo alla funzione eulero, ma non esegue ancora la funzione. Per eseguirla usiamo invece le parentesi, esattamente come per una normale chiamata:
passo(&sist, &punto);Il ciclo che realizza la simulazione non deve più sapere quale algoritmo sia stato scelto:
double t = 0.0;
Sistema sistema;
Punto punto;
/* inizializzazione */
while(t < t_max) {
printf("%g %g %g\n", t, punto.x, punto.v);
passo(&sistema, &punto);
t += sistema.dt;
}A ogni iterazione passo chiamerà eulero, eulero_cromer oppure velocity_verlet, a seconda della scelta iniziale. La decisione viene quindi presa a runtime, mentre il programma è in esecuzione.
Il vantaggio non consiste soltanto nel risparmiare qualche riga. Tutto ciò che non dipende dall’algoritmo, come condizioni iniziali, ciclo temporale, scrittura
dei risultati e calcolo dell’errore, rimane in un unico punto del programma. Per aggiungere un nuovo integratore sarà sufficiente scrivere una funzione con la stessa firma e assegnarne l’indirizzo a passo.
Un puntatore a funzione svolge dunque, per le funzioni, un ruolo simile a quello che un normale puntatore svolge per i dati: permette di memorizzare e passare ad altre parti del programma non un valore su cui operare, ma l’operazione da eseguire.
5Il moto del pendolo semplice¶
Torniamo ora al pendolo semplice introdotto nell’eq. (5). A differenza dell’oscillatore armonico, il pendolo è un sistema non lineare e rappresenta quindi un banco di prova più realistico per gli algoritmi di integrazione numerica discussi finora. Introducendo la velocità angolare , l’equazione del moto può essere riscritta come il sistema di due equazioni del primo ordine
In assenza di attrito l’energia meccanica
è una costante del moto. Possiamo quindi utilizzare la conservazione dell’energia come ulteriore strumento per confrontare la qualità dei diversi algoritmi.
Figure 8:Dall’alto verso il basso, i tre pannelli mostrano la posizione , la velocità e l’energia meccanica in funzione del tempo per il pendolo semplice simulato utilizzando i metodi di Eulero, Eulero-Cromer e Velocity Verlet con parametri m/s, m, Kg, rad, rad/s e passo di integrazione s.
La Figura 8 mostra una differenza qualitativa già osservata nel caso dell’oscillatore armonico. Il metodo di Eulero introduce una deriva sistematica: come abbiamo visto, a ogni oscillazione il pendolo acquista una piccola quantità di energia numerica e la traiettoria si allontana progressivamente da quella fisica.
Diversamente, ed esattamente come per l’oscillatore armonico, per Eulero-Cromer e Velocity Verletl’energia numerica non è esattamente costante, ma compie piccole oscillazioni attorno al valore corretto senza mostrare una deriva sistematica (nella figura questo si vede solo per Eulero-Cromer). Come già discusso, questa proprietà è legata alla simpletticità dei due algoritmi: la discretizzazione modifica leggermente la dinamica del sistema, ma ne preserva la struttura geometrica nello spazio delle fasi. In effetti, per un integratore simplettico di ordine ci aspettiamo genericamente un andamento del tipo[13]
dove è una funzione limitata. L’ampiezza delle oscillazioni dell’energia deve quindi scalare come
Possiamo verificare numericamente questa previsione. Una possibile misura dell’errore energetico è il massimo scostamento dal valore iniziale,
In alternativa possiamo considerare le fluttuazioni dell’energia attorno al suo valor medio,
Le due quantità non sono identiche: è sensibile anche a un eventuale spostamento del valor medio rispetto a , mentre misura soltanto le fluttuazioni attorno alla media. Tuttavia, quando l’errore energetico è limitato e oscillatorio e la forma delle oscillazioni non cambia qualitativamente al variare di , entrambe presentano lo stesso comportamento asintotico,
Figure 9:Massimo scostamento dell’energia dal valore iniziale, (in alto) e fluttuazioni dell’energia attorno al valor medio, (in basso), in funzione del passo temporale per il pendolo semplice. Le due rette continue indicano le dipendenze e . I parametri delle simulazioni sono gli stessi della Figura 8.
I risultati della Figura 9 sono coerenti con quanto ci aspettiamo dall’ordine dei metodi. Eulero-Cromer è un algoritmo simplettico del primo ordine e le oscillazioni dell’energia scalano quindi come . Velocity Verlet è invece simplettico e del secondo ordine, e infatti le oscillazioni dell’energia scalano come .
Il metodo di Eulero mostra anch’esso, a tempo finale fissato e nel regime di piccoli , un errore energetico che diminuisce linearmente con . Questo fatto, però, non implica che Eulero sia simplettico: la differenza fondamentale è il comportamento temporale dell’errore. Nel caso di Eulero l’errore tende ad accumularsi producendo una deriva dell’energia, mentre per Eulero-Cromer e Velocity Verlet esso rimane limitato e oscillatorio. La sola dipendenza di da non è quindi sufficiente a stabilire se un algoritmo sia simplettico. La caratteristica distintiva dei metodi simplettici è la combinazione delle due proprietà: l’ampiezza dell’errore energetico diminuisce con l’ordine previsto e, soprattutto, non presenta una deriva secolare a tempi lunghi.
Ordine del metodo e ordine dell’errore sull’energia
È importante non identificare in generale l’ordine di un metodo con l’esponente osservato misurando una singola quantità fisica come l’energia.
Un metodo di ordine garantisce che, a tempo fisico fissato, l’errore globale sullo stato del sistema (posizione e velocità) sia dell’ordine di . Una particolare osservabile può però essere meno sensibile al termine dominante dell’errore e mostrare una convergenza apparentemente più rapida.
Un esempio si incontra con il metodo Runge-Kutta del secondo ordine, che introdurremo più avanti. RK2 è un metodo del secondo ordine e non è simplettico; tuttavia, per l’oscillatore armonico l’errore dominante è principalmente un errore di fase, che non modifica l’energia. Il primo contributo alla deriva energetica compare quindi a un ordine superiore e, su un intervallo temporale fissato, si trova
Questo comportamento può rimanere visibile anche per il pendolo in determinati regimi, ma non rappresenta una proprietà generale di tutti i sistemi non lineari.
Analogamente, il metodo Runge-Kutta del quarto ordine ha un errore globale ma, per motivi simili, per il pendolo si trova .
Il messaggio importante è quindi che l’ordine del metodo non coincide necessariamente con l’oordine dell’errore di un particolare osservabile, mentre per un integratore simplettico di ordine l’ampiezza delle oscillazioni energetiche è solitamente dell’ordine di .
6Runge-Kutta¶
Molti dei problemi complessi da risolvere con metodi numerici non riguardano sistemi in cui l’energia si conserva. La simpletticità non è quindi sempre un requisito necessario. Vediamo subito un esempio.
6.1L’oscillatore armonico smorzato¶
Un oggetto che si muove lentamente in un fluido viscoso è sottoposto, sotto opportune condizioni, a una forza di attrito proporzionale e opposta alla sua velocità[14]. Nel caso di un oscillatore armonico, la dinamica del sistema è descritta dalla seguente equazione differenziale:
dove è il coefficiente di attrito, che determina l’intensità della dissipazione di energia. Introducendo la quantità
si possono distinguere tre diversi regimi dinamici, a seconda del segno di .
: smorzamento sottocritico. La soluzione può essere scritta nella forma
dove e sono costanti determinate dalle condizioni iniziali. Il sistema oscilla con pulsazione , mentre l’ampiezza delle oscillazioni si riduce esponenzialmente nel tempo.
: smorzamento critico. La soluzione ha la forma
Il sistema non oscilla e ritorna all’equilibrio nel minor tempo possibile senza oltrepassarlo.
: smorzamento sovracritico. Ponendo
la soluzione può essere scritta come una combinazione di due esponenziali decrescenti e non presenta oscillazioni:
Per , in tutti e tre i regimi si ha e per . L’energia meccanica
decresce infatti secondo
e tende a zero a tempi lunghi.
Figure 10:Come in Figura 7 per per l’oscillatore armonico smorzato. Il tempo totale di simulazione è , mentre i parametri utilizzati sono , (e quindi ), , , .
La Figura 10 mostra l’errore globale che i tre algoritmi che abbiamo studiato compiono nel caso di un oscillatore armonico smorzato. Qualitativamente, i risultati sono simili a quelli ottenuti per l’oscillatore armonico non smorzato e mostrati nella Figura 7, con una differenza importante: il metodo di Velocity Verlet, pur risultando quantitativamente più accurato dei metodi di Eulero e di Eulero-Cromer, presenta in questo caso un errore globale che scala come .
Questa perdita di accuratezza è dovuta al fatto che la formulazione standard del Velocity Verlet è costruita per sistemi nei quali l’accelerazione dipende dalla posizione, ma non dalla velocità. Nell’oscillatore smorzato, invece,
e il calcolo dell’accelerazione al passo successivo richiede quindi anche una stima della nuova velocità. Se il metodo viene applicato senza modificarne la struttura per trattare esplicitamente questa dipendenza, l’accelerazione viene valutata utilizzando una velocità non ancora aggiornata in modo pienamente consistente. L’errore introdotto da questa approssimazione è di ordine per ogni singolo passo e si accumula nel corso dell’integrazione, producendo un errore globale di ordine . Il Velocity Verlet standard perde pertanto, in presenza di forze dipendenti dalla velocità, la convergenza del secondo ordine che possiede per sistemi conservativi con accelerazione dipendente dalla sola posizione.
Per ovviare a questo problema si può estendere il metodo Velocity Verlet al caso in cui la forza dipenda esplicitamente dalla velocità, oppure utilizzare un metodo non simplettico che non soffre di questi problemi. In questa sezione presenteremo il più famoso di questi metodi, noto come Runge-Kutta, discutendo la sua implementazione al secondo (RK2) e al quarto (RK4) ordine. Per rendere la trattazione della stabilità e dell’accuratezza direttamente confrontabile con quella dei capitoli precedenti, presenteremo inizialmente i nuovi metodi utilizzando l’oscillatore armonico non smorzato, eq. (29). L’oscillatore armonico smorzato sarà invece il sistema modello che utilizzeremo alla fine della sezione per confrontare gli algoritmi già esaminati con i due RK che introdurremo qui.
6.2Metodo Runge-Kutta del secondo ordine (RK2)¶
I metodi introdotti finora approssimano il valore medio della derivata nell’intervallo utilizzando informazioni disponibili agli estremi dell’intervallo stesso. Possiamo però ottenere una stima più accurata osservando che, per una funzione sufficientemente regolare, il valore della derivata nel punto medio dell’intervallo costituisce spesso una buona approssimazione del suo valore medio.
L’idea alla base del metodo Runge-Kutta del secondo ordine consiste quindi nello stimare la derivata nel punto medio dell’intervallo e utilizzarla per aggiornare la soluzione. Partendo dal valore noto , si esegue innanzitutto un mezzo passo con il metodo di Eulero:
Questa quantità fornisce una stima della soluzione al tempo . Possiamo quindi calcolare una nuova stima della derivata nel punto medio:
Infine, utilizziamo questa derivata per avanzare di un passo completo:
Ricordando che, nel nostro caso, , applicando il metodo al nostro sistema di equazioni otteniamo
dove è l’accelerazione calcolata nel punto medio dell’intervallo, ottenuta integrando di mezzo passo e . Calcolare due volte l’accelerazione (in e in ) è il prezzo computazionale che si paga per migliorare l’accuratezza rispetto al metodo di Eulero. Questo è un prezzo che molte volte (ma non necessariamente sempre) vale la pena di pagare.
6.2.1Stabilità ed accuratezza¶
Considerando che, per l’oscillatore armonico, , le equazioni di aggiornamento (125) si possono scrivere come
che permette di scrivere la matrice di propagazione per il metodo di RK2:
il cui determinante vale
Poiché il determinante della matrice è sempre maggiore di uno, il metodo non conserva l’energia e quindi non può essere simplettico, e le traiettorie generate nello spazio delle fasi spiraleggiano verso l’esterno. Come fatto in precedenza, utilizziamo il polinomio caratteristico per calcolare gli autovalori, che valgono
I due autovalori sono sempre complessi coniugati e quindi hanno lo stesso modulo, che vale : il metodo è, come quello di Eulero, incondizionatamente instabile, ma in questo caso la differenza tra il modulo al quadrato degli autovalori e 1 è più piccola ( vs. )[15]. Di conseguenza, l’instabilità diventa evidente dopo tempi di integrazione molto più lunghi rispetto ad Eulero.
Per quanto riguarda l’accuratezza, notiamo prima di tutto che l’aggiornamento delle posizioni, cioè la prima delle equazioni (125), è lo stesso del metodo di Velocity Verlet, eq. (92). Di conseguenza, i due algoritmi condividono l’errore di troncamento locale per le posizioni, che va come . Per le velocità applichiamo lo stesso procedimento visto per Velocity Verlet espandendo per ottenere
Se ora sottraiamo questa quantità da quella teorica, eq. (99), e assumiamo come al solito che al tempo lo stato numerico coincida con quello esatto (, , ), otteniamo:
e quindi anche per la velocità l’errore locale dell’algoritmo RK2 è di ordine .
Quindi, sia per le posizioni che per le velocità l’errore globale scala come : RK2 è un algoritmo del secondo ordine nel tempo. La figura che mostra questo andamento è mostrata e discussa più sotto.
6.3Metodo Runge-Kutta del quarto ordine (RK4)¶
Il metodo RK2 migliora l’accuratezza dell’integrazione utilizzando una stima della derivata nel punto medio dell’intervallo. Possiamo però ottenere una stima ancora più accurata del valore medio della derivata combinando informazioni provenienti da più punti dell’intervallo stesso.
L’idea alla base del metodo Runge-Kutta del quarto ordine consiste nel costruire una successione di stime della derivata e combinarle opportunamente per ottenere una migliore approssimazione del valore medio di tra e .
Si definiscono innanzitutto quattro stime della derivata:
che rappresenta la derivata all’inizio dell’intervallo,
che fornisce una prima stima della derivata nel punto medio,
che costituisce una stima migliorata della derivata nel punto medio,
e infine
che rappresenta una stima della derivata alla fine dell’intervallo.
Queste quattro quantità, pesate opportunamente, si possono combinare per ottenere una stima della derivata media:
dove i vanno scelti in modo da minimizzare l’errore. Si può dimostrare (vedi box sotto per una derivazione semplificata) che fissando e il metodo fornisce una stima particolarmente accurata della derivata media nell’intervallo, raggiungendo un’accuratezza molto superiore rispetto a Eulero ed RK2, al costo di quattro valutazioni della funzione per ogni passo temporale.
L’aggiornamento della soluzione assume pertanto la forma
Se siamo interessato a un sistema dinamico unidimensionale come l’oscillatore armonico, il metodo RK4 va applicato simultaneamente alle due variabili e . In altre parole, a ogni passo temporale dobbiamo costruire quattro stime sia per la derivata della posizione, cioè la velocità, sia per la derivata della velocità, cioè l’accelerazione.
Partendo dallo stato noto al tempo , definiamo innanzitutto
Queste sono le derivate valutate all’inizio dell’intervallo. Usiamo poi queste quantità per stimare lo stato del sistema a metà passo:
e calcoliamo le derivate in questo punto intermedio:
Ripetiamo ora la stessa procedura, ma usando per ottenere una stima migliorata dello stato a metà passo:
da cui
Infine, usiamo per stimare lo stato alla fine dell’intervallo:
e calcoliamo l’ultima coppia di derivate:
L’aggiornamento completo si ottiene quindi combinando le quattro stime con gli stessi pesi già ricavati per il caso generale:
Nota Bene: per l’oscillatore armonico l’accelerazione dipende solo dalla velocità,e quindi le quantità si ottengono semplicemente valutando nei diversi punti intermedi costruiti dall’algoritmo. Nel caso più generale (ad esempio quello dell’oscillatore smorzato o forzato), l’accelerazione dipende anche da e da , e quindi è importante aggiornare correttamente entrambe le variabili nei passi intermedi.
6.3.1Stabilità ed accuratezza¶
Come per gli algoritmi visti, studiamo la stabilità di RK4 applicandolo all’oscillatore armonico. Poiché in questo caso i calcoli sono piuttosto lunghi e non fondamentali per quello che ci interessa, li riporto in un box più in basso. Discutiamo invece i risultati principali:
Il determinante della matrice di propagazione è
Gli autovalori invece valgono
Poiché gli autovalori sono complessi coniugati, il loro modulo è uguale al determinante:
Questo risultato mostra che, come RK2, anche RK4 non è simplettico. Tuttavia la deviazione da 1 compare soltanto a partire dall’ordine , cioè è molto piccola per passi temporali sufficientemente piccoli. A differenza di RK2, l’espressione del determinante ha un termine negativo e uno positivo, e quindi può cambiare segno a seconda dei parametri. La condizione per cui la dinamica non esplode è , cioè
che, per , equivale a
da cui otteniamo la condizione di stabilità
Il metodo RK4 è dunque condizionatamente stabile per l’oscillatore armonico. Il limite di stabilità è meno restrittivo di quello trovato per Eulero-Cromer e Velocity Verlet, eq. (78), ma questo non significa che RK4 sia sempre preferibile. Infatti, RK4 non è simplettico: anche quando è stabile, non conserva esattamente l’area nello spazio delle fasi. Per il determinante è leggermente minore di 1, quindi lo schema introduce una piccola dissipazione numerica: le orbite nello spazio delle fasi tendono a spiraleggiare molto lentamente verso l’interno. Per , invece, il determinante (e quindi il modulo di entrambi gli autovalori) diventa maggiore di 1 e la soluzione numerica diverge.
Discutiamo ora l’accuratezza dell’algoritmo. Come viene dimostrato formalmente nel box sotto, ma si può anche inferire notando che i coefficienti numerici sono scelti in maniera da eguagliare al quarto ordine lo sviluppo di Taylor della soluzione analitica, l’errore di troncamento locale è . L’errore globale accumulato scala quindi come
Figure 11:Come in Figura 6, con, in aggiunta, l’errore ottenuto applicando gli algoritmi Runge-Kutta del secondo e quarto ordine.
La figura Figura 11 mostra l’andamento degli errori globali di RK2 ed RK4 insieme a quelli degli altri algoritmi visti finora in questa sezione. È evidente che l’errore di RK4 diminuisce molto più rapidamente di quello di RK2 o, equivalentemente, Velocity Verlet, quando si riduce il passo temporale. Il grafico mostra addirittura una saturazione per piccoli valori di : per questi valori del passo di integrazione, l’errore dovuto alla precisione numerica delle operazioni in virgola mobile comincia a dominare l’errore totale, e l’errore di troncamento diventa trascurabile.
Riassumendo, RK4 è molto accurato su intervalli di tempo finiti, ma non rispetta esattamente la struttura geometrica dei sistemi conservativi, e la mancata simpletticità può produrre una lenta deriva artificiale dell’energia. Per simulazioni molto lunghe di sistemi conservativi, un metodo simplettico come Velocity Verlet può quindi produrre un comportamento qualitativamente migliore, anche se l’ordine formale di accuratezza è più basso.
6.4Confronto tra algoritmi: l’oscillatore armonico smorzato¶
Figure 12:Come in Figura 11 per l’oscillatore armonico smorzato. Il tempo totale di simulazione è , mentre i parametri utilizzati sono , (e quindi ), , , .
La Figura 12 mostra l’andamento degli errori globali ottenuti risolvendo numericamente l’equazione dell’oscillatore armonico smorzato, eq. (113), con i diversi metodi di integrazione introdotti in questa sezione. Si vede bene il comportamento di scala dell’errore globale degli algoritmi RK2 e RK4 rimanga invariato anche in sistemi in cui, come in questo caso, la forza dipende esplicitamente dalla velocità. Questo fa sì che RK4 (che è marginalmente più complesso da implementare e, per un calcolatore, da eseguire) sia il metodo più comune per risolvere equazioni (o sistemi di equazioni) differenziali che non richiedono la proprietà di simpletticità.
7Dal moto unidimensionale ai sistemi planetari¶
Finora abbiamo applicato gli algoritmi di integrazione numerica a sistemi con un solo grado di libertà, descritti da una posizione e da una velocità . Gli stessi metodi possono però essere utilizzati per studiare sistemi multidimensionali formati da molti corpi. In questo caso lo stato del sistema non è più rappresentato da una singola coppia , ma dall’insieme delle posizioni e delle velocità di tutti i corpi.
Come esempio consideriamo un sistema formato da una stella di massa , mantenuta fissa nell’origine, e da pianeti vincolati a muoversi sullo stesso piano. Il pianeta -esimo ha massa , posizione
e velocità
L’ipotesi di moto planare è appropriata quando le posizioni e le velocità iniziali appartengono allo stesso piano. Le forze gravitazionali rimangono infatti contenute in quel piano e non possono generare una componente del moto perpendicolare a esso.
Mantenere la stella immobile è invece un’approssimazione: anche la stella dovrebbe accelerare per effetto dell’attrazione dei pianeti. L’approssimazione è ragionevole quando la massa della stella è molto maggiore della somma delle masse planetarie, o in altre parole
così che il suo spostamento risulti molto piccolo rispetto alle dimensioni delle orbite. Ad esempio, per il sistema solare la somma delle masse dei pianeti è kg (con Giove che da solo rappresenta più del 70% della massa planetaria), mentre la massa del Sole vale kg. Quindi, se volessimo simulare il nostro sistema solare, dove la massa planetaria è appena lo 0.134% di quella solare, la relazione (197) sarebbe valida.
7.1Pianeti non interagenti¶
Consideriamo inizialmente pianeti che interagiscono con la stella, ma non tra loro. La forza gravitazionale esercitata dalla stella sul pianeta è
dove è la costante di gravitazione universale e
è la distanza del pianeta dalla stella. Applicando la seconda legge di Newton si ottiene
La massa del pianeta non compare nell’accelerazione: a parità di posizione e velocità iniziali, tutti i corpi seguono la stessa traiettoria indipendentemente dalla loro massa.
In componenti, le equazioni del moto sono
Per ogni pianeta dobbiamo quindi integrare quattro equazioni differenziali del primo ordine. Un sistema di pianeti è descritto complessivamente da variabili dinamiche. Nel caso non interagente, tuttavia, le equazioni relative a pianeti diversi sono indipendenti: stiamo semplicemente risolvendo problemi di Keplero separati.
7.2Orbite circolari¶
Un caso particolarmente semplice è quello di un pianeta in orbita circolare di raggio . In questo caso l’accelerazione gravitazionale deve coincidere con l’accelerazione centripeta:
La velocità necessaria per ottenere un’orbita circolare è quindi
Se inizialmente il pianeta si trova nel punto , una possibile condizione iniziale è pertanto
Il periodo dell’orbita vale
da cui segue la terza legge di Keplero,
Le orbite circolari sono particolarmente utili per verificare un programma: raggio, velocità ed energia devono rimanere costanti, a meno dei piccoli errori introdotti dall’integrazione numerica.
7.3Energia e momento angolare¶
Poiché la forza gravitazionale è conservativa, a ciascun pianeta possiamo associare l’energia meccanica
In assenza di interazioni fra i pianeti, l’energia di ciascun pianeta si conserva separatamente. Si conserva quindi anche l’energia totale
Per un’orbita circolare, sostituendo , si ottiene
L’energia è negativa perché il pianeta è gravitazionalmente legato alla stella. Più in generale, valori negativi dell’energia corrispondono a orbite legate, mentre un corpo con energia non negativa può allontanarsi indefinitamente dalla stella.
Il momento lineare dei pianeti non si conserva, perché la stella fissa esercita su di essi una forza esterna. Per descrivere un sistema completamente isolato sarebbe necessario lasciare libera anche la stella e integrare la sua equazione del moto. D’altro canto, la forza gravitazionale è una forza centrale: è sempre parallela a e non esercita quindi alcun momento rispetto all’origine. Di conseguenza, il momento angolare del pianeta si conserva. Nel moto planare l’unica componente non nulla è quella perpendicolare al piano:
Nel caso non interagente si conserva separatamente ogni , e di conseguenza anche
Energia e momento angolare forniscono due strumenti fondamentali per valutare la qualità dell’integrazione numerica. Una deriva sistematica di queste quantità può segnalare un passo temporale troppo grande o un algoritmo poco adatto allo studio di sistemi conservativi.
7.4Interazioni fra i pianeti¶
In un sistema planetario reale ogni pianeta è attratto non soltanto dalla stella, ma anche dagli altri pianeti. La forza esercitata dal pianeta sul pianeta è
L’accelerazione complessiva del pianeta diventa quindi
Le equazioni dei diversi pianeti sono ora accoppiate: per conoscere l’accelerazione di un pianeta dobbiamo conoscere simultaneamente le posizioni di tutti gli altri. In generale il problema non può più essere scomposto in orbite kepleriane indipendenti e non possiede una soluzione analitica.
L’interazione può produrre precessioni, scambi di energia e momento angolare, risonanze orbitali e, in opportune condizioni, dinamiche caotiche o espulsioni dal sistema. Anche quando le interazioni sono deboli, i loro effetti possono accumularsi su tempi molto lunghi e modificare sensibilmente le orbite.
L’energia meccanica totale deve ora includere anche l’energia potenziale associata a ogni coppia di pianeti:
La condizione garantisce che ogni coppia venga contata una sola volta. Se sommassimo su tutti gli indici distinti , conteremmo infatti due volte la stessa interazione: una come coppia e una come coppia .
In presenza di interazioni non si conserva più l’energia di ciascun pianeta: i pianeti possono scambiarsi energia. Si conserva invece l’energia totale del sistema. Analogamente, il momento angolare di ogni singolo pianeta può variare, mentre si conserva il momento angolare totale
7.5Integrazione numerica di un sistema di molti corpi¶
Gli algoritmi introdotti in questo capitolo possono essere applicati senza modifiche concettuali. Posizioni e velocità non sono più singoli numeri, ma insiemi di array (o array di strutture). Poiché le forze gravitazionali dipendono dalle posizioni ma non dalle velocità, Velocity Verlet è una scelta naturale: è un metodo del secondo ordine e, essendo simplettico, descrive generalmente meglio l’evoluzione a lungo termine dei sistemi conservativi.
È però essenziale che, a ogni fase dell’algoritmo, le accelerazioni siano calcolate usando posizioni riferite allo stesso istante. Non possiamo aggiornare completamente il primo pianeta e utilizzare subito la sua nuova posizione per calcolare l’accelerazione del secondo: in questo modo il risultato dipenderebbe arbitrariamente dall’ordine con cui i pianeti sono memorizzati.
Un passo di Velocity Verlet deve quindi essere organizzato collettivamente:
si calcolano le accelerazioni di tutti i pianeti;
si aggiornano tutte le posizioni;
si calcolano le nuove accelerazioni usando tutte le posizioni aggiornate;
si aggiornano tutte le velocità.
Nel caso non interagente, il calcolo delle accelerazioni richiede un numero di operazioni proporzionale a . Includendo le interazioni, per ciascuno degli pianeti dobbiamo sommare il contributo degli altri : il costo di un passo cresce quindi come . Questa differenza diventa fondamentale nelle simulazioni con un numero molto grande di corpi.
Rimane infine un problema pratico. Finora il numero delle variabili era stabilito direttamente nel codice. In un programma generale per il moto planetario, invece, vogliamo scegliere il numero di pianeti durante l’esecuzione. Dobbiamo pertanto riservare una quantità di memoria che dipenda da un valore non noto al momento della compilazione. Nella prossima sezione vedremo come farlo in C mediante l’allocazione dinamica della memoria.
7.6C: Allocazione dinamica della memoria con malloc e free¶
Nel programma appena descritto vogliamo scegliere il numero di pianeti durante l’esecuzione, per esempio prendendolo come argomento dalla riga di comando. Non possiamo quindi dichiarare un array di dimensione fissata nel codice sorgente:
Pianeta pianeti[2];Una dichiarazione di questo tipo permette di simulare esattamente due pianeti. Se volessimo modificarne il numero dovremmo cambiare il codice e compilarlo nuovamente.
Possiamo invece riservare la memoria necessaria durante l’esecuzione mediante la funzione malloc, dichiarata nell’header stdlib.h. Supponiamo che la struttura che rappresenta un pianeta sia stata definita nel modo seguente:
typedef struct {
double m;
double x, y;
double vx, vy;
} Pianeta;Dopo aver determinato il numero di pianeti, possiamo dichiarare un puntatore e allocare dinamicamente un array:
int N = atoi(argv[1]); /* primo argomento da riga di comando */
Pianeta *pianeti = malloc(N * sizeof(Pianeta));La funzione malloc riceve come argomento il numero di byte da allocare. L’espressione sizeof(Pianeta) (o, equivalententemente, sizeof *pianeti), restituisce il numero di byte necessari per memorizzare un oggetto di tipo Pianeta. Il prodotto
N * sizeof(Pianeta)è quindi la quantità di memoria necessaria per contenere pianeti, espressa in byte[16].
Il valore restituito da malloc è l’indirizzo iniziale della regione di memoria allocata. Possiamo quindi utilizzare pianeti come un normale array:
pianeti[0].m = 1.0e-3;
pianeti[0].x = 1.0;
pianeti[0].y = 0.0;
pianeti[0].vx = 0.0;
pianeti[0].vy = 1.0;o, più in generale,
for (int i = 0; i < N; i++) {
/* inizializzazione del pianeta i-esimo */
}In generale, gli array dinamici sottostanno alle stesse regole degli array normali:
la memoria è contigua:
pianeti[0],pianeti[1], ,pianeti[N - 1]identificano strutturePianetaconsecutive.accedere a zone di memoria non allocate (cioè che non fanno parte dei byte richiesti tramite
malloc) dà luogo a undefined behaviour che, nel migliore dei casi farà crashare il programma, mentre potrebbe più insidiosamente dare luogo a comportamenti non riproducibili del codice.vale l’aritmetica dei puntatori:
pianeti[i + 1]e*(pianeti + i + 1)si riferiscono allo stesso elemento.
Notate inoltre che in C non è necessario convertire esplicitamente il valore restituito da malloc. È quindi possibile scrivere sia
Pianeta *pianeti = malloc(N * sizeof(Pianeta));che
Pianeta *pianeti = (Pianeta *) malloc(N * sizeof(Pianeta));La conversione da void *, il tipo restituito da malloc, a un altro tipo di puntatore avviene infatti automaticamente.
7.6.1Inizializzazione della memoria e calloc¶
malloc riserva la memoria, ma non ne inizializza il contenuto. Immediatamente dopo l’allocazione, i campi delle strutture contengono valori indeterminati:
Pianeta *pianeti = malloc(N * sizeof(Pianeta));
/* pianeti[0].x non possiede ancora un valore utilizzabile */Prima di leggere un campo dobbiamo quindi assegnargli esplicitamente un valore. Nel programma sui pianeti questo avverrà naturalmente leggendo le masse e le condizioni iniziali dalla riga di comando o da un file. Se vogliamo invece che la memoria venga inizialmente azzerata, possiamo usare calloc al posto di malloc:
Pianeta *pianeti = calloc(N, sizeof(Pianeta));A differenza di malloc, calloc riceve separatamente il numero di elementi e la dimensione di ciascun elemento, e inizializza a zero tutti i byte della memoria allocata. Per esempio, può essere comodo inizializzare a zero due array destinati a contenere le componenti delle accelerazioni:
double *ax = calloc(N, sizeof(double));
double *ay = calloc(N, sizeof(double));Anche il risultato di calloc deve essere controllato:
if(ax == NULL || ay == NULL) {
fprintf(stderr, "Impossibile allocare la memoria\n");
exit(1);
}D’altro canto, se dobbiamo comunque assegnare immediatamente un valore a tutti gli elementi, come nel caso delle condizioni iniziali dei pianeti, tanto vale usare malloc.
7.6.2Liberare la memoria con free¶
La memoria allocata dinamicamente rimane riservata finché non viene liberata esplicitamente. Quando un array non serve più, dobbiamo restituire la memoria al sistema mediante free:
free(pianeti);Nel nostro programma, alla fine della simulazione scriveremo per esempio:
free(ax);
free(ay);
free(pianeti);
return 0;Ogni regione ottenuta mediante malloc o calloc deve essere liberata esattamente una volta. Se perdiamo il puntatore senza aver chiamato free, la memoria rimane occupata ma non è più accessibile: si verifica quello che viene chiamato un memory leak, cioè una perdita di memoria.
Dopo free, il puntatore continua formalmente a contenere il vecchio indirizzo, ma la memoria corrispondente è stata “restituita” al sistema operativo, e quindi non appartiene più al programma. Non possiamo quindi usarlo:
free(pianeti);
/* Errore: la memoria è già stata liberata, probabile segmentation fault */
pianeti[0].x = 1.0;Non bisogna neppure chiamare free due volte sullo stesso puntatore:
free(pianeti);
free(pianeti); /* Errore */7.6.3Il ciclo di vita dell’array¶
L’utilizzo della memoria dinamica segue quindi quattro passaggi:
determiniamo durante l’esecuzione il numero di oggetti da memorizzare;
allochiamo la memoria con
mallococalloc;controlliamo il risultato e utilizziamo l’array;
liberiamo la memoria con
freequando non serve più.
Nel caso dei pianeti, lo schema complessivo assume la forma
int N = atoi(argv[1]);
Pianeta *pianeti = malloc(N * sizeof(Pianeta));
if(pianeti == NULL) {
fprintf(stderr, "Impossibile allocare la memoria\n");
return 1;
}
/* Lettura delle condizioni iniziali */
/* Integrazione delle equazioni del moto */
/* Analisi o scrittura dei risultati */
free(pianeti);
return 0;La memoria dinamica ci permette così di scrivere un solo programma capace di simulare un numero arbitrario di pianeti, limitato soltanto dalla memoria disponibile, senza conoscere quando scriviamo o compiliamo il codice.
E non solo: le equazioni differenziali appaiono in praticamente ogni ambito scientifico, o comunque in cui analisi e modelli quantitativi sono possibili.
È possibile rendere questa definizione, che qui sembra piuttosto generica, formale e non ambigua.
Che deriva dalla formula di bisezione .
Altri esempi di metodi comunemente utilizzati per risolvere sistemi di equazioni differenziali (spesso alle derivate parziali) sono gli elementi finiti (finite element methods) e la risoluzione in spazio di Fourier tramite Fast Fourier transform (FFT).
Se avete un occhio attento potete notare qualche discrepanza tra la posizione teorica e quella ottenuta con in prossimità di massimi e minimi
in effetti questo vale per qualunque sistema dinamico unidimensionale
In spazi delle fasi a più dimensioni (sistemi con gradi di libertà, dove lo spazio delle fasi ha dimensione ), la simpletticità è una condizione molto più restrittiva della semplice conservazione del volume. Un algoritmo simplettico deve conservare non solo il volume totale (), ma anche le proiezioni delle aree orientate su tutte le coppie di piani coordinati coniugati , dove è il momento coniugato a .
con quasi conservazione si intende quella proprietà per cui l’energia meccanica di un sistema fluttua intorno a un valore costante. Quando integriamo numericamente delle equazioni differenziali non possiamo sperare di fare meglio.
Questo risultato si può ottenere immediatamente ricordando che .
Si veda il box più in basso sul perché possiamo farlo.
Questa proprietà deriva dal fatto che, per sufficientemente piccolo, una dinamica simplettica può essere interpretata come quella generata da un sistema con un’Hamiltoniana (non sapete cos’è? Lo saprete presto grazie al corso di Meccanica Analitica) leggermente modificata rispetto a quella originale,
Questa proprietà, che può essere formalizzata tramite la cosiddetta backward error analysis, implica che l’errore sull’energia fisica rimanga generalmente limitato anche per tempi di integrazione molto lunghi.
Il regime in cui la forza di attrito è proporzionale alla velocità si può quantificare introducendo il numero di Reynolds , dove è la densità del fluido, la sua viscosità dinamica, la velocità caratteristica e la dimensione caratteristica dell’oggetto. La legge di attrito lineare è valida per . Questa condizione si realizza tipicamente per oggetti molto piccoli, velocità ridotte o fluidi con elevata viscosità cinematica. A numeri di Reynolds elevati, in molti regimi il contributo dominante alla resistenza del fluido è invece approssimativamente proporzionale al quadrato della velocità.
Se , come si dovrebbe sempre avere.
Avremmo potuto scrivere anche
malloc(N * sizeof *pianeti);. In questo casosizeofavrebbe restituito il numero di byte necessari per memorizzare un oggetto del tipo a cui puntapianeti.