Equazioni Differenziali Ordinarie(ODE) - unina.stidue.netunina.stidue.net/Calcolo Numerico...

Calcolo Numerico II - CdL Specialistica Ing.Informatica- A.d’Alessio Equazioni Differenziali Ordinarie

1

Le equazioni differenziali sono uno dei modelli di base nella risoluzione di problemi di meccanica, meteorologia, chimica, astronomia, scienze economiche, etc. e descrivono fenomeni in movimento o in evoluzione nel tempo. Un’equazione differenziale descrive come varia una funzione, soluzione del problema, esprimendo il legame tra essa e le sue derivate. In questo capitolo non sono trattati gli aspetti teorici relativi alle equazioni differenziali, ma soltanto gli aspetti relativi alla loro risoluzione numerica. Ricordiamo alcune definizioni e concetti di base. Il problema differenziale del primo ordine di Cauchy (problema a valori iniziali) è: y’=f(t,y(t)) , t∈ [t0,T] e descrive fenomeni di evoluzione nel tempo. Se esso ammette una soluzione, ammette una famiglia di soluzioni. Una soluzione particolare y(t), t∈ [t0,T] si ottiene assegnando la condizione iniziale y(t0)=y0 Ad esempio le equazioni y' =e-t e y’=-y hanno rispettivamente le infinite soluzioni (fig1a,1b): y(t) = C-e-t e y(t)=Ce-t con C costante arbitraria.

Fig.1a Fig.1b

Esempio1 . Decadimento radioattivo. L’età dei reperti archeologici si determina sfruttando la proprietà di alcuni elementi (isotopi) di essere instabili, ossia di decadere in isotopi di altri elementi mediante l’emissione di nuclei di elio, elettroni o fotoni (elementi radioattivi). In particolare il carbonio14C decade in carbonio12C. Tale processo è descritto dall’equazione: N’(t)= -kN(t) Una famiglia di soluzioni è N(t)=Ce-k(t-to). Imponendo la condizione iniziale N(t0) = N0 si ha l’unica soluzione è N(t)=N0e-k(t-to).

Equazioni Differenziali Ordinarie(ODE)

0 0.5 1 1.5 2 2.5 3-5

-4

-3

-2

-1

0

1

2

3dy/dx=exp(-x),y(x)=c-exp(-x)

0 0.5 1 1.5 2 2.5 3-5

-4

-3

-2

-1

0

1

2

3dy/dx=exp(-x),y(x)=c-exp(-x)

0 0.5 1 1.5 2 2.5 3-4

-3

-2

-1

0

1

2

3

4dy/dx=-y,y(x)=c*exp(-x)

0 0.5 1 1.5 2 2.5 3-4

-3

-2

-1

0

1

2

3

4dy/dx=-y,y(x)=c*exp(-x)


2

Sistemi ed equazioni di ordine superiore Molti fenomeni sono descritti da sistemi di equazioni differenziali. I metodi relativi ad un’equazione del primo ordine si estendono semplicemente ai sistemi del primo ordine. Gli stessi metodi si possono utilizzare per risolvere equazioni di ordine p>1(con derivate di ordine p). Si dimostra che una ODE di ordine p si può trasformare in un sistema di p equazioni del primo ordine. Consideriamo un’equazione di ordine p : y(p) (t) = F(t,y,y’,y’’,..,y(p-1)) con condizioni iniziali : y(t0)=y0, y’(t0)=y1, ...., y(p-1)(t0)=yp-1 Posto : z1 (t)=y(t) , z2 (t)=y’(t),...., zp(t)=y(p-1)(t) l’equazione si trasforma nel sistema di m equazioni differenziali :

z’1=z2(t) z’2=z3(t) ..... z’p=F(t,z1,z2,...,zp)

con condizioni iniziali: z1(t0)=y0, z2(t0)=y1, ...., zp (t0)=yp. Ad esempio un oscillatore armonico è descritto dall’equazione del secondo ordine: y’’(t) = -y(t) Posto z(t) = (y(t),y’(t)) si ha che : z’(t) = ( y’(t),-y(t)) da cui il sistema di due equazioni in due incognite:

z’1(t) = z2(t) z’2(t)= -z1(t).

Esempio 2. Sviluppo di una popolazione. I modelli di sviluppo di una popolazione(persone, animali,batteri,…) si basano sull’ipotesi che la velocità di crescita è proporzionale in ogni istante alla popolazione esistente, se non scarseggia il nutrimento e la popolazione varia in funzione delle nascite e delle morti. Se a(t),b(t) = tasso di nascita e morte,si è osservato che il tasso di nascita diminuisce al crescere della popolazione (diminuzione di spazio,di risorse di cibo,..),per cui si pone l’ipotesi che: b(t)=b0-b1(t), b0,b1>0 e a(t)=a0<b0 e la dinamica della popolazione è descritta dal modello logistico di crescita: p’(t)=Kp(pmax-p) t>0 p(0)=p0 con K=b1 , pmax=(b0-a0)/b1.


3

Un’applicazione di tale modello è lo studio della popolazione batterica in un lago, che consente di stimare quando la popolazione supera una soglia di tolleranza e prendere provvedimenti per depurare le acque. Esempio 3. Equazioni di Lotka-Volterra(problema predatore -preda). Consiste nello studio della dinamica di una popolazione P costituita da due specie in relazione predatore-preda: la preda si procura cibo sufficiente, ma è uccisa dal predatore(volpi-conigli,antibiotici-batteri,…). Sia: r(t)= numero prede al tempo t f(t)= numero predatori al tempo t Il modello matematico che descrive tale problema è un sistema di due equazioni differenziali non lineari(non risolvibile analiticamente): r'(t)= 2r-αrf f’(t)=-f+αrf r(t0)=r0 f(t0)=f0 con α quantità positiva costante, α =0 rappresenta la crescita della popolazione in assenza di predatori. Se α >0 i predatori incontrano le prede con una probabilità proporzionale a rf (probabilità di incontro),si suppone che i predatori diminuiscano per cause naturali in assenza di prede ed aumentino in presenza di prede. Le soluzioni del sistema sono periodiche, con un periodo che dipende dal valore iniziale.

Risoluzione numerica di una ODE La maggior parte dei problemi di ODE che si incontrano nella pratica non sono risolvibili analiticamente o la soluzione non è esprimibile in forma tale da poter essere facilmente calcolata. Sono quindi necessarie tecniche numeriche per la risoluzione di una ODE. La prima risale a Newton, per lo studio del moto di una cometa, ed i calcoli erano effettuati a mano.

Condizionamento di un problema a valori iniziali Prima di risolvere numericamente una ODE bisogna analizzare il condizionamento, ossia il comportamento della soluzione in presenza di una perturbazione del valore iniziale y(t0)=y0+ε.

0 1 2 3 4 5 6 7 80

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

t

p(t)

equazione della logisticasol.per k=pmax=1 e diversi valori di p0


4

Tale perturbazione comporta la scelta di una soluzione diversa nella famiglia di soluzioni. Se la distanza tra la curva soluzione vera e quella perturbata diminuisce il problema è bencondizionato (fig 1b). Se tale distanza cresce al crescere di t il problema è malcondizionato, come visualizzato in figura (y’(t)=y(t), y(t)=Cet).

