(x) = f (x, y(x), y (x),K, y x - Dipartimento di Scienze di...

Generalità

Problema: soluzione di una equazione differenziale alle derivate ordinarie di ordine n:

soggetta alle n condizioni iniziali:

( )( ) ( ) ( )( ) ( )( )),,,,( 11 xyxyxyxfxy nn −= K

( )( ) 1,,1,00 −== nkyxy kk K

Ovvero ricercare la soluzione di un sistema di nequazioni differenziali ordinarie del primo ordine, mediante opportune tecniche numeriche

Problema di Cauchy

Problema di Cauchy ad una sola equazione:

( ) ( )( )( )

=

=

00

,'

yxy

xyxfxy

condizione iniziale

Esempio: crescita (limitata) di una popolazione:

( ) ( ) ( )

( )

( )

( )

−−=

=ΨΨ+Ψ

=⇒

=

−= +

1/ln

e)()(1

)(

0

1d

d

00 yBK

tt

tBty

yyBtytCy

tty KCt

0 5 10 15 20 250

2

4

6

8

10

12

Problema di Cauchy: esempio

Esempio: crescita (limitata) di una popolazione:

1105.0 0 === yBC>> B=10;

>> C=1/2;

>> y0=1;

>> k=-log(B/y0-1);

>> t=linspace(0,25,100);

>> Psi=exp(C*t+k);

>> y=B*Psi./(1.+Psi);

>> plot(t,y)

Problema di Cauchy: esistenza ed unicitàIn generale, f(x,y(x)) può avere una espressione complicata o tale per cui la soluzione del problema di Cauchy non èesprimibile in forma chiusa, in tal caso si impiega in metodo numerico per ricavarne una soluzione approssimata.

Prima di cercarne una soluzione approssimata bisognerebbe chiedersi se la soluzione esiste e se è unica!!!

Definizione: una funzione f(x,y) si dice lipschitziana in yuniformemente rispetto ad x in D⊂ℜ2, se esiste una costante L>0 tale che:

DyxyxyyLyxfyxf ⊂∀−≤− ),(,),(),(),( 212121

Nota: condizione sufficiente affinché una funzione f(x,y) sia lipschitziana in D, è che fy(x,y) esista e sia limitata in D.

Problema di Cauchy: esistenza ed unicitàTeorema 1: sia f(x,y) definita e continua in Q:

e sia lipschitziana in y uniformemente rispetto a x in Q, allora esiste un’unica soluzione y(x)∈C1[x0,x0+β] del problema di Cauchy in I=[x0,x0+β], con β=min[a,b/M], M=maxQ|f(x,y)|.

( ) [ ] [ ]{ }bybyyaxxxyxQ +−∈+∈ℜ∈= 00002 ,,,,

Teorema 2: sia f(x,y) definita e continua in S:

e sia lipschitziana in y uniformemente rispetto a x in S, allora esiste un’unica soluzione y(x)∈C1[x0,x0+β] del problema di Cauchy in I=[x0,x0+β].

( ) [ ]{ }ℜ∈+∈ℜ∈= yxxxyxS ,,, 002 β

Problema di Cauchy: condizionamentoPer poter risolvere numericamente il problema di Cauchy ènecessario che il problema sia ben-posto.Definizione: il problema di Cauchy

è ben-posto se, dette y(x;y0) e y(x;y0+δ) le soluzioni con con-dizioni iniziali y(x0)=y0 e y(x0)=y0+δ rispettivamente, si ha:

( ) ( )( )xyxfxy ,' = condizione iniziale

( ) ( ) ],[;; 00 baxyxyyxy ∈∀<+− εδ

dove ε>0 è una prefissata tolleranza, purché δ=δ(ε) sia sufficientemente piccolo. (Esercizio consigliato [GL] 6.20)Nota: se f(x,y(x)) soddisfa le condizioni del Teorema 2allora il problema di Cauchy associato è ben-posto.

( ) 00 yxy =con

Metodo di Eulero esplicito

Problema di Cauchy:( ) ( )( )( )

=

+∈=

00

00 ],[,'

yxy

xxxxyxfxy β

Condizione iniziale

[ ][ ]fhxyxRxyxhfxy

fhxyxRxhyxyxy;);(,))(,()(

;);(,)(')()(

00000

00001

++==++=

),( 0001 yxhfyy +=

Discretizzazione di I: nhniihxxi /,,00 β==+= L

Sviluppo in serie di Taylor in x=x0:

Soluzione approssimata per y(xi):

Trascuro il resto R[x1,y(x1);h;f].

Sviluppo in serie di Taylor in x=x1:


[ ][ ]fhxyxRxyxhfxy

fhxyxRxhyxyxy;);(,))(,()(

;);(,)(')()(

11111

11112

++==++=

),( 1112 yxhfyy +=Soluzione approssimata per y(x2):

Il valore esatto di y(x1) non è disponibile, si usa il valore approssimato y1, ed inoltre si trascura il resto R[x1,y(x1);h;f].

Noto il generico valore approssimato xi.

Sviluppo in serie di Taylor in x=xi:


[ ][ ]fhxyxRxyxhfxy

fhxyxRxhyxyhxyxy

iiiii

iiiiii

;);(,))(,()(;);(,)(')()()( 1

++==++=+=+

),(1 iiii yxhfyy +=+Soluzione approssimata per y(xi):

Uso del valore approssimato yi≈y(xi).Si trascura il resto R[xi,y(xi);h;f].

Algoritmo: il metodo di Eulero esplicito fornisce la succes-sione di approssimazioni:


1,,0),(1 −=+=+ niyxhfyy iiii K

Alla soluzione è associato l’errore globale di troncamento:

dovuto alla linearizzazione (sostituzione delle y(x) con la retta tangente), ossia all’errore di troncamento locale e l’errore dovuto all’accumulo degli errori di troncamento locali (ci si muove lungo la retta tangente in yi con coefficiente angolare f(xi,yi), che è una approssimazione della retta tangente in y(xi)con coefficiente angolare f(xi,y(xi))).

iii yxy −= )(ε

Pb. di Cauchy:

Metodo di Eulero esplicito: esempio

( ) ( )( )

=

∈=

10

]1,0[2'

y

xxyxy

Condizione iniziale ( ) xexy 2=

Soluzione

Metodo di Eulero:

=

−=+=+=+

1

1,,02),(

0

1

y

niyhyyxhfyy iiiiii K

Discretizzazioni:

2005.0,,0)2

101.0,,0)1

