Dispense del corso di Laboratorio di Metodi Numerici per le...

Dispense del corso

Laboratorio di Metodi Numerici per le

Equazioni Differenziali

Dott. Marco Caliari

a.a. 2008/09

Questi appunti non hanno nessuna pretesa di completezza. Sono solo al-cune note ed esercizi che affiancano il corso di Metodi Numerici per le Equa-zioni Differenziali. Sono inoltre da considerarsi in perenne “under revision”e pertanto possono contenere discrepanze, inesattezze o errori.

Indice

0 Preliminari 5

0.1 Metodi iterativi per sistemi di equazioni lineari . . . . . . . . 5

0.2 Metodi di Richardson . . . . . . . . . . . . . . . . . . . . . . . 5

0.2.1 Metodo del gradiente precondizionato . . . . . . . . . . 6

0.2.2 Metodo del gradiente coniugato precondizionato . . . . 7

0.2.3 Test d’arresto . . . . . . . . . . . . . . . . . . . . . . . 8

0.3 Memorizzazione di matrici sparse . . . . . . . . . . . . . . . . 9

0.3.1 Alcuni comandi per matrici sparse . . . . . . . . . . . . 10

0.4 Metodo di Newton per sistemi di equazioni non lineari . . . . 10

0.4.1 Metodo di Newton modificato . . . . . . . . . . . . . . 11

0.5 Esponenziale di matrice . . . . . . . . . . . . . . . . . . . . . 12

0.5.1 Formula delle variazioni delle costanti . . . . . . . . . 12

0.5.2 Calcolo di exp(A) . . . . . . . . . . . . . . . . . . . . . 13

0.6 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

I ODEs 17

1 Equazioni differenziali ordinarie 18

1.1 theta-metodo . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.1.1 Caso lineare . . . . . . . . . . . . . . . . . . . . . . . . 19

1.1.2 Verifica della correttezza dell’implementazione . . . . . 19

1.2 Metodi di Runge–Kutta embedded . . . . . . . . . . . . . . . . 20

1.2.1 Passo di integrazione variabile . . . . . . . . . . . . . . 22

1.3 Integratori esponenziali . . . . . . . . . . . . . . . . . . . . . . 22

1.4 Metodi di splitting esponenziale . . . . . . . . . . . . . . . . . 24

1.4.1 Strang’s splitting . . . . . . . . . . . . . . . . . . . . . 24

1.5 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3

4 INDICE

II BVPs 27

2 Problemi ai limiti 282.1 Differenze finite . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.1.1 Condizioni di Dirichlet omogenee . . . . . . . . . . . . 282.1.2 Condizioni di Dirichlet non omogenee . . . . . . . . . . 292.1.3 Condizioni di Neumann . . . . . . . . . . . . . . . . . 29

2.2 Collocazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.2.1 Polinomi ortogonali . . . . . . . . . . . . . . . . . . . . 312.2.2 Collocazione negli zeri di polinomi ortogonali . . . . . . 33

2.3 Norme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.4 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

III PDEs 36

3 Equazioni iperboliche 373.1 Equazione del trasporto . . . . . . . . . . . . . . . . . . . . . 37

3.1.1 Metodi espliciti . . . . . . . . . . . . . . . . . . . . . . 383.1.2 Metodi impliciti . . . . . . . . . . . . . . . . . . . . . . 383.1.3 Verifica dell’ordine . . . . . . . . . . . . . . . . . . . . 383.1.4 Estensione . . . . . . . . . . . . . . . . . . . . . . . . . 393.1.5 Condizioni al bordo . . . . . . . . . . . . . . . . . . . . 39

3.2 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

4 Equazioni paraboliche 414.1 Equazione di diffusione-convezione . . . . . . . . . . . . . . . . 41

4.1.1 Condizioni al bordo . . . . . . . . . . . . . . . . . . . . 414.1.2 Problema bidimensionale . . . . . . . . . . . . . . . . . 42

4.2 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5 Metodo degli elementi finiti 445.1 Il problema di Poisson . . . . . . . . . . . . . . . . . . . . . . 445.2 Equazione del calore . . . . . . . . . . . . . . . . . . . . . . . 465.3 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

Capitolo 0

Preliminari

0.1 Metodi iterativi per sistemi di equazioni

lineari

I metodi iterativi per la soluzione del sistema lineare

Ax = b (1)

si basano sull’idea di calcolare la soluzione come limite di una successione divettori

x = liml→∞

x(l) .

Una strategia generale per costruire la successione x(l)l e basata sullo split-ting A = P − M , ove P e non singolare. Assegnato x(1), il termine x(l+1) ecalcolato ricorsivamente come

Px(l+1) = Mx(l) + b, l ≥ 1 (2)

Posto e(l) = x − x(l), si ha

e(l) = Bel−1, B = P−1M = I − P−1A ,

ove B e chiamata matrice di iterazione.

Lemma 1. Si ha liml→∞ e(l) = 0 per ogni e(1) se e solo se liml→∞ Bl = 0,cioe se e solo se ρ(B) < 1.

0.2 Metodi di Richardson

Indicato con r(l) il residuo

r(l) = b − Ax(l) = Ax − Ax(l) = A(x − x(l)) = Ae(l) ,

5

6 CAPITOLO 0. PRELIMINARI

il metodo (2) puo essere riscritto come

P (x(l+1) − x(l)) = r(l) . (3)

In questo contesto, P viene chiamata matrice di precondizionamento o pre-condizionatore di A e viene scelta in modo che la matrice di iterazioneB = I − P−1A abbia un raggio spettrale minore di 1 e la risoluzione di(3) sia “facile”.

Una generalizzazione dello schema (3) e il metodo di Richardson: datox(1), x(l+1) e calcolato ricorsivamente come

P (x(l+1) − x(l)) = αr(l) ,

ove α e un opportuno parametro di accelerazione. Dati x(1) e r(1) = b−Ax(1),l’algoritmo per calcolare x(l+1) e

Pz(l) = r(l)

x(l+1) = x(l) + αz(l)

r(l+1) = r(l) − αAz(l)

(4)

Il costo di un’iterazione e dato essenzialmente dalla risoluzione di un sistemalineare Pz(l) = r(l) facile e dal prodotto matrice-vettore Az(l). Tali metodirisulteranno particolarmente vantaggiosi per matrici sparse, in cui il numerodi elementi diversi da zero e O(N) piuttosto che O(N2) (e dunque il costodi un prodotto matrice-vettore e O(N)), se l’ordine della matrice e N .

Il calcolo del residuo r(l+1) = r(l) − αAz(l) (invece di r(l+1) = b − Ax(l+1))permette di ridurre la propagazione, attraverso il prodotto matrice-vettore,degli errori, in quanto il vettore z(l), contrariamente a x(l+1), diminuisce inmodulo al crescere di l.

0.2.1 Metodo del gradiente precondizionato

Siano A e P simmetriche e definite positive. Il metodo di Richardson puoessere generalizzato con una scelta dinamica del parametro di accelerazione,prendendo α = αl in modo tale che

∥

∥x(l+1) − x∥

∥

A, ‖y‖A =

√

yTAy

sia minima. Si ha

∥

∥x(l+1) − x∥

∥

2

A= (x(l) + αlz

(l) − x)TA(x(l) + αlz(l) − x) =

= α2l z

(l)TAz(l) − 2αlz(l)TA(x − x(l)) + (x(l) − x)TA(x(l) − x)

0.2. METODI DI RICHARDSON 7

e dunque il minimo e dato dalla scelta

αl =z(l)Tr(l)

z(l)TAz(l).

Il metodo ottenuto si chiama metodo del gradiente precondizionato. Dati x(1)

e r(1), l’algoritmo per calcolare x(l+1) e

Pz(l) = r(l)

αl =z(l)Tr(l)

z(l)TAz(l)

x(l+1) = x(l) + αlz(l)

r(l+1) = r(l) − αlAz(l)

(5)

Nel caso si scelga P = I, si ottiene il metodo del gradiente.

0.2.2 Metodo del gradiente coniugato precondizionato

Siano A e P simmetriche e definite positive. Il metodo del gradiente coniugatoprecondizionato e una generalizzazione del metodo di Richardson in cui

x(l+1) = x(l) + αlp(l)

ove i p(l)l sono coniugati, cioe soddisfano

p(i)TAp(j) = 0, i 6= j

Per soddisfare questa proprieta e necessaria l’introduzione di un ulterioreparametro βl. Dati x(1), r(1), Pz(1) = r(1) e p(1) = z(1), l’algoritmo percalcolare x(l+1) e

αl =z(l)Tr(l)

p(l)TAp(l)

x(l+1) = x(l) + αlp(l)

r(l+1) = r(l) − αlAp(l)

Pz(l+1) = r(l+1)

βl+1 =z(l+1)Tr(l+1)

z(l)Tr(l)

p(l+1) = z(l+1) + βl+1p(l)

(6)

Teorema 1. Il metodo del gradiente coniugato applicato ad una matrice diordine N converge in al piu N iterazioni (in aritmetica esatta).