Per analizzare il condizionamento di una ODE consideriamo il problema : 1) y’(t)=λy (t) t∈[0,T] y(0)=1 la cui soluzione è y(t)=eλt. Perturbando il valore iniziale yε(0)=1+ε, si ottiene la soluzione yε(t)=(1+ε) eλt. Da cui : yε(t)-y(t) = ε eλt t∈[0,T] e quindi, per λ< 0 l’errore diminuisce, mentre per λ>0 l’errore si amplifica di un fattore < eλT. Pertanto il problema 1) è: bencondizionato se λ< 0 malcondizionato se λ> 0 Per il problema generale: 2) y’(t)=f(t,y(t)) t∈[t0,T] y(t0)=y0 si dimostra che, se f(t,y) e la derivata parziale fy ( Jacobiano dell’equazione ed indicata con J) sono continue allora il problema è: bencondizionato se fy < 0 per t∈[t0,T] malcondizionato se fy > 0 per t∈[t0,T] Ciò comporta la possibilità di avere intervalli in cui il problema è bencondizionato e intervalli in cui è malcondizionato. Esempio 4. L’equazione y'(t)= -10(t-1)y ha soluzioni y(t)=Ce-5(t-1)(t-1) . Poiché fy = -10(t-1) il problema è benconzionato per t>1 e malcondizionato per t<1, come si vede nella figura seguente.


5

Osserviamo che un problema malcondizionato non si può risolvere numericamente con una accuratezza accettabile, ma un problema bencondizionato non assicura una risoluzione numerica “facile”.

Discretizzazione di un problema a valori iniziali I metodi che saranno analizzati si basano sulla discretizzazione del problema, ossia:

si discretizza l’intervallo [t0,T] in un insieme di nodi o punti di griglia t0 =t0<x1<…<tn=T, tk= t0+kh, k=0,...,n, h=(T-t0)/n passo

si discretizza la ODE, ossia si sostituisce al problema continuo un problema discreto

si determinano i valori yk, approssimazioni dei valori y(xk).

Nota la soluzione in tali punti,sarà possibile ricostruire un’approssimazione su tutto l’intervallo [t0,T] tramite un’interpolazione polinomiale o con funzioni splines. A partire da y0, in ogni nodo la soluzione è approssimata da un valore yk ottenuto in funzione di valori precedenti. I metodi in cui il valore attuale dipende solo dal precedente sono detti one-step Si costruisce, a partire dal valore iniziale, l’insieme di punti yk cercando di seguire la traiettoria della curva soluzione:

Ovviamente yk≠y(tk), non appartiene alla curva soluzione passante per t0, ma ad un’altra curva, è soluzione esatta di un problema con stessa equazione e diversa condizione iniziale. Ad ogni passo una perturbazione, (errore di troncamento e di round-off), produce il passaggio da una

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

1

2

3

4

5

6sol. per c=1,2,3,4,5,6,alfa=5

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

1

2

3

4

5

6sol. per c=1,2,3,4,5,6,alfa=5

(t0,y0)

(ti,yi)

(tn,yn)


6

curva ad un’altra, per cui, se il problema è malcondizionato, i valori calcolati si allontanano da quelli veri e l’errore si amplifica, se il problema è bencondizionato l’errore non si amplifica.

Metodo di Eulero esplicito

Il metodo di Eulero esplicito è il metodo di discretizzazione più semplice ma non è utilizzato in pratica perché poco accurato ed efficiente. Esso è però utile per comprendere i fenomeni (convergenza,stabilità,..) connessi ad un metodo numerico per la risoluzione di una ODE. Si considerano i primi due termini della serie di Taylor di y(t): y(t+h) = y(t)+h y’(t) +y’’(z)h2/2 = y(t)+h f(t,y(t)) +y’’(z)h2/2 z∈(t,t+h) e si trascura il resto. In questo modo si ottiene una discretizzazione del modello matematico, ossia il problema continuo 2) si trasforma nel problema discreto : yk+1= yk+h f(tk,yk) k=0,1,…,n-1. Ciò equivale ad approssimare la derivata con il rapporto incrementale: y’(tk) = f(tk,yk)≈ [y(tk+h)-y(tk)]/h= (yk+1- yk)/h A partire da y0 si calcola y1 , stima di y(t1), con la formula: y1= y0 + f(t0,y0) a partire da y1 si calcola y2 , stima di y(t2), con la formula: y2= y1 + f(t1,y1)

e così via fino a T. Tale metodo è detto esplicito o in avanti perchè la soluzione numerica dipende solo dal valore calcolato precedentemente e perché si ottiene discretizzando la derivata con un rapporto incrementale in avanti. Geometricamente yi+1 è sulla retta tangente in (t i,yi) la curva soluzione passante per (ti,yi).


7

In Matlab una function che implementa il metodo di Eulero esplicito è: function[t,u]=eulero(t0,T,h,f,y0) %risolve l’equazione differenziale y’=f(t,y(t)) con Eulero in [to,T] con passo h e valore iniziale y0 %in output si ha il vettore t dei punti e u della soluzione approssimata nei punti t t=[t0:h:T];u(1)=y0; for i=1:length(t)-1 y=u; u=[u,y(i)+h*f(t(i),y(i))]; end La funzione f si può definire tramite un function file e si richiama con un function handle, o tramite una anonimous function . Esempio 5. y’(t)=-10(t-1)y(t), t∈[0,2], con y(0)=e-5 e soluzione y(t)=e-5(t-1)^2. Definiamo f(t,y(t)) tramite la function( con due parametri di input): function x=f(t,y) x = -10.*(t-1).*y ; Chiamiamo la function eulero : >>[t,u]=eulero(0,2,0.1,@f,exp(-5)); Grafico dei punti, della soluzione approssimata interpolando con una spline e della soluzione vera: >>x=0:0.1:2;y=c*exp(-5*(x-1).^2); >>e=spline(t,u,x); %spline cubica interpolante not-a-knot >>plot(t,u,'o',x,y,'r',x,e,'g')

ti ti+1

y(ti+1)

y(ti)yi+1


8

La risoluzione di un sistema di equazioni con il metodo di Eulero esplicito (e con gli altri metodi che vedremo in seguito) si ottiene in Matlab semplicemente vettorizzando le formule. In Matlab una function che implementa il metodo di Eulero esplicito relativamente ad un sistema è function [t,y] = eulsist(fun,t0,tf,y0,h) % Risolve con il metodo di Eulero un sistema di equazioni differenziali definito da y'=fun(t,y(t)) %in[T0 ,TF] con passo h e y0 vettore di valori iniziali % output : vettore t dei punti , y matrice con numero di colonne=dimensione del % sistema,colonne= soluzioni approssimate nei punti t % La funzione fun deve ritornare un vettore della stessa lunghezza di y0 t=[t0:h:tf]; y(1,:)=y0; for k=1:length(t)-1 k1=fun(t(k),y(k,:)); y(k+1,:)=y(k,:)+h*k1; end Risolviamo ad esempio il problema predatore preda (volpi-conigli), con a=0.01, r0=300, f0=150 . Definiamo f tramite il function file:

function f=LotVol(t,y) f=[2*y(1)-0.001*y(1)*y(2) -y(2)+0.01*y(1)*y(2)]; Risolviamo il sistema : >> y0=[300; 150]; >> [t,y]=eulsist(@LotVol,0,50,y0,0.005); >> plot(t,y(:,1),t,y(:,2)); >> legend('conigli','volpi');xlabel('tempo')

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

sol. vera

Eulero h=0.1

0 5 10 15 20 25 300

100

200

300

400

500

600

coniglivolpi


9

I problemi da affrontare in generale nella risoluzione numerica di una ODE sono:

Stabilità : lo studio di come gli errori introdotti nei calcoli precedenti si amplificano nella soluzione calcolata yn

Convergenza : Un metodo numerico è convergente se, fissato tk, la soluzione numerica tende a quella vera al tendere all’infinito del numero di suddivisioni dell’intervallo, ossia se

limh→0| y(tk)-yk | =0 ∀k=0,..n Ciò vuol dire che si può ottenere nel risultato l’accuratezza voluta, scegliendo un h opportunamente piccolo. Errore e stabilità della soluzione numerica di una ODE Per studiare stabilità e convergenza bisogna quindi ottenere una stima dell’errore globale in ogni punto. Nella soluzione calcolata di un problema di ODE si sommano due tipi di errori:

Errore di discretizzazione o troncamento : dipende dal metodo numerico scelto per discretizzare la ODE

