METODO DI EULERO ESPLICITO

METODO DI EULERO ESPLICITO
u0 dato
un+1 = un + hf (tn , un ) 0 ≤ n ≤ Nh − 1
(1)
Scrivere una function
[tn,un]=eulero esp(odefun,tspan,y0,Nh) con
INPUT:
odefun: espressione della f
tspan=[t0,T]: vettore di 2 componenti: istante iniziale e finale
dell’intervallo
y0: valore scalare: la condizione iniziale
Nh: numero (intero) di passi temporali (Nh `e tale che T = tNh ).
OUTPUT:
tn: vettore colonna contenente gli istanti temporali da t0 a tNh .
un: vettore colonna contenente la soluzione numerica negli istanti
temporali tn .
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
1
N.B. esprimere f in funzione di due variabili t e y con y non
dipendente da t.
ES: f (t, y ) = t − y (t): diventa f=@(t,y)[t-y]
N.B. Generare il vettore dei tn con il comando linspace:
tn=linspace(tspan(1),tspan(2),Nh+1);
Inizializzare il vettore un della stessa dimensione di tn
e memorizzare nella prima componente di un il valore di y0 .
Con un ciclo costruire sequenzialmente le componenti del vettore
un
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
2
Esercizio 1 (esode1.m)
Si consideri il problema di Cauchy
′
y = t − y t ∈ (−1, 5]
y (−1) = 1.
(2)
Scrivere un m-file che:
1) richiami la function eulero esp.m con h = 0.001 (calcolare
Nh = (T − t0)/h e prendere la parte intera (usare fix, ceil,
round));
2) rappresenti graficamente la soluzione numerica e quella esatta
(fuori dalla function eulero esp.m).
La soluzione esatta `e: y (t) = t − 1 + 3e −(t+1) .
3) calcoli l’errore
eh =
max |y (tn ) − un |
n=0,...,Nh
(fuori dalla function eulero esp.m)
Sol. Si ottiene err=5.541e-4
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
3
I grafici della soluzione numerica e di quella esatta sono:
Zoom
4.5
sol numerica
sol esatta
sol numerica
sol esatta
0.24
4
0.22
3.5
0.2
3
0.18
2.5
0.16
2
0.14
1.5
0.12
1
0.1
0.5
0
−1
0
1
2
3
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
4
0.08
5
−0.5
0
0.5
1
4
Esercizio 2 (esode2.m)
Scrivere un m-file che, con un ciclo su h risolva il problema di
Cauchy dell’Es. 1. Si prendano questi valori di h:
h = 1, .8, .5, .1, .05, .01, .005, .001.
Aprire una finestra grafica e plottare la soluzione esatta.
Per ogni valore di h:
- richiamare la function eulero esp per calcolare [tn,un],
- plottare la soluzione numerica (con un colore diverso) sulla stessa
figura sui cui `e stata disegnata la soluzione esatta,
- calcolare l’errore
eh =
max |y (tn ) − un |
n=0,...,Nh
e memorizzare tali errori in un vettore.
Ricordare che la sol. esatta `e: y (t) = t − 1 + 3e −(t+1) .
Fuori dal ciclo su h, rappresentare in scala loglog gli errori eh
trovati in funzione degli h utilizzati e verificare che Eulero esplicito
`e un metodo convergente del primo ordine, cio`e che l’errore eh `e
proporzionale ad h(cio`e esiste C > 0 t.c. eh = C · h) quando
h → 0.
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
5
Plottando su uno stesso grafico la soluzione esatta e la soluzione
numerica dei problemi dell’esercizio 2 per alcuni valori di h si ha:
5
4
sol esatta
h=1
h=0.8
h=0.5
h=0.1
0.6
0.4
3
0.2
0
2
−0.2
1
−0.4
−0.6
0
−0.8
−1
−1
0
1
2
3
4
5
−1
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
sol esatta
h=1
h=0.8
h=0.5
h=0.1
h=0.05
h=0.01
h=0.005
0.8
1
Tanto pi`
u h `e piccolo, tanto pi`
u la soluzione numerica `e vicina alla
soluzione esatta (CONVERGENZA).
Per h → 0, maxn |un − y (tn )| → 0.
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
6
L’errore in funzione di h `e:
1
10
Errore Eulero
h
0
10
−1
errore
10
−2
10
−3
10
−4
10
−3
10
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
−2
−1
10
10
0
10
h
7
METODO DI EULERO IMPLICITO
u0 dato
un+1 = un + hf (tn+1 , un+1 ), 0 ≤ n ≤ Nh − 1
(3)
[tn,un]=eulero_imp(odefun,tspan,y0,Nh)
INPUT:
odefun: espressione della f
tspan=[t0,T]: vettore con istante iniziale e finale
dell’intervallo
y0: valore scalare: la condizione iniziale
Nh: numero intero di passi temporali (Nh `e tale che T = tNh ).
OUTPUT:
tn: vettore colonna contenente gli istanti temporali da t0 a
tNh .
un: vettore colonna contenente la soluzione numerica negli
istanti temporali tn .
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
8
Ad ogni passo tn dobbiamo risolvere l’equazione non lineare con
incognita un+1 :
un+1 = un +hf (tn+1 , un+1 ) ⇔ F (un+1 ) = un+1 −un +hf (tn+1 , un+1 ) = 0
Attenzione: la funzione F cambia ad ogni passo tn perch`e dipende
da un e da tn+1 .
L’equazione non lineare F (un+1 ) = 0 pu`o essere risolta con secanti.
Per semplicit`a poniamo w = un+1 .
Secanti per risolvere F (w ) = 0 `e:

 w (0) , w (1) dati
F (w (k) ) − F (w (k−1) )
 w (k+1) = w (k) −
w (k) − w (k−1)
k≥1
Utilizziamo la function secant
(www.ing.unibs.it/gervasio/CS/matlab)
[zero,res,niter]=secant(fun,w0,w1,tol,nmax)
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
9
Cosa prendiamo come w0 e w1?
Al primo passo temporale poniamo w (1) = u0 e w (0) = u0 + h.
Per i successivi passi temporali poniamo w (1) = un e w (0) = un−1 .
n=1
F=@(w)w-un(n)-h*odefun(tn(n+1),w);
un(n+1)=secant(F,un(n)+h,un(n),tol,nmax);
for n=2:Nh
F=@(w)w-un(n)-h*odefun(tn(n+1),w);
un(n+1)=secant(F,un(n-1),un(n),tol,nmax);
end
Esercizio 3. Ripetere lo stesso lavoro svolto nell’esercizio 2 per
testare la correttezza della function eulero imp.m.
Si riesce a dimostrare che il metodo di Eulero implicito ha ordine di
convergenza 1.
Verificare numericamente questo risultato teorico.
c
Paola
Gervasio - Calcolo Scientifico - 2014/15
10

Download Report