Dimostrazione. La dimostrazione (omessa) si basa essenzialmente sul fattoche p(1), . . . , p(N) sono vettori linearmente indipendenti e non ce ne possonoessere piu di N .

Per questo motivo, tale metodo e detto semiiterativo.

Stima dell’errore

Vale la seguente stima dell’errore:

∥

∥e(l)∥

∥

A≤ 2

(

√

cond2(P−1A) − 1√

cond2(P−1A) + 1

)l−1∥

∥e(1)∥

∥

A

dalle quale si osserva che

• la stima d’errore decresce in ogni caso, poiche il numeratore e piupiccolo del denominatore;

• in particolare, nel caso P = I;

• tanto piu e piccolo il numero di condizionamento di P−1A, tanto piu ilmetodo ha convergenza veloce;

• nel caso limite di P = A, si ha∥

∥e(l)∥

∥

A≤ 0.

0.2.3 Test d’arresto

Un primo stimatore e costituito dal residuo: si arresta cioe il metodo iterativoquando

∥

∥r(l)∥

∥ ≤ tol · ‖b‖ .

Infatti, dalla precedente si ricava

∥

∥e(l)∥

∥

‖x‖ ≤ tol · cond(A) .

Una modifica consiste in

∥

∥r(l)∥

∥ ≤ tol ·∥

∥r(1)∥

∥ , (7)

che coincide con il precedente nel caso in cui come x(1) venga scelto il vettoredi zeri.

0.3. MEMORIZZAZIONE DI MATRICI SPARSE 9

0.3 Memorizzazione di matrici sparse

Sia A una matrice sparsa di ordine N con m elementi diversi da zero. Esistonomolti formati di memorizzazione di matrici sparse. Quello usato da GNUOctave e il Compressed Column Storage (CCS). Consiste di tre array: unprimo, data, di lunghezza m contenente gli elementi diversi da zero dellamatrice, ordinati prima per colonna e poi per riga; un secondo, ridx, dilunghezza m contenente gli indici di riga degli elementi di data; ed un terzo,cidx, di lunghezza N+1, il cui primo elemento e 0 e l’elemento i+1-esimo e ilnumero totale di elementi diversi da zero nelle prime i colonne della matrice.Per esempio, alla matrice

A =

1 0 0 00 2 3 04 0 5 60 0 0 7

corrispondono i vettori

data = [1, 4, 2, 3, 5, 6, 7]

ridx = [1, 3, 2, 2, 3, 3, 4]

cidx = [0, 2, 3, 5, 7]

Dato un vettore x, il prodotto matrice-vettore y = Ax e implementatodall’algoritmo in Tabella 1.

function y = ccsvect(data,ridx,cidx,x)

%

% y = ccsvect(data,ridx,cidx,x)

%

n = length(x);

y = zeros(size(x));

for j = 1:n

for i = cidx(j)+1:cidx(j+1)

y(ridx(i)) = y(ridx(i))+data(i)*x(j);

end

end

Tabella 1: Prodotto matrice-vettore in formato CCS.

In GNU Octave, il formato CCS e l’implementazione del prodotto matri-ce-vettore sono automaticamente usati dalla function sparse e dall’operato-re *, rispettivamente.


0.3.1 Alcuni comandi per matrici sparse

• Il comando speye(N) genera la matrice identita di ordine N .

• Il comando spdiags(v,0,N,N), ove v e un vettore colonna, genera lamatrice diagonale di ordine n avente v in diagonale. Se la dimensionedi v e minore di n, la diagonale viene riempita con zeri posti dopo ilvettore v. Se invece la dimensione di v e maggiore di N , vengono usatesolo le prime N componenti di v.

Sia V la matrice

V =

v11 v12 v13

v21 v22 v23...

......

vN1 vN2 vN3

Il comando spdiags(V,-1:1,N,N) genera la matrice

v12 v23 0 0 . . . 0v11 v22 v33 0 . . . 0

0 v21 v32. . . . . .

......

. . . . . . . . . . . ....

0 . . . 0 vN−2 1 vN−1 2 vN3

0 . . . . . . 0 vN−1 1 vN2

0.4 Metodo di Newton per sistemi di equa-

zioni non lineari

Consideriamo il sistema di equazioni non lineari

f1(x1, x2, . . . , xN) = 0

f2(x1, x2, . . . , xN) = 0

...

fN(x1, x2, . . . , xN) = 0

che puo essere riscritto, in forma compatta,

f(x) = 0 .

Dato x(1), il metodo di Newton per calcolare x(l+1) e

J (l)δx(l) = −f(x(l))

x(l+1) = x(l) + δx(l)(8)

0.4. METODO DI NEWTON PER SISTEMI DI EQUAZIONI NON LINEARI11

ove J (l) e la matrice Jacobiana, definita da

J(l)ij =

∂fi(x(l))

∂x(l)j

. (9)

Il criterio d’arresto solitamente usato e

∥

∥δx(l)∥

∥ ≤ tol .

function [x,iter,stimaerrore] = newton(sistema,jacobiano,x,tol,maxit,varargin)

%

% [x,iter,stimaerrore] = newton(sistema,jacobiano,x,tol,maxit,varargin)

%

iter = 0;

stimaerrore = -feval(jacobiano,x,varargin:)\feval(sistema,x,varargin:);

while (norm(stimaerrore) > tol) & (iter < maxit)

x = x+stimaerrore;

iter = iter+1;

stimaerrore = -feval(jacobiano,x,varargin:)\feval(sistema,x,varargin:);

end

stimaerrore = norm(stimaerrore);

if (stimaerrore > tol)

warning(’Impossibile raggiungere la tolleranza richiesta’)

warning(’entro il numero massimo di iterazioni consentito.’)

end

Tabella 2: Metodo di Newton

0.4.1 Metodo di Newton modificato

Il metodo di Newton (8) richiede il calcolo della matrice Jacobiana e la sua“inversione” ad ogni passo k. Questo potrebbe essere troppo oneroso. Unastrategia per ridurre il costo computazionale e usare sempre la stessa matriceJacobiana J (1), oppure aggiornarla solo dopo un certo numero di iterazioni.In tal modo, per esempio, e possibile usare la stessa fattorizzazione L(l)U (l)

per piu iterazioni successive.


0.5 Esponenziale di matrice

Data una matrice quadrata A ∈ RN×N , si definisce

exp(A) =∞

∑

j=0

Aj

j!.

Tale serie converge per qualunque matrice A, essendo A un operatore linearetra spazi di Banach e avendo la serie esponenziale raggio di convergenza ∞.Se A e B sono permutabili (cioe AB = BA), allora

exp(A + B) = exp(A) exp(B) .

0.5.1 Formula delle variazioni delle costanti

Data l’equazione differenziale

u′(t) = au(t) + b(u(t)), t > 0

u(t0) = u0

(10)

la soluzione puo essere scritta analiticamente mediante la formula delle va-riazioni delle costanti

u(t) = e(t−t0)au0 +

∫ t

t0

e(t−τ)ab(u(τ))dτ . (11)

Infatti, si ha

u′(t) = ae(t−t0)au0 + a

∫ t

t0

e(t−τ)ab(u(τ))dτ + e(t−t)ab(u(t)) = au(t) + b(u(t)) .

Si osservi che

∫ t

t0

e(t−τ)adτ = −1

a

∫ t

t0

−ae(t−τ)adτ = −1

ae(t−τ)a

∣

∣

∣

∣

t

t0

=

= −1

a

(

1 − e(t−t0)a)

= (t − t0)e(t−t0)a − 1

(t − t0)a=

= (t − t0)ϕ1((t − t0)a) ,

ove

ϕ1(z) =ez − 1

z=

∞∑

j=0

zj

(j + 1)!. (12)

0.5. ESPONENZIALE DI MATRICE 13

Consideriamo ora un sistema differenziale

u′(t) = Au(t) + b(u(t)), u ∈ RN×1, A ∈ R

N×N

u(t0) = u0

(13)

Ancora, la soluzione esplicita puo essere scritta come

u(t) = exp((t − t0)A)u0 +

∫ t

t0

exp((t − τ)A)b(u(τ))dτ .

Consideriamo l’approssimazione b(u(τ)) ≈ b(u(t0)) = b(u0). Allora

u(t) ≈ exp((t − t0)A)u0 +

∫ t

t0

exp((t − τ)A)b(u0)dτ =

= exp((t − t0)A)u0 + (t − t0)ϕ1((t − t0)A)b(u0) .

(14)

0.5.2 Calcolo di exp(A)

Come per la risoluzione di sistemi lineari, non esiste il modo per calcolareexp(A), ma diversi modi, ognuno adatto a particolari situazioni.

Matrici piene, di modeste dimensioni

Questi metodi si applicano, in pratica, a quelle matrici per le quali si usanoi metodi diretti per la risoluzione di sistemi lineari.

Decomposizione spettrale Se la matrice e diagonalizzabile, cioe A =V DV −1, allora exp(A) = V exp(D)V −1, ove exp(D) e la matrice diagonalecon elementi ed1 , ed2 , . . . , edN . Basta infatti osservare che

