ODE
PROBLEMA DI CAUCHY IN 1-D
Sia f : I x RR, I  R.
Cerchiamo y: I R con y derivabile in I:
 y'  f ( x, y( x))

 y( x0 )  y0
 x I
(1)
Se x  t tale problema è detto PROBLEMA AI VALORI
INIZIALI.
Supponiamo che siano verificate le condizioni di esistenza
e unicità della soluzione.
Teorema di esistenza e unicità
Sia G  n+1 un dominio e f : G  n una funzione continua
che soddisfi la condizione di Lipschitz:
f ( x, u)  f ( x, v)  L u  v
 (x,u),(x,v)  G e qualche costante L > 0.
  (x0,u0)  G  [x0-a, x0+a] con a > 0 tale che il problema
u '  f ( x, u )

 u ( x0 )  u0
abbia soluzione unica in tale intervallo.
SUCCESSIONE DI APPROSSIMAZIONI DELLA y
Esistono metodi numerici ad un passo o a più passi.
DEF: Un metodo numerico si dice ad un passo se  n  0,
yn+1 dipende solo da yn.
METODI AD UN PASSO
Sviluppo in serie di Taylor di y(x) attorno xi.
Supponendo che y(x) sia sufficientemente regolare
h2
y ( xi 1 )  y ( xi )  hy ' ( xi ) 
y ' ' ( xi )  ...
2
tronchiamo al k-esimo termine :
yi 1  yi  hTk ( xi , yi ; h)
con
h
h k 1 ( k )
Tk ( xi , yi ; h)  y ' ( xi )  y ' ' ( xi )  ... 
y ( xi )
2
k!
Poiché è richiesto il calcolo delle derivate non è conveniente.
E’ meglio usare metodi ad un passo che utilizzano
l’informazione al passo precedente per calcolare la soluzione
al passo successivo.
yi 1  yi  h ( xi , yi ; h)
METODO DI EULERO
Ponendo k = 1 si ottiene :
yi 1  yi  hf ( xi , yi )
che è il metodo di Eulero in avanti o esplicito. Oppure
yi 1  yi  hf ( xi 1 , yi 1 )
che è il metodo di Eulero all’indietro o implicito.
DEF:
Un metodo si dice esplicito se yi+1 dipende solo dai valori ai
passi precedenti.
Un metodo si dice implicito se yi+1 dipende da se stessa
attraverso f.
Questi ultimi richiedono la risoluzione di un problema non
lineare se f non è lineare in y.
METODO DEI TRAPEZI
h
yi 1  yi  [ f i  f i 1 ]
2
dove si è posto f i  f ( xi , yi ) . La suddetta formula è stata
ricavata integrando la (1) ed applicando il metodo del trapezio.
METODO DI HEUN
h
yi 1  yi  [ f i  f ( xi 1 , yi  hf ( xi , yi ))]
2
La formula appena esposta è stata ricavata applicando il metodo
del trapezio ed utilizzando il metodo di Eulero in avanti per
calcolare le yi+1.
ANALISI DEI METODI AD UN PASSO
Indicando, come prima, con y(xi+1) la soluzione in xi+1 e con yi+1
la soluzione approssimata in xi+1 si ha :
yi 1  yi  h ( xi , yi ; h)
y ( xi 1 )  y ( xi )  h ( xi , y ( xi ); h)   i 1
dove  è l’errore al passo i+1. Riscriviamolo come
 i 1  h i 1 (h)
La quantità τi+1(h) è detta errore di troncamento locale.
Definiamo, invece, errore di troncamento globale la quantità
 (h)  max  i 1 (h)
0i  N 1
τ dipende dalla soluzione del problema di Cauchy.
La funzione incremento  caratterizza completamente il metodo
ad un passo ed è tale che
lim  ( yi , y ( xi ); h)  y ( xi )  f ( xi , y ( xi ))
h 0
y ( xi 1 )  y ( xi )
 y ' ( xi ), si ha :