Errore di round-off : dipende dal sistema aritmetico floating-point utilizzato. L’errore di troncamento dipende da due fattori: dallo schema di discretizzazione scelto, che introduce un errore nei singoli passi, e dall’accumulo e propagazione di tali errori nel valore corrente calcolato, ossia si definiscono: errore locale e globale di troncamento. Errore locale di troncamento : eT

k+1= y(tk+1)-yk+1, dove y(tk)=yk è l’errore introdotto in un singolo passo di applicazione dello schema, e si valuta supponendo che i valori calcolati al passo k siano esatti. E’ cioè una stima della bontà dello schema numerico. Errore globale di troncamento: Ek+1= y(tk+1)-yk+1, dove y(tk)≠yk È l’errore commesso nell’approssimare il valore esatto y(tk) con yk, a cui contribuiscono gli errori locali dei singoli passi e la loro propagazione (l’errore globale non è semplicemente la somma degli errori locali).

Ciò che si vuole in definitiva è un errore globale piccolo, ma ciò che si riesce a stimare e controllare in generale è l’errore locale. Un altro concetto fondamentale riguardante l’accuratezza è l’ordine di un metodo.

y(ti)

yi+1

ti ti+1

y(ti+1)

Errore locale

t0 t1 t2

Errore locale

Errore globale


10

Un metodo è di ordine p se esiste una costante M∈R tale che: eT

k+1≅ Mhp+1= O(hp+1) con M che dipende da f e dalle sue derivate, dall’intervallo [t0,T] ma è indipendente da h. In genere ciò implica che l’errore globale è O(hp), ossia , dimezzando il passo, l’errore globale diminuisce di un fattore pari a 2p. Errore locale e globale del metodo di Eulero esplicito Sia y’’(t) continua in [t0,T]. Supponendo y(tk)=yk il metodo di Eulero Esplicito è : yk+1= yk+h f(tk,yk) = y(tk)+h f(tk,y(tk)) Sviluppando y(tk+1) in formula di Taylor : y(tk+1) = y(tk)+h f(tk,y(tk)) +y’’(z)h2/2 z∈(tk,tk+1) Si ha : eT

k+1= y(tk+1)-yk+1= y’’(z)h2/2= O(h2)

quindi il metodo è di ordine p= 1 (O(hp+1) = O(h2) ) e si ha che: limh→0 eT

k+1= 0. Un metodo per il quale vale tale relazione è detto consistente, cioè lo schema numerico è una “buona” approssimazione del modello continuo. Per valutare l’errore globale bisogna considerare y(tk)≠yk, k>0 Si ha: y(tk+1) = y(tk)+h f(tk,y(tk)) +y’’(z)h2/2 z∈(tk,tk+1) e yk+1= yk+h f(tk,yk) con y(tk)≠yk, k>0 Sottraendo membro a membro si ottiene :

Ek+1= y(tk+1)-yk+1= y(tk)-yk+ h [f(tk,y(tk))-f(tk,yk)] + y’’(z)h2/2= Ek + h [f(tk,y(tk))-f(tk,yk)] + eT

k+1

Applicando il teorema del valor medio di Lagrange, si ha f(tk,y(tk))-f(tk,yk) = fy(tk,ξ)[ y(tk)-yk]= JEk ξ∈ (y(tk)-yk) da cui si la la seguente relazione tra errore locale e globale:

Ek+1 = (1+hJ)Ek + eTk+1 = errore di propagazione + errore locale

Si ha che (1+hJ)Ek rappresenta l’accumulo degli errori precedenti, mentre (1+hJ) rappresenta il fattore di amplificazione dell’errore globale Ek nel passare da tk ad tk+1 e quindi se:

| 1+hJ| < 1 ⇒ l’errore diminuisce ⇔ metodo di Eulero esplicito stabile |1+hJ| > 1 ⇒ l’errore aumenta ⇔ metodo di Eulero esplicito instabile

Se il problema è malcondizionato si ha fy = J>0 ⇒ | 1+hJ| > 1 sempre ⇒ Eulero instabile


11

Se invece il problema è bencondizionato si ha: fy = J<0 ⇒ | 1+hJ| < 1 ⇔ h < -2/J ⇒ Eulero stabile ⇔ 0 < h <2/|J| Quindi, se il problema è bencondizionato, il metodo di Eulero esplicito è condizionatamente stabile. Esempio 6. Consideriamo il problema: y’(t)=-10y(t) t∈[0,1] y(0)=1 si ha J= -10, per cui il metodo di Eulero esplicito è stabile se h< 2/10=0.2. Nella figura vi sono i risultati per h=0.15(stabile) e h=0.25 (instabile).

In generale un metodo si dice (assolutamente) stabile per un certo passo h ed un dato problema differenziale, se una perturbazione in un valore yk non si amplifica nei valori successivi. Poiché in tale definizione vi è dipendenza dal problema, per evitare ciò (e per i metodi più complessi di Eulero), si studia l’assoluta stabilità rispetto al problema test:

y’(t)=λy(t) t∈[0,T] y(0)=y0

che ha soluzione y(t)=y0eλt. Se il metodo risulta stabile per il problema test quando: a<hλ<b allora risulterà stabile, per una ODE generale y’=f(t,y), quando: a<hfy=hJ<b Se non si hanno informazioni su y o fy la determinazione di h non è agevole. Conviene seguire, come vedremo in seguito, una procedura adattativa per il calcolo automatico dell’ h ottimale. Si definisce la :

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-4

-3

-2

-1

0

1

2

3

4

5

6

sol vera h=0.25

h=0.15

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-4

-3

-2

-1

0

1

2

3

4

5

6

sol vera h=0.25

h=0.15


12

Regione di assoluta stabilità : insieme dei valori h>0 e λ tali che una perturbazione in yk non si amplifica nei valori successivi. Per il metodo di Eulero esplicito la regione di assoluta stabilità è data dai valori di h e |λ| tali che h|λ|<2. Osserviamo che più è grande λ più deve essere piccolo h per garantire la stabilità. Convergenza del metodo di Eulero esplicito Per verificare la convergenza bisogna dimostrare che: limh→0Ek = limk→∞Ek= 0 Poiché si ha: Ek+1 = (1+hJ)E k + eT

k+1 il metodo converge se:

eTk+1 tende a zero con h: consistenza

e l’errore di propagazione non “esplode” per h che diventa piccolo, ossia se la regione di assoluta stabilità contiene l’origine: stabilità. Ciò significa che esiste un h’ tale che,per ogni h, 0<h<h’, l’errore si mantiene limitato. Per il metodo di Eulero esplicito si dimostra il seguente: Teorema di convergenza. Sia |fy|< L<∞ in [t0,T] ed y(t) con derivata seconda continua e tale che: |y’’|<K<∞, posto s=tk-t0,si ha : |Ek| < (esL-1) [E0 + hK]/2L Se E0→0, per tk fissato, si ha che yk→y(tk), ossia: |Ek| = | yk-y(tk)| = O(h) ∀k=0,..n Il metodo converge, ossia si può ottenere l’accuratezza voluta scegliendo opportunamente h. Fissata cioè una tolleranza TOL, esiste un h abbastanza piccolo tale che |Ek|<TOL. Il metodo di Eulero esplicito è di ordine 1, ossia, dimezzando h, l’errore diminuisce di un fattore pari a 1/2. Consideriamo l’equazione dell’esempio 2 in [0,20] con K=0.001, pmax=1000,p0=1. Ad un passo più piccolo corrisponde una soluzione approssimata migliore (l’errore si dimezza con h).


13