A2 = (V DV −1)2 = (V DV −1)(V DV −1) = V D2V −1

e scrivere exp(A) come serie di Taylor. La decomposizione spettrale di unamatrice costa, in generale, O(N3). Si ottiene in GNU Octave con il comandoeig.

Approssimazione razionale di Pade Si considera un’approssimazionerazionale della funzione esponenziale

ez ≈ a1zp−1 + a2z

p−2 + . . . + ap

b1zq−1 + b2zq−2 + . . . + bq

,


ove bq = 1 per convenzione. Essa e chiamata diagonale quando p = q. Sipuo dimostrare che le approssimazioni diagonali sono le piu efficienti. Fissatoil grado di approssimazione, si sviluppa in serie di Taylor la funzione espo-nenziale e si fanno coincidere quanti piu coefficienti possibile. Per esempio,fissiamo p = q = 2. Si ha allora

(

1 + z +z2

2+

z3

6+ . . .

)

(b1z + 1) = a1z + a2

b1z + 1 + b1z2 + z +

z2

2+ o(z2) = a1z + a2

da cui

1 = a2

b1 + 1 = a1

b1 +1

2= 0

L’approssimazione di Pade si estende banalmente al caso matriciale. Consi-derando sempre il caso p = q = 2, si ha

exp(A) ≈ B = (b1A + I)−1(a1A + a2I) ,

cioe B e soluzione del sistema lineare (b1A + I)B = a1A + a2I.L’approssimazione di Pade e accurata solo quando |z| < 1/2 (o, nel caso

matriciale, ‖A‖2 < 1/2). Per la funzione esponenziale esiste una tecnica,chiamata scaling and squaring che permette di aggirare il problema. Si usainfatti la proprieta

ez =(

ez/2)2

=(

ez/2j)2j

.

Se |z| > 1/2, allora |z|/2j < 1/2 per j > log2(|z|) + 1. Si calcola dunquel’approssimazione di Pade di ez/2j

e poi si eleva al quadrato j volte. Per lafunzione ϕ1 vale

ϕ1(z) =1

2(ez/2 + 1)ϕ1

(z

2

)

Anche l’approssimazione di Pade matriciale ha costo O(N3). In GNU Octavesi usa una variante di questa tecnica nel comando expm.

Matrici sparse, di grandi dimensioni

I metodi visti nel paragrafo precedente ignorano l’eventuale sparsita dellematrici. Inoltre, negli integratori esponenziali, non e mai richiesto di calcolareesplicitamente funzioni di matrice, ma solo funzioni di matrice applicate avettori, cioe exp(A)v (e l’analoga differenza tra calcolare A−1 e A−1v). Sipossono allora usare dei metodi iterativi.

0.5. ESPONENZIALE DI MATRICE 15

Metodo di Krylov Mediante la tecnica di Arnoldi e possibile, tramiteprodotti matrice-vettore, decomporre A in A ≈ VmHmV T

m , ove Vm ∈ Rn×m,

V Tm Vm = I, Vme1 = v e Hm e matrice di Hessenberg di ordine m (con m ≪

N). Allora AVm ≈ VmHm e

exp(A)v ≈ Vm exp(Hm)e1 .

Il calcolo di exp(Hm) e fatto mediante l’approssimazione di Pade. Il costodella tecnica di Arnoldi e O(Nm2) se A e matrice sparsa. E necessario inoltrememorizzare la matrice Vm.

Interpolazione su nodi di Leja Se il polinomio pm(z) interpola ez neinodi ξ0, ξ1, . . . , ξm, allora pm(A)v e una approssimazione di exp(A)v. E unabuona approssimazione se i nodi sono buoni (non equispaziati, per esempio)e se sono contenuti nell’involucro convesso dello spettro di A. E difficilestimare a priori il grado di interpolazione m necessario. E conveniente usarela formula di interpolazione di Newton

pm−1(z) = d1 +d2(z−ξ1)+d3(z−ξ1)(z−ξ2)+ · . . . ·+dm(z−ξ1) · · · (z−ξm−1)

ove dii sono le differenze divise. Tale formula si puo scrivere, nel casomatriciale,

pm−1(A)v = pm−2v + dmwm, wm =

(

m−1∏

i=1

(A − ξiI)

)

v = (A − ξm−1)wm−1

Dunque, la complessita e O(Nm) e richiesta la memorizzazione di un solovettore w.

Quali nodi usare? I nodi di Chebyshev, molto buoni per l’interpolazione,non possono essere usati, in quanto non permettono un uso efficiente dellaformula di interpolazione di Newton (cambiano tutti al cambiare del grado).I nodi di Leja sono distribuiti asintoticamente come i nodi di Chebyshev e,dati i primi m − 1, ξm e il nodo per cui

m−1∏

i=1

|ξm − ξi| = maxξ∈[a,b]

m−1∏

i=1

|ξ − ξi| ,

ove l’intervallo [a, b] e in relazione con lo spettro di A, per esempio [a, b] =σ(A) ∩ y = 0. Il primo nodo coincide, solitamente, con l’estremo dell’in-tervallo [a, b] di modulo massimo. E chiaro che l’insieme dei primi m nodidi Leja coincide con l’unione di ξm con l’insieme dei primi m − 1 nodi diLeja.


0.6 Esercizi

1. Implemetare le functions [data,ridx,cidx] = full2ccs(A) e [A] =

ccs2full(data,ridx,cidx).

2. Data una matrice A memorizzata in formato CCS e un vettore x,implementare la seguente formula: y = ATx.

3. Implementare il metodo del gradiente precondizionato (5) con una func-tion [x,iter] = pg(A,b,tol,maxit,P,x0), con test d’arresto dato da(7). Si testi il metodo, per esempio, sulle matrici di Toeplitz di ordine nA = toeplitz(sparse([4,1,zeros(1,N-2)])). Si testino come pre-condizionatori la matrice identita e la matrice formata dalla diagonaledi A (precondizionatore diagonale).

4. Si risolvano 6 sistemi lineari con le matrici di Hilbert di ordine N =4, 6, 8, 10, 12, 14 (hilb(N)) e termine noto scelto in modo che la solu-zione esatta sia il vettore [1, 1, . . . , 1]T usando il comando \ di GNUOctave, il metodo del gradiente precondizionato e il metodo del gra-diente coniugato precondizionato, implementato dalla function di GNUOctave pcg (si consulti il manuale relativo a questa function). Per que-sti ultimi due, si usi una tolleranza pari a 10−6, un numero massimo diiterazioni pari a 2000, il precondizionatore diagonale e un vettore ini-ziale x(1) di zeri. Si riporti, per ogni n, il numero di condizionamentodella matrice, l’errore in norma Euclidea rispetto alla soluzione esattae il numero di iterazioni dei metodi iterativi.

5. Risolvere il sistema non lineare

f1(x1, x2) = x21 + x2

2 = 1

f2(x1, x2) = sin(πx1/2) + x32

con il metodo di Newton (8). Si usi una tolleranza pari a 10−6, unnumero massimo di iterazioni pari a 150 e un vettore iniziale x(1) =[1, 1]T.

6. Si risolva lo stesso sistema non lineare usando sempre la matrice Jaco-biana relativa al primo passo e aggiornando la matrice Jacobiana ognir iterazioni, ove r e il piu piccolo numero di iterazioni che permettedi ottenere la soluzione con la tolleranza richiesta calcolando solo duevolte la matrice Jacobiana.

7. Si risolva lo stesso sistema non lineare usando la function fsolve diGNU Octave.

Parte I

ODEs

17

Capitolo 1

Equazioni differenziali ordinarie

1.1 θ-metodo

Consideriamo il sistema di ODEs

y1(t) = f1(y1(t), y2(t), . . . , yN(t), t)

y2(t) = f2(y1(t), y2(t), . . . , yN(t), t)

...

yN(t) = fN(y1(t), y2(t), . . . , yN(t), t)

con dato iniziale

y1(t0) = y10

y2(t0) = y20

...

yN(t0) = yN 0

che puo essere riscritto, in forma compatta,

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

y(t0) = y0

(1.1)

Notiamo come il sistema non autonomo (1.1) puo essere ricondotto ad unoautonomo

y(t) = f(y(t), yN+1(t))

yN+1(t) = 1

y(t0) = y0

yN+1(t0) = t0

(1.2)

ponendo yN+1(t) = t.

18

1.1. THETA-METODO 19

Il θ-metodo per il sistema (1.2) si scrive

yn+1 − yn

k= (1 − θ)f(yn) + θf(yn+1) (1.3)

ove tn+1 = tn +k e yn ≈ y(tn). Chiaramente, il θ metodo si riduce al metododi Eulero esplicito per θ = 0, al metodo di Eulero implicito per θ = 1 eal metodo di Crank–Nicolson per θ = 1/2. Nel caso implicito (θ 6= 0), adogni passo si deve risolvere un sistema di equazioni non lineari F (x) = 0,x = yn+1, ove

F (x) = x − kθf(x) − yn − k(1 − θ)f(yn).

La matrice Jacobiana associata e

Jij(x) = I − kθ∂fi(x)

