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

Dispense del corso

Laboratorio di Metodi Numerici per le

Equazioni Differenziali

Dott. Marco Caliari

a.a. 2007/08

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 . . . . . . . . . . . . . . . . 8

0.3.1 Alcuni comandi per matrici sparse . . . . . . . . . . . . 9

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

0.4.1 Metodo di Newton modificato . . . . . . . . . . . . . . 10

0.5 Esponenziale di matrice . . . . . . . . . . . . . . . . . . . . . 11

0.5.1 Formula delle variazioni delle costanti . . . . . . . . . 11

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

0.6 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

I ODEs 17

1 Equazioni differenziali ordinarie 18

1.1 θ-metodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

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

1.2.1 Passo di integrazione variabile . . . . . . . . . . . . . . 21

1.2.2 Ordine dei metodi con passo variabile . . . . . . . . . . 21

1.3 Integratori esponenziali . . . . . . . . . . . . . . . . . . . . . . 21

1.4 Metodi di splitting esponenziale . . . . . . . . . . . . . . . . . 23

1.4.1 Strang’s splitting . . . . . . . . . . . . . . . . . . . . . 23

1.5 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3

4 INDICE

II BVPs 26

2 Problemi ai limiti 272.1 Differenze finite . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.1.1 Condizioni di Dirichlet omogenee . . . . . . . . . . . . 272.1.2 Condizioni di Dirichlet non omogenee . . . . . . . . . . 282.1.3 Condizioni di Neumann . . . . . . . . . . . . . . . . . 282.1.4 Norme . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.2 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

III PDEs 31

3 Equazioni iperboliche 323.1 Equazione del trasporto . . . . . . . . . . . . . . . . . . . . . 32

3.1.1 Metodi espliciti . . . . . . . . . . . . . . . . . . . . . . 333.1.2 Metodi impliciti . . . . . . . . . . . . . . . . . . . . . . 333.1.3 Estensione . . . . . . . . . . . . . . . . . . . . . . . . . 333.1.4 Condizioni al bordo . . . . . . . . . . . . . . . . . . . . 34

3.2 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

4 Equazioni paraboliche 354.1 Equazione di diffusione-convezione . . . . . . . . . . . . . . . . 35

4.1.1 Condizioni al bordo . . . . . . . . . . . . . . . . . . . . 354.1.2 Problema bidimensionale . . . . . . . . . . . . . . . . . 36

4.2 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

5 Metodo degli elementi finiti 385.1 Il problema di Poisson . . . . . . . . . . . . . . . . . . . . . . 385.2 Equazione del calore . . . . . . . . . . . . . . . . . . . . . . . 395.3 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

6 Metodi spettrali 416.1 Trasformata di Fourier discreta . . . . . . . . . . . . . . . . . 41

6.1.1 Note . . . . . . . . . . . . . . . . . . . . . . . . . . . . 426.2 Trasformata di Hermite . . . . . . . . . . . . . . . . . . . . . . 43

6.2.1 Note . . . . . . . . . . . . . . . . . . . . . . . . . . . . 446.3 Equazione di Schrodinger . . . . . . . . . . . . . . . . . . . . . 45

6.3.1 L’equazione di Schrodinger e la trasformata di Fourier . 456.3.2 L’equazione di Schrodinger e la trasformata di Hermite 46

6.4 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 − N , ove P e non singolare. Assegnato x(0), il termine x(l+1) ecalcolato ricorsivamente come

Px(l+1) = Nx(l) + b, l ≥ 0 (2)

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

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

ove B e chiamata matrice di iterazione.

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

0.2 Metodi di Richardson

Indicato con r(l) il residuo

r(l) = b− Ax(l) = Ax− Ax(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(0), x(l+1) e calcolato ricorsivamente come

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

ove α e un opportuno parametro di accelerazione. Dati x(0) e r(0) = b−Ax(0),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 .

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)

e dunque il minimo e dato dalla scelta

αl =z(l)Tr(l)

z(l)TAz(l).