Esempio 7. Consideriamo l’equazione: y’(t)=-2y(t) tє[0,4] y(0) =1 soluzione vera e-2t applichiamo il metodo di Eulero dimezzando progressivamente h, e calcoliamo poi l’errore nell’ultimo punto: %test eulero convergenza h=0.5;num=[];err=[]; for j=1:10 [t,u]=eulero(0,4,h,@feul,1); num=[num;length(t)]; %errore nell’ultimo punto er=abs(exp(-2*t(end))-u(end)); err=[err;er]; h=h/2; end Osserviamo il numero eccessivo di punti per ottenere circa 6 cifre corrette. Nei risultati ottenuti finora non si è considerato l’errore di round-off che interviene nei calcoli. Se l’errore di round-off ad ogni passo è ≅ ε, allora l’errore totale di round-off sarà: Er.o.=nε = ε (T-t0)/h ⇒ errore totale=O(Mh+ε/h) Ciò significa che un h troppo piccolo può condurre ad un accumulo tale di errori di round-off per cui la soluzione calcolata può divergere. Esiste un h ottimale(difficile da valutare) al di sotto del quale l’errore di round-off dà il contributo dominante all’errore totale:h ott= minh(Mh+ε/h).

Per quanto riguarda la complessità computazionale,si misura in base al numero di valutazioni di funzione. Il costo totale nel calcolo di yend a partire da y0 è il costo del singolo passo per il numero di passi, e quindi dipende da h, ossia dall’ordine del metodo e dalla regione di stabilità. Si vuole h il più grande possibile per ridurre l’errore di round-off e il costo computazionale ma h deve essere tale da:

0.01 0.51 1.01 1.51 2.01 2.510

5

10

h

erro

re

errore totale

troncam . round -off

h ott. 0.01 0.51 1.01 1.51 2.01 2.510

5

10

h

erro

re

errore totale


h ott. 0.01 0.51 1.01 1.51 2.01 2.510

5

10

h

erro

re

errore totale


h ott. 0.01 0.51 1.01 1.51 2.01 2.510

5

10

h

erro

re

errore totale


h ott.

Errore relativo Numero di punti3.3546e-004 93.2020e-004 172.3501e-004 331.4114e-004 657.7044e-005 1294.0206e-005 2572.0532e-005 5131.0374e-005 10255.2143e-006 20492.6140e-006 4097


14

soddisfare l’accuratezza richiesta garantire la stabilità

Il metodo di Eulero Esplicito è: Convergente di ordine 1 Condizionatamente stabile : h|λ|<2

Si può ottenere nel risultato l’accuratezza voluta scegliendo opportunamente h, ma l’ordine è basso e ciò comporta un numero eccessivo di punti per ottenere una buona accuratezza(esempio 7), inoltre la scelta del passo h è vincolata dalla stabilità, e spesso h risulta impraticabile. Vi è quindi la necessità di metodi che consentono un h più grande, ossia

• metodi con una regione di stabilità più ampia • metodi più accurati (ordine>1)

Metodo di Eulero Implicito (BE) In maniera analoga al metodo di Eulero esplicito, si può discretizzare la derivata con il rapporto incrementale all’indietro: y’(tk+1) = f(tk+1,yk+1)≈ [y(tk+1)-y(tk+1-h)]/h= (yk+1- yk)/h ottenendo il metodo di Eulero implicito (BE): yk+1= yk+h f(tk+1,yk+1) k=0,1,…,n-1. Tale metodo è detto implicito poichè la soluzione numerica dipende, attraverso f(tk+1,yk+1), anche da se stessa. Consideriamo l’equazione: y’(t)=-y3 y(0)=1 Applichiamo il primo passo di BE: y1 = y0 +hf(t1,y1)=1-0.5y1

3

otteniamo la seguente equazione non lineare in y1 : y1 -1-0.5y1

3 = 0 e così via. Analizziamo la stabilità del metodo BE utilizzando l’equazione test:

y’(t)=λy(t) t∈[0,T] y(0)=y0

che ha soluzione y(t)=y0eλt. Si ha: yk+1= yk+h λyk+1 k=0,1,…,n-1. Da cui: yk+1= yk[1/(1-h λ)] Il fattore di amplificazione da un passo all’altro è 1/(1-hλ). Se λ<1, ossia il problema è bencondizionato, si ha che: 1/(1-hλ)<1 sempre


15