∂xj

.

1.1.1 Caso lineare

Il caso piu frequente e quello lineare autonomo

y(t) = Ay(t)

y(t0) = y0

con passo di integrazione k costante. In tal caso, il metodo si scrive

(I − kθA)yn+1 = yn + k(1 − θ)Ayn

Nel caso implicito, si tratta dunque di risolvere un sistema lineare di matriceI − kθA ad ogni passo. Pertanto, per problemi di piccola dimensione, econveniente precalcolare la fattorizzazione LU della matrice. Altrimenti, sipuo considerare un metodo iterativo, ove si scelga come vettore iniziale peril calcolo di yn+1 la soluzione al passo precedente yn.

1.1.2 Verifica della correttezza dell’implementazione

Supponiamo di aver implementato un metodo per la soluzione del sistemadifferenziale

y(t) = f(y(t))

y(t0) = y0

e di volerne testare la corretta implementazione. L’idea e quella di creareuna soluzione artificiale x(t), inserirla nell’equazione e calcolarne il residuo

x(t) − f(x(t)) = g(t)

20 CAPITOLO 1. EQUAZIONI DIFFERENZIALI ORDINARIE

A questo punto, si risolve il sistema differenziale

y(t) = f(y(t)) + g(t)

y(t0) = x(t0)

fino ad un tempo T fissato, con due discretizzazioni di passo costante k1 =(T − t0)/n1 e k2 = (T − t0)/n2, rispettivamente. Si avranno errori finalierr1 = ‖yn1 − x(T )‖ = Ckp

1 e err2 = ‖yn2 − x(T )‖ = Ckp2. Si ha dunque

err2

err1

=

(

k2

k1

)p

,

da cui

log err2 − log err1 = p(log k2 − log k1) = −p(log n2 − log n1) .

Dunque, rappresentando in un grafico logaritmico-logaritmico l’errore in di-pendenza dal numero di passi, la pendenza della retta corrisponde all’ordinedel metodo, cambiato di segno. Tale verifica e valida anche nel caso di passinon costanti.

Falsa superconvergenza

Supponiamo che la soluzione esatta di un problema differenziale sia y = 1 eche, discretizzandolo con n = 1, 2, 4, 8 passi si ottengano le soluzioni y1 = 17,y2 = 9, y4 = 5 e y8 = 3. Supponiamo poi di averne calcolato una soluzionedi riferimento con n = 16 passi y = y16 = 2. In Figura 1.1 il grafico deglierrori che si ottiene.

1.2 Metodi di Runge–Kutta embedded

Supponiamo di avere un metodo di Runge–Kutta esplicito di ordine p − 1 ilcui tableau e riportato nella Tabella 1.1 e un altro metodo di Runge–Kutta diordine p il cui tableau e riportato nella Tabella 1.2. E chiaro che, dopo avercostruito il primo metodo, con una sola nuova valutazione della funzione f sipuo costruire il secondo metodo. Una tale coppia di metodi si dice embeddede si scrive di solito un unico tableau, come nella Tabella 1.3. Consideriamoil sistema differenziale

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

y(tn) = yn(p)

1.2. METODI DI RUNGE–KUTTA EMBEDDED 21

110

1 10errore

n

w.r.t riferimentow.r.t. esatta

Figura 1.1: Ordine dei metodi

0c2 a21

c3 a31 a32...

......

. . .

cs as1 as2 . . . as s−1

b1 b2 . . . bs−1 bs

Tabella 1.1: Metodo di Runge–Kutta di ordine p − 1.

0c2 a21

c3 a31 a32...

......

. . .

cs as1 as2 . . . as s−1

cs+1 as+11 as+12 . . . as+1 s−1 as+1 s

b1 b2 . . . bs−1 bs bs+1

Tabella 1.2: Metodo di Runge–Kutta di ordine p.

ove yn(p) e l’approssimazione di y(tn) ottenuta con il metodo Runge–Kutta di

ordine p. Si ha allora∥

∥

∥yn+1(p) − yn+1

(p−1)

∥

∥

∥ =∥

∥

∥yn+1(p) − y(tn+1) + y(tn+1) − yn+1

(p−1)

∥

∥

∥ / Ckpn+1 , (1.4)

ove kn+1 = tn+1 − tn e il passo di integrazione e Ckpn+1 e l’errore di tronca-

mento locale del metodo di ordine p − 1. Se si vuole controllare tale errore


0c2 a21

c3 a31 a32...

......

. . .

cs+1 as+11 as+12 . . . as+1 s

b1 b2 . . . bs

b1 b2 . . . bs bs+1

Tabella 1.3: Metodi di Runge–Kutta embedded di ordine p − 1 e p.

con una tolleranza tol, si puo allora imporre, ad ogni passo, che∥

∥

∥yn+1(p) − yn+1

(p−1)

∥

∥

∥ ≤ tol ,

rifiutando yn+1(p) nel caso la disuguaglianza non sia soddisfatta e calcolando un

nuovo passo di integrazione minore del precedente. Siccome p e l’ordine delmetodo piu accurato, ci si aspetta che l’errore globale sia controllato dallastessa tolleranza.

1.2.1 Passo di integrazione variabile

Per calcolare il successivo passo di integrazione kn+2, imponiamo che valga

Ckpn+2 = tol

e, ricavando 1/C da (1.4), si ha

kn+2 =

tol∥

∥

∥yn+1(p) − yn+1

(p−1)

∥

∥

∥

1/p

· kn+1 .

Per evitare che il passo di integrazione cambi troppo bruscamente, si puoadottare una correzione del tipo

kn+2 = min

2, max

0.6, 0.9 ·

tol∥

∥

∥yn+1(p) − yn+1

(p−1)

∥

∥

∥

1/p

· kn+1 .

1.3 Integratori esponenziali

Il metodo Eulero esponenziale per la risoluzione di (13) e

yn+1 = exp(kA)yn + kϕ1(kA)b(yn) (1.5)

1.3. INTEGRATORI ESPONENZIALI 23

ove yn ≈ y(tn), tn+1 = tn + k.

Proposizione 1. Per il metodo di Eulero esponenziale (1.5), se y(tn) = yn,si ha

‖y(tn+1) − yn+1‖ = O(k2)

e

y(tn+1) = yn+1, se b(y(t)) = b(y0) ≡ b .

Dunque il metodo e esatto per problemi lineari autonomi a coefficienti co-stanti, di ordine 1 altrimenti.

Dimostrazione. Si ha

yn+1 = exp(kA)yn +

∫ tn+1

tn

exp((tn+1 − τ)A)g(tn)dτ ,

ove si e posto g(t) = b(y(t)). Per la formula di variazioni delle costanti (11)

y(tn+1) = exp(hA)yn +

∫ tn+1

tn

exp((tn+1 − τ)A)g(τ)dτ =

= exp(kA)yn +

∫ tn+1

tn

exp((tn+1 − τ)A)(g(tn) + g′(τn)(τ − tn))dτ =

= yn+1 +

∫ tn+1

tn

exp((tn+1 − τ)A)g′(τn)(τ − tn))dτ =

= yn+1 + k2ϕ2(kA)g′(τn) ,

ove

ϕ2(z) =ez − 1 − z

z2=

∞∑

j=0

zj

(j + 2)!(1.6)

Proposizione 2. Per un problema lineare, non autonomo

y′(t) = Ay(t) + b(t), t > 0

y(t0) = y0

il metodo esponenziale—punto medio

yn+1 = exp(kA)yn + kϕ1(kA)b(tn + k/2)

e di ordine 2.


Dimostrazione. Procedendo come sopra, si arriva a

y(tn+1) = yn+1 +

∫ tn+1

tn

exp((tn+1 − τ)A)g′(τn + k/2)(τ − (tn + h/2))dτ =

= yn+1 +

∫ tn+1

tn

exp((tn+1 − τ)A)g′(τn + k/2)(τ − tn − k/2)dτ =

= yn+1 + (k2ϕ2(kA) − k/2kϕ1(kA))g′(τn + k/2) =

= yn+1 +

(

k2I

2+

k3A

6+ o(k3) − k2I

2− k3A

2+ o(k3)

)

g′(τn + k/2) .

1.4 Metodi di splitting esponenziale

Consideriamo il seguente problema differenziale lineare

y′(t) = Ay(t) + By(t)

y(0) = y0

Puo essere conveniente eseguire uno splitting e risolvere i due problemi

v′(t) = Av(t)

w′(t) = Bw(t)

e poi combinare opportunamente le due soluzioni.

1.4.1 Strang’s splitting

La seguente combinazione di tre stadi

exp(kA/2) exp(kB) exp(kA/2)y0

si chiama Strang’s splitting e vale

exp(kA/2) exp(kB) exp(kA/2)y0 = exp(k(A + B))y0 + O(k3) .

Per vederlo, basta riscrivere(

I +kA

2+

1

2

(

kA

2

)2

+ O(k3)

)

·(

I + kB +1

2(kB)2 + O(k3)

)

·

·(

I +kA

2+

1

2

(

kA

2

)2

+ O(k3)

)

=

=

