equazioni differenziali Soluzione di sistemi di -...
Transcript of equazioni differenziali Soluzione di sistemi di -...
![Page 1: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/1.jpg)
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
![Page 2: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/2.jpg)
Laboratorio di Calcolo B 7
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 =+++′′+′+ hh txtatxtatxtatxta n
n
![Page 3: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/3.jpg)
Laboratorio di Calcolo B 8
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
![Page 4: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/4.jpg)
Laboratorio di Calcolo B 9
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.
�
�
+−+
++−′′+−′+=
nn xxxfn
xxxfxxxfxfxf
))((!
1
))((21))(()()(
00)(
200000
![Page 5: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/5.jpg)
Laboratorio di Calcolo B 10
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)
![Page 6: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/6.jpg)
Laboratorio di Calcolo B 11
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 ∆+=∆
Primo Ordine )()0()0()( 2tOtxxtx ∆+∆′+=∆
Si assume che l’accelerazione sia nulla (cioè che la velocità resti costante)
![Page 7: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/7.jpg)
Laboratorio di Calcolo B 12
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 ∆+=∆
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)
![Page 8: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/8.jpg)
Laboratorio di Calcolo B 13
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(2
2)0( xtxt
x
![Page 9: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/9.jpg)
Laboratorio di Calcolo B 14
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
![Page 10: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/10.jpg)
Laboratorio di Calcolo B 15
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
![Page 11: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/11.jpg)
Laboratorio di Calcolo B 16
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
∆+++++=∆
∆+∆=
∆+∆=
∆+∆=
∆=
![Page 12: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/12.jpg)
Laboratorio di Calcolo B 17
Soluzione di equazioni Soluzione di equazioni differenziali: riassuntodifferenziali: riassunto
x(0)
È nota la posizione iniziale (condizioni al contorno)
![Page 13: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/13.jpg)
Laboratorio di Calcolo B 18
Soluzione di equazioni Soluzione di equazioni differenziali: riassuntodifferenziali: riassunto
x(0)Si usa l’equazione differenziale per calcolare la derivata prima nel punto iniziale
x’(0) = f(x(0),0)
![Page 14: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/14.jpg)
Laboratorio di Calcolo B 19
Soluzione di equazioni Soluzione di equazioni differenziali: riassuntodifferenziali: riassunto
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
![Page 15: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/15.jpg)
Laboratorio di Calcolo B 20
Soluzione di equazioni Soluzione di equazioni differenziali: riassuntodifferenziali: riassunto
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)
![Page 16: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/16.jpg)
Laboratorio di Calcolo B 21
Soluzione di equazioni Soluzione di equazioni differenziali: riassuntodifferenziali: riassunto
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
![Page 17: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/17.jpg)
Laboratorio di Calcolo B 22
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.
![Page 18: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/18.jpg)
Laboratorio di Calcolo B 23
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
),(),()()()()()()()()(
![Page 19: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/19.jpg)
Laboratorio di Calcolo B 24
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 pianeta.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.
![Page 20: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/20.jpg)
Laboratorio di Calcolo B 25
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)( 11�
�
��
�
�
�
� ′′=′′
( )
=′
=′
ttxtxtvtvftv
tvtx
NNii
ii
),(,),(),(,),()(
)()(
11�
�
��
�
�
�
�
��
![Page 21: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/21.jpg)
Laboratorio di Calcolo B 26
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
∆++=∆∆++=∆
∆++∆=
+∆=
∆=∆=
���
�
��
�
�
�
�
�
�
�
�
�
�
�
�
�
��
�
�
�
��
�
![Page 22: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/22.jpg)
Laboratorio di Calcolo B 27
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
∆++=∆∆++=∆
∆++∆=
+∆=
∆=∆=
���
�
��
�
�
�
�
�
�
�
�
�
�
�
�
�
��
�
�
�
��
�
Passo 1Si calcola il valore
degli x’i(0)Banale: x’i(0) = vi(0)
![Page 23: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/23.jpg)
Laboratorio di Calcolo B 28
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
∆++=∆∆++=∆
∆++∆=
+∆=
∆=∆=
���
�
��
�
�
�
�
�
�
�
�
�
�
�
�
�
��
�
�
�
��
�
Passo 2Si calcola il valore dei
v’i(0) usando l’equazione
![Page 24: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/24.jpg)
Laboratorio di Calcolo B 29
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
∆++=∆∆++=∆
∆++∆=
+∆=
∆=∆=
���
�
��
�
�
�
�
�
�
�
�
�
�
�
�
�
��
�
�
�
��
�
Passo 3Si estrapolano xi e vi al primo ordine a ∆∆∆∆t/2 e si
calcola x’i(∆∆∆∆t/2)
![Page 25: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/25.jpg)
Laboratorio di Calcolo B 30
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
∆++=∆∆++=∆
∆++∆=
+∆=
∆=∆=
���
�
��
�
�
�
�
�
�
�
�
�
�
�
�
�
��
�
�
�
��
�
Passo 4Si calcola v’i(∆∆∆∆t/2)
![Page 26: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/26.jpg)
Laboratorio di Calcolo B 31
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
∆++=∆∆++=∆
∆++∆=
+∆=
∆=∆=
���
�
��
�
�
�
�
�
�
�
�
�
�
�
�
�
��
�
�
�
��
�
Passo 5Si calcola xi(∆∆∆∆t) al
secondo ordine in ∆∆∆∆t
![Page 27: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/27.jpg)
Laboratorio di Calcolo B 32
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
∆++=∆∆++=∆
∆++∆=
+∆=
∆=∆=
���
�
��
�
�
�
�
�
�
�
�
�
�
�
�
�
��
�
�
�
��
�
Passo 6Si calcola vi(∆∆∆∆t) al
secondo ordine in ∆∆∆∆t
![Page 28: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/28.jpg)
Laboratorio di Calcolo B 33
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.
![Page 29: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/29.jpg)
Laboratorio di Calcolo B 34
Schema a blocchiSchema a blocchiInput 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’accelerazione del corpo i-mo
Grafico dei punti
IN: N,D,t, xi(t), vi(t)
OUT: grafico
t += ∆tt > Tmax ? End
Start
t = Tmin
![Page 30: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/30.jpg)
Laboratorio di Calcolo B 35
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 = (double*)malloc(N*D*sizeof(double));
Loopint n,d;for (n=0; n<N; n++) {
for (d=0; d<D; d++) {vect[n*D + d] = ....
}}
![Page 31: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/31.jpg)
Laboratorio di Calcolo B 36
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 = (double**)malloc(N*sizeof(double*));for (n=0; n<N; n++) {
vect[n] = (double *)malloc(D*sizeof(double));}
Loopint n,d;for (n=0; n<N; n++) {
for (d=0; d<D; d++) {vect[n] [d] = ....
}}
![Page 32: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/32.jpg)
Laboratorio di Calcolo B 37
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
Componenti x, y,zdel corpo 2
Componenti x, y,zdel corpo 3
...
Soluzione 1: vettore di DxN componenti
![Page 33: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/33.jpg)
Laboratorio di Calcolo B 38
Manipolazione dei vettoriManipolazione dei vettori
vect[0]vect[1]vect[2]vect[3]
vect[0][0]vect[0][1]vect[0][2]
Componenti x, y,zdel corpo 1
...
vect[1][0]vect[1][1]vect[1][2]
Componenti x, y,zdel corpo 2
vect[2][0]vect[2][1]vect[2][2]
Componenti x, y,zdel corpo 3
Soluzione 2: matrice [N][D]
![Page 34: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/34.jpg)
Laboratorio di Calcolo B 39
Visualizzazione delle traiettorieVisualizzazione delle traiettoriePer visualizzare le traiettorie in due dimensioni si può utilizzare la funzioneR_Graph(int n, double *x, double *y, ‘ ‘)che gia conoscete. Se ∆t è piccolo e la dimensione del marker pure (R_SetMarkerSize) questa funzione disegna una linea continua lungo la traiettoria.Allo scopo di fare spettacolo, più che di capire il problema fisico, si possono utilizzare
Queste funzioni identificano tramite il parametro “obj”un particolare grafico e ne consentono la modifica in tempi successivi. In questo modo potete implementare il movimento di palline in tre dimensioni per rappresentare i punti.
R_GraphObj(char *name, int n, double *x, double *y, ‘M’)R_3DGraphObj(char *obj, int n, double *x, double *y, double *z, ‘M’)
![Page 35: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/35.jpg)
Laboratorio di Calcolo B 40
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.
![Page 36: equazioni differenziali Soluzione di sistemi di - ge.infn.itparodifa/labc/aa0102/Trasp/EqDiff.pdf · equazioni differenziali Potreste avere l’impressione di non sapere nulla sulle](https://reader031.fdocumenti.com/reader031/viewer/2022020121/5c6573d009d3f2916e8cc564/html5/thumbnails/36.jpg)
Laboratorio di Calcolo B 41
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’accelerazione del 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