e quindi il metodo è incondizionatamente stabile, ossia è stabile per ogni valore di h. Ciò vuol dire che l’accuratezza richiesta è l’unica condizione sulla scelta di h. Per quanto riguarda la convergenza valgono gli stessi risultati del metodo esplicito, si dimostra,sotto opportune ipotesi, che , se E0→0, per tk fissato, |Ek| = | yk-y(tk)| = O(h) Ossia converge ed è di ordine 1. I metodi impliciti in generale sono più costosi di quelli espliciti poiché richiedono ad ogni passo, se f non è lineare rispetto ad y, la risoluzione di un’equazione non lineare, ma hanno proprietà di stabilità ottimali e per alcuni problemi (problemi stiff) sono gli unici computazionalmente praticabili. Per risolvere l’equazione non lineare si può usare la function matlab fzero: x = fzero(@fun,x0), che tenta di determinare l’approssimazione dello zero di fun in prossimità del valore iniziale x0. Se fun è una funzione parametrica(a più parametri di input) la chiamata è: x = fzero(@(x) fun(x,t,..),x0) dove @(x) identifica la variabile rispetto a cui determinare lo zero. Nel caso di Eulero implicito come valore iniziale si utilizza l’ultimo valore disponibile della soluzione approssimata. Analogamente per la risoluzione di un sistema di equazioni differenziali si usa la function matlab fsolve che risolve un sistema di equazioni non lineari. In Matlab una function che implementa il metodo di Eulero implicito è: function[t,u]=eulback(t0,T,h,f,y0) %risolve un’equazione differenziale y’=f(t,y(t)) con BE in [to,T],con passo h e valore iniziale y0 % output: vettore t dei punti e u della soluzione approssimata,chiama fzero per la risoluzione %dell’eq.non lineare(@(b) identifica in fzero la variabile %rispetto a cui determinare lo zero t=[t0:h:T];u(1)=y0; for i=2:length(t) tx=t(i);x=u(i-1); sol=fzero(@(b) b-x-h.*f(tx,b),x); u=[u,sol]; end Risolviamo con Eulero Esplicito ed Implicito l’equazione : y'(t)=-2y(t) tє[0,10] y(0)=0 che ha come soluzione y(t)=e-2t. La condizione di stabilità per Eulero esplicito impone h<2/│λ│=1, mentre per Eulero Implicito non c’è nessuna limitazione su h. Definiamo la function che definisce f(t,y): function x=f(t,y) x=-2*y; e chiamiamo le function eulero ed eulback con passo h=0.9(fig. a sinistra) e h= 1.1(fig. a destra) :


16

Si osserva che, mentre Eulero Implicito con entrambi i valori di h è stabile , anche se ovviamente per il valore minore di h si ha un’approssimazione migliore, Eulero Esplicito per h=1.1 è instabile e quindi la soluzione calcolata “esplode” , mentre per h=0.9 è stabile ma troppo prossimo alla soglia di stabilità, per cui l’errore di propagazione fa sì che l’errore globale sia molto grande. Per avere una buona soluzione con Eulero Esplicito bisogna prendere un h più piccolo, ossia più punti di discretizzazione, ad esempio h=0.4. Nella figura seguente vi sono le soluzioni calcolate con eulero con h=0.4 e eulback con h=0.9.

Infine effettuiamo un confronto sui tempi tra i due metodi: %script file tempi che effettua un confronto sui tempi temp=cputime; [t,u]=eulero(0,1,0.001,@f,0); temp=cputime-temp; temp1=cputime; [t,u]=eulback(0,1,0.001,@f,0); temp1=cputime-temp1; >> temp temp = 0.015625000000000 >> temp1 temp1 = 0.562500000000000

0 1 2 3 4 5 6 7 8 9 10-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

sol.veraBEEE

h=0.9

0 1 2 3 4 5 6 7 8 9 10-6

-4

-2

0

2

4

6

sol.veraBEEE

h=1.1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

euleroeulbacksoluzione


17

Come era prevedibile BE è più costoso di Eulero Esplicito, ma, come vedremo, per alcuni problemi risulta l’unico tipo di metodo praticabile.

Metodi Runge-Kutta Poiché il metodo di Eulero si ottiene considerando i primi due termini della serie di Taylor, si possono ottenere metodi di ordine superiore considerando più termini della serie: y(t+h) = y(t)+h y’(t) +y’’(t)h2/2+...+ y(p)(t)hp/p! + O(hp+1) In tal modo si ottiene un metodo di ordine p, ossia, dimezzando il passo, l’errore diminuisce di un fattore pari a 1/2p, il che equivale ad ottenere la stessa accuratezza del metodo di Eulero con un h più grande. Il problema di tali metodi da un punto di vista computazionale, che li rende impraticabili, è la valutazione delle derivate di f. I metodi Runge-Kutta evitano la necessità di valutare le derivate di f e nel contempo hanno la stessa accuratezza di un metodo di Taylor, risultano semplici da implementare e sono per tali motivi tra i metodi più utilizzati per risolvere un problema a valori iniziali. Si basano sull’idea di approssimare la derivata, noto yk, con una combinazione lineare F(tk,yk,h) di valori di f(t,y) in opportuni punti di [tk,tk+h] ossia: yk+1 = yk+h F(tk,yk,h) i coefficienti di tale combinazione lineare ed i punti sono determinati in modo da avere un errore locale eT

k+1= O(hp+1). Ciò si ottiene imponendo che lo sviluppo in serie di Taylor di yk+1 coincida con lo sviluppo di Taylor di y(tk+1), di punto iniziale tk, fino ai termini di ordine O(hp). Runge-Kutta di ordine 2 RK2(metodo di Heun). Il metodo consiste in: yk+1 = yk+ h/2[f(tk,yk)+f(tk+h,yk+ hf(tk,yk))] tale che eT

k+1= O(h3). Esempio 8.Consideriamo il problema: y’(t)=-y(t)+2cos(t) t∈[0,10] y(0)=1 (soluzione y(t)=sin(t)+cos(t)). In Figura vi sono le soluzioni ottenute con Eulero e RK2 con h=0.5.

0 1 2 3 4 5 6 7 8 9 10-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

Eulero

RK2

sol.vera

0 1 2 3 4 5 6 7 8 9 10-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

Eulero

RK2

sol.vera


18

Il metodo RK2 ha ordine maggiore di Eulero ma richiede due valutazioni di f. Confrontando il numero di cifre esatte ≅ -log10|yn-y(xn)| si ha:

numero cifre esatte h Eulero RK2 0.5 0.05

0.77481.85

1.23683.2813

Ossia, dividendo h per 10 Eulero guadagna ≅ 1 cifra e RK2 ≅ 2 cifre. Il metodo più utilizzato, che realizza un compromesso accettabile tra numero di valutazioni di f ed accuratezza, è Runge-Kutta a 4 stadi (RK4), di ordine 4, ossia con eT

k+1= O(h5).Il metodo è: yk+1 = yk+(1/6)[/K1+2K2 +2K3 +K4] con: K1 = h f(tk,yk) K2 = h f(tk+h/2,yk+K1/2) K3 = h f(tk+h/2,yk+K2/2) K4 = h f(tk+h,yk+K3)

Per studiare la stabilità ci si riferisce al problema test e si ottiene:

RK2 RK4

h|λ|<2 h|λ|<2.7853

La regione di stabilità tende quindi ad ampliarsi al crescere dell’ordine. Per la convergenza dei metodi Runge-Kutta, si può dimostrare un teorema analogo a quello relativo al metodo di Eulero, ossia, sotto opportune ipotesi di regolarità su f, un metodo RK di ordine p converge. Confrontiamo il metodo di Eulero esplicito, RK2 ed RK4, relativamente all’esempio 7, a parità di costo, ossia di numero di valutazioni di f.

Da tale confronto si evidenzia che l’ordine più elevato bilancia il maggior numero di valutazioni richieste, nel senso che, a parità di costo, si ottiene una maggiore accuratezza. Analogamente al metodo di Eulero, anche per i metodi Runge-Kutta, l’errore totale cresce a causa dell’errore di round-off. Il valore di hott dipende dall’ordine del metodo, ed aumenta all’aumentare dell’ordine. Appare quindi evidente il vantaggio nell’utilizzo di formule di ordine più elevato. Stima dell’errore e del passo Un software efficiente per la risoluzione di una ODE deve fornire una soluzione numerica “accurata” (con un’accuratezza TOL prefissata in input) ed efficiente (costo computazionale

t Eu., h=0.125 RK2, h=0.25 RK4, h=0.5 Sol.vera 4 -1.45975 -1.39362 -1.40983 -1.41044 6 0.66338 0.68096 0.68085 0.68075 8 0.90760 0.82697 0.84315 0.84385 10 -1.41878 -1.36923 -1.38261 -1.38309


19

accettabile). E’ quindi necessario ottenere una stima dell’errore, per verificare la bontà della soluzione calcolata rispetto alla tolleranza richiesta. Le formule analitiche dell’errore locale e globale non sono utilizzabili nella pratica; una tecnica attuabile consiste nello stimare l’errore locale come differenza dei risultati ottenuti con due metodi di ordine diverso. Le coppie di formule più usate sono del tipo RK23 e RK45. Tale tecnica è efficiente se i metodi usano valutazioni di funzione in comune, come ad esempio le formule di Runge-Kutta Fehlberg (RKF), tra cui RKF45, che usa 6 valutazioni di funzione per ottenere una formula di ordine 5, di tali valutazioni poi 4 sono usate per ottenere una formula di ordine 4 .L’errore locale è stimato dalla differenza delle approssimazioni ottenute con le due formule: eT

k+1= |zk+1-yk+1| . Il numero di valutazioni di funzione per ogni passo è quindi 6. Più precisamente le formule RKF45 sono:

yk+1=yk+[ (25/216)K1+(1408/2565)K3+(2197/4104)K4-(1/5)K5] zk+1=zk+[(16/135)K1+(6656/12825)K3+(28561/56430)K4-(9/50)K5+(2/55)K6] eT

k+1≅|(1/360)K1-(128/4275)K3-(2197/75240)K4+(1/50)K5+(2/55)K6| dove: K1 = h f(tk,yk) K2 = h f(tk+h/4,yk+K1/4) K3 = h f(tk+3h/8, yk+(3K1+9K2)/32)) K4 = h f(tk+12h/13, yk+(1932K1-7200K2+7296K3)/2197) K5= h f(tk+h, yk+439K1/216-8K2+3680K3/513-845K4/4104) K6= h f(tk+h/2, yk-8K1/27+2K2-3544K3/2565+1859K4/4104-11K5/40) Un problema fondamentale è la scelta del passo, che, come si è visto, dipende dall’ordine del metodo, dalla relativa regione di stabilità e dall’accuratezza richiesta, ma con un passo troppo piccolo si ha un costo computazionale elevato ed un eccessivo accumulo dell’errore di round-off. Si vuole cioè determinare h il più grande possibile che consente di ottenere la soluzione entro la tolleranza fissata. Spesso un passo costante non è conveniente, in quanto, in certi problemi, l’errore è tale che un h buono in una zona dell’intervallo non lo è altrove. E’ opportuno quindi avere la possibilità di variare h in maniera adattativa, così da renderlo sufficientemente piccolo solo dove è necessario. Scelto un passo h iniziale ed un intervallo realistico di variabilità per hk, si valuta una stima dell’errore e si modifica hk così da soddisfare il test: |eT

k|< TOL hk errore per unit step In pratica, se l’errore è molto più piccolo del richiesto, si può aumentare hk. Se poi il test non è verificato si diminuisce hk. Per calcolare un hk ottimale in maniera automatica, si supponga di usare un metodo di ordine p, ossia: |eT

k|≅M hk p+1 Se |eT