Pertanto, poiché lim
h 0
h
lim  i ( h)  0 da cui lim  (h)  0.
h 0
h 0
che dà la consistenza del metodo numerico con il problema di
Cauchy.
Un metodo di dice consistente quando lim  ( h)  0.
h 0
p
Un metodo si dice consistente di ordine p se  (h)   (h ) , h0.
ES.:Dimostrare la consistenza dei metodi di Eulero ed Heun.
ZERO - STABILITA’
Un metodo numerico del tipo
yi 1  yi  h ( xi , yi ; h)
si dice zero-stabile se
h0   C  0 : Zi  yi  C
h  [0, h0 ]
dove Zi e yi sono le soluzioni di :
Z i 1  Z i  h[ ( xi , Z i ; h)   i 1 ]

Z 0  y0   0

con  i  
 yi 1  yi  h ( xi , yi ; h)

y ( x0 )  y0

Tale stabilità riguarda il comportamento del metodo numerico
quando h0. Essa assicura che il metodo sia poco sensibile
alle piccole perturbazioni.
DEF.: Un metodo si dice convergente se i
y( xi )  yi  C (h)
Un metodo si dice convergente di ordine p se i
y( xi )  yi  Ch p
TEOREMA DI CONVERGENZA
 y '  f ( x, y ( x))
Ip:-y(x) sia soluzione di 
 y( x0 )  y0
-yi+1 sia soluzione approssimata, yi 1  yi  h ( xi , yi ; h)
- sia Lipschitziana nella 2a variabile
-  ( x, u; h)   ( x, v; h)  L u  v x  [a, b]
0  h  b  a,  u, v  , h 
-sia   max  i , y( xi 1 )  y( xi )  h ( xi , y( xi ); h)  h i
ba
n
l j  y( x j )  y j 
Ts:
Dim.

L
[l
L ( x j  x0 )
 1]
j  1,..., N
li 1  y ( xi 1 )  yi 1 
 y ( xi )  h ( xi , y ( xi ); h)  h i  yi  h ( xi , yi ; h) 
 li  h[ ( xi , y ( xi ); h)   ( xi , yi ; h)]  h i
li 1  (1  hL) li  h  l1  (1  hL) l0  h
l2  (1  hL) l1  h  (1  hL) 2 l0  h [1  (1  hL)]
...
n
(
1

hL
)
1
n
ln  (1  hL) l0  h
hL
n
n
(
1


)

l
ma
con   0, l0  0 . Da cui la tesi.
Poiché la condizione di consistenza è 0 per h0, si ha:
TEOREMA: Un metodo ad un passo è convergente se e solo
se è consistente.
Analizziamo più in dettaglio l’errore del metodo di Eulero
esplicito. Sia:
li 1  y ( xi 1 )  yi
l’errore globale, cioè la differenza tra la soluzione analitica e
quella approssimata, e sia:
yi*1  yi  hf ( xi , yi )
la soluzione ottenuta con un passo del metodo di Eulero
esplicito.
INTERPRETAZIONE GEOMETRICA
DEGLI ERRORI PER METODI AD UN
PASSO
Vediamo come è possibile controllare gli errori nei metodi ad un
passo.
Sia y(x) la soluzione di:
y’ = f(x,y)
y (x0) = y0
e sia yi una approssimazione ad y(xi) in qualche xi.
y(xi)-yi rappresenta l’errore globale che è quello che vogliamo
tenere limitato da una certa accuratezza e che è difficile da
stimare.
Ciò che vogliamo fare è controllare l’errore globale
controllando l’errore locale.
• INSERIRE DISEGNO
Vediamo cos’e l’errore locale.
yi 1  yi  h ( xi , y i , h)
ottenuto con il metodo.
u(x) è una curva integrale di y’=f(x,y) che passa per yi, quindi
soddisfa
u’ = f(x, u)
u(xi) = yi
L’errore locale è quindi: u(xi+1)- yi+1
Esso quindi ci dice quanto bene può essere seguita la curva u(x)
con un passo.
Pertanto l’errore globale e locale sono correlati da:
y(xi+1)-yi+1=[y(xi+1)-u(xi+1)]+[u(xi+1)-yi]
L’errore globale ha quindi due componenti:
i. y(xi+1)-u(xi+1) misura di quanto distano le due curve
integrali y(x) e u(x) (dipende quindi dalla ODE ed è legata
alla “Stabilità” del problema);
ii.
u(xi+1)-yi misura quanto bene il metodo risolve u’=f(x,u).
u(xi)=yi, è legato al metodo e può essere reso piccolo
aumentando l’accuratezza del metodo (decrescendo il
passo h oppure aumentando l’ordine del metodo).
Ricaviamo u(xi+1)-yi:
u(xi+1) = u(xi)+hΦ(xi, u(xi), h)+ht
ma u(xi)=y => yi+1 = u(xi)+hΦ(xi, yi, h)+ht
=> u(xi+1)-yi+1= ht
Supponiamo di usare due metodi, uno di ordine p e uno di ordine
q, e vediamo di stimare l’errore locale:
Partiamo per entrambi da (xi, yi):
yi+1 = yi+hФ1(xi, yi, h) p
y’i+1 = yi+hФ2(xi, yi, h) q
p<q
Si ha:
u(xi+1)=u(xi)+hФ1(xi, u(xi); h)+ht
u(xi+1)=u(xi)+hФ2(xi, u(xi); h) +ht’
yi+1-y’i+1=ht – ht’ ovvero
yi+1-y’i+1=ht+O(hq+1) =>
ht ~ yi+1-y’i+1 e quindi:
u(xi+1)-yi+1≈ y’i+1-yi+1
ASSOLUTA STABILITA’
Tale tipo di stabilità riguarda la propagazione degli errori dei
passi precedenti.
Un metodo è assolutamente stabile se, per h fissato, yi è
limitato per xi  .
Consideriamo il Problema Test:
 y '  y (t )

 y (0)  1
