Meccanica Razionale per Allievi Ingegneri Aerospaziali Esercitazioni numeriche

Meccanica Razionale per Allievi Ingegneri “Aerospaziali” Esercitazioni numeriche 2008 Equazioni Differenziali e “Ritratti in Fase” per alcuni problemi della Meccanica - Uso di Matlab

Obiettivi del Corso Il corso propone agli allievi alcuni strumenti concettuali e fornisce un primo training per lo studio di problemi della Meccanica Razionale con i metodi del calcolo computerizzato, con riferimento a sistemi dinamici didattici. In particolare presenta alcuni integratori per ODE, offre una introduzione alla dinamica di sistemi non lineari e di sistemi dotati di cicli limite e propone agli allievi l’uso critico di metodi numerici per tracciare i ritratti in fase di alcuni dispositivi classici particolarmente rappresentativi. Il corso offre altresì metodi introduttivi per la soluzione di problemi alle derivate parziali (problemi di Poisson, problemi di Neumann, equazioni iperboliche).

Questioni e date: • Primi Elementi del linguaggio (oct. ) 6, 13 • Integratori ODE: Eulero, Heun, Runge. 20 • Valori iniziali e valori al contorno. Precisione: sviluppi di Taylor per sin, cos, esponenziale. Piano delle fasi. Oscillatori, pendolo. Periodo del pendolo. 27 • Ripasso integratori num. Crescita di popolazioni, prelievo, logistica, Lotka, Van Der Pol, Keplero. Due palle. Matrici, sist. lineari. (nov. ) 3 • Pendolo forzato, bipendolo. Eq. di Lagrange. 10 • Modi di oscillazione: autovettori-autovalori 17 • PDE paraboliche e PDE ellittiche. 24

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. • Attenzione: MATLAB consente con grande facilità sia calcoli immediati, sia la costruzione di “m-files” che possono contenere o programmi o funzioni, da chiamare per nome e da eseguire a chiamata.

Primi programmi • Programma uno: stampare un titolo, poi nome e matricola dell’ allievo e l’intestazione ‘Laboratorio di meccanica …’ (si può usare l’ istruzione “disp”) • Programma due: scrivere su un file le coordinate di sette punti equidistanti su una circonferenza (usando le funzioni sin e cos) : • t = 0: 2*pi/7: 2*pi; x=cos (t); y=sin (t) • Le variabili di Matlab sono intrinsecamente vettori, di lunghezza adeguata all’uso dato dal contesto. La prima istruzione è un ciclo “for” implicito che descrive 8 valori, da 0 a 2 pigreco a passi di 1/7 di giro; le successive due istruzioni generano ascisse ed ordinate dei punti richiesti.

Poligoni e stelle • I punti trovati sono vertici di un poligono disegnabile con il comando “ plot(x, y) ”. Il numero di lati può essere dato da input “ npunti = input(‘ nvertici = ‘) ”. • Programma tre : come due, ma ordinando i punti in modo che siano i vertici di una stella, una volta congiunti successivamente con corde opportune. I medesimi risultati potrebbero essere ottenuti mediante l’uso del software “gnuplot” in ambiente Linux o mediante applicazioni Excel in ambiente Window. • Programma 4: come 2 e 3 , ma assegnando da input il numero di vertici ( e/o il numero di giri , anche non chiusi).

Integratori per ODE e per PDE • Esplorando questioni matematiche e applicazioni molto diverse, darò uno strumento generale e facilmente adattabile, come un “cacciavite a punte intercambiabili”. • Faremo “integratori numerici per Equazioni Differenziali Ordinarie (ODE) ”, con valutabili stabilità e precisione. • Un integratore “del primo ordine” diviene, con semplici modifiche successive, di ordine 2 -3 e 4 -5 , adeguato ad uno svariato numero di problemi della ingegneria (e della biologia) e utilizzabile anche per equazioni differenziali PDE di tipo iperbolico e di tipo parabolico • Vedremo anche esempi di PDE di tipo ellittico, parabolico e iperbolico ed un esempio di trattamento di matrici sparse.

Equazioni differenziali ordinarie, metodo di Eulero in linguaggio Matlab. %equazioni : dx/dt=-y; dy/dt=x , %condizioni iniziali x 0=1, y 0=0 disp(‘Laboratorio di Meccanica Razionale 2008’); disp(‘allievo Mario Bianchi, matr 543007’); disp(‘Integratore ODE –Eulero. Moti Armonici’); t=0; tf=10: npti = 100; delt=tf/npti % calcoli di inizio ciclo x 0 =1; y 0 =0; x=x 0 ; y=y 0 %valori iniziali e loro assegnazione for i=1: npti ; xx(i)=x; yy(i)=y; tt(i)=t; %storage variabili deltx=-y*delt; delty=x *delt; % calcolo degli incrementi t=t+delt; x=x+deltx; y=y+delty; % loro assegnazione end; plot (tt, xx), plot(tt, yy), plot(xx, yy)

Euler-ode-matlab Altra implementazione Matlab del corpo programma • • • a=[0, -1; 1, 0]; %matrice per calcolo delle derivate dt=0. 01; tf=10 %assegnazione dati ciclo y=[1; 0]; yy=y % valori iniziali e prima risposta for t=0: dt: tf; r=a*y; %calcolo delle derivate con la matrice a y=y+r*dt %incremento della funzione vettoriale y yy=[yy, y] %storage dellesuccessive risposte end; plot(yy(1, : ), yy(2, : )) %graficazione

Osservazioni: dettagli del linguaggio Matlab • Notare la assegnazione della matrice a, che serve per calcolare il vettore derivate r corrispondente a ciascun valore del vettore y. • Notare la assegnazione dei successivi valori y alla matrice risultati yy. • Notare l’uso di due punti per sintetizzare tutto il range di un indice.

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 P 1 = P(0+Δt) = P 0+ ΔP. Il vettore ΔP è ortogonale ad OP (prodotto scalare dei due: xdx/dt+ vdv/dt = xv – vx = 0) ed è tangente alla traiettoria. Ripetere la costruzione ottenendo P 2 = P(0 + 2Δt). • Nel 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 (vedi più avanti).

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 dt 2. • Una approssimazione migliore si otterrebbe sostituendo al segmento di retta un arco di parabola, che oltre a partire da P(t) con tangente d. P(t), giungesse a P(t+Δt ) con tangente d. P(t+Δt). • Un risultato simile può essere ottenuto se si sa calcolare la tangente nel punto P(t+1/2Δt).

Predictor-corrector • Detto ΔP* l’incremento valutato con la tangente in P, calcoliamo la tangente in P+ΔP*. Il segmento di parabola che cerchiamo descrive un ΔP che è il valor medio tra quello calcolato usando la tangente valutata in P e questo calcolato con la tangente in P+ΔP*. • 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.

Errore nell’integrazione numerica Con il metodo “di Eulero”, a partire da P=(x 0, y 0) viene tracciato un segmento PP’ortogonale ad OP, quindi un successivo segmento ortogonale al nuovo raggio OP+PP’, e il nuovo raggio, OP’, è più lungo del raggio iniziale. Ripetendo l’esecuzione del programma per valori diversi del parametro “npunti”, si constata che l’esito della integrazione dipende dal valore di “npunti” e l’errore diminuisce fortemente aumentando tale valore, cioè riducendo “delt”. Pur essendo l’equazione differenziale soddisfatta nei punti iniziali, e, col nuovo valore R, nei punti di arrivo, R non resta invariante. Le traiettorie appaiono iscritte in spirali che si aprono e non in cerchi di raggio costante. Ciononostante il Metodo di Eulero è uno strumento prezioso anche dal punto di vista teorico, giacché fornisce una dimostrazione “costruttiva” del Teorema di Cauchy sulle Equazioni Differenziali.

