Metodi elementari per l’integrazione numericadi equazioni differenziali ordinarie

Marco Sansottera

Introduzione

Consideriamo un’equazione differenziale ordinaria autonoma

x(t) = f (x(t)) , x : I → R ,

con dato iniziale x(0) = x0 .

Analitica vs. NumericaGenericamente, non e possibile determinare una soluzione analitica, masolamente una soluzione approssimata numericamente.

AttenzioneQuanto e affidabile una soluzione numerica?

Quanto bene riesco a descrivere l’evoluzione del sistema?

Approccio numerico

Approssimare il flusso determinato dalla soluzione con una mappa:

si parte dal dato iniziale x0;

si sceglie una opportuna discretizzazione temporale;

si cerca una formula per approssimare la soluzione, e.g., usando lasoluzione e la sua derivata nei punti precedentemente calcolati.

I cosiddeti metodi one-step sono quelli che sfruttano la soluzione e la suaderivata nell’ultimo punto precedentemente calcolato.

L’idea e quella di sfruttare, in modo opportuno, lo sviluppo di Taylor

x(t + τ) = x(t) + τ x(t) +τ 2

2x(t) + . . . ,

o la formulazione integrale

x(t + τ)− x(t) =

∫ τ

0

x(t + h)dh =

∫ τ

0

f (t + h)dh ,

dove l’integrale viene sostituito da una opportuna formula di quadratura.

Metodi elementari one-step

x = f (x)

Metodi elementari one-step

Eulero Eulero implicitoxn+1 = xn + τ f (xn) xn+1 = xn + τ f (xn+1)

Trapezi Punto medioxn+1 = xn + τ

2

(f (xn) + f (xn+1)

)xn+1 = xn + τ f

( xn+xn+1

2

)Eulero ed Eulero implicito: metodi di ordine 1, l’errore locale e O(τ 2).

Trapezi e Punto medio: metodi di ordine 2, l’errore locale e O(τ 3).

Esercizi

Oscillatore armonico — 1Si consideri un oscillatore armonico. Al variare dei dati iniziali:determinare numericamente la soluzione e tenere traccia dell’errorecommesso rispetto alla soluzione analitica e di quello relativo all’energia.Discutere la dinamica osservata.

Lennard-Jones — 1Si consideri il sistema

q = v ,

v = −φ′(q) , φ(q) = q−12 − 2q−6 .

Al variare dei dati iniziali: determinare numericamente la soluzione etenere traccia dell’errore commesso nell’energia. Discutere la dinamicaosservata e confrontarla con uno studio qualitativo “analitico”.

Esercizi

Lennard-Jones — 1 (stima errore)

Stimare l’errore del metodo di Eulero sfruttando lo sviluppo di Taylordella soluzione.E immediato ottenere una stima del tipo

‖en+1‖ ≤ (1 + τa(tn))‖en‖+ τ 2b(tn) +O(τ 2) ,

con

a(tn) = ‖Jf (x(tn))‖ , b(tn) = max[tn,tn+1]

1

2‖(x(τ))‖ .

Nell’esempio di Lennard-Jones, verificare che l’errore in energia aumentabruscamente in corrispondenza di particolari condizioni e discutere talecomportamento.Si consiglia un grafico con i valori di a(tn), b(tn) e dell’errore in energia.

High-order methods: Runge-Kutta

RK2

X1 = xn +τ

2f (xn)

xn+1 = xn + τ f (X1)

RK4

X1 = xn

X2 = xn +τ

2f (X1)

X3 = xn +τ

2f (X2)

X4 = xn + τ f (X3)

xn+1 = xn +τ

6

(f (X1) + f (X2) + f (X3) + f (X4)

)

Leap-frog

x = −∇V (x) = f (x)

