Calcolo Numerico 2

Calcolo numerico II

Indice

• Introduzione alle equazioni differenziali 3
• Metodi ad un passo 4
• Metodi di Eulero 4
• Codici d’implementazione 11
• Costruzioni di altri metodi generici 13
• Codici d’implementazione 16
• Metodi di Runge-Kutta 18
• Codici d’implementazione 20
• Metodi di Estrapolazione 23
• Stime a posteriori 25
• Codici d’implementazione con passo adattativo 28
• Equazioni alle differenze 33
• Introduzione ai Metodi multistep 37
• Metodi a più passi lineari 39
• Codici d’implementazione 43
• Metodi di Adams 45
• Metodi predictor-corrector 47
• Codici d’implementazione 48
• Problema ai valori al contorno 50
• Metodi risolutivi 54
• Metodo Shooting 54
• Metodo delle differenze finite 55
• Codici d’implementazione 65
• Metodo delle differenze finite per problemi non lineari 73
• Codici d’implementazione 77
• Metodo degli elementi finiti 79

1 Introduzione alle equazioni differenziali

Le equazioni differenziali sono equazioni, la cui incognita non è uno scalare o un vettore, ma una funzione ~x = ~x(t).

Una generica equazione differenziale appare nella forma G(~x, x^(k)₀, . . . , ~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 ~x^(k) = f(x, x^(k−1), 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:

x⁰ = f(x, t)

x⁰⁰ = f(x, x⁰, t)

Definizione 1.4. (Soluzione)

Sia f : E_n^k × R → Rⁿ, allora x : I → Rⁿ è soluzione su I di un’equazione differenziale del primo o secondo ordine, a seconda che m sia 2 o 3, se:

• X è derivabile k volte.
• ∀x ∈ I, (x, t) oppure (x, x⁰, t) ∈ E.
• ∀x ∈ 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:

x^(k) = f(x, x⁰, . . . , x^(k−1), t)

x(t₀) = x₀

x⁰(t₀) = x₁

. . .

x^(k−1)(t₀) = x_k−1

Teorema 1.6. (Teorema di Peano)

Sia f : Ω ⊂ R^nk × R → Rⁿ con f ∈ C⁰(Ω) e (x, t) ∈ Ω, allora il PC ha almeno una soluzione.

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

Sia f : Ω ⊂ R^nk × J → Rⁿ con J intervallo, f ∈ C⁰(Ω) e lipschitziana rispetto a (x, x⁰, . . . , x^(k−1)), uniformemente rispetto a t su I con J compatto, allora ∀(x₀, x₁, . . . , x_k−1, t) ∈ J il Problema di Cauchy ha 1! soluzione su J.

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

Sia f : Ω ⊂ R^nk+1 → Rⁿ con f ∈ C¹(Ω) e (x₀, x₁, . . . , x_k−1, t) ∈ Ω, allora ∃ U(t₀) t.c. il Problema di Cauchy ha 1! soluzione su U(t₀).

2 Metodi ad un passo

2.1 Metodi di Eulero

Lemma 2.1. (Lemma di Gronwall)

Sia 0 ≤ v(t) ≤ A + B ∫_t0^t v(s)ds con t ∈ (t₀, t₀ + T), A, B, T ≥ 0. Allora v(t) ≤ Ae^BT.

Nota. Se fosse v(t) = A + B ∫_t0^t v(s)ds avrei v⁰(t) = Bv(t), v(t₀) = A ⇒ v(t) = Ae^B(t−t₀).

Dimostrazione. Sia z(t) = A + B ∫_t0^t v(s)ds. Segue che:

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

z(t₀) = A

⇒ z⁰(t) − Bz(t) ≤ 0

⇒ z⁰(t)e^{−B(t−t₀)} − Bz(t)e^{−B(t−t₀)} ≤ 0

⇒ (z(t)e^−Bt)⁰ ≤ 0

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

z(t)e^{−B(t−t₀)} ≤ A

z(t) ≤ Ae^B(t−t₀)

v(t) ≤ z(t) ≤ Ae^B(t−t₀) ≤ Ae^BT.

Lemma 2.2. (Disuguaglianze e identità)

• 2ab ≤ ηa² + b²/η, ∀η > 0 poiché si ha η(a − b/η)² ≥ 0.
• a(a − b) = 1/2(a² + (a − b)² − b²).
• Se a ≤ b + c, allora a² ≤ b² + c² (disuguaglianza triangolare).

Lemma 2.3.

Da a|e_n+1| ≤ a|e_n| + b|σ_n| si ottiene in modo ricorsivo:

|e_n+1| ≤ |e₀|(a/a)ⁿ⁺¹ + b/a ∑_i=0ⁿ (a/a)ⁿ⁻ⁱ max(|σ_n|). (2.1.1)

Teorema 2.4.

Dato un PVI:

y⁰ = f(t, y)

y(t₀) = y₀

Se considero il PVI:

ỹ⁰ = f(t, ỹ) + σ(t)

ỹ(t₀) = y₀ + σ₀

Con f(x, y) lipschitziana rispetto a y, uniformemente rispetto a t, allora questo problema non è ben posto.

Dimostrazione. Se considero:

|y(t) − ỹ(t)| = |∫_t0^t (f(s, y(s)) − f(s, ỹ(s)) + σ(s))ds|

≤ ∫_t0^t |f(s, y(s)) − f(s, ỹ(s))|ds + ∫_t0^t |σ(s)|ds, con t ∈ (t₀, t₀ + T), T ∈ R.

= L ∫_t0^t |y(s) − ỹ(s)|ds + ∫_t0^t |σ(s)|ds.

Per il lemma di Gronwall:

|y(t) − ỹ(t)| ≤ (|σ₀| + ∫_t0^t |σ(s)|ds) · e^LT.

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:

y₀ · y₁ = f₀ · f₁ ≤ 0.

(y₀)² ≤ 0.

Teorema 2.6.

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

Dimostrazione. Dato:

δ̃(t) = ỹ⁰(t) − y⁰(t) = f(t, ỹ(t)) − f(t, y(t)) + σ(t)

δ(t₀) = 0

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

⇒ δ̃(t)δ(t) ≤ σ(t)δ(t)

⇒ (δ̃(t))² ≤ 2σ(t)δ(t).

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

1/2 δ²(t) − 1/2 δ²(t₀) ≤ ∫_t0^t σ(s)δ(s)ds.

δ²(t) ≤ 2 ∫_t0^t |σ(s)| · |δ(s)|ds ≤ 2 max_x∈(t0,t)|δ(s)| ∫_t0^t|σ(s)|ds.

A questo punto ho che 2 max_x∈(t0,t)|δ(s)| ∫_t0^t|σ(s)|ds è monotona crescente.

⇒ max_x∈(t0,t)|δ(t)|² ≤ 2 max_x∈(t0,t)|δ(s)| ∫_t0^t|σ(s)|ds.

Pongo a = max_x∈(t0,t)|δ(t)|, b = ∫_t0^t|σ(s)|ds e quindi ho a² ≤ 2ab.

Inoltre 2ab ≤ a²/2 + 2b².

2 max_x∈(t0,t)|δ(t)|² ≤ a² + 4b².

Cioè:

max_x∈(t0,t)|δ(t)|² ≤ 2 + 4 ∫_t0^t|σ(s)|ds.

E quindi:

max_x∈(t0,t)|δ(t)| ≤ 2 + 2 √ ∫_t0^t|σ(s)|ds. (2.1.2)

Definizione 2.7. (Metodo di Eulero esplicito)

Dato il PC:

y⁰ = f(t, y(t)), t ∈ [t₀, t₀ + T]

y(t₀) = y₀

Assegniamo ad ogni punto nel piano una velocità. Noi sappiamo che f(t₀, y₀) = y⁰(t₀), quindi, in U(t₀, y₀) se considero la retta tangente y₀ + f(t₀, y₀)(t − t₀) al posto della funzione si avrà un errore o(h) e l’idea è di riconsiderare il PC in t₀ + h = t₁, chiamiamo quindi Y₁ = y₀ + hf(t₀, y₀) e ricorsivamente Y_n = Y_n−1 + hf(t_n, Y_n−1). Questo metodo è chiamato metodo di Eulero esplicito.

Nota. Quando si utilizza il metodo di Eulero esplicito, come da qualsiasi altro metodo, si commettono due tipi di errori differenti.

• Il primo errore lo si commette muovendosi sulla retta tangente al posto che sulla funzione in un intorno della soluzione.
• 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))/h ' y⁰(t) = f(t, y(t)).

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