Errore per il moto circolare e metodo di ordine più alto. P’ O P p

• Predictor corrector ( RINORMALIZZAZIONE: Una prima terapia può essere “rinormalizzare”le soluzioni, cioè obbligare R a mantenere il valore iniziale, dividendo ad ogni passo x ed y per un fattore opportuno, R. ) • All’inizio del passo la derivata ha la direzione definita dalla equazione, ortogonale ad OP, al termine del passo la nuova direzione della tangente è ortogonale non più ad OP, ma ad OP’=OP+PP’. • Per descrivere un tratto di curva più somigliante alla circonferenza si può optare per un tratto di parabola, di cui siano dati il punto iniziale e due tangenti: all’inizio ed al termine e si ripete il calcolo degli incrementi utilizzando il valor medio delle derivate nei punti iniziale e finale del segmento di traiettoria. L’ordine dell’errore diviene proporzionale a dt 3 : se dt è 0. 1, dt 3 è 0. 001. • Il metodo che si ottiene è detto di Heun o con altra terminologia Predictor Evaluator Corrector Evaluator…: PECE

Pseudo linguaggio per Heun • % ndiv è il numero di suddivisioni intermedie tra due punti successivi del grafico, e diventerà strategico per il metodo delle “sezioni di Poincaré, per eq. dipendenti da t. • • • Float tempo, tmax, delt, x, v, xa, va, rx, rv Integer i, j, k, npti, ndiv, nru Leggi tmax, npti, ndiv, nru tmax=10, npti=50, ndiv=1, nru=3 delt=0. 5*(tmax/npti)/ndiv Leggi x 0, v 0 ; Tempo=0 ; x=x 0; v=v 0 (Oscar: ) rx=v; rv=-x For i=1 to npti; Registra tempo, x, v, rx, rv; For j=1 to ndiv; xa=x+rx*delt; va=v+rv*delt; For k=1 to nru; x=xa+rx*delt; v=va+rv*delt; (Oscar: ) rx=v; rv=-x; Next k; Next j; Next i; FINE

Metodo Heun con successive correzioni Il metodo valuta la posizione dopo l’intervallo di tempo dt/2 e memorizza tale posizione nelle variabili xa, ya; quindi valuta la posizione dopo un successivo intervallo dt/2 (variabili xx, yy), valuta ivi le derivate ed impiega tali derivate (valutate alla fine del passo) per calcolare più corretti incrementi, sempre relativi al secondo semintervallo dt/2, da aggiungere alla posizione memorizzata in xa, ya. La procedura valuta la traiettoria come fosse approssimata da una parabola anziché da un segmento. La correzione può essere ulteriormente migliorata, conforme il valore di nru. L’ulteriore variabile ndiv misura il numero di “passini” da valutare tra successivi punti che saranno effettivamente tracciati nel grafico. Si esplori il metodo per ndiv=1, nru=1 , 2, 3 , 4, poi nru=2 e ndiv=2, 8, 16.

Metodo di Heun in Matlab disp('Laboratorio di Meccanica Razionale 2008'); disp('allievo Mario Bianchi, matr 543007'); disp('Integratore ODE (metodo di Heun) per due variabili'); x 0 =1; y 0 =0; x=x 0 ; y=y 0 ; t=0; tf=10; nru=3; npti = 100; delt=(tf/npti)/2; for i=1: npti ; xx(i)=x; yy(i)=y; tt(i)=t; deltx=-y*delt; delty=x*delt; t=t+delt; xa=x+deltx; ya=y+delty; for j =: nru; x=xa+deltx; y=ya+delty; deltx=-y*delt; delty=x *delt; end ; end; plot (tt, xx), plot(tt, yy), plot(xx, yy)

Heun per più variabili disp('ODE, metodo di Heun per più variabili'); nvar=2 ; y 0 =[0, 1]; ya=[0, 0]; r=[0, 0]; y=y 0 ; t=0; tf=10; nru=3; npunti = 100; delt=(tf/npunti)/2; for i=1: npunti ; xx(i)=x; yy(i)=y; tt(i)=t; r=oscar(y); t=t+delt; ya=y+r*delt; for j =: nru; y=ya+r*delt; r=oscar(y) ; end; Function[r] =oscar(y); r=[0, -1; 1, 0]*y; end;

Heun-ode-matlab • • a=[0, -1; 1, 0]; % matrice per calcolare le derivate di sin cos dt=0. 2; tf=10; nru=2; % parametri del calcolo y=[1; 0]; yy=y; %valori iniziali for t=0: dt: tf; % ciclo nel tempo r=a*y; dy=r*dt/2; ya=y+dy; for iru=1: nru; y=ya+dy; r=a*y; dy=r*dt/2; end; yy=[yy, y]; end; plot(yy(1, : ), yy(2, : ))

Funzioni e parametri per gli Integratori ODE L’integratore ODE chiede di definire una funzione che calcoli le derivate delle incognite da integrare, a partire dai valori delle variabili e del tempo. Vanno inoltre forniti i valori iniziali ed eventuali parametri per controllare la precisione del calcolo IL Matlab offre programmi potenti ed efficienti per l’integrazione di equazioni differenziali con metodi di diverso ordine (2 3 e 4 -5) e per problemi che richiedono particolare attenzione all crescita dell’errore ed alle diverse velocità di evoluzione delle varie componenti (stiffness). I programmi di integrazione che costruiamo in questo corso fanno comprendere il valore ed il significato dei principali criteri di scelta tra vari metodi.

Invarianti • Di fondamentale importanza nel nostro contesto sono certe funzioni, dipendenti dalle variabili e dai parametri del problema, che mantengono costante il loro valore durante il moto, mentre le variabili da cui dipendono assumono successivi valori. Nel caso di sin e cos sappiamo che sin 2(t)+cos 2(t) =1 per qualsiasi valore di t. • Da dx/dt = -y , dy/dt = x , moltiplicando ambo i termini della prima per x ed ambo i termini della seconda per y e sommando si ottiene d(x 2 + y 2)/dt=2( -yx+xy) = 0, cioè la grandezza x 2 + y 2 dovrebbe restare invariante. È infatti R 2, raggio della circonferenza al quadrato. • L’eventualità che nel calcolo che si sta svolgendo R 2 od R non resti costante indica la possibile presenza di errori.

Metodi simplettici Se nell’integratore si impiega per l’incremento di y il nuovo valore della variabile x, anziché il vecchio come fatto or ora, si ottiene un risultato di grande inattesa precisione numerica : xxnuovo=xxvecchio dt*yyvecchio; yynuovo= yyvecchio+dt*xxnuovo . È come aver valutato la velocità xxnuovo al tempo t+dt usando la posizione yyvecchio al tempo t+dt/2, e poi la nuova posizione yynuovo al tempo t+3/2 dt mediante la velocità al tempo t+dt: l’ordine dell’errore va da dt 1 a dt 3. . Si perde la simultaneità dei valori, che è recuperabile poi con opportune interpolazioni tra risultati consecutivi. Il metodo dato dalla formula “preserva il volume nello spazio delle fasi”, realizzando trasformazioni intrinsecamente più precise anche a passi maggiori. Varianti multipasso del metodo consentono precisioni anche molto maggiori.

Predictor corrector polinomiali • Anziché la tangente all’origine, il “predittore” può essere il polinomio che passa per gli ultimi tre o quattro punti già calcolati, il “correttore” sarà un altro polinomio che usa i punti “vecchi”e la tangente calcolata nel punto “predetto”. • Nel punto iniziale il polinomio “predittore” dovrà essere calcolato con una opportuna serie di Taylor oppure valutando i primi punti con il metodo di Runge Kutta che presenteremo.