(

I + k(A + B) +1

2(k(A + B))2 + O(k3)

)

.

1.5. ESERCIZI 25

Ovviamente si puo generalizzare al caso non lineare, combinando allo stessomodo le soluzioni ottenute con un qualunque metodo numerico; inoltre, sipossono scambiare i ruoli di A e B.

1.5 Esercizi

1. Si calcoli y(1), ove y(t) = Ay(t), y(0) = [1, . . . , 1]T, con A data da A =

100*toeplitz(sparse([-2,1,zeros(1,N-2)])), n = 10, usando il θ-metodo con θ = 0, 1/2, 1 e diversi passi temporali k = 2−3, 2−4, . . . , 2−8.Si confrontino i risultati con la soluzione di riferimento ottenuta usandoθ = 1/2 e k = 2−10, mettendo in evidenza l’ordine del metodo usato.Si provi anche il valore θ = 2/3, discutendo i risultati ottenuti.

2. Si risolva il sistema di ODEs

A(t) = −2a(t)A(t)

a(t) = A(t)2 + Ω(t)2 − a(t)2 − 1

Ω(t) = −2(a(t) + A(t))Ω(t)

(1.7)

con dato iniziale

A(0) = 0.5

a(0) = 2

Ω(0) = 10

fino ad un tempo finale T = 15, producendo un grafico della quantitaE(t) = (A(t)2 + a(t)2 + Ω(t)2 + 1)/(2A(t)).

3. Si implementi il seguente metodo di Runge–Kutta di ordine 2 (chiama-to metodo di Eulero Modificato) per un sistema non autonomo (1.1)con una function y = rk2(fun,y0,k,time), ove time e un vettore

012

12

0 1

bidimensionale contenente il tempo iniziale e quello finale.

4. Si implementi il seguente metodo di Runge–Kutta di ordine 3 per un si-stema non autonomo (1.1) con una function y = rk3(fun,y0,k,time),ove time e un vettore bidimensionale contenente il tempo iniziale equello finale.


012

12

1 -1 216

23

16

5. Si implementi il metodo RK23, Runge–Kutta di ordine 3 con metodoembedded Runge–Kutta di ordine 2, con passo variabile. Si risolva ilsistema differenziale (1.7), producendo il grafico della quantita E(t) finoad un tempo finale T = 15. Si metta in evidenza l’ordine del metodocon un grafico logaritmico-logaritmico.

6. Si implementi il metodo RK45 il cui tableau e riportato nella Tabella 1.4e lo si testi sul sistema differenziale (1.7).

014

14

38

332

932

1213

19322197

−72002197

72962197

1 439216

−8 3680513

− 8454104

12

− 827

2 −35442565

18594104

−1140

25216

0 14082565

21974104

−15

16135

0 665612825

2856156430

− 950

255

Tabella 1.4: Metodi di Runge–Kutta embedded di ordine 4 e 5.

Parte II

BVPs

27

Capitolo 2

Problemi ai limiti

2.1 Differenze finite

Consideriamo il seguente problema ai limiti

y′′(x) = f(y(x), y′(x), x), x ∈ (a, b) (2.1)

2.1.1 Condizioni di Dirichlet omogenee

Se vengono prescritti i valori al bordo y(a) = y(b) = 0, conviene usareuna discretizzazione interna. Si introducono cioe i nodi xi = a + hi, i =1, . . . , N − 2 ove h = (b − a)/(N − 1). Le matrici di discretizzazione alsecondo ordine della derivata prima e seconda diventano allora

1

2h

0 1 0 0 · · · 0−1 0 1 0 · · · 0

0. . . . . . . . . . . .

......

. . . . . . . . . . . . 00 · · · 0 −1 0 10 · · · 0 0 −1 0

e1

h2

−2 1 0 0 · · · 01 −2 1 0 · · · 0

0. . . . . . . . . . . .

......

. . . . . . . . . . . . 00 · · · 0 1 −2 10 · · · 0 0 1 −2

In GNU Octave e possibile generare la prima matrice con il comandotoeplitz(sparse([0,-1,zeros(1,N-2)]),sparse([0,1,zeros(1,N-2)]))

e la seconda con il comando toeplitz(sparse([-2,1,zeros(1,N-2)])) (ameno dei fattori di scala). I vettori [−1, 0, 1]/(2h) e [1,−2, 1]/h2 vengonochiamati stencils della discretizzazione.

28

2.1. DIFFERENZE FINITE 29

2.1.2 Condizioni di Dirichlet non omogenee

Se vengono prescritti i valori y(a) = ya 6= 0 o y(b) = yb 6= 0, conviene ancorausare una discretizzazione interna, ma si devono correggere i valori ai bordi.Per esempio, indicate con y1, y2, . . . , yN−2 le incognite, la discretizzazione alsecondo ordine della derivata prima diventa

y′(x1)y′(x2)

...

...

...y′(xN−2)

≈ 1

2h

0 1 0 0 · · · 0−1 0 1 0 · · · 0

0. . . . . . . . . . . .

......

. . . . . . . . . . . . 00 · · · 0 −1 0 10 · · · 0 0 −1 0

y1

y2.........

yN−2

+1

2h

−ya

0......0yb

Alternativamente, si possono introdurre anche i nodi al bordo. Si ha alloraxi = a + h(i − 1), i = 1, . . . , N . Poi si discretizza il problema senza tenerconto dei valori al bordo: per esempio, la discretizzazione del problema ailimiti

y′′(x) = 0, x ∈ (a, b)

y(a) = ya

y(b) = yb

diventa

1

h2

−2 1 0 0 · · · 01 −2 1 0 · · · 0

0. . . . . . . . . . . .

......

. . . . . . . . . . . . 00 · · · 0 1 −2 10 · · · 0 0 1 −2

y1

y2......

yN−1

yN

=

00......00

Infine, si correggono le equazioni relative ai nodi al bordo

1

h2

1 0 0 0 · · · 01 −2 1 0 · · · 0

0. . . . . . . . . . . .

......

. . . . . . . . . . . . 00 · · · 0 1 −2 10 · · · 0 0 0 1

y1

y2......

yN−1

yN

=1

h2

ya

0......0yb

2.1.3 Condizioni di Neumann

Se vengono prescritti i valori della derivata prima y′(a) = y′a o y′(b) = y′

b,e necessario introdurre almeno l’incognita in corrispondenza del nodo a o

30 CAPITOLO 2. PROBLEMI AI LIMITI

del nodo b ed approssimare la derivata prima con uno stencil non simme-trico, eventualmente con ordine di approssimazione minore. Per esempio, ladiscretizzazione del problema ai limiti

y′′(x) = 0, x ∈ (a, b)

y′(a) = y′a

y(b) = yb

potrebbe essere

1

h2

−h h 0 0 · · · 01 −2 1 0 · · · 0

0. . . . . . . . . . . .

......

. . . . . . . . . . . . 00 · · · 0 1 −2 10 · · · 0 0 1 −2

y1

y2

y3......

yN−1

+1

h2

−h2y′a

0......0yb

= 0

o, meglio,

1

h2

−32h 2h −1

2h 0 · · · 0

1 −2 1 0 · · · 0

0. . . . . . . . . . . .

......

. . . . . . . . . . . . 00 · · · 0 1 −2 10 · · · 0 0 1 −2

y1

y2

y3......

yN−1

+1

h2

−h2y′a

0......0yb

= 0

Infatti lo stencil all’indietro [−3, 4,−1]/(2h) produce una approssimazionedel secondo ordine in h della derivata prima.

2.2 Collocazione

Consideriamo l’approssimazione

y(x) ≈ y(x) =N

∑

j=1

yjφj(x)

per un problema differenziale (lineare)

D(y(x)) = b(x)

2.2. COLLOCAZIONE 31

ove φj(x)j e un insieme di funzioni ortonormali in un intervallo [a, b] ⊂R ∪ −∞, +∞

∫ b

a

φj(x)φi(x)dx = δji

Si ha alloraN

∑

j=1

yjD(φj(x)) = b(x) (2.2)

A questo punto, collochiamo l’equazione (2.2) in un insieme di nodi xiNi=1

per ottenereN

∑

j=1

yjD(φj(xi)) = b(xi)

ossiaAy = b, A = (aij) = (D(φj(xi))), bi = b(xi)

2.2.1 Polinomi ortogonali

Dato un intervallo [a, b] ed una funzione peso w(x) non negativa, e possibilecostruire una famiglia di polinomi ortonormali Lj(x)∞j=1, Lj(x) di gradoj − 1, che soddisfa

∫ b

a

Lj(x)Li(x)w(x)dx = δji

Tali polinomi soddisfano una relazione di ricorrenza a tre termini, del tipo

L1(x) = k dato

L2(x) = mx + q dato

Lj+1(x) = aj(x)Lj(x) + bj(x)Lj−1(x) j ≥ 2

Per alcuni intervalli e pesi, si conoscono esplicitamente i polinomi ortogo-nali (generalmente sono riportati in una forma non ortonormale). Vediamoqualche esempio (in forma ortonormale):

