Equazioni DifferenzialiOrdinarie

Metodi Multi-step

Studenti:Cutuli SebastianoCantone Salvatore

Seminario

Metodi Matematici per l’Ottimizzazione

A.A. 2011/2012

Indice

• 1 Introduzione

• 2 I Metodi Multi-Step

• 2.1 Metodi impliciti• 2.1.1 Metodo di Adams-Bashford

• 2.2 Metodi espliciti• 2.2.1 Metodo di Adams-Moulton

Introduzione

Si definisce equazione differenziale un’equazione che coinvolge una o più derivate di una funzione incognita. Se tutte le derivate sono calcolate rispetto ad una sola variabile indipendente, l’equazione si dirà equazione differenziale ordinaria (ODE).

Quando sono presenti derivate rispetto a più variabili indipendenti, avremo invece una equazione differenziale alle derivate parziali (PDE).

Una equazione differenziale avrà ordine n, se n è l’ordine massimo delle derivate che vi compaiono.

La forma generale di una ODE di ordine n è:

Dove y=y(x) è la funzione cercata.

Una funzione è detta soluzione di una ODE se essa riduce l’equazione ad una identità quando viene sostituita nell’equazione.

0),...,,,( )( nyyyxF

Introduzione

ODE di ordine 1

La forma generale di una equazione differenziale ordinaria di primo ordine è:

dove y=y(x) è la funzione incognita.Per semplicità assumiamo che l’equazione sia scritta nella forma:

con la condizione y=y0 quando x=x0:

0),,( yyxF

),( yxfy

00 )( yxy

Esempio…Determinare le primitive di una funzione f(x) significa risolvere

y’(x) = f(x)

dove l'incognita è la funzione y(x).

Questa equazione è un semplice esempio di equazione differenziale.

Ad esempio se:

y’(x) = 2x

le soluzioni che si ottengono integrando 2x sono

y(x) = x2+ c

ossia le soluzioni sono incognite e ciascuna è individuata da un diverso valore della costante reale c, ognuna detta soluzione particolare.

La costante c può essere determinata imponendo un'ulteriore condizione.

Ad esempio se vogliamo che y(1) = 3 allora c = 2 e y(x) =x2 + 2.

ODE di ordine 1L’interpretazione geometrica della soluzione è data dal fatto che ogni soluzione particolare è rappresentata da una curva nel piano (x,y) e, preso un punto P(x,y) arbitrario , possiamo calcolare f(x,y) che deve essere uguale alla pendenza della tangente alla curva desiderata nel punto P.

Campo di direzione

),(tan yxf

Problema di Cauchy

L’obiettivo è capire come evolve un sistema a partire da alcuni dati iniziali: ciò si traduce imponendo che la soluzione del problema differenziale passi per un determinato punto del piano x,y.

Si ottiene allora il cosiddetto problema ai valori iniziali o problema di Cauchy:

Sia con valore iniziale.

trovare una funzione con y derivabile in

tale che:

IIxf ,:Iy :

00)(

))(,(

yxy

Ixxyxfy

I0y

Problema di Cauchy

Alcune proprietà che devono essere verificate dal problema di Cauchy sono:

• Esistenza: vedere se il problema ammette esistenza della soluzione.

• Unicità: vedere se la soluzione è unica.

Teorema di esistenza ed unicità

Sia un dominio ed una funzione

continua che soddisfi la condizione di Lipschitz:

Comunque si scelgano e una costante segue che:

ha soluzione unica in tale intervallo.

1 nD nDf :

vuLvxfuxf ),(),(

Dvxux ),(),,(0L

00

0000 )(

),(:0],,[,),(

uxu

uxfuxxDux

Metodi Multistep

Per risolvere le equazioni differenziali si ricorre a metodi numerici che possono essere ad un passo o a più passi, poiché trovare una soluzione analitica è un’operazione piuttosto difficile.

Un metodo numerico per risolvere Equazioni Differenziali Lineari inizia da un punto iniziale e poi tramite un passo computazionale (metodo) trova il punto di soluzione successiva.

Il processo continua con i passi successivi per mappare la soluzione.

Metodi Multistep

• I Metodi a Passo Singolo (come ad esempio il Metodo di Eulero) si riferiscono ad un solo punto precedente e dei suoi derivati per determinare il valore corrente.

• I Metodi a più Passi o Multistep ricercano una soluzione più accurata possibile, mantenendo e utilizzando le informazioni provenienti da passaggi precedenti, piuttosto che abbandonarle. Di conseguenza, i metodi multistep si riferiscono a diversi punti precedenti ed ai loro valori derivati.

Metodi Multistep Impliciti ed Espliciti

• Un metodo si dice esplicito se yi+1 dipende solo dai valori dei passi precedenti.

• Un metodo si dice implicito se yi+1 dipende anche da se stesso, attraverso la funzione f.

Questi ultimi metodi richiedono la soluzione di una equazione non lineare se f è non lineare in y.

Un’equazione differenziale è detta lineare se la funzione f è lineare nella funzione incognita e nelle derivate, ovvero è della forma:

0)()(')( xcyxbyxa

Instabilità dei metodi one-step

Nei metodi One-Step di tipo esplicito, dove ad ogni passo viene calcolato la soluzione nel punto ti+1,

utilizzano in modo opportuno i valori già noti nei nodi precedenti, sono schemi in generale poco costosi e facili da implementare, ma presentano spesso dei fenomeni di instabilità numerica in particolare quando la soluzione esatta del problema che si sta approssimando presenta dei forti gradienti (ricordiamo che il gradiente di una funzione è definito come il vettore che ha per componenti cartesiane le derivate parziali della funzione. Il gradiente rappresenta quindi la direzione di massimo incremento di una funzione di n variabili).

Instabilità dei metodi one-step

Per ottenere una soluzione non oscillatoria si dovrebbe utilizzare un passo h così piccolo da rendere proibitivo l'uso di metodi espliciti. In questi casi è indispensabile ricorrere a metodi impliciti, cioè metodi in cui la funzione f deve essere valutata anche in yi+1.

Vediamo allora come possono essere ricavati dei metodi di questo tipo.

Dal teorema fondamentale del calcolo integrale e dalla definizione di equazione differenziale otteniamo:

1

1 ))(,(n

n

nn

t

t

dttytfyy

Instabilità dei metodi one-step

Se adesso usiamo degli schemi numerici per approssimare l’integrale della slide precedente otteniamo degli schemi per approssimare il problema di Cauchy. Se per esempio approssimiamo l'integrale con l'area del rettangolo di base h ed altezza f(tn; un) otteniamo il Metodo di Eulero Esplicito:

mentre se usiamo l'area del rettangolo di base h ed altezza f(tn+1; un+1) otteniamo il Metodo di Eulero Implicito:

),(1 nnnn uthfuu

),( 111 nnnn uthfuu

Instabilità dei metodi one-step

Utilizzando la stessa tecnica si possono ottenere anche metodi così detti a più passi, cioè in cui il valore di un+1 dipende non solo da un un, ma

anche dagli altri valori { un-1, un-2,.…., un-q } calcolati

precedentemente.

Sostituendo nell’integrale

l'integrale del polinomio che interpola f nei nodi

{ tn-q, tn-q+1,….,tn }, si ottengono i così detti Metodi di

Adams-Bashforth che sono di tipo esplicito, mentre sostituendo l'integrale del polinomio che interpola f nei nodi { tn-q, tn-q+1,….,tn, tn+1 } si ottengono i così

detti Metodi di Adams-Moulton che sono di tipo implicito.

1

1 ))(,(n

n

nn

t

t

dttytfyy

Schema generale

L'equazione principale dei Metodi Multistep si ottiene dall'equazione

supponendo che l'intervallo di integrazione [ t0; T ] sia suddiviso in un certo numero di intervallini (N) di ampiezza h (passo di discretizzazione ((T-t0)/N)

ed integrando l'equazione da tn-j a tn+k:

Applichiamo una formula di quadratura per il calcolo dell'integrale.

1

1 ))(,(n

n

nn

t

t

dttytfyy

kn

jn

jnkn

t

tdttytftyty ))(,()()(

Schema generale

Per il calcolo dell'integrale si applicano le formule di Quadratura di Newton-Cotes.

Utilizzando q+1 punti { tn-q,….,tn } costruiamo il polinomio di Lagrange ed integriamolo nell'intervallo [tn-j ; tn+k].

dove Li(x) rappresentano i polinomi elementari di Lagrange, cioè:

q

i

iininq xLytfxp0

)(),()(

q

ill lnin

lni

xx

xxxL

0

)(

Schema generaleSostituendo quindi il polinomio nell'integrale, si

ottiene la seguente relazione fra i valori della yi

dove, per compattezza di notazione si è posto

e dove i coefficienti βqi si ottengono per integrazione dei polinomi elementari di Lagrange

A seconda dei diversi valori di k, j e q si ottengono differenti metodi multistep.

q

i

inqijnkn fhyy0

),( lll ytff

dxli

lxt

tdttL

h

k

j

q

ill

kn

jn

iqi

0

)(1

Schema generale

In particolare, per k = 1 e j = 0 si ottengono i cosiddetti metodi di Adams-Bashforth

(di tipo esplicito)

Per k = 0 e j = 1 otteniamo i cosiddetti metodi di Adams-Moulton (di tipo implicito)

o meglio

q

i

inqinn fhyy0

1

q

i

inqinn fhyy0

1

q

iinqinn fhyy

011

Metodo di adams-bashforth

Questi metodi vengono detti espliciti poiché per calcolare il passo n-esimo si utilizza una combinazione lineare dei primi q+1 passi.

con

Otteniamo diversi metodi di tale formula se a partire da questa poniamo p=1,2,3,…

Se difatti poniamo p=1 (Adams-Bashforth ad 1 passo) ritorniamo ai metodi one-step perché l’equazione diventa:

q

i

inqinn fhyy0

1 dxli

lxq

ill

qi

1

0 0

11 nnn hfyy

Metodo di adams-bashforth

Se poniamo p=2 (Adams-Bashforth a 2 passi) :

Ovvero, per trovare il passo (n+1)-esimo vengono utilizzati i valori trovati al passo n-2 ed n-1.

Oppure ancora se poniamo p=3 (Adams-Bashforth a 3 passi) si ottiene:

]3[2

211 nnnn ffh

yy

]51623[12

3211 nnnnn fffh

yy

Ode & matlab

Risoluzione dell’equazione differenziale:

con soluzione:

- Intervallo di approssimazione [t0,tf]=[0,5]

- N° approssimazioni (passi) n = 10

- Step-Size : h=(tf-t0)/n=1/2

Poiché tratteremo un metodo esplicito a 2 passi, ed il valore iniziale non è sufficiente per eseguire lo script di risoluzione, utilizzeremo il metodo di Eulero esplicito per calcolare il valore y(2)

1)0(

)5(5'

y

tseneyy t

)5cos( te t

Ode & matlab

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

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1.2Adams-Bashforth 2 num.passi=10

AB2

f(x)

Ode & matlab

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

-2

-1

0

1

2

3Adams-Bashforth 3 num.passi=10

AB3

f(x)

Ode & matlab

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

-10

-5

0

5

10

15

20

25Adams-Bashforth 4 num.passi=10

AB4

f(x)

Ode & matlabconsiderazioni metodo a.b.

Come è possibile notare l’errore cresce soprattutto con il metodo di A.B.4, il quale non converge alla soluzione.

1 2 3 4 5 6 7 8 9 100

5

10

15

20

25Andamento Errore

ERRORE AB2

ERRORE AB4ERRORE AB3

Ode & matlabconsiderazioni metodo a.b.

Per ottenere la convergenza abbiamo usato le seguenti accorgenze:

• Passo di integrazione ¼ per tutti i metodi

Ottenendo :

• Andamento errore diminuito rispetto ai valori precedentemente rilevati

Ode & matlabconsiderazioni metodo a.b.

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

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1Adams-Bashforth 2 num.passi=20

AB2

f(x)

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

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1Adams-Bashforth 3 num.passi=20

AB3

f(x)

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

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1Adams-Bashforth 4 num.passi=20

AB4

f(x)

0 2 4 6 8 10 12 14 16 18 200

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Andamento errore con num.passi 20

ERRORE AB2

ERRORE AB4ERRORE AB3

Risultati sperimentaliSono stati trovati i seguenti risultati sperimentali (errori relativi):

NERRORE

BASHFORTH 2

ERRORE

BASHFORTH 3

ERRORE

BASHFORTH 4

10 8,936618 314,5462 3.170,614

20 0,108160 0,02121 0,50884

40 0,006005 0,010660 0,00280

80 0,000699 0,00122 0,00137

Osservazioni: Aumentando il n° di passi (dunque dimezzando man mano h), la soluzione approssimata si avvicina sempre più alla soluzione esatta, indice del fatto che il metodo numerico converge.Inoltre è facile osservare che la funzione che si avvicina più di tutte alla soluzione reale è quella di ordine 2.

METODO DI ADAMS-MOULTONSono molte le analogie tra il metodo di Adams-Bashforth ed il metodo di Adams-Moulton, infatti sostituendo all'integrale della funzione f(t; y(t)) nell’equazione:

l'integrale del polinomio interpolatore di Lagrange (con grado p+1 per ottenere il metodo Adams-Moulton a p+1 passi) che interpola la f nei nodi, si ottengono i così detti metodi di Adams-Moulton che fanno parte della famiglia dei metodi impliciti. Pertanto il polinomio interpolatore utilizza non solo i punti tn-p, tn-p+1,….., tn, analogamente ad Adams-Bashforth, metodo esplicito, ma anche tn+1.

kn

jn

jnkn

t

tdttytftyty ))(,()()(

Illustrazione metodo a.moulton

• Nella figura apprezziamo

che il metodo di Adams-

Moulton tiene conto del

passo tn+1 rispetto a

Bashforth, pertanto tale

metodo implicito

dovrebbe avere

un’accuratezza

maggiore.

Formule generali Adams-Moulton

• Dopo tutte le considerazioni del caso la formula generale del metodo Moulton si presenta così:

• Con

• Tuttavia preferiamo scrivere la formula secondo questa forma:

q

i

inqinn fhyy0

1 1

0

1

q

i

inqinn fhyy

dxli

lxq

q

ill

i

0

1 0

METODO DI ADAMS-MOULTONA partire dalla formula precedente iteriamo p sulla

sommatoria e vediamo come si presenta la

formula per p=1,2,3,4 e 5

• p=1 yn= yn-1 + hfn (ovvero Eulero - implicito)

• p=2 yn+1= yn + h/2 [fn+1 + fn ]

• p=3 yn+2= yn+1 + h (5/12 fn+2 – 2/3 fn+1 -1/12 fn ]

• p=4 yn+3= yn+2 + h (3/8 fn+3 + 19/24 fn+2 -5/24 fn+1 +1/24 fn]

• p=5 yn+4= yn+3 + h/720 (251 fn+4 + 646 fn+3 - 264 fn+2 + 106 fn+1 -19 fn ]

• …..

METODO DI ADAMS-MOULTONPoiché si tratta di un metodo implicito dove ad ogni

passo

è richiesta la soluzione di un equazione non lineare, è

necessario trovare un metodo per risolvere le equazioni

non lineari: ad esempio nel metodo di Eulero Implicito la

soluzione è data dalla formula:

yn+1 = yn + h*f(tn+1,yn+1)

dunque yn+1 è la soluzione dell’equazione ad ogni passo:

yn+1 - yn - h*f(tn+1,yn+1) = 0

Se f ha una forma semplice si può risolvere facilmente yn+1

altrimenti è necessario ricorrere ad un metodo per

la risoluzione di equazioni non lineari.

Funzione fzero di matlabMatlab calcola la soluzione del problema f(x)=0 con la function fzero.In fzero sono accoppiati un metodo a convergenza veloce (interpolazione quadratica inversa) insieme ad un metodo di tipo bisezione, per determinare un intervallo in cui f cambi segno, e per evitare che il metodo diverga.L’espressione è x = fzero(funzione; x0) dove “funzione” è una stringa di caratteri che contiene il nome della funzione di cui si sta cercando uno zero, ed x0 è un valore di partenza dal quale cominciare la ricerca.Per risolvere l’equazione yn+1 - yn - h*f(tn+1,yn+1) = 0è stato utilizzato il comando

y(n+1) = fzero |@(y_inc) y_inc - y(n) -h*f (t(n+1),y_inc), y(n)|

dove ‘@(y_inc) y_inc - h*f (t(n+1) – y(n)’ indica che si vuole

trovare la soluzione dell’equazione y_inc -h*f (t(n+1),y_inc)-y(n)=0

con incognita y_inc; il parametro y(n) indica che verrà scelto questo valore(il valore precedente) come stima iniziale della soluzione.

Ode & matlab

Risoluzione dell’equazione differenziale:

con soluzione:

- Intervallo di approssimazione [t0,tf]=[0,5]

- N° approssimazioni (passi) n = 10

- Step-Size : h=(tf-t0)/n=1/2

Poiché Adams-Moulton di 2° ordine è un metodo a 2

passi, utilizzeremo il metodo di Eulero implicito per

calcolare il valore y(2) e analogamente y(3), y(4) e

y(5) per i metodi di ordine 3°,4° e 5°.

1)0(

)5(5'

y

tseneyy t

)5cos( te t

Risultati sperimentali

Risultati sperimentali

Risultati sperimentaliSono stati trovati i seguenti risultati sperimentali (errori relativi):

N ERRORE

MOULTON 2

ERRORE MOULTON 3

ERRORE MOULTON 4

ERRORE MOULTON 5

10 0,160915 1,192922 2,143834 1,470348

20 0,044363 0,026204 0,113504 0,114283

40 0,003998 0,0055349 0,026335 0,0054211

80 0,000555 0,001929 0,00007787 0,007616

160 0,00066595 0,01657 0,0422435 0,02716

Osservazioni: Aumentando il n° di passi (dunque dimezzando man mano h), la soluzione approssimata si avvicina sempre più alla soluzione esatta, indice del fatto che il metodo numerico converge.Inoltre è facile osservare che la funzione che si avvicina più di tutte alla soluzione reale è quella di ordine 2.