Soluzione di sistemi di equazioni differenziali -...

Laboratorio di Calcolo B 6

Soluzione di sistemi di Soluzione di sistemi di equazioni differenzialiequazioni differenziali

Potreste avere l’impressione di non sapere nulla sulle equazioni differenziali, e di non averne mai incontrata una. In realtà questa impressione è sbagliata, perché la legge di Newton

F = masi può scrivere nella forma

F(x(t),t) = mx’’(t)che mostra esplicitamente come si tratta in sostanza di una equazione differenziale del secondo ordine.Nel seguito ci riferiremo in modo esplicito alla risoluzione delle equazioni del moto, anche se i metodi illustrati sono di validità generale


Equazioni differenziali lineari Equazioni differenziali lineari Si tratta di equazioni del tipo:

dove x(t) è la funzione incognita e le an(t) sono funzioni note. In generale, un’equazione differenziale non possiede una singola soluzione (ovvero una singola funzione x(t) che la soddisfi), ma una famiglia di soluzioni. La soluzione appropriata per un dato problema è quella che rispetta le condizioni al contorno. Nel caso si stia studiando una evoluzione temporale, le condizioni al contorno sono date come i valori della funzione e delle sue derivate fino alla (n-1)-esima ad un dato istante (ad esempio per t=0).

