Calcolo Numerico 2

Calcolo Numerico II

1

Indice

1 Introduzione alle equazioni differenziali 3

2 Metodi ad un passo 4

2.1 Metodi di Eulero . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1.1 Codici d’implementazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.2 Costruzioni di altri metodi generici . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2.1 Codici d’implementazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3 Metodi di Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.3.1 Codici d’implementazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.4 Metodi di Estrapolazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.5 Stime a posteriori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.5.1 Codici d’implementazione con passo adattativo . . . . . . . . . . . . . . . . 28

2.6 Equazioni alle differenze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3 Introduzione ai Metodi multistep 37

3.1 Metodi a più passi lineari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.1.1 Codici d’implementazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.2 Metodi di Adams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.3 Metodi predictor-corrector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.3.1 Codici d’implementazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4 Problema ai valori al contorno 50

5 Metodi risolutivi 54

5.1 Metodo Shooting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.2 Metodo delle differenze finite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

5.2.1 Codici d’implementazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.3 Metodo delle differenze finite per problemi non lineari . . . . . . . . . . . . . . . . 73

5.3.1 Codici d’implementazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5.4 Metodo degli elementi finiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

2

1 Introduzione alle equazioni differenziali

Le equazioni differenziali sono equazioni, la cui incognita non è uno scalare o un vettore, ma una

~t

funzione ~x = ~x ( ). ~

~t, (k)

0

Una generica equazione differenziale appare nella forma G( ~x, x , . . . , ~x ) = 0

Definizione 1.1. (Equazione ordinaria) ∈

Si tratta di una equazione differenziale la cui incognita è ~x = ~x (t) con t R

Definizione 1.2. (Forma ordinaria normale di ordine k)

Si tratta di una equazione differenziale ordinaria, che può venir riscritta nel seguente modo

~ ~

(k) (k−1)

0 00

~x = f (x, x , x , . . . , ~x , t)

Nota.

D’ora in avanti non facciamo più distinzione tra ~x e x per poter alleggerire le notazioni.

Definizione 1.3. (Equazioni del primo e del secondo ordine)

Ci riferiamo a queste equazioni:

0

x = f (x, t)

00 0

x = f (x, x , t)

Definizione 1.4. (Soluzione)

nk n n

⊂ × → ⊂ →

Sia f : E , allora x : I è soluzione su I di un equazione differenziale

R R R R R

del primo o secondo ordine (a seconda che m sia 2 o 3) se

1. x è derivabile k volte 0

∀x ∈ ∈

2. I, (x, t) oppure (x, x , t) E

∀x ∈

3. I l’equazione differenziale è soddisfatta

Definizione 1.5. (Problema di Cauchy (PC) - Problema ai valori iniziali (PVI))

Si tratta di una particolare equazione differenziale scritta nella forma

0

(k) (k−1)

 x = f (x, x , . . . , x , t)





 x(t ) = x

 0 0





 0

x (t ) = x

0 1

..



 .







 (k−1)

 x (t ) = x

0 k−1

Teorema 1.6. (Teorema di Peano)

nk n 0

⊂ × → ∈ ∈

Sia f : Ω C (Ω) con (x, t) Ω allora il PC ha almeno una soluzione.

R R R

Teorema 1.7. (Teorema di esistenza e unicità globale del PC)

nk n 0

⊂ × → ⊂ ∈

Sia f : Ω con J intervallo f C (Ω) e lipschitziana rispetto a

J R R R nk

0 (k−1) × ∀I ⊂

(x, x , . . . , x ), uniformemente rispetto a t su I con J compatto, allora

R

nk

∀(x ∈ ∗

, x , . . . , x , t) J il Problema di Cauchy ha 1! soluzione su J

R

0 1 (k−1)

Teorema 1.8. (Teorema di esistenza e unicità locale del PC)

nk+1 n 1

⊂ → ∈ ∈

Sia f : Ω con f C (Ω) e (x , x , . . . , x , t) Ω, allora

