Soluzione Numerica di Equazioni Differenziali …riccardo.broglia/Didattica/AN2003...Metodo di...
36
Soluzione Numerica di Equazioni Differenziali Ordinarie
Transcript of Soluzione Numerica di Equazioni Differenziali …riccardo.broglia/Didattica/AN2003...Metodo di...
Soluzione Numerica di Equazioni Differenziali
Ordinarie
2
Generalità
Problema: soluzione di una equazione differenziale alle derivate ordinarie di ordine n:
( )( ) ( ) ( )( ) ( )( )),,,,( 11 xyxyxyxfxy nn −= K
soggetta alle n condizioni iniziali:( )( ) 1,,1,00 −== nkyxy kk K
Ovvero ricercare la soluzione di un sistema di nequazioni differenziali ordinarie del primo ordine, mediante opportune tecniche numeriche
3
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
4
Problema di Cauchy: esempio
Esempio: crescita (limitata) di una popolazione:
1105.0 0 === yBC
0 5 10 15 20 250
2
4
6
8
10
12
>> 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)
5
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.
6
Problema di Cauchy: esistenza ed unicitàTeorema 1: sia f(x,y) definita e continua in Q:
( ) [ ] [ ]{ }bybyyaxxxyxQ +−∈+∈ℜ∈= 00002 ,,,,
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)|.
Teorema 2: sia f(x,y) definita e continua in S:
( ) [ ]{ }ℜ∈+∈ℜ∈= yxxxyxS ,,, 002 β
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+β].
7
Problema di Cauchy: condizionamentoPer poter risolvere numericamente il problema di Cauchy ènecessario che il problema sia ben-posto.Definizione: il problema di Cauchy
( ) ( )( )xyxfxy ,' = condizione iniziale( ) 00 yxy =con
è 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:
( ) ( ) ],[;; 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.
8
Metodo di Eulero esplicito
( ) ( )( )( )
=
+∈=
00
00 ],[,'
yxy
xxxxyxfxy β
Condizione inizialeProblema di Cauchy:
Discretizzazione di I: nhniihxxi /,,00 β==+= L
Sviluppo in serie di Taylor in x=x0:[ ]
[ ]fhxyxRxyxhfxyfhxyxRxhyxyxy
;);(,))(,()(;);(,)(')()(
00000
00001
++==++=
Trascuro il resto R[x1,y(x1);h;f].
),( 0001 yxhfyy +=Soluzione approssimata per y(xi):
9
Metodo di Eulero esplicito
Sviluppo in serie di Taylor in x=x1:
[ ][ ]fhxyxRxyxhfxy
fhxyxRxhyxyxy;);(,))(,()(
;);(,)(')()(
11111
11112
++==++=
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].
),( 1112 yxhfyy +=Soluzione approssimata per y(x2):
10
Metodo di Eulero esplicito
Noto il generico valore approssimato xi.
Sviluppo in serie di Taylor in x=xi:
[ ][ ]fhxyxRxyxhfxy
fhxyxRxhyxyhxyxy
iiiii
iiiiii
;);(,))(,()(;);(,)(')()()( 1
++==++=+=+
Uso del valore approssimato yi≈y(xi).Si trascura il resto R[xi,y(xi);h;f].
),(1 iiii yxhfyy +=+Soluzione approssimata per y(xi):
11
Metodo di Eulero esplicito
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:
iii yxy −= )(ε
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))).
12
Metodo di Eulero esplicito: esempio
( ) ( )( )
=
∈=
10
]1,0[2'
y
xxyxy
Condizione iniziale ( ) xexy 2=
SoluzionePb. di Cauchy:
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
13
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
14
Metodo numerici: generalità
( ) ( )( )( )
=
+∈=
00
00 ],[,'
yxy
xxxxyxfxy β
Condizione inizialeProblema di Cauchy:
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)('
15
Metodo numerici: generalità
( )∫∫++
=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
+=− ∑=
−+−++
Ovvero, in forma sintetica:
( ) ( )fhxyxRfhxyxyxyxGxy iiikikiiki ;);(,;);(,),(),(,)( 1 += −+++ K
16
Metodo numerici: generalità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)).
17
Metodo numerici: generalità
( ) ( )( )( )
=
+∈=
00
00 ],[,'
yxy
xxxxyxfxy β
Condizione inizialeProblema di Cauchy:
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
ξ
ξ
18
Metodo numerici: generalità( )
( )!1)(
!)(,)()(
)(1
1
)1(
+++=+ +
=
−
∑ mfh
kxyxfhxyhxy
mm
m
k
kk ξ
Trascurando l�errore di troncamento ed utilizzando i valori approssimati calcolati, per x=xi si ha:
∑=
−−
+ +=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φ
19
Metodo numerici: generalità
( ) knifhyyyxGy ikikiiki −== −+++ ,,0;;,,,, 1 KK
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:
( ) ( ) 1,,0;;,;;,1 −=+==+ nifhyxhyfhyxGy iiiiii Kφ
rappresenta un generico metodo one-step esplicito.
20
Metodo numerici: consistenza, ordineUn metodo numerico deriva dalla discretizzazione di una relazione esatta del tipo:
( )( )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).
21
Metodo numerici: consistenza, ordine( ) ( )fhxyxyxyxGxyfhxyxR ikikiikiii ;);(,),(),(,)(;);(, 1 K−+++ −=
Un metodo numerico è detto consistente se:
( ) 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:
( ) )(;);(, 1+= pii hOfhxyxR
è detto ordine di accuratezza del metodo.
22
Metodi numerici: errore globale di tron.Il valore nel nodo xi fornito dal metodo numerico è valutato impiegando valori approssimati già calcolati:
( )fhyyyxGy ikikiiki ;;,,,, 1 K−+++ =
La differenza tra il valore esatto ed il valore calcolato:
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.
23
Metodi numerici: convergenza
Un metodo numerico è detto convergente se:
0maxlim00
=≤≤→ inihε β=nhcon
ossia se la soluzione numerica tende, al diminuire del passo di discretizzazione, alla soluzione analitica del problema di Cauchy.
24
Metodi numerici: stabilitàRicordando che:
( ) ( )[ ]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
25
Metodi numerici espliciti one-step
( )fhyxhyy iiii ;;,1 φ+=+Un metodo esplicito one-step:
è ottenuto dalla discretizzazione della relazione esatta:
( ) ( )fhxyxRfhxyxhxyhxy ;);(,;);(,)()( ++=+ φ
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
==→φ
26
Metodi espliciti one-step: EuleroIl metodo esplicito di Eulero può essere ottenuto dalla formula di Taylor arrestata al primo ordine:
],[)(''2
))(,()()(2
hxxyhxyxhfxyhxy +∈++=+ ξξ
=
−=+=+
)(
1,,0),(
00
1
xyy
niyxhfyy iiii KAlgoritmo
Condizione iniziale
))(,(););(,( xyxffhxyx =φ Il metodo è consistente
)(''2
););(,(2
ξyhfhxyxR = Il metodo è al primo ordine
Si può dimostrare che il metodo è convergente
27
Metodi espliciti one-step: Runge-Kutta( )
( )!1)(
!)(,)()(
)(1
1
)1(
+++=+ +
=
−
∑ mfh
kxyxfhxyhxy
mm
m
k
kk ξ
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.
28
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
ii
• Esercizio consigliato [GL] 6.6
Il metodo di Heun è un metodo esplicito ad un passo, a due stadi.
29
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
• Esercizio consigliato [GL] 6.7
Si dimostra che il metodo di Runge-Kutta classico èaccurato al quarto ordine, consistente e convergente.
30
Metodi one-step espliciti: esempio( ) ( )( )
=
∈−=
20
]2,0['
y
xxxyxy
Condizione iniziale ( ) 1++= xexy x
SoluzionePb. di Cauchy:
( )
++=
++=
++=
=
++++=+
34
23
12
1
43211
,2
,2
2,
2
),(
)22(6
hkyhxfk
hkyhxfk
hkyhxfk
yxfk
kkkkhyy
ii
ii
ii
ii
ii
Runge-Kutta classicoMetodo di Eulero
),(1 iiii yxhfyy +=+
Metodo di Heun
++=
=
++=+
),(
),(
)(2
12
1
211
hkyhxfk
yxfk
kkhyy
ii
ii
ii
31
Metodi one-step espliciti: esempio
>> 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)
32
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
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
33
Metodi numerici espliciti one-step
Per un metodo esplicito one-step, in assenza di errori di arrotondamento si dimostra che:
LehC
xxLp
pi
i 1)( 0 −≤
−
ε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.
34
Metodi numerici espliciti one-stepDetto ηi+1 l�errore di arrotondamento che si produce nel calcolo di yi+1 ad ogni passo, si può scrivere:
11 );;,( ++ ++= iiiii fhyxhyy ηφ
se |ηj|≤η,∀j, per l�errore globale di troncamento si ha:
+
−≤
−
hhC
Le p
p
xxL
i
i ηε 1)( 0
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
35
Metodi numerici espliciti one-step
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
36
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.