(Y_n+1 − Y_n)/h ' f(t_n, Y_n).

Dal quale si scrive subito Y_n+1 = Y_n + hf(t_n, Y_n).

Di seguito diamo le definizioni dei seguenti errori.

Definizione 2.8. (Errore locale di troncamento)

Si tratta del primo errore, il muoversi sulla retta tangente, e riscalando per h si ha:

δ_n(h) = 1/h(y(t_n+1) − [y(t_n) + hf(t_n, y(t_n))]).

Esso è quindi l’errore in una iterazione supponendo di avere la soluzione esatta, nel caso particolare usando Eulero esplicito.

Definizione 2.9. (Consistenza)

Diciamo che un metodo è consistente se δ_n(h) → 0 per h → 0, ∀n.

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

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

δ_n(h) = 1/h[y(t_n) + hy⁰(t_n) + h²/2 y⁰⁰(t_n) + o(h²) − y(t_n) − hy⁰(t_n)] = h/2 y⁰⁰(t_n) + o(h) = O(h).

E quindi δ_n(h) = O(h).

Definizione 2.10. (Errore globale di troncamento)

Si tratta dell’errore:

e_n = y(t_n) − Y_n.

Definizione 2.11. (Metodo convergente)

Diciamo che un metodo è un metodo convergente se e_n → 0, per h → 0.