R R 0 1 (k−1)

∃ U(t U(t

) t.c. il Problema di Cauchy ha 1! soluzione su )

0 0

3

2 Metodi ad un passo

2.1 Metodi di Eulero

Lemma 2.1. (Lemma di Gronwall)

t

Z + BT

∈ ∈ ≤

≤ ≤ v(s)ds con t (t , t + T ), A, B, T . Allora v(t) Ae

Sia 0 v(t) A + B R

0 0

t

0

Nota. t

Z 0 B(t−t )

⇒

Se fosse v(t) = A + B v(s)ds avrei v (t) = Bv(t), v(t ) = A v(t) = Ae 0

0

t 0

Dimostrazione. t

Z v(s)ds. Segue che

Sia z(t) = A + B t 0 ( 0 ≤

z (t) Bz(t)

0 ≤ ⇒

z (t) = Bv(t) Bz(t) ≤

z(t ) A

0

0

⇒ − ≤

z (t) Bz(t) 0

0 −B(t−t −B(t−t

) )

− ≤

z (t)e Bz(t)e 0

0 0

−Bt 0 ≤

(z(t)e ) 0

Si ha che quindi la funzione è monotona decrescente e di conseguenza

−B(t−t ) ≤

z(t)e A

0

B(t−t )

≤

z(t) Ae 0 B(t−t ) BT

≤ ≤ ≤

v(t) z(t) Ae Ae

0

Lemma 2.2. (disuguaglianze e identità)

2 2 2

≤ ∀η − ≥

1. 2ab ηa + b /η, > 0 poiché si ha η(a b/η) 0

1 2 2 2

− −

2. a(a b) = (a + (a b) + b )

2

2 2 2

≤ ≤

3. Se a b + c , allora a b + c (disuguaglianza triangolare)

Lemma 2.3.

|e | ≤ | |

Da a|e + b|σ si ottiene in modo ricorsivo

n+1 n n n −

a 1

X

n i n

|e | ≤ |e | |σ |) ≤ |e | |)

