Metodo di Calcolo Numerico per Equazioni differenziali Ordinarie

Approssimiamo la soluzione in:

Suddividiamo in N intervalli mediante i punti:

Una volta conosciuta la soluzione in questi punti potremo approssimare la soluzione in tutto l’intervallo (magari con una interpolazione polinomiale).

Per adesso considereremo una suddivisione uniforme:

],[ 0 Tt

Tttt N ,...,, 10

N

tThhnttn

00

00 )( ),(' ytyytfy

Considereremo un problema ai valori iniziali (di Cauchy):

La f definisce una campo di pendenze.

Metodi a un passo

In generale, studiaremo strategie che rientrano nella categoria dei cosiddetti metodi a un passo; il nome deriva dal fatto che per calcolare la soluzione numerica al tempo tn+1 é sufficiente conoscere la soluzione numerica al tempo n:

;h,f,ytΦhyy nnnn 1

Caratterizza un metodo specifico. Rappresenta una approssimazione numerica della media della funzione f tra gli istanti n, n+1.

),(1 htd;h,f,ytΦh

yynnn

nn

Errore locale di discretizzazione

(di troncamento)

Funzione Incrementale

Metodi a un passo

Consistenza: Un metodo a un passo si dice consistente nell’intervallo di integrazione se d(t,h) é infinitesimo per h tendente a zero. Piú precisamente esiste una funzione d(h) tale che:

Inoltre un metodo a un passo si dirá di ordine di consistenza p se:

Convergenza: se la Φ é lipschitziana rispetto a y si ha convergenza. Richiede anche la stabilitá….vedi propagazione dell’errore….

0)(lim)(),(0

hdhdhtdh

)()( phOhd

2

1

)2

1