1)
Sol:
y(t) = et
se Re(λ)<0 => lim |y(t)| =0
t->0
DEF: Un metodo numerico è assolutamente stabile se yi,
soluzione di 1) è tale che:
yi  0 per ti  
2)
Poiché yi è funzione di hλ si ha:
Regione di Assoluta Stabilità ≡ {z = hλ є C vera la 2)}
METODO DI EULERO IN AVANTI
Applichiamo il metodo alla 1):
la 2) è vera se:
yi+1=yi+hλyi
y0=1 =>
yi=(1+hλ)i
|1+h λ|<1
R.A per Eulero in Avanti:
hλ є C- 0<h<(2/|λ|)
C-={z є C : Re(z)<0}
METODO DI EULERO INDIETRO
1
yi 
(1  h )i
Per ogni hλ non appartiene a {zєC: |z-1|<1}
METODI DI RUNGE-KUTTA
Tutti i metodi ad un passo possono essere dedotti, come già detto,
dallo sviluppo in serie di Taylor:
yi+1=yi+hTk(xi, yi; h)
dove:
Tk(xi, yi, h)=y’(xi)+h/2y’’(xi)+…+h^(k-1)/k! h(k)(xi)
Il calcolo delle derivate di f può essere oneroso. D’altronde, i
metodi visti precedentemente sono di basso ordine. Un buon
compromesso tra la semplicità dei metodi di basso ordine e la
serie di Taylor troncata ad un alto ordine è dato dai metodi di
Runge-Kutta.
Rispetto ai metodi Multi-Step, che vedremo più avanti, si ha lo
svantaggio che occorrono molte valutazioni della f per
raggiungere la stessa accuratezza. L’idea dei metodi di R-K è di
costruire formule del tipo:
yi+1= yi+hΦ(xi, yi, h)
Con Φ coincidente con Tk per un certo numero di termini senza
l’utilizzo esplicito delle derivate. Per un certo metodi di ordine K:
Φ(xi, yi, h)=A1f(θ1, γ1)+…+Akf(θk, γk)
Per il metodo di Eulero, che può essere interpretato come R-K del
1° ordine, il punto (θ1, γ1) ≡ (x0, y0)
Nei R-K del 2° ordine si hanno i punti:
(xi, yi), (xi+αh, yi+αhf(xi, yi)) =>
yi 1  yi  h[ A1 f ( xi , yi )  A2 f ( xi  h, yi  hf ( xi , yi ))]
Espandiamo f(xi+αh, yi+αhf(xi, yi)) attorno ad (xi,yi):
f ( xi  h, yi  hf ( xi , yi )) 
 f ( xi , yi )  hf x ( xi , yi )  hf y ( xi , yi ) f ( xi , yi )  O ( h 2 )
 yi 1  yi h[ A1 f ( xi , yi )  A2 ( f ( xi , yi )  f x ( xi , yi )  hf y ( xi , yi ) f ( xi yi ))]