Polinomi di Legendre: sono ortonormali nell’intervallo [−1, 1] rispetto al-la funzione peso w(x) = 1. Soddisfano la relazione di ricorrenza

L1(x) =1√2

L2(x) =

√

3

2x

Lj+1(x) =

√

(2j + 1)(2j − 1)xLj(x) − (j − 1)√

2j+12j−3

Lj−1(x)

j, j ≥ 2


Polinomi di Chebyshev: sono ortonormali nell’intervallo [−1, 1] rispettoalla funzione peso w(x) = 1√

1−x2. Soddisfano la relazione di ricorrenza

T1(x) =1√π

T2(x) =

√

2

πx

T3(x) = 2xT2(x) −√

2T1(x)

Tj+1(x) = 2xTj(x) − Tj−1(x), j ≥ 3

Polinomi di Hermite: sono ortonormali nell’intervallo [−∞, +∞] rispettoalla funzione peso w(x) = e−x2

. Soddisfano la relazione di ricorrenza

H1(x) =14√

π

H2(x) =

√2x

4√

π

Hj+1(x) =

√2xHj(x) −√

j − 1Hj−1(x)√j

, j ≥ 2

Dalle relazioni di ricorrenza dei polinomi ortogonali e possibile ricavare lerelazioni di ricorrenza tra le derivate prime, seconde, etc.

Ci sono tre risultati fondamentali per ogni famiglia di polinomi ortonor-mali Lj(x)j:

1. gli zeri xjmj

m=1 del polinomio Lj+1(x) sono interni all’intervallo diortogonalita;

2. gli zeri del polinomio Lj+1(x) sono distinti dagli zeri di qualunque altropolinomio della stessa famiglia;

3. gli zeri del polinomio Lj+1(x) permettono, assieme agli opportuni pesidi quadratura wj

mjm=1, di calcolare esattamente integrali del tipo

∫ b

a

p2j−1(x)w(x)dx =

j∑

m=1

p2j−1(xm)wm

con p2j−1(x) polinomio di grado 2j − 1. In particolare,

∫ b

a

Li(x)2w(x)dx =

j∑

m=1

Li(xm)2wm = 1, i = 1, 2, . . . , j

2.3. NORME 33

2.2.2 Collocazione negli zeri di polinomi ortogonali

Data una famiglia di polinomi Lj(x)j ortonormali rispetto alla funzionepeso w(x), la famiglia di funzioni φj(x) = Lj(x)w(x)1/2j risulta essere

ortonormale. E di facile costruzione, se si conosce la relazione di ricorrenzatra i polinomi. E ha un insieme di nodi di collocazione candidato: se siconsidera un’approssimazione che coinvolge le prime N funzioni, gli zeri delpolinomio LN+1 non annullano le funzioni φj(x)N

j=1 (oltre che risultareottimali ai fini della quadratura).

Collocazione negli zeri dei polinomi di Legendre

La collocazione negli zeri dei polinomi di Legendre risulta particolarmentefacile da implementare. Come visto, si arriva a dover risolvere un sistemalineare del tipo

Ay = b, A = (aij) = (D(Lj(xi))), bi = b(xi)

L’unico problema e dato dalle condizioni al bordo: come imporle? Sarebbecomodo poter disporre dei nodi di collocazione x1 = a e xj = b. A quelpunto, la prima e l’ultima riga della matrice potrebbero essere modificatein modo da riprodurre le condizioni al bordo. Pertanto, si preferisce usare inodi di Legendre–Lobatto, definiti come gli zeri del polinomio (1 − x2)L′

j(x).

2.3 Norme

Data una funzione y(x) e due discretizzazioni su nodi equispaziati [y1, y2, . . . , yN ] ≈[y(x1), y(x2), . . . , y(xN)] e [y1, y2, . . . , yM ] ≈ [y(x1), y(x2), . . . , y(xM)], xii ⊂[a, b], xii ⊂ [a, b], con N 6= M , non ha molto senso confrontare gli errori‖[y(x1) − y1, y(x2) − y2, . . . , y(xN) − yN ]‖2 e ‖[y(x1) − y1, y(x2) − y2, . . . , y(xM) − yM ]‖2.Si preferisce usare la norma infinito, oppure la norma ‖[v1, v2, . . . , vN ]‖L2 =

‖v‖2

√

b−aN

, che risulta essere una approssimazione mediante quadratura con

formula dei rettangoli della norma in L2.

2.4 Esercizi

1. Si risolva il problema ai limiti

y′′(x) = y(x) + x, x ∈ (0, 1)

y(0) = 0

y(1) = 0

(2.3)


usando il metodo delle differenze finite del secondo ordine. Sapendoche la soluzione esatta e y(x) = (ex − e−x)/(e − e−1) − x, si mostriinoltre l’ordine del metodo mediante un grafico logaritmico-logaritmicodell’errore in norma infinito.

2. Si calcoli l’ordine di discretizzazione associata allo stencil centrato [1,−8, 0, 8,−1]/(12h)per la derivata prima e associata allo stencil centrato [−1, 16,−30, 16,−1]/(12h2)per la derivata seconda.

3. Si risolva il problema ai limiti (2.3) usando il metodo delle differenzefinite associate agli stencils dell’esercizio precedente, usando la corre-zione al secondo ordine per il primo e per l’ultimo nodo. Si mostriinoltre l’ordine del metodo mediante un grafico logaritmico-logaritmicodell’errore in norma infinito.


y′′(x) + y′(x) + y(x) − cos(x) = 0, x ∈ (0, π/2)

y(0) = 0

y(π/2) = 1

usando il metodo delle differenze finite del secondo ordine. Sapendo chela soluzione esatta e y(x) = sin(x), si mostri inoltre l’ordine del me-todo mediante un grafico logaritmico-logaritmico dell’errore in normainfinito.


y′′(x) + y′(x) + y(x) − cos(x) = 0, x ∈ (0, π/2)

y′(0) = 1

y(π/2) = 1

usando il metodo delle differenze finite del secondo ordine. Sapendo chela soluzione esatta e y(x) = sin(x), si mostri inoltre l’ordine del me-todo mediante un grafico logaritmico-logaritmico dell’errore in normainfinito.


y′′(x) = cos(y(x)), x ∈ (0, 1)

y(0) = 0

y(1) = 0

usando il metodo delle differenze finite del secondo ordine.

2.4. ESERCIZI 35

7. Si ricavi la relazione di ricorrenza dei polinomi ortonormali nell’inter-vallo [−∞,∞] rispetto alla funzione peso w(x) = e−α2x2

8. Noti gli zeri dei polinomi di Legendre e i pesi di quadratura della ri-spettiva formula gaussiana, si ricavino i nodi e i pesi di una formulagaussiana nell’intervallo [a, b] rispetto al peso w(x) = 1.

9. Si risolva il problema ai limiti (2.3) usando il metodo di collocazio-ne con polinomi di Legendre. Gli N nodi di collocazione in [a, b] e lavalutazione dei polinomi di Legendre e delle loro derivate sono dati dal-la function [L,x,L1,L2] = legendrepolynomials(N,a,b). Si mostriinoltre l’ordine del metodo mediante un grafico logaritmico-logaritmicodell’errore in norma infinito.

Parte III

PDEs

36

Capitolo 3

Equazioni iperboliche

3.1 Equazione del trasporto

Consideriamo il seguente problema iperbolico scalare

∂y

∂t(x, t) + c

∂y

∂x(x, t) = 0 (x, t) ∈ (a, b) × (0, +∞)

y(x, 0) = y0(x) x ∈ (a, b)

y(a, t) = ϕ(t) t ∈ (0, +∞)

(3.1)

ove c > 0 (dunque x = a e punto di inflow). Introduciamo una discretizza-zione spaziale a passo h, xj = a + jh, j = 1, . . . , N − 1, h = (b − a)/(N − 1)e temporale a passo k, tn = nk. Sia ora yn

j l’approssimazione della soluzioney(xj, tn) e λ = k/h. Un qualsiasi metodo alle differenze finite puo esserescritto nella forma

yn+1j = yn

j − λ(Fmj+1/2 − Fm

j−1/2) (3.2)

con Fmj+1/2 = F (ym

j , ymj+1), ove F (·, ·) e detto flusso numerico ed e, in generale,

un’approssimazione del flusso esatto cy(xj+1/2, tm) associato all’equazione(3.1). Se m = n, il metodo sara esplicito, se invece m = n + 1 il metodo saraimplicito. La (3.2) puo essere interpretata come una discretizzazione in xj etm dell’equazione

d

dtyj(t) +

1

hcy(xj+1/2, t) − cy(xj−1/2, t) = 0

che deriva da∂y

∂t(x, t) + c

∂y

∂x(x, t) = 0

37

38 CAPITOLO 3. EQUAZIONI IPERBOLICHE

integrando nello spazio tra xj−1/2 e xj+1/2 e ponendo

yj(t) =1

h

∫ xj+1/2

xj−1/2

y(x, t)dx

La scelta del flusso numerico caratterizza il metodo.

3.1.1 Metodi espliciti

Citiamo i seguenti metodi espliciti, indicando per ognuno l’ordine dell’erroredi troncamento:

Eulero in avanti/centrato (O(k + h2)):

Fmj+1/2 =

1

2c(yn

j+1 + ynj )

Lax–Friedrichs (O(h2/k + k + h)):

Fmj+1/2 =

1

2[c(yn

j+1 + ynj ) − λ−1(yn

j+1 − ynj )]

Lax–Wendroff (O(k2 + h2 + h2k)):

Fmj+1/2 =

1

2[c(yn

j+1 + ynj ) − λc2(yn

j+1 − ynj )]

upwind (o Eulero in avanti/decentrato) (O(k + h)):

Fmj+1/2 =

1

2[c(yn

j+1 + ynj ) − |c|(yn

j+1 − ynj )]

3.1.2 Metodi impliciti

Eulero all’indietro/centrato (O(k + h2)):

Fmj+1/2 =

1

2c(yn+1

j+1 + yn+1j )

3.1.3 Verifica dell’ordine

Per verificare l’ordine di un metodo di ordine p nel tempo e q nello spazio, epossibile scegliere un passo di discretizzazione temporale k proporzionale allapotenza q/p del passo di discretizzazione spaziale. In questo modo, l’errore

globale risulta essere O(kp + hq) = O(hqpp + hq) = O(hq).

3.2. ESERCIZI 39

3.1.4 Estensione

I metodi descritti sopra si estendono, senza modifiche, al caso c < 0, ove perole condizioni al bordo in (3.1) vanno imposte in x = b. Si estendono inoltrea problemi a coefficienti non costanti in forma conservativa

∂y

∂t(x, t) +

∂(c(x, t)y(x, t))

∂x= 0 ,

ove Fmj+1/2 = F (xj+1/2, y

mj , ym

j+1). Se c(x, t) > 0, le condizioni al bordo di

inflow sono in x = a, se c(x, t) < 0 le condizioni al bordo di inflow sono inx = b. Negli altri casi, essendo c(x, t) la velocita di propagazione, si possonooriginare shocks.

3.1.5 Condizioni al bordo

Nel punto di outflow non e possibile usare una discretizzazione della derivatache utilizzi punti esterni al dominio. Si possono allora usare discretizzazioninon simmetriche, eventualmente con ordine di approssimazione minore (vedi§ 2.1.3). Oppure, si possono estrapolare i valori mancanti.

3.2 Esercizi

1. Data l’approssimazione del secondo ordine della derivata y′(xN−1) ≈(yN − yN−2)/(2h), si ricavi per estrapolazione lineare il valore yN dayN−2 e yN−1 e lo si inserisca nell’approssimazione. Determinare l’ordinedell’approssimazione risultante.

2. Data l’approssimazione del secondo ordine della derivata seconda y′′(xN−1) ≈(yN −2yN−1 + yN−2)/(2h), si ricavi per estrapolazione quadratica il va-lore yN da yN−3, yN−2 e yN−1 e lo si inserisca nell’approssimazione.Determinare l’ordine dell’approssimazione risultante.

3. Data l’approssimazione del secondo ordine della derivata y′(xN−1) ≈(yN −yN−2)/(2h), si ricavi per estrapolazione quadratica il valore yN dayN−3, yN−2 e yN−1 e lo si inserisca nell’approssimazione. Si verifichi chel’approssimazione risultante e y′(xN−1) ≈ (3yN−1−4yN−2+yN−3)/(2h).Determinare l’ordine dell’approssimazione risultante.

4. Data l’approssimazione del secondo ordine della derivata seconda y′′(xN−1) ≈(yN −2yN−1 +yN−2)/h

2, si ricavi per estrapolazione cubica il valore yN

40 CAPITOLO 3. EQUAZIONI IPERBOLICHE

da yN−4, yN−3, yN−2 e yN−1 e lo si inserisca nell’approssimazione. Si ve-rifichi che l’approssimazione risultante e y′′(xN−1) ≈ (2yN−1 − 5yN−2 +4yN−3 − yN−4)/h

2.

5. Dato il problema di trasporto

∂y

∂t(x, t) + c

∂y

∂x(x, t) = 0 (x, t) ∈ (a, b) × (0, +∞)

y(x, 0) = y0(x) x ∈ (a, b)

y(a, t) = 0 t ∈ (0, +∞)

(3.3)

con c > 0 e y0(x) = 0 per x ≤ a, si dimostri che la soluzione e

y(x, t) = y0(x − ct) .

6. Si calcoli numericamente la soluzione dell’esercizio precedente al tempot = 1, con a = 0, b = π e y0(x) = 0 per x ≤ 0 e y0(x) = sin(x)50 per0 < x ≤ π usando i metodi descritti in § 3.1.1 e in § 3.1.2. Si mostrinogli ordini dei metodi.

Capitolo 4

Equazioni paraboliche

4.1 Equazione di diffusione-convezione

Consideriamo la seguente PDE

∂

∂tu(x, t) + c

∂

∂xu(x, t) = d

∂2

∂x2u(x, t), (x, t) ∈ (a, b) × (0, +∞)

u(x, 0) = u0(x) x ∈ (a, b)

u(a, t) = 0 t ∈ (0, +∞)

ux(b, t) = 0 t ∈ (0, +∞)

con c, d > 0 (in generale, c puo assumere qualunque segno, ma d deve esserepositivo). Si chiama equazione di diffusione-convezione, in questo caso concondizioni al bordo miste Dirichlet–Neumann. Si puo allora discretizzare condifferenze finite nello spazio, ottenendo un sistema di ODE

u′(t) = Au(t)

u(0) = u0

(4.1)

e integrare quest’ultimo nel tempo (metodo delle linee). Detto h il passo didiscretizzazione spaziale, bisogna assicurarsi che il numero di Peclet di griglia,definito come h|c|/(2d) sia minore di 1, per evitare problemi di stabilita.

4.1.1 Condizioni al bordo

Consideriamo una discretizzazione a N punti, con x1 = a e xN = b, con inco-gnite u1, u2, . . . , uN . La prima equazione del sistema di ODE riguarda u1(t)e sara della forma u′

1(t) = . . .. Siccome si vuole imporre una condizione diDirichlet costante, e sufficiente modifcare opportunamente la prima riga della

41

42 CAPITOLO 4. EQUAZIONI PARABOLICHE

matrice A in modo da ottenere l’equazione u′1(t) = 0. Per quanto riguarda