0.2. METODI DI RICHARDSON 7

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

e r(0), 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(l)TAp(l−1) = 0 .

Per soddisfare questa proprieta e necessaria l’introduzione di un ulterioreparametro βl. Dati x(0), r(0), Pz(0) = r(0) e p(0) = z(0), 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(0), . . . , p(N−1) sono vettori linearmente indipendenti e non ce ne possonoessere piu di N .

Per questo motivo, tale metodo e detto semiiterativo.


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 · L(A) .

Una modifica consiste in∥

∥r(l)∥

∥ ≤ tol ·∥

∥r(0)∥

∥ , (7)

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

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.

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. MEMORIZZAZIONE DI MATRICI SPARSE 9

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.

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(0), 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)

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 .

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 (0), 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 11

function [x its errest] = newton(fun,fun1,x0,tol,maxits,varargin)

%

% [x its errest] = newton(fun,fun1,x0,tol,maxits,varargin)

%

its = 0;

errest = -feval(fun1,x0,varargin:)\feval(fun,x0,varargin:);

while (norm(errest) > tol) & (its <= maxits)

its = its+1;

x = x0+errest;

errest = -feval(fun1,x,varargin:)\feval(fun,x,varargin:);

x0 = x;

end

if (its > maxits)

warning(’Raggiunto il numero massimo di iterazioni.’)

end

Tabella 2: Metodo di Newton

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)

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.

0.5. ESPONENZIALE DI MATRICE 13

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 ≈ a0 + a1z + . . .+ apzp

b0 + b1z + . . .+ bqzq,

ove b0 = 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 = 1. Si ha allora

(

1 + z +z2

2+z3

6+ . . .

)

(b0 + b1z) = a0 + a1z

b0 + b1z + b0z + b1z2 + b0

z2

2+ o(z2) = a0 + a1z

da cui

b0 = 1

b0 = a0

b1 + b0 = a1

b1 +b02

= 0

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

exp(A) ≈ B = (b0I + b1A)−1(a0I + a1A) ,

cioe B e soluzione del sistema lineare (b0I + b1A)B = a0I + a1A.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

(z

2

)(z

4ϕ1

(z

2

)

+ 1)

.

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.

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

Tm , ove Vm ∈ R

N×m,V T

mVm = 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(z) = d0 + d1(z− ξ0) + d2(z− ξ1)(z− ξ2) + · . . . ·+dm(z− ξ1) · · · (z− ξm−1)

0.6. ESERCIZI 15

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

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