da cui:
ma:
h
 ( xi , yi ; h )  ( A1  A2 ) f ( xi , yi )  [ f x  f y f ]
2
h
Tk ( xi , yi ; h )  f ( xi , yi )  [ f x  f y f ]
2
Quindi, perché si abbia Φ = T2 si deve avere:
A1+A2= 1
1
αA2 = 2
che danno luogo ad una famiglia di metodi R-K del 2° ordine.
I più noti di tali metodi sono quelli di Eulero modificato, di Heun
e di Raltson.
Eulero modificato: α =
1
2
A2 = 1 A1 = 0
h
h
y i1  y i  hf(x i  , y i  f(x i , yi ))
2
2
Che è equivalente a calcolare y
h
i
2
y
i
1
2
Calcolare la pendenza:
con Eulero:
h
 yi  f ( xi , yi )
2
y'
1
i
2
 f (x
1
i
2
,y
1
i
2
)
ed usarla per tutto l’intervallo yi 1  y1  hf ( x
1
i
2
Disegno
,y
1
i
2
)
METODI DI RUNGE-KUTTA A PASSO
VARIABILE
Poiché tali metodi sono ad un passo e’ semplice rendere tale
passo adattativo, cioè tale da ridurre l’errore.
Per ridurre l’errore e’ necessario poterlo stimare. Ciò può essere
fatto in due modi:
1. stesso metodo con due passi diversi (h, 2h)
2. due metodi di ordine diverso ma con lo stesso numero di
stadi
 yn 1  y ( xn 1 )   ( yn )h p 1  (h p  2 )
caso 1)
Metodo di ordine p. Partendo dal dato esatto: y(xn) = yn l’errore
locale sia minore di ε. Si ha:
 yn 1  y ( xn 1 )   ( yn )h p 1  (h p  2 )
dove Φ(yn) è una funzione incognita.
Stesso calcolo con passo 2h a partire da xn-1
 yn^1  y ( xn 1 )  ( yn )( 2h) p 1  (h p  2 )
Sottraendo:
(2 p 1  1) h p 1( yn )  yn 1  yn^1  ( h p  2 ) 
 y ( xn 1 )  yn 1
yn 1  yn^1
~

p 1
2 1
Se |ξ| < ε si prosegue altrimenti si dimezza il passo.
In generale, il passo raddoppia se


2 p 1
caso 2)
come già visto, usando 2 schemi di ordine p є p+1 la differenza tra
le soluzioni approssimate dà una stima dell’errore di troncamento
locale per lo schema di ordine inferiore.
Metodo di Heun, α=1, A1=A2=½
è usato per rendere esplicito il metodo dei trapezi
h
yi 1  yi  [ f ( xi , yi )  hf ( xi  h. yi  f ( xi , yi ))]
2
Metodo di Ralston, α=3/4, A1=⅓, A2=⅔
tale metodo da il minimo errore di troncamento
yi+1=yi+(h/3)(k1+2k2)
k1=f(xi, yi)
k2=f(xi+(3/4)h, yi+(3/4)hk1)
Metodi R-K espliciti generali
yi+1= yi+h
m
c  k
j 1
i
j
j 1
k1=f(xi, yi), kj=f(xi+ αjh, yi+h   jl  kl )
l 1
per j=2,…n
Il metodo più noto è quello del 4° ordine:
h
yi 1  yi  ( k1  2k2  2k3  k4 )
6
k1  f ( xi , yi )
h
h
k2  f ( xi  , yi  k1 )
2
2
h
h
k3  f ( xi  , yi  k2 )
2
2
k4  f ( xi  h, yi  h  k3 )
Metodi di ordine maggiore non sono convenienti poiché
richiedono un numero troppo grande di valutazioni della f.
METODI MULTISTEP
Integriamo la ODE tra tn-j e tn+k
y (t n  k )  y (t n  j ) 
t n k
 f (t , y(t ))dt
t n j
e applichiamo una quadratica di Newton-Cotes poiché
supponiamo la suddivisione uniforme dell’intervallo.
Utilizziamo q+1 punti: tn-q, tn-q+1, …, tn
e costruiamo il polinomio di Lagrange, integrando poi in [tn-j, tn+k]
pq ( x)  i 0 f (t n i , yn i ) Li ( x)
q
x  xn  l
Li ( x)  
l  0 x n  i  xn  l
q
l i
Integrando il polinomio si ha:
q
yn  k  yn  j  h  qi f n i
i 0
f l  f (tl , yl )
k
q
1
xl
 qi   Li (t )dt   
