Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. differenziali di primo grado

Problema soluzione numerica dell’equazione differenziale:

x′(t) = f(x(t), t)

dato x0 = x(t0).

Soluzione iterativa utilizzando l’equazione differenziale e lo sviluppo inserie di Taylor: relazione tra posizioni a distanza temporale h:xn+1 ↔ xn (dove xn+1 ≡ x(tn + h) e xn ≡ x(tn))

1 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. differenziali di primo grado

Problema soluzione numerica dell’equazione differenziale:

x′(t) = f(x(t), t)

dato x0 = x(t0).

Soluzione iterativa utilizzando l’equazione differenziale e lo sviluppo inserie di Taylor: relazione tra posizioni a distanza temporale h:xn+1 ↔ xn (dove xn+1 ≡ x(tn + h) e xn ≡ x(tn))

Metodo di Eulero

k1 = hf(xn, tn) derivata nel punto iniziale

xn+1 = xn + k1 +O(h2) estrapolazione a tn + h

2 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. differenziali di primo grado

Problema soluzione numerica dell’equazione differenziale:

x′(t) = f(x(t), t)

dato x0 = x(t0).

Soluzione iterativa utilizzando l’equazione differenziale e lo sviluppo inserie di Taylor: relazione tra posizioni a distanza temporale h:xn+1 ↔ xn (dove xn+1 ≡ x(tn + h) e xn ≡ x(tn))

Metodo di Runge Kutta al second’ordine

k1 = hf(xn, tn) derivata nel punto iniziale

k2 = hf(xn +k1

2, tn +

h

2) estrap. al punto intermedio e derivata

xn+1 = xn + k2 +O(h3) estrapolazione a tn + h

3 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. differenziali di primo grado

Problema soluzione numerica dell’equazione differenziale:

x′(t) = f(x(t), t)

dato x0 = x(t0).

Soluzione iterativa utilizzando l’equazione differenziale e lo sviluppo inserie di Taylor: relazione tra posizioni a distanza temporale h:xn+1 ↔ xn (dove xn+1 ≡ x(tn + h) e xn ≡ x(tn))

Metodo di Runge Kutta al quart’ordine

k1 = hf(xn, tn) k2 = hf

(xn +

k1

2, tn +

h

2

)k3 = hf

(xn +

k2

2, tn +

h

2

)k4 = hf(xn + k3, tn + h)

xn+1 = xn +k1

6+

k2

3+

k3

3+

k4

6+O(h5)

4 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. diff. di ordini superiori

Una generica equazione differenziale del secondo ordine e del tipo

x′′(t) = f(t, x(t), x

′(t))

Si puo trasformare in un sistema di due equazioni del primoordine ponendo

x′(t) = v(t)