(

m−1∏

i=0

(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, ξm e il nodo per cui

m−1∏

i=0

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

m−1∏

i=0

|ξ − ξ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’inter-vallo [a, b] di modulo massimo. E chiaro che l’insieme dei primi m + 1 nodidi Leja coincide con l’unione di ξm con l’insieme dei primi m nodi di Leja.

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 ordineNA = 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 GNU


Octave 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(0) 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(0) =[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.2. METODI DI RUNGE–KUTTA EMBEDDED 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 = yk+1, ove

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

La matrice Jacobiana associata e

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

∂xj

.

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 di

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.

ordine 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. Indicata conyn+1

(p−1) l’approssimazione di y(tn+1) con il metodo Runge–Kutta di ordine p−1

e con yn+1(p) l’approssimazione di y(tn+1) con il metodo Runge–Kutta di ordine

p, si ha

∥

∥

∥yn+1(p) − yn+1

(p−1)

∥

∥

∥ =∥

∥

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

(p−1)

∥

∥

∥ / Ckp , (1.4)

20 CAPITOLO 1. EQUAZIONI DIFFERENZIALI ORDINARIE

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.

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.

ove k = tn+1 − tn e il passo di integrazione e Ckp e l’errore di troncamen-to locale del metodo di ordine p − 1. Infatti, si sta risolvendo il sistemadifferenziale

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

y(tn) = yn(p)

con entrambi i metodi. Se si vuole controllare tale errore con una tolleranzatol, 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.3. INTEGRATORI ESPONENZIALI 21

1.2.1 Passo di integrazione variabile

Indicato con kold l’attuale passo di integrazione e con knew il nuovo passod’integrazione, usando (1.4), si ha

knew =

tol∥

∥

∥yn+1(p) − yn+1

(p−1)

∥

∥

∥

1/p

· kold .

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

knew = min

2,max

0.6, 0.9 ·

tol∥

∥

∥yn+1(p) − yn+1

(p−1)

∥

∥

∥

1/p

· kold .

1.2.2 Ordine dei metodi con passo variabile

Supponiamo che l’integrazione nell’intervallo di tempo [t0, T ] sia condotta daun metodo di ordine p in n1 passi di lunghezza media k1 = (T − t0)/n1, conun errore finale err1 = Ckp

1 e in n2 passi di lunghezza media k2 = (T −t0)/n2,con un errore finale err2 = Ckp

2. 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.

1.3 Integratori esponenziali

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

un+1 = exp(kA)un + hϕ1(kA)b(un) (1.5)

ove un ≈ u(tn), tn+1 = tn + k.


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

‖u(tn+1) − un+1‖ = O(k2)

e

u(tn+1) = un+1, se b(u(t)) = b(u0) ≡ b .

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

Dimostrazione. Si ha

un+1 = exp(kA)un +

∫ tn+1

tn

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

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

u(tn+1) = exp(hA)un +

∫ tn+1

tn

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

= exp(kA)un +

∫ tn+1

tn

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

= un+1 +

∫ tn+1

tn

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

= un+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 non autonomo

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

u(t0) = u0

il metodo esponenziale—punto medio

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

e di ordine 2.

1.4. METODI DI SPLITTING ESPONENZIALE 23

Dimostrazione. Procedendo come sopra, si arriva a

u(tn+1) = un+1 +

∫ tn+1

tn

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

= un+1 +

∫ tn+1

tn

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

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

= un+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

u′(t) = Au(t) +Bu(t)

u(0) = u0

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)u0

si chiama Strang’s splitting e vale

exp(kA/2) exp(kB) exp(kA/2)u0 = exp(k(A+B))u0 + 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)

)

.


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)/A(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.

1.5. ESERCIZI 25

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

26

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. . . . . . . . . . . . 0

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

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

e1

h2

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

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

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

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

In GNU Octave e possibile generare la prima matrice con il comando toeplitz(sparse([0-1 zeros(1,N-2)]),sparse([0 1 zeros(1,N-2)])) e la seconda con il co-mando toeplitz(sparse([-2 1 zeros(1,N-2)])) (a meno dei fattori discala). I vettori [−1, 0, 1]/(2h) e [1,−2, 1]/h2 vengono chiamati stencils delladiscretizzazione.

27

28 CAPITOLO 2. PROBLEMI AI LIMITI

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. . . . . . . . . . . . 0

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

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

y1

y2.........

yN−2

+1

2h

−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 y0 = y(a) o yN−1 = y(b) ed ap-prossimare la derivata prima con uno stencil non simmetrico, eventualmentecon ordine di approssimazione minore. Per esempio, la discretizzazione delproblema ai limiti

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

y′(a) = y′ay(b) = yb

potrebbe essere

1

h2

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

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

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

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

y0

y1

y2......

yN−2

+1

h2

−h2y′a0......0yb

= 0

2.1.4 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)], conN 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 /√N , che risulta essere

una approssimazione mediante quadratura trapezoidale della norma in L2.

2.2. ESERCIZI 29

2.2 Esercizi

1. Si risolva il problema ai limiti

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

y(0) = 0

y(1) = 0

(2.2)

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 [1,−8, 0, 8,−1]/(12h)per la derivata prima e associata allo stencil [−1, 16,−30, 16,−1]/(12h2)per la derivata seconda.

3. Si risolva il problema ai limiti (2.2) 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.

30 CAPITOLO 2. PROBLEMI AI LIMITI


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

y(0) = 0

y(1) = 0

usando il metodo delle differenze finite del secondo ordine.