dx
q t n j
i l
 j l 0
t n k
l i
k, j  q determinano vari metodi multistep
k=1 j=0  ADAMS- BASHFORTH
k=0 j=1  ADAMS- MOULTON
k=1 j=1  NYSTRÖM
espliciti
impliciti
Metodi Adams-Bashforth (espliciti)
Sono basati sulla quadratura interpolatoria dell’integrale:
k=1, j=0
q
yn1  yn  hi 0  qi f ni
Se q=0  EULERO ESPLICITO : yn+1=yn+hfn
Metodi Adams-Moulton (impliciti)
Sono basati sulla quadratura interpolatoria dell’integrale:
q
k=0, j=1
yn  yn1  hi 0  qi f ni
preferibile riscriverlo come
yn1  yn  hi 0  qi f ni 1
q
Se q=0  EULERO IMPLICITO: yn+1=yn+hfn+1
h
y

y

[ f n  f n 1 ]
Se q=1  CRANK-NICHOLSON (TRAPEZI): n 1
n
2
Metodi di Nyström
k=1, j=1
yn1  yn1  hi 0  qi f ni
q
Se q=0  METODO DEL PUNTO MEDIO
yn1  yn1  2hf (tn , yn )
METODI PREDICTOR-CORRECTOR
Risolvendo un problema di Cauchy non lineare con uno schema
implicito è richiesto, ad ogni passo, la risoluzione di un’equazione
non lineare. Si possono usare: metodi di Punto Fisso, metodo di
Newton, … Ciò richiederà un guess iniziale vicino alla soluzione
sia per problemi di convergenza, sia per diminuire il numero di
iterazioni. Ciò può essere ottenuto usando in coppia un metodo
Ciò può essere ottenuto usando in coppia un metodo esplicito
(predictor) che fornisce un buon dato iniziale per il metodo
implicito (corrector) che è generalmente più stabile.
Un esempio di tale metodo è quello di Heun
h
yn 1  yn  [ f (t n , yn )  f (t n 1 , yn  hf n )]
2
in cui il predictor è Eulero in avanti e il corrector è il metodo di
Crank-Nicholson
P.
~
yn 1  yn  hf (t n , yn )
C.
h
yn 1  yn  [ f (t n , yn )  f (t n 1 , ~
yn 1 )
2
L’ordine di convergenza è q se p ha ordine q-1 e c ha ordine q.
Generalmente si usano i metodi di Adams in coppia (2-3, 3-4) per
ottenere PC di ordine pari a quello del corrector.
Metodi BDF (Backward Differentiating Formula)
Famiglia di schemi complementari a quelli di Adams. Lì si usa una
quadratura per approssimare l’integrale, nei BDF si approssima la y’.
Se si hanno q+1 punti e si conosce un’approssimazione della
soluzione nei punti n-q+1, …, n+1 si può determinare una pq la
cui derivata interpola la y’.
Calcoliamo la derivata in uno dei nodi tk
p’ (t k )  f(t k , yk )
Se k=n, il metodo è esplicito; se k=n+1, implicito. In generale:
q
l

 i yni 1  hf n
i 0
q
l

 i yni 1  hf n1
i 0
ESPLICITO
l  n  q  1,..., n  1
IMPLICITO
dove i coefficienti sono dati dalle derivate del polinomio di
Lagrange.
Es.:
q=1 Eulero avanti
q=2 Punto Medio
q=3 Instabile
LMM (LINEAR MULTISTEP METHODS)
Una generalizzazione dei metodi multistep che include i metodi di
Adams e i metodi BDF, è data dalla famiglia dei metodi multistep
lineari.
Un metodo multistep lineare ha la forma:
q
 y
i 0
i
q
n i
 h  i f n1
i 0
CONVERGENZA
L’analisi della convergenza è più complicata poiché:
1. La sol. approssimata è influenzata pure dagli errori nei valori
di partenza:
lj = yj - y(xj) j=0,..,k-1
Tali valori si dicono CONSISTENTI se
lim y j (h)  y ( x j )  0
h 0
j=0,..,k-1.
2. I metodi M-S possono essere instabili. Per mostrare ciò,
vediamo un esempio:
 y'  y

 y  
Soluzione: y(x)=ex. Tale problema è detto Problema Test.
Per >0 il problema è instabile. Analizziamo il comportamento di
qualche metodo multistep nel caso <0.
Consideriamo il metodo del Punto Medio
yn  2  yn  2hf n 1
che ha l’equazione, alle diff. associata:
x 2  2hx  1  0
le cui soluzioni sono:
r1  h  1  h 2 2
r2  h  1  h 2 2
che ha sol. generale
yn  1r1n   2 r2n
Ricaviamo 1 e 2 :
 y 0  1  1   2

h
y

e
 1r1   2 r2
 1
 1 
e h  r2
 1  O(h 2 2 )
2 1  h 2 2
r1  e h
2 
 O(h 33 )
2 1  h 2 2
per h0,
per >0,
per <0,
11, 20
| r1|>| r2|>0
termine dominante: 1 r1n
0<r1<1, r2<-1 termine dominante: 2 r2n
Pertanto per <0 la soluzione diverge da quella vera.
Questo perché la ODE ha una sola soluzione mentre l’equazione
alle differenze di ordine K ha k soluzioni di cui una corrisponde
alla vera soluzione per aversi convergenza è quindi necessario
che le altre soluzione rimangano limitate. Analizziamo quindi il
comportamento delle equazioni alle differenze relativamente al
problema della stabilità.
DEF: L’equazione alle differenze
k 1
Z n  k   am z n  m  0
n  0,1,...
m 0
con coefficienti a0,…,ak-1 costanti è detta stabile se tutte le sue
soluzioni sono limitate.
Per cercare delle condizioni facilmente verificabili per stabilire la
convergenza di un metodo MS partiamo dall’errore locale di
discretizzazione
k
k
h    j y ( x  hj )  h  j y ' ( x  hj )
j 0
j 0
Abbiamo visto che il metodo è consistente se 0 per h0.
E’ detto di ordine p se:
h  O(h p 1 )
Se y(x) è sufficientemente differenziabile si può esprimere h come
h  c0 y ( x)  c1hy ' ( x)  ...  c p h p y ( p ) ( x)  ...
Infatti espandendo y(x+hj) e y’(x+hj) attorno x si ha
(hj ) 2
y ( x  hj )  y ( x)  hjy' ( x)  y ' ' ( x)
 ...
2
(hj ) 2
y ' ( x  hj )  y ' ( x)  hjy' ' ( x)  y ' ' ' ( x)
 ...
2
che dà h se poniamo
c0   0  1  ...   n
c1  1  2 2  ...  k k  (  0  1  ...   k )
...
1
1
n
n
cn  (1  2  2  ...  k  k ) 
( 1  2 n 1  2  ...  k n 1  k )
n!
(n  1)!
...
Se c0=c1=…..=cp=0, cp+10, il metodo è di ordine p.
Vediamo le proprietà di un metodo convergente.
Se il metodo multistep converge, c0 deve essere nullo.
Sia dato il problema
 y'  0

 y (0)  1
sol: y(x)=1
k
   j yn  j  0
j 0
Fissiamo x e verifichiamo n ed h:
x  (n  k )h  x0

0
Supponiamo il metodo convergente (non alla soluzione y=1):
yn  k  y ( x ) per h0
y n  s  y ( x )   j ( h)
0jk k fissato
lim 
h 0
j
( h)  0
0
  j y ( x )    j j ( h )  0
k
k
j 0
j 0
k
k
j 0
j 0
0    j y ( x )  y ( x )  j
 c0  0
Dimostriamo ora che un metodo convergente alla soluzione ha
ordine almeno 1.
Sia dato il problema:
 y'  1

 y (0)  0
sol: y(x)=x
k
k
j 0
j 0
   j yn  j  h  j
Una soluzione è data da :
k
yi  ihM con M 

j 0
k
j
 j
j 0
j
Se x  ( n  k ) h
yn  k  (n  k )hM
E poiché la soluzione è y(x)=x :
k
k
j 0
j 0
M  1    j  j j  c1  0
Un metodo che è almeno di ordine 1 è detto consistente allora una
condizione necessaria per la convergenza è la consistenza ma essa
non è sufficiente. Solo se anche la condizione della radice è
soddisfatta allora si ha convergenza.
Infatti, se il metodo è convergente, lo è pure per il problema
 y'  0

 y (0)  0
sol.: y(x)=0
k
   j yn  j  0
j 0
che è soddisfatta da:
ym  h(ri ) m
dove ri è soluzione del polinomio caratteristico.
Poiché si abbia convergenza, si deve avere:
h0
yn  k  y ( x ) ma yn+k=h(ri)n+k
yn  k  y ( x )  0  ri  1
Se ri non è uno zero semplice, ma ha molteplicità m:
y j  h qj (ri ) j
j  0,1,... q  m  1
Per j = n+k:
yn  k  h(n  k ) q (ri ) n  k ma
h( n  k )  x
n
yn  k  0  ri  1
La condizione della radice è:
1) ri  1 se ri è uno zero semplice del polinomio caratteristico
2) ri  1 se ri non è uno zero semplice del polinomio caratteristico
Per un metodo consistente, il polinomio caratteristico ha una
radice r1=1 detta radice principale. Infatti, in tal caso
k
c0  0   j  p(1)  0
j 0
I metodi
yn  2  yn  2hf n 1
h
yn  2  yn  [ f n  4 f n 1  f n  2 ]
3
hanno p(x) = x2-1 e quindi soddisfano il criterio della radice.
Sono consistenti e quindi convergenti. Eppure non sono buoni da
usare in pratica. Abbiamo già visto che per i MS non basta la sola
convergenza poiché le equazioni alle differenze hanno soluzioni in
più rispetto alla ODE. Tali soluzioni, dette Parasitiche, devono
rimanere piccole rispetto alla radice principale e ciò porta al
concetto di stabilità relativa.
Applicando il metodo al problema test:
 y '  y

 y ( 0)  1