Teorema 2.12.

Il metodo di Eulero esplicito è convergente.

Dimostrazione. Possiamo riscrivere il PC, o PVI, in forma discreta; si ha, scrivendo y_n al posto di y(t_n) e σ_n al posto di σ_n(h):

y_n+1 = y_n + hf(t_n, y_n) + hσ_n e Y_n+1 = Y_n + hf(t_n, Y_n).

Se poniamo e_n+1 = y_n+1 − Y_n+1 allora otteniamo:

e_n+1 = e_n + h[f(t_n, y_n) − f(t_n, Y_n)] + hσ_n

⇒ |e_n+1| ≤ |e_n| + h|L(y_n − Y_n)| + h|σ_n| = (1 + hL)|e_n| + h|σ_n|.

⇒ |e_n+1| ≤ (1 + hL)(1 + hL)|e_n−1| + h|σ_n| + h|σ_n| ≤ . . . ≤ (1 + hL)ⁿ⁺¹|e₀| + h ∑_j=0ⁿ(1 + hL)^n−j|σ_j|.

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

|e_n+1| ≤ hσ ∑_j=0ⁿ(1 + hL)^j = hσ ((1+hL)ⁿ⁺¹ − 1)/hL ≤ σ/L e^LT.

Poiché σ = Mh, M = max(y²), si ha che il metodo è convergente di ordine 1.

Definizione 2.13. (Problemi di Stabilità)

Diciamo che un metodo, quando esso viene applicato ad un problema dissipativo lineare, cioè nella forma y⁽ⁿ⁾ = λ ∑_i y⁽ⁱ⁾, ha un problema di stabilità o che è instabile, se scelto il passo finito h, dal punto di vista teorico esso è infinitesimale, ma dal punto di vista pratico, su un calcolatore, no, la soluzione Y_n non tende a 0 quando n tende a infinito ma vi è un errore che viene propagato al 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:

y⁰ = λy, Re(λ) < 0

y(0) = 1

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

Y_n+1 = Y_n + hλY_n = (1 + hλ)Y_n = . . . = (1 + hλ)ⁿY₀ = (1 + hλ)ⁿ.

Se h > 0 è fissato, si ha che se |1 + hλ| > 1 allora per n → infinito Y_n → infinito al contrario di y(t_n). Questo è un problema di stabilità, perché il metodo di Eulero esplicito non converge alla soluzione poiché h è troppo grande, nello specifico serve |1 + hλ| < 1 ⇒ −2 < Re(hλ) < 0 ⇒ h < −2/Re(λ).

Nota. Il caso descritto nell’esempio si può generalizzare facilmente; data y : R → R^k e l’equazione:

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

y(0) = y₀

Allora si deve scegliere il passo h nel seguente modo.

Si ha che ∃ X t.c. XAX⁻¹ = Δ diagonale con λ_i autov. di A.

A = X⁻¹ΔX ⇒ y⁰ = Ay ⇒ X⁻¹y⁰ = Δ(X⁻¹y).

Ponendo z = X⁻¹y si ha z⁰ = Δz, si ha quindi, su ogni singola componente, una equazione del tipo z⁰_i = f(z_i).

Possiamo quindi applicare il metodo di Eulero esplicito, e per evitare problemi di stabilità dovremmo scegliere h t.c. |1 + hλ_i| < 1, ∀i = 1 . . . k, cioè h < 2/|max|λ_i|.

Nel caso invece λ fosse a sua volta in funzione di t, sarebbe allora necessario scegliere h tale che h ≤ 2/max(λ(t)).

Definizione 2.16. (Regione di assoluta stabilità)

Si tratta di un sottoinsieme di C i cui elementi z = h · λ non causano problemi di stabilità al 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_n+1 − Y_n)/h ' f(t_n+1, Y_n+1)

Otteniamo il metodo iterattivo Y_n+1 = Y_n + hf(t_n+1, Y_n+1).

Che viene chiamato implicito perché la soluzione Y_n+1 dipende da Y_n+1 stessa. A questo punto si possono usare diversi metodi per risolvere il metodo implicito, come ad esempio porre:

Y_n+1 − Y_n − hf(t_n+1, Y_n+1) = g(Y_n+1).

Usando il metodo di bisezione, che però è ”lento”, oppure il metodo di Newton, notando che y⁰_n+1 − y⁰_n. Ma se si avesse y⁰ = 1 o se non fossimo abbastanza vicini al punto, il metodo di Newton potrebbe non convergere.

Un’altra strada è quella di porre Y_n+1 = Φ(Y_n+1) = Y_n + hf(t_n+1, Y_n+1), valutare quindi il problema come problema di punto fisso.

Si nota subito che Φ⁰(S) = hf_y(t_n+1, S) e che per la convergenza di questo metodo serve h < 1/f_y perché Φ sia una contrazione, condizione sufficiente, ma non necessaria.

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 è consistente, infatti si ha:

hσ_n(h) = y(t_n+1) − y(t_n) − h f(t_n+1, y(t_n+1)).

hσ_n(h) = y_n + hy⁰_n + h²/2 y⁰⁰_n + . . . − y_n − h(y⁰_n + hy⁰⁰_n + . . .) = y_n.

Convergenza

Si può dimostrare che il metodo di Eulero implicito converge con ordine 1.

Sia y(t_n+1) = y(t_n) + hf(t_n+1, y(t_n+1)) + hσ_n+1 e Y_n+1 = Y_n + hf(t_n+1, Y_n+1), posto allora e_n+1 = y(t_n+1) − Y_n+1 si ha:

e_n+1 = e_n + h[f(t_n+1, y_n+1) − f(t_n+1, Y_n+1)] + hσ_n+1.