7. Si risolva il problema ai limiti (2.2) 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

31

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 − λ(hmj+1/2 − hm

j−1/2) (3.2)

con hmj+1/2 = h(ym

j , ymj+1), ove h(·, ·) 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 saraimlicito. La (3.2) puo essere interpretata come una discretizzazione in xj etm dell’equazione

d

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

che deriva da∂y

∂t(x, t) + c

∂y

∂x(x, t) = 0

integrando nello spazio tra xj−1/2 e xj+1/2 e ponendo Yj(t) =∫ xj+1/2

xj−1/2y(x, t)dx.

La scelta del flusso numerico caratterizza il metodo.

32

3.1. EQUAZIONE DEL TRASPORTO 33

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)):

hmj+1/2 =

1

2c(yn

j+1 + ynj )

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

hmj+1/2 =

1

2[c(yn

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

j+1 − ynj )]

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

hmj+1/2 =

1

2[c(yn

j+1 + ynj ) − λc2(yn

j+1 − ynj )]

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

hmj+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)):

hmj+1/2 =

1

2c(yn+1

j+1 + yn+1j )

3.1.3 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

∂u

∂t(x, t) +

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

∂x= 0 ,

ove hmj+1/2 = h(xj+1/2, y

mj , y

mj+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.

34 CAPITOLO 3. EQUAZIONI IPERBOLICHE

3.1.4 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 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.

3. 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

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.

4. 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) .

5. 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. Si chiama equazione di diffusione-convezione, in questo casocon condizioni al bordo miste Dirichlet–Neumann. Si puo allora discretizzaredifferenze 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 x0 = a e xN−1 = b, coninconite u0, u1, . . . , uN−1. La prima equazione del sistema di ODE riguardau0(t) e sara della forma u′0(t) = . . .. Siccome si vuole imporre una condizionedi Dirichlet costante, e sufficiente modifcare opportunamente la prima rigadella matrice A in modo da ottenere l’equazione u′0(t) = 0. Per quanto riguar-da la condizione di Neumann, si puo pensare di introdurre la variabile fittizia

35

36 CAPITOLO 4. EQUAZIONI PARABOLICHE