k|> TOL hk bisogna determinare un nuovo passo, hk=q hk , tale che: |eT

k|≅TOL hk ≅ TOLq hk Poiche: eT

k ≅ M hkp+1=Mqp+1 hk p+1 ≅ qp+1 eT

k si ha: TOLq hk ≅ qp+1 eT

k ⇔ q = [TOL hk / eTk]1/p


20

Tale stima del passo è utilizzata per calcolare il nuovo passo e la nuova approssimazione, se il test è verificato, altrimenti si riduce ancora h e si calcola una nuova stima dell’errore. La caratteristica fondamentale di un metodo a passo variabile è che automaticamente genera il passo adatto alla stabilità. Algoritmo RKF45 t = t0, y = y0, hmax = (T-t)/a, hmin = (T-t)/b, h = (T-t)/c while (t<T and h>hmin) if t+h>T then h=T-t endif calcola Ki , i=1,...,6 calcola R=eT

k if R < TOLh then t=t+h calcola y nuova approssimazione endif if R≠0 then q= 0.84[TOL h/ R]1/4 h=min(hmax, qh) endif endwhile if(t<T) ‘probabile singolarità’ endif Nell’algoritmo precedente bisogna definire le costanti a,b,c, (ad esempio potrebbero essere a = 5, b = 20000, c=100). Si è corretta la stima di q con 0.84 per renderla meno conservativa e nell’ultimo if si tiene conto della possibilità che h diventi minore di hmin, il che significa che ci può essere una singolarità o l’equazione ha particolari caratteristiche (è stiff) e va trattata con metodi opportuni. Nelle figure seguenti vi è la soluzione di y’=-30y con ode45 del matlab(RK45) con TOL=10-4 e TOL=10-7.

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1RKF45 Tol=10 -7

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1RKF45 Tol=0.0001


21

Problemi stiff Una equazione (o un sistema) differenziale si dice stiff se la sua soluzione è composta da un termine che non varia molto con t (stazionario) e da un termine che, in un intervallo limitato di tempo, ha una grande variazione per poi decadere rapidamente a zero (transiente). Il termine stiff è stato introdotto nel 1952 a proposito dello studio di un sistema di molle delle quali alcune molto più rigide (“stiff”) di altre. Le equazioni stiff descrivono un sistema fisico caratterizzato da costanti di tempo (velocità di decadimento) molto diverse tra loro e in cui la scala temporale è piccola se paragonata all’intervallo di integrazione. Equazioni e problemi stiff si devono risolvere in molti problemi fisici e chimici (analisi di circuiti, problemi chimici, calcoli relativi a reattori nucleari,…). Esempio 9. Il seguente problema (bencondizionato) y' = -100(y-sint)+cost y(0) = 1 ha come soluzione : y(t)=sint+e-100t costituita dal termine stazionario sint e dal termine transiente e-100t . In generale i metodi espliciti approssimano bene il termine stazionario, ma l’errore associato al termine transiente è dominante e provoca effetti disastrosi. Nella zona relativa al transiente, l’accuratezza richiede un passo piccolo, poi si potrebbe utilizzare un passo più grande, relativo alla componente stazionaria, ma, per assicurare la stabilità, il passo h è limitato proprio dalla componente transiente ad assumere valori molto piccoli. Infatti in generale un problema stiff è tale che fy << 0, il che comporta un h molto piccolo per garantire la stabilità, e quindi una complessità e conseguentemente un tempo di calcolo eccessivi. La componente transiente quindi, anche se non dà più contributo alla soluzione quasi subito, influenza però in modo decisivo l’intervallo di stabilità. Non si può dare una definizione di ODE stiff, non è una proprietà dell’equazione, ma un comportamento che si evidenzia quando l’ODE è risolta numericamente e che dipende dall’equazione, dall’intervallo e dal metodo numerico. Una ODE stiff non è “difficile” da risolvere se si sceglie il metodo adatto. Nell’equazione dell’esempio per assicurare la stabilità deve essere h<2/100 per cui l’equazione risulta stiff in un intervallo abbastanza grande rispetto ad h. In pratica, quando si cerca di risolvere un’equazione stiff l’instabilità generata da un h non sufficientemente piccolo genera un grande errore locale, da cui la richiesta di un h sempre più piccolo, che indica in genere che il problema è stiff. Spesso, per un’equazione stiff, il passo h che assicura la stabilità è computazionalmente impraticabile perché troppo piccolo e conviene usare un metodo implicito incondizionatamente stabile, ossia senza limitazioni su h (dipendente solo dall’accuratezza richiesta), anche se ad ogni passo può richiedere la soluzione di un’equazione non lineare. Nelle figure seguenti sono rappresentate, insieme alla soluzione, le approssimazioni calcolate con Eulero implicito ed esplicito a parità di passo h=0.02 relativamente all’equazione dell’esempio 9 in [0,0.5].


22

Si vede che Eulero esplicito, a differenza di BE, che approssima bene la soluzione, ha un andamento altamente oscillatorio. Osserviamo infine che tutti i metodi descritti si applicano ad un sistema di equazioni differenziali, semplicemente vettorializzando le formule risolutive. Equazioni Differenziali in Matlab Il Matlab possiede numerose funzioni (ODEsuite) con interfaccia omogenea che risolvono un sistema di ODE, anche relative a sistemi stiff. Le function presenti in Matlab sono: per equazioni non stiff:

ODE45 metodo RK45 a passo variabile ODE23 metodo RK23 a passo variabile

per equazioni stiff i metodi impliciti:

ODE15s metodo multistep ODE23s metodo ad un passo

Una definizione pratica di stiffness: se un problema è “difficile” (tempo di calcolo eccessivo o warning su h troppo piccolo) per ode45 e “facile” per ode15s è stiff. Osserviamo che se un problema è difficile per entrambi, non vuol dire che è stiff, ma che bisogna usare un metodo più sofisticato di ode45(ode113). In generale per risolvere un’equazione differenziale (o un sistema ) in Matlab bisogna:

Scrivere una function con due parametri di ingresso t ed y che valuta la funzione f(t,y) Applicare un solutore adatto al problema Effettuare un grafico dei risultati

0 0.1 0,2 0,3 0,4 0,50

0.2

0.4

0.6

0.8

1metodo BE

0 0.1 0.2 0.3 0.4 0.5-2

-1

0

1

2metodo Eulero Esplicito

0 0.1 0,2 0,3 0,4 0,50

0.2

0.4

0.6

0.8

1metodo BE

0 0.1 0.2 0.3 0.4 0.5-2

-1

0

1

2

0 0.1 0,2 0,3 0,4 0,50

0.2

0.4

0.6

0.8

1metodo BE

0 0.1 0.2 0.3 0.4 0.5-2

-1

0

1

2metodo Eulero Esplicito


23

Il comando più semplice per risolvere un’equazione differenziale è: [t,y]=solver(fun,tspan,y0)

solver nome del solutore che si vuole utilizzare input

fun handle di function con due parametri obbligatori : t scalare e y vettore colonna che valuta f(t,y)

tspan vettore contenente gli estremi dell’intervallo di integrazione, se =[t0<t1<..<tm] si ha la soluzione in punti fissati dall’utente

y0 vettore di valori iniziali

output t vettore colonna contenente i punti generati dalla discretizzazione y array contenente la soluzione, la riga i è la componente di y all’ istante t(i)

Ad esempio: Richiamiamo l’ode solver con un handle @ alla function, >> [t,z]=ode45(@fun,[0 10],1); >>plot(t,z,'*') >>title('sol eq.diff ') function yp=fun(t,y) yp=-y+2.*cos(t) ;

Le function dell’ODEsuite prevedono argomenti addizionali nella chiamata: [t,y,t1,y1]=solver(fun,tspan,y0,options,p1,p2,…)

options permette di modificare alcuni parametri di default di solver p1,p2,… parametri opzionali che solver può passare a fun e alle funzioni specificate in options