2222)2(

1111)1(

=⇔===

=⇔===

nhniihxnhniihx

i

i

K

K

Metodo di Eulero esplicito: esempio

0 0.2 0.4 0.6 0.8 11

2

3

4

5

6

7

8

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

1.4

Soluzione Errore

( ) xexy 2=

1.01 =h

05.02 =h

1.01 =h

05.02 =h

• Esercizio consigliato [GL] 6.1

Problema di Cauchy:

Metodo numerici: generalità

( ) ( )( )( )

=

+∈=

00

00 ],[,'

yxy

xxxxyxfxy β

Condizione iniziale

Discretizzazione di I (non necessariamente uniforme):

nhniihxxi /,,00 β==+= L

La forma discreta del problema differenziale si può ottenere tramite una integrazione di entrambe i membri tra due nodi qualsiasi xi, xi+k:

( )∫∫++

=ki

i

ki

i

x

x

x

x

xxyxfxxy d)(,d)('

Utilizzando una formula di quadratura basata su N≤k nodi compresi tra xi e xi+k:


( ) ( )fhxyxRxyxfaxyxy ii

N

rlkilkiriki rr

;);(,)(,)()(1

+=− ∑=

−+−++

( ) ( )fhxyxRfhxyxyxyxGxy iiikikiiki ;);(,;);(,),(),(,)( 1 += −+++ K

Ovvero, in forma sintetica:

( )∫∫++

=ki

i

ki

i

x

x

x

x

xxyxfxxy d)(,d)('

Trascurando l’errore di troncamento ed utilizzando i valori approssimati calcolati:


( ) knifhyyyxGy ikikiiki −== −+++ ,,0;;,,,, 1 KK

Od analogamente:

( ) knifhyyyxhycy ikikii

k

jjijki −=+= −++

−

=++ ∑ ,,0;;,,,, 1

1

0KKφ

La funzione φ caratterizza il metodo ed è detta funzione di incremento; essendo una combinazione lineare di valori assunti dalla f(x,y(x)) ne conserva tutte le proprietà, in particolare è lipschitziana se lo è la f(x,y(x)).

Problema di Cauchy:


( ) ( )( )( )

=

+∈=

00

00 ],[,'

yxy

xxxxyxfxy β

Condizione iniziale

Discretizzazione di I (non necessariamente uniforme):

nhniihxxi /,,00 β==+= L

La forma discreta del problema differenziale si può ottenere tramite uno sviluppo in serie di Taylor della funzione y(x):

( )( )

( )!1)(

!)(,)(

!1)(

!)()()(

)(1

1

)1(

)1(1

1

)(

+++=

=+

++=+

+

=

−

++

=

∑

∑

mfh

kxyxfhxy

myh

kxyhxyhxy

mm

m

k

kk

mm

m

k

kk

ξ

ξ

Trascurando l’errore di troncamento ed utilizzando i valori approssimati calcolati, per x=xi si ha:

Metodo numerici: generalità( )

( )!1)(

!)(,)()(

)(1

1

)1(

+++=+ +

=

−

∑ mfh

kxyxfhxyhxy

mm

m

k

kk ξ

∑=

−−

+ +=m

k

iik

kii k

yxfhhyy1

)1(1

1 !),(

La derivata di ordine k della f(x,y(x)) è valutata ricorsivamente:

=

+= −−

))(,())(,(

))(,())(,())(,())(,()0(

)1()1()(

xyxfxyxf

xyxfxyxfxyxfxyxf ky

kx

k

Sostituendo e riordinando si ottiene:

( ) knifhyyyxhycy ikikii

k

jjijki −=+= −++

−

=++ ∑ ,,0;;,,,, 1

1

0KKφ

Un metodo numerico è detto esplicito o implicito a seconda se yi+k non compaia o compaia nella funzione G.Un metodo è detto one-step o ad un passo se nella funzione G compare solo la soluzione in yi, cioè se k=1.Per k>1 si ha un metodo multi-step o a più passi, i valori yi+jper j=0,…,k-1, richiesti per la valutazione di yi+k sono detti valori di innesco.

In particolare, la relazione:


( ) knifhyyyxGy ikikiiki −== −+++ ,,0;;,,,, 1 KK

( ) ( ) 1,,0;;,;;,1 −=+==+ nifhyxhyfhyxGy iiiiii Kφ

rappresenta un generico metodo one-step esplicito.

Un metodo numerico deriva dalla discretizzazione di una relazione esatta del tipo:

Metodo numerici: consistenza, ordine

( )( )fhxyxR

fhxyhkxykhxyxGkhxy;);(,

;);(,),)1((),(,)(+

+−++=+ K

Il termine R(x,y(x);h;f) è detto errore di troncamento locale, e per x=xi diventa:

( ) ( )fhxyxyxyxGxyfhxyxR ikikiikiii ;);(,),(),(,)(;);(, 1 K−+++ −=

E rappresenta la differenza tra il valore esatto della soluzioney(xi+k) ed il valore approssimato fornito dallo schema numeri-co descritto dalla funzione G operante con i valori esatti y(xi).

Un metodo numerico è detto consistente se:

Metodo numerici: consistenza, ordine( ) ( )fhxyxyxyxGxyfhxyxR ikikiikiii ;);(,),(),(,)(;);(, 1 K−+++ −=

( ) 0;);(,lim0

=→ h

fhxyxR iih

Per un metodo essere consistente significa che la forma discreta tende alla formulazione continua se il passo della discretizzazione tende a zero (e per la soluzione????).Il massimo intero positivo p per il quale:

è detto ordine di accuratezza del metodo.

( ) )(;);(, 1+= pii hOfhxyxR

Il valore nel nodo xi fornito dal metodo numerico è valutato impiegando valori approssimati già calcolati:

Metodi numerici: errore globale di tron.

La differenza tra il valore esatto ed il valore calcolato:

( )fhyyyxGy ikikiiki ;;,,,, 1 K−+++ =

iii yxy −= )(ε

è detta errore globale di troncamento. Si può scrivere:

( ) ( )[ ]fhyxyxGfhxyxyxGfhxyxRyxy

ikiiikii

iikikiki

;;,,;);(,),(,););(,()(

KK ++

+++

−++=−=ε

L’errore globale di troncamento è somma dell’errore locale di troncamento e di un errore che dipende dall’uso dei valori ap-prossimati yi anziché dei valori esatti nella valutazione di G.

Un metodo numerico è detto convergente se:

Metodi numerici: convergenza

con0maxlim00

=≤≤→ inihε

ossia se la soluzione numerica tende, al diminuire del passo di discretizzazione, alla soluzione analitica del problema di Cauchy.

β=nh

Ricordando che:Metodi numerici: stabilità

( ) ( )[ ]fhyxyxGfhxyxyxGfhxyxRyxy

ikiiikii

iikikiki

;;,,;);(,),(,););(,()(

KK ++

+++

−++=−=ε

si vede che la consistenza (errore di troncamento locale tendente a zero) non assicura la convergenza della soluzione (errore globale di troncamento tendente a zero), ma è anche indispensabile che la propagazione degli errori dovuti all’uso dei valori approssimati (accumulo) non vengano amplificati, cioè che il metodo sia stabile.Si dimostra che (teorema di Lax):

stabilità + consistenza = convergenza

Metodi numerici espliciti one-step

( ) ( )fhxyxRfhxyxhxyhxy ;);(,;);(,)()( ++=+ φ

è ottenuto dalla discretizzazione della relazione esatta:

( )fhyxhyy iiii ;;,1 φ+=+Un metodo esplicito one-step:

e per l’errore locale di troncamento si ha:

( ) ( )fhxyxh

xyhxyh

fhxyxR ;);(,)()(;);(, φ−−+=

Quindi un metodo esplicito one-step è consistente se e solo se:

( ) ))(,()(';);(,lim0

xyxfxyfhxyxh

==→φ

Metodi espliciti one-step: Eulero

],[)(''2

))(,()()(2

hxxyhxyxhfxyhxy +∈++=+ ξξ

Il metodo esplicito di Eulero può essere ottenuto dalla formula di Taylor arrestata al primo ordine:

=

−=+=+

)(

1,,0),(

00

1

xyy

niyxhfyy iiii KAlgoritmo

))(,(););(,( xyxffhxyx =φ

Condizione iniziale

)(''2

););(,(2

ξyhfhxyxR =

Il metodo è consistente

Il metodo è al primo ordine

Si può dimostrare che il metodo è convergente

Metodi espliciti one-step: Runge-Kutta

Metodi numerici con ordine superiore al primo possono ottenersi considerando più termini nell’espansione in serie di Taylor; il metodo risulta laborioso e richiede una certa regolarità per la funzione f(x,y(x)).

Metodi Runge-Kutta espliciti: vengono ricavati esprimendo la funzione incremento φ come una combinazione lineare di opportuni valori assunti dalla funzione f(x,y(x)) nell’intervallo [x,x+h], sotto la condizione che il metodo sia consistente e che abbia l’ordine di accuratezza p desiderato.

( )( )!1

)(!

)(,)()()(

1

1

)1(

+++=+ +

=

−

∑ mfh

kxyxfhxyhxy

mm

m

k

kk ξ

Metodi Runge-Kutta: HeunMetodo di Heun: è un metodo Runge-Kutta esplicito accurato al secondo ordine:

[ ]

=

−=++++=+

)(

1,,0)),(,(),(2

00

1

xyy

niyxhfyhxfyxfhyy iiiiiiii K

Si dimostra che il metodo di Heun è accurato al secondo ordine, consistente e convergente. Forma alternativa:

++=

=

++=+

),(

),(

)(2

12

1

211

hkyhxfk

yxfk

kkhyy

ii

ii

iiIl metodo di Heun è un metodo esplicito ad un passo, a due stadi.


Metodo Runge-Kutta classicoMetodo di Runge-Kutta classico: metodo accurato al quarto ordine, esplicito ad un passo a quattro stadi:

( )

++=

++=

++=

=

++++=+

34

23

12

1

43211

,2

,2

2,

2

),(

)22(6

hkyhxfk

hkyhxfk

hkyhxfk

yxfk

kkkkhyy

ii

ii

ii

ii

ii

Si dimostra che il metodo di Runge-Kutta classico èaccurato al quarto ordine, consistente e convergente.


Pb. di Cauchy:

Metodi one-step espliciti: esempio( ) ( )( )

=

∈−=

20

]2,0['

y

xxxyxy

Condizione iniziale ( ) 1++= xexy x

Soluzione

Metodo di Eulero

),(1 iiii yxhfyy +=+

Metodo di Heun

++=

=

++=+

),(

),(

)(2

12

1

211

hkyhxfk

yxfk

kkhyy

ii

ii

ii

( )

++=

++=

++=

=

++++=+

34

23

12

1

43211

,2

,2

2,

2

),(

)22(6

hkyhxfk

hkyhxfk

hkyhxfk

yxfk

kkkkhyy

ii

ii

ii

ii

ii

Runge-Kutta classico

>> t0=0;t=2;y0=2;n=10;

>> f='y-t';

>> text=linspace(t0,t,n+1);

>> ext=exp(text)+text+1;

>> [Teu,Ueu]=EULESP(t0,t,n,f,y0);

>> [Trk2,Urk2]=RK2(t0,t,n,f,y0);

>> [Trk4,Urk4]=RK4(t0,t,n,f,y0);

>> plot(Teu,Ueu,Trk2,Urk2,Trk4,Urk4,text,ext)

>> plot(Teu,ext-Ueu,Trk2,ext-Urk2,Trk4,ext-Urk4)

Metodi one-step espliciti: esempio

0 0.5 1 1.5 20

0.2

0.4

0.6

0.8

1

1.2

1.4

Metodi one-step espliciti: esempio

0 0.5 1 1.5 20

2

4

6

8

10

12

Soluzione Errore

Eulero

Heun

R.K. ClassicoEulero

Heun

R.K. Classico

Soluzione esatta


LehC

xxLp

pi

i 1)( 0 −≤

−

ε

Per un metodo esplicito one-step, in assenza di errori di arrotondamento si dimostra che:

L: costante di Lipschitz di f(x,y(x))

Cp: costante dipendente dal metodo

Teorema: sia φ(x,y(x);h)∈C0(S) lipschitziana in y, allora un metodo one-step esplicito è convergente se e solo se èconsistente. Inoltre se il metodo è di ordine p, si ha:

piii khyxy ≤−= )(ε

dove k è una costante indipendente da i e da h.


+

−≤

−

hhC

Le p

p

xxL

i

i ηε 1)( 0

Detto ηi+1 l’errore di arrotondamento che si produce nel calcolo di yi+1 ad ogni passo, si può scrivere:

se |ηj|≤η,∀j, per l’errore globale di troncamento si ha:

11 );;,( ++ ++= iiiii fhyxhyy ηφ

In particolare per il metodo di Eulero:

+

−≤

−

hhM

Le xxL

i

i ηε2

1)( 0

[ ])(''max

00 ,xyM

xxx β+∈=con:

• Esercizi consigliati [GL] 6.2, 6.4


Metodo di Eulero

hhott

Errore globale

Errore di arrotondamento

Errore di troncamento

• Per h<hott=√(2η/M) l’errore di troncamento è uguale a quello di arrotondamento e la maggiorazione dell’errore globale ha un minimo

• Per h<hott predomina l’errore di arrotondamento

• Per h>hott predomina l’errore di troncamento

Formula del punto centrale

L’algoritmo del metodo mid-point o del punto centrale è:

0),(2 112 ≥+= +++ iyxhfyy iiii

con y1, valore di innesco. Per l’errore di troncamento locale, nell’ipotesi che y∈C3(I), si ha:

],[)('''3 2

3

+∈= iii xxyhR ξξ

Quindi si tratta di un metodo esplicito a due passi del secondo ordine, che richiede, ad ogni passo, la una sola valutazione della funzione f, a differenza del metodo di Heun che necessita di due valutazioni funzionali.

(x) = f (x, y(x), y (x),K, y x - Dipartimento di Scienze di...

Documents

Transcript of (x) = f (x, y(x), y (x),K, y x - Dipartimento di Scienze di...