Meccanica Razionale per Allievi
Ingegneri “Aerospaziali”
Esercitazioni numeriche 2006
Realizzazione dei “Ritratti in Fase”
per alcuni problemi classici della
Meccanica
Obiettivi del Corso
Scopo del presente corso è fornire agli allievi
strumenti concettuali ed un primo training per lo
studio di problemi della Meccanica Razionale con i
metodi del calcolo computerizzato. In particolare si
fornisce una introduzione alla dinamica di sistemi
non lineari e di sistemi dotati di cicli limite e si
propone agli allievi l’uso critico di integratori
numerici per tracciare i ritratti in fase di alcuni
dispositivi classici particolarmente rappresentativi.
Programmazione “ben fatta”
• Un programma di calcolo “ben fatto” ha un
punto di inizio e una condizione di termine
ben definiti
• L’ obbiettivo del programma è definito ed
enunciato (espressamente)
• Sono assegnati i dati di ingresso e sono ben
definiti i risultati attesi
Primi programmi
• Programma uno: stampare un titolo con nome e
matricola dell’ allievo
• Programma due: scrivere su un file le coordinate di
nove punti equidistanti su una circonferenza (usare
le funzioni sin e cos)
• Programma tre : come due, ma ordinando i punti in
modo che siano i vertici di una stella, una volta
congiunti successivamente con corde opportune.
P.es. mediante l’uso del software “gnuplot” in
ambiente Linux o mediante applicazioni Excel in
ambiente Window .
Equazioni Differenziali Ordinarie
• Metodi numerici per ”Problemi di Caucy” (vedi il corso di
Analisi2).
• Vogliamo considerare equazioni differenziali:
• - lineari: particella libera, moto per gravità, crescita
Malthusiana, decadimento radioattivo, oscillatore e
repulsore armonico;
• - non lineari: pendolo, crescita di popolazioni con risorse
limitate, moto kepleriano, bipendolo;
• -dotate di ciclo limite: orologio meccanico, problema di
Lotka-Volterra, problema di Van Der Pol;
• Con esiti caotici: pendolo forzato, mappa standard
O.D.E Lineari. Part.Libera
• Particella libera: d2x/dt2 = 0 ;
Per la soluzione introduciamo la variabile v,
velocità, definita da v = dx/dt;
• l’equazione del secondo ordine (d2x/dt2 = 0) diviene
una coppia di equazioni del primo ordine:
dv/dt = 0 ; dx/dt = v ; con le condizioni iniziali
x(0) = x0 , v(0) = v0 .
• Vedi più avanti le traiettorie rettilinee parallele nel
piano delle fasi.
Moto per gravità
• d2x/dt2 = g. Come sopra introduciamo la variabile v,
velocità, definita da v = dx/dt;
• l’equazione del secondo ordine diviene una coppia
di equazioni del primo ordine: dv/dt = g; dx/dt = v ;
con le condizioni iniziali x(0) = x0 , v(0) = v0 .
• Più avanti sarà richiesto di tracciare anche per
questo caso le traiettorie nel piano delle fasi.
Oscillatore e Repulsore Armonico
L’equazione d2x/dt2 = - ω2x , con ω2 = k /m , detta
dell’oscillatore armonico, descrive il comportamento
di un punto materiale di massa m soggetto al richiamo di una forza elastica di costante k. L’equazione è
rappresentativa del comportamento di sistemi fisici
in prossimità della posizione di equilibrio stabile.
Se il coefficiente del termine lineare fosse invece
positivo, si avrebbe il “repulsore armonico”,
rappresentativo del comportamento di sistemi fisici
in prossimità degli equilibri instabili. Ne è esempio
il moto di un punto in un campo centrifugo, p.es. il
moto relativo di un punto materiale su una giostra.
Crescita Malthusiana e
decadimento nucleare
• I fenomeni sono descritti dalla equazione differenziale
ordinaria lineare del primo ordine dy/dt = A y , con A
positivo nel caso della crescita di popolazioni viventi. A
rappresenta il numero di nuovi nati nella unità di tempo
(vedi anche il modello di crescita del capitale bancario).
• Per A negativo il modello può rappresentare il decadimento
radioattivo o lo svuotamento di una spugna (o di uno strato
di terreno) impregnati di acqua (o l’evoluzione di una
popolazione senile in cui A è il numero di decessi nella unità
di tempo).
• Gli esiti numerici dei due problemi possono essere
sorprendenti in termini di precisione di calcolo e di
dominanza degli errori numerici sulla vera soluzione.
Sviluppo di Taylor per il problema di
Cauchy: crescita Malthusiana.
• Gli allievi hanno studiato ampiamente la teoria delle serie e
delle successioni nei corsi precedenti. Il Calcolo Numerico
richiede effettive utilizzazioni di tale teoria. Invitiamo a
rileggere nei testi di analisi la teoria dello sviluppo in serie
di Taylor-McLaurin per le funzioni più comuni.
• Come primo esempio si consideri l’equazione
dy/dt = A y , rappresentativa per numerosi problemi
concreti, dagli interessi bancari fino alla crescita delle
popolazioni biologiche.
• Nel nostro esempio, semplificato e standardizzato, si ponga
la condizione iniziale y(0) = 1 ed anche il coefficiente A sia
costante (ed uguale ad 1). Vogliamo capire cosa significa il
termine “ordine di precisione del solutore numerico”.
Il grafico dell’ esponenziale
• Si cerca y(T).
• Cerco uno sviluppo di Taylor, quindi ho bisogno delle
successive derivate di Y, calcolate nel punto iniziale.
• A è costante, da dy/dt = A y, si ha: d2y/dt2 = A dy/dt = A2
y; d3y/dt3 = A d2y/dt2 = A2 dy/dt = A3 y . Quindi:
y(T)= y(0) ( 1+ AT + ½ (AT)2+1/6 (AT)3+ 1/24…).
Con le consuete osservazioni per la convergenza della
serie (prodotto AT < 1) e per la valutazione del resto.
Accorgimenti per il calcolo
Una volta scelto il numero di termini della serie, il calcolo va
organizzato evitando di calcolare grandi fattoriali con i loro
inversi e potenze successive di AT, ma raccogliendo i fattori
comuni così da ridurre l’errore di troncamento che si
verifica quando si sommano termini via via più piccoli
y=yo (1+AT(1+1/2 AT(1+1/3 AT(1+1/4 AT(1+1/5 AT(…).
Come esercitazione richiesta si costruisca una procedura che
chieda in input l’ordine massimo e realizzi il prodotto dei
termini annidiati a partire dal più interno. In uscita la
procedura fornirà, oltre al risultato approssimato voluto,
anche una stima del primo termine trascurato.
Esempi : passo passo
• Per valutare y a tempi t > R (raggio di convergenza), è
possibile valutarla per t < R e quindi scegliere il nuovo
valore di y come nuovo punto iniziale per lo sviluppo in
serie di Taylor. L’espressione dello sviluppo in serie resta
la medesima ed è sufficiente aggiornare i valori della
funzione e delle derivate.
• La scelta del valore di y al termine del primo passo come
valore iniziale per un passo successivo illustra la filosofia
generale dei metodi numerici che intendiamo adottare.
• Attenzione: gli errori commessi in ciascun passo vanno
sommati in valore assoluto per ottenere la stima dell’errore
complessivo.
DA SVOLGERE PER L’ESAME
• Può convenire accorciare il passo elementare, aumentando
il numero dei passi, e simultaneamente riducendo il
numero dei termini dello sviluppo in serie. La scelta
dipende dal problema particolare ed è conveniente quando
diviene grande l’onere del computo delle successive
derivate.
• La semplicità dell’esempio suggerito consente almeno in
questi primi casi di valutare la soluzione dell’equazione
proposta dapprima con passi piuttosto lunghi, prossimi ad
R, e numerosi termini della serie di Taylor, poi con serie
di Taylor più brevi e con un maggiore numero di brevi
passi di calcolo.
Si può scegliere A=1,T=0.5. Sviluppo fino alla 4 o alla 6
potenza, poi, per confronto T =1, T=0.25,T=0.125
Esempi : seni e coseni
• Consideriamo ora la coppia di equazioni:
dx/dt= -y ; dy/dt= x ; con x(0)= 1 ; y(0)=0.
• È immediato calcolare le derivate successive delle
due espressioni e quindi ottenere
d2x/dt2= - dy/dt= - x
d2y/dt2= dx/dt= - y
d3x/dt3= - dx/dt= y
d3y/dt3= - dy/dt= -x
eccetera.
• Svolgendo quindi lo sviluppo in serie di Taylor:
x(T)= x(0) – y(0) T – ½ x(0) T2 + 1/6 y(0) T3 ……
y(T)= y(0) + x(0) T – ½ y(0) T2 – 1/6 x(0) T3 ……
Esempi : seni & coseni, segue
• E quindi si possono riorganizzare i termini:
x(T)=x(0) (1-1/2 T2+1/(2x3x4)T4–1/(2x3x4x5x6) T6 ….)–
y(0)T(1-1/(2x3)T2 +1/(2x3x4x5)T4–1/5040 T6….)
y(T)=y(0) (1-1/2 T2 + 1/24 T4 – 1/720 T6 ….) +
+ x(0)T(1-1/6 T2 +1/120 T4 – 1/5040 T6 ….)
• Anche ora è opportuno trasformare le somme di termini
rapidamente decrescenti in prodotti come p.es
x(T)=x(0) (1-(1/2 T2(1-1/(3x4)T2(1+1/(5x6) T2 (….)–
y(0)T(1-1/(2x3)T2(1+(1/(4x5)T2(1-1/(6x7)T2 (….)
etc... e quindi valutare se la scelta di compiere numerosi
passi di piccola ampiezza può essere complessivamente
vantaggiosa in termini di precisione dei risultati.
Seni & coseni, integrali primi
Risulta particolarmente utile, nello studio dei sistemi
dinamici, considerare gli eventuali “Integrali Primi” o
“Costanti del moto” ammessi dal sistema.
Per la coppia di equazioni dx/dt= -y ; dy/dt= x ; si può
notare, moltiplicando la prima per x e la seconda per y e
sommando i prodotti,
x dx/dt= -x y ; y dy/dt= y x ; d(x2)/dt +d(y2)/dt= 0
da cui si deduce che la somma x2+y2 deve restare costante
ed uguale al valore iniziale (in questo caso x(0)2+y(0)2=1).
Poter valutare nel corso della integrazione numerica se e
come varia la costante del moto consente di misurare
direttamente l’errore numerico che si commette.
Metodi numerici generali
• Gli argomenti presentati oggi servono per costruire
un linguaggio didattico comune.
• Dovrebbero essere tutti a voi completamente noti e
serviranno nelle prossime lezioni per costruire
strumenti di lavoro più sofisticati ed efficaci.
• Ma dovrebbero altresì essere utilizzabili per
verificare e per scegliere metodi risolutivi adeguati
ai problemi che si porranno alla vostra attenzione.
DA PORTARE ALL’ESAME
Per casa:
• Organizzare il calcolo della funzione esponenziale
per passi successivi, impiegando sviluppi in serie
di Taylor a 4, a 6 ed a 10 termini,
• Con passi temporali di ampiezza prossima ad R
fino a 4R , poi con passi di R/2 e di R/8
• Illustrare brevemente i risultati in termini di
precisione e rapidità di calcolo.
• Organizzare similmente il calcolo di seni e coseni
fino a raggiungere 4 π e 8 π .
Oscillatore armonico 1
• L’esempio sviluppato per seni e coseni si trasferisce
immediatamente alla descrizione del sistema
dinamico chiamato “oscillatore armonico”
• La realizzazione fisica elementare è rappresentata da
un punto materiale dotato di massa, soggetto a forza
centrale elastica
• Inizialmente lo immaginiamo vincolato su una retta.
Si potrà constatare che il moto deve essere piano e
che opportune condizioni iniziali consentono il moto
rettilineo.
Oscillatore Armonico 2
• L’equazione del moto può essere scritta :
m d2x/dt2 = - k x
• Con m massa, k costante elastica, origine delle
coordinate nel centro della forza elastica e
coordinate iniziali compatibili con il moto rettilineo
che si intende descrivere.
Velocità; equazioni del primo
ordine.
• Introduco la variabile v =dx/dt
• Il sistema di O.D.E. diviene
• m dx/dt = m v
• m dv/dt = -k x
• Una prima trasformazione usualmente operata è
dividere entrambe le equazioni per il fattore m
(massa: costante positiva).
Normalizzazione
• Si ottiene così:
dx/dt = v
dv/dt = - k/m x
Anche k/m è una costante positiva e ciò permette di
introdurre il nuovo parametro ω con ω2 = k/m.
La seconda equazione diviene: dv/dt= -ω2x. Se
cambio l’unità di misura dei tempi tale equazione
diviene dv/dt= - x dove in luogo dei secondi
utilizzo come unità dei tempi la “pulsazione” ω.
Condizioni iniziali
• L’ equazione dunque può scriversi :
d2x/dt2 = -x
• Oppure
dx/dt = v
dv/dt = -x
Alla equazione vanno aggiunte le condizioni iniziali,
che consentono di selezionare una soluzione
particolare nella famiglia (a due parametri) di tutte
le soluzioni possibili.
Serie di Taylor per l’oscillatore
• Scelgo come come condizioni iniziali i valori di x e
di v = dx/dt al tempo iniziale, t0 e li chiamo x0 e v0 .
Posso subito esprimere la soluzione nella forma
dello sviluppo in serie di Taylor.
• x = x0 + dx/dt (t- t0) + ½ d2x/dt2 (t- t0)2 +
+ (½ *1/3) d3x/dt3 (t- t0)3 + . . .
con le successive derivate espresse dalla equazione
originale: d2x/dt2 = dv/dt = -x;
d3x/dt3=- dx/dt =-v e calcolate nel tempo t0
Oscillatore, seni, coseni.
• Si ottiene :x(t) = x0 + v0 (t- t0) - ½ x0 (t- t0)2 +
-(1/6) v0(t- t0)3 + 1/24 . x0 (t- t0)4 + . . .
Cioè x(t) = x0 ( 1- (t- t0)2 (1/2-1/12 (t- t0)2 (….)))
+ v0 (t- t0) ( 1- 1/6 (t- t0)2 (…)))
• Si riconoscono nelle espressioni che moltiplicano x0
e v0 le espressioni di Taylor per cos(ωt) e sin(ωt),
con ω = 1 per la particolare scelta effettuata per
l’unità di misura del tempo.
• Si noti che abbiamo scelto anziché la forma di
somma di termini (rapidamente) decrescenti, la
forma di prodotto di fattori sempre più prossimi ad
1, per ragioni di precisione numerica complessiva.
DA PORTARE ALL’ESAME
Il programma richiesto
• Le espressioni sopra individuate permettono di
valutare x(t) e v(t) per successivi valori di t. Si
suggerisce di scegliere l’intervallo 0 <= t <=10 ed in
tale intervallo 100 punti.
• Si costruisca un file testo con 101 record di terne
t,x(t),v(t) e, mediante il software gnuplot, si
rappresentino su un grafico le curve x(t) ed y(t)
sovrapposte e su un altro grafico la curva ottenuta
unendo tutte le coppie (x,t) di ciascun tempo
successivo.
Metodo di Eulero
A partire dal punto P(x(0),v(0)) = P, si valuti il vettore
ΔP = (dx/dt*Δt, dv/dt*Δt), mediante le espressioni
dx/dt = v ; dv/dt = - x Si valuti ora il nuovo punto P1 =
P(0+Δt) = P0+ ΔP. Il vettore ΔP è ortogonale ad OP ( il
prodotto scalare dei due è x dx+v dv = xv – vx = 0 ) ed
è tangente alla traiettoria. Ripetere la costruzione
ottenendo P2 = P(0 + 2Δt).
Nel programma (programma 5) si registrino sul file
uscita le successive coppie di coordinate relative ai
successivi valori t+ nΔt e si tracci la traiettoria così
ottenuta, confrontandola con quella ottenuta mediante
lo sviluppo di Taylor.
Eulero con correttore
• Nel metodo di Eulero l’arco di traiettoria percorso
nell’intervallo t — t + Δt viene approssimato con un
segmento rettilineo tangente alla traiettoria in t .
È una approssimazione del primo ordine e l’errore
vale circa derivata seconda per dt2 .
• Una approssimazione migliore si otterrebbe
sostituendo al segmento di retta un arco di parabola,
che oltre a partire da P(t) con tangente dP(t),
giungesse a P(t+Δt ) con tangente dP(t+Δt).
• Tale risultato potrebbe essere ottenuto se si potesse
calcolare la tangente nel punto P(t+1/2Δt)
Predictor-corrector
• Proviamo ad utilizzare come valore approssimato
della tangente nel punto P(t+Δt) il valore dato dalla
derivata calcolata in P+ΔP. Si valuti un nuovo valore
di ΔP come valor medio di quello calcolato al tempo
t e di quello valutato nella approssimazione al primo
ordine. L’errore che si commette ora è del terzo
ordine e consente di ricalcolare la traiettoria con
errore ancora minore.
• Il metodo complessivo può essere definito un
“predictor-corrector-evaluation-corrector”,
esemplare di metodi più elaborati che si studieranno
in seguito.
Implementazione
•
•
•
•
•
•
•
•
•
•
•
STAMPA DI TITOLO ED AUTORE
Float tempo,tmax,delt,x,v,u,w,xa,va
Integer i,j,k,npunti,ndiv,nru
Leggi tmax, npunti,ndiv,nru
delt=0.5*(tmax/npunti)/ndiv
Leggi x0,v0 ;Tempo=0 ; x=x0; v=v0
(Oscar: ) dx=v;dv=-x
For i=1 to npunti ;Registra tempo,x,v,dx,dv
For j=1 to ndiv ;xa=x+dx*delt;va=v+dv*delt;
For k=1 to nru ;x=xa+dx*delt;v=va+dv*delt
(Oscar:) dx=v;dv=-x;
Next k; Next j; Next i;
FINE
Suggerimenti per l’esecuzione
• xa e va sono variabili ausiliarie. Si userà la
medesima struttura di programma anche per un altro
metodo.
• Si scelga inizialmente tmax =10, npunti =100, ndiv
=1, 2, 4, 8, 16.
• Nru =1 realizza il metodo di Eulero; = 2 Eulero con
correttore; =3 con ulteriore correzione
• Si chiede di confrontare le traiettorie ottenute con
nru=1 e con ndiv diversi e le traiettorie con ndiv=1
ed nru diversi e scegliere e motivare la scelta
combinata preferibile.
Serie di Taylor (complementi)
Si chiede di trovare sull’intervallo [to, to+T] la soluzione
dell’equazione differenziale: dy/dt = G ( y , t ) con la
condizione iniziale y(to) = yo. G è analitica nel punto
(yo,to). Vale lo sviluppo y(t)= (1/i!) y(i) (to) (t-to)i
Se [t-to] è più grande del raggio di convergenza della
serie , l’errore di troncamento non tende verso zero al
crescere del numero dei termini. Ma si può
suddividere l’intervallo [to to+T] in una catena di
intervalli parziali definiti dai successivi punti tj .
Supponendo ora di avere già trovato nel punto tj la y
e tutte le derivate necessarie, posso scrivere
l’approssimazione di Taylor per t compreso tra tj e tj+1 :
y(t)  zj(t) = (1/i!)y(i)(t-to)i .
Formule di integrazione del
2do ordine che usano le sole derivate prime
Nel caso si intenda arrestare lo sviluppo di Taylor al termine di
grado 1, si ottiene yj+1 =zj(tj+1) = yj + h G(yj,tj). È la formula che
descrive il “metodo di Eulero”, impiegato per la dimostrazione
costruttiva della esistenza ed unicità della soluzione y(t) al
problema di Cauchy proposto.
Considero l’uguaglianza y(t+h)=y(t)+G(y,t+) d

Se l’integrale è sostituito dal prodotto h·G(y(t),t), l’errore è O(h2)
.Una approssimazione dell’ordine di O(h3) si ottiene valutando
l’integrale con la formula dei Trapezi: ½ h (G(y(t),t) + G(y(t+h),t+h)
Ma non è ancora noto il valore di y al tempo t+h.Ne conosco però
una approssimazione: y(t+h)= y(t)+h G(y,t) +O(h2)
L’integrale può essere ottenuto come prodotto di h per G(y,t) con
y valutato in un opportuno punto interno all’intervallo h di
integrazione.
Errore O(h2) ed O(h3)
• Una prima approssimazione di y è:
y* = yi + h G(yi,ti) = yi+1 + O(h2)
tale approssimazione consente di valutare :
yi+1= yi + ½ h (G(yi,ti)+G(y*,ti+1)). Che errore si
commette utilizzando y* anziché yi+1 ?
• L’errore è espresso da:
½ h (G(yi+1,ti)-G(y*,ti+1))= ½ hGy(yM,ti)(y*-yi+1))
dove la derivata parziale di G va valutata in un
opportuno yM , intermedio tra y* e yi+1.. L’errore
quindi è O(h3).
Altre formule di 2do ordine
• Un’altra formula al 2do ordine è:
y* = yi + ½ h G(yi,ti)
yi+1 =yi + h(Gy*,ti + ½ h)
• Nei casi in cui le derivate terze sono piccole,
una formula del secondo ordine con errore
minimo è data da:
• k1=h G(y,t) ; k2=h G(y+2/3 k1,t+2/3 h) ;
Δy= ¼ ( k1+3 k2)
Formule del 4to ordine
• È molto diffusa la formula del 4to ordine:
k1=hG(y,t);
k2=hG(y+ ½ k1,t+ ½ h);
k3=h(G(y+ ½ k2, t+ ½ h);
k4=h(G(y+k3, t+h);
yi+1= yi + 1/6 (k1+2k2+2k3+k4)
Formule parametriche
(complementi)
• Tale formula del 4to ordine è dedotta per s=1 da
una formula nel parametro libero s:
k1=hG(y,t);
k2=hG(y+ ½ k1,t+ ½ h);
k3=h(G(y+ ½k1+1/(2s)(k2-k1), t+ ½ h);
k4=h(G(y+k2 + s (k3-k2), t+h);
yi+1= yi + 1/6 (k1+4k2+2s(k3-k2)+k4)
Altre formule (complementi)
La formula Runge – Kutta – Gill che
suggeriamo impiega il valore s=(0.5)
e la vedremo esplicitamente più avanti.
Un’altra formula al 4to ordine è:
k1=hG(y,t);
k2=hG(y+ 1/3 k1,t+ 1/3 h);
k3=h(G(y+k2-1/3 k1, t+2/3 h);
k4=h(G(y+k3+k1-k2, t+h);
yi+1= yi + 1/8 (k1+3k2+3k3+k4).
Metodi di Runge – Kutta di
ordine > 4 (complementi)
Negli esempi visti si sono costruiti i vettori ki di incrementi
attribuibili ad y, per valutare diverse determinazioni di G(y,t).
I metodi visti si generalizzano proponendo i vettor α e p e la
matrice β costanti, di dimensione q e q x q-1, mediante i quali si
valuterà il vettore k degli incrementi da assegnare ad y, per
calcolare G in punti diversi del dominio, in modo simile a come
si sono calcolate le k1 e k2 nei metodi del secondo ordine.
Posto qm=q-1
k1= h G(y,t);
k2= h G(y+b21 k1, t+a2 h);
k3= h G(y+b31 k1 + b32 k2 , t+a3 h );
kq= h G(y+bq1 k1+…+ b q qm k qm , t+aq h).
y(t+h)  z(h) = y(t)+
i pi ki
Metodi di Runge – Kutta
espliciti (complementi)
La matrice beta sarà triangolare inferiore, giacché si valutano i
successivi elementi di k solo mediante gli elementi già calcolati.
(è questo il significato dell`aggettivo “espliciti “).
Le componenti di alfa, p e beta vanno scelte in modo da
realizzare la precisione prevista.
A tal fine si definisce le scarto E(h) tra la soluzione cercata y(t+h)
ed il valore fornito dall’algoritmo, z(h). Tale E può essere
sviluppata in serie di Taylor e, per metodi di ordine s, E deve
annullarsi ad h=0 con le sue prime s derivate
Per metodi del primo ordine (q=1) si ha z(h)=y + p1 h G(y,t);
E(h) = y(t+h) - z(h) = y(t+h)-y(t)- p1h G(y,t). Si può scrivere la
serie di Taylor per E in funzione di h e si vede che p1 deve
valere 1.
Il piano delle fasi
Le immagini dei “ritratti in Fase” sono tratte dagli appunti delle lezioni di G.Benettin per gli
allievi ingegneri dell’Informazione
• Uno strumento prezioso per lo studio di problemi
della Fisica Matematica è l’uso del “piano delle
fasi”, cioè la rappresentazione delle curve
‘posizione-velocità’ nel piano cartesiano x-y. Nel
caso dell’oscillatore armonico tali curve sono ellissi
concentriche, o, con scelta opportuna della scala,
cerchi, con il raggio determinato dal valore
dell’invariante x2+v2 . Le condizioni iniziali con x0
=0 sono rappresentate da punti sulla retta verticale,
quelle con v0 = 0 sulla retta orizzontale.
Alcuni ritratti in fase
Verso di percorrenza dell’orbita
Nei sistemi meccanici il verso con cui il punto rappresentativo percorre le orbite è orario: un punto in x=0 con v>0
parte dal ramo positivo dell’asse verticale ed invade il primo
quadrante man mano che x aumenta (perché v è >0) e v
diminuisce (allorché dv/dt è ~-x, quindi < 0). Poi quando
v=0, x raggiunge il massimo positivo ed il punto entra nel
quarto quadrante, con x >0, decrescente e v che diviene <0.
Nella immagine sono illustrate le ellissi dell’oscillatore
armonico, così percorse, le iperboli del repulsore armonico,
per le quali vale la medesima osservazione, le rette della
particella libera, percorse da sinistra a destra per v >0 e da
destra verso sinistra per v<0. L’asse delle x (v=0) non è una
traiettoria ma un insieme di punti di equilibrio isolati.
Moto per gravità
• L’equazione è d2x/dt2 = -mg (con g = 9.81)
• Le orbite sono parabole ad asse orizzontale.
• Presentare come esercizio un esempio di moto per
gravità del punto libero, con condizioni iniziali
diverse ( p. es. velocità iniziali positive di valori
diversi e posizioni iniziali positive con v=0 )
Equazioni differenziali ordinarie
non lineari
• Pendolo, crescita con risorse limitate, moto
kepleriano, bipendolo
• Equazioni con ciclo limite: Van der Pol,
Lotka Volterra
• Equazioni con esiti caotici: pendolo forzato,
mappa standard
Moto del pendolo
• Si distinguono tre tipi di moto assai diversi:
• Le librazioni, od oscillazioni, rappresentate da
ellissi, (con centro in (0, 0) ed in (2nπ , 0), n intero
positivo o negativo), ellissi che assomigliano sempre
più a cerchi per oscillazioni di piccola ampiezza.
• Una curva limite, rappresentata da v = 0, x = ± π,
punto isolato in equilibrio instabile, e da due
traiettorie che partono e terminano quanto si vuole
vicino al punto isolato, passano per x = 0, v = ± 2 e
vengono percorse in tempo infinito!!!
• Le rotazioni, orarie p.es. passanti per x = 0, v >2 ed
antiorarie, p.es. passanti per x = 0, v < –2
Ritratto in fase: oscillatore
armonico,buca di potenziale
Buca con colle, doppia buca
Ritratto in fase per il pendolo,
relazione con il potenziale
Tempi di percorrenza delle
traiettorie
• Mediante i programmi di integrazione definita già
noti dal corso di Analisi Numerica si possono
valutare i tempi di percorrenza delle traiettorie .
• Si suggerisce di valutare il tempo necessario per il
pendolo per compiere l’oscillazione completa
(periodo) in funzione della semiampiezza o della
energia del pendolo, p.es. graficando il periodo T per
oscillazioni di 2 gradi, 5° , 15° 30° 90° e fino a 179
gradi e valutando v dalla espressione v2+2mgh=cost.
e T dalla T=4  0 xmax ( 1/v)dx
Periodo del pendolo
Struttura dell’algoritmo per
grafici nel piano delle fasi
• L’algoritmo per il tracciamento dei grafici nel piano delle
fasi prevede le routines fondamentali:
1)Scelta del problema differenziale da risolvere e delle scale
grafiche di rappresentazione
2)Scelta del metodo di integrazione numerica e della
tolleranza accettabile per gli eventuali integrali primi
3)Assegnazione dei parametri del calcolo: tempo iniziale e
finale dell’integrazione, numero di punti da rappresentare
nell’intervalli, n di ulteriori suddivisioni per il raffinamento
della precisione.
4)Assegnazione dei valori iniziali
5)Ripetizione del calcolo con altri valori iniziali
6)Ripetizione con altri parametri o con altri metodi.
Struttura dei dati
• Un primo blocco di variabili contiene gli estremi
dell’intervallo di integrazione, l’ampiezza del passo di
calcolo e gli estremi dei cicli di ripetizione dei passi
elementari e delle rappresentazioni dei successivi punti
• Il blocco delle variabili risultato della integrazione conterrà
per esempio il tempo t, le variabili x e y, le velocità u e v ,
assieme alle relative derivate ed a variabili ausiliarie da
impiegare nel Runge Kutta. Le variabili possono essere
organizzate in vettori: vettore y, di componenti t,x,v ed
indice da 0 a 2, vettore r , di componenti dt/dt=1, dx/dt =v,
dv/dt e vettore ya contenente le variabili ausiliarie. Oppure
con indici da 0 a 4 per il bipendolo e per il problema di
Keplero in coord. rettangolari.
Il passo base di calcolo
• Per ogni problema va assegnata una
particolare funzione vettoriale che assegna le
derivate delle funzioni incognite valutate al
tempo attuale
• Per l’oscillatore armonico: dx/dt= -y;dy/dt=x;
• Per il pendolo:
dx/dt=v;dv/dt=-sin(x)
• Per il problema di Lotka-Volterra:
dx/dt = α x - βx y
dy/dt = - γ y +δ x y
Oscillatore armonico smorzato e
amplificato
L’equazione per l’oscillatore armonico con effetti di
smorzamento o di amplificazione è:
d2x/dt2 + β dx/dt +x =0;
dove se β è una costante positiva si descrive un moto
smorzato, se β è una costante negativa si descrive un
moto amplificato (provare a graficare le mappe con
valori diversi di β )usando la consueta
trasformazione in due variabili.
Oscillatore di Van der Pol:
• L’equazione di Van der Pol è:
d2x/dt2 + µ(x) dx/dt +x =0
• Dove µ(x) è una funzione che assume valori negativi
per x “piccolo” e valori positivi per x “grande”,
descrivendo il comportamento di alcuni dispositivi
elettronici impiegati nei circuiti amplificatori ed
esemplarmente la valvola “triodo”, L’equazione
descrive il comportamento di oscillatori elettronici
disegnati per produrre onde di varie forme ed in
particolare onde “quadre”.
Van der Pol
• Rappresentiamo la funzione µ(x) con l’espressione
µ0( x2 /a2– 1) che è negativa se x è compreso tra -a e
+a e positiva fuori di tale intervallo. L’equazione di
Van der Pol è : d2x/dt2+β ( x2– 1) dx/dt + x = 0,
ottenuta con un fattore di scale nelle x, evidenziando
β, parametro adimensionale significativo.
• Si può scrivere: d/dt ( dx/dt + β ( x3/3 –x ) ) =-x,
cioè dy/dt = -x/ β, posto β y= dx/dt + β ( x3/3 –x) ,
da scriversi come: dx/dt = β (y- γ (x)). L’equazione
di secondo grado è così trasformata in una coppia di
equazioni di primo grado. Si è posto: γ(x)= x3/3 –x.
Come d’uso faremo attenzione al caso β >> 1.
Van der Pol osservazioni
• Sulla cubica γ, dx/dt è = 0 e la traiettoria ha tangente
verticale, con dy/dt di segno opposto a x; fuori dalla
cubica γ, la dx/dt non è più zero, ma è fortemente
condizionata dalla presenza del parametro β >> 1.
Lontano dalla cubica la traiettoria è pressoché
orizzontale, essendo dx/dt moltiplicato per β. In
prossimità della cubica dy/dt diviene rilevante,
negativo se x>0, e la traiettoria si piega. Discesa
sotto la γ, la traiettoria scende, perché dy/dt è
negativo, ma punta verso la γ, perché anche dx/dt è
negativo, senza tuttavia mai poterla attraversare
perché la pendenza dx/dt, che è negativa, diminuisce
ancor più in vicinanza di γ.
Van der Pol
Crescita con risorse limitate
• Si desidera talora descrivere l’evoluzione di
popolazioni che risentono di effetti di
sovraffollamento. Tali effetti possono essere descritti
dall’equazione: dx/dt = α x -ζ x2. Esistono due
“punti critici”, valori di x per i quali la crescita è
zero: x=0 ed x= α/ζ. Per tali valori di x il grafico ha
tangente orizzontale. La pendenza del grafico è
massima per x= α / (2ζ).
• A tale valore conviene attestare la dimensione della
popolazione se si intende utilizzare economicamente
la produzione e massimizzare il prelievo, che potrà
essere uguale alla crescita corrispondente a tale
dimensione, cioè α2 / (2ζ).
Problema di Lotka-Volterra
• Si intende descrivere l’interazione tra due popolazioni
viventi: le prede x ed i predatori y. Le prede dispongono di
risorse alimentari illimitate e si riproducono con tasso α.
I predatori, in assenza di prede, muoiono con tasso γ.
In assenza di predatori le prede crescono in modo
esponenziale: dx/dt =αx; in assenza di prede i predatori
decadono: dy/dt =- γy. L’interazione tra le due popolazioni
è descritta dai parameri β e δ : nell’unità di tempo in
presenza di y  0 la x diminuisce per il termine –β xy
mentre in presenza di x  0 la y aumenta del termine δ xy
Le equazioni sono dunque:
dx/dt = α x – β x y ;
dy/dt = - γ y + δ x y
Le equazioni di Lotka-Volterra
• dx/dt = α x – β x y ;
dy/dt = - γ y + δ x y
• Le derivate si annullano nei “punti critici” o “di
stazionarietà” : (x=0,y=0) e (y= α/ β , x= γ/δ)
• Sono l’origine degli assi ed un punto, intersezione di due
rette, orizzontale e verticale, tagliate ortogonalmente da
tutte le traiettorie
• In certe condizioni le prede possono risentire effetti di
“affollamento”, che descriviamo con il termine – ζ x2 e che
produce un altro punto critico in (y=0,x= α / ζ ).
• dx/dt = α x – β x y – ζ x2
• dy/dt = - γ y + δ x y
Problema di Lotka-Volterra
Lotka Volterra modificato dal
termine di affollamento
Altri problemi
• Pendolo smorzato
Pendolo forzato
bipendolo
• Bipendolo: Piano verticale, origine fissa O e due punti, Q e
P, con O Q =const,QP = const; d(OQ )/dt = v, d(QP)/dt
=w, quindi d(OP)/dt =v + w. Le equazioni cardinali vanno
proiettate sulle direzioni di v e di w, così da annullare
l’effetto delle reazioni vincolari, che restano non calcolate
fino alla determinazione delle leggi del moto.Energia
cinetica e potenziale hanno le espressioni:
• W= ½ mQ v2 + ½ mP(v+w)2,
V= g(mQ hQ + mP hP),
W= ½(mQ +mP) v2 + ½mP (w2 + 2 v w) =
½ (mQ +mP) θ’ 2 + ½ mP (φ’2 +2 θ’ φ’ cos (θ-φ)).
Posti ad uno i moduli OQ e QP ed indicando con θ’, φ’
le derivate temporali delle coordinate angolari.
• h Q= - cos θ; hP = hQ - cos φ;
V= - g( (m + m ) cos θ+ m cos φ)
Grafica bipendolo
Equazioni di Lagrange
• OQ=OQ(θ); QP=QP(φ);
vQ=d(OQ)/dt= (OQ)/ θ dθ/dt = d(OQ)/ θ θ’
• vP =d(OP)/dt= (OP)/ θ θ’ + (OP)/ dφ φ’
• Nota che (OP)/ φ è uguale a vP/ φ’
• Si cerca una espressione per aP•(OP)/θ .
Si valuti vP • (OP)/θ che si può calcolare come
vP •  vP/  θ’ e quindi ½ (v2P) /  θ’.
• Lagrange suggerisce di valutare aP • (OP)/θ come
d/dt((v2P/2)/θ’)- (v2P/2)/  θ
Equazioni di Lagrange
• Le Equazioni di Lagrange possono essere ottenute
come equazioni cardinali proiettate sul sistema di
vettori v e w, ortogonali alle reazioni vincolari, che
definiscono gli “spostamenti virtuali”:
• d/dt (W/(θ’) - W/  θ = -  V/ d θ
• d/dt (W/(φ’) - W/  φ = -  V/  φ
• (mP + mq) d2θ/dt2 + mq cos (θ-φ) d2φ /dt2 =
= -g (mP + mq) sin θ+ mq (d φ / dt) 2 sin (θ – φ)
• mq cos (θ-φ) d2θ/dt2 +
mq d2φ /dt2 =
= -g mq sin φ - mq (d θ /dt) 2 sin (θ – φ)
Forma matriciale
•  (mP +mq) , mqcos(θ-φ)  d2 θ/dt2 
|
||
|=
•  mqcos(θ-φ)
, mq d2 φ /dt2
•
=
•
- g(mP +mq)sin θ + mq (d φ /dt)2 sin (θ – φ) 
|
|
 -g mq sin φ
- mq (d θ /dt)2 sin (θ – φ) 
Soluzione (regola di Cramer)
Il sistema espresso dalla relazione matriciale:
A B x =  C
D E y =  F
si risolve con le posizioni:
den=AE-BD; numx=CE-BF;numy=AF-CD;
e quindi:
x=numx/den; y=numy/den
Soluzioni particolari
• Una volta individuate le formule per esprimere le
derivate prime delle velocità, il sistema di derivate può
essere scritto all’interno del ciclo di integrazione.
• Vanno assegnate le condizioni iniziali su posizioni e
velocità.
• Poste le velocità iniziali uguali a zero suggerisco due
tipi di condizioni iniziali per le posizioni:
• 1) la sola  oppure la sola  diverse da zero; p.es.
uguale ad 1 o a valori diversi più o meno piccoli
• 2)  = 2  e  = -  2  e  con valori diversi più o
meno piccoli
Cond.
iniz.
Tipo 1
Modi di oscillazione
Le condizioni iniziali di tipo1 presentano inizialmente
eccitata solo una delle coordinate, ma nel moto la
reazione vincolare mobilizza immediatamente anche
l’altra coordinata e le due coordinate evolvono in
modo apparentemente indipendente e scoordinato.
Nel piano [ , ] le traiettorie riempiono con “figure
di Lissajous” quadrilateri centrati in O, con assi propri
obliqui agli assi coordinati e poggiati sul valore della
variabile assegnata.
Cond.
iniz.
Tipo2
Modi normali (ortogonali)
Le condizioni iniziali di tipo 2 evolvono in modo da
mantenere la proporzione iniziale tra le due variabili,
che inoltre si muovono con la medesima frequenza.
Nel piano [ , ] le traiettorie occupano sottili rettangoli
centrati in O ed orientati nel 1° e 3° oppure nel 2° e
4° quadrante.
Le condizioni di tipo 2 eccitano i “modi normali”,
indipendenti. Le condizioni generiche, e quindi anche
le 1, vanno lette come combinazioni di moti tipo 2.
AUTOVETTORI
Il sistema A B x =  λ x
D E y =  λ y
determina particolari vettori (x,y) che, moltiplicati
dalla matrice, sono amplificati o ridotti dal fattore λ,
ma si mantengono della medesima direzione ed
hanno il nome di “autovettori”. Sono possibili due
determinazioni di λ , che prendono il nome di
“autovalori” . I vettori (x,y ) corrispondenti sono gli
“autovettori” associati a ciascun λ .
AUTOVALORI
La soluzione della equazione matriciale richiede che
det (A-λ) B  =0
D (E- λ) 
Cioè ((A-λ)(E-λ) –BD) = 0; λ2 – (A+E) λ +AE –BD =0
( (A+E) è la “traccia”, AE –BD è il “determinante”)
Ciascun autovalore permette di calcolare una possibile
determinazione del quoziente x/y o y/x) e quindi un
autovettore , per il quale è indeterminata la attuale
lunghezza.
Le “condizioni iniziali di tipo 2” corrispondono ad
autovettori della matrice del sistema lineare proposto.
Forze elastiche su punti allineati
• Si considerino ora i punti O, Q, P allineati. O è fisso,
Q ha la coordinata x, P ha la coordinata y entrambe
relative ad O. Su Q agisce una forza elastica k1 di
centro O ed una forza elastica k2 agisce tra Q e P.
• Si scrivano le equazioni del moto, per le due
funzioni incognite x ed y .
Forze elastiche su punti allineati (gr)
Due punti con forze elastiche.
•
•
•
•
•
FQ = -k1 x + k2 (y-x)
FP = - k2 (y-x)
mQ d2x/dt2 = - (k1 + k2 ) x + k2 y
mP d2y/dt2 =
k 2 x - k2 y
Poste uguali ad 1 , per semplicità, le masse e le costanti
elastiche, l’equazione agli autovalori è (–2 –λ)(–1–λ) –1 = 0
Cioè λ2 + 3λ + 1 = 0, con soluzioni (-3 ± sqrt(5))/2,
-0.381966 e -2.618 . Introdotti tali valori nella equazione
matriciale si ottiene x = 0.618 y e x = - 1.618 y.
Questi sono gli autovettori che danno luogo ai “modi
normali” di oscillazione.
Moto Kepleriano
In coordinate rettangolari x ed y:
• r2 = x2 + y2 ; r3 = r2 sqrt(r2)
• dx/dt=v ; dv/dt=-Kx/(r3);
• dy/dt=w; dw/dt=-Ky/(r3)
In coordinate centrali r e θ : posto OP = ru, u è il versore
radiale e v ad esso ortogonale è il versore trasverso.
La velocità d(ru)/dt ha componenti
dr/dt u+ r(du/dt) = dr/dt u+ r (dθ/dt) v
e l’accelerazione
(d2 r/dt2 - r (d θ /dt) 2) u + (r (d2θ /dt2) +2 (dr/dt) (dθ/dt)) v=
(d2 r/dt2 - r (d θ /dt) 2) u + 1/r (d(r 2(dθ/dt)) /dt) v.
Moto Kepleriano 2
Conforme la legge di Newton l’accelerazione è radiale:
a= (-K/r2)u . L’annullarsi di d(r 2(dθ/dt)) /dt
conduce all’integrale primo r2ω = cost; con ω = dθ/dt , e
dθ/dt = cost/r2 ; dr/dt = dr/dθ dθ/dt = - cost d(1/r)/dθ.
d2r/dt2 = dθ/dt d(dr/dt) / dθ = - cost dθ/dt d2 (1/r)/dθ2
= - cost2 (1/r2) d2 (1/r) / dθ2 = - K/r2 .
Appare come se la variabile naturale fosse anziché r la y=1/r.
si può infatti scrivere
-cost2(1/r2) +d2 (1/r)/dθ2 = y2 +d2 y/dθ2 = cost2 /K, che è
analoga alla equazione del pendolo, con un termine non
omogeneo.
Moti centrali coord,radiali
Moti caotici
• Moti caotici
Mappa standard :
• (x,v) -> Φ (x,v) ->(x+v,v+εsin(x))
il comportamento caotico è regolato dal
valore del piccolo parametro ε.
Moti caotici
Runge Kutta del secondo
ordine
• Un Runge Kutta del secondo ordine nella
incognita (vettoriale) y può essere
programmato prevedendo la variabile
(vettoriale) ausiliaria ya ( e l’ulteriore variabile
yo destinata a contenere i valori iniziali) . La
variabile (vettoriale) r contiene, dopo
ciascuna invocazione della routine “derivate”,
i valori delle derivate delle componenti di y,
calcolate dai valori delle y attuali.
Metodo di Heun
• Si valuti dunque r, derivata di y; si valuti un
incremento parziale ½ h r e si memorizzi ya=y
+ ½ h r. Un valore provvisorio per y al tempo
t+h si ottiene sommando ancora l’incremento
parziale: y=ya+ ½ h r. Tale valore è quello che
si otterrebbe con il metodo di Eulero (preciso
entro O(h2)). Si invoca nuovamente la routine
per il calcolo delle derivate e si ottengono i
nuovi valori delle componenti di r che offrono un
valore di y preciso entro O(h3) y=ya+ ½ h r.
Una ulteriore correzione si ottiene, una volta
valutata r nel più preciso valore più recente di y,
ripetendo ancora il calcolo y=ya+ ½ h r.
Precisione del 4to ordine
punti da rappresentare e passi di integrazione
• L’intervallo complessivo di calcolo 0-T va
suddiviso in intervalli minori, tra uno e l’altro dei
punti da rappresentare in uscita , e quindi, se
necessario, in sottointervalli che consentano di
restare in un conveniente raggio di convergenza per gli sviluppi formali. L’effettivo passo elementare di calcolo sarà caratterizzato dall’incremento temporale : step = T/(npunti*ndiv)
Precisione del
to
4
ordine
cicli interni di calcolo
Tale passo elementare comprende un ciclo
(passoru) ripetuto quattro volte (indice iru) che
inizia con la chiamata a deriv e contiene un ciclo
interno (incre), da ripetere per ciascuna variabile
del sistema dinamico (indice ivar).
Tale ciclo interno è il nocciolo dell’integratore e
definisce due variabili interne, rr ed aum, cui va
consentita la massima precisione disponibile.
Runge Kutta Gill
• q = (0.5); aa = [1/2,1-q,1+q,1/6] ;
• bb=[2,1,1,2]; cc=[1/2,1-q,1+q,1/2];
• Iniziotraccia ;Valin;
(la routine assegna i valori iniziali e determina il valore logico FINE)
• : Ripeti per ipunti = 1 to npunti
• Ripeti per idiv= 1 to ndiv; Ripeti per iru = 1 to 4; deriv
• Ripeti per ivar = 1 to nvar(questo è il ciclo interno “incre”)
rr=r(ivar)*step;
aum=aa(iru)*(rr-bb(iru)*ya(ivar))
y(ivar)=y(ivar)+aum
ya(ivar) =ya(ivar) +3*aum –cc(iru)*rr
Fineivar ; Fineiru ; Fineidiv ; segnapunto;
Fineipunti; Valin; se non fine vai a iniziotraccia ; FIN
procedure
• Costanti vettoriali : sono aa,bb,cc; “vettori” sono y, che
contiene le variabili del problema, r che contiene le derivate,
ya che contiene variabili ausiliarie, corrispondenti alle y , y0
che contiene i valori iniziali, chiesti dalla procedura valin ed
assegnati dalla medesima alle y , ad ogni inizio di traiettoria.
• Valin : assegna i valori iniziali per il giro successivo e
stabilisce il valore logico della variabile “fine”
• Deriv: assegna le espressioni delle derivate
• Incre: assegna gli incrementi alle variabili
• Passoru: valuta gli incrementi parziali secondo il metodo
Runge-Kutta
• Segnapunto: trasferisce al file output ( o ai vettori da plottare
i valori delle variabili integrate
• Iniziotraccia inizializza una nuova tavola per lo stesso
problema, modificando npunti,ndiv.
Runge Kutta Gill
Le formule Runge – Kutta - Gill realizzano la
precisione del 4to ordine. Si definiscono in
primo luogo la costante q = (0.5) e tre vettori di
costanti: aa = [1/2,1-q,1+q,1/6], bb=[2,1,1,2] e
cc=[1/2,1-q,1+q,1/2].
Tali costanti vengono impiegate nel ciclo interno,
di indice iru che va da 1 a 4 per esplorare le quattro determinazioni di G che realizzano il processo. Ciascuna iterazione del ciclo iru inizia con la
chiamata alla funzione deriv e viene applicata,
con l’indice ivar, a ciascuna variabile dinamica.
Integratori passo passo Runge-Kutta
• L’intervallo complessivo di calcolo 0 - T va
suddiviso in intervalli minori, tra l’uno e l’altro
dei punti da rappresentare in uscita , e quindi,
se necessario, in sottointervalli che consentano
di restare in un conveniente raggio di
convergenza per gli sviluppi formali. L’effettivo
passo elementare di calcolo sarà:
step = T/(npunti*ndiv).
Ciclo interno Runge-Kutta Gill
• Tale passo elementare comprende un ciclo
(passoru) ripetuto quattro volte (indice iru) che
inizia con la chiamata a deriv e contiene un
ciclo interno (incre), da ripetere per ciascuna
variabile del sistema dinamico (indice ivar).
• Il ciclo interno di indice iru è il nocciolo
dell’integratore e definisce due variabili interne,
rr ed aum, cui va consentita la massima
precisione disponibile.
Metodi simplettici
Runge-K,Multipasso,problemi
stiff