la condizione di Neumann, si puo pensare di introdurre la variabile fittiziauN+1(t). Trascurando la dipendenza temporale, la discretizzazione dell’ope-ratore di convezione per l’incognita uN prima dell’imposizione delle condizionial bordo e ux(xN) ≈ (uN+1 − uN−1)/(2h) e la discretizzazione dell’operatoredi diffusione per la medesima e uxx(xN) ≈ ((uN+1 − 2uN + uN−1)/h

2. Unmodo per imporre la condizione di Neumann nulla e imporre che il valoredi uN+1 sia uguale al valore di uN−1. Le discretizzazioni diventano alloraux(xN) ≈ 0 e uxx(xN) ≈ ((−2uN + 2uN−1)/h

2.

4.1.2 Problema bidimensionale

Consideriamo la PDE

∂

∂tu(x, y, t)+c1

∂

∂xu(x, y, t) + c2

∂

∂yu(x, y, t) =

= d1∂2

∂x2u(x, y, t) + d2

∂2

∂y2u(x, y, t),

(x, y, t) ∈ (a1, b1) × (a2, b2) × (0, +∞)

u(x, y, 0) = u0(x, y) (x, y) ∈ [a1, b1] × [a2, b2]

u(a1, y, t) = 0 t ∈ (0, +∞)

u(x, a2, t) = 0 t ∈ (0, +∞)

ux(b1, y, t) = 0 t ∈ (0, +∞)

uy(x, b2, t) = 0 t ∈ (0, +∞)

con c, d > 0. Se indichiamo con A1 e A2 le matrici di discretizzazione deglioperatori di derivazione in x e in y rispettivamente (con le condizioni albordo), allora il sistema differenziale che ne risulta sara

u′(t) = (I ⊗ A1 + A2 ⊗ I)u(t)

u(0) = u0

Le incognite u1(t), u2(t), . . . si intendono numerate da sinistra a destra e dalbasso all’alto.

4.2 Esercizi

1. Si risolva la PDE

∂

∂tu(x, t) =

∂2

∂x2u(x, t) + (1 + π2)et sin(πx), (x, t) ∈ (0, 1) × (0, +∞)

u(x, 0) = sin(πx) x ∈ (0, 1)

u(0, t) = u(1, t) = 0 t ∈ (0, +∞)

4.2. ESERCIZI 43

usando il metodo delle linee. Si confrontino gli integratori esponenzialiEulero esponenziale e esponenziale—punto medio e il metodo Euleroimplicito per l’integrazione fino al tempo T = 1. Si mostrino gli ordinidegli integratori temporali, sapendo che la soluzione esatta e u(x, t) =et sin(πx).

2. Si risolva l’equazione di diffusione-convezione

∂

∂tu(x, y, t)+c1

∂

∂xu(x, y, t) + c2

∂

∂yu(x, y, t) =

= d1∂2

∂x2u(x, y, t) + d2

∂2

∂y2u(x, y, t),

(x, y, t) ∈ (0, π)2 × (0, +∞)

u(x, y, 0) = (sin(x) sin(y))6 (x, y) ∈ [0, π]2

u(0, y, t) = 0 t ∈ (0, +∞)

u(π, y, t) = 0 t ∈ (0, +∞)

u(x, 0, t) = 0 t ∈ (0, +∞)

u(x, π, t) = 0 t ∈ (0, +∞)

con d1 = 5, d2 = 1, c1 = −1 e c2 = −10. Si usino differenze finitecentrate del secondo ordine nello spazio e il θ-metodo, RK23 e Euleroesponenziale nel tempo.

3. Si risolva l’equazione di diffusione-convezione dell’esercizio precedente,con le condizioni al bordo

ux(0, y, t) = 0 t ∈ (0, +∞)

u(π, y, t) = 0 t ∈ (0, +∞)

uy(x, 0, t) = 0 t ∈ (0, +∞)

u(x, π, t) = 0 t ∈ (0, +∞)

Capitolo 5

Metodo degli elementi finiti

5.1 Il problema di Poisson

Vogliamo risolvere il seguente problema differenziale

− d2

dx2u(x) = f(x), x ∈ (0, 1)

u(0) = u(1) = 0

con il metodo degli elementi finiti (lineari). Introduciamo una discretizzazio-ne dell’intervallo [0, 1] a passo variabile, come in figura.

φi+1

xN

h1 hN−1

xixixi−1 xi+1x1 x2

φiφ1

hi−1 hi

Si ha allora, per 1 < i < N ,

φi(x) =

x − xi−1

hi−1

, xi−1 ≤ x ≤ xi

xi+1 − x

hi

, xi ≤ x ≤ xi+1

0, altrimenti

e

d

dxφi(x) =

1

hi−1

, xi−1 ≤ x ≤ xi

− 1

hi

, xi ≤ x ≤ xi+1

0, altrimenti

44

5.1. IL PROBLEMA DI POISSON 45

mentre, per i = 1,

φ1(x) =

x2 − x

h1

, x1 ≤ x ≤ x2

0, altrimenti

e

d

dxφ1(x) =

− 1

h1

, x1 ≤ x ≤ x2

0, altrimenti

e, per i = N ,

φN(x) =

x − xN−1

hN−1

, xN−1 ≤ x ≤ xN

0, altrimenti

e

d

dxφN(x) =

1

hN−1

, xN−1 ≤ x ≤ xN

0, altrimenti

Dunque, per 1 < i < N ,

aii = a(φi, φi) =

∫ xi

xi−hi−1

(

1

hi−1

)2

dx +

∫ xi+hi

xi

(

− 1

hi

)2

dx =1

hi−1

+1

hi

ai i+1 = ai+1 i = a(φi, φi+1) =

∫ xi+hi

xi

− 1

hi

· 1

hi

dx = − 1

hi

Per i = 1 e i = N , si ha invece

a11 =

∫ x1+h1

x1

(

− 1

h1

)2

dx =1

h1

aNN =

∫ xN

xN−hN−1

(

− 1

hN−1

)2

dx =1

hN−1

Per quanto riguarda il calcolo di 〈f, φi〉 si puo ricorrere alla formula di quadra-tura del trapezio che risulta essere sufficientemente accurata. Si ha dunque,per 1 < i < N ,

fi = 〈f, φi〉 =

∫ xi

xi−1

f(x)x − xi−1

hi−1

dx +

∫ xi+1

xi

f(x)xi+1 − x

hi

dx ≈

≈ f(xi)

2hi−1 +

f(xi)

2hi =

f(xi)

2(hi−1 + hi) .

46 CAPITOLO 5. METODO DEGLI ELEMENTI FINITI

Per i = 1 e i = N si ha invece

f1 = 〈f, φ1〉 =

∫ x2

x1

f(x)x2 − x

h1

dx ≈ f(x1)

2h1

fN = 〈f, φN〉 =

∫ xN

xN−1

f(x)x − xN−1

hN−1

dx ≈ f(xN)

2hN−1

A questo punto si risolve il sistema lineare Au = f , dopo aver opportunamen-te modificato la matrice e il termine noto per imporre le condizioni al bordodi Dirichlet. Nel caso di condizioni di Neumann omogenee, non e necessariaalcuna modifica. Infatti, la forma debole del problema di Poisson e

∫ 1

0

du

dx(x)

d

dxφi(x)dx − du

dx(x)φi(x)

∣

∣

∣

1

0=

∫ 1

0

f(x)φi(x)dx, i = 1, . . . , N

Il secondo termine e naturalmente zero per 1 < i < N . Non considerarloneanche per i = 1 e i = N (cioe porlo a zero) significa richiedere du

dx(x1) =

dudx

(xN) = 0, visto che φ1(x1) = φN(xN) = 1. Nel caso di condizioni diNeumann non omogenee (per esempio in 0), basta modificare la prima rigadel sistema secondo l’equazione

∫ 1

0

du

dx(x)

d

dxφ1(x)dx = −du

dx(0) +

∫ 1

0

f(x)φ1(x)dx

5.2 Equazione del calore

La discretizzazione mediante elementi finiti della PDE

∂

∂tu(x, t) =

∂2

∂x2u(x, t), (x, t) ∈ (0, 1) × (0, +∞)

u(x, 0) = u0(x) x ∈ (0, 1)

u(0, t) = u(1, t) = 0 t ∈ (0, +∞)

porta al sistema di ODEs

Pu′(t) = Au(t)

u(0) = u0

ove, per 1 < i < N ,

pii =

∫ xi+1

xi−1

φi(x)φi(x)dx =hi−1 + hi

3

pi i+1 =

∫ xi+1

xi

φi(x)φi+1(x)dx =hi

6

(5.1)

5.3. ESERCIZI 47

mentre, per i = 1 e i = N ,

p11 =

∫ x2

x1

φ1(x)φ1(x)dx =h1

3

pNN =

∫ xN

xN−1

φN(x)φN(x)dx =hN−1

3

(5.2)

Le condizioni di Dirichlet omogenee per un nodo xi si impongono sostituendola riga corrispondente di P con zeri e 1 in diagonale e la riga corrispondentedi A con zeri.

La matrice P si chiama matrice di massa. Anche usando un metodoesplicito per la risoluzione del sistema differenziale, e necessaria l’inversionedella matrice di massa. Per tale motivo, si puo ricorrere alla tecnica del masslumping che consiste nel rendere diagonale la matrice P sostituendo ogni suariga con una riga di zeri e la somma degli elementi originali in diagonale. Talemodifica e equivalente all’approssimazione degli integrali in (5.1) mediante laformula del trapezio e dunque non riduce l’accuratezza del metodo. Infatti,la matrice P

(−1)L A (PL la matrice di massa con lumping) risulta uguale alla

matrice che si ottiene discretizzando con differenze finite centrate del secondoordine.

5.3 Esercizi

1. Si risolva l’equazione differenziale

− d2

dx2u(x) = ex, x ∈ (0, 1)

u(0) = 0

u(1) = 0

mediante il metodo degli elementi finiti e si confronti la soluzioneottenuta con quella esatta.


− d2

dx2u(x) = ex, x ∈ (0, 1)

u′(0) = 0

u(1) = 0


48 CAPITOLO 5. METODO DEGLI ELEMENTI FINITI


− d2

dx2u(x) = ex, x ∈ (0, 1)

u′(0) = 1

u(1) = 0



− d

dx

(

(1 + x)d

dxu(x)

)

= 1, x ∈ (0, 1)

u(0) = u(1) = 0

usando il metodo degli elementi finiti lineari, con griglia non equispa-ziata. Si usi la formula dei trapezi per il calcolo degli integrali.


∂

∂tu(x, t) =

∂2

∂x2u(x, t), (x, t) ∈ (0, 1) × (0, +∞)

u(x, 0) = sin(2πx) x ∈ (0, 1)

u(0, t) = u(1, t) = 0 t ∈ (0, +∞)

usando il metodo degli elementi finiti lineari, con griglia non equispa-ziata. Una volta ottenuto il sistema di ODE, lo si integri fino al tempot = 1 applicando uno dei metodi conosciuti per la sua risoluzione. Sitesti anche la tecnica del mass lumping.


∂

∂tu(x, t) =

∂2

∂x2u(x, t), (x, t) ∈ (0, 1) × (0, +∞)

u(x, 0) = sin(2πx) x ∈ (0, 1)

u(0, t) = 0 t ∈ (0, +∞)

ux(1, t) = 2π t ∈ (0, +∞)

usando il metodo degli elementi finiti lineari, con griglia non equispa-ziata. Una volta ottenuto il sistema di ODE, lo si integri fino al tempot = 1 applicando il metodo di Eulero esponenziale.

Dispense del corso di Laboratorio di Metodi Numerici per le...

Documents

Transcript of Dispense del corso di Laboratorio di Metodi Numerici per le...