0)()()()()()()()( )(210 =+++′′+′+ KK txtatxtatxtatxta n

n


Il metodo di Il metodo di EuleroEuleroConsideriamo un’ equazione del primo ordine:

cui è associata la condizione iniziale x(0) = x0.L’equazione medesima ci fornisce il valore della derivata prima (la velocità) al tempo t=0:

Se assumiamo che in un breve intervallo ∆t la velocità resti costante potremo scrivere, applicando la formula del moto rettilineo uniforme:

Questa relazione, nota la posizione iniziale x(0) ci fornisce x(∆t) al primo ordine. La procedura può essere iterata per ricavare x(t) a t generico. La precisione è modesta.

( )ttxftbtxtatx ),()()()()( =+=′

( )0),0()0( xfx =′

txfxtxxtx ∆+=∆′+=∆ )0,()0()0()( 00


La serie di La serie di TaylorTaylorPer chiarire che cose si intende con “ricavare x(t) al primo ordine” dobbiamo introdurre la serie di Taylor. Si può dimostrare che se una funzione f(x) e tutte le sue derivate sono continue in un intervallo vale la relazione

dove x e x0 sono punti generici dell’intervallo. In altre parole la serie di Taylor ci garantisce che la conoscenza del valore di una funzione e di tutte le sue derivate in un punto x0 è sufficiente a determinare il valore della funzione in un punto generico x.

K

K

+−+

++−′′+−′+=

nn xxxfn

xxxfxxxfxfxf

))((!1

))((21))(()()(

00)(

200000


Precisione degli sviluppi di Precisione degli sviluppi di TaylorTaylorSe applichiamo la serie di Taylor al nostro problema di estrapolare il moto di un corpo note la posizione e le sue derivate calcolate al tempo t=0, possiamo ottenere precisioni differenti troncando la serie in diverse posizioni.

Ordine Zero )()0()( tOxtx ∆+=∆

Si assume che la velocità sia nulla (cioè che la posizione resti costante)




Primo Ordine )()0()0()( 2tOtxxtx ∆+∆′+=∆

Si assume che l’accelerazione sia nulla (cioè che la velocità resti costante)




Primo Ordine )()0()0()( 2tOtxxtx ∆+∆′+=∆

Secondo Ordine )()0(21)0()0()( 32 tOtxtxxtx ∆+∆′′+∆′+=∆

Si assume che la derivata dell’accelerazione sia nulla (cioè che l’accelerazione resti costante)


Il metodo di Il metodo di RungeRunge--KuttaKuttaPer ottenere una precisione maggiore occorre eseguire uno sviluppo al secondo ordine:

Abbiamo però il problema di ricavare x’’(0).Sviluppiamo al prim’ordine x’(∆t/2):

da cui:

2)0(21)0()0()( txtxxtx ∆′′+∆′+=∆

2)0()0(

2txxtx ∆′′+′=

∆′

′−

∆′

∆=′′ )0(

22)0( xtxt

x


Il metodo di Il metodo di RungeRunge--KuttaKuttaSe sostituiamo nell’espressione al secondo ordine di x(∆t) otteniamo:

Il valore della velocità nel punto ∆t/2 può essere ricavato calcolando la posizione al prim’ordine con il metodo di Eulero ed utilizzando l’equazione differenziale:

)(2

)0()( 3tOttxxtx ∆+∆

∆′+=∆

∆

∆=

∆′

∆′+=

∆

2,

22

2)0()0(

2ttxftx

txxtx


Il metodo di Il metodo di RungeRunge--KuttaKuttaRiassumendo possiamo dire che il metodo di Runge-Kutta di ordine due consiste nell’eseguire una estrapolazione del primo ordine da x(0) a x(∆t/2), nel valutare la derivata x’(∆t/2) e nell’utilizzarla per ottenere una stima di x(∆t) esatta al secondo ordine.Vediamo la sequenza di calcolo:

)()0()(2,2

)0(

)0),0((

32

12

1

tOkxtx

tkxftk

xftk

∆++=∆

∆

+∆=

∆=Calcolo della derivata nel punto iniziale

Estrapolazione al primo ordine in ∆t/2 e calcolo della derivata

Estrapolazione al secondo ordine in ∆t


RungeRunge--KuttaKutta di ordine quattrodi ordine quattroEsiste un metodo di Runge-Kutta esatto al quart’ordine, del quale riportiamo solo la sequenza di calcolo, ma che è sempre conveniente usare:

)(6336

)0()(

),)0((2,

2)0(

2,

2)0(

)0),0((

54321

34

23

12

1

tOkkkkxtx

tkxftk

tkxftk

tkxftk

xftk

∆+++++=∆

∆+∆=

∆

+∆=

∆

+∆=

∆=


Soluzione di equazioni Soluzione di equazioni differenziali: riassuntodifferenziali: riassunto

x(0)

È nota la posizione iniziale (condizioni al contorno)



x(0)Si usa l’equazione differenziale per calcolare la derivata prima nel punto iniziale

x’(0) = f(x(0),0)



x(0) Si assume che la derivata prima resti costante, e si calcola lo spostamento nell’intervallo ∆t (al primo ordine)

x’(0)

x(∆t) = x(0) + x’(0)∆t



x(0)

Alternativamente ci si ferma a metà strada e si calcola la derivata prima in ∆t/2 con l’eq. differenziale

x’(0)

x(∆t) + O(∆t2)

x(∆t/2) = x(0) + x’(0)∆t/2

x’(∆t/2) = f(x(∆t/2),∆t/2)



x(0)

Sfruttando il fatto che

possiamo ricavare il valore della derivata seconda nel punto iniziale e quindi determinare la posizione a ∆t al secondo ordine

x’(0)

x(∆t) + O(∆t2)

x(∆t/2)

x’(∆t/2)

( )[ ])0(2/)/2()0( xtxtx ′−∆′∆=′′

x(∆t) = x(0) + x’(0)∆t + x’’(0)∆t2/2 +O(∆t3)


Sistemi di equazioni Sistemi di equazioni differenzialidifferenziali

Ci sono tre motivi per cui ci dobbiamo occupare di sistemi di equazioni differenziali e non di equazioni singole:

1) i problemi di dinamica del corpo rigido sono problemi vettoriali, per cui avremo tre equazioni per ciascuna delle componenti cartesiane;

2) i problemi interessanti da un punto di vista fisico coinvolgono più di un corpo, e quindi avremo tre equazioni vettoriali per ciascun corpo coinvolto;