si ha
sol: y(x)=ex
ym   1e mh  2(r2 ) m  ....  k(rk ) m   (h p 1 )
Soluzione del Problema Test:
y ( xm )  y0 e mh
Quindi ym è una buona approssimazione di y(xm) se:
1) 1y0, i0
i=2,…,k
2) ri<<ex
i=2,…,k
•
•
E’ soddisfatto se i valori di partenza sono buoni
E’ relativo alla stabilità relativa
Un metodo MS si dirà relativamente stabile se:
|r1|>|ri|
i=2,..,k
L’intervallo di stabilità relativa è il più grande intervallo (,),
 tale che il metodo è R.S.  h  (,). Se  è grande, h
dovrà essere piccolo. Con tale tipo di stabilità si controlla l’errore
relativo.Infatti:
m
m
 k  rk 
 2  r2 
ym   1 (r1 ) [1     .....    ]
 1  r1 
 1  r1 
m
Stabilità assoluta
Spesso è importante fare un’analisi di stabilità tenendo il passo h
fissato, e ciò permette di controllare l’errore assoluto, un metodo è
assolutamente stabile se gli errori ai passi precedenti non
aumentano.
Tale concetto si applica anche ai metodi onestep, come abbiamo
già visto applicando il metodo M-S
k 1
k
j 0
j 0
yn k   j yn j  h  j f ( xn j , yn j )
al problema test
 y ' (t )  y (t )

 y (0)  1
t>0
Si ha:
k 1
k
j 0
j 0
yn  k    j yn  j  h   j yn  j
k
 (
j 0
j
 h j ) yn  j  0
Per x=xn, si ha:
k
 (
j 0
j
 h j ) y ( xn  j )  h n
j
 h j )ln  j  h n
k
e sottraendo
 (
j 0
E quindi gli errori soddisfano un’equazione alle differenze le cui
soluzioni sono:
lm  1 (r1 )  ....   k (rk ) 
m
m
h
k
h   j
j 0
Diremo che un metodo M-S soddisfa la condizione assoluta delle
radici se esiste h0>0:
|rj(h)|<1
j=0,..,k hh0
Pertanto C.N.S. affinché un metodo M-S sia assolutamente stabile,
ovvero che
|yn|0 per tn 
è che esso soddisfi la condizione assoluta delle radici.
L’assoluta stabilità implica la zero stabilità, mentre il viceversa
non è vero.
Scarica

ODE