si ottiene {x′(t) = v(t)

v′(t) = a(t, x(t), v(t))

che e un sistema del primo ordine nelle variabili x(t) e v(t).Ovviamente lo stesso “trucco” funziona per tutti gli ordini.

5 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. diff. per N corpi in Ddimensioni

La soluzione di un sistema di equazioni non presenta particolariproblemi.

Lo sviluppo riportato di seguito (per RK al secondo ordine) richiede che

le operazioni siano rigidamente sequenziali: si deve quindi fare

attenzione ad eseguire ogni passo del calcolo su tutte le equazioni

prima di passare al passo successivo.

I

{~k1i = h~vi(tn)~w1i = h~fi(~v1(tn), ..., ~x1(tn), ..., tn)

II

{~k2i = h

(~vi(tn) + ~w1i

2

)~w2i = h~fi(~v1(tn) + ~w11

2, ..., ~x1(tn) +

~k11

2, ..., tn + h

2)

III

{~xi(tn + h) = ~xi(tn) + ~k2i +O(h3)~vi(tn + h) = ~vi(tn) + ~w2i +O(h3)

6 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. diff. per N corpi in Ddimensioni

La soluzione di un sistema di equazioni non presenta particolariproblemi.

Lo sviluppo riportato di seguito (per RK al secondo ordine) richiede che

le operazioni siano rigidamente sequenziali: si deve quindi fare

attenzione ad eseguire ogni passo del calcolo su tutte le equazioni

prima di passare al passo successivo.

I

{~k1i = h~vi(tn)~w1i = h~fi(~v1(tn), ..., ~x1(tn), ..., tn)

II

{~k2i = h

(~vi(tn) + ~w1i

2

)~w2i = h~fi(~v1(tn) + ~w11

2, ..., ~x1(tn) +

~k11

2, ..., tn + h

2)

III

{~xi(tn + h) = ~xi(tn) + ~k2i +O(h3)~vi(tn + h) = ~vi(tn) + ~w2i +O(h3)

I: calcolo di ~x′i(tn) (banale: ~x

′i(tn) = ~v(tn))

calcolo di ~v′i(tn) usando l’equazione differenziale

7 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. diff. per N corpi in Ddimensioni

La soluzione di un sistema di equazioni non presenta particolariproblemi.

Lo sviluppo riportato di seguito (per RK al secondo ordine) richiede che

le operazioni siano rigidamente sequenziali: si deve quindi fare

attenzione ad eseguire ogni passo del calcolo su tutte le equazioni

prima di passare al passo successivo.

I

{~k1i = h~vi(tn)~w1i = h~fi(~v1(tn), ..., ~x1(tn), ..., tn)

II

{~k2i = h

(~vi(tn) + ~w1i

2

)~w2i = h~fi(~v1(tn) + ~w11

2, ..., ~x1(tn) +

~k11

2, ..., tn + h

2)

III

{~xi(tn + h) = ~xi(tn) + ~k2i +O(h3)~vi(tn + h) = ~vi(tn) + ~w2i +O(h3)

II: si estrapolano ~xi e ~vi al primo ordine a tn + h/2

e si calcolano ~x′i(tn + h/2) e ~v

′i(tn + h/2).

8 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Riassunto: soluzione equ. diff. per N corpi in Ddimensioni

La soluzione di un sistema di equazioni non presenta particolariproblemi.

Lo sviluppo riportato di seguito (per RK al secondo ordine) richiede che

le operazioni siano rigidamente sequenziali: si deve quindi fare

attenzione ad eseguire ogni passo del calcolo su tutte le equazioni

prima di passare al passo successivo.

I

{~k1i = h~vi(tn)~w1i = h~fi(~v1(tn), ..., ~x1(tn), ..., tn)

II

{~k2i = h

(~vi(tn) + ~w1i

2

)~w2i = h~fi(~v1(tn) + ~w11

2, ..., ~x1(tn) +

~k11

2, ..., tn + h

2)

III

{~xi(tn + h) = ~xi(tn) + ~k2i +O(h3)~vi(tn + h) = ~vi(tn) + ~w2i +O(h3)

III: si calcolano ~xi(tn + h) e ~vi(tn + h)

9 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Schema (un corpo una dimensione) (I)

10 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Schema (un corpo una dimensione) (II)

void cinematica(double &t, double h, double &x,double &v);

contenuto:

k1 = hv(tn)

w1 = hf(tn, x(tn), v(tn))

k2 = h(v(tn) + w1/2)

w2 = hf(tn + h/2, x(tn) + k1/2, v(tn) + w1/2)

x(tn + h) = x(tn) + k2 +O(h3)

v(tn + h) = v(tn) + w2 +O(h3)

double dinamica(double t, double x, double v);restituisce: f(tn, x(tn), v(tn)) (accelerazione)

11 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Schema (un corpo una dimensione) (III)

Come disegnare i punti ? Siccome i punti vengono aggiunti uno allavolta e utile utilizzare il costruttore default per TGraph:

TGraph::TGraph()

e poi, quando sono disponibili, aggiungere ad uno ad uno i punti con

void TGraph::SetPoint(int i, double x, double y)

Ad esempio:TGraph *gr = new TGraph;

...

gr->SetPoint(0,x0,y0);

gr->Draw(...);

for (...){

...

gr->SetPoint(i++,x,y);

...

gPad->Modified();

gPad->Update();

}

Disegna una serie di punti

TGraph *gr = new TGraph;

...

gr->SetPoint(0,x0,y0);

gr->Draw(...);

for (...){

...

gr->SetPoint(0,x,y);

...

gPad->Modified();

gPad->Update();

}

Aggiorna le coordinate del punto (→punto in movimento, animazione)

12 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Esercitazione in laboratorio

1 Risolvere numericamente con il metodo di RK al second’ordine

d2x

dt2+ β

dx

dt+ ω2x = 0

confrontandola con la soluzione analitica

x(t) = Ae−β2 tcos(ω1t + φ)

dove ω21 = ω2 − β2/4. Provate i casi β = 0 (oscillatore armonico)

e β = 1 (oscillatore smorzaro) entrambi con ω = 4.

2 Impostate il metodo di RK al quarto ordine e ripetete il punto 1

3 Affrontate uno di questi due problemi (facoltativo):

applicate un metodo per il controllo dell’erroreaffrontate il problema della dinamica di molti corpiapplicando il metodo RK al sistema Terra-Sole.

13 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Valutazione dell’errore (I)

E chiaro dalla trattazione generale sulla soluzione numerica delleequazioni differenziali, che la precisione e tanto migliore quanto piupiccolo e il passo h.

Il problema che si pone e quello di scegliere un valore sensato perl’incremento temporale.

Conservazione dell’energiaUna tecnica possibile consiste nel calcolare, ad ogni passo, l’energia delsistema. Tale energia dovrebbe rimanere costante. Se si individuanodelle variazioni (al di la di una certa soglia predeterminata) significa chele precisione e insufficiente.In questo caso si puo intervenire riducendo il passo.

14 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Valutazione dell’errore (II)

Raddoppio del passoMetodo generale (vale per sistemi conservativi e non).Assumiamo di utilizzare un metodo di Runge-Kutta di ordine n.Per ogni passo dell’iterazione si fa un passo 2h e, indipentemente, duepassi di lunghezza h. Denotiamo con y(t + 2h) la soluzione esatta econ y1 e y2 le soluzioni approssimate ottenute rispettivamente con unsolo passo pari a 2h o con due passi pari ciascuno a h.

y(t + h) = y1 + (2h)n+1α +O(hn+2)

y(t + h) = y2 + 2(h)n+1α +O(hn+2)

(dove α e un numero dell’ordine di y(n+1)/(n + 1)!)Sostituendo e indicando ∆ = y2 − y1

y(t + h) = y2 +∆

2n − 1+O(hn+2)

si ottiene il valore “esatto” dell’errore fino all’ordine O(hn+1).

Tale errore puo essere mantenuto al di sotto di un certa soglia

prefissata diminuendo h.15 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Schema (con controllo dell’errore)

16 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Sistemi di equazioni differenziali (I)

Lo studio della dinamica di N corpi in D dimensioni comporta lasoluzione di un sistema di N×D equazioni (×2 per la separazione dellevariabili x e v).

Due le soluzioni possibili

vettore 2×N×DPoco naturale: appesantisce la funzione dinamica ma rendecinematica assolutamente generale (facilmente estendibile a gradisuperiori al secondo)

2 matrici N×DMolto “fisico”: rende ben leggibile dinamica ma complicacinematica (il comportamente di x e v non e simmetrico).

17 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Sistemi di equazioni differenziali (II)

Allocazione:

double *x;

x = new double[N*D*2];

Loop:

for (int n=0;n<N;n++){

for (int d=0;d<D;d++){

x[2*n*D+d ] = ... // pos.

x[(2*n+1)*D+d] = ... // vel.

}

}

x[0] comp. x corpo 1

x[1] comp. y corpo 1

x[2] comp. z corpo 1

x[3] comp. vx corpo 1

x[4] comp. vy corpo 1

x[5] comp. vz corpo 1

x[6] comp. x corpo 2

x[7] comp. y corpo 2

x[8] comp. z corpo 2

....

Prototipi

void dinamica(double t, double *x, double *dxdt);(il vettore dxdt contiene le derivate di x)

void cinematica(double &t, double h, double *x);

Codice cinematica (runge-kutta al second’ordine):

~k1 = h~f(~x(tn), tn) ~k2 = h~f(~x(tn) +~k1

2, tn +

h

2)

~x(tn + h) = ~x(tn) + ~k2 +O(h3)

Assumo qui e nel seguito, per semplicita di notazione, che N, D e il vettore dellemasse mass[N] siano variabili globali.18 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Sistemi di equazioni differenziali (III)

Allocazione:double *x,*v;

x = new double*[N];

v = new double*[N];

for (int n=0;n<N;n++){

x[n] = new double[D];

v[n] = new double[D];

}

Loop:for (int n=0;n<N;n++){

for (int d=0;d<D;d++){

x[n][d] = ... // pos.

v[n][d] = ... // vel.

}

}

x[0][0] comp. x corpo 1

x[0][1] comp. y corpo 1

x[0][2] comp. z corpo 1

v[0][0] comp. vx corpo 1

v[0][1] comp. vy corpo 1

v[0][2] comp. vz corpo 1

x[1][0] comp. x corpo 2

x[1][1] comp. y corpo 2

x[1][2] comp. z corpo 2

....

Prototipi

double dinamica(int i, int d, double t, double **x, double **v);(dinamica ritorna la componente d-esima dell’accelerazione del corpo i-esimo)

void cinematica(double &t, double h, double **x, double **v);

Codice cinematica (runge-kutta al second’ordine):~k1i = h ~vi(tn) ~w1i = h ~fi( ~v1(tn), ..., ~x1(tn), ..., tn)

~k2i = h(

~vi(tn) +~w1i2

)~w2i = h ~fi( ~v1(tn) +

~w112

, ..., ~x1(tn) +~k112

, ..., tn + h2

)

~xi(tn + h) = ~xi(tn) + ~k2i + O(h3) ~vi(tn + h) = ~vi(tn) + ~w2i + O(h3)

19 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Suggerimenti

CinematicaCinematica deve, all’inizio, allocare dinamicamente i vettori necessari per ilmetodo di runge-kutta (k, w e gli eventuali vettori di coordinate e velocitatemporanee) e distruggerli alla fine (il metodo non e ottimale per la velocitadi esecuzione ma evita “pasticci” con le variabili globali).

DinamicaE utile schematizzare la funzione dinamica come somma delle interazione congli N-1 corpi piu un eventuale contributo di forza esterna.

double dinamica(int i, int d, double t, double **x, double **v);

double fInter(int i, int j, int d, double t, double **x, double **v);

double fExter(int i, int d, double t, double **x, double **v);

double dinamica(int i, int d, double t, double **x, double **v){

double ftot=0.0;

for (int j=0;i<N;j++){

if (i!=j) ftot+=fInter(i,j,d,x,v,t);

}

ftot+=fExter(i,d,x,v,t);

return ftot/mass[i];

}

20 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Disegno 3D

Assi:

TView::TView();

void TView::Draw();

Grafico

TPolyMarker3D::TPolyMarker3D();

void TPolyMarker3D::SetPoint(int n, double x, double y, double z)

(analogo a TGraph)

21 / 24

Riassunto

Eq. diff.:un corpo

Esercitazione

Controlloerrore

Sistemi diequazioni

Struttura a classe

Pratico per controllare creazione, distruzione e generalizzarne l’utilizzoa problemi diversi.class RungeKutta{

public:

RungeKutta(int N, int D, int rkord);

// costruttore che alloca i vettori e definisce l’ordine di RK

~RungeKutta(); // distruttore (dealloca i vettori);

void Cinematica();

void SetMass(int n, double val);

double GetMass(int n);

void SetPos(int n, int d, double val);

double GetPos(int n, int d);

void SetVel(int n, int d, double val);

double GetVel(int n, int d);

void SetT(double t0);

double GetT();

void SetStep(double h);

double GetStep();

void SetInter(double (*fInter)(..)){m_fInter = fInter;}

void SetExter(double (*fExter)(..)){m_fExter = fExter;}

...

private:

int m_N,m_D;

double m_t,m_h;

double **m_x,**m_v,*m_mass;

double **m_k1,**m_w1,...;

double (*m_fInter)(...);

double (*m_fExter)(...);

double m_dinamica(...);

...

};22 / 24