t1 parametro opzionale generato dalla function events y1 parametro opzionale generato dalla function events

Per definire le options si usa il comando odeset. Digitando >>odeset si ha l’elenco dei valori di default ed il nome delle variabili che si possono modificare, tra cui

0 2 4 6 8 10 12 14-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

sol.verasol. calcolata

0 2 4 6 8 10 12 14-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

sol.verasol. calcolata


24

RelTol 1.e-3 tolleranza errore relativo AbsTol 1.e-6 tolleranza errore assoluto MaxStep |tspan|/10 valore max di h InitialStep calcolato passo iniziale scelto OutputFcn function controllo output Events function Il test di accuratezza utilizzato dai solver è: e(i) <= max(RelTol*abs(y(i)),AbsTol) Spesso per problemi reali di una certa complessità i valori di RelTol,AbsTol di default sono troppo grandi, ossia bisogna richiedere più cifre corrette. Per modificare il valore di default delle variabili di odeset: options = odeset(‘nome1’,val1,’nome2’,val2,…) Ciò vuol dire che la tolleranza relativa controlla l’errore, a meno che y(i) non sia troppo “piccolo”. In teoria il valore di default dà circa tre cifre significative corrette, ma può accadere che, se ad esempio l’intervallo temporale è grande, l’accumulo dell’errore di round-off porta ad una soluzione con una accuratezza minore e quindi non affidabile. Esempio 10. Consideriamo l’equazione non lineare del II ordine di Van Der Pool, che rappresenta un circuito RLC in cui la resistenza è sostituita da un elemento attivo (tubo vuoto o semiconduttore): y’’(t)-c(1-y2)y’+y=0 y(0)=2 y’(0)=0 dove c rappresenta lo smorzamento del sistema. La soluzione è un ciclo limite stabile. L’equazione diventa il sistema : x’(t)=y(t) y’(t)=c(1-x2)y-x con x(0)=2, y(0)=0 Fissiamo c=1 = smorzamento debole e creiamo l’ode file : function f=vdpol(t,y) f=[y(2);(1-y(1)^2)*y(2)-y(1)]; Risolviamo il sistema:

variabile default descrizione