a + b (a a + b max(|σ (2.1.1)

n+1 0 n−i 0 n

−

a 1

Teorema 2.4.

( (

0 0

y = f (t, y) ỹ = f (t, ỹ) + (t)

Dato un PVI se considero il PVI, con f (x, y) lipschitziana

g

y(t ) = y ỹ(t ) = y +

0 0 0 0 0

rispetto a y, uniformemente rispetto a t allora questo problema non è ben posto.

Dimostrazione.

Se considero t

Z

|y(t) − − − −

ỹ(t)| = (f (s, y(s)) f (s, ỹ(s)) (s)ds) 0

t

0

t t

Z Z +

≤ |f − |(s)|ds | | ∈ ∈

(s, y(s)) f (s, ỹ(s))|ds + + con t (t , t + T ), T R

0 0 0

t t

0 0

t t

Z Z

· |y(s) −

| | |(s)|ds + L ỹ(s)| ds

= +

0 t

t {z }

|

0 0 v(s)

| {z }

A 4

⇒ Per il lemma di Gronwall

t LT

R

|y(t) − ≤ | |(s)|ds) ·

ỹ(t)| (| + e

0 t

0

Definizione 2.5. (Problemi dissipativi)

Diciamo che f è un problema dissipativo se la funzione f è monotona non crescente nella variabile

y.

Nota.

Se f è monotona decrescente, f (0) = 0 si ha

0

y = f

0 · · ≤

y y = f y 0

0

1 2 ≤

(y ) 0

2

Teorema 2.6.

Dato un problema dissipativo, la maggiorazione dell’errore è minore di un problema non dissipa-

tivo.

Dimostrazione.

Dato

( 0 0 0

− −

δ (t) = ỹ (t) y (t) = f (t, ỹ(t)) f (t, y(t)) + (t)

δ(t ) =

0 0

0

⇒ − · −

δ̃ (t)δ(t) = [f (t, ỹ(t)) f (t, y(t))] [ỹ(t) y(t)] +(t)δ(t)

| {z }

≤0

0

⇒ ≤

δ̃ (t)δ(t) (t)δ(t)

0

12 2

⇒ ≤

( δ̃ (t)) (t)δ(t)

Posso allora integrare sia a destra, sia a sinistra, mantenendo la disequazione.

t

Z

1

1 2 2

− ≤

δ (t) δ (t ) (s)δ(s)ds

0

2 2 t

| {z } 0

20

t t t

Z Z Z

2 20 20 20

≤ ≤ |(s)| · |δ(s)| ≤ |δ(s)| |(s)|ds

δ (t) +2 (s)δ(s)ds +2 ds + 2 max

x∈(t ,t)

t t t

0

0 0 0

t

Z

20 |δ(s)| |(s)|ds

A questo punto ho che + 2 max è monotona crescente

x∈(t ,t) t

0 0 t

Z

2 20

⇒ |δ ≤ |δ(s)| |(s)|ds

max (t)| + 2 max

x∈(t ,t) x∈(t ,t) t

0 0 0

t

Z 2 20

|δ(t)|, |(s)|ds ≤

Pongo a = max b = e quindi ho a + 2ab

x∈(t ,t) t

0 0 2

2 a

a 2 2 20 2

⇒ ≤

≤ +

+ 2b a + 2b

Inoltre 2ab 2 2 2

t

2

|δ Z

max (t)|

x∈(t ,t)

2 20 0

|δ ≤ |(s)|ds

max (t)| + +2

2

x∈(t ,t) t

0 0

cioè 2

t

Z

2 20

|δ ≤ |(s)|ds

max (t)| 2 + 4

x∈(t ,t) t

0 0

e quindi √ t

Z

|δ(t)| ≤ |(s)|ds

max 2 + 2 (2.1.2)

0

x∈(t ,t) t

0 0

5

Definizione 2.7. (Metodo di Eulero esplicito)

( 0 ∈

y = f (t, y(t)), t [t , t + T ]

0 0

Dato il PC y(t ) = y

0 0 0

assegniamo ad ogni punto nel piano una velocità. Noi sappiamo che f (t , y ) = y (t ), quindi,

0 0 0

−

in U (t , y ) se considero la retta tangente y + f (t , y )(t t ) al posto della funzione si avrà un

0 0 0 0 0 0

errore o(h) e l’idea è di riconsiderare il PC in t + h = t , chiamiamo quindi Y = y + hf (t , y ) e

0 1 1 0 0 0

ricorsivamente Y = Y + hf (t , Y ). Questo metodo è chiamato metodo di Eulero esplicito.

n n−1 n n−1

Nota.

Quando si utilizza il metodo di Eulero esplicito, (come da qualsiasi altro metodo) si commettono

due tipi di errori differenti.

1. Il primo errore lo si commette muovendosi sulla retta tangente al posto che sulla funzione

in un intorno della soluzione

2. Il secondo errore è che consideriamo al passo α il valore f (t , Y ) che non è tangente alla

α α

prima curva!

Un modo alternativo per costruire il metodo di Eulero esplicito è notare che

−

y(t + h) y(t)

0 + '

y (t) = f (t, y(t)) h

·

Ponendo allora y(x + n h) = Y si ha

n −

Y Y

n+1 n

'

f (t , Y )

n n h

Dal quale si scrive subito Y = Y + hf (t , Y )

n+1 n n n

Di seguito diamo le definizioni dei seguenti errori.

Definizione 2.8. (Errore locale di troncamento) 1 si ha

Si tratta del primo errore (il muoversi sulla retta tangente) e riscalando per h

1 −

δ (h) = (y(t ) [y(t ) + hf (t , y(t ))])

n n+1 n n n

h

Esso è quindi l’errore in una iterazione supponendo di avere la soluzione esatta, nel caso partico-

lare usando eulero esplicito.

Definizione 2.9. (consistenza) −−−→ ∀n

Diciamo che un metodo è consistente se δ (h) 0

n h→0 α α

→ ·

Diciamo che un metodo è di ordine α se per h 0 si ha δ (h) = c h = O(h )

n

Nota.

Il metodo di Eulero è consistente con ordine 1, infatti applicando Taylor si ha:

2 2

0 00 0 00

h h 2 2

· · · −

hδ (h) = y(t ) + hy (t ) + y (t ) + y(t ) + hy (t ) = y (t ) + o(h ) = O(h )

n n n n n n n

2 2

1

E quindi δ (h) = O(h )

n

Definizione 2.10. (Errore globale di troncamento)

Si tratta dell’errore:

−

e = y(t ) Y

n n n 6

Definizione 2.11. (Metodo convergente) → ∀n →

Diciamo che un metodo è un metodo convergente se e 0, per h 0

n

Teorema 2.12.

Il metodo di Eulero esplicito è convergente.

Dimostrazione.

Possiamo riscrivere il PC (o PVI) in forma discreta; si ha, scrivendo y al posto di y(t ) e σ al

n n n

posto di σ (h)y = y + hf (t , y) + hσ e Y = Y + hf (t , Y ).

n n+1 n n n n+1 n n n

−

Se poniamo e = y Y allora otteniamo

n+1 n+1 n+1

− ⇒ |e | ≤ |e | − |

e = e + h[f (t , y ) f (t , Y )] + hσ + h|L(y Y )| + h|σ = (1 +

n+1 n n n n n n n+1 n n n n

| |

hL)|e + h|σ

n+1 n n

X

n+1 j

⇒ |e | ≤ | | ≤ · · · ≤ |e | |σ |

+h (1 + hL)

(1 + hL)|e + h|σ (1 + hL)

n+1 n n 0 n−j

{z }

| j=0

=0

|σ

Sia allora σ = max (h)| si ha che

n n+1 −1

(1+hL) σ σ σ σ

n+1 hL n+1 hL(n+1) LT

j

P ≤ ≤ ≤

|e | ≤ (1 + hL) (e ) = e e

hσ (1 + hL) = hσ

n+1 1+hL−1 L L L L

00

h

Poiché σ = )) si ha che il metodo è convergente di ordine 1.

M (M = max(y

2 2

2

Definizione 2.13. (Problemi di Stabilità)

Diciamo che un metodo, quando esso viene applicato ad un problema dissipativo lineare, cioè

(n) (i)

P

nella forma y = λ y , ha un problema di stabilità o che è instabile, se scelto il passo finito

i

h (dal punto di vista teorico esso è infinitesimale, ma dal punto di vista pratico, su un calcolatore,

∞,

no) la soluzione Y non tende a 0 quando n tende a ma vi è un errore che viene propagato al

n

crescere di n (e quindi di t ).

n

Un metodo che non dà problemi di stabilità sotto alcune ipotesi viene chiamato assolutamente

stabile.

Nota. Ovviamente è possibile definire e controllare la stabilità su equazioni differenziali non

lineari, ma la sua definizione sarebbe diversa da quella appena data e l’instabilità difficile da

controllare.

Definizione 2.14. (A-stabilità)

Diciamo che un metodo è A-stabile se la sua regione di stabilità non dipende dal passo temporale.

Esempio 2.15. (Problema modello con Eulero esplicito)

Dato il problema modello: ( 0

y = λy, Re(λ) < 0

y(0) = 1

λt

Banalmente si ha che y(t) = e , se però applichiamo il metodo di Eulero esplicito si ha: Y =

n+1

n n

· · ·

Y + hλY = (1 + hλ)Y = = (1 + hλ) Y = (1 + hλ) . Se h > 0 è fissato, si ha che se

n n n 0

|1 → ∞

+ hλ| > 1 allora per n Y 0 al contrario di y(t ). Questo è un problema di stabilità,

9

n n

perchè il metodo di Eulero esplicito non converge alla soluzione poiché h è troppo grande, nello

−2

|1 ⇒ −2 ⇒

specifico serve + hλ| < 1 < Re(hλ) < 0 h < Re(λ)

Nota. k k

→

Il caso descritto nell’esempio si può generalizzare facilmente; data y : e l’equazione

R R

( 0 ∈

y (t) = Ay(t), A M (k) simm., def. negativa

y(0) = y 0 7

Allora si deve scegliere il passo h nel seguente modo.

−1

∃

Si ha che X t.c XAX = ∆ diagonale con λ autov. di A.

i

−1 0 −1 −1 0 −1

⇒ ⇒ ⇒

A = X∆X y = Ay = X∆x y (X y ) = ∆(X y)

−1

Ponendo z = X y si ha z = ∆z, si ha quindi, su ogni singola componente, una equazione del

i

0

tipo z = f (z )

i i

i

Possiamo quindi applicare il metodo di Eulero esplicito, e per evitare problemi di stabilità dovre-

2

|1 | ∀i

mo scegliere h t.c. + hλ < 1 = 1 . . . k, cioè h <

i |

max|λ

i

Nel caso invece λ fosse a sua volta in funzione di t, sarebbe allora necessario scegliere h tale che

2

≤

h max(λ(t))

Definizione 2.16. (Regione di assoluta stabilità) ·

Si tratta di un sottoinsieme di i cui elementi z = h λ non causano problemi di stabilità al

C

metodo utilizzato.

Definizione 2.17. (Metodo di Eulero implicito)

Per ottenere il metodo di Eulero esplicito abbiamo considerato l’approssimazione della derivata

prima da sinistra, se invece consideriamo la seguente approssimazione

−

Y Y

n+1 n

'

f (t , Y )

n+1 n+1 h

Otteniamo il metodo iterattivo Y = Y + hf (t , Y )

n+1 n n+1 n+1

Che viene chiamato implicito perchè la soluzione Y dipende da Y stessa. A questo punto

n+1 n+1

si possono usare diversi metodi per risolvere il metodo implicito, come ad esempio porre

− −

Y Y hf (t , Y ) = g(Y )

n+1 n n+1 n+1 n+1

Usando il metodo di bisezione, che però è ”lento”, oppure il metodo di Newton, notando che

−Y

Y

0 0

'

n+1 n

y . Ma se si avesse y 1 o se non fossimo abbastanza vicini al punto, il metodo di

h

Newton potrebbe non convergere.

Un’altra strada è quella di porre Y = Φ(Y ) = Y + hf (t , Y ) valutare quindi il pro-

n+1 n+1 n n+1 n+1

blema come problema di punto fisso.

0 0

Si nota subito che Φ (S) = hf (t , S) e che per la convergenza di questo metodo serve

n+1

y 1 perché Φ sia una contrazione

(condizione sufficiente, ma non necessaria) h < 0

f

y

Nota.

In questo caso il vincolo da porre ad h, in confronto al metodo di Newton implicito, è due volte

maggiore, serve un passo h minore in confronto, la condizione è però sufficiente, non necessaria.

Consistenza

Si ha che il metodo di eulero implicito e’ consistente, infatti si ha

− −

hσ (h) = y(t ) y(t ) h f (t , y(t )

)

n n+1 n n+1 n+1

| {z }

0

y (t )

n+1

2 2

0 00 0 00 00

h h

· · · − − −

hσ (h) = y + hy + y + y h(y + hy + . . . ) = y

n n n

n n n n n

2 2

Convergenza

Si puo’ dimostrare che il metodo di Eulero implicito converge con ordine 1.

Sia y(t ) = y(t ) + hf (t , y(t )) + hσ e Y = Y + hf (t , Y ), posto allora

n+1 n n+1 n+1 n n+1 n n+1 n+1

8

−

e = y(t ) Y si ha

n+1 n+1 n+1 − ⇒ |e | ≤ |e | − |

e = e + h[f (t , y ) f (t , Y )] + hσ