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 = - k /m x , detta dell’oscillatore

armonico, descrive il comportamento di un punto materiale soggetto al richiamo di una forza elastica. 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 ottenered2x/dt2= - dy/dt= - x d2y/dt2= dx/dt= - yd3x/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.esx(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 primiRisulta 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 deve restare costante la somma x2+y2, ed uguale al valore iniziale ( in questo caso x(0)2+y(0)2=1

Valutare nel corso della integrazione numerica tale espressione consente di avere una indicazione dell’errore numerico complessivamente commesso.

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’ESAMEPer 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 pigreco e 8 pigreco

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

• mdx/dt=mv

• mdv/dt=-kx

• 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 xAnche 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 = -xAlla 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 le funzioni circolari

• 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

.• 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’ESAMEIl 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 a 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, confron-tandola 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 t j . Supponendo ora di avere già trovato nel punto t j la y e tutte le derivate necessarie, posso scrivere l’approssimazione di Taylor per t compreso tra t j e tj+1 : y(t) zj(t) = (1/i!)y(i)(t-to)i .

Formule di integrazione del 2do ordine

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 fasiLe 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’orbitaNei sistemi meccanici il verso con cui il punto rappresen-

tativo 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 dal punto isolato v=0, x=±π, 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 iniziali5)Ripetizione del calcolo con altri valori iniziali6)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 -1 ed 1 e positiva fuori di tale intervallo. L’equazione diviene :

d2x/dt2 + β ( x2– 1) dx/dt + x = 0dove β è il parametro adimensionale significativo.

Come d’uso poniamo attenzione al caso β >> 1. Le equazioni sono: dx/dt = β (y- γ (x)); con γ(x)= x3/3 -x

dy/dt = - 1/ β x;

Van der Pol osservazioni

• Sulla cubica gamma, dx/dt è = 0 e la traiettoria ha tangente verticale, con dy/dt di segno opposto a x; fori dalla cubica gamma dx/dt non è = 0 , ma è fortemente condizionata dalla presenza del para-metro β >> 1. Lontano dalla cubica la traiettoria è pressoché orizzontale, essendo dx/dt moltiplicato per beta. In prossimità della cubica dy/dt diviene rilevante, negativo se x>0, e la traiettoria si piega. Discesa sotto la gamma, la traiettoria scende, perché dy/dt è negativo, ma punta verso la gamma, perché anche dx/dt è negativo, senza tuttavia mai poterla attraversare perché in vicinanza di gamma dx/dt diminuisce .

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ζ).

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 δ xyLe 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: due punti, P e Q, con !OP!=const ed !

OP! = const; d(OP)/dt = v, d(QP)/dt =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.

• W= ½ v2 + ½ mq(v+w)2,V= g(mP hP +mq hq) W= ½(mP +mq) v2 + ½mq (w2 + 2 v w) =½ ((mP +mq) θ’ 2 + mq (φ’2 +θ’φ’cos (θ-φ)), posti ad 1 v e w

• h P=cosθ; hq= hP+cos φ; V= g((mP+mq)cosθ+ mq cos φ)

Equazioni di Lagrange• Le Equazioni di Lagrange sono le 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 + mqcos(θ-φ) d2φ /dt2 == -g(mP+mq) sinθ+ mq (dφ/dt) 2 sin (θ – φ)

• mqcos(θ-φ) 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 = Fsi 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

Modi di oscillazione• Le condizioni iniziali di tipo 1 presentano inizial-

mente 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.

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.

Moto KeplerianoIn 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 vettore radiale e v = d(ru)/dt . Conforme la legge di Newton l’accelerazione è radiale: dv/dt= (-K/r2). Si ha l’integrale primo r2ω = cost; con ω = dθ/dt .

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 convergen-za per gli sviluppi formali. L’effettivo passo ele-mentare di calcolo sarà caratterizzato dall’incre-mento temporale : step = T/(npunti*ndiv)

Precisione del 4to ordinecicli 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 ; FINE

procedure• Costanti strutturate: sono aa,bb,cc• 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 quat-tro determinazioni di G che realizzano il proces-so. 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 rap-presentare in uscita , e quindi, se necessario, in sotto-intervalli 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). 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.

Metodi simplettici

Runge-K,Multipasso,problemi stiff