Metodi di ordine superiore • I metodi qui presentati hanno il grande pregio di impiegare la sola derivata prima nel punto di calcolo, oppure ai due estremi di un intervallo, od anche nel solo punto medio dell’intervallo. Al confronto lo sviluppo in polinomi di Taylor richiede la conoscenza di più derivate successive tutte da valutare all’inizio dell’intervallo. Proponiamo di confrontare le precisioni ottenibili con polinomi di Taylor di vario ordine e con la suddivisione del dominio della soluzione in più sottointervalli.

Integratori passo • 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 mantenere gli sviluppi formali in un conveniente raggio di convergenza. L’effettivo passo elementare di calcolo sarà: step = T/(npunti*ndiv).

Metodi Runge-Kutta: derivate prime in punti diversi anziché derivate di alto ordine nello stesso punto • Il passo elementare comprende un ciclo (passo-ru) 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). • Variante Gill: Il ciclo di indice iru , nocciolo dell’integratore, definisce due variabili interne, rr ed aum, cui va consentita la massima precisione disponibile, anziché i quattro vettori di dimensione nvar del Runge classico. Le derivate prime delle diverse variabili sono calcolate in punti diversi della traiettoria, corrispondenti ad opportuni valori della variabile indipendente t. La formula Runge-Kutta tradizionale sceglie t 1=0, t 2= ½ dt, t 3= ½ dt e t 4= dt.

Derivate prime in punti diversi • Il metodo tradizionale di Runge e Kutta valuta la derivata prima in quattro punti opportuni, poi usa una combinazione dei valori trovati per valutare l'incremento di y r 4 r 3 r 2 r 1 t 0+dt/2