3) l’equazione del moto è una equazione differenziale del secondo ordine, mentre noi abbiamo un metodo di soluzione per le equazioni del primo ordine.


Equazioni del secondo ordineEquazioni del secondo ordineUna equazione lineare del secondo ordine è del tipo:

Possiamo trasformare questa equazione in un sistema di due equazioni del primo ordine; basta porre

e si ottiene:

che è un sistema del primo ordine nelle variabili x(t) e v(t). Naturalmente, siccome ci sono due equazioni saranno necessarie due condizioni iniziali: x(0) e v(0).Il trucchetto funziona anche per gli ordini superiori.

( )ttxtxftctxtbtxtatx ),(),()()()()()()( ′=++′=′′

)()( tvtx =′

( )

=++=′=′

ttxtvftctxtbtvtatvtvtx

),(),()()()()()()()()(


Esercitazione in laboratorioEsercitazione in laboratorioIn laboratorio dovrete realizzare un programma per la soluzione di un generico problema del moto di N corpi in D dimensioni (con D = 1, 2 o 3). Questo programma potrà essere applicato a diversi problemi fisici interessanti.A scopo di esempio, e per verificare che tutto funzioni bene, la procedura di calcolo dovrà essere applicata al problema della determinazione, in due dimensioni, della traiettoria di un satellite attorno ad un pianeta assunto fermo.Dovrete anche realizzare il grafico della traiettoria del satellite.Nel seguito ci concentreremo nell’analisi dettagliata di questo problema, e specializzeremo un poco (ma non molto) la generica procedura valida per qualunque sistema di equazioni differenziali lineari.


Il calcolo della traiettoriaIl calcolo della traiettoriaIl problema che dobbiamo risolvere é:

che, come abbiamo visto, va riscritto nella forma:

Se abbiamo N corpi (i=1...N) e risolviamo il problema in tre dimensioni si tratta di un sistema di 2x3xN equazioni del primo ordine.Dobbiamo osservare che la quasi totalità dell’ interdipendenza tra le equazioni del sistema risiede nel termine di forza: la componente j-ma della forza che agisce sul punto i è funzione di tutte le componenti delle posizioni e delle velocità di tutti gli Ncorpi.

( )ttxtxtxtxFm

tx NNii

i ),(,),(),(,),(1)( 11r

Krr

Krrr ′′=′′

( )

=′

=′

ttxtxtvtvftv

tvtx

NNii

ii

),(,),(),(,),()(

)()(

11r

Krr

Krrr

rr


Soluzione del sistemaSoluzione del sistemaVisto che siamo capaci di risolvere una singola equazione differenziale, la soluzione di un sistema di molte equazioni non presenta particolari problemi. Di fatto il metodo che usiamo è rigidamente sequenziale, per cui si deve solo fare attenzione ad eseguire ogni passo del calcolo su tutte le equazioni prima di passare al passo successivo. Proviamo con il metodo di Runge-Kutta di ordine due.

)()0()()()0()(2,,

2)0(,,

2)0(

2)0(

)0),0(,,)0(),0(,),0(()0(

32

32

111

1112

12

1111

tOwvtvtOkxtx

tkxwvftwwvtk

xxvvftwvtk

iiiiii

iii

ii

NNiiii

∆++=∆∆++=∆

∆++∆=

+∆=

∆=∆=

rrrrrr

K

rr

Kr

rrrr

rr

rK

rrK

rrrrr


Soluzione del sistemaSoluzione del sistema

)()0()()()0()(2,,

2)0(,,

2)0(

2)0(

)0),0(,,)0(),0(,),0(()0(

32

32

111

1112

12

1111

tOwvtvtOkxtx

tkxwvftwwvtk

xxvvftwvtk

iiiiii

iii

ii

NNiiii

∆++=∆∆++=∆

∆++∆=

+∆=

∆=∆=

rrrrrr

K

rr

Kr

rrrr

rr

rK

rrK

rrrrr

Passo 1Si calcola il valore

degli x’i(0)Banale: x’i(0) = vi(0)



)()0()()()0()(2,,

2)0(,,

2)0(

2)0(

)0),0(,,)0(),0(,),0(()0(

32

32

111

1112

12

1111

tOwvtvtOkxtx

tkxwvftwwvtk

xxvvftwvtk

iiiiii

iii

ii

NNiiii

∆++=∆∆++=∆

∆++∆=

+∆=

∆=∆=

rrrrrr

K

rr

Kr

rrrr

rr

rK

rrK

rrrrr

Passo 2Si calcola il valore dei

v’i(0) usando l’equazione



)()0()()()0()(2,,

2)0(,,

2)0(

2)0(

)0),0(,,)0(),0(,),0(()0(

32

32

111

1112

12

1111

tOwvtvtOkxtx

tkxwvftwwvtk

xxvvftwvtk

iiiiii

iii

ii

NNiiii

∆++=∆∆++=∆

∆++∆=

+∆=

∆=∆=

rrrrrr

K

rr

Kr

rrrr

rr

rK

rrK

rrrrr

Passo 3Si estrapolano xi e vi al primo ordine a ∆t/2 e si

calcola x’i(∆t/2)



)()0()()()0()(2,,

2)0(,,

2)0(

2)0(

)0),0(,,)0(),0(,),0(()0(

32

32

111

1112

12

1111

tOwvtvtOkxtx

tkxwvftwwvtk

xxvvftwvtk

iiiiii

iii

ii

NNiiii

∆++=∆∆++=∆

∆++∆=

+∆=

∆=∆=

rrrrrr

K

rr

Kr

rrrr

rr

rK

rrK

rrrrr

Passo 4Si calcola v’i(∆t/2)



)()0()()()0()(2,,

2)0(,,

2)0(

2)0(

)0),0(,,)0(),0(,),0(()0(

32

32

111

1112

12

1111

tOwvtvtOkxtx

tkxwvftwwvtk

xxvvftwvtk

iiiiii

iii

ii

NNiiii

∆++=∆∆++=∆

∆++∆=

+∆=

∆=∆=

rrrrrr

K

rr

Kr

rrrr

rr

rK

rrK

rrrrr

Passo 5Si calcola xi(∆t) al

secondo ordine in ∆t



)()0()()()0()(2,,

2)0(,,

2)0(

2)0(

)0),0(,,)0(),0(,),0(()0(

32

32

111

1112

12

1111

tOwvtvtOkxtx

tkxwvftwwvtk

xxvvftwvtk

iiiiii

iii

ii

NNiiii

∆++=∆∆++=∆

∆++∆=

+∆=

∆=∆=

rrrrrr

K

rr

Kr

rrrr

rr

rK

rrK

rrrrr

Passo 6Si calcola vi(∆t) al

secondo ordine in ∆t


Soluzione del sistemaSoluzione del sistemaIn conclusione, il calcolo della traiettoria di un sistema di N corpi in D dimensioni si riduce all’esecuzione iterativa di una sequenza di calcoli:

• La sequenza non dipende dallo specifico problema in esame, in quanto dipende solo dalla cinematica.

• La sequenza richiede il calcolo delle forze agenti su ogni corpo. Questa è la sola parte dove emergono le peculiarità della dinamica del sistema

• La sequenza va eseguita in parallelo su DxNequazioni. Di conseguenza è necessario manipolare vettori o matrici a DxN componenti.


Schema a blocchiSchema a blocchi

Input dei dati

IN: Nome del fileOUT: N,D

Tmin, Tmax, ∆txi(Tmin), vi(Tmin)

Cinematica

IN: N, D,t, ∆t, xi(t), vi(t)funzione Dinamica

OUT: xi(t+ ∆t), vi(t+∆t)

Dinamica

IN: i, dt, xi(t), vi(t)

OUT: componente d dell’accelerazionedel corpo i-mo

Grafico dei punti

IN: N,D,t, xi(t), vi(t)

OUT: grafico

t += ∆tt > Tmax ? End

Start

t = Tmin


Manipolazione dei vettoriManipolazione dei vettoriIl programma deve manipolare vettori di DxNcomponenti; sia D che N possono variare, per cui è necessario allocare dinamicamente. Ci sono due soluzioni possibili:Soluzione 1: vettore di DxN componenti

Allocazionedouble *vect;vect = new double[N*D];

Loopint n,d;for (n=0; n<N; n++) {

for (d=0; d<D; d++) {vect[n*D + d] = ....

}}


Manipolazione dei vettoriManipolazione dei vettoriIl programma deve manipolare vettori di DxNcomponenti; sia D che N possono variare, per cui è necessario allocare dinamicamente. Ci sono due soluzioni possibili:Soluzione 2: matrice [N][D]

Allocazionedouble **vect;int n;vect = new double*[N];for (n=0; n<N; n++) {

vect[n] = new double[D];}

Loopint n,d;for (n=0; n<N; n++) {

for (d=0; d<D; d++) {vect[n] [d] = ....

}}


Manipolazione dei vettoriManipolazione dei vettori

vect[0]vect[1]vect[2]vect[3]vect[4]vect[5]vect[6]vect[7]vect[8]

Componenti x, y,zdel corpo 1



...

Soluzione 1: vettore di DxN componenti


Manipolazione dei vettoriManipolazione dei vettori

vect[0]vect[1]vect[2]vect[3]

vect[0][0]vect[0][1]vect[0][2]


...





Soluzione 2: matrice [N][D]


Visualizzazione delle traiettorieVisualizzazione delle traiettoriePer visualizzare le traiettorie in due dimensioni si può utilizzare un oggetto di tipo TGraph:

TGraph(int n, double *x, double *y)che gia conoscete. Se ∆t è piccolo e la dimensione del marker pure (metodo SetMarkerSize) questa funzione disegna una linea continua lungo la traiettoria.Allo scopo di fare spettacolo, più che di capire il problema fisico, potete utilizzare il metodo

che consente lo spostamento dei punti, e quindi realizza una sorta di animazione del moto dei corpi lungo la traiettoria calcolata. Ricordate che, dopo la chiamata a SetPoint dovete chiamare il metodo Update di TCanvaso TPad per aggiornare il grafico. Se poi siete appassionati di effetti speciali potete creare un sistema di assi 3D con un oggetto di tipo TView e poi un grafico 3D con TPolyMarker3D.

SetPoint(int n, double *x, double *y)


Verifica della precisione di Verifica della precisione di calcolocalcolo

È chiaro, dalla trattazione generale sulla soluzione numerica delle equazioni differenziali, che la precisione è tanto migliore quanto più piccolo è ∆t. Il problema che si pone è quello di scegliere un valoresensato per l’incremento temporale.Una tecnica possibile consiste nel calcolare, ad ogni passo, l’energia totale del sistema. Tale energia dovrebbe rimanere costante. Se si individuano delle variazioni significa che la precisione è insufficiente (o che il sistema non è conservativo). In questo caso si può intervenire riducendo ∆t.


Schema a blocchi completoSchema a blocchi completo

Cinematica

IN: N, D,t, ∆t, xi(t), vi(t)funzione Dinamica

OUT: xi(t+ ∆t), vi(t+∆t)

Dinamica

IN: i, d, N, Dt, xi(t), vi(t)

OUT: componente d dell’accelerazionedel corpo i-mo

Grafico dei punti

IN: N,D,t, xi(t), vi(t)

OUT: grafico

t += ∆tt > Tmax ? End

Energia

IN: N, D∆t, t, xi(t), vi(t)

OUT: ∆t

Input dei dati

IN: Nome del fileOUT: N,D

Tmin, Tmax, ∆txi(Tmin), vi(Tmin)

Start

t = Tmin

Soluzione di sistemi di equazioni differenziali -...

Documents

Transcript of Soluzione di sistemi di equazioni differenziali -...