10 ) ,(

hy''(ξy)y(t

,....,N- nytf hyy

nn

nnnn

1. Metodo di Eulero Esplicito

)(

)()()()()(

dt

tdyhtyhty

h

tyhty

dt

tdy

Metodi a un passo

Sviluppo di Taylor sino al primo ordine.

Metodi a un passo

)(

)()()()()(

dt

tdyhhtyty

h

htyty

dt

tdy

10 ) ,( 111 ,....,N- nytf hyy nnnn

2. Metodo di Eulero Implicito

Ha lo stesso grado di accuratezza del metodo esplicito peró richiede la risoluzione di una equazione non lineare (comunque in alcuni casi i metodi impliciti possono presentare dei vantaggi).

Metodi a un passo

fffyffy

ytfy

) O(h)y''(yh)y'(yhy)y(t

ytyt

'''

),('2

1 30

2000

Se arrestiamo lo sviluppo di Taylor al secondo ordine, e calcoliamo la derivata seconda di y, si ottiene: (metodo del secondo ordine)

1.0

2

1 21

,N-,..n

), yf(ty

),yf(t

t

),yf(th), yhf(tyy nn

nnnnnnnn

3. Metodo basato sullo Sviluppo di Taylor arresto al 2 ordine:

Metodi a un passo

Coordinate di P:

),(2

22/1

nnnp

nnp

ytfh

yy

httx

Per costruire un metodo del secondo ordine (come il precedente) senza dover calcolare e valutare le derivate di f, si puó ragionare in questo modo: htt nn 1

),(22/

),(2/ nnnp

npnn

np ytfh

yyh

yyytf

h

yym

nt

)( nty

1nt

h

2/1nt

2/h

Ppy

Soluzione reale

1.0

))(221

,N-,..n

, ytfh

, yh

f(thyy nnnnnn

4. Metodo di Eulero Modificato:

In pratica, abbiamo approssimato la soluzione reale nell’ intervallo con la retta tangente in tn per stimare il valore della funzione al tempo tn+h/2 (punto P); infine abbiamo approssimato l’incremento della soluzione in attraverso la pendenza in P (si cerca di migliorare la stima della retta secante tra due instanti di integrazione).

Questo metodo ha una accuratezza del secondo ordine.

2/, htt nn

htt nn ,

Metodi a un passo

Metodi a un passo

Come abbiamo appena visto per aumentare l’accuratezza il metodo di Eulero esplicito basta migliorare l’approssimazione della secante tra due instanti di integrazione: un’altra idea potrebbe essere l’utilizzo di una media tra le pendenze ai tempi tn e tn+h.

1.0

))(),2

11

,N-,..n

, ythfh, yf(t yf(thyy nnnnnnnn

5. Metodo di Heun (differente del Metodo dei Trapezi):

Per valutare la funzione al passo tn+h si é utilizzato una passo dell’Eulero esplicito si é considarato un’incremento lineare con pendenza definita al tempo tn. Differisce dal metodo dei Trapezi poiché quest’ultimo é implicito (vedi dopo).

Metodi a un passo: Metodi Runge-Kutta

khyhtfk

,N-,.. nkhyy

P

jjijnini

M

iiinn

1

11

),(

1.0

6. Metodo Runge-kutta a M stadi (livelli):

Metodo Esplicito

Metodo Imsplicito

Metodo Semi-imsplicito

1iP

MP

iP

Generalizzazione delle osservazioni utilizzate per arrivare all’ Eulero modificato e alla formula di Heun.

Metodi a un passo: Metodi Runge-Kutta

La denominazione di esplicito, implicito, o semi-implicito dipende dalla minore o maggiore facilitá nel derivare i vari ki : in un caso si ricaveranno in cascata, mentre nell’altro si dovrá risolver un sistema di equazioni. Possiamo giungere ai metodi Runge-kutta in altro modo; infatti in generale possiamo utlizzare varie formule di quadratura per ottenere metodi giá conosciuti:

dttytfyyn

n

t

tnn

1

))(,(1

),( nn ytfh

),( 11 nn ytfh

),(),(2 11 nnnn ytfytfh

Eulero esplicito

Eulero implicito

7. Metodo dei Trapezi

(ricavato dalla omonima formula di quadratura)

Dunque sfruttando delle formule di quadratura che utilizzano i punti:

Miht in ,...,1

Dunque sfruttando altre formule di quadratura che utilizzano i nodi :

Miht in ,...,1

Metodi a un passo: Metodi Runge-Kutta

che dipendono dai corrispondenti valori incogniti:

)( hty in

Possiamo ottenere formule del tipo:

))(( ))(,(1

1 nnin

M

iininn tyyhtyhtfhyy

Il problema viene risolto approssimando a loro volta l’integrale seguente con una fomula di quadratura:

P

jinninij

ht

tninn

htyhtfh

dttytftyhtyin

n

11

1

))(,(

))(,()()(

Metodi a un passo: Metodi Runge-Kutta

La formula di quadratura in questo caso potrá utlizzare tutti o solo alcuni nodi (diaciamo P):

Chiaramente se P=i-1 sará un metodo implicito; inoltre P al massimo sará uguale a M. Se imponiamo che la formula di quadratura sia esatta almeno per funzioni costanti per ogni i troviamo la relazione:

MiP

jiji ,...,1

1

Unendo i varii procendimenti di quadratura troviamo le formule di Runge-Kutta:

Metodi a un passo: Metodi Runge-Kutta

YhtfhyY

,N-,.. nYhtfhyy

P

jjinijni

M

iiininn

1

11

),(

1.0),(

Queste formule possono essere rappresentate in forma sintentica con una tabella di coefficienti (consideriamo il caso piú generale P=M):

11

MM

M1

1M

1 M……………

……………

……………

……

…

……

…

1

MVettore W dei Pesi

Matrice dei coefficienti B

Vettore delle ascisse

Se la matrice B é triangolare inferiore il metodo sará esplicito e gli Yi si calcolano facilmente in cascata. In questo caso la condizione sugli alfa impone: . Se include anche la diagonale sará semiesplicito e la soluzione sará ricorsiva, mentre con B piena il metodo si dice implicito e richiede la risoluzione di un sistema non lineare.

Metodi a un passo: Metodi Runge-Kutta

01

0 0

1Eulero Esplicito

1 1

1Eulero Implicito

Trapezi

0 0

1/2

01 1/2 1/2

1/2

A un livello esplicito

A un livello implicito

A due livelli semi-implicito

Metodi a un passo: Metodi Runge-Kutta

Heun

0 0

1/2

0

1 1 0

1/2

Eulero Modificato

0 0

1

0

1/2 1/2 0

0

Metodi a un passo: Metodi Runge-Kutta

Il massimo ordine di accuratezza p(M) raggiungibile con un metodo a M livelli varia in questo modo:

M p(M)

1

2

3

4

5

6

7

8

1

2

3

4

4

5

5

6

Metodi espliciti Metodi impliciti

p(M)=2M

Metodi a un passo: Metodi Runge-Kutta

khyhtfk

ytfk

kkhyy

nn

nn

nn

),(

),(

12

1

22111

Troviamo i coeff. per i metodi a due stadi espliciti imponendo un certo grado di accuratezza:

hOy

ff

t

fhfhyy

hOy

ff

t

fhf hO

y

fk

t

fhf k

nn )()(

)()(

32

2211

2212

Sviluppando k2 secondo Taylor arrestando al primo ordine:

Qui risulta al secondo…

Metodi a un passo: Metodi Runge-Kutta

hOy

ff

t

fhfhy y nn )()( 3

22

211

hOy

ff

t

fhfhy y nn )( 32

1

Confrontando la formula precedente con il reale sviluppo di Taylor della formala arrestato al secondo ordine, troviamo le condizioni sui coefficienti:

I coeff. dovranno compiere:

2

12

1

1

2

2

21

Queste condizioni sono rispettate da Eulero Modificato e da Heun. Inoltre si vede α=β come avevamo giá detto.

Analisi dell’errore: Metodi Runge-Kutta

L’errore di troncamento locale (errore locale di discretizzazione) indica l’errore all’integrare un passo tra due instanti di tempo…ora questo errore si propagherá al passo successivo sommandosi al seguente errore di integrazione.

),(1 htd;h,f,ytΦh

yynnn

nn

Errore locale di discretizzazione

E’ utile definire ora l’errore di troncamento locale al passo n-esimo:

)()()(),( 1 nnnnn t,ytΦhtytyht Valori esatti !!

In pratica l’errore che commetto integrando dal passo n a n+1, supponendo perfettamente conosciuta la soluzione al tempo n.

Analisi dell’errore: Metodi Runge-Kutta

E’ chiaro che all’errore di trocamento locale al passo n, dovremo sommare l’errore commesso dalle precedenti integrazioni:

Errore accumulato totale

),()( htEEEtyyE nPRLOCPRnnn

Soluzione esatta

Errore di propagazione

Errore di troncamento locale

Stabilitá

Consistenza

))(,()(1 nnnnPR tythtyyE

Valore ottenuto dalle varie integrazioni

Valore esatto

Analisi dell’errore: Metodi Runge-Kutta

In generale per definire un metodo convergente si richiede che l’errore di troncamento locale tenda a zero al decrescere il passo di integrazione h (consistenza) e che l’errore di propagazione non si amplifichi passo dopo passo (stabilitá) .

Convergenza = Consistenza + Stabilitá

Se h é infinitesimo lo é anche l’errore di troncamento

Gli errori non si amplificano al propagarsi

Diciamo che la consistenza é una condizione statica, che suppone la convergenza al diminuire il passo di integrazione h (ovvero che la nostra approssimazione migliori con un passo h piú piccolo).

La stabilitá controlla la dinamica del nostro modello, in modo tale che errori successivi non portino a approssimazioni assolutamente erronee.

Analisi dell’errore: Metodi Runge-Kutta

Stabilitá: Metodi Runge-Kutta

In generale e’ difficile analizzare la stabilitá di un metodo, per questo ci limitiamo a una classe particolare di equazioni differenziali test:

i

yy

yytfy

0)0(

),('

La cui soluzione generale é :

))sin()(cos()( 0 titeyty t

Noi conisidereremo ovviamente le soluzioni stabili con alfa<0.

Stabilitá: Metodi Runge-Kutta

Applichiamo ora per esempio il metodo di Eulero esplicito:

nnnnnnn yhyhyytfhyy )1(),(1

Questa equazione alle differenze é stabile se:

-2 -1

hα

hβ

i

11 h Regione di assoluta stabilitá per il metodo di Eulero esplicito: Il passo h deve essere sufficientemente piccolo (dato un λ).

Stabilitá: Metodi Runge-Kutta

Altre regioni di assolutá stabilitá:

11

1

hEulero Implicito

12/1

2/1

h

h Trapezi

Stabilitá: Metodi Runge-Kutta

Se le regioni di assoluta stabilitá contengono il semipiano α<0 allora il metodo si dice incondizionatamente stabile o assolutamente stabile poiché risulta stabile per tutti i λ della equazione test stabili, e per ogni passo h.

I Metodi impliciti risultano migliori se si analizza la stabilitá.

Pur essendo l’equazione test un caso particolare puó servire per studiare almeno localmente equazioni piú generali. Infatti intorno a un punto (tn,yn) possiamo linearizzare rispetto a y:

y

ytf

yyy

ytfytfy

),(

),(),('