Formule del 4 to ordine • La formula del 4 to ordine tradizionale è: k 1=h. G(y, t); k 2=h. G(y+ ½ k 1, t+ ½ h); k 3=h(G(y+ ½ k 2, t+ ½ h); k 4=h(G(y+k 3, t+h); yi+1= yi + 1/6 (k 1+2 k 2+2 k 3+k 4)

Runge Kutta Gill Le formule della variante Gill realizzano la precisione del 4 to ordine con errore minimo e buona stabilità. Si assegna q = sqrt(0. 5), costante 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.

Metodi Runge – Kutta di ordine > 4 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 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 k 1 e k 2 nei metodi del secondo ordine. Posto qm=q-1 k 1= h G(y, t); k 2= h G(y+β 21 k 1, t+ 2 h); k 3= h G(y+β 31 k 1 + β 32 k 2 , t+ 3 h . . . . kq= h G(y+βq 1 k 1+…+ β q qm k qm , t+ q h). y(t+h) z(h) = y(t)+ i pi ki

Metodi di Runge – Kutta espliciti Nei metodi “Espliciti” 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 “). Per taluni problemi detti problemi “stiff”, con autovalori delle soluzioni molto diversi, con alcune variabili che si comportano come esponenziali decrescenti a tempo di caduta molto più breve di altre, possono essere opportuni metodi “impliciti”come ha ampiamente illustrato lo studioso Gear. Le componenti di alfa, p e beta vanno scelte in modo da realizzare l’ordine di precisione previsto. A tal fine si definisce lo 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 + p 1 h G(y, t); E(h) = y(t+h) - z(h) = y(t+h)-y(t)- p 1 h G(y, t). Si può scrivere la serie di Taylor per E in funzione di h e si vede che p 1 deve valere 1.

Metodi impliciti Per gli incrementi nell'equazione dy/dt=Ay si confrontino le espressioni: y(t+dt)=y(t)+dt * dy(t)/dt; y(t+dt)=y(t)+dt * dy(t+dt)/dt nella seconda la derivata è calcolata nel punto tuttora incognito di ascissa t+dt. A volte l'espressione particolarmente semplice, p. es. lineare, della derivata consente di calcolare direttamente “in forma implicita” il valore y(t+dt) soddisfacente l'equazione agli incrementi finiti. Si verifichi il caso dell'esponenziale decrescente.

Esponenziale decrescente metodo implicito del primo ordine • Da: y(t+dt)=y(t)+dt*A*y(t+dt) si ottiene: y(t+dt)(1 -dt*A) = y(t), cioè: y(t+dt)=y(t)/(1 - dt*A) il metodo può già essere implementato dagli allievi per scrivere una variante opportuna del metodo di Heun. . • Esiste una famiglia di metodi Runge-Kutta impliciti e i. I metodi di Gear per equazioni “stiff” sviluppano tale idea.

Runge-Kutta-Gill. Pseudolinguaggio q = sqrt(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 ya=0, i valori iniziali y; e la var. logica “fine” , per terminare il programma) • Ripeti per ipti = 1 to npti; segnapunto; • Ripeti per idiv= 1 to ndiv; Ripeti per iru = 1 to 4; deriv ( nvar è l’ordine del sistema ODE cioè il numero di equazioni differenziali del primo ordine) • Ripeti per ivar = 1 to nvar; 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 ; Fineipti; Valin; se non “fine” vai a iniziotraccia ; FINE

Procedure, Variabili, Costanti vettoriali : sono aa, bb, cc; “vettori” sono y, variabili del problema, r, derivate, ya, variabili ausiliarie, corrispondenti alle y , y 0 valori iniziali, chiesti dalla procedura valin ed assegnati 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 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.

RKGill in Matlab. Preliminari disp('Laboratorio di Meccanica Razionale 2008'); disp('allievo Mario Bianchi, matr 543007'); disp('Integratore RKGill (metodo di IV ordine ) '); q = sqrt(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]; % costanti del metodo npunti=60; ndiv=1; tf=10; step=(tf/npunti)/ndiv; % init = inizializzazione dei parametri nvar=2; ya=[0, 0]; y 0=[1, 0]; y=y 0 ; t=0; % valin = valori iniziali della soluzione cercata

RKGill: Passo di calcolo • for i=1: npunti ; xx(i)=y(1); yy(i)=y(2); tt(i)=t; for idiv=1: ndiv ; for iru = 1: 4; r=oscar(y); % calcolo derivate for ivar = 1: nvar; % 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; • end ; end; t=t+step; end; plot(xx, yy) ; hold

passo RKGill Matlab. Variante. • for i=1: npunti ; xx(i)=y(1); yy(i)=y(2); tt(i)=t; for idiv=1: ndiv ; • for iru=1: 4; a=aa(iru); b=bb(iru); c=cc(iru) ; r=oscar(y); %( il ciclo ivar=1: nvar è sviluppato internamente da Matlab) rr=r*step; aum=a*(rr-b*ya); y=y+aum; ya =ya + 3*aum - c*rr; • end; t=t+step; end; plot(xx, yy) ; hold

ODE condizioni al contorno • Come semplicissimo esempio di problema al contorno da risolvere con metodi ODE, esemplare per risolvere problemi reali assai più difficili, propongo all’allievo di studiare il moto di un proiettile che deve raggiungere un bersaglio assegnato ( esempio elementare di più complesse tecniche di “shouting”). • Si assegnino la posizione del bersaglio e condizioni iniziali di prima ipotesi. • Si valutino le correzioni successive da dare alle condizioni iniziali perché il proiettile raggiunga il bersaglio.

Precisione delle soluzioni numeriche, sviluppo di Taylor • Per comprendere il significato della espressione “Precisione del quarto ordine”, proveremo a valutare le funzioni più semplici ed accessibili mediante i rispettivi sviluppi in serie già studiati nei precedenti corsi di matematica, utilizzando alcuni accorgimenti per ridurre errori numerici nelle operazioni di somma con numeri molto piccoli.

Crescita e decadimento esponenziali 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 esponenziale. • 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-Mc. Laurin 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) per l'equazione dy/dt = A y, con A per ora costante, • 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: d 2 y/dt 2 = A dy/dt = A 2 y; d 3 y/dt 3 = A d 2 y/dt 2 = A 2 dy/dt = A 3 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, in base al valore stimato del resto della serie che verrà trascurato nel calcolo effettivo, il calcolo va organizzato evitando di valutare grandi fattoriali con i loro inversi e potenze successive sempre più piccole di AT, ma raccogliendo i fattori comuni per ridurre l’errore di troncamento che si verifica se si sommano termini via più piccoli y=yo(1+AT(1+1/2 AT(1+1/3 AT(1+1/4 AT(1+1/5 AT(…)))))).

Errore relativo, resto della serie, errore totale accumulato • Come esercitazione (obbligatoria) si costruisca una procedura che chieda in input l’ordine massimo della ridotta da calcolare e realizzi il prodotto dei termini annidati a partire dal più interno. • In uscita la procedura dovrà fornire, oltre al risultato approssimato voluto, anche una stima del primo termine trascurato e del resto. • Si confronti tale errore con il valore della funzione nel caso di esponenziali crescenti e nel caso di esponenziali decrescenti.

Metodi passo Per valutare y a tempi t > R (raggio di convergenza), è possibile scegliere tau < R e quindi usare il valore y(tau) 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 di integrazione che adottiamo. Attenzione: gli errori commessi in ciascun passo vanno sommati in valore assoluto per ottenere la stima

Passo di calcolo e numero di passi • 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 e 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 d 2 x/dt 2= - dy/dt= - x d 2 y/dt 2= dx/dt= - y d 3 x/dt 3= - dx/dt= y d 3 y/dt 3= - dy/dt= -x eccetera. • Svolgendo quindi lo sviluppo in serie di Taylor: x(T)= x(0) – y(0) T – ½ x(0) T 2 + 1/6 y(0) T 3 …… y(T)= y(0) + x(0) T – ½ y(0) T 2 – 1/6 x(0) T 3 ……

Esempi : seni & coseni, segue E quindi si possono riorganizzare i termini: x(T)=x(0) (1 -1/2 T 2+1/(2 x 3 x 4)T 4– 1/(2 x 3 x 4 x 5 x 6) T 6 …. )– y(0)T(1 -1/(2 x 3)T 2 +1/(2 x 3 x 4 x 5)T 4– 1/5040 T 6…. ) y(T)=y(0) (1 -1/2 T 2 + 1/24 T 4 – 1/720 T 6 …. ) + + x(0)T(1 -1/6 T 2 +1/120 T 4 – 1/5040 T 6 …. ) Anche ora è opportuno trasformare le somme di termini rapidamente decrescenti in prodotti come p. es. x(T)=x(0) (1 -(1/2 T 2(1 -1/(3 x 4)T 2(1+1/(5 x 6) T 2 (…. )– y(0)T(1 -1/(2 x 3)T 2(1+(1/(4 x 5)T 2(1 -1/(6 x 7)T 2 (…. ) 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(x 2)/dt +d(y 2)/dt= 0 da cui si deduce che la somma x 2+y 2 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 costruiscono un linguaggio per la matematica numerica ed in particolare per le ODE. • Dovrebbero essere tutti completamente noti e servono per costruire strumenti di lavoro assai sofisticati ed efficaci. • Sono altresì utilizzabili per verificare e per scegliere metodi risolutivi adeguati ai problemi che si pongono all’attenzione dello studioso.

Esercizi • 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 4 R , 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 d 2 x/dt 2 = - 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= -ω2 x. 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 : d 2 x/dt 2 = -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 condizioni iniziali i valori x e v, al tempo iniziale, t 0 e li chiamo x 0 e v 0. Posso esprimere la soluzione nella forma dello sviluppo in serie di Taylor : x = x 0 + dx/dt (t- t 0) + ½ d 2 x/dt 2 (t- t 0)2 + + (½ *1/3) d 3 x/dt 3 (t- t 0)3 + . . . con le successive derivate espresse dalla equazione originale: d 2 x/dt 2 = dv/dt = -x ; d 3 x/dt 3=- dx/dt =-v , calcolate nel tempo t 0.

Oscillatore, seni, coseni. Si ottiene : x(t) = x 0 + v 0 (t- t 0) - ½ x 0 (t- t 0)2 + -(1/6) v 0(t- t 0)3 + 1/24 . x 0 (t- t 0)4 +. . . Cioè x(t) = x 0 ( 1 – (t – t 0)2 (1/2 – 1/12 (t – t 0)2 (…. ))) + v 0 (t – t 0) ( 1 – 1/6 (t – t 0)2 (…))) Si riconoscono nelle espressioni che moltiplicano x 0 e v 0 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.

Oscillatore e serie di Taylor 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. Si realizzi il programma Matlab che produca i grafici richiesti e si conf rontino con quelli dei programmi per ODE.

Equazioni Differenziali Ordinarie • Metodi numerici per ”Problemi di Caucy” (vedi il corso di Analisi 2 e di Calcolo Numerico). • 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

Altri problemi differenziali rilevanti. Popolazioni problemi lineari e non. • Popolazioni : Crescita lineare : dy/dt= a y • Crescita con risorse limitate : cercare gli equilibri e il punto di crescita massima : dy/dt= a y –by 2 • Competizione tra popolazioni : cercare gli equilibri dx/dt= a x – b x y – cx 2 ; dy/dt= d y – e x y – fy 2 • Prede e predatori (Equazioni di Volterra e di Lotka): dx/dt= a x–b xy ; dy/dt= –d y+ e xy. Anche per queste cercare i punti di equilibrio e valutare l’andamento dei cicli. La presenza di cicli chiusi fa sospettare l’esistenza di una funzione potenziale.

Moto Kepleriano • Il moto di un pianeta attorno al sole è descritto dalle equazioni : • d 2 P/dt 2=-OP/R 3 con x(0)=x 0, y(0)=y 0 e R 2=x*x+y*y; R 3= R 2*sqrt(R 2); • quindi: d 2 x/dt 2=-x/R 3 ; d 2 y/dt 2=-y/ R 3 • Definite le velocità: v=dx/dt e w= dy/dt si ottiene un sistema di 4 eq. diff. in 4 incognite: • dx/dt=v; dy/dt=w; • dv/dt=-x/R 3; dw/dt=-y/R 3 da integrare con opportuni valori iniziali.

Keplero In coordinate rettangolari x ed y: • r 2 = x 2 + y 2 ; r 3 = r 2 sqrt(r 2) • dx/dt=v ; dv/dt=-Kx/(r 3); • dy/dt=w; dw/dt=-Ky/(r 3) 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 (d 2 r/dt 2 - r ( d θ /dt)2) u + (r(d 2θ /dt 2) +2 (dr/dt) (dθ/dt)) v = (d 2 r/dt 2 - r ( d θ /dt)2) u + 1/r (d(r 2(dθ/dt)) /dt) v.

Altri moti centrali Moti centrali sono anche quelli prodotti su un punto da forze elastiche e in generale da forze che passano per un punto fisso, p. es. proiezione sul piano orizzontale del moto di un punto soggetto a gravità vincolato ad un cono rovesciato o vincolato su una forma per budino (liscia, ovviamente): quest’ultimo è un modello proposto per le forze atomiche sugli elettroni. Il metodo di calcolo già realizzato per sin e cos può essere esteso immediatamente al problema con quattro equazioni del primo ordine, tracciando poi diverse ellissi in corrispondenza di diversi valori iniziali.

Riduzione del numero di variabili Il moto centrale si presta ad illustrare un metodo fondamentale della Meccanica razionale: l’uso degli integrali primi o costanti del moto per ridurre il numero delle variabili. Già nel caso di sin e cos, si può descrivere l’evoluzione del punto rappresentativo nel piano (x, y) come moto a velocità angolare uniforme su una circonferenza di raggio costante: R=R 0, d /dt=u, = 0+ u t . Il problema fisico descritto dalle variabili cartesiane x ed y, che si evolvono nel tempo in modo “complicato”, può essere descritto dalle nuove variabili R costante e lineare in t.

Momento angolare Nel caso dei moti centrali la forza è diretta verso il centro, pertanto si annulla il prodotto OP F, si annulla la grandezza OP d 2 OP/dt 2 che è derivata di L= OP d. OP/dt =c 3 r 2 d /dt (verificare!!! : è domanda di esame). La costanza di L impone che OP e v. P siano sempre sullo stesso piano, ortogonale ad L. Possono dunque essere usate coordinate polari con u e w, versori ortogonali e w = du /d , u = - dw /d . OP=ru; d. OP/dt = dr/dt u+r du/dt= u dr/dt + w r d /dt; dr/dt è la velocità radiale, r du/dt = w r d /dt è la velocità trasversa. L’energia cinetica è T= ½ m ((dr/dt )2 + r 2 ( d /dt)2 ). In astronomia si osservano direttamente e d /dt, mentre r va calcolata con operazioni elaborate.

Moti centrali: Accelerazioni d 2 OP/dt 2 = d 2 r/dt 2 u-r (d /dt)2 u +2 dr/dt d /dt w = =(d 2 r/dt 2 -r (d /dt )2 ) u + 2 w(dr/dt d /dt +r d 2 /dt 2 ). Il termine in u è l’accelerazione radiale, il termine in w è l’accelerazione trasversa, che è nulla per le leggi di Keplero: moltiplicando la accelerazione trasversa per r si vede che è costante ||N||=r 2 d /dt =k, con N= OP d. OP/dt. Se ne deduce d /dt = k/ r 2, da usare nell’espressione di T. Resta costante anche l’energia totale 1/2 (d. OP/dt )2 + V(r ). Si ottiene un sistema di equazioni differenziali del secondo ordine nelle variabili radiali r e dr/dt, con i dati iniziali r 0 , 0, k , da rappresentare nel piano delle fasi (r, dr/dt). e d /dt vengono dedotte da r e dr/dt e mediante si può tradurre il grafico nelle consuete coordinate cartesiane x ed y. Il moto è descritto dalla funzione di Lagrange L* = ½(dr/dt)2–½ k 2/r 2 –V( r), con V( r) =1/r , potenziale. Anziché 1/r il potenziale potrebbe essere costante (punto soggetto a gravità su di un cono rovesciato) oppure = – r 2, forza elastica.

Equazioni di Lagrange Posto v = dr/dt, velocità radiale, la Funzione di Lagrange è L=½v 2–½ k 2/r 2 –V( r) ; il piano delle fasi è ( r, v ). Il potenziale kepleriano è V= 1/r. La funzione L diviene L*=½v 2–V*( r) con V*= ½ k 2/r 2 - 1/r, che esprime in r ed r 2 la medesima dinamica. Si studino anche i casi dei potenziali (senza *) V= - r 2, e V=cost. , potenziale elastico e potenziale di gravità su un cono rovesciato. L’Equazione di Lagrange è : d/dt( L/ v)- L/ r =0; da cui: dv/dt - k 2/r 3 +1/r 2 = 0 con dr/dt=vr. Integrando la coppia di ODE si ottengono le curve da rappresentare nel piano delle fasi (r, vr). Si usi inoltre l’informazione d /dt = k/ r 2 : il sistema delle tre equazioni si integra con i metodi studiati per le O. D. E. , ottenendo r(t) e (t) con cui si possono tracciare i grafici cartesiani delle traiettorie.

Moto Kepleriano 2 Conforme la legge di Newton l’accelerazione è radiale: a= (-K/r 2)u. L’annullarsi di d(r 2(dθ/dt)) /dt conduce all’integrale primo r 2ω = cost; con ω = dθ/dt , e dθ/dt = cost/r 2 ; dr/dt = dr/dθ dθ/dt = - cost d(1/r)/dθ. d 2 r/dt 2 = dθ/dt d(dr/dt) / dθ = - cost dθ/dt d 2 (1/r)/dθ 2 = - cost 2 (1/r 2) d 2 (1/r) / dθ 2 = - K/r 2 . Appare come se la variabile naturale fosse anziché r la y=1/r. si può infatti scrivere : -cost 2(1/r 2) +d 2 (1/r)/dθ 2 = = cost 2 y 2 +d 2 y/dθ 2 = cost 2 /K, che è analoga alla equazione del pendolo, con un termine non omogeneo.

Moti centrali coord, radiali

Il piano delle fasi • 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 x 2+v 2. Le condizioni iniziali con x 0 = 0 sono rappresentate da punti sulla retta verticale, quelle con v 0 = 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 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. Ciò non vale per altri problemi. In Lotka, il verso è antiorario.

Equazioni differenziali ordinarie non lineari • • • Esempi proposti: Pendolo, crescita con risorse limitate, Lotka Volterra Moto kepleriano, bipendolo. Equazioni con ciclo limite: Van der Pol. 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 (2 nπ , 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 v 2+2 mgh= cost. e periodo T dalla T=4 0 xmax (1/v)dx.

Periodo del pendolo • L’equazione è: dv/dt=-sin(x) ; moltiplico ambo i membri per dx=vdt : dv/dt vdt= -sin(x) dx. Integrando i due membri tra estremi corrispondenti si ha ½ v 2 –½v 02=cos (x)-cos(x 0). • Posto ad es. x 0 = /2 e v 0=0, ½ v 2 =cos (x), v= dx/dt= (2 cos (x)), e il tempo di caduta da x 0= /2 a x=0 è T 1/4 = dt = (dt/dx) dx= (1/v)dx, integrale definito da valutare con i metodi dei trapezi o delle parabole o dei punti di Gauss. • Il periodo della oscillazione completa è 4 T 1/4.

Periodo del pendolo

Pendolo forzato. Sezione di Poincaré Sul consueto pendolo agisca anziché la forza di gravità costante una forza periodica del tipo: g*=g+b cos( t) con = (g/l): L’equazione è allora dv/dt=-sen( x) + cos ( t); con, p. es. , : =2 / 20. Il consueto piano delle fasi ( , v) non è adeguato a rappresentare il moto, che ha un carattere tridimensionale, dovendosi considerare la variabile t, da cui dipende la forzante. Per tali situazioni il matematico Poincaré ha suggerito di rappresentare sul piano ( , v) le “sezioni” della traiettoria tridimensionale tagliate dai successivi piani ai tempi t = 2 n /. I punti così ottenuti si dispongono in figure di cui è possibile discutere svariate proprietà. Per realizzare le figure con i programmi già realizzati diviene ora particolarmente opportuno l’uso del parametro ndiv.

Forzante dipendente dal tempo • Per valutare correttamente la forzante che dipende dal tempo in modo esplicito è utile introdurre anche t tra le variabili da integrare, scrivendo un’equazione differenziale con derivata r = 1 costante e val. iniziale y 0 = 0. • Il sistema che proponiamo per descrivere il pendolo forzato ha quindi tre equazioni in tre incognite: tempo t, velocità v, posizione x, con derivata del tempo=1 costante, derivata della posizione = v, velocità, derivata della velocità = a, accelerazione, con a= -sin( x) + cos( t). Suggerisco di sperimentare = 1, = 2 /6 ed con valori diversi da 0. 1 a 0. 5. • Si utilizzi il parametro ndiv, numero di passi tra un punto e l’altro, che sia multiplo di 6, così da segnare un punto ogni giro completo della forzante e non rappresentare gli esiti dei passi intermedi.

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 = a 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.

Le m-function per il calcolo delle derivate • Per ogni problema va assegnata una particolare funzione vettoriale che assegna le derivate delle funzioni incognite valutate al tempo attuale e che può essere memorizzata come m-file • 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 function r = lotka(y); r(1)= y(1)*α(1 - y(2)*β/α); r(2)= -y(2)*γ(1 - y(1)*δ/γ);

Oscillatore armonico smorzato e amplificato L’equazione per l’oscillatore armonico con effetti di smorzamento o di amplificazione è: d 2 x/dt 2 + β dx/dt +x = 0; da scrivere come sistema di due equazioni: dv/dt=- β dx/dt -x; dx/dt=v , 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 è: d 2 x/dt 2 + µ(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( x 2 /a 2– 1) che è negativa se x è compreso tra -a e +a e positiva fuori di tale intervallo. L’equazione di Van der Pol è : d 2 x/dt 2+β ( x 2– 1) dx/dt + x = 0, ottenuta con un fattore di scala nelle x, evidenziando β, parametro adimensionale significativo. • Si può scrivere: d/dt ( dx/dt + β ( x 3/3 –x ) ) =-x, cioè dy/dt = -x/ β, posto β y= dx/dt + β ( x 3/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)= x 3/3 –x. Come d’uso faremo attenzione al caso β >> 1.

Van der Pol osservazioni • dx/dt è = 0 sulla cubica γ: y= x 3/3 –x, ivi 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 β. 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: figure

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 - x 2. 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 – ζ x 2 e che produce un altro punto critico in (y=0, x= α / ζ ). I termini con τ rappresentano il prelievo della pesca. Si osservi che in presenza di prelievo il punto di stazionarietà diviene y= (α-τ)/ β , x= (γ+τ)/δ dx/dt = α x – β x y – ζ x 2 - τx • dy/dt = - γ y + δ x y – τy

Problema di Lotka-Volterra

Lotka Volterra modificato dal termine di affollamento

Altri problemi • Pendolo smorzato

O. D. E Lineari. Part. Libera • Particella libera: d 2 x/dt 2 = 0 ; Per la soluzione introduciamo la variabile v, velocità, definita da v = dx/dt; • l’equazione del secondo ordine (d 2 x/dt 2 = 0) diviene una coppia di equazioni del primo ordine: dv/dt = 0 ; dx/dt = v ; con le condizioni iniziali x(0) = x 0 , v(0) = v 0. • Vedi più avanti le traiettorie rettilinee parallele nel piano delle fasi.

Moto per gravità • d 2 x/dt 2 = 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) = x 0 , v(0) = v 0. • Sarà richiesto di tracciare anche per questo caso le traiettorie nel piano delle fasi.

Oscillatore e Repulsore Armonico L’equazione d 2 x/dt 2 = - ω2 x , 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.

Errore O(h 2) ed O(h 3) Una prima approssimazione di y è: y* = yi + h G(yi, ti) = yi+1 + O(h 2) 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))= ½ h. Gy(y. M, ti)(y*-yi+1)) dove la derivata parziale di G va valutata in un opportuno y. M , intermedio tra y* e yi+1. L’errore quindi è O(h 3).

Altre formule di 2 do ordine • Un’altra formula al 2 do 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: • k 1=h G(y, t) ; k 2=h G(y+2/3 k 1, t+2/3 h) ; Δy= ¼ ( k 1+3 k 2)

Formule del 4 to ordine • È molto diffusa la formula del 4 to ordine: k 1=h. G(y, t); k 2=h. G(y+ ½ k 1, t+ ½ h); k 3=h(G(y+ ½ k 2, t+ ½ h); k 4=h(G(y+k 3, t+h); yi+1= yi + 1/6 (k 1+2 k 2+2 k 3+k 4)

Formule parametriche (complementi) • Tale formula del 4 to ordine è dedotta per s=1 da una formula nel parametro libero s: k 1=h. G(y, t); k 2=h. G(y+ ½ k 1, t+ ½ h); k 3=h(G(y+ ½k 1+1/(2 s)(k 2 -k 1), t+ ½ h); k 4=h(G(y+k 2 + s (k 3 -k 2), t+h); yi+1= yi + 1/6 (k 1+4 k 2+2 s(k 3 -k 2)+k 4)

Altre formule (complementi) La formula Runge – Kutta – Gill che suggeriamo impiega il valore s= (0. 5). Un’altra formula al 4 to ordine è: k 1=h. G(y, t); k 2=h. G(y+ 1/3 k 1, t+ 1/3 h); k 3=h(G(y+k 2 -1/3 k 1, t+2/3 h); k 4=h(G(y+k 3+k 1 -k 2, t+h); yi+1= yi + 1/8 (k 1+3 k 2+3 k 3+k 4).

Precisione del 4 to 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. Nella variante che lascia il ciclo ivar al Matlab, anche rr ed aum assumono la dimensione nvar.

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 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 , y 0 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.

Forze elastiche su punti allineati

Due punti con forze elastiche. FQ = -k 1 x + k 2 (y-x) ; FP = - k 2 (y-x) ; da cui: m. Q d 2 x/dt 2=- (k 1 + k 2 ) x + k 2 y; m. Pd 2 y/dt 2= k 2 x - k 2 y. Poniamo, per questo esempio, le costanti m e k uguali ad 1. Il sistema di 4 equazioni del primo ordine diviene : dv/dt = -2 x + y; dw/dt = x – y; dx/dt = v; dy/dt=w nelle 4 incognite v, w, x, y, e per esso sono utilizzabili i metodi ODE proposti, in particolare RKGill. Vanno rappresentate le traiettorie sul piano x, y, per diversi valori delle coordinate iniziali x 0 e y 0. Si potranno porre v 0 e w 0 entrambe = 0.

Autovalori Si studia se, come accade nel moto armonico in una sola variabile, può accadere che accelerazione e spostamento siano tra loro proporzionali : dv/dt= x, dw/dt= y. Esplicitando si ha : x=-2 x+y; y=x-y, cioè (-2 - )x +y = 0; x+(-1 - )y = 0 È necessario che il discriminante del sistema di equazioni sia 0 : (– 2 –λ)(– 1–λ) – 1 = 0 ; λ 2 + 3λ + 1 = 0, soluzioni ½ (-3 ± 5) : λ 1= - 0. 381966 ; λ 2 = - 2. 618, gli autovalori. Introdotti gli autovalori nell’equazione matriciale si ottiene x 1= 0. 618 y 1 e x 2 = - 1. 618 y 2. Le coppie (x 1 , y 1) , (x 2 , y 2) sono gli “autovettori” che danno luogo ai “modi normali” di oscillazione: in ciascuno dei modi normali x e y cambia-no simultaneamente mantenendo il loro rapporto costante.

Valori iniziali Nel “modo” 1, cioè lungo la retta che passa per O e pe (x 1 , y 1), x ed y hanno lo stesso segno (retta nei quadranti 1 e 3), nel “modo” 2 , (x 2 , y 2) hanno segno opposto (retta nei quadranti 2 e 4). Le due rette vanno pensate come possibili nuovi assi di riferimento. Qualunque altro vettore iniziale ossia qualunque quaterna x 0, y 0, v 0, w 0 può essere rappresentata come (x 1 , y 1) + (x 2 , y 2) per le posizioni e (v 1 , w 1)+ (v , w 2) per le velocità e i coefficienti , , , restano costanti, mentre modo 1 e modo 2 evolvono ciascuno con la propria frequenza, dando luogo ad evoluzioni anche complicate.

Esercizio: più punti in un box. • Sviluppare un analogo esempio, in cui i due punti A e B, legati da tre molle sono racchiusi tra due punti fissi, O e Q. • Sviluppare un esempio con tre punti A, B, C legati da quattro molle tra i punti fissi O, Q • Nella fisica dei solidi i modi di oscillazione di un continuo modellizzato come limite per infiniti punti di questi esempi si chiamano “fononi”.

Modello di continuo elastico • L’esempio, che può essere ampliato a più punti connessi da forze elastiche, ed è di fatto un modello lineare per continui elastici e per reti anche multiconnesse, si presta ad illustrare diversi fenomeni e introduce vari metodi di calcolo, tra cui le “differenze finite” e “gli elementi finiti”. In particolare qui evidenziamo i “modi normali delle oscillazioni”, ossia oscillazioni che coinvolgono insieme tutti i punti del sistema, ma con frequenze diverse per ogni diverso “modo”. In ciascun modo normale tutti i punti si muovono simultaneamente di moto armonico, con la medesima frequenza e mantenendo sempre proporzionali tra loro i valori delle coordinate. Ogni altro movimento è rappresentabile combinazione lineare di modi normali, selezionati dalle condizioni iniziali.

Bipendolo Si consideri il dispositivo in figura: O è un punto fisso, OQ e QP hanno lunghezza costante l=1 e direzioni misurate da (t) e da O (t), i punti Q e P sono dotati di massa m=1 e sono soggetti a Q gravità. Si valuti Energia Cinetica e Potenziale del sistema. P P

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= ½ m. Q v 2 + ½ m. P(v+w)2, V= g(m. Q h. Q + m. P h. P), W= ½(m. Q +m. P) v 2 + ½m. P (w 2 + 2 v w) = ½ (m. Q +m. P) θ’ 2 + ½ m. P (φ’ 2 +2 θ’ φ’ cos (θ-φ)). Posti ad uno i moduli OQ e QP ed indicando con θ’, φ’ le derivate temporali delle coordinate angolari. • h Q= - cos θ; h. P = h. Q - cos φ; V= - g( (m. Q + m. P) cos θ + m. P cos φ)

Grafica bipendolo

Equazioni di Lagrange • OQ=OQ(θ); QP=QP(φ); v. Q=d(OQ)/dt= (OQ)/ θ • dθ/dt = (OQ)/ θ θ’ • v. P =d(OP)/dt= (OP)/ θ θ’ + (OP)/ φ φ’ • Nota che (OP)/ φ è uguale a v. P/ φ’ • Si cerca una espressione per a. P • (OP)/ θ. Si valuti v. P • (OP)/ θ che si può calcolare come v. P • v. P/ θ’ e quindi ½ (v 2 P) / θ’. [θ’=dθ/dt] • Lagrange suggerisce di valutare a. P • (OP)/ θ come d/dt( (v 2 P/2)/ θ’)- (v 2 P/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” e sono: d/dt ( W/( θ’) - W/ θ = - V/ θ d/dt ( W/( φ’) - W/ φ = - V/ φ nel presente esercizio: • (m. P + mq) d 2θ/dt 2 + mq cos (θ-φ) d 2φ /dt 2 = = -g (m. P + mq) sin θ+ mq (d φ / dt) 2 sin (θ – φ) • mq cos (θ-φ) d 2θ/dt 2 + mq d 2φ /dt 2 = = -g mq sin φ - mq (d θ /dt) 2 sin (θ – φ)

Forma matriciale (m. P +mq) , mqcos(θ-φ) d 2 θ/dt 2 | | = mqcos(θ-φ) , mq d 2 φ /dt 2 - g(m. P +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. ugual 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 tipo 1 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. Tipo 2

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

movie • Sperimentare la istruzione Matlab per rappresentare sullo schermo il movimento studiato, disegnando le successive configurazioni del sistema (p. es. il bipendolo o le tre palle nel box) in successive frames, da presentare il sequenza.

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

P. D. E. • La sigla PDE significa: Equazioni Differenziali alle Derivate Parziali. Molti importanti problemi studiati dalla Fisica Matematica e dalla ingegneria sono regolati da vari tipi di P. D. E. . • Ci limitiamo a proporre qualche esempio per illustrare metodi numerici di soluzione per problemi di equilibrio elastico, equilibrio termodinamico, evoluzione della temperatura, propagazione ondosa in fluidi ed in corpi elastici. • Sovente per la soluzione di effettivi problemi ingegneristici vanno costruite matrici di grandi dimensioni, con pochi coefficienti su ogni riga (“matrici sparse”), da utilizzare ripetutamente (“metodi iterativi”).

P. D. E. ellittiche • Equilibrio termico in un continuo • Distribuzione di sforzi in una struttura reticolare od in un continuo • Distribuzione delle pressioni nel suolo • Correnti e tensioni in reti elettriche • Pressioni e portate in reti idrauliche in regime stazionario

P. D. E. paraboliche • Evoluzione delle temperature in un continuo • Evoluzione di concentrazioni o di pressioni in un mezzo dispersivo

P. D. E. iperboliche • Propagazione di onde in fluidi od in corpi elastici : onde sonore , determinazione di fratture nei solidi; onde sismiche, determinazione delle strutture del sottosuolo

Distribuzione di temperature

Condizioni al contorno • problema 1) Lati orizzontali: si supponga che sul lato superiore T=100, sul lato inferiore T=0; lati verticali: T=100 sul lato sinistro, T=0 sul lato destro. • problema 2) : Si suppongano condizioni di simmetria sui lati orizzontali (T sui punti del lato inferiore uguali alle corrispondenti del lato superione), e sui verticali: T=100 lato sinistro, T=0 destro. • All'equilibrio sui nodi interni la temperatura in ciascun nodo è la media dei nodi adiacenti; • Per l'evoluzione nel tempo, in ciascun nodo l'incremento della temperatura nel tempo dt è valutato per differenze finite come somma delle derivate seconde spaziali.

Evoluzione della temperatura • L'evoluzione nel tempo del campo di temperature è regolata da un'equazione PDE di tipo parabolico: T/dt= 2 T/ x 2+ 2 T/ y 2. • Un possibile metodo numerico per la ricerca della evoluzione temporale, dopo aver definito un opportuno reticolo di nodi, di dimensione n 2, propone un sistema di n 2 equazioni algebriche: per ciascun nodo la derivata temporale della temperatura è uguale alla somma delle due derivate seconde nello spazio, derivate da approssimare con le “differenze finite” delle temperature nel nodo e nei quattro nodi adiacenti.

Differenze finite per le derivate seconde • La derivata prima nello spazio è approssimata da : (Ti+1 –Ti)/(xi+1 -xi); la derivata seconda nello spazio sarà approssimata da: 2((Ti+1 –Ti)/(xi+1 -xi) - (Ti –Ti-1)/(xi-xi-1)) /(xi+2 -xi-1) ; Cioè, uniformando le delta x : 4( - 2 Ti +Ti-1 +Ti+1)/ (xi+2 -xi-1)2 poiché al denominatore delta x è raddoppiato,

Soluzione stazionaria • La soluzione stazionaria è regolata da una PDE di tipo ellittico: è nulla la somma delle derivate parziali seconde in x e in y della temperatura: 2 T/ x 2+ 2 T/ y 2 = 0. • Una soluzione numerica alle “differenze finite” di tale equazione nel reticolo di nodi definito più sopra propone di calcolare la temperatura in ciascun nodo come media delle temperature nei nodi adiacenti.

Matrice sparsa • Risulta un sistema di n 2 equazioni algebriche lineari, descritto da una matrice di n 2 righe ed n 2 colonne. In ciascuna riga solo 5 termini sono diversi da zero: i termini sulla diagonale hanno valore 4 , mentre i due termini adiacenti alla diagonale, assieme ai due termini distanti n dalla diagonale hanno valore -1. • Risolvere tale sistema di eq. lineari simultanee col metodo generale di Gauss, richiede n 6/3 moltiplicazioni, la maggior parte delle quali con fattori =0. Ma se la matrice è “ben condizionata” le equazioni possono essere risolte “in successione”, ripetendo iterativamente la ricerca della soluzione. • Qui proponiamo un metodo iterativo “economico”.

Soluzione iterativa per l'equilibrio a=zeros(10) ; dim=10; niter=50; for j=1: dim ; a(1, j)=100; a(j, 1)=100; end; %condizioni ai margini for k=1: niter; for j=2: dim-1; for i=2: dim-1; a(i, j)=(a(i-1, j)+a(i+1, j) +a(i, j-1)+a(i, j+1))/4; end; a

Soluzioni esplicita ed implicita Per lo studio della evoluzione nel tempo le nuove temperature al tempo t+ t nel punto (x, y) possono essere ottenute esplicitamente dalle temperature nel punto corrispondente e nei punti adiacenti al tempo t (soluzione esplicita diretta). Attenzione: l'intervallo temporale t deve essere minore di 1/2 ( x 2+ y 2) per evitare sgradevoli instabilità numeriche. Si può anche impostare un sistema di n 2 ODE. In alternativa si scrive un sistema di equazioni lineari nelle variabili al tempo t+ t e con termini noti la medesime variabili al tempo t (soluzione implicita).

Soluzione esplicita e soluzione iterativa Per la ricerca della distribuzione di equilibrio (PDE parabolica) , studiare dapprima casi con pochi punti (2 x 2 o 4 x 4), risolverli iterativamente a mano, quindi costruire il risolutore iterativo Matlab e verificare il numero di iterazioni necessario per ridurre il residuo a 1/1000. • Quindi costruire la matrice per il calcolo diretto delle temperature e verificare che i risultati diretto ed iterativo coincidano. • Si osservi che la matrice ha solo 5 elementi diversi da zero in ogni riga. Un notevole risparmio nella occupazione di memoria si ottiene utilizzando la modalità di calcolo delle “matrici sparse”, di cui abbiamo suggerito un primo semplice esempio.

Struttura dati per Matrici Sparse L’informazione contenuta in una matrice sparsa può essere rappresentata dal “vettore delle adiacenze” che contiene, per ogni riga, gli indici ed i valori di ogni elemento diverso da zero. Per il problema del calore sopra descritto, costruire il vettore delle adiacenze e studiare come impiegarlo per realizzare la soluzione iterativa del problema. Il metodo diviene strategico per la soluzione di problemi di “calcolo delle strutture”con il metodo degli “elementi finiti”, in cui il numero di nodi può divenire molto grande.

Vettore delle adiacenze • Per ogni matrice “sparsa” si può costruire una coppia di vettori che contiene, per ogni riga della matrice l’indice ed il valore degli elementi diversi da zero, assieme ad un altro vettore che contiene per ciascuna riga l’indice del primo elemento registrato. • Strutture analoghe sono impiegabili per descrivere i “grafi” e per descrivere problemi da risolvere col metodo degli “elementi finiti”.

Un problema iperbolico: ottimizzare un carburatore Un carburatore è schematizzato come un tubo uniforme suddiviso in 9 sezioni che identificano 8 intervalli. In ciascuna sezione definisco la velocità dell’aria, in ogni intervallo definisco la pressione. . La velocità aumenta se è positiva la differenza di pressione tra monte e valle, la pressione aumenta se è positiva la differenza delle velocità a monte ed a valle. Le condizioni al contorno possono essere ad es. pressione costante nella prima cella a monte e pressione variabile assegnata all’ultima sezione a valle. Si può costruire un sistema di ODE di ordine 16, per il quale è adeguato il metodo RKGill già studiato.

Modello di tubo vibrante Nel volume che precede la prima sezione suppongo p costante (aria atmosferica), nel volume che segue l’ultima suppongo p variabile col coseno (pistone che si muove); nelle sezioni sono definite le velocità dell’aria, nei riquadri le pressioni. Valgono le relazioni v/ t= - p/ x, p/ t= - v/ x. Se derivo la prima relazione per x e la seconda per t ottengo la formulazione delle eq. iperboliche come differenze di derivate seconde nel tempo e nello spazio per v e similmente per p: v/ t - v/ x =0; p/ t - p/ x =0 . P=cost P=cos(t)

Un vettore di ODE per il problema iperbolico • Sono 8 celle in cui sono incognite le pressioni e 9 sezioni in cui sono incognite le velocità, in tutto il sistema è di 17 equazioni in 17 incognite, con la pressione dopo la sezione finale variabile assegnata, dipendente dal tempo e quella prima della sezione iniziale costante (aria atmosferica). Anche il tempo può essere considerato variabile da integrare, con valore iniziale 0 e derivata 1, per semplificare la gestione di integratori a passo variabile, così le equazioni sono 18 o 19 se si vuole tenere nel medesimo vettore anche il valore della pressione costante.

Struttura dati • Una possibile struttura dati per il problema del tubo vibrante è : y(1)=tempo, con derivata r(1)=1 e valore iniziale y 0(1)=0, y(2)=p 0, con valore=1, costante e derivata r(2)=0, y(3)= velocità alla sezione 1, con r(3)= p 0 -p 1=y(2)-y(4) e valore iniziale 0, y(4) = pressione nel box 1, con derivata r(4)= v 1 v 2=y(3)-y(5) e valore iniziale=1, fino alla pressione finale che ha valore iniziale=1 e derivata =-0. 5 sin (t) = -0. 5 sin(y(1)). • Una opportuna funzione vettoriale costruisce tutte le r in funzione delle y, e sarà invocata all’interno del consueto integratore Runge Kutta Gill, o di altri opportuni integratori offerti da Matlab.

Metodo “leap frog” per PDE iperboliche Ricordando il “metodo simplettico”, la differenza tra valori di p in celle contigue definisce esplicitamente l'incremento del valore di v nella sezione intermedia e similmente la differenza tra valori di v in sezioni successive definisce esplicitamente l'incremento del valore di p nella compresa, a condizione di intendere che le v iniziali siano state definite al tempo t= -1/2 dt e le v al tempo dt/2 siano calcolate mediante le p al tempo t=0. Le p al tempo t=dt saranno calcolate usando le v al tempo t=dt/2. Ciò consente un veloce calcolo esplicito.

• • • Indice lezioni Presentazione del corso. Matlab. Eulero per sin e cos. Moti centrali, popolazioni, Van Der Pol. Metodi simplettici. Heun. Predictor corrector. RKGill. Serie di Taylor. Errori di approssimazione e di troncamento. Passi piccoli vs. ordine alto. Forze elastiche su due punti. Punti e molle in un box. Autovettori. Bipendolo. Eq. Di Lagrange. Animazioni. Pendolo forzato. Sezione di Poincaré. (Test intermedio). Problemi ellittici e parabolici. Temperature. Differenze finite, matrici “sparse”. Metodi diretti e metodi iterativi. Problemi iperbolici. Tubi vibranti. Carburatore. Ripassi e conversazioni.