>> [t Effe>> p

Speche Ciò esedellesepaNel >> p

Si p funcf=[y( Rich

CE

t,y]=ode45(

ttuiamo il gplot(t,y(:,1),

esso è più

agisce comche si vedempio, se lee fasi è l’oarati). piano delle

plot(y(:,1),y(

uò generali

ction f=vdpo(2);c*(1-y(1

hiamiamo l’

Calcolo NumEquazioni D

(@vdpol,[0

grafico com'-',t,y(:,2),'--

utile visualme paramee è una proe variabili dorbita vista

e fasi (piano(:,2) ,2,0,'*')

izzare la fu

ol(t,y,c) )^2)*y(2)-y

ode solver

0 2-3

-2

-1

0

1

2

3

merico II - CDifferenziali

20],[2;0]);

ponente pe-')

izzare una etro, ossia foiezione stadipendenti rdall’alto, e

o y,y’) evide)

nction pass

(1)];

con c com

4 6 8

CdL SpeciaOrdinarie

er compone

variabile dfare il graficatica a due rappresenta

e ciò ha un

enziando la

sando dei p

me ultimo pa

10 12 14

listica Ing.In

ente:

dipendente co della sodimensioniano le coo

n significato

a condizione

parametri d

arametro:

16 18 20

y

y'

nformatica-

dall’altra, cluzione nel della traietrdinate di u

o fisico più

e iniziale:

dal solver, a

- A.d’Alessi

con la varial piano dellttoria assocun satellite interessan

ad esempio

o

abile indipee fasi (pia

ciata al siste in orbita,

nte dei due

o c:

endente no x,y). ema(ad il piano

e grafici


26

>> [t,y]=ode45(@vdpol,[0 20],[2;0],[ ],1); parametri non presenti Consideriamo ora il caso di c=1000 = smorzamento forte in [0,200], il comando: >> [t,y]=ode45(@vdpol,[0 1000],[2;0],[ ],1000); richiede molto tempo e quindi si può dedurre che il sistema è stiff ossia che y e y’ hanno grandi variazioni in tempi molto diversi. E’ opportuno allora arrestare i calcoli e ricorrere ad una function adatta a tale tipo di problema: >>[t,y]=ode15s(@vdpol,[0 200],[2;0],[ ],1000); Effettuando il grafico in scala semilogaritmica delle due componenti si evidenziano chiaramente le diverse scale temporali caratteristiche di ciascuna componente, da cui la natura stiff del problema: >>semilogy(t,abs(y(:,1)),'-s',t,abs(y(:,2)),'--p'),axis([0 1000 10^-6 10^4])

Esempio 11. Pendolo semplice senza attrito. Consideriamo il moto di un pendolo non lineare di massa m e lunghezza L. Questo problema ad esempio è alla base della dinamica delle macchine automatizzate o del braccio di un robot.

0 100 200 300 400 500 600 700 800 900 100010

-6

10-4

10-2

100

102

104

y

y'

O

θ L

P g


27

Si sposta P dalla posizione iniziale di un angolo θ0 e lo si lascia andare con velocità y0. P si muove lungo un arco del cerchio di centro O, oscillando intorno alla verticale. Il sistema che descrive l’evoluzione dell’angolo θ formato dal pendolo con la verticale è: x'(t)=y y’(t)=(-g/L)sinx dove: x = posizione (=θ) y = velocità angolare g = gravità =9.81 , L = lunghezza (=1) x(0)=θ0 , y(0)=0 (velocità iniziale) Risolviamo il sistema in un intervallo piccolo: >> x0=[0.5 0]; >>[t,x]=ode45(@pend,[0 5],x0); >> plot(t,x(:,1),'r',t,x(:,2),'g') La funzione che definisce il sistema è :

function zprime=pend(t,z) z1=z(2); z2=-9.81*sin(z(1)); zprime=[z1;z2]; Osserviamo che la function deve avere entrambi i parametri di input. Si ottiene :

Il punto P si muove lungo un arco oscillando intorno alla posizione verticale, poichè la velocità iniziale è al di sotto di quella limite (=2√ ), oltre la quale il pendolo compie rotazioni complete. Poiché il pendolo oscilla nel piano delle fasi la traiettoria è chiusa, infatti si ha : >> plot(x(:,1),x(:,2),'LineWidth',2)

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-2

-1.5

-1

-0.5

0

0.5

1

1.5

2


28

Risolviamo il problema nell’intervallo [0,50] con la stessa velocità iniziale e θ0 = π-0.01, ossia il pendolo si trova inizialmente in posizione quasi verticale con velocità nulla, anche in questo caso il pendolo oscilla e la traiettoria nel piano delle fasi è chiusa. Si ottiene invece : >>[t,x]=ode45(@pend,[0 50,[pi-0.01 0]); >> length(t) ans = 709 >> plot(x(:,1),x(:,2),'LineWidth',2) La soluzione calcolata numericamente non è fisicamente accettabile, poiché l’errore è controllato dalla tolleranza di default, in questo caso troppo grande. Modifichiamo i valori di default per la tolleranza usando la struttura options: >> options=odeset('RelTol',1e-6,'AbsTol',1e-7); >> [t,x]=ode45(@pend,[0 50],[pi-0.01 0],options); >> length(t) ans = 2869 >> plot(x(:,1),x(:,2),'-c','LineWidth',2)

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

-10 0 10 20 30 40 50 60 70 80-8

-6

-4

-2

0

2

4

6

8

-4 -3 -2 -1 0 1 2 3 4-8

-6

-4

-2

0

2

4

6

8


29

La soluzione è corretta, e la richiesta di una maggiore accuratezza corrisponde ovviamente ad un aumento del numero di passi. Si pone quindi il problema di stabilire se la soluzione calcolata ha un’accuratezza accettabile, che corrisponde a stabilire un “buon” valore per la tolleranza. La tecnica migliore è quella sperimentale. Si calcola la soluzione con la tolleranza di default RelTol e poi con RelTol/10, si fa il grafico delle curve, se si sovrappongono vuol dire che la soluzione è accettabile, altrimenti si divide la tolleranza fino a che le curve non sono sufficientemente vicine. Il programma seguente effettua questo test. %test tolleranza relativa per il pendolo semplice reltol=10^-3; color='r'; hold on for k=1:2 opts=odeset('RelTol',reltol); [t,y]=ode45(@pend,[0 50],[pi-0.01 0],opts); plot(t,y(:,1),'color',color,'LineWidth',2) reltol=reltol/10; color='g'; end hold off Si ottiene:

La soluzione non è accettabile. Diminuiamo RelTol fino a 10-6 e si ottiene una soluzione accettabile:

0 5 10 15 20 25 30 35 40 45 50-10

0

10

20

30

40

50

60

70

80

RelTol=10-3

RelTol=10-4


30

Finora abbiamo supposto di conoscere l’intervallo [t0,T] in cui è definito il problema, ma in molte applicazioni la determinazione di T è un aspetto fondamentale del problema. Esempio 11. Consideriamo il problema di un corpo che cade soggetto alla forza di gravità e alla resistenza dell’aria. Si vuole calcolare la traiettoria del corpo, e quando il corpo tocca terra. In tal caso si devono arrestare i calcoli. E’ possibile risolvere questo problema definendo, tramite options, una function events di cui determinare lo zero. Il problema dell’esempio è descritto dall’equazione differenziale: y’’(t) = -1+y’2(t) con condizioni iniziali y(0)=1(posizione),y’(0)=0 (velocità). che diventa il sistema: y’1=y2 y’2=-1+y2

2 con y1(0)=1, y2(0)=0. Per quale valore di t y(t)=0? Creiamo la function: function [stop,isterminal,direction] = g(t,y)

stop funzione di cui calcolare lo zero isterminal 1 se si vuole che il solver termini quando stop è 0, 0 altrimenti direction [ ] default (determina tutti gli zeri)

e passiamola, con un handle alla function, al solver tramite options: options=odeset(‘events’,@g) La chiamata al solver diventa: [t,Y,TE,YE] = solver(odefun,tspan,y0,options)

TE tempo in cui la function stop vale zero YE valore della soluzione in TE

Riprendiamo l’esempio, e definiamo la funzione f(t,y): function yd=f(t,y) yd=[y(2);-1+y(2)^2];

0 5 10 15 20 25 30 35 40 45 50-4

-3

-2

-1

0

1

2

3

4


31

Definiamo la function events g: function [ystop,isterminal,direction] = g(t,y) ystop=y(1) ; % funzione di cui determinare lo zero isterminal= 1; direction = [ ] ; A questo punto possiamo calcolare la traiettoria del corpo fino a che non tocca terra: >>options=odeset(‘events’,@g); >>y0=[1;0]; >>[t,y,tfinal,yfinal]=ode45(@f,[0 inf],y0,options); >>tfinal >>yfinal >>plot(t,y(:,1),[0 tfinal],[1 0],’o’) >>axis([-.1 tfinal+.1 -.1 1.1]) >>text(tfinal-0.4,0,[‘tfinal’ = num2str(tfinal)]) Si ottiene: tfinal = 1.6575 yfinal = -0.0000 -0.9300

Esempio 12. Equazioni di Lorentz. Sono equazioni relative al problema delle previsioni metereologiche, che sono difficili ed incerte non per la complessità del problema ma per la natura delle equazioni che lo modellizzano. Tali equazioni appartengono ad una classe di equazioni dette del caos , sono un modello della turbolenza atmosferica all’interno di un accumulo di nubi. In esse un cambiamento piccolo nei valori iniziali produce soluzioni molto

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

0

0.2

0.4

0.6

0.8

1

y

corpo che cade

tfinal =1.6575


32

diverse, per cui l’accuratezza della soluzione dipende in maniera cruciale dalla bontà con cui si misurano le condizioni iniziali. Il meteorologo Lorentz scoprì questo fenomeno nel 1961. Le equazioni seguenti sono un modello semplificato del caos, quelle originali sono molto più complesse. x’(t)= 10(y-x) y’(t)= -xz+28x-y z’(t)=xy-8z/3 con valori iniziali: x(0)=-2,y(0)=-3.5,z(0)=21 Definiamo la function : function f=lorenz(t,x) f=zeros(3,1); f(1)=10*(x(2)-x(1)); f(2)=-x(1)*x(3)+28*x(1)-x(2); f(3)=x(1)*x(2)-8*x(3)/3; Risolviamo il sistema (la soluzione x è una matrice di tre colonne): >>x0=[-2 -3.5 21]; % valori inz. >>[t,x]=ode45(@lorenz,[0 10],x0); >> plot(t,x)

Il grafico nello spazio è >>plot3(x(:,1),x(:,2),x(:,3))

0 1 2 3 4 5 6 7 8 9 10-30

-20

-10

0

10

20

30

40

50equazioni di Lorentz


33

Valutiamo l’effetto di un piccolo cambiamento dei valori iniziali su y(t): >> hold on >> plot(t,x(:,2),'g') >> x0=[-2.04 -3.5 21]; >>[t,x]=ode45(@lorenz,[0 10],x0); >> plot(t,x(:,2),'r')

I due grafici coincidono fino a 1.5 poi si differenziano sempre più. Il grafico nel piano delle fasi è detto grafico butterfly del caos:

-20-10

010

20

-40-20

0

20400

10

20

30

40

50

0 1 2 3 4 5 6 7 8 9 10-25

-20

-15

-10

-5

0

5

10

15

20

25


34

>> plot(x(:,1),x(:,3))

>> plot(x(:,1),x(:,2))

La function seguente produce un’ animazione della soluzione: function muovilorentz(tspan,x0) %anima la soluzione delle equazioni di Lorentz del caos %come una cometa %tspan:intervallo di tempo

-20 -15 -10 -5 0 5 10 15 205

10

15

20

25

30

35

40

45phase plan

-20 -15 -10 -5 0 5 10 15 20-25

-20

-15

-10

-5

0

5

10

15

20

25phase plan


35

%x0 : vettore condizioni iniziali close all;cla; [t,x]=ode45(@lorentz,tspan,x0); %colore sfondo set(gca,'color','black'); %crea un handle al plot del punto punto=plot3(x(1,1),x(1,2),x(1,3),'o',... 'EraseMode','xor','MarkerEdgeColor','y','MarkerFaceColor','g','MarkerSize',10); hold on %crea un handle al plot della linea scia=plot3(x(1,1),x(1,2),x(1,3),'b-',... 'EraseMode','none','LineWidth',2); axis([min(x(:,1)) max(x(:,1)) min(x(:,2)) max(x(:,2)) min(x(:,3)) max(x(:,3)) ]); hold on for j=2:length(t) set(gca,'color','black'); %specifica I dati da rappresentare nei due plot set(punto,'XData',x(j,1),'YData',x(j,2),'ZData',x(j,3)) set(scia,'XData',x(j-1:j,1),'YData',x(j-1:j,2),'ZData',x(j-1:j,3)) %esegue immediatamente il comando grafico drawnow %rallenta l’animazione pause(0.01) end hold off %sottofunzione che definisce il sistema function f=lorentz(t,x) f=zeros(3,1);f(1)=10*(x(2)-x(1)); f(2)=-x(1)*x(3)+28*x(1)-x(2); f(3)=x(1)*x(2)-8*x(3)/3; Una chiamata è: muovilorentz([0 20],[5 -30 -5])

Equazioni Differenziali Ordinarie(ODE) - unina.stidue.netunina.stidue.net/Calcolo Numerico...

Documents

Transcript of Equazioni Differenziali Ordinarie(ODE) - unina.stidue.netunina.stidue.net/Calcolo Numerico...