uN(t). Trascurando la dipendenza temporale, la discretizzazione dell’opera-tore di convezione per l’incognita uN−1 prima dell’imposizione delle condizio-ni al bordo e ux(xN−1) ≈ (uN−uN−2)/(2h) e la discretizzazione dell’operatoredi diffusione per la medesima e uxx(xN−1) ≈ ((uN − 2uN−1 + uN−2)/h

2. Unmodo per imporre la condizione di Neumann nulla e imporre che il valo-re di uN sia uguale al valore di uN−2. Le discretizzazioni diventano alloraux(xN−1) ≈ 0 e uxx(xN−1) ≈ ((2uN−1 + 2uN−2)/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 37

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, π] × [0, π] × (0,+∞)

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

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

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

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

uy(x, a2, 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.

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 φi+1

x1x0 xN

h1 hi+1 hN−1hi

xixixi−1 xi+1

Si ha allora

φi(x) =

x− xi−1

hi

, xi−1 ≤ x ≤ xi

xi+1 − x

hi+1

, xi ≤ x ≤ xi+1

0, altrimenti

e

d

dxφi(x) =

1

hi

, xi−1 ≤ x ≤ xi

−1

hi+1

, xi ≤ x ≤ xi+1

0, altrimenti

38

5.2. EQUAZIONE DEL CALORE 39

Dunque, per 0 < i < N ,

aii = a(φi, φi) =

∫ xi

xi−hi

(

1

hi

)2

dx+

∫ xi+hi+1

xi

( −1

hi+1

)2

dx =1

hi

+1

hi+1

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

∫ xi+hi+1

xi

−1

hi+1

· 1

hi+1

dx =−1

hi+1

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

a00 =

∫ x0+h1

x0

(−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 diquadratura del trapezio che risulta essere sufficientemente accurata. Si hadunque

fi = 〈f, φi〉 =

∫ xi

xi−1

f(x)x− xi−1

hi

dx+

∫ xi+1

xi

f(x)xi+1 − x

hi+1

dx ≈

≈ f(xi)

2hi +

f(xi)

2hi+1 =

f(xi)

2(hi + hi+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.

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

40 CAPITOLO 5. METODO DEGLI ELEMENTI FINITI

ove, per 0 < i < N ,

pii =

∫ xi+1

xi−1

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

3

pi i+1 =

∫ xi+1

xi

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

6

(5.1)

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

− 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.

2. Si risolva la PDE

∂

∂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.

Capitolo 6

Metodi spettrali

6.1 Trasformata di Fourier discreta

Sia u ∈ L2([a, b] → C) periodica (u(a) = u(b)) e N > 0 pari e fissato. Allora

u(x) =∞∑

m=−∞

umφFm(x) (6.1)

ove

φFm(x) =

ei(m−1−N/2)2π(x−a)/(b−a)

√b− a

.

E facile vedere che

〈φFj (x), φF

m(x)〉L2 =

∫ b

a

φFj (x)φF

m(x)dx = δjm

e dunque

um = 〈u, φFm〉L2 =

∫ b

a

u(x)e−i(m−1−N/2)2π(x−a)/(b−a)

√b− a

dx =

=√b− a

∫ 1

0

u((b− a)y + a)e−i(m−1−N/2)2πydy ≈

≈√b− a

N

N∑

n=1

(

u(xn)eiNπyn)

e−i(m−1)2πyn = um

ove x = (b − a)y + a, yn = (n − 1)/N e xn = (b − a)yn + a. Il vettoreN · [u1, u2, . . . , uN ]/

√b− a puo essere scritto come prodotto matrice-vettore

F [u(x1)eiNπy1 , u(x2)e

iNπy2 , . . . , u(xN)eiNπyN ]t, ove

F = (fmn), fmn = e−i(m−1)2πyn .

41

42 CAPITOLO 6. METODI SPETTRALI

Alternativamente, si puo usare la Fast Fourier Transform. Il comando diGNU Octave fft applicato al vettore [u(x1)e

iNπy1 , u(x2)eiNπy2 , . . . , u(xN)eiNπyN ]

produce il vettore N ·[u1, u2, . . . , uN ]/√b− a, cosı come il comando fftshift

applicato al risultato del comando fft applicato al vettore [u(x1), u(x2), . . . , u(xN)],come risulta chiaro dal codice riportato in Tabella 6.1.

clear all

N = 128;

a = -5;

b = 7;

y = linspace(0,1,N+1).’;

y = y(1:N);

x = (b-a)*y+a;

u = rand(N,1)+i*rand(N,1);

F = exp(-i*[0:N-1].’*2*pi*y.’);

uhat1 = F*(u.*exp(i*N*pi*y));

uhat2 = fft(u.*exp(i*N*pi*y));

uhat3 = fftshift(fft(u));

norm(uhat1-uhat2)

norm(uhat2-uhat3)

norm(uhat1-uhat3)

Tabella 6.1: Trasformata di Fourier discreta

Dunque

u(xn) ≈N∑

m=1

ukei(m−1−N/2)2π(xn−a)/(b−a)

√b− a

=N√b− a

1

N

(

N∑

m=1

ukei(m−1)2πyn

)

e−iNπyn = ûn .

Il comando ifft applicato al vettore [u1, u2, . . . , uN ] produce il vettore√b− a·

[û1eiNπy1 , û2e

iNπy2 , . . . , ûNeiNπyN ]/N , mentre, se applicato al risultato del

comando fftshift applicato al vettore [u1, u2, . . . , uN ], produce il vetto-re

√b− a · [û1, û2, . . . , ûN ]/N , come risulta chiaro dal codice riportato in

Tabella 6.2.

6.1.1 Note

La Fast Fourier Transform di un vettore di lunghezza N ha costo O(N logN),mentre il prodotto matrice-vettore O(N2). Tali costi sono pero asintotici e

6.2. TRASFORMATA DI HERMITE 43

clear all

N = 128;

a = -5;

b = 7;

y = linspace(0,1,N+1).’;

y = y(1:N);

x = (b-a)*y+a;

uhat = rand(N,1)+i*rand(N,1);

F = exp(-i*[0:N-1].’*2*pi*y.’);

IF = F’;

uhathat1 = (IF*uhat).*exp(-i*N*pi*y)/N;

uhathat2 = ifft(uhat).*exp(-i*N*pi*y);

uhathat3 = ifft(fftshift(uhat));

norm(uhathat1-uhathat2)



Tabella 6.2: Antitrasformata di Fourier discreta

nascondono i fattori costanti. Inoltre, GNU Octave puo far uso di implemen-tazioni ottimizzate di algebra lineare (come, ad esempio, le librerie ATLAS).In pratica, dunque, esiste un N0 sotto il quale conviene usare il prodottomatrice-vettore e sopra il quale la FFT.

La 6.1 si deve intendere

limM→+∞

∥

∥

∥

∥

∥

u(x) −M∑

m=−M

umφF(x)

∥

∥

∥

∥

∥

L2

= 0 .

L’identita di Parseval assicura che

‖u‖L2 =

(∫ b

a

|u(x)|2dx)1/2

≈ ‖u‖2 .

6.2 Trasformata di Hermite

Sia u ∈ L2(R → C). Allora

u(x) =∞∑

m=1

umφHm(x) (6.2)

oveφH

m(x) = Hm−1(x)e−x2/2


avendo indicato con Hm−1(x) il polinomio di Hermite normalizzato di gradom − 1 (Hm−1(x)e

−x2/2 e la corrispondente funzione di Hermite). Si puodimostrare che

∫

∞

−∞

φHj (x)φH

m(x)dx = δjm

e dunque

um = 〈u, φHm〉L2 =

∫

∞

−∞

u(x)Hm−1(x)e−x2/2dx ≈

≈N∑

n=1

u(xn)ex2n/2Hm−1(xn)wn = um .

essendo gli xn e i wn i nodi e i pesi di quadratura della formula di Gauss–Hermite.

Dunque

u(xn) ≈N∑

m=1

umHm−1(xn)e−x2n/2 = ˆun .

I polinomi di Hermite si costruiscono mediante la seguente relazione di ricor-renza

H0(x) =14√π, H1(x) =

√2x

4√π

Hm+1(x) =

√2xHm(x) −√

mHm−1(x)√m+ 1

, m > 0

6.2.1 Note

Non esiste la Fast Hermite Transform. L’unico modo per calcolare i coeffi-cienti e usare il prodotto matrice-vettore.

Vale la seguente:

− d2

dx2

(

Hm−1(x)e−x2/2

)

+x2Hm−1(x)e−x2/2 = (2m−1)Hm−1(x)e

−x2/2 . (6.3)

La 6.2 si deve intendere

limM→+∞

∥

∥

∥

∥

∥

u(x) −M∑

m=1

umφH(x)

∥

∥

∥

∥

∥

L2

= 0 .

L’identita di Parseval assicura che

‖u‖L2 =

(∫

∞

−∞

|u(x)|2dx)1/2

= ‖u‖2 .

6.3. EQUAZIONE DI SCHRODINGER 45

6.3 Equazione di Schrodinger

L’equazione di Schrodinger (in forma normalizzata) si scrive

i∂

∂tu(x, t) = −1

2

∂2

∂x2u(x, t) +

1

2x2u(x, t), (x, t) ∈ R × (0, T )

u(x, 0) = u0(x), ‖u0‖L2 = 1 x ∈ R

ove u(x, t) ∈ C (di solito l’incognita nell’equazione di Schrodinger si indicacon ψ). Si puo dimostrare che ‖u(·, t)‖L2 = 1, e cio dice, in qualche modo,che

limx→∞

u(x, t) = 0 .

Si parla in questo caso di condizioni al bordo asintotiche.

6.3.1 L’equazione di Schrodinger e la trasformata diFourier

Prima di tutto restringiamo il dominio a [a, b] ed eseguiamo lo splitting

i∂

∂tv(x, t) = −1

2

∂2

∂x2v(x, t)

i∂

∂tw(x, t) =

1

2x2w(x, t)

Nella prima equazione, decomponiamo

v(x, t) ≈N∑

m=1

vm(t)φFm(x) ,

moltiplichiamo per una funzione di base φFj (x) ed integriamo in [a, b]. Si

ottiene il sistema

i

v′1(t)v′2(t)

...v′N−1(t)v′N(t)

= −1

2

λF1 0 0 . . . 00 λF

2 0 . . . 0...

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

0 . . . 0 λFN−1 0

0 . . . 0 0 λFN

v1(t)v2(t)

...vN−1(t)vN(t)

con

λFm =

(

i(m− 1 −N/2)2π

b− a

)2

.


A questo punto, mediante l’antitrasformata di Fourier, possiamo approssima-re la soluzione v(x, t) nei punti xn. Risulta dunque possibile adesso risolvereil sistema

i

w1(t)′

w2(t)′

...wN−1(t)

′

wN(t)′

=1

2

x21 0 0 . . . 00 x2

2 0 . . . 0...

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

0 . . . 0 x2N−1 0

0 . . . 0 0 x2N

w1(t)w2(t)

...wN−1(t)wN(t)

ove wm(t) = w(xm, t).Notiamo che un potenziale V (x) diverso da quello armonico x2/2 non

comporterebbe nessuna modifica essenziale.Le condizioni al bordo non sono ben definite (sono periodiche se il poten-

ziale V (x) soddisfa V (a) = V (b)).

6.3.2 L’equazione di Schrodinger e la trasformata diHermite

Decomponiamo

u(x, t) ≈N∑

m=1

um(t)φHm(x) ,

moltiplichiamo per una funzione di base φHj (x) ed integriamo in R. Per la

proprieta (6.3) si ottiene il sistema

i

u′1(t)u′2(t)

...u′N−1(t)u′N(t)

=1

2

λH1 0 0 . . . 00 λH

2 0 . . . 0...

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

0 . . . 0 λHN−1 0

0 . . . 0 0 λHN

u1(t)u2(t)

...uN−1(t)uN(t)

conλH

m = 2m− 1 .

Un potenziale V (x) diverso da quello armonico forza lo splitting dell’equa-zione in

i∂

∂tv(x, t) = −1

2

∂2

∂x2v(x, t) +

1

2x2v(x, t)

i∂

∂tw(x, t) =

(

V (x) − 1

2x2

)

w(x, t)

In questo caso, le condizioni al bordo sono di tipo asintotico.

6.4. ESERCIZI 47

6.4 Esercizi

• Mostrare che se u0(x) = e−x2/2/ 4√π, allora la soluzione dell’equazione

di Schrodinger e u(x, t) = e−x2/2−it/2/ 4√π.

• Dimostrare che la soluzione dell’equazione di Schrodinger u(x, t) sod-disfa ‖u(·, t)‖L2 = 1 (Sugg.: calcolare d

dt

∫

R|u(x, t)|2dx).

• Risolvere l’equazione di Schrodinger con il metodo di splitting esponen-ziale, usando le differenze finite centrate del secondo ordine e condizionidi Dirichlet nulle per il primo problema. Mostrare la convergenza qua-dratica nel tempo e nello spazio, prendendo come condizione inizialeu0(x) = e−x2/2/ 4

√π.

• Risolvere l’equazione di Schrodinger con il metodo di splitting esponen-ziale e la trasformata di Fourier. Mostrare la convergenza quadraticanel tempo e spettrale nello spazio, prendendo come condizione inizialeu0(x) = e−x2/2/ 4

√π.

• Risolvere l’equazione di Schrodinger con il metodo di splitting esponen-ziale e la trasformata di Hermite, prendendo come potenziale V (x) =x2/3.

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

Documents

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