{x = v ,

v = f (x) .

LF2

xn+1 = xn + τvn+1/2

vn+3/2 = vn+1/2 + τ f (xn+1)

AttenzionePer inizializzare l’algoritmo devo calcolare v1/2, posso usare un mezzopasso di Eulero o qualche schema piu sofisticato.

Verlet

x = −∇V (x) = f (x)

{x = v ,

v = f (x) .

Verlet (velocita)

vn+1/2 = vn +τ

2f (xn)

xn+1 = xn + τvn+1/2

vn+1 = vn+1/2 +τ

2f (xn+1)

Verlet (posizioni)

xn+1/2 = xn +τ

2vn

vn+1 = vn + τ f (xn+1/2)

xn+1 = xn+1/2 +τ

2vn+1

Esercizi

Oscillatore armonico — 2Si consideri un oscillatore armonico. Al variare dei dati iniziali:determinare numericamente la soluzione e tenere traccia dell’errorecommesso rispetto alla soluzione analitica e di quello relativo all’energia.Discutere le differenze tra i metodi Runge-Kutta e quelli simplettici(Leap-frog e Verlet).

Lennard-Jones — 2Si consideri il sistema

q = v ,

v = −φ′(q) , φ(q) = q−12 − 2q−6 .

Al variare dei dati iniziali: determinare numericamente la soluzione etenere traccia dell’errore commesso nell’energia. Discutere la dinamicaosservata e confrontarla con uno studio qualitativo “analitico”. Discuterele differenze tra i metodi Runge-Kutta e quelli simplettici (Leap-frog eVerlet).

Contopoulos — Henon-Heiles

Un modello di oscillatori interagenti

H =ω0

2(q20 + p20) +

ω1

2(q21 + p21) + q20q1 −

1

3q31

Consideriamo ω0 e ω1 positive:

il potenziale ha un minimo locale in (0, 0);

il potenziale ha tre punti di sella

(0, ω1)

(±√ω20 + 2ω0ω1

2,−ω0

2

)che corrispondono ai valori di energia

E1 =ω31

6E2 =

ω30

24+ω20ω1

8

Contopoulos — Henon-Heiles

Sempre per ω0 e ω1 positive, e interessante restringersi a quei valoridell’energia per cui la dinamica e confinata per tutti i tempi!

Defininamo l’energia di fuga o energia critica

Ec = min

(ω31

6,ω30

24+ω20ω1

8

)

Per 0 < E < Ec la dinamica e confinata.

Per una spiegazione dettagliata del modello di Henon-Heiles, si vedahttp://www.mat.unimi.it/users/antonio/metmod/Note_6.pdf

Esercizi

Sezione Poincare modello Henon-HeilesAl variare dell’energia, studiare la sezione di Poincare in corrispondenzaall’iperpiano q0 = 0 proiettando sul piano q1 p1.

Per i dati iniziali: fissare il valore dell’energia E0, scegliere arbitrariamente(q1, p1) (in modo sensato), porre q0 = 0 (in modo da essere sullasezione!) e determinare p0 > 0 in modo che l’energia iniziale sia E0.

Se 0 < E < Ec , riportare anche la “curva critica” data dall’equazione

ω1

2(q21 + p21)− 1

3q31 = E0 .

Studiare in dettaglio i casi

ω0 = 1, ω1 =√5−12 ;

ω0 = 1, ω1 = 1;

ω0 = 1, ω1 = 2;

ω0 = 1, ω1 = −√5−12 .

Esercizi

Il problema di Keplero nel piano

Consideriamo l’Hamiltoniana

H =1

2(p21 + p22)− 1

r, r =

√q21 + q22 .

con i seguenti dati iniziali

q1(0) = 1− e , q2(0) = 0 , p1(0) = 0 , p2(0) =

√1 + e

1− e,

per eccentricita 0 ≤ e < 1.

Integrare il problema con lo schema di Verlet e discutere ilcomportamento per tempi “lunghi”.

Discutere i risultati ottenuti (suggerimento, tutte gli integrali primi sonoben conservati durante l’integrazione?).

Metodi splitting

Consideriamo un’Hamiltoniana

H = H1 + H2

con H1 e H2 integrabili.

AttenzioneTipicamente la somma di due Hamiltoniane integrabili non e integrabile.

L’dea dei metodi “splitting” e quella di trattare separatamente i flussiΦH1(t) e ΦH1(t), componendoli in modo opportuno. Ad esempio

ΦH1(t/2) ◦ ΦH2(t) ◦ ΦH1(t/2)

Il flusso di un’Hamiltoniana e ovviamente una trasformazione simplettica,quindi i metodi “splitting” sono simplettici per costruzione.

Contopoulos — Henon-Heiles

Consideriamo nuovamente l’Hamiltoniana

H =ω0

2(q20 + p20) +

ω1

2(q21 + p21) + q20q1 −

1

3q31

Il metodo leapfrog e sostanzialmente dato dallo splitting

H1 =ω0

2p20 +

ω1

2p21 , H2 =

ω0

2q20 +

ω1

2q21 + q20q1 −

1

3q31

EsercizioConfrontare i risultati ottenuti con lo schema leapfrog e quelli dati dalseguente splitting

H1 =ω0

2(q20 + p20) +

ω1

2(q21 + p21) , H2 = q20q1 −

1

3q31

(suggerimento, integrare esattamente gli oscillatori armonici!)

Approssimazione di Landau-Teller

Obiettivo di questa parte e di riprodurre parte dei risultati dell’articolohttp://www.mat.unimi.it/users/sansotte/pdf_files/hamsys/

BenCarSem-1993.pdf

Consideriamo un’Hamiltoniana

H =p2x2

+p2y2

+ω2

2y2 + V (r) , r = x − y .

Consideriamo il potenziale

V (r) = V0d

re−r/d

Vogliamo quantificare lo scambio di energia tra i gradi di liberta legatiall’energia vibrazionale della molecola e quelli traslazionali, in funzionedella frequenza ω (si veda figura 4 a pagina 187 dell’articolo).

Approssimazione di Landau-Teller

Esercizio 1/2

H =p2x2

+p2y2

+ω2

2y2 + V0

d

re−r/d , r = x − y .

Fissiamo V0 = 1, ω = 1, d = 1 e x0 = 20d (provare anche x0 = 70d).

Diamo come energie iniziali

Ex(0) =p2x2

= 2 , Ey (0) = ωI = 1 .

py =√ω√

2I cosφ , y =1√ω

√2I sinφ .

Quando x torna al punto di partenza valutiamo ∆Ex(φ).

Scegliamo fasi iniziali equispaziate e calcoliamo le componenti di Fourier

Ek =1

N

N−1∑l=0

∆E(2πl

N

)e−2iπkl/N

Approssimazione di Landau-Teller

Esercizio 2/2

considerare lo “splitting” tipo leapfrog (posizioni-momenti);

considerare lo “splitting”

H1 =p2x2

+p2y2

+ω2

2y2 , H2 = V0

d

re−r/d ;

riportare in un grafico il decadimento dei moduli delle componenti diFourier |E0|, |E1|, |E2| ed |E3|, al variare della frequenza ω;

riportare in un grafico il comportamento dello scambio di energia, alvariare della fase iniziale;

Integratori simplettici per Hamiltoniane quasi-integrabili

Obiettivo di questa parte e di riprodurre parte dei risultati dell’articolohttp://www.mat.unimi.it/users/sansotte/pdf_files/hamsys/

LasRob-2000.pdf

Consideriamo un’Hamiltoniana

H = A + εB

con A e B Hamiltoniane integrabili.

L’idea e di sfruttare il piccolo parametro ε per migliorare la precisionedegli schemi di integrazione di tipo splitting.

Integratori simplettici per Hamiltoniane quasi-integrabili

SABA2n = ec1τLAed1τLεB . . . ednτLεB ecn+1τLAednτLεB . . . ed1τLεB ec1τLA

SABA2n+1 = ec1τLAed1τLεB . . . ecn+1τLAedn+1τLεB ecn+1τLA . . . ed1τLεB ec1τLA

SBAB2n = ed1τLεB ec2τLA . . . ecn+1τLAedn+1τLεB ecn+1τLA . . . ec2τLAed1τLεB

SBAB2n+1 = ed1τLεB ec2τLA . . . edn+1τLεB ecn+2τLAedn+1τLεB . . . ec2τLAed1τLεB

SABA3

c1 =1

2−√

15

10, c2 =

√15

10, d1 =

5

18, d2 =

4

9.

SBAB3

c2 =1

2−√

5

10, c3 =

√5

5, d1 =

1

12, d2 =

5

12.

Correttore

Gli schemi di integrazione SABAn e SBABn hanno una precisioneO(τ 2nε+ τ 2ε).

E possibile migliorare ulteriormente questi schemi introducendo uncorrettore che elimini il termine O(τ 2ε).

SnC = e−(τ/2)cnL{{A,B},B}Sne−(τ/2)cnL{{A,B},B}

c3 =54− 13

√5

648(SABA3) c3 =

13− 5√

5

288(SBAB3)

EsercizioImplementare il metodo SBAB3C per il modello di Henon-Heiles.

ESERCIZI DI

INTEGRAZIONE

NUMERICA

Dinamica molecolare (Lennard-Jones) 1/2

Consideriamo N particelle identiche che interagiscono a coppie con unpotenziale

φ(r) = ε

(( rr

)12− 2( rr

)6), r = ‖q1 − q2‖ .

Fissiamo ε = 0.4 e r = 1. Consideriamo N = 400 particelle sui nodi di unreticolo cartesiano (di lato unitario) 20× 20 e perturbiamo le posizioniiniziali in modo random uniforme di modulo massimo 0.01.

SE le posizioni non perturbate corrispondessero ad un equilibrio,dovremmo osservare oscillazioni random attorno ai valori di equilibrio.

Simulare la dinamica e discutere i risultati.

Dinamica molecolare (Lennard-Jones) 2/2

Definiamo l’intorno Ii = {j 6= i : ‖qj − qi‖ < ∆}, ad esempio fissando∆ = 1.3r . Per ogni coppia di particelle nell’intorno Ii , definiamo l’angoloθj,i,k con vertice in qi e lati (qj − qi ), (qk − qi ) e per m ≥ 1 calcoliamo

µm =N∑i=1

∑j,k∈Ii

cos(mθjik)

Studiare l’evoluzione di µ4 e µ6 e discutere i risultati.

Sistema SolareConsideriamo un modello di sistema solare costituito da 10 corpi: Sole,Mercurio, Venere, Terra, Marte, Giove, Saturno, Urano, Nettuno ePlutone.

Nel sistema baricentrico, l’Hamiltoniana assume la forma

H =1

2

10∑i=1

‖pi‖2

mi− G

9∑i=1

10∑j=i+1

mimj

‖qi − qj‖

I dati iniziali, generati tramite l’interfacciahttps://ssd.jpl.nasa.gov/horizons.cgi

sono contenuti nella cartella compressahttp://www.mat.unimi.it/users/sansotte/hamsys/solar_

system_2440400-50.tgz

Simulare la dinamica e discutere i risultati.