Capitolo 1 Equazioni alle Derivate...

Capitolo 1

Equazioni alle Derivate Parziali

1.1 Introduzione

Un’Equazione alle Derivate Parziali (in breve PDE, acronimo dei terminiinglesi Partial Differential Equation) e un’equazione che mette in relazioneuna funzione incognita dipendente da due (o piu) variabili indipendenti allesue derivate parziali rispetto a queste variabili. La necessita di utilizzaretali equazioni sta nel fatto che queste consentono di descrivere in modo piuaccurato determinati fenomeni naturali. Infatti quando si cerca di descriverefenomeni dipendenti da diverse variabili indipendenti (piu comunemente po-sizione e tempo) allora e necessario utilizzare un modello differenziale allederivate parziali. Un esempio di PDE e il seguente

a∂2u

∂x2+ 2b

∂2u

∂x∂y+ c

∂2u

∂y2+ f

(

x, y, u,∂u

∂x,∂u

∂y

)

= 0

dove anche a, b, c possono essere funzioni di x, y, u e delle derivate parzialiprime di u. Generalmente le derivate parziali del secondo ordine possonoessere indicate anche in forma piu compatta:

uxx =∂2u

∂x2, uxy =

∂2u

∂y∂x, uyy =

∂2u

∂y2

e, analogamente quelle del primo ordine.

ux =∂u

∂x, uy =

∂u

∂y

1

CAPITOLO 1. EQUAZIONI ALLE DERIVATE PARZIALI 2

cosicche la forma dell’equazione appena vista diviene:

auxx + 2buxy + cuyy + f (x, y, u, ux, uy) = 0.

Le equazioni alle derivate parziali del secondo ordine:

a(x, y, u)uxx + 2b(x, y, u)uxy + c(x, y, u)uyy + f (x, y, u, ux, uy) = 0

sono le piu diffuse. Un’equazione alle derivate parziali si dice di ordine p sep e il massimo ordine di derivata che vi compare.Generalmente la scelta delle variabili indipendenti dipende dal problema fisi-co: infatti le variabili x, y, z indicano una posizione nello spazio, mentre lavaribile t indica il tempo. Talvolta anche le variabili x1, x2, x3 indicano unaposizione nello spazio. Considerando quindi le due equazioni

uxx + uyy + f(x, y, u) = 0, uxx + utt + f(x, t, u) = 0

esse sono matematicamente equivalenti ma fisicamente no, perche la primadescrive un fenomeno indipendente dal tempo (cioe stazionario) che riguardaun dominio bidimensionale (la posizione e descritta dalle variabili (x, y))mentre nel secondo caso il fenomeno descritto evolve nel tempo in un dominiomonodimensionale.Nell’equazione del secondo ordine

a(x, y, u)uxx + 2b(x, y, u)uxy + c(x, y, u)uyy + f (x, y, u, ux, uy) = 0

non compare la derivata uyx perche, in ipotesi di continuita, applicando ilTeorema di Schwarz, le derivate parziali miste sono uguali, cioe:

uxy = uyx.

Va anche precisato che nelle equazioni piu diffuse non e presente la deriva-ta uxy, perche talvolta non ha significato fisico mentre in altri casi con unopportuno cambiamento di variabile essa puo diventare nulla.

1.1.1 Operatori Differenziali

Spesso le equazioni alle derivate parziali sono rappresentate in forma piucompatta utilizzando determinati operatori differenziali, tra i quali:


1. Il Gradiente di u(x, y, t):

gradu(x, y, t) = ∇u(x, y, t) =

ux(x, y, t)uy(x, y, t)ut(x, y, t)

2. La Divergenza di una funzione vettoriale u(x, y, t) = (u1, u2, u3):

divu(x, y, t) =∂u1

∂x+∂u2

∂y+∂u3

∂t.

3. Il Laplaciano di u(x, y, t):

∆u = ∇2u(x, y, t) = div(grad(u(x, y, t)) =∂2u

∂x2+∂2u

∂y2+∂2u

∂t2.

4. Il Curl (o Rotore) della funzione vettoriale (o del campo vettoriale)u = (u1, u2, u3):

Curl(u) = gradu× u =

∂u3

∂y−∂u2

∂z

∂u1

∂z−∂u3

∂x

∂u2

∂x−∂u1

∂y

che trova applicazione nella terza equazione di Maxwell (la cosiddettaLegge di Faraday-Neumann-Lenz) che stabilisce che il rotore del cam-po elettrico e uguale e opposto alla derivata dell’intensita del campomagnetico e che, in generale, consente di descrivere il comportamentodi un campo vettoriale a ruotare rispetto ad un punto.

Descriviamo ora i piu noti esempi di equazioni del secondo ordine.

Esempio 1.1.1 L’equazione d’onda

∂2u

∂x2+∂2u

∂y2=

1

c2∂2u

∂t2


descrive, in funzione della posizione e del tempo, lo spostamento, rispettoal punto di equilibrio, di una corda vibrante. L’equazione descrive anche ilcampo elettrico o magnetico in un’onda elettromagnetica oppure l’intensita dicorrente oppure il potenziale lungo una linea di trasmissione. La quantita ce la velocita di propagazione dell’onda.

Esempio 1.1.2 L’equazione di diffusione

κ

(

∂2u

∂x2+∂2u

∂y2

)

=∂u

∂t

descrive la temperatura in una regione che non contiene sorgenti di calore, e siapplica anche alla diffusione di un composto chimico in un mezzo permeabile(liquido, mezzo poroso) avente concentrazione u(x, y, t). La costante κ vienedetta diffusivita.

Esempio 1.1.3 L’equazione di Laplace

∇2u = 0 ⇔∂2u

∂2x+∂2u

∂2y= 0

puo essere ottenuta dall’equazione di diffusione ponendo a zero la deriva-ta rispetto al tempo e descrive la distribuzione di temperatura in regimestazionario di un solido privo di sorgenti di calore. L’equazione di Laplacedescrive anche il potenziale elettrostatico in una regione priva di carica elet-trica. Si applica anche al flusso di un fluido incompromibile in una regionesenza vortici, sorgenti o scarichi.

Esempio 1.1.4 L’equazione di Poisson

∇2u = ρ(x, y) ⇔∂2u

∂2x+∂2u

∂2y= ρ(x, y)

descrive la stessa situazione dell’equazione di Laplace ma in una regione incui c’e carica elettrica oppure una sorgente di calore o di fluido. La funzioneρ(x, y) si chiama densita di sorgente e dipende anche da costanti fisiche. Peresempio nell’equazione di Poisson:

∂2u

∂2x+∂2u

∂2y= −

ρ(x, y)

ε

dove u(x, y) rappresenta il potenziale elettrostatico di una regione dello spazio,ρ(x, y) rappresenta la densita della carica elettrica e ε e la permittivita dellasostanza.


Esempio 1.1.5 Studiando il fenomeno della filtrazione dell’acqua nel sotto-suolo supponiamo che Ω sia una regione dello spazio occupata da un cosiddet-to mezzo poroso (terra o argilla per esempio), che abbia conduttivita idraulicaK. Indicando con φ(x, y, z) il livello dell’acqua, e con

q = (qx, qy, qz)

la velocita di filtrazione, allora applicando la legge di Darcy risulta che talevelocita e proporzionale alla variazione del livello dell’acqua:

q = −K

(

∂φ

∂x,∂φ

∂y,∂φ

∂z

)

e inoltre, per la proprieta di conservazione della massa, la divergenza di qdeve essere nulla:

divq = 0.

Applicando la definizione di divergenza segue:

divq =∂qx∂x

+∂qy∂y

+∂qz∂z

=

= −K

(

∂2φ

∂x2+∂2φ

∂y2+∂2φ

∂z2

)

= 0.

Da cui segue che la funzione φ soddisfa l’equazione di Laplace:

∆φ = 0.

La principale motivazione che spinge a risolvere numericamente le equazionialle derivate parziali sta nel fatto che non esistono tecniche analitiche generaliper ottenere la soluzione anche se per alcuni tipi di equazioni (soprattuttolineari) puo essere scritta esplicitamente sotto forma di serie di Fourier. Tut-tavia la convergenza di tali sviluppi in serie e lenta quindi per ottenere unabuona approssimazione e richiesto il calcolo di un numero di termini parti-colarmente elevato. Per altri tipi di equazioni si puo scrivere l’espressionedella soluzione teorica sotto forma di integrali che spesso non possono esserecalcolati esattemente ma solo approssimati numericamente.Considerando un’equazione alle derivate parziali del secondo ordine essa puoessere di tipo Lineare:

a(x, y)uxx+2b(x, y)uxy+c(x, y)uyy+d(x, y)ux+e(x, y)uy+f(x, y)u+g(x, y) = 0


oppure Quasi-Lineare:

a(x, y, u, ux, uy)uxx + 2b(x, y, u, ux, uy)uxy + c(x, y, u, ux, uy)uyy+

+f(x, y, u, ux, uy) = 0

oppure ancora Semi-Lineare:

a(x, y)uxx + 2b(x, y)uxy + c(x, y)uyy + f(x, y, u, ux, uy) = 0.

1.2 Classificazione delle equazioni alle deri-

vate parziali

Consideriamo l’equazione alle derivate parziali lineare

auxx + 2buxy + cuyy + dux + euy + fu+ g = 0, (1.1)

con (x, y) ∈ Ω ⊂ R2 e tale che a2 + b2 + c2 6= 0 per ogni (x, y) ∈ Ω. La

classificazione delle equazioni alle derivate parziali avviene in base al segnoassunto dalla quantita

∆ = b2 − ac.

Infatti un’equazione alle derivate parziali del secondo ordine si dice:

1. iperbolica se ∆ > 0,

2. ellittica se ∆ < 0,

3. parabolica se ∆ = 0.

Tale classificazione (per la verita un po’ superata) dipende solo dall’analogiaformale tra la (1.1) e l’equazione completa di una conica

ax2 + 2bxy + cy2 + dx+ ey + f = 0,

che, in funzione del valore assunto dal discriminante ∆, rappresenta

1. un’iperbole se b2 − ac > 0,

2. un’ellisse se b2 − ac < 0,

3. una parabola se b2 − ac = 0.


Si tratta di una classificazione non univoca perche se i coefficienti a, b, c dipen-dono da x e y allora il tipo di equazione dipende dal dominio dell’equazione.Per esempio, considerando l’equazione del secondo ordine:

y(x2 + 1)uxx + (x2 − 1)uyy + 3x+ y = 0

con (x, y) ∈ Ω ⊂ R2 risulta

∆ = −y(x2 − 1)(x2 + 1)

e la situazione e schematizzata nel seguente grafico, in cui lungo le retterosse l’equazione e parabolica, nella zona verde e iperbolica, mentre nellazona arancio e ellittica.

x = −1 x = 1y

x

Ovviamente questo tipo di classificazione presenta anche un altro problema,in quanto consente di riconoscere solo equazioni alle derivate parziali delsecondo ordine. Un modo indubbiamente migliore, e sicuramente univoco,per classificare un’equazione alle derivate parziali e quello di farlo in base


al tipo di fenomeno che essa descrive. Le equazioni alle derivate parzialipossono essere divise sommariamente in due tipi: equazioni stazionarie, incui tutte le variabili sono spaziali, ed equazioni di evoluzione, le quali pre-sentano una derivazione sia rispetto allo spazio che rispetto al tempo. Leequazioni di evoluzione rappresentano modelli matematici di fenomeni fisiciche subiscono variazioni nel tempo e sono molto importanti nella descrizionedei fenomeni ondulatori, termodinamici, nei processi di diffusione e nella di-namica delle popolazioni. Le equazioni alle derivate parziali, come dettoin precedenza, vengono anche classificate in ellittiche, paraboliche ed iper-boliche. Le equazioni ellittiche sono di tipo stazionario mentre le parabolichee le iperboliche sono equazioni di evoluzione. Le equazioni di evoluzionepossono essere viste come delle equazioni differenziali ordinarie senza vari-abili spaziali, infatti uno dei metodi piu efficaci per risolvere le equazioni dievoluzione e quello di approssimarle con un sistema di equazioni differenzialiordinarie. Questo tipo di classificazione puo essere ulteriormente raffinatoconsiderando altre classi di fenomeni fisici, per esempio i cosiddetti Problemidi avvezione che descrivono il moto di masse gassose (e che sono particolar-mente utilizzati nella descrizione di fenomeni meteorologici).Volendo classificare le equazioni introdotte in precedenza si puo osservare che

• L’equazione di Laplace

∂2u

∂x2+∂2u

∂y2= 0.

e Stazionaria ed Ellittica (∆ = −1);

• L’equazione di Poisson

∂2u

∂x2+∂2u

∂y2= ρ(x, y)

e Stazionaria ed Ellittica (∆ = −1);

• L’equazione d’onda

∂2u

∂x2− µ

∂2u

∂t2= 0, µ > 0

e di Evoluzione ed Iperbolica (∆ = µ);


• L’equazione di diffusione

κ∂2u

∂x2−∂u

∂t= 0,

e di Evoluzione e Parabolica (∆ = 0);

• L’equazione di Burgers

∂u

∂t+ u

∂u

∂x= 0, (1.2)

e di Evoluzione, ed ha numerose applicazioni in idrodinamica e gasdi-namica e nello studio delle onde acustiche.

Anche se l’ultimo tipo di equazione non rientra nella classificazione vista inprecedenza (infatti non e del secondo ordine), si puo verificare che e di tipoIperbolico, in quanto puo essere riscritta nel seguente modo:

∂2u

∂t2− u2∂

2u

∂x2+∂u

∂t

∂u

∂x− u

(

∂u

∂x

)2

= 0. (1.3)

Infatti derivando (1.2) rispetto a x e a t si ottiene:

∂2u

∂x∂t+

(

∂u

∂x

)2

+ u∂2u

∂x2= 0

e∂2u

∂t2+∂u

∂t

∂u

∂x+ u

∂2u

∂t∂x= 0.

Ricavando la derivata mista dalla prima equazione e sostituendo il valorenella seconda si ottiene la relazione (1.3). In questo caso ∆ = u2, da cuisegue che l’equazione e appunto iperbolica (purche u 6= 0).

Capitolo 2

Derivazione numerica

2.1 Introduzione

I metodi per la risoluzione numerica di equazioni alle derivate parziali chesaranno descritti nei prossimi capitoli sono detti Metodi alle differenze finitee sono basati sull’approssimazione discreta delle derivate parziali che com-paiono nell’equazione. In questo capitolo affronteremo quindi il problemarelativo all’approssimazione delle derivate prima e seconda di una funzionein un punto del dominio utilizzando opportune combinazioni lineari tra ivalori assunti dalla funzione in tale punto e in altri punti ad esso adiacenti.Nei successivi paragrafi considereremo il caso di una funzione in una singolavariabile. Per semplicita supponiamo che tale funzione, f(t), sia continuae differenziabile in un intervallo [a, b] fino ad un opportuno ordine k. In-izialmente considereremo il caso in cui l’intervallo [a, b] sia stato suddivisoin un insieme di sottointervalli di uguale ampiezza (griglia uniforme), perpoi considerare il caso in cui la funzione sia nota in un insieme di punti nonequidistanti (griglia non uniforme). Nei capitoli successivi sara descritto l’u-so di tali approssimazioni per le derivate parziali di funzioni definite su uninsieme Ω ⊂ R

2 (oppure Ω ⊂ R3).

2.2 Approssimazione discreta delle derivate

Come detto in precedenza supponiamo che f ∈ Ck([a, b]) e suddividiamol’intervallo di variabilita di t in sottointervalli di ampiezza h. Consideriamotre punti consecutivi appartenenti a tale reticolazione, rispettivamente ti−1,

10

CAPITOLO 2. DERIVAZIONE NUMERICA 11

ti e ti+1 tali cheti−1 = ti − h, ti+1 = ti + h.

Scriviamo lo sviluppo in serie di Taylor di f(ti+1) prendendo come puntoiniziale ti:

f(ti+1) = f(ti) + hf ′(ti) +h2

2f ′′(ti) +

h3

6f ′′′(ti) +

h4

24f iv(ξi), ξi ∈ [ti, ti+1]

e procediamo in modo analogo per f(ti−1):

f(ti−1) = f(ti) − hf ′(ti) +h2

2f ′′(ti) −

h3

6f ′′′(ti) +

h4

24f iv(ηi), ηi ∈ [ti−1, ti].

Sommiamo ora le due espressioni

f(ti+1) + f(ti−1) = 2f(ti) + h2f ′′(ti) +h4

24

[

f iv(ξi) + f iv(ηi)]

ricavando

f ′′(ti) =f(ti+1) − 2f(ti) + f(ti−1)

h2−h2

24

[

f iv(ξi) + f iv(ηi)]

e, trascurando l’ultimo termine, l’approssimazione della derivata seconda e:

f ′′(ti) ≃f(ti+1) − 2f(ti) + f(ti−1)

h2(2.1)

mentre si puo provare che l’errore vale:

E(f ′′(ti)) = −h2

12f iv(ξ), ξ ∈ [ti−1, ti+1].

Nel seguente grafico viene evidenziata l’interpretazione geometrica della for-mula appena ricavata.


ti ti+1ti−1

y = f(t)

Infatti l’approssimazione appena trovata coincide con il valore della deriva-ta seconda della parabola passante per i punti (ti−1, f(ti−1)), (ti, f(ti)) e(ti+1, f(ti+1)).Infatti scrivendo l’equazione di tale parabola come:

p(t) = a(t− ti)(t− ti−1) + b(t− ti−1) + c

risultac = f(ti−1)

b =f(ti) − f(ti−1)

h

a =f(ti+1) − 2f(ti) + f(ti−1)

2h2

e la proprieta segue poiche:

p′′(t) = 2a =f(ti+1) − 2f(ti) + f(ti−1)

h2.

Poniamoci il problema di approssimare derivata prima e procediamo nellostesso modo cioe scrivendo le serie di Taylor per f(ti+1) e f(ti−1) :

f(ti+1) = f(ti) + hf ′(ti) +h2

2f ′′(ti) +

h3

6f ′′′(σi), σi ∈ [ti, ti+1]


f(ti−1) = f(ti) − hf ′(ti) +h2

2f ′′(ti) −

h3

6f ′′′(µi), µi ∈ [ti−1, ti]

e questa volta sottraiamo la seconda dalla prima:

f(ti+1) − f(ti−1) = 2hf ′(ti) +h3

6[f ′′′(σi) + f ′′′(µi)]

ottenendo

f ′(ti) =f(ti+1) − f(ti−1)

2h−h2

12[f ′′′(σi) + f ′′′(µi)]

e, trascurando l’ultimo termine, l’approssimazione della derivata prima e:

f ′(ti) ≃f(ti+1) − f(ti−1)

2h(2.2)

mentre si puo provare che l’errore vale:

E(f ′(ti)) = −h2

6f ′′′(δ), δ ∈ [ti−1, ti+1].

ti ti+1ti−1

y = f(t)

m =f(ti+1) − f(ti−1)

2h

La formula (2.2) prende il nome di formula alle differenze centrali. Osservi-amo che sia per questa che per l’approssimazione numerica per la derivata


seconda l’errore dipende da h2, sono formule cioe del secondo ordine. Vedi-amo ora altre due approssimazioni per la derivata prima. Infatti possiamoanche scrivere:

f(ti+1) = f(ti) + hf ′(ti) +h2

2f ′′(ξi), ξi ∈ [ti, ti+1]

da cui si ricava immediatamente la formula alle differenze in avanti:

f ′(ti) ≃f(ti+1) − f(ti)

h(2.3)

con errore

E(f ′(ti)) = −h

2f ′′(ξi).

ti ti+1ti−1

y = f(t)

m =f(ti+1) − f(ti)

h

Analogamente si ricava

f(ti−1) = f(ti) − hf ′(ti) +h2

2f ′′(µi), µi ∈ [ti−1, ti]

da cui si ricava immediatamente la formula alle differenze all’indietro:

f ′(ti) ≃f(ti) − f(ti−1)

h(2.4)


con errore

E(f ′(ti)) = −h

2f ′′(µi).

Queste due formule hanno ordine 1, quindi sono meno precise rispetto al-la formula alle differenze centrali, tuttavia hanno il pregio di poter essereapplicate quando la derivata prima e discontinua in ti.

ti ti+1ti−1

y = f(t)

m =f(ti) − f(ti−1)

h

2.3 Approssimazioni di ordine superiore e di

derivate di grado superiore

Per determinare approssimazioni di ordine superiore per le derivate prima eseconda di una funzione di variabile reale e necessario aumentare il numerodi punti che sono coinvolti. A puro titolo di esempio calcoliamo un’approssi-mazione per la derivata seconda in ti che coinvolge due a punti a destra e duea sinistra (quindi in tutto 5 punti). Il modo di procedere risulta ovviamentemolto piu complesso e sicuramente meno intuitivo rispetto a quanto vistofinora. Si vogliono determinare i coefficienti della seguente relazione

f ′′(ti) ≃ αf(ti−2) + βf(ti−1) + γf(ti) + δf(ti+1) + εf(ti+2) (2.5)


in modo tale che l’approssimazione a sinistra abbia l’ordine p voluto, cioerisulti:

f ′′(ti) = αf(ti−2) + βf(ti−1) + γf(ti) + δf(ti+1) + εf(ti+2) + chpf (p)(ηi)

con ηi ∈]ti−2, ti+2[. Scriviamo ora gli sviluppi in serie di Taylor rispetto alpunto ti dei valori che compaiono in (2.5):

f(ti±2) = f(ti) ± 2hf ′(ti) + 2h2f ′′(ti) ±4h3

3f ′′′(ti) +

2h4

3f iv(ti) ±

4h5

15fv(σi)

(2.6)

f(ti±1) = f(ti)±hf′(ti)+

h2

2f ′′(ti)±

h3

6f ′′′(ti)+

h4

24f iv(ti)±

h5

120fv(µi). (2.7)

Sostituendo (2.6) e (2.7) in (2.5) e raccogliendo i termini con il medesimoordine di derivata si deve imporre che il risultato abbia nulli i coefficientidei termini differenziali di ordine 0,1,3 e 4, mentre quello di ordine 2 deveessere uguale a 1. In questo modo si arriva ad ottenere un sistema di cinqueequazioni nelle cinque incognite:

α +β +γ +δ +ε = 0

−2hα −hβ +hδ +2hε = 0

2h2α +h2

2β +h2

2δ +2h2ε = 1

−43h3α −h3

6β +h3

6δ +4

3h3ε = 0

23h4α +h4

24β +h4

24δ +2

3h4ε = 0

che ammette come soluzione

α = ε = −1

12h2, γ = −

5

2h2, β = δ =

4

3h2

che consente di ottenere la seguente approssimazione

f ′′(ti) ≃1

h2

[

−1

12f(ti−2) +

4

3f(ti−1) −

5

2f(ti) +

4

3f(ti+1) −

1

12f(ti+2)

]

(2.8)con un errore

E(f ′′(ti)) = ch4f (vi)(ξi), c ∈ R,


e pertanto di ordine 4, purche f(t) sia una funzione di classe Cvi([a, b]).Anche per la derivata prima e possibile ricavare approssimazioni discrete diordine maggiore coinvolgendo un numero di punti superiore. Consideriamo iseguenti sviluppi in serie di Taylor:

f(ti+1) = f(ti) + hf ′(ti) +h2

2f ′′(ti) +

h3

6f ′′′(ηi) (2.9)

f(ti+2) = f(ti) + 2hf ′(ti) + 2h2f ′′(ti) +4

3h3f ′′′(σi) (2.10)

Moltiplicando (2.9) per 4 e sottraendo la (2.10) si ottiene

4f(ti+1) − f(ti+2) = 3f(ti) + 2hf ′(ti) +2

3h3f ′′′(ηi) −

4

3h3f ′′′(σi),

che consente di ottenere la seguente approssimazione

f ′(ti) ≃4f(ti+1) − 3f(ti) − f(ti+2)

2h(2.11)

con errore

E(f ′′(ti)) =2

3h2 [f ′′′(ηi) − 2f ′′′(σi)] ,

da cui segue che la formula ha ordine 2.Procedendo come per la (2.8) e prendendo 2 punti a sinistra di ti e uno adestra, si puo ottenere, per la derivata prima, la seguente approssimazionedi ordine 3:

f ′(ti) ≃f(ti−2) − 6f(ti−1) + 3f(ti) + 2f(ti+1)

6h, (2.12)

e la seguente di ordine 4:

f ′(ti) ≃f(ti−2) − 8f(ti−1) + 8f(ti+1) − f(ti+2)

12h.

Come esempio si riporta nella seguente tabella l’errore che si commetteapprossimando il valore della derivata prima della funzione

f(x) = log x

in x = 4, utilizzando le formule introdotte in questo paragrafo e per valori dih decrescenti. Si osservi come la diminuzione dell’errore rispetto al valore dih sia molto piu accentuata per le formule di ordine piu elevato.


h (2.2) (2.3) (2.4) (2.11) (2.12)1 5.41E − 3 2.69E − 2 3.77E − 2 6.45E − 3 3.46E − 3

0.5 1.31E − 3 1.44E − 2 1.71E − 2 2.01E − 3 3.09E − 40.25 3.26E − 4 7.50E − 3 8.15E − 3 5.69E − 4 3.40E − 50.1 5.21E − 5 3.07E − 3 3.18E − 3 9.86E − 5 2.04E − 60.01 5.21E − 7 3.18E − 3 3.13E − 4 1.04E − 6 1.96E − 90.005 1.30E − 7 1.56E − 4 1.56E − 4 2.60E − 7 2.45E − 10

Nella seguente tabella sono riportati gli errori commessi nell’approssimazionedella derivata seconda per la stessa funzione f(x) = log x, in x = 4.

h (2.1) (2.8)

1 2.0385E − 3 4.2214E − 40.5 4.9343E − 4 2.1603E − 50.25 1.2239E − 4 1.2904E − 60.1 1.9539E − 5 3.2629E − 80.01 1.9531E − 7 4.9546E − 120.005 4.8825E − 8 1.4567E − 13

Approssimazione di derivate di ordine superiore

Per motivi di completezza citiamo al termine di questo paragrafo alcuneapprossimazioni di derivate di ordine superiore al secondo. In particolare perla derivata terza possiamo ottenere la seguente approssimazione, di ordine 2:

f ′′′(ti) ≃f(ti+2) − 2f(ti+1) + 2f(ti−1) − f(ti−2)

2h3.

Per quello che riguarda la derivata quarta la seguente approssimazione

f iv(ti) ≃f(ti+2) − 4f(ti+1) + 6f(ti) − 4f(ti−1) + f(ti−2)

h4

ha ordine 2. Bisogna considerare che, in realta, tali approssimazioni ven-gono utilizzate piuttosto raramente poiche le equazioni alle derivate parzialiche consideremo nel seguito, e che sono quelle incontrate piu spesso nelleapplicazioni, sono al piu del secondo ordine.


2.4 Approssimazione delle derivate su griglie

non uniformi

Su griglie non uniformi l’approssimazione della derivata prima e piuttostosemplice. Infatti supponendo di conoscere i valori (ti − h1, f(ti − h1)) e(ti + h2, f(ti + h2)) la derivata prima in ti puo essere approssimata con ilcoefficiente angolare della retta che passa per tali punti:

f ′(ti) ≃f(ti + h2) − f(ti − h1)

h1 + h2

Consideriamo, a puro titolo di esempio, il problema di approssimare la deriva-ta seconda della funzione f(t) nel punto di ascissa ti ma considerando i valoridella funzione nei punti non equidistanti ti−h e ti+sh, con 0 < s < 1. Svilup-piamo in serie di Taylor la funzione nel punto ti + sh prendendo come puntoiniziale ti

f(ti + sh) = f(ti) + shf ′(ti) +s2h2

2f ′′(ti) +

s3h3

6f ′′′(ξi), ξi ∈ [ti, ti + sh]

e procediamo in modo analogo per f(ti−1):

f(ti−1) = f(ti) − hf ′(ti) +h2

2f ′′(ti) −

h3

6f ′′′(ηi), ηi ∈ [ti−1, ti].

Moltiplichiamo per s quest’ultima espressione

sf(ti−1) = sf(ti) − shf ′(ti) +sh2

2f ′′(ti) −

sh3

6f ′′′(ηi)

e sommiamola con quella di f(ti + sh):

f(ti+sh)+sf(ti−1) = f(ti)(1+s)+h2

2f ′′(ti)s(1+s)+

sh3

6

[

s2f ′′′(ξi) − f ′′′(ηi)]

ricavando

f ′′(ti) = 2f(ti + sh) − f(ti)(1 + s) + sf(ti−1)

sh2(1 + s)+

h

3(1 + s)

[

f ′′′(ηi) − s2f ′′′(ξi)]

e, trascurando l’ultimo termine, l’approssimazione della derivata seconda e:

f ′′(ti) ≃ 2f(ti + sh) − f(ti)(1 + s) + sf(ti−1)

sh2(1 + s)(2.13)


mentre l’errore vale:

E(f ′′(ti)) =h

3(1 + s)

[

f ′′′(ηi) − s2f ′′′(ξi)]

.

Volendo ottenere un’approssimazione piu accurata si puo utilizzare la se-guente formula, della quale si omette la dimostrazione, che utilizza l’ulteriorevalore f(ti−2) ma e di ordine 2:

f ′′(ti) ≃1

h2

[

s− 1

s+ 2f(ti−2) +

2(2 − s)

s+ 1f(ti−1) −

3 − s

sf(ti) +

6f(ti + sh)

s(s+ 1)(s+ 2)

]

.

Capitolo 3

Equazioni ellittiche

3.1 L’equazione di Laplace

Vediamo ora di descrivere una tecnica per la risoluzione numerica della piuelementare equazione ellittica lineare, l’equazione di Laplace:

uxx + uyy = 0, (3.1)

ovvero

u ≡∂2u

∂x2+∂2u

∂y2= 0.

Se una funzione u(x, y) e di classe C2 in un determinato sottoinsieme Ω ⊆ R2

ed e una soluzione di (3.1) nello stesso Ω allora prende il nome di funzionearmonica. L a principale proprieta di queste funzioni e enunciata nel seguenteteorema.

Teorema 3.1.1 (Principio del massimo) Sia Ω una regione limitata esemplicemente connessa e sia Γ la sua frontiera. Sia Ω = Ω∪Γ. Se u(x, y) earmonica su Ω e continua su Ω, allora u(x, y) assume il suo valore massimosu Γ.

L’equazione di Laplace puo essere associata ad un problema di Dirichletquando, assegnata una funzione f(x, y) di classe C2(Γ) si cerca una funzioneu(x, y) tale che:

1. u(x, y) e continua su Ω ∪ Γ;

2. u(x, y) = f(x, y) per ogni (x, y) ∈ Γ;

21

CAPITOLO 3. EQUAZIONI ELLITTICHE 22

3. u(x, y) e armonica in Ω.

In alternativa si puo imporre la cosiddetta condizione di Neumann in cui, alposto della condizione 2., si impone che sia

∂u

∂n= f(x, y)

cioe sia assegnata la derivata normale di u(x, y) rispetto alla curva Γ. Ricor-diamo che se nT = (nx, ny), e il vettore normale allora

∂u

∂n= nx

∂u

∂x+ ny

∂u

∂y.

Consideriamo ora la risoluzione dell’equazione di Laplace prendendo Ω ugualeal rettangolo [a, b] × [c, d], con b > a e d > c. In questo caso un metodoe quello di approssimare l’operatore differenziale dopo avere discretizzatol’insieme Ω. Innanzitutto si suddivide l’intervallo [a, b] in N parti uguali,ognuna di ampiezza

h =b− a

Nper poi porre x0 = a, e

xi = xi−1 + k = a+ ih i = 1, 2, . . . , N.

Lo stesso procedimento si segue per le ordinate suddividendo l’intervallo [c, d]in M parti uguali, ciascuna di ampiezza

k =d− c

M

per poi porre y0 = c, e

yj = yj−1 + h = c+ jk j = 1, 2, . . . ,M.

In questo modo si ottiene l’insieme discreto di punti del piano

RN+1,M+1 =

(xi, yj) ∈ R2|xi = a+ ih, i = 0, N, yj = c+ jk, j = 0,M

.

La risoluzione numerica del problema di Dirichlet associato consiste nell’ap-prossimare opportunamente la funzione u(x, y) nei punti appartenenti all’in-sieme RN+1,M+1, tenendo presente che la soluzione e nota sul perimetro delrettangolo [a, b] × [c, d]. L’insieme discretizzato e evidenziato nella seguentefigura.


y

xa b

c

d

O

3.2 Il metodo a 5 punti per l’equazione di

Laplace

L’idea alla base del metodo a 5 punti e quella di approssimare le derivateparziali seconde nei punti del reticolo RN+1,M+1 e imporre che tali approssi-mazioni soddisfano l’equazione di Laplace. Poniamo innanzitutto

ui,j ≃ u(xi, yj), i = 0, 1, . . . , N, j = 0, 1, . . . ,M

e, per approssimare la derivata parziale seconda uxx(xi, yj), consideriamo iseguenti 3 punti del reticolo (xi−1, yj), (xi, yj) e (xi+1, yj) e, applicando laformula (2.1) supponendo costante il valore yj, risulta:

∂2u

∂x2(xi, yj) ≃

ui+1,j − 2ui,j + ui−1,j

h2.

Analogamente per approssimare uyy(xi, yj) consideriamo i seguenti 3 pun-ti del reticolo (xi, yj−1), (xi, yj) e (xi, yj+1) e, applicando la formula (2.1),risulta:

∂2u

∂y2(xi, yj) ≃

ui,j+1 − 2ui,j + ui,j−1

k2.


Adesso possiamo imporre che queste approssimazioni soddisfano l’equazionedi Laplace

ui+1,j − 2ui,j + ui−1,j

h2+ui,j+1 − 2ui,j + ui,j−1

k2= 0

(ui+1,j − 2ui,j + ui−1,j)k2 + (ui,j+1 − 2ui,j + ui,j−1)h

2 = 0.

Il metodo a cinque punti assume quindi la forma finale

h2ui,j−1 + k2ui−1,j − 2(h2 + k2)ui,j + k2ui+1,j + h2ui,j+1 = 0

con i = 1, . . . , N − 1, e j = 1, . . . ,M − 1.Tenendo presente che la funzione u(x, y) e nota sul bordo del rettangoloalcune delle approssimazioni non devono essere calcolate, infatti:

u0,j = u(x0, yj) = u(a, yj) = f(a, yj), j = 0, . . . ,Mui,0 = u(xi, y0) = u(xi, c) = f(xi, c), i = 0, . . . , NuN,j = u(xN , yj) = u(b, yj) = f(b, yj), j = 0, . . . ,Mui,M = u(xi, yM) = u(xi, d) = f(xi, d), i = 0, . . . , N.

Lo schema numerico e sintetizzato nel seguente stencil:

−2(h2 + k2) k2k2

h2

h2


3.2.1 Ordinamento delle incognite

Le relazioni che legano le incognite ui,j formano un sistema lineare la cuistruttura dipende dal modo con cui vengono ordinate tali incognite. Unprimo modo di riordinare le incognite ui,j e quello di porre

uT = (u1,1, u2,1, . . . , uN−1,1, u1,2, . . . , uN−1,2, . . . , u1,M−1, . . . , uN−1,M−1)

ed e schematizzato dal seguente grafico.

1 2 3 4 5 6 7 8

9 10 11 12 13 14 15 16

17 18 19 20 21 22 23 24

25 26 27 28 29 30 31 32

Tale ordinamento prende il nome di Ordinamento naturale (o lessicografico).Si inizia quindi dal punto (x1, y1), che viene numerato con 1, e si procede versodestra. Appena terminata la riga si passa a quella superiore, si considerail punto (x1, y2) e cosı via. La prima equazione (cioe quella per la primaincognita) si ottiene quindi per i = j = 1:

k2u1,0 + h2u0,1 − 2(h2 + k2)u1,1 + h2u2,1 + k2u1,2 = 0

equivalente a

−2(h2 + k2)u1,1 + h2u2,1 + k2u1,2 = −k2u1,0 − h2u0,1,

in quanto ogni volta che una delle relazioni coinvolge un punto sulla fron-tiera dell’insieme tale valore, essenco noto, contribuisce al termine noto del


sistema.La seconda equazione si ottiene per i = 2 e j = 1:

k2u2,0 + h2u1,1 − 2(h2 + k2)u2,1 + h2u3,1 + k2u2,2 = 0

equivalente a

h2u1,1 − 2(h2 + k2)u2,1 + h2u3,1 + k2u2,2 = −k2u2,0.

L’equazione relativa alla generica incognita di posto (i, j) ha al piu 5 coef-ficienti diversi da 0: 3 relativi alle 3 incognite numerate consecutivamente(i − 1, j), (i, j) e (i + 1, j), uno relativo ad una precedente, (i, j − 1), e unoad una successiva (i, j+ 1), numerate rispettivamente N − 1 prima ed N − 1dopo. Per ottenere le approssimazioni e quindi necessario risolvere un sistemalineare

Au = b

che ha la seguente struttura tridiagonale a blocchi

A =

T JJ T J

. . . . . . . . .

J T JJ T

dove J = k2IN−1, essendo IN−1 e la matrice identita di ordine N − 1, e T ela seguente matrice tridiagonale di dimensione N − 1:

T =

−2(h2 + k2) h2

h2 −2(h2 + k2) h2

. . . . . . . . .

h2 −2(h2 + k2) h2

h2 −2(h2 + k2)

.

Nella Figura 3.1 e riportata la struttura della matrice nel caso della grigliariportata come esempio in precedenza. Il sistema lineare da risolvere hauna struttura molto sparsa: se indichiamo con n la sua dimensione (n =(M − 1)(N − 1)) poco meno di 5n elementi sono diversi da zero su n2 ele-menti della matrice dei coefficienti. Questo significa che il sistema puo essere


risolto in modo piu efficiente utilizzando i cosiddetti metodi iterativi (vedereCapitolo 7) e non diretti (tipo fattorizzazione LU) che distruggerebbero lastruttura sparsa della matrice. Come esempio di questa affermazione nellaFigura 3.2 riportiamo la struttura della matrice triangolare superiore U dellafattorizzazione LU di A, come si puo osservare ha ben 231 elementi diversida zero (da confrontare con i 136 di A).

0 5 10 15 20 25 30

0

5

10

15

20

25

30

nz = 136

Ordinamento Lessicografico

Figura 3.1: Struttura della matrice dei coefficienti per l’ordinamentoLessicografico

L’ordinamento lessicografico non e l’unico modo di numerare le incognite delsistema lineare in oggetto, altri sono:

1. Ordinamento Cuthill-McKee;

2. Ordinamento Red-Black;

3. Ordinamento Multicolore.


0 5 10 15 20 25 30

0

5

10

15

20

25

30

nz = 231

Ordinamento Lessicografico−Fattore U

Figura 3.2: Struttura del fattore triangolare U per la matrice dei coefficientidefinita dall’ordinamento Lessicografico

Nell’ordinamento Cuthill-McKee, proposto nel 1969, si ordinano le incognitepartendo da un punto fissato e numerando i punti adiacenti in modo tale daprivilegiare quelli che si trovano lungo una determinata direzione (di solitoquella diagonale); se non ce ne sono si prende un altro nodo adiacente (o quel-lo a destra oppure quello in alto) e si prosegue. La tecnica e schematizzatanel seguente grafico.


Ordinamento di Cuthill-McKee

1 2 6 7 14 15 22 23

3 5 8 13 16 21 24 29

4 9 12 17 20 25 28 30

10 11 18 19 26 27 31 32

L’obiettivo di questo ordinamento e quello di diminuire l’ampiezza di bandadella matrice. Infatti in Figura 3.3 e visualizzata la struttura della matrice Acon tale ordinamento mentre la Figura 3.4 mostra la struttura della matriceU triangolare superiore. Osserviamo come il numero di elementi diversi dazero sia piuttosto basso (148) rispetto al caso precedente.Altri due ordinamenti, cui accenniamo per motivi di completezza prevedonoche i nodi vengano suddivisi (colorati) in due (o piu) tipi. Nell’ordinamen-to Red-Black sono divisi in rossi e neri, in modo tale che due nodi adiacentiabbiano colori diversi. Finita la colorazione i nodi sono numerati usando l’or-dinamento lessicografico applicato ad un colore per volta, come si evidenzianel seguente schema.


0 5 10 15 20 25 30

0

5

10

15

20

25

30

nz = 136

Ordinamento Cuthill−McKee

Figura 3.3: Struttura della matrice dei coefficienti per l’ordinamento Cuthill-McKee

Ordinamento di Red-Black

1 17 2 18 3 19 4 20

21 5 22 6 23 7 24 8

9 25 10 26 11 27 12 28

29 13 30 14 31 15 32 16


0 5 10 15 20 25 30

0

5

10

15

20

25

30

nz = 148

Ordinamento Cuthill−McKee−Fattore U

Figura 3.4: Struttura del fattore triangolare U per la matrice dei coefficientidefinita dall’ordinamento Cuthill-McKee

In Figura 3.5 e visualizzata la struttura della matrice A con tale ordinamentomentre la Figura 3.6 mostra la struttura della matrice U triangolare superi-ore. Osserviamo come il numero di elementi diversi da zero (170) sia piu altorispetto all’ordinamento Cuthill-McKee ma inferiore rispetto all’ordinamentolessicografico.L’ordinamento multicolore e uguale al Red-Black ma si usano piu colori, soli-tamente 4 o 6, e comunque un numero pari. I nodi sono ordinati colorandoliin sequenza uno di ogni colore diverso, mentre per la riga successiva si partedalla coppia di colori successiva, in modo tale che su ogni colonna ci sianosolo due colori. Nella figura seguente vediamo un esempio con 4 colori.


0 5 10 15 20 25 30

0

5

10

15

20

25

30

nz = 136

Ordinamento Red−Black

Figura 3.5: Struttura della matrice dei coefficienti per l’ordinamento Red-Black

Ordinamento Multicolore con 4 colori

1 9 17 25 2 10 18 26

19 27 3 11 20 28 4 12

5 13 21 29 6 14 22 30

23 31 7 15 24 32 8 16


0 5 10 15 20 25 30

0

5

10

15

20

25

30

nz = 170

Ordinamento Red−Black−Fattore U

Figura 3.6: Struttura del fattore triangolare U per la matrice dei coefficientidefinita dall’ordinamento Red-Black

Nel caso in cui i colori siano sei il primo colore della seconda riga coincidecon il quarto della sequenza, cosicche sulla stessa colonna siano accoppiatisempre il primo ed il quarto, il secondo con il quinto ed il terzo con il sesto.Se sono otto allora i colori accoppiati sono il primo con il quinto, il secondocon il sesto e cosı via.In Figura 3.7 e visualizzata la struttura della matrice A con tale ordinamentomentre la Figura 3.8 mostra la struttura della matrice U triangolare superi-ore. Osserviamo come il numero di elementi diversi da zero (191) sia piu altorispetto agli ordinamenti Cuthill-McKee e Red-Black ma inferiore rispettoall’ordinamento lessicografico.Se il dominio e quadrato, oppure si puo scegliere h = k, l’espressione delmetodo a 5 punti si semplifica:

ui,j−1 + ui−1,j − 4ui,j + ui+1,j + ui,j+1 = 0


0 5 10 15 20 25 30

0

5

10

15

20

25

30

nz = 136

Ordinamento Multicolor

Figura 3.7: Struttura della matrice dei coefficienti per l’ordinamento a 4colori

e gli elementi della matrice A sono indipendenti da h e k:

A =

T IN−1

IN−1 T IN−1

. . . . . . . . .

IN−1 T IN−1

IN−1 T


0 5 10 15 20 25 30

0

5

10

15

20

25

30

nz = 191

Ordinamento Multicolor−Fattore U

Figura 3.8: Struttura del fattore triangolare U per la matrice dei coefficientidefinita dall’ordinamento Multicolore con 4 colori

dove IN−1 e la matrice identita di ordine N − 1 e T e la seguente matricetridiagonale di dimensione N − 1:

T =

−4 11 −4 1

. . . . . . . . .

1 −4 11 −4

.

Esercizio. Ordinare le incognite dell’equazione di Laplace con condizione diDirichlet definita nel dominio in figura utilizzando gli ordinamenti lessicogra-fico e Red-Black. Supponendo di risolvere l’equazione applicando il metodoa 5 punti schematizzare, a fianco del dominio, la struttura della matrice deicoefficienti del relativo sistema lineare.


Ordinamento Lessicografico

1 2 3 4 5 6

7 8 9 10 11 12

13 14

15 16

17 18


1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Ordinamento Red-Black

1 10 2 11 3 12

13 4 14 5 15 6

7 16

17 8

9 18


1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

3.3 Convergenza del metodo a 5 punti

L’esistenza della soluzione numerica fornita dal metodo a 5 punti e con-seguenza dell’unicita della soluzione del sistema lineare definito dal metodo.Piu complessa risulta l’analisi dell’errore, cioe la valutazione dell’affidabilitadella soluzione numerica rispetto a quella teorica. Infatti e intuitivo che l’in-troduzione delle approssimazioni delle derivate parziali e causa di errore nellasoluzione numerica. Supponendo, per semplicita che h = k, allora si dice cheun metodo numerico e convergente se

limh→0

maxi,j

|u(xi, yj) − uij| = 0,

cioe se la massima differenza tra la soluzione numerica e quella teorica tendea zero al tendere a zero del passo di discretizzazione. Definiamo il cosiddet-to errore di troncamento locale come l’errore che si ottiene sostituendo lasoluzione teorica al posto di quella numerica nello schema:

τh(xi, yj) =u(xi, yj−1) + u(xi−1, yj) − 4u(xi, yj) + u(xi+1, yj) + u(xi, yj+1)

h2.

Condizione necessaria per la convergenza e che il metodo sia consistente,condizione che si verifica quando l’errore di troncamento locale tende a zero


quando h→ 0. In questo caso, considerando che la funzione u(x, y) soddisfal’equazione di Laplace possiamo scrivere l’errore di troncamento locale come

τh(xi, yj) =

[

u(xi−1, yj) − 2u(xi, yj) + u(xi+1, yj)

h2−∂2u(xi, yj)

∂x2+

+u(xi, yj−1) − 2u(xi, yj) + u(xi, yj+1)

h2−∂2u(xi, yj)

∂y2

]

.

da cui segue la consistenza del metodo poiche l’errore nell’approssimazionedelle derivate tende a zero quando h→ 0. In realta vale il seguente risultatoda cui segue la convergenza dello schema.

Teorema 3.3.1 Sia u ∈ C4(Ω), allora esiste una costante C ∈ R tale che

maxi,j

|u(xi, yj) − uij| ≤ CMh2

dove M e il massimo valore assoluto assunto dalle derivate quarte di u in Ω.

Esempio 3.3.1 Applichiamo il metodo a 5 punti all’equazione di Laplace

uyy + uxx = 0

definita sul rettangolo [0, 4π] × [0, 2] e con condizioni iniziali:

u(x, 0) = 2 sinx, 0 ≤ x ≤ 4πu(x, 2) = (e2 + e−2) sinx, 0 ≤ x ≤ 4πu(0, y) = 0, 0 ≤ y ≤ 2u(4π, y) = 0, 0 ≤ y ≤ 2

che ammette come soluzione teorica u(x, t) = sinx(ey + e−y).Nella Figura 3.9, viene riportato l’andamento della soluzione numerica pren-dendo N = M = 20, (il numero di punti della discretizzazione e pari, circa,a 1600).


02

46

810

1214

0

0.5

1

1.5

2

−10

−5

0

5

10

x

y

Figura 3.9: Soluzione dell’equazione di Laplace ottenuta con il metodo a 5punti.

3.4 Equazione di Laplace su domini non po-

ligonali

L’uso delle differenze divise funziona bene quando il dominio Ω e un rettan-golo, oppure un quadrato o un poligono che puo essere scomposto come ununione di quadrati o rettangoli. Quando invece il contorno Γ del dominio diintegrazione Ω e un poligono oppure una curva regolare a tratti l’approssi-mazione delle derivate parziali e piuttosto problematica. Vediamo prima uncaso particolare.


3.4.1 Il caso del cerchio unitario

Se il dominio Ω coincide con un cerchio e quindi Γ e una circonferenza, allorae possibile ricondursi ad un dominio rettangolare. Infatti, supponiamo, persemplicita, che sia

Γ =

(x, y) ∈ R2 : x2 + y2 = 1

mentre Ω e, ovviamente, il cerchio

Ω =

(x, y) ∈ R2 : x2 + y2 < 1

.

In questo caso il dominio puo essere trasformato in un rettangolo cambiandole coordinate cartesiane in polari:

x = ρ cos θ, y = ρ sin θ.

In coordinate polari il dominio diventa

Ω =

(ρ, θ) ∈ R2 : 0 ≤ ρ < 1, 0 ≤ θ < 2π

mentreΓ =

(ρ, θ) ∈ R2 : ρ = 1

.

Per evitare confusione poniamo

v(ρ, θ) = u(ρ cos θ, ρ sin θ).

Nella seguente figura e riportata la corrispondenza tra il dominio circolare incoordinate cartesiane ed il rettangolo in coordinate polari. Ad ogni circon-ferenza di raggio r < 1 viene associato un segmento parallelo all’asse θ, diequazione ρ = r e 0 ≤ θ < 2π. Ad ogni raggio della circonferenza unitariaviene invece associato un segmento parallelo all’asse ρ, di equazione θ = θ0,ed 0 ≤ ρ ≤ 1, (se θ0 e l’angolo che il raggio forma con l’asse delle ascisseorientato nel verso positivo).


ρx

θy

La condizione al contorno e assegnata ovviamente su Γ cioe:

u(x, y) = u(cos θ, sin θ) = f(cos θ, sin θ).

L’equazione di Laplace in coordinate polari diventa

∂2v

∂ρ2+

1

ρ2

∂2v

∂θ2+

1

ρ

∂v

∂ρ= 0.

Il problema piu delicato riguarda la trasformazione delle condizioni inizialiin coordinate polari. Il caso piu semplice e il segmento ρ = 1 che coincidecon la circonferenza Γ in coordinate cartesiane, quindi

v(1, θ) = f(cos θ, sin θ), 0 ≤ θ ≤ 2π.

I segmenti che si ottengono per 0 < ρ < 1 e θ = 0 e 0 < ρ < 1 e θ = 2π sonoesattamente lo stesso segmento nel cerchio originale. Il valore della funzionesu tale segmento non e noto percio e necessario assegnare una cosiddettacondizione di periodicita:

v(r, 0) = v(ρ, 2π), 0 < ρ < 1.


L’ultimo segmento da considerare e quello che si ottiene assegnando il valoreρ = 0, 0 ≤ θ ≤ 2π. L’intero segmento corrisponde ad un unico punto, cioel’origine del piano cartesiano. Quindi la funzione v deve essere costantelungo tale linea e, per esprimere tale condizione imponiamo la condizione diNeumann:

∂v

∂θ(0, θ) = 0.

Una volta determinate le condizioni sulla frontiera del dominio nel piano po-lare e possibile risolvere numericamente l’equazione di Laplace in coordinatepolari usando le approssimazioni alle differenze divise per le derivate parziali.

3.4.2 Il caso di un qualsiasi dominio irregolare

Analizziamo ora il caso in cui la frontiera del dominio Ω sia una curva chiusae regolare senza una particolare forma. In questo caso si considera un ret-tangolo [a, b]× [c, d] tale da contenere sia Ω che Γ e si discretizza tale regionecome gia visto in precedenza, definendo l’insieme discreto:

RN+1,M+1 = (xi, yj) ∈ [a, b] × [c, d]|i = 0, . . . , N, j = 0, . . . ,M .

y

xa b

c

d

Ω

Γ

O


L’insieme dei punti discreti che appartengono sia al rettangolo [a, b] × [c, d]che al dominio Ω si denota con

Ωhk = RN+1,M+1 ∩ Ω.

Tale insieme viene detto insieme dei punti interni. Ogni punto interno (xi, yj)ha quattro punti vicini, cioe (xi±1, yj) e (xi, yj±1).Un punto vicino ad un punto interno che non appartiene a Ωhk si dice puntodi confine. L’insieme dei punti di confine si indica con Γhk.I punti angolari sono invece i punti dell’insieme RN+1,M+1 che non sono neinterni ne di confine, ma risultano essere vertici di una cella che contienealmeno un punto interno.

rs

rsbc

bc

bc

rs

bc

bc

ut

ut

ut

rs

ut

ut

rs

Γ

Ω

Punti di confiners

Punti angolariut

Punti internibc

I punti angolari non hanno alcun ruolo nella soluzione numerica di equazioniellittiche. L’ipotesi che si deve fare sulla griglia che si considera e che isegmenti che congiungono punti interni devono essere interamente contenutinel dominio Ω. Si evita cioe il caso evidenziato dal secondo grafico:

(xi+1, yj)(xi, yj)(xi, yj) (xi+1, yj)


Tale situazione puo essere evitata scegliendo opportunamente il passo didiscretizzazione oppure effettuando un opportuno cambio di variabile (peresempio una rotazione degli assi). Il problema dei domini irregolari sorgequando nelle approssimazioni delle derivate seconde intervengono valori del-la funzione calcolati nei punti di confine. Per ovviare a tale inconveniente cisono diverse tecniche, due sono le piu usate:

1. come valori nei punti di confine si utilizzano gli stessi valori assuntidalla condizione al contorno;

2. si utilizza il valore della condizione al contorno nel punto sulla curva Γpiu vicino al punto interno.

Nel primo caso, considerando la seguente situazione:

bc bc

bcΓ

(xi, yj)

(xi, yj+1)

(xi+1, yj)

si pone:

u(xi, yj+1) = f(xi, yj+1), u(xi+1, yj) = f(xi+1, yj).

Tale assegnazione rappresenta, dal punto di vista matematico, una forzatura,poiche in realta non e noto neanche se la funzione f(x, y) sia calcolabile intali punti. Nel secondo caso, evidenziato dal seguente grafico, si utilizza ilpunto nel quale si conosce il valore della soluzione.


bc bcbc

Γ

xi − h xi xi + hxi + sh

Applicando la formula (2.13) si otterrebbe l’approssimazione:

uxx(xi + sh, yj) ≃ 2u(xi + sh, yj) − (1 + s)u(xi, yj) + su(xi−1, yj)

sh2(1 + s).

Tale approssimazione risulta essere del primo ordine quindi meno precisarispetto all’approssimazione della derivata seconda su punti equidistanti.Questo risultato ha come conseguenza un’approssimazione numerica menoprecisa in prossimita del contorno del dominio di integrazione.

3.5 Metodi di ordine elevato per equazioni

ellittiche

Consideriamo per semplicita l’equazione ellittica di Poisson

∂2u

∂x2+∂2u

∂y2= g(x, y)

definita sul quadrato Ω = [0, 1] × [0, 1].Il metodo a 5 punti ovviamente puo essere applicato anche all’equazione diPoisson. La matrice dei coefficienti e la stessa che abbiamo visto nel casodell’equazione di Laplace, cambia il termine noto del sistema lineare che euguale ai valori della funzione g(x, y) calcolata nei punti del dominio discreto.Avendo definito l’equazione su un quadrato possiamo scegliere un numero dipunti uguale sia per la variabile x che per y (e quindi h = k), ottenendo unamatrice dei coefficienti indipendente dal valore del passo di discretizzazione.


Il metodo a 5 punti ha ordine 2 perche l’errore di discretizzazione dipendeda h2.Si possono ottenere metodi di ordine superiore cercando di approssimare inmodo piu accurato le derivate parziali seconde.Utilizzando l’approssimazione (2.8) per le derivate parziali seconde porta alseguente metodo

−1

12ui,j−2 +

4

3ui,j−1 −

1

12ui−2,j +

4

3ui−1,j − 5ui,j +

4

3ui+1,j+

−1

12ui+2,j +

4

3ui,j+1 −

1

12ui,j+2 = h2g(xi, yj)

Il metodo puo essere descritto dal seguente stencil:

−5

4

3

4

3

4

3

4

3

− 1

12

− 1

12

− 1

12− 1

12uij = h2g(xi, yj)

Il metodo e piu preciso del metodo a 5 punti ma ci sono problemi quandoviene applicato in prossimita della frontiera in quanto non e disponibile unnumero di punti sufficiente.Il sistema lineare

Au = b


che si deve risolvere ha una struttura piu complessa di quello ottenuto con ilmetodo a 5 punti. La matrice dei coefficienti e la seguente

A =

T0 J0 OJ1 T J1 J2

J2 J1 T J1 J2

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

J2 J1 T J1 J2

J2 J1 T J1

O J0 T0

dove, utilizzando l’ordinamento lessicografico, la matrice T ha una strutturaa 5 diagonali (pentadiagonale), la matrice T0 e tridiagonale, mentre le matriciJ0, J1 e J2 sono matrici diagonali.

3.6 Il Metodo a 9 punti

Utilizzando una diversa approssimazione delle derivate parziali seconde dellafunzione u(x, y) si puo derivare uno schema numerico che utilizza nove punticollocati intorno al nodo di riferimento. Senza descrivere in dettaglio taletecnica, lo stencil del metodo applicato con lo stesso passo di discretizzazionerispetto alle variabili x e y, e il seguente:

- 10

3

2

3

2

3

2

3

2

3

1

6

1

6

1

6

1

6

uij = h2g(xi, yj)

Il metodo e preciso quanto il metodo a 5 punti (l’errore dipende da h2) mal’uso di maggiori informazioni (l’approssimazione coinvolge appunto 9 punti)


consente di ottenere in pratica risultati migliori;Il sistema lineare

Au = b

che si deve risolvere ha una struttura piu complessa di quello ottenuto con ilmetodo a 5 punti. La matrice dei coefficienti e la seguente

A =

T JJ T J

J T J. . . . . . . . .

J T JJ T J

J T

dove, utilizzando l’ordinamento lessicografico, le matrici T e J sono entrambetridiagonali.

Capitolo 4

Equazioni paraboliche

4.1 Equazioni di Evoluzione

Le equazioni di evoluzione sono, in modo molto semplicistico, le equazionialle derivate parziali che descrivono fenomeni dipendenti dal tempo e checontengono termini con derivate rispetto al tempo. Si tratta sostanzialmentedelle equazioni differenziali ordinarie alle quali sono aggiunte le variabilispaziali. Numericamente vengono risolte considerando prima una discretiz-zazione dell’intervallo temporale di osservazione del fenomeno, e, successi-vamente, vengono poste in relazione le approssimazioni in istanti di tempoconsecutivi. La conoscenza delle approssimazioni numeriche nei punti al tem-po tn consente di calcolare simultaneamente quelle in tn+1. I metodi possonoessere di due tipi:

1. Metodi Espliciti: partendo dalla conoscenza delle approssimazioni nu-meriche in un istante di tempo e possibile calcolare esplicitamente quelleall’istante successivo;

2. Metodi Impliciti: lo schema numerico mette simultaneamente in re-lazione le approssimazioni in due (o piu) istanti di tempo consecutivie, per calcolare l’approssimazione nell’istante di tempo successivo, enecessario risolvere un sistema lineare.

I metodi numerici descritti nel capitolo precedente per le equazioni ellit-tiche appartengono a questa seconda classe di metodi, in quanto, poiche lasoluzione e indipendente dal tempo, le approssimazioni sono calcolate simul-taneamente in tutti i punti del dominio dicreto.

50

CAPITOLO 4. EQUAZIONI PARABOLICHE 51

Una buona parte delle equazioni di evoluzione in una dimensione possonoessere descritte attraverso la cosiddetta equazione conservativa del flusso:

∂u

∂t= −

∂F (u)

∂x(4.1)

dove u ed F possono anche essere vettori, e, in alcuni casi, F puo dipenderenon solo da u ma anche dalle derivate parziali di u. La funzione (o il vet-tore) F e detto flusso conservativo. Le equazioni di evoluzione corrispondono(grossolanamente) a quelle di tipo parabolico e iperbolico. In questo capitolosaranno descritti alcuni metodi numerici per equazioni paraboliche, inizian-do da quella che viene considerata come il prototipo appartenente ad essa,l’equazione del calore.

4.2 Problemi di diffusione

I problemi di diffusione sono caratterizzati da un’equazione del tipo

∂u

∂t=

∂

∂x

[

σ∂u

∂x

]

(4.2)

dove σ e il cosiddetto coefficiente di diffusione (o anche dissipazione). Con-siderando l’equazione (4.1) questa si ottiene scegliendo

F (u) = −σ∂u

∂x.

Si deve supporre che sia σ > 0 altrimenti l’equazione (4.2) ha soluzioni fisi-camente instabili: la soluzione tenderebbe a concentrarsi e non a diffondersi.Se il coefficiente di diffusione e costante allora l’equazione (4.2) diventa

ut(x, t) = σuxx(x, t), (4.3)

che e nota come Equazione del Calore, supponendo t ≥ 0 e 0 ≤ x ≤ L,con L costante reale positiva. Supponiamo inoltre che sia stata assegnata lacondizione iniziale all’istante t = 0:

u(x, 0) = f(x), 0 ≤ x ≤ L

e le condizioni al contorno in x = 0 e x = L:

u(0, t) = g1(t) t ≥ 0u(L, t) = g2(t) t ≥ 0


in cui le funzioni f(x), g1(t) e g2(t) soddisfano le cosiddette condizioni diomogeneita:

f(0) = g1(0), f(L) = g2(0). (4.4)

t

x

O f(x)

g1(t) g2(t)

L

L’equazione del calore descrive fenomeni di diffusione termodinamica, infattiu(x, t) puo rappresentare, per esempio, l’evoluzione nel tempo della densitadi calore di una sbarra termoconduttrice di lunghezza L e spessore trascur-abile che viene sottoposta ad una certa temperatura iniziale (la condizioneal contorno f(x)) e tale che la temperatura agli estremi sia fissata dalle con-dizioni poste dalle funzioni g1(t) e g2(t). Se dopo l’istante iniziale la sorgentedi calore viene tolta allora la densita di calore u(x, t) del punto della barradi ascissa x al tempo t soddisfa l’equazione del calore.L’equazione di conducibilita del calore puo essere generalizzata in molti modi:

1. aggiungendo una variabile spaziale (cioe u = u(x, y, t)):

ut(x, y, t) = σ(uxx(x, y, t) + uyy(x, y, t));

2. aggiungendo un termine forzante f(x, t):

ut(x, t) = σuxx(x, t) + f(x, t);


3. aggiungendo un coefficiente di diffusione variabile a(x):

∂u

∂t=

∂

∂x

[

a(x)∂u

∂x

]

con a(x) funzione differenziabile e tale che 0 < a(x) <∞ per x ∈ [0, L];

4. prendendo x appartenente ad un intervallo arbitrario della retta reale esostituendo le condizioni al contorno con la condizione che la funzioneu(x, t) abbia quadrato integrabile, cioe:

∫ +∞

−∞

[u(x, t)]2 dx <∞, t ≥ 0,

condizione che si verifica considerando un’altra equazione paraboli-ca, l’equazione di Schrodinger, equazione alla base della meccanicaquantistica:

ι~∂u

∂t= −

~2

2m

∂2u

∂x2+ U(x)u

in cui m e la massa della particella, ι e l’unita immaginaria, ~ e lacostante di Plank, U(x) e l’energia potenziale della particella e u(x, t)e la funzione d’onda. In questo caso si deve richiedere che, oltre allecondizioni iniziali venga soddisfatta anche la seguene relazione

∫ +∞

−∞

|u(x, t)|2dx = 1,

in quanto il modulo della soluzione rappresenta una densita di proba-bilita.

Per l’equazione parabolica e possibile definire due tipi di problemi. Il primoe il problema ai valori iniziali, in cui si tratta di trovare una funzione u(x, t),definita e continua per x ∈ R e t ≥ 0, che soddisfi l’equazione alle derivateparziali per x ∈ R e t > 0 e la condizione iniziale u(x, 0) = f(x), x ∈ R,come schematizzato nella seguente figura.


t

xOu(x, 0) = f(x)

Il secondo e il problema ai valori al contorno, in cui, assegnata una costanteL > 0, si deve trovare una funzione u(x, t), definita e continua per 0 ≤ x ≤ Le t ≥ 0, che soddisfi l’equazione alle derivate parziali per 0 < x < L e t > 0e le condizioni iniziali:

u(x, 0) = f(x) 0 ≤ x ≤ Lu(0, t) = g1(t) t ≥ 0u(L, t) = g2(t) t ≥ 0.

che, a loro volta, devono soddisfare le condizioni di omogeneita (4.4).

t

xO u(x, 0) = f(x)

u(0, t) = g1(t)

u(L, t) = g2(t)

L


4.3 Il metodo di Eulero Esplicito per l’equa-

zione del calore

Supponiamo ora di dover risolvere l’equazione del calore

ut = σuxx, 0 ≤ x ≤ L, t ≥ 0,

con condizione iniziale

u(x, 0) = f(x), 0 ≤ x ≤ L (4.5)

e condizioni al contorno

u(0, t) = g1(t), u(L, t) = g2(t), t ≥ 0. (4.6)

Poiche t ≥ 0 e non e ovviamente possibile calcolare la soluzione all’infini-to si sostituisce t ≥ 0 con t ∈ [0, Tmax]. La costante Tmax e generalmentedeterminata dalla fisica del fenomeno in osservazione. L’equazione vieneintegrata numericamente negli istanti di tempo tn ≤ Tmax. La risoluzionenumerica di questa equazione richiede la sua discretizzazione sia rispetto altempo che rispetto allo spazio. Il problema discreto e ottenuto mediante l’ap-prossimazione delle derivate parziali attraverso differenze finite. Scegliamoun intero positivo Nx e definiamo nella striscia

(x, t) : x ∈ [0, L], t ≥ 0

una griglia rettangolare formata dai punti di coordinate (xj, tn) tale che

xj = j∆x, j = 0, 1, . . . , Nx + 1, tn = n∆t, n ≥ 0,

con

∆x =L

(Nx + 1), ∆t =

Tmax

Nt

.

Il valore ∆t rappresenta l’intervallo di tempo tra due approssimazioni suc-cessive. I punti (xj, tn) del dominio discreto sono di tre tipi, evidenziati nelseguente grafico.


t

x

O x1 x2 . . . xNx L

t1

t2

...

...

Punti al contorno

Punti iniziali

Punti interni

L’approssimazione di u(xj, tn) e denotata con unj . Utilizzando l’operatoredifferenza centrale possiamo approssimare la derivata parziale seconda nelmodo gia visto in precedenza:

∂2u

∂x2(xj, tn) ≃

unj+1 − 2unj + unj−1

(∆x)2

mentre utilizzando l’operatore differenza in avanti per la derivata temporale,si ottiene

∂u

∂t(xj, tn) ≃

un+1j − unj

∆t.

Sostituendo in (4.3) si ha

un+1j − unj

∆t= σ


(∆x)2

un+1j = unj +

σ∆t

(∆x)2

(


)

e infine si ha come risultato il seguente Metodo di Eulero Esplicito:

un+1j = αunj−1 + (1 − 2α)unj + αunj+1 (4.7)


dove

α =σ∆t

(∆x)2

e una costante che prende il nome di Numero di Courant. Partendo da n = 0la formula (4.7) consente di determinare esplicitamente le approssimazioniu1j , j = 1, . . . , Nx, pertanto si tratta di un metodo di tipo esplicito. In modo

analogo, conoscendo le approssimazioni al livello n, consente di calcolarequelle al livello n+ 1. La condizione iniziale (4.5) fornisce i valori

u0j = f(j∆x), j = 0, 1, 2, . . . , Nx + 1.

Quando j = 1 e j = Nx la formula (4.7) utilizza le condizioni sulla frontiera(4.6):

un0 = g1(tn), unNx= g2(tn).

4.3.1 Condizione di Stabilita per il metodo di EuleroEsplicito

Come e noto i metodi semplici non sono sempre quelli migliori, infatti consi-deriamo che la soluzione teorica di un’equazione parabolica possiede la stessaproprieta di massimo-minimo delle funzioni armoniche (cioe il massimo edil minimo della soluzione teorica si trovano sulla frontiera del dominio). Equindi auspicabile che anche la soluzione numerica possieda tale proprieta.Per analizzare tale proprieta consideriamo il caso in cui il dominio discretiz-zato ha un solo nodo interno, prendendo quindi ∆x = L/2, e come condizioneiniziale la funzione che assume valore ε > 0 nel nodo interno e 0 nei due puntial contorno. Il dominio discreto e rappresentato nella seguente figura.

t

x

O L/2 L

t1

t2

t3

...

...


Quindiu0

0 = u02 = 0, u0

1 = ε.

Applicando lo schema di Eulero Esplicito si ottiene:

u11 = (1 − 2α)u0

1 = (1 − 2α)ε

e, ad un generico istante di tempo tn:

un1 = (1 − 2α)nε.

Ora, in base alla proprieta di massimo-minimo, deve essere

0 ≤ un1 ≤ ε ⇒ 0 ≤ (1 − 2α)n ≤ 1

e quindi0 ≤ 1 − 2α ≤ 1.

Poiche α > 0 la seconda disuguaglianza e sicuramente verificata, quindi deveessere

1 ≥ 2α ⇒ α ≤1

2.

Tale disuguaglianza, supponendo σ = 1, implica la seguente restrizione suipassi di discretizzazione

∆t ≤(∆x)2

2.

Se ∆x ≃ 10−2 deve essere∆t ≤ 0.5 · 10−4

quindi se e necessario integrare l’equazione su un intervallo di tempo moltogrande si deve usare un valore ∆t molto piccolo e di conseguenza il numerodi passi temporali e incredibilmente grande.

4.4 Analisi di Stabilita di von Neumann

Come dimostra il semplice esempio descritto nel paragrafo precedente la sosti-tuzione delle derivate parziali con approssimazioni discrete non garantisce chele soluzioni numeriche siano affidabili, a meno di porre determinati vincoli suipassi di dicretizzazione. Per studiare questo fenomeno e necessario introdurreuna tecnica, detta analisi di stabilita di von Neumann. Tale analisi e di tipo


locale, nel senso che si suppone che i coefficienti dell’equazione che definisceil metodo (detta equazione alle differenze) siano costanti rispetto al tempo eallo spazio (o cambiano di poco). Si suppone inoltre che le soluzioni di taliequazioni siano della forma:

unj = ξneικj∆x (4.8)

dove κ e un numero reale, detto numero d’onda, e ξ e un numero complessodipendente da κ. Il numero ξ, ed in particolare le sue potenze, indicano ladipendenza della soluzione rispetto al tempo. Infatti la soluzione numericaunj cresce esponenzialmente se |ξ(κ)| > 1 per un certo valore di κ. Per questoil numero ξ e detto fattore di amplificazione per il numero d’onda κ. Lasemplice analisi della stabilita vista nel precedente paragrafo corrisponde,in realta, allo studio dei valori che puo assumere il numero di Courant inmodo tale che la soluzione numerica si mantenga limitata. In questo caso cioequivale a richiedere che il fattore di amplificazione abbia modulo minore di(o uguale a) 1. Per applicare l’analisi di von Neumann al metodo di EuleroEsplicito sostituiamo l’espressione (4.8) nello schema numerico (4.7):

ξn+1eικj∆x = ξneικj∆x + α[

ξneικ(j+1)∆x − 2ξneικj∆x + ξneικ(j−1)∆x]

.

Dividendo per ξneικj∆x si ottiene

ξ = 1 + α[

eικ∆x − 2 + e−ικ∆x]

= 1 + α [2 cos(κ∆x) − 2] .

Poiche

cos(κ∆x) = 1 − 2 sin2

(

κ∆x

2

)

⇒ cos(κ∆x) − 1 = −2 sin2

(

κ∆x

2

)

quindi

ξ = 1 − 4α sin2

(

κ∆x

2

)

.

Poiche ξ e un numero reale devono essere soddisfatte simultaneamente ledisequazioni:

−1 ≤ ξ ≤ 1.

In particolare la seconda risulta sempre soddisfatta, mentre dalla prima segueche

2α sin2

(

κ∆x

2

)

≤ 1

che e soddisfatta, come visto prima, se α ≤ 1/2.


4.5 θ−Metodi

Il metodo di Eulero esplicito fa parte di una pia ampia classe di scheminumerici, detti θ−Metodi perche dipendono da un parametro reale θ. Laderivata prima viene approssimata utilizzando la formula alle differenze inavanti:

ut (xj, tn) ≃un+1j − unj

∆t,

mentre la derivata parziale seconda rispetto a x viene approssimata attraversouna combinazione lineare tra i valori in due istanti di tempo consecutivi:

uxx (xj, tn) ≃ [(1 − θ)uxx(xj, tn) + θuxx(xj, tn+1)] ,

con θ ∈ R (in particolare θ ∈ [0, 1]). Approssimando le derivate seconde nelmodo consueto si ottiene la formula:

un+1j − unj

∆t= σ

[

(1 − θ)unj+1 − 2unj + unj−1

(∆x)2+ θ

un+1j+1 − 2un+1

j + un+1j−1

(∆x)2

]

,

da cui, indicato con α il numero di Courant, si ricava

un+1j = unj + (1 − θ)α

[


]

+ θα[

un+1j+1 − 2un+1

j + un+1j−1

]

.

Un metodo di questa classe puo essere rappresentato attraverso il seguentestencil, in cui i punti celesti rappresentano quelli coinvolti nell’approssi-mazione delle derivate seconde mentre i rossi quelli nell’approssimazione delladerivata prima.

(xj, tn)(xj−1, tn) (xj+1, tn)

(xj, tn+1) (xj+1, tn+1)(xj−1, tn+1)


Se θ = 0 il metodo coincide con il Metodo di Eulero Esplicito:

un+1j = unj + α

[


]


(xj, tn+1)Stencil per il Metodo di Eulero Esplicito

Se θ = 0.5 il metodo viene detto Metodo di Crank-Nicolson:

un+1j = unj +

α

2

[


]

+α

2

[

un+1j+1 − 2un+1

j + un+1j−1

]

,

ovvero

−α

2un+1j−1 + (1 + α)un+1

j −α

2un+1j+1 =

α

2unj−1 + (1 − α)unj +

α

2unj+1.


(xj, tn+1) (xj+1, tn+1)(xj−1, tn+1)Stencil per il Metodo di Crank-Nicolson


Se θ = 1 il metodo coincide con il Metodo alle Differenze Centrali Implicito:

un+1j = unj + α

[

un+1j+1 − 2un+1

j + un+1j−1

]

,

ovvero−αun+1

j−1 + (1 + 2α)un+1j − αun+1

j+1 = unj .

(xj, tn)

(xj, tn+1) (xj+1, tn+1)(xj−1, tn+1)Stencil per il Metodo alle Differenze Centrali Implicito

I metodi di Crank-Nicolson e alle differenze centrali implicito sono entrambidi tipo implicito, poiche ad ogni istante consentono di determinare simul-taneamente tutti i valori delle approssimazioni numeriche all’istante tn+1 ri-solvendo un sistema lineare a struttura tridiagonale (l’espressione del metodocoinvolge solo 3 incognite al livello n+ 1).Applicando l’analisi di stabilita di von Neumann si scopre che se θ ≥ 0.5allora la stabilita e incondizionata (non c’e alcuna restrizione sui passi diintegrazione), mentre se 0 ≤ θ < 0.5 la condizione di stabilita e la seguente

α ≤1

2(1 − 2θ).

Esempio 4.5.1 Applichiamo i metodi descritti alla seguente equazione delcalore:

ut − uxx = 0

con condizioni iniziali:

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


u(0, t) = u(1, t) = 0, t ≥ 0

che ammette come soluzione teorica u(x, t) = e−π2t sin πx.

Nelle Figure 4.1, 4.2 e 4.3, e stato applicato il metodo di Eulero esplicitonell’intervallo temporale [0, 0.4] con il numero di Courant che e, rispettiva-mente, inferiore, uguale o maggiore di 0.5. Nelle prime due figure si osservache la soluzione numerica concorda con quella teorica, anche se il passo didiscretizzazione e molto piccolo (si osservi la griglia temporale molto fitta).Nella terza figura la soluzione numerica presenta un andamento fortementeoscillatorio che non soddisfa la proprieta del massimo. Nelle Figure 4.4 e 4.5sono riportate le soluzioni numeriche fornite, rispettivamente, dai metodi alledifferenze centrali implicito e di Crank-Nicolson. Si osserva facilmente chele soluzioni numeriche concordano con quella teorica indipendentemente dalvalore del numero di Courant riportato nei grafici.

00.2

0.40.6

0.81

0

0.1

0.2

0.3

0.40

0.2

0.4

0.6

0.8

1

x

Metodo di Eulero Esplicito−α=0.25

t

Figura 4.1: Soluzione dell’equazione del calore ottenuta con il metodo diEulero esplicito con α = 0.25.


00.2

0.40.6

0.81

0

0.1

0.2

0.3

0.40

0.2

0.4

0.6

0.8

1

x

Metodo di Eulero Esplicito−α=0.5

t

Figura 4.2: Soluzione dell’equazione del calore ottenuta con il metodo diEulero esplicito con α = 0.5.

4.6 L’equazione del calore in due dimensioni

Aggiungendo una seconda variabile spaziale l’equazione del calore diventa

ut = σ(uxx + uyy), (x, y, t) ∈ D × R+

dove D e un sottoinsieme di R2, con condizioni al contorno

u(x, y, 0) = u0(x, y), (x, y) ∈ D

u(x, y, t) = g(x, y, t), (x, y) ∈ ∂D, t ≥ 0.

Tale equazione descrive, per esempio, la temperatura di una piastra bidimen-sionale la cui forma coincide con quella dell’insieme D e altezza trascurabile(rispetto alle altre dimensioni), che viene a contatto con una sorgente dicalore (descritta dalla funzione u0(x, y)) all’istante t = 0, e il cui bordo (cioe∂D) viene tenuto a temperatura costante.Nei paragrafi seguenti sono descritti alcuni semplici metodi numerici per larisoluzione numerica di tale equazione.


00.2

0.40.6

0.81

0

0.02

0.04

0.06

0.08

0.1−0.2

0

0.2

0.4

0.6

0.8

1

1.2

x

Metodo di Eulero Esplicito−α=2

t

Figura 4.3: Soluzione dell’equazione del calore ottenuta con il metodo diEulero esplicito con α = 2.

Il metodo di Eulero Esplicito in due dimensioni

Supponiamo che il dominio spaziale D sia il rettangolo [0, L] × [0,M ] cheviene discretizzato definendo, al solito, i valori

xi = i∆x, yj = j∆y,

mentre anche il tempo viene suddiviso nel solito modo

tn = n∆t, n = 0, 1, . . .

Siauni,j ≃ u(xi, yj, tn),

il metodo di Eulero Esplicito e la generalizzazione di quello gia incontrato

un+1i,j − uni,j

∆t= σ

[

uni+1,j − 2uni,j + uni−1,j

(∆x)2+uni,j+1 − 2uni,j + uni,j−1

(∆y)2

]

.


00.2

0.40.6

0.81

0

0.5

1

1.5

20

0.2

0.4

0.6

0.8

1

x

Metodo alle differenze centrali implicito−α=44

t

Figura 4.4: Soluzione dell’equazione del calore ottenuta con il metodo alledifferenze divise implicito.

Livello n

Livello n+ 1

ii− 1 i+ 1j − 1

j

j + 1


00.2

0.40.6

0.81

0

0.5

1

1.5

2−0.2

0

0.2

0.4

0.6

0.8

1

x

Metodo di Crank−Nicolson implicito−α=88

t

Figura 4.5: Soluzione dell’equazione del calore ottenuta con il metodo diCrank-Nicolson.

La condizione di stabilita diventa:

∆t

(∆x)2+

∆t

(∆y)2≤

1

2.

Osserviamo che tale condizione e ben piu restrittiva di quella vista in prece-denza, infatti supponendo di utilizzare lo stesso passo di discretizzazione siaper x che per y la condizione di stabilita diventa:

∆t

(∆x)2≤

1

4.

Il Metodo di Crank-Nicolson

Esattamente come nel caso unidimensionale e possibile definire la classe deiθ−metodi. Nel caso del metodo di Crank-Nicolson la derivata seconda nelpunto della griglia di coordinate (xi, yj, tn) viene approssimato con il valor


medio tra due istanti di tempo successivi:

ut(xi, yj, tn) ≃ σ1

2[uxx(xi, yj, tn+1) + uyy(xi, yj, tn+1)]

+σ1

2[uxx(xi, yj, tn) + uyy(xi, yj, tn)]

Il metodo assume la seguente espressione

un+1i,j − uni,j

∆t=σ

2

[

un+1i+1,j − 2un+1

i,j + un+1i−1,j

(∆x)2+un+1i,j+1 − 2un+1

i,j + un+1i,j−1

(∆y)2

]

+σ

2

[

uni+1,j − 2uni,j + uni−1,j

(∆x)2+uni,j+1 − 2uni,j + uni,j−1

(∆y)2

]

.

Livello n

Livello n+ 1

ii− 1 i+ 1j − 1

j

j + 1

Il Metodo alle Differenze Centrali Implicito

Il metodo alle differenze centrali implicito consiste nell’approssimare al livellotemporale n+ 1 le derivate parziali seconde:

un+1i,j − uni,j

∆t= σ

[

un+1i+1,j − 2un+1

i,j + un+1i−1,j

(∆x)2+un+1i,j+1 − 2un+1

i,j + un+1i,j−1

(∆y)2

]

.


Ad ogni passo si deve risolvere un sistema lineare con struttura pentadiag-onale. Il metodo e incondizionatamente stabile (non ci sono restrizioni sulpasso temporale di discretizzazione).

Livello n

Livello n+ 1

ii− 1 i+ 1j − 1

j

j + 1

Il Metodo delle Direzioni Alternate

Si effettuano due mezzi passi temporali, uno in ogni direzione:

un+1/2i,j − uni,j

∆t/2= σ

[

Lxxun+1/2i,j + Lyyu

ni,j

]

(4.9)

un+1i,j − u

n+1/2i,j

∆t/2= σ

[

Lxxun+1/2i,j + Lyyu

n+1i,j

]

(4.10)

dove Lxx ed Lyy indicano le approssimazioni delle derivate seconde rispettoad x e ad y. Esplicitando la relazione (4.9) si ottiene

un+1/2i,j − uni,j

∆t/2= σ

[

un+1/2i−1,j − 2u

n+1/2i,j + u

n+1/2i+1,j

(∆x)2+uni,j−1 − 2uni,j + uni,j+1

(∆y)2

]

,

da cui, posto α = σ∆t/(2(∆x)2) e β = σ∆t/(2(∆y)2) si ottiene la relazionerelativa al primo mezzo passo temporale:

−αun+1/2i−1,j + (1 + 2α)u

n+1/2i,j − αu

n+1/2i+1,j = βuni,j−1 + (1 − 2β)uni,j + βuni,j+1.


In modo simile si ricava l’espressione relativa al secondo mezzo passo tempo-rale (4.10):

−βun+1i,j−1 + (1 + 2β)un+1

i,j − βun+1i,j+1 = αu

n+1/2i−1,j + (1 − 2α)u

n+1/2i,j + αu

n+1/2i+1,j .

Il metodo e implicito ma ad ogni passo temporale si devono risolvere due siste-mi lineari che hanno struttura tridiagonale, ed inoltre e incondizionatamentestabile.

Livello n

Livello n+ 1/2

Livello n+ 1

ii− 1 i+ 1j − 1

j

j + 1

Lo Schema di Yanenko

Si effettuano due discretizzazioni temporali con passo dimezzato ognunarispetto ad una direzione diversa (la prima verso x la seconda verso y).

un+1/2i,j − uni,j

∆t/2= σLxxu

n+1/2i,j (4.11)

un+1i,j − u

n+1/2i,j

∆t/2= σLyyu

n+1i,j (4.12)

Esplicitando la relazione (4.11) si ottiene

un+1/2i,j − uni,j

∆t/2= σ

un+1/2i−1,j − 2u

n+1/2i,j + u

n+1/2i+1,j

(∆x)2


da cui, posto α = σ∆t/(2(∆x)2) si ottiene la relazione relativa al primomezzo passo temporale:

−αun+1/2i−1,j + (1 + 2α)u

n+1/2i,j − αu

n+1/2i+1,j = uni,j.

Posto β = σ∆t/(2(∆y)2) ed esplicitando la relazione (4.12) si ottiene l’e-spressione eaplicita del secondo mezzo passo temporale:

−βun+1i,j−1 + (1 + 2β)un+1

i,j + βun+1i,j+1 = u

n+1/2i,j .

Ad ogni passo temporale si devono risolvere due sistemi lineari che hannostruttura tridiagonale, ed inoltre il metodo e incondizionatamente stabile.

Il Metodo di Hopscotch

E uno schema che e composto da due parti, una esplicita e l’altra implicita.

un+1i,j − uni,j

∆t= σ

[

Lxxuni,j + Lyyu

ni,j

]

, se i+ j + n pari;

un+1i,j − uni,j

∆t= σ

[

Lxxun+1i,j + Lyyu

n+1i,j

]

, se i+ j + n dispari.

Se n e pari e si pone σα = ∆t/(∆x)2 e β = σ∆t/(∆x)2, allora si ottiene laseguente espressione se i+ j e pari:

un+1i,j = βuni,j−1 + αuni−1,j + (1 − 2α− 2β)uni,j + αuni+1,j + βuni−1,j,

e la seguente se i+ j e dispari:

un+1i,j = uni,j +

1

1 + 2α+ 2β

[

βun+1i,j−1 + αun+1

i−1,j + αun+1i+1,j + βun+1

i−1,j

]

.

Quando n e dispari allora la prima formula vale se i+ j e dispari e la secondase i + j e pari. Il metodo e poco costoso, in quanto e esplicito, poiche nellaseconda relazione vengono utilizzate le approssimazioni al passo n + 1 chesono gia state calcolate. Inoltre il metodo e incondizionatamente stabile.

Esempio 4.6.1 Applichiamo il metodo di Eulero esplicito al’equazione delcalore in due dimensioni

ut = uyy + uxx


definita sul quadrato [0, 1] × [0, 1] e con condizioni iniziali:

u(x, y, 0) =

300 (x− 0.5)2 + (y − 0.5)2 ≤ 0.2

0 altrimenti

u(x, y, t) = 0, x = 0, x = 1, 0 ≤ y ≤ 1, t ≥ 0y = 0, y = 1, 0 ≤ x ≤ 1, t ≥ 0

Nella Figure 4.6, 4.7 e 4.8 viene riportato l’andamento della condizione in-iziale e della soluzione numerica in una griglia composta da 100 punti circae in due istanti di tempo successivi (t pari a circa otto millesimi di secondo emezzo secondo rispettivamente). Si puo osservare che una condizione inizialediscontinua dia luogo ad una soluzione che invece e molto regolare. I valori∆x, ∆y e ∆t sono stati scelti in modo tale che venga soddisfatta la condizionedi stabilita (in particolare e stato posto ∆x = ∆y = 10−2 e ∆t = 2.5 · 10−5.

00.2

0.40.6

0.81

0

0.2

0.4

0.6

0.8

10

50

100

150

200

250

300

xy

Figura 4.6: Condizione iniziale per l’equazione del calore in due dimensioni.


00.2

0.40.6

0.81

0

0.2

0.4

0.6

0.8

10

50

100

150

200

250

300

xy

Figura 4.7: Soluzione numerica per l’equazione del calore dopo circa 8millisecondi.

00.2

0.40.6

0.81

0

0.2

0.4

0.6

0.8

10

50

100

150

200

250

300

xy

Figura 4.8: Soluzione numerica per l’equazione del calore dopo circa mezzosecondo.

Capitolo 5

Equazioni iperboliche

5.1 L’equazione d’onda

Le equazioni iperboliche rappresentano probabilmente la classe che descrive ilpiu ampio numero di fenomeni in diversi campi della fisica (fluidodinamica,acustica, elettromagnetismo e cosı via). Le equazioni di Eulero per fluidicomprimibili, le equazioni di Einstein per la relativita generale, sono esempidi equazioni iperboliche. L’equazione del secondo ordine di tipo iperbolicopiu nota e sicuramente l’equazione d’onda:

utt(x, t) − c2uxx(x, t) = 0, c 6= 0, (5.1)

infatti∆ = c2 > 0.

L’equazione (5.1) ammette una formulazione come equazione del primo or-dine:

ut(x, t) + cux(x, t) = 0, c ∈ R. (5.2)

Infatti derivando l’equazione (5.2) rispetto al tempo

∂2u

∂t2+ c

∂2u

∂t∂x= 0,

e rispetto allo spazio∂2u

∂t∂x+ c

∂2u

∂x2= 0, (5.3)

e, sostituendo la derivata mista in (5.3) si ricava appunto l’equazione (5.1). Inentrambe le equazioni c rappresenta una costante prefissata, cioe la velocita di

74

CAPITOLO 5. EQUAZIONI IPERBOLICHE 75

propagazione dell’onda (per un’onda sonora che si propaga nell’aria il valorec e pari a circa 340 m/sec, in genere il valore dipende dal tipo di onda). Lafunzione u(x, t) esprime l’ampiezza dell’onda, una misura della sua intensitain funzione della posizione x al tempo t. Per un’onda sonora nell’aria u(x, t)esprime la pressione dell’aria in diversi punti dello spazio, per una cordavibrante esprime lo spostamento fisico della corda rispetto alla posizione diriposo. Se il dominio spaziale e bidimensionale (l’onda si propaga nel pianoxy) l’equazione (5.1) puo essere conseguentemente modificata:

utt(x, t) = c2(uxx(x, t) + uyy(x, t)), (x, y) ∈ D ⊆ R2.

Un esempio di tale equazione si ha con il modo della pelle di un tamburocircolare tesa rigidamente. Ovviamente nulla vieta di inserire un’ulteriorevariabile spaziale per descrivere fenomeni quali la propagazione di onde inun mezzo omogeneo elastico, e il caso per esempio delle onde sismiche dellaTerra oppure le onde ultrasoniche usate per rilevare difetti nei materiali. Unamodifica dell’equazione d’onda consiste nel fatto che, piu realisticamente, lavelocita puo dipendere dall’ampiezza, il che conduce all’equazione non lineare

utt(x, t) = c(u)2(uxx(x, t) + uyy(x, t)), (x, y) ∈ D ⊆ R2.

All’equazione d’onda del secondo ordine (5.1) e possibile associare due tipiproblemi. Il primo e il il problema ai valori iniziali (o Problema di Cauchy)in cui si deve determinare una funzione u(x, t), definita e continua per x ∈ R

e t ≥ 0, che soddisfi l’equazione alle derivate parziali per x ∈ R e t > 0 e lecondizioni iniziali:

u(x, 0) = f1(x) x ∈ R

ut(x, 0) = f2(x) x ∈ R,

come schematizzato nella seguente figura.


t

xOu(x, 0) = f1(x) ut(x, 0) = f2(x)

Il secondo e il problema ai valori iniziali e al contorno, in cui, assegnatauna costante L > 0, si deve trovare una funzione u(x, t), definita e continuaper 0 ≤ x ≤ L e t ≥ 0, che soddisfi l’equazione alle derivate parziali per0 < x < L e t > 0 e le condizioni iniziali:

u(x, 0) = f1(x) 0 ≤ x ≤ Lut(x, 0) = f2(x) 0 ≤ x ≤ Lu(0, t) = g1(t) t ≥ 0u(L, t) = g2(t) t ≥ 0.

t

x

O

u(x, 0) = f1(x), ut(x, 0) = f2(x)

u(0, t) = g1(t)

u(L, t) = g2(t)

L


Le funzioni che definiscono le condizioni al contorno devono soddisfare lecondizioni di omogeneita agli angoli del dominio:

f1(0) = g1(0), f1(L) = g2(0), f2(0) = g′1(0), f2(L) = g′2(0).

t

x

O f1(x), f2(x)

g1(t) g2(t)

L

La risoluzione per via analitica del problema di Cauchy e possibile suppo-nendo per semplicita c = 1 ed effettuando il seguente cambio di variabile:

ξ = x+ t, ψ = x− t

ovvero

x =ξ + ψ

2, t =

ξ − ψ

2

e definendo la funzione

U(ξ, ψ) = u(x(ξ, ψ), t(ξ, ψ)) = u

(

1

2(ξ + ψ),

1

2(ξ − ψ)

)

.

Si osserva innanzitutto che∂2U

∂ξ∂ψ= 0. (5.4)

Infatti∂U

∂ξ=∂u

∂x

∂x

∂ξ+∂u

∂t

∂t

∂ξ=

1

2

∂u

∂x+

1

2

∂u

∂t


e, calcolando la derivata parziale seconda:

∂2U

∂ξ∂ψ=

1

2

∂

∂ψ

[

∂u

∂x+∂u

∂t

]

=

=1

2

[

∂x

∂ψ

(

∂2u

∂x2+

∂2u

∂x∂t

)

+∂t

∂ψ

(

∂2u

∂t∂x+∂2u

∂t2

)]

=

=1

2

[

1

2

∂2u

∂x2+

1

2

∂2u

∂t∂x−

1

2

∂2u

∂x∂t−

1

2

∂2u

∂t2

]

= 0.

L’uguaglianza a zero deriva dall’ipotesi che la funzione u(x, t) soddisfa l’e-quazione d’onda e dall’uguaglianza delle derivate parziali miste.Poiche Uξψ = 0 possiamo considerare la derivata Uξ come funzione della solavariabile ξ quindi integrando (5.4) rispetto a ψ si ottiene:

Uξ = F1(ξ)

e, integrando nuovamente rispetto a ξ:

U(ξ, ψ) =

∫ ξ

0

F1(z)dz +G2(ψ),

dove F1 e G2 sono due funzioni arbitarie differenziabili. Posto

G1(ξ) =

∫ ξ

0

F1(z)dz

risultaU(ξ, ψ) = G1(ξ) +G2(ψ).

Tornando alle variabili x e t si ha che la soluzione deve essere:

u(x, t) = G1(x+ t) +G2(x− t).

Sostituendo le condizioni iniziali risulta:

u(x, 0) = G1(x) +G2(x) = f1(x)ut(x, 0) = G′

1(x) −G′2(x) = f2(x).

e, differenziando la prima equazione:

G′1(x) +G′

2(x) = f ′1(x)


si ricava:

G′1(x) =

1

2[f ′

1(x) + f2(x)]

G′2(x) =

1

2[f ′

1(x) − f2(x)] ,

da cui, integrando rispetto a x, risulta:

G1(x) =1

2

[

f1(x) +

∫ x

0

f2(z)dz

]

G2(x) =1

2

[

f1(x) −

∫ x

0

f2(z)dz

]

.

Sostituendo tali formule nell’espressioni di u(x, t) si ottiene:

u(x, t) =1

2

[

f1(x+ t) +

∫ x+t

0

f2(z)dz

]

+1

2

[

f1(x− t) −

∫ x−t

0

f2(z)dz

]

=

=1

2

[

f1(x+ t) + f1(x− t) +

∫ x+t

x−t

f2(z)dz

]

che prende il nome di Formula di D’Alembert. Da tale formula segue chela funzione u(x, t) e determinata univocamente in base alla conoscenza dellefunzioni f1 ed f2 tra i punti (x − t, 0) e (x + t, 0). L’intervallo [x − t, x + t]viene detto intervallo di dipendenza del punto (x, t). Il dominio interno altriangolo di vertici (x, t), (x − t, 0) e (x + t, 0) ed evidenziata nella figuraseguente si chiama regione di dipendenza.


(x, t)

(x− t, 0) (x+ t, 0)

t

xO

Le rette congiungenti i punti (x, t) e (x − t, 0), (x, t) e (x + t, 0) sono dettecaratteristiche dell’equazione d’onda in (x, t).Osserviamo che, nel caso in cui il problema sia ai valori al contorno, la formuladi D’Alembert puo essere usata per calcolare la soluzione solo nel triangolodi vertici (0, 0), (L, 0) e (L/2, L/2).

(L/2, L/2)

(L, 0)

t

x

(0, 0)


Il metodo delle caratteristiche

Considerando, in generale, l’equazione iperbolica

∂u

∂t+ c

∂u

∂x= 0

le caratteristiche sono le rette di equazione

x = x0 ± ct

Considerando la caratteristica x(t) = x0 + ct, lungo questa la funzioneu(x(t), t) rimane costante, infatti calcolando la derivata rispetto a t,

du

dt(x(t), t) =

∂u

∂t+∂u

∂x

dx

dt= 0.

Volendo calcolare la soluzione in un punto (x, t) si individua l’equazione dellacaratteristica

x = x0 + ct

e, sapendo che

u(x, t) = u(x0 + ct, t) = u(x0, 0) = f1(x0),

poiche il valore della funzione lungo la caratteristica e indipendente dal tem-po. Se la condizione iniziale e discontinua (per esempio presenta un salto)allora anche la discontinuita si propaga alla medesima velocita dell’onda,come si evince dal seguente grafico.

x

ut

t = 0t = t1

t = t2

t = t3


5.2 L’equazione d’onda del primo ordine

I metodi che sono stati descritti nel Capitolo 4 possono essere applicati, in mo-do piu o meno naturale, anche alle equazioni d’onda. In particolare il Metododi Crank-Nicolson fornisce buoni risultati anche in questo caso, soprattuttonel caso di problema ai valori iniziali e al contorno. E opportuno evitaredi ripetere la descrizione di metodi gia discussi e considerare metodi diversiche possono essere applicati esclusivamente ad equazioni iperboliche, anchedi tipo non lineare, come sara descritto in seguito. Nei seguenti paragrafisaranno descritti alcuni metodi che possono essere applicati all’equazione delprimo ordine (5.2).

5.2.1 Il metodo di Lax-Friedrichs

Un primo modo per risolvere numericamente l’equazione d’onda del primoordine potrebbe essere quello di approssimare la derivata temporale con laformula alle differenze in avanti e la derivata spaziale con quella alle differenzecentrali, ottenendo la seguente formula:

un+1j = unj −

c∆t

2∆x(unj+1 − unj−1). (5.5)

Tale metodo risulta, in verita, piuttosto instabile (basta applicare l’analisi divon Neumann e verificare che il fattore di amplificazione e sempre, in modulo,maggiore di 1). Il metodo di Lax-Friedrichs consiste nello stabilizzare ilmetodo (5.5) sostituendo il valore l’approssimazione unj con il valor medio:

unj −→unj+1 + unj−1

2

ottenendo, appunto, l’espressione del metodo di Lax-Friedrichs:

un+1j =

1

2(unj+1 + unj−1) −

c∆t

2∆x(unj+1 − unj−1).

5.2.2 Il metodo di Lax-Wendroff

Il metodo di Lax-Wendroff, esplicito, risolve numericamente l’equazione (5.2)partendo dall’espansione in serie di Taylor della funzione u(xj, tn+∆t) rispet-to alla variabile temporale e prendendo (xj, tn) come punto iniziale:

u(xj, tn + ∆t) ≃ u(xj, tn) + ∆tut(xj, tn) +(∆t)2

2utt(xj, tn).


Applicando la relazione (5.2):

u(xj, tn + ∆t) ≃ u(xj, tn) − c∆tux(xj, tn) +(∆t)2

2utt(xj, tn).

Poiche la funzione u(x, t) soddisfa l’equazione (5.1) si ha

u(xj, tn + ∆t) ≃ u(xj, tn) − c∆tux(xj, tn) +(c∆t)2

2uxx(xj, tn). (5.6)

Le derivate spaziali vengono approssimate usando la consueta formula del sec-ondo ordine mentre per la derivata prima si applica la formula alle differenzecentrali.

un+1j = unj −

c∆t

2∆x

(

unj+1 − unj−1

)

+(c∆t)2

2(∆x)2

(


)

.

Posto

α =c∆t

∆x

si ottiene lo schema

un+1j = unj −

α

2

(

unj+1 − unj−1

)

+α2

2

(


)

=α

2(1 + α)unj−1 + (1 − α2)unj −

α

2(1 − α)unj+1.

5.2.3 Condizione di Courant, Friedrichs e Lewy

Consideriamo ora la risoluzione numerica dell’equazione d’onda (5.2) con-siderando per ora c = 1. L’approssimazione un+1

j dipende, sia nei metodiimpliciti che espliciti, da approssimazioni al livello precedente n, in parti-colare da unj , u

nj−1 e unj+1. A loro volta tali approssimazioni dipendono da

altre al livello n − 1, in particolare da un−1j±k , con k = −2, . . . , 2, cosı via. In

questo modo procedendo a ritroso e possibile definire una specie di dominiodi dipendenza discreto che composto dalle approssimazioni, dal livello 0 allivello n, necessarie al calcolo di un+1

j .


(xj, tn+1)

t

xO

Appare ovvio che tale dominio discreto debba avere necessariamente unlegame con quello continuo che abbiamo definito in precedenza. Se il do-minio continuo includesse quello discreto questo vorrebbe dire che l’approssi-mazione un+1

j e stata ottenuta considerando solo una parte dei valori da cuidipende il valore teorico u(xj, tn+1), sicuramente tale approssimazione nu-merica non puo essere un valore affidabile. Al contrario se il dominio discre-to contiene quello continuo significa che la soluzione numerica ha utilizzatoeffettivamente tutti i dati necessari (e anche altri).

(xj, tn+1)

t

xO

Dal punto di vista matematico si deve richiedere che tale situazione si verifi-chi, imponendo opportune condizioni sui passi di discretizzazione spaziale e


temporale. Infatti e necessario richiedere che la retta caratteristica passanteper (xj, tn+1) intersechi l’intervallo di dipendenza del metodo numerico.

(xj, tn+1)

h h(xj−1, tn) (xj+1, tn)

(xj, tn)

∆t

La condizione viene verificata se la retta di colore blu tratteggiata ha uncoefficiente angolare inferiore rispetto a quello della retta caratteristica (chein questo caso vale 1), cioe se

∆t

∆x≤ 1. (5.7)

La relazione (5.7) prende il nome di Condizione di Courant, Friedrichs e Lewy(spesso indicata come condizione CFL). Se la velocita c 6= 1 allora il rapporto∆t/∆x deve essere inferiore al coefficiente angolare delle rette caratteristiche,che hanno equazione

x = x0 + ct

e quindi∆t

∆x≤

1

c⇒

c∆t

∆x≤ 1.

5.3 L’equazione d’onda del secondo ordine

Nei successivi paragrafi saranno descritti alcuni metodi numerici per risolverel’equazione del secondo ordine (5.1) con velocita costante c 6= 0. Si supponeche siano assegnate le funzioni f1(x) ed f2(x) tali che

u(x, 0) = f1(x), ut(x, 0) = f2(x), 0 ≤ x ≤ L.


Un metodo esplicito

Innanzitutto si definisce la griglia suddividendo l’intervallo [0, L] in sottoint-ervalli di ampiezza ∆x = L/(N + 1) e definendo gli istanti di tempo multiplidi un valore ∆t:

xj = j∆x, j = 0, 1, 2, . . . , N + 1, tn = n∆t, n = 0, 1, 2, . . . .

Le derivate parziali seconde sono approssimate nel modo consueto:

uxx(xj, tn) ≃unj+1 − 2unj + unj−1

(∆x)2, utt(xj, tn) ≃

un+1j − 2unj + un−1

j

(∆t)2


j

(∆t)2− c2


(∆x)2= 0


j = c2(∆t)2

(∆x)2(unj+1 − 2unj + unj−1)

Posto

α =(c∆t)2

(∆x)2

e ricavato un+1j :

un+1j = 2unj − un−1

j + α(unj+1 − 2unj + unj−1)

un+1j = αunj−1 + 2(1 − α)unj + αunj+1 − un−1

j , j = 1, . . . , N.

I primi valori che e possibile calcolare sono u2j , che necessitano della conoscen-

za di u1j , poiche i valori u0

j sono forniti dalla condizione iniziale

u0j = u(xj, 0) = f1(xj).

Il problema e ora quello di approssimare la soluzione nei punti (xj,∆t), cioeconoscere i valori:

u1j ≃ u(xj,∆t), j = 1, . . . , N.

Per questo motivo si utilizza l’espansione in serie di Taylor:

u(xj,∆t) ≃ u(xj, 0) + ∆t ut(xj, 0) +(∆t)2

2utt(xj, 0). (5.8)


Poiche la funzione u(x, t) soddisfa l’equazione d’onda, allora possiamo sosti-tuire utt(x, t) con c2uxx(x, t), e le condizioni iniziali per u(x, t) e ut(x, t):

u(xj,∆t) ≃ f1(xj) + ∆tf2(xj) +(c∆t)2

2uxx(xj, 0).

L’ultimo termine della serie viene approssimato come al solito:

uxx(xj, 0) ≃u(xj+1, 0) − 2u(xj, 0) + u(xj−1, 0)

(∆x)2=

=f1(xj+1) − 2f1(xj) + f1(xj−1)

(∆x)2

cosicche si ottiene la seguente approssimazione:

ui,1 ≃ f1(xj) + ∆tf2(xj) +(c∆t)2

2

f1(xj+1) − 2f1(xj) + f1(xj−1)

(∆x)2.

Un metodo implicito

Per risolvere l’equazione d’onda si puo discretizzare in modo diverso la deriva-ta seconda di tipo spaziale:

uxx(xj, tn) ≃1

2[uxx(xj, tn+1) + uxx(xj, tn−1)]

utt(xj, tn) ≃un+1j − 2unj + un−1

j

(∆t)2

uxx(xj, tn+1) ≃un+1j+1 − 2un+1

j + un+1j−1

(∆x)2

uxx(xj, tn−1) ≃un−1j+1 − 2un−1

j + un−1j−1

(∆x)2

uxx(xj, tn) ≃1

2

[

un+1j+1 − 2un+1

j + un+1j−1

(∆x)2+un−1j+1 − 2un−1

j + un−1j−1

(∆x)2

]

.

Sostituendo le approssimazioni nell’equazione alle derivate parziali si ottiene:


j

(∆t)2= c2

[

un+1j+1 − 2un+1

j + un+1j−1

2(∆x)2+un−1j+1 − 2un−1

j + un−1j−1

2(∆x)2

]

.


Ponendo

α =(c∆t)2

2(∆x)2

si giunge all’espressione finale del metodo

−α un+1j−1 + (2α+ 1)un+1

j − αun+1j+1 = αun−1

j−1 + (2α− 1)un−1j + αun−1

j+1 + 2unj .

Se la soluzione numerica e nota ai livelli tn e tn−1 allora, utilizzando le con-dizioni al contorno, la formulazione del metodo costituisce un sistema linearedi N−1 equazioni nelle N−1 incognite un+1

j , j = 1, N−1. Il sistema ha unastruttura tridiagonale a predominanza diagonale quindi ammette un’unicasoluzione. Per calcolare la soluzione al livello t1 si puo utilizzare lo stessometodo di approssimazione (5.8) visto per il metodo esplicito.

5.4 Equazioni Iperboliche non Lineari

La risoluzione numerica di equazioni iperboliche non lineari e decisamentepiu complessa rispetto al caso lineare. Il problema e che spesso le non linear-ita inducono delle discontinuita nelle soluzioni teoriche anche in presenza dicondizioni iniziali regolari. In questo paragrafo considereremo soltanto unaclasse di equazioni iperboliche non lineari e descriveremo brevemente alcuneproblematiche connesse alla loro risoluzione numerica e contemporaneamentealcuni semplici metodi numerici senza avere la pretesa di esaurire un argo-mento che e molto complesso e che richiederebbe una trattazione molto piudettagliata. Consideriamo la seguente equazione, gia incontrata nel capitoloprecedente (equazione (4.1)),

ut(x, t) + [F (u)]x = 0 x ∈ R, t ≥ 0u(x, 0) = u0(x) x ∈ R

che viene detta in forma conservativa. Nell’ipotesi che la funzione F (u)sia sufficientemente regolare rispetto a u e a x allora l’equazione puo essereriscritta nella forma non conservativa:

ut(x, t) + F ′(u)ux = 0 x ∈ R, t ≥ 0u(x, 0) = u0(x) x ∈ R.

Spesso alla condizione iniziale si aggiunge una condizione di periodicita:

u(a, t) = u(b, t), t ≥ 0,


cosicche l’equazioni viene risolta quando x ∈ [a, b]. L’esempio piu noto diequazione di questo tipo e probabilmente l’equazione di Burgers:

F (u) =1

2u2

in cuiF ′(u) = u

e quindi l’equazione in forma non conservativa e

ut(x, t) + u(x, t)ux(x, t) = 0

cioe la velocita dell’onda coincide con il valore della funzione stessa. In questocaso le caratteristiche x(t) sono le soluzioni dell’equazione differenziale

x′(t) = u(x, t)

infatti

d

dtu(x(t), t) = ux(x, t)

dx

dt+ ut(x, t) = ux(x, t)u(x, t) + ut(x, t) = 0.

Poiche lungo le caratteristiche la funzione u(x, t) e costante allora questedevono essere necessariamente delle rette, infatti partendo dalla condizioneiniziale u0(x0) = u(x0, 0), esse si determinano sostituendo tale valore nell’e-quazione differenziale

x′(t) = u(x0, 0).

Le rette caratteristiche hanno equazione

x(t) = x0 + u(x0, 0)t = x0 + u0(x0)t

e coefficiente angolare 1/u0(x0), quindi non sono parallele, e possono in-tersecarsi. Nel punto in cui avviene tale intersezione la funzione u(x, t) ediscontinua, poiche presenta un salto. In questi casi si introduce il concet-to di soluzione debole dell’equazione differenziale per indicare una soluzioneche puo essere discontinua, oppure non differenziabile in qualche punto deldominio, in contrapposizione con le cosiddette soluzioni classiche. Se duecaratteristiche si intersecano in un punto verso la direzione crescente di t al-lora la discontinuita (o la soluzione debole) viene detta shock. Un altro casoin cui le caratteristiche possono dare luogo a discontinuita si verifica qualorale due rette partono dai due lati di una discontinuita e si allontanano da ques-ta. In questo caso la discontinuita prende il nome di onda di rarefazione. Lasituazione e schematizzata nel seguente grafico in cui le caratteristiche celestiformano un’onda di rarefazione, mentre quelle verdi formano uno shock.


x

uShockRarefazione

u0(x)

Per risolvere numericamente equazioni iperboliche non lineari, l’intervallo[a, b] di variabilita di x viene suddiviso in sottointervalli [xj, xj+1] di ugualeampiezza ∆x (nell’eventualita in cui l’intervallo sia infinito allora si fissanogli estremi in modo opportuno), si definisce una decomposizione detta dualedi intervalli:

Ij =]

xj− 1

2

, xj+ 1

2

[

, xj+ 1

2

= xj +∆x

2.

Integrando la forma conservativa dell’equazione iperbolica tra tn e tn+1 rispet-to al tempo e tra xj− 1

2

e xj+ 1

2

rispetto allo spazio si ottiene

∫ xj+1

2

xj− 1

2

dx

∫ tn+1

tn

[ut(x, t) + [F (u)]x] dt = 0

∫ xj+1

2

xj− 1

2

dx

∫ tn+1

tn

ut(x, t)dt+

∫ tn+1

tn

dt

∫ xj+1

2

xj− 1

2

[F (u)]xdt = 0

∫ xj+1

2

xj− 1

2

[u(x, tn+1) − u(x, tn)] dx+

+

∫ tn+1

tn

[

F(

u(

xj+ 1

2

, t))

− F(

u(

xj− 1

2

, t))]

dt = 0.


Definiamo l’approssimazione della media di u(x, tn) sull’intervallo Ij:

unj =1

∆x

∫ xj+1

2

xj− 1

2

u(x, tn)dx.

Sia

Hnj+1/2 =

1

∆t

∫ tn+1

tn

F(

u(

xj+ 1

2

, t))

dt

la media sull’intervallo temporale [tn, tn+1] del flusso in xj+ 1

2

. Posto

α =∆t

∆x

allora la relazione precedente puo essere riscritta come

un+1j = unj − α[Hn

j+1/2 −Hnj−1/2]

Da tale relazione si puo dedurre il seguente schema numerico

un+1j = unj − α[Φn

j+1/2 − Φnj−1/2]

in cui Φnj+1/2 = Φ(unj , u

n+1j ) e il cosiddetto flusso numerico che approssima

Hnj+1/2, mentre unj indica l’approssimazione del valore u(xj, tn). In generale

ogni schema numerico puo essere scritto sotto forma

un+1j = G(unj−l, . . . , u

nj−1, u

nj , u

nj+1, . . . , u

nj+l) (5.9)

in cui G e una determinata funzione.

Definizione 5.4.1 Uno schema numerico del tipo (5.9) e detto consistentese la funzione di flusso numerica coincide con il vero flusso F nel caso difunzioni costanti, cioe:

Φ(u, u) = F (u), ∀u.

Definizione 5.4.2 Si dice che uno schema numerico e detto monotono se lafunzione G aumenta in maniera monotona rispetto ad ognuna delle variabili.In questo caso un tale schema e limitato, cioe:

∃C > 0 : supj,n

|unj | ≤ C.


La condizione di Courant-Friedrichs-Lewy

Per l’equazione iperbolica non lineare la condizione di Courant-Friedrichs-Lewy diventa

|F ′(u)|∆t

∆x≤ 1,

disequazione che e sicuramente soddisfatta se

∆t

∆x≤

1

max |F ′(u)|. (5.10)

Il passo di discretizzazione temporale deve soddisfare la seguente disequazio-ne

∆t ≤∆x

max |F ′(u)|.

Il problema di tale vincolo e che e molto difficile stimare il massimo di F ′(u),funzione che spesso e non lineare. Poiche si puo ragionevolmente ipotizzareche la funzione soddisfa la proprieta di massimo-minimo (dipende dal feno-meno ondulatorio descritto) allora si puo stimare che il massimo di u(x, t)coincida con il massimo della condizione iniziale:

max |F ′(u)| = max |u(x, t)| = maxx

|u(x, 0)|.

La condizione (5.10) diventa quindi

∆t

∆x≤

1

maxx |u(x, 0)|. (5.11)

Lo Schema di Lax-Friedrichs

La funzione di flusso numerica ha la seguente espressione:

Φnj+1/2 =

1

2

[

F nj+1 − F n

j

]

−1

2α

[

unj+1 − unj]

con F nj = F (unj ).

Lo schema numerico si scrive quindi

un+1j =

1

2

[

unj+1 + unj−1

]

−α

2

[

F nj+1 − F n

j−1

]

.

Si tratta di uno schema del primo ordine che e monotono se

|F ′(unj )| ≤∆x

∆t=

1

α, ∀j, n.


Lo Schema di Godunov

La soluzione viene approssimata utilizzando una funzione costante a trat-ti. Il problema viene affrontato come se fosse un insieme di problemi concondizione iniziale che presenta diversi salti. In particolare si cerca, ad ungenerico istante di tempo tn la funzione u∗(x, t), con t ∈ [tn, tn+1], soluzionedel problema

ut + [F (u)]x = 0 x ∈ R, t ≥ 0u(x, tn) = unj xj−1/2 ≤ x ≤ xj+1/2.

Se ∆t e sufficientemente piccolo allora non ci sono interazioni tra questiproblemi locali. Si ha quindi

u∗(x, t) = u

(

x− xj+1/2

t− tn;unj , u

nj+1

)

dove x ∈ [xj, xj+1] e t ∈ [tn, tn+1].Il valore un+1

j e ottenuto come media di u∗(x, tn+1) su [xj−1/2, xj+1/2] :

un+1j =

1

∆x

∫ xj+1

2

xj− 1

2

u∗(x, tn+1)dx.

Lo schema numerico di Godunov diventa

un+1j = unj −

∆t

∆x

[

F (u(0;unj , un+1j ) − F (u(0;unj−1, u

nj )]

.

In generale solitamente lo schema di Godunov non si usa esattamente neitermini visti ma viene semplificato scegliendo il valore del flusso numericouguale a quello del flusso dell’equazione valutato in un’approssimazione in xjo xj+1. In questo caso si sceglie il flusso numerico ponendo:

Φnj+1/2 =

F (unj ) seF (unj+1) − F (unj )

unj+1 − unj≥ 0

F (unj+1) seF (unj+1) − F (unj )

unj+1 − unj< 0.

e quindi il metodo numerico ha la seguente espressione

un+1j = unj − α

[

Φnj+1/2 − Φn

j−1/2

]

.


Esempio 5.4.1 Consideriamo l’equazione di Burgers con condizione iniziale

u(x, 0) = 1 + 0.5 sin(2πx), 0 ≤ x ≤ 1

per t ∈ [0, 3]. Nelle figure seguenti sono riportate le soluzioni numericheottenute con i metodi di Lax-Friedrichs e di Godunov. Il metodo di Godunovconsente di evidenziare meglio la formazione di uno shock. I valori di ∆xe ∆t sono stati scelti in modo tale che venga soddisfatta la condizione diCourant, Friedrichs e Lewy (5.11), considerando che

maxx∈[0,1]

|u(x, 0)| = 1.5.

Figura 5.1: Soluzione delll’equazione di Burgers ottenuta con il metodo diLax-Friedrichs.


Figura 5.2: Soluzione dell’equazione di Burgers ottenuta con il metodo diGodunov.

Capitolo 6

Il Metodo agli Elementi Finiti

6.1 La formula di Green

Il metodo agli elementi finiti puo essere descritto efficacemente consideran-do la soluzione della semplice equazione ellittica di Poisson in un dominiobidimensionale. Sia assegnata quindi la seguente equazione

∂2u

∂x2+∂2u

∂y2= −f(x, y) ⇔ ∆u = −f

definita su un dominio Ω ⊂ R2 avente come frontiera la curva Γ, e su cui sono

assegnate le condizioni di Dirichlet alla frontiera. Uno strumento essenzialeper comprendere il metodo agli elementi finiti e la formula di Green. SiaΩ un insieme aperto avente come frontiera la curva chiusa e regolare Γ.Supponiamo inoltre che sia assegnata una funzione vettoriale

v(x, y) = [v1(x, y), v2(x, y)]T

continua e differenziabile in Ω. Il teorema della divergenza in spazi bidimen-sionali stabilisce che

∫

Ω

div v(x, y)dxdy =

∫

Γ

vT (x, y)nds,

ricordando che la divergenza della funzione vettoriale v e

div v(x, y) =∂v1

∂x+∂v2

∂y.

96

CAPITOLO 6. IL METODO AGLI ELEMENTI FINITI 97

Nell’integrale a secondo membro compare il prodotto scalare tra i vettori ve n che rappresenta la direzione normale alla curva Γ nel punto in questioneorientata secondo il verso uscente dalla curva stessa. Per derivare la formuladi Green si deve considerare una funzione scalare v ed una vettoriale w.Applicando le consuete regole di derivazione si ottiene

∇(vw) = (∇v)Tw + v∇w

che consente di esprimere (∇v)Tw come

(∇v)Tw = −v∇w + ∇(vw).

Integrando tale uguaglianza su Ω e applicando il teorema della divergenza siottiene

∫

Ω

(∇v)Tw dxdy = −

∫

Ω

v∇w dxdy +

∫

Ω

∇(vw) dxdy

= −

∫

Ω

v∇w dxdy +

∫

Γ

vwTn ds.

(6.1)

Tale uguaglianza puo essere considerata come una generalizzazione dellaformula di integrazione per parti. Il teorema di Green si ottiene da (6.1)sostituendo semplicemente la funzione vettoriale w con il gradiente di unafunzione scalare u(x, y) definita su Ω:

∫

Ω

∇vT∇u dxdy = −

∫

Ω

v∇(∇u) dxdy +

∫

Γ

v∇uTn ds.

Osservando che il gradiente del gradiente di u e il Laplaciano della funzione,e che la funzione ∇uTn, detta derivata normale, e indicata con

∇uTn =∂u

∂n,

allora si ottiene la formula di Green:∫

Ω

∇vT∇u dxdy = −

∫

Ω

v∆u dxdy +

∫

Γ

v∂u

∂nds.


6.2 Il metodo agli elementi finiti

Per risolvere numericamente il problema (cioe approssimare la funzione inco-gnita u) e necessario trasformare tale equazione in un sistema che possa essererisolto numericamente. I modi per approssimare u sono molti, comunqueil primo requisito e che tali approssimazioni appartengano ad uno spaziovettoriale di dimensione piccola. Inoltre ci sono altri requisiti. Per esempio, edifficile approssimare numericamente polinomi di grado elevato. Per ricavareil sistema di equazioni si usa comunemente la cosiddetta formulazione deboledel problema.Definiamo

a(u, v) ≡

∫

Ω

∇uT∇v dxdy =

∫

Ω

(

∂u

∂x

∂v

∂x+∂u

∂y

∂v

∂y

)

dxdy

e

(f, v) ≡

∫

Ω

f(x, y)v(x, y) dxdy.

Una proprieta immediata del funzionale a e che esso e bilineare, cioe e linearerispetto ad entrambi gli argomenti:

a(µu1 + µ2u2, v) = µ1a(u1, v) + µ2a(u2, v), ∀µ1, µ2 ∈ R

a(u, λ1v1 + λ2v2) = λ1a(u, v1) + λ2a(u, v2), ∀λ1, λ2 ∈ R.

Notiamo che (u, v) indica il prodotto scalare tra u e v nello spazio L2(Ω), cioel’insieme delle funzioni definite su Ω e aventi quadrato integrabile, infatti

(u, v) ≡

∫

Ω

u(x, y)v(x, y) dxdy.

Se v(x, y) e una funzione tale v(x, y) = 0 se (x, y) ∈ Γ, allora la formula diGreen applicata alla soluzione dell’equazione di Poisson diviene:

a(u, v) = −(∆u, v).

La formulazione debole del problema iniziale consiste nel selezionare unsottospazio V ⊂ L2 e quindi nel definire il seguente problema:

Trovare u ∈ V tale che a(u, v) = (f, v), per ogni v ∈ V .

Per capire quali possano essere le scelte per lo spazio V e bene osservare chela formulazione debole del problema coinvolge solo il prodotto scalare tra i


gradienti di u e v e quindi lo spazio V puo essere l’insieme delle funzioniderivabili e con derivata prima continua. Questo insieme e definito comeH1(Ω). Considerando anche le condizioni al contorno allora le funzioni inV devono essere nulle su Γ cosicche lo spazio vettoriale viene indicato conH1

0 (Ω).Il metodo agli elementi finiti consiste nell’approssimare il problema debolecon uno a dimensione finita ottenuto sostituendo lo spazio V con un sot-tospazio di funzioni che sono definite come polinomi di grado basso su piccolipezzi (detti appunto elementi) del dominio originale. Si considera ora laregione Ω del piano che viene suddivisa in triangoli (indipendentemente dal-la forma della curva che e la frontiera del dominio), come mostrato nellaseguente figura in cui un dominio a forma di ellisse e stato suddiviso in trian-goli. E chiaro che quando il numero di triangoli e molto elevato allora il loroinsieme ricoprira quasi perfettamente il dominio continuo qualunque formaquesto abbia.

•

•

•

•

•

•

•••

•

•

•

• • •

•

•

•

•

•

•

•

•

•

Quindi il dominio e approssimato dall’unione Ωh di m triangoli Ki:

Ωh =m⋃

i=1

Ki.

L’unica restrizione da porre e che nessun vertice di un triangolo appartengaal lato di un altro triangolo, cioe i triangoli possono condividere solo interilati e i relativi vertici. La dimensione della griglia h e definita come

h = maxi=1,...,m

diam(Ki)


dove diam(Ki), diametro del triangolo Ki, e la lunghezza del lato piu lungo.Lo spazio a dimensione finita Vh e definito come lo spazio di tutte le funzioniche sono lineari a tratti e continue nella regione Ωh e che sono zero sulcontorno Γ. Quindi

Vh =

ϕ∣

∣ ϕ|Ωhcontinua, ϕ|Γ = 0, ϕ|Kj

lineare per ogni j

,

dove ϕ|X indica la restrizione della funzione ϕ al sottoinsieme X. Se xj,j = 1, . . . , n, sono i nodi della triangolazione, allora una funzione ϕj ∈ Vhpuo essere associata ad ogni nodo, in modo tale che la famiglia di funzioniϕj(x) soddisfa le seguenti condizioni:

ϕj(xi) = δij =

1 se i = j0 se i 6= j.

Queste condizioni definiscono univocamente le funzioni ϕj e inoltre questeformano una base dello spazio Vh, quindi ogni funzione ϕ ∈ Vh puo essereespressa nella forma

ϕ(x) =n∑

i=1

ξiϕi(x).

L’approssimazione agli elementi finiti consiste nello scrivere la cosiddettacondizione di Galerkin per funzioni in Vh definendo cosı il seguente problema:

Trovare u ∈ Vh tale che a(u, v) = (f, v), per ogni v ∈ Vh.

Poiche u ∈ Vh ci sono n gradi di liberta. Poiche la funzione a(·, ·) e linearerispetto ai suoi argomenti allora e necessario imporre solo le condizioni:

a(u, ϕi) = (f, ϕi), i = 1, . . . , n.

Scrivendo u in termini della base ϕi come

u =n∑

i=1

ξiϕi(x).

e sostituendo nel problema di minimo si ottiene il sistema lineare

n∑

j=1

αijξj = βi


doveαij = a(ϕi, ϕj), βi = a(f, ϕi).

A questo punto si deve risolvere il sistema

Ax = b

in cui gli elementi della matrice A sono i valori αij, mentre quelli del vettore bsono i valori βj. Inoltre la matrice A e simmetrica e definita positiva, infattie evidente che

∫

Ω

∇ϕi∇ϕjdxdy =

∫

Ω

∇ϕj∇ϕidxdy

da cui segue che αij = αji. Per vedere che A e definita positiva prima no-tiamo che a(u, u) ≥ 0, per ogni funzione u. Se fosse a(ϕ, ϕ) = 0 per unafunzione appartenente a Vh, allora il gradiente di ϕ dovrebbe essere nulloquasi ovunque in Ωh. Poiche ϕ e lineare in ogni triangolo e continua in Ω,essa deve essere costante nello stesso insieme. Inoltre poiche e nulla sul con-torno dell’insieme deve essere identicamente nulla su Ω. Il risultato segueapplicando la relazione

(Aξ, ξ) = a(ϕ, ϕ)

con

ϕ =n∑

i=1

ξiϕi,

che e vera per ogni vettore ξ = [ξ1, . . . , ξn]T .

Un’importante osservazione e che la matrice A e molto sparsa. Infatti l’ele-mento αij e diverso da zero solo quando i supporti delle due funzioni dellabase (cioe i domini dove esse sono definite) hanno un triangolo in comune, o,equivalentemente, quando i nodi i e j sono vertici di uno stesso triangolo.Assegnato un nodo i il coefficiente αij e diverso da zero quando il nodo je uno dei vertici di un triangolo adiacente i. In pratica la matrice A vienecostruita sommando i contributi di tutti i triangoli applicando la formula

a(ϕi, ϕj) =∑

k

aK(ϕi, ϕj)

in cui la somma viene fatta su tutti i triangoli e

aK(ϕi, ϕj) =

∫

K

∇ϕi∇ϕjdx.


Notiamo che aK(ϕi, ϕj) e zero a meno che i nodi i e j sono entrambi verticidi K. Quindi un triangolo contrinuisce con valori diversi da zero ai 3 verticinella suddetta forma. La matrice 3 × 3:

AK =

aK(ϕi, ϕi) aK(ϕi, ϕj) aK(ϕi, ϕk)aK(ϕj, ϕi) aK(ϕj, ϕj) aK(ϕj, ϕk)aK(ϕk, ϕi) aK(ϕk, ϕj) aK(ϕk, ϕk)

associata al triangoloK(i, j, k) con vertici i, j, k e detta matrice degli elementidi stiffness. Per formare la matrice A e necessario sommare tutti i contributiaK(ϕk, ϕm) in posizione k,m della matrice. Questo procedimento viene dettoprocesso di assemblaggio:

A =N∑

κ=1

A[κ]

in cui N e il numero degli elementi e A[κ] e una matrice che ha solo 9 elementidiversi da zero pur essendo di dimensione uguale al numero dei nodi.

Esempio 6.2.1 Vediamo come esempio un dominio triangolare suddiviso inaltri 4 triangoli piu piccoli e che presenta 6 nodi. Nelle seguenti figure sonostati evidenziate le strutture delle matrici A[κ], con κ = 1, 2, 3, 4. Il numero dimatrici da assemblare coincide con il numero di elementi (triangoli), men-tre ciascuna di tali matrici presenta elementi diversi da zero solo in queglielementi di posto (i, j) tali che i nodi xi e xj sono vertici di tale triangolo.Nell’ultimo grafico viene evidenziata la struttura della matrice assemblata.Appare evidente come tale struttura dipende ovviamente dal modo con cuisono stati ordinati i nodi e che ci sono alcuni elementi che sono ottenuti dalcontributo di una singola matrice A[κ] (nel caso in cui un nodo e vertice diun solo triangolo), altri che invece hanno il contributo di tutte le matrici: ilnodo 5, per esempio, e vertice di tutti i triangoli.


6

4 5

1 2 3

1

2

43

Figura 6.1: Dominio per l’esempio 6.2.1.

6

4 5

1 2 3

1

2

43

A[1] =

Figura 6.2: Matrice A[1] relativa al triangolo 1.


6

4 5

1 2 3

2

1

43

A[2] =


6

4 5

1 2 3

1

43

2

A[3] =



6

4 5

1 2 3

1

4

2

3

A[4] =


6

4 5

1 2 3

A =

1

4

2

3

Figura 6.6: Struttura della matrice assemblata relativa al dominio a sinistra.


Condizioni al contorno particolarmente semplici non creano problemi peril metodo agli elementi finiti. Infatti il modo migliore per manipolare lecondizioni di Dirichlet e quello di includere i valori al contorno come incognitee modificare il sistema lineare incorporando tali valori. Cosı ogni equazioneassociata ad un punto sulla frontiera e sostituita da un’equazione del tipou(xi, yj) = f(xi, yj). Questo implica che all’interno del sistema alla matricedei coefficienti viene incorporata una piccola matrice identita. Nel caso dellecondizioni di Neumann la formula di Green fornisce le equazioni

∫

Ω

∇uT∇ϕjxdy = −

∫

Ω

fϕjdxdy +

∫

Γ

ϕj∂u

∂nds

che coinvolge appunto i valori della derivata normale forniti dalla condizionedi Neumann. Se tali valori sono noti solo in un insieme discreto di punti sullafrontiera allora si possono approssimare gli altri valori utilizzando l’interpo-lazione lineare (cioe approssimando l’integrale con la formula dei trapezi)o con il valor medio della funzione (cioe approssimando l’integrale con laformula del punto di mezzo). Se le condizioni al contorno sono tutte di Neu-mann allora e possibile che la matrice dei coefficienti del sistema possa esseresingolare, in tal caso si deve eliminare una delle equazioni del sistema oppurerisolvere il sistema tenendo conto della singolarita della matrice.

Generazione e raffinamento della griglia

Generare una triangolazione per gli elementi finiti puo essere fatto in modopiuttosto semplice partendo da una prima triangolazione e poi effettuandouna raffinamento in modo uniforme oppure in determinate aree del dominio.Il modo piu semplice di effettuare il raffinamento consiste nel prendere i trepunti medi dei lati di un triangolo e creando altri quattro triangoli che sosti-tuiscono quello di partenza. Ovviamente tale processo puo essere ripetuto an-che piu volte. Nella seguente figura si vede l’applicazione di tale raffinamentoal dominio visto nel precedente paragrafo.


6

4 5

1 27 8

9 10 11 12

13

14 15

3

1

2

3 4

5

6

8

12

13

7

910

14

11

15 16

In questo esempio particolare la matrice dei coefficienti, una volta riassem-blata, presenta la seguente struttura:

A =

Un vantaggio di questo raffinamento e che consente di preservare gli angolidella triangolazione originaria, e questo e importante perche gli angoli di una


buona triangolazione devono soddisfare determinati limiti. Tuttavia l’usoindiscriminato di strategie di raffinamento uniforme puo portare a qualcheproblema. Inoltre e preferibile raffinare la griglia nelle zone dove la soluzionedel problema presenta grandi variazioni, piuttosto che in quelle dove ha unandamento piuttosto uniforme. Anche l’introduzione dei nodi nei punti medidei lati andrebbe limitata a casi particolari, sarebbe opportuno aggiungerequalche punto e unirlo ai vertici dei triangoli gia creati.Il metodo degli elementi finiti, che abbiamo brevemente descritto in questocapitolo, solo nei tratti essenziali, ha dimostrato, negli ultimi anni, unaenorme versatilita, al punto da essere il metodo probabilmente piu utiliz-zato nella risoluzione numerica di equazioni alle derivate parziali. Il van-taggio indiscusso e la possibilita di applicazione a domini di qualsiasi tipo,sia in due che in tre dimensioni, problema che invece hanno i metodi alledifferenze finite. Un secondo, indiscutibile, vantaggio e l’uso di strategiedi griglia variabile, cioe la possibilita di utilizzare un numero di elementivariabile a seconda del comportamento della soluzione, che lo rende partico-larmente flessibile. Dall’altro lato va consiederata anche una certa difficoltanella necessita di riformulare il problema differenziale in termini variazionali.

Capitolo 7

Metodi iterativi per sistemisparsi

7.1 Il Metodo di Jacobi

Come abbiamo visto nel capitolo dedicato alle equazioni ellittiche i sistemilineari che derivano dalla discretizzazione di equazioni alle derivate parzialihanno grandi dimensioni e sono sparsi, cioe una buona parte degli elementidella matrice dei coefficienti sono nulli. Tuttavia quando si applicano i metodidiretti a tali sistemi succede che le matrici perdono la struttura di sparsita,cioe molti elementi nulli diventano diversi da zero e inoltre si ha il problemadi gestire matrici di grosse dimensioni, il che puo causare un notevole degradodelle prestazioni dei metodi usati. Per questi motivi si introduce una nuovaclasse di metodi, detti metodi iterativi. Supponiamo di dover risolvere ilsistema Ax = b, con A matrice non singolare, b 6= 0. Assumiamo inoltre chegli elementi diagonali della matrice aii, i = 1, . . . , n, siano diversi da 0. Lai−esima equazione del sistema si scrive

ai1x1 + ai2x2 + · · · + ainxn = bi

e, isolando xi risulta:

xi =

(

bi −

n∑

j=1,j 6=i

aijxj

)

1

aiii = 1, . . . , n.

Queste n equazioni sono del tutto equivalenti al sistema di partenza tut-tavia la loro forma suggerisce particolari procedimenti iterativi per cercare la

109

CAPITOLO 7. METODI ITERATIVI PER SISTEMI SPARSI 110

soluzione. A partire da un’approssimazione iniziale x(0) = (x(0)1 , x

(0)2 , . . . , x

(0)n )

si calcola la successione di vettori x(k) ponendo

x(k+1)i =

(

bi −n∑

j=1,j 6=i

aijx(k)j

)

1

aiii = 1, . . . , n; (7.1)

con k = 0, 1, 2, . . . .La generica componente i−esima del vettore al passo k + 1 e calcolata permezzo di tutte le componenti del vettore al passo k eccetto la i−esima.Questo procedimento iterativo prende il nome di metodo di Jacobi.

7.2 Il Metodo di Gauss-Seidel

Una variante del metodo di Jacobi si ottiene osservando che, quando si calcolax

(k+1)i si possono utilizzare le approssimazioni x

(k+1)j , con j = 1, . . . , i − 1,

ottenendo

x(k+1)i =

1

aii

(

bi −

i−1∑

j=1

aijx(k+1)j −

n∑

j=i+1

aijx(k)j

)

. (7.2)

Si ottiene in questo modo il classico Metodo di Gauss-Seidel, dove la compo-nente i−esima al passo k + 1 e calcolata per mezzo delle componenti dallaprima alla i−1-esima al passo k+1 e dalla i+1-esima alla n−esima al passok. Si deve osservare che entrambi i metodi appena introdotti non utilizzanola componente i−esima al passo k. Per questo si introduce una nuova vari-ante che coinvolge tale valore a partire da un parametro ω 6= 0. Si proponelo schema

x(k+1)i =

ω

aii

[

bi −

i−1∑

j=1

aijx(k+1)j −

n∑

j=i+1

aijx(k)j

]

+ (1 − ω)x(k)i . (7.3)

Questa classe di metodi prende il nome di Metodi di Rilassamento. Si osservafacilmente che se si pone ω = 1 il metodo di Rilassamento coincide con ilmetodo di Gauss-Seidel.Tutti i metodi appena descritti tengono in conto il carattere sparso dellamatrice dei coefficienti poiche si puo evitare di eseguire prodotti del tipo aijxjquando aij e nullo. Per decidere quando fermare il calcolo delle iterazioni


si puo pensare di fissare a priori una tolleranza ε e prendere x(k+1) comeapprossimazione della soluzione quando risulta

‖x(k+1) − x(k)‖

‖x(k+1)‖< ε

per una fissata norma vettoriale. Di solito si sceglie come approssimazioneiniziale x(0) il vettore nullo.

7.3 Il Metodo del Gradiente Coniugato

Consideriamo il problema di minimizzare la funzione

Φ(x) =1

2xTAx − xTb

dove b ∈ Rn e A ∈ R

n×n e una matrice simmetrica e definita positiva. Lasoluzione del problema e:

x = A−1b

quindi ha senso applicare i metodi di minimizzazione di Φ(x) per risolvere ilsistema lineare

Ax = b.

Un primo metodo per minimizzare Φ(x) e il metodo di steepest descent (tra-ducibile in italiano come metodo di piu ripida discesa). Esso si basa sul-l’osservazione che nel punto xc ∈ R

n la funzione Φ decresce piu rapidamentenella direzione del gradiente con il segno cambiato:

−∇Φ(xc) = b − Axc.

Poniamorc = b − Axc

il vettore detto residuo di xc. Se il residuo e diverso da zero cerchiamo ditrovare il parametro α > 0 tale che:

Φ(xc + αrc) < Φ(xc).


Vediamo ora come calcolare il valore di α.

Φ(xc + αrc) =1

2(xc + αrc)

TA(xc + αrc) − (xc + αrc)Tb

=1

2(xTc Axc + αrTc Axc + αxTc Arc + α2rTc Arc) − xTc b − αrTc b =

=1

2

[

α2rTc Arc + 2α(

rTc Axc − rTc b)]

+1

2xTc Axc − xTc b

=1

2

[

α2rTc Arc − 2αrTc rc]

+1

2xTc Axc − xTc b.

Dovendo minimizzare la funzione Φ(xc + αrc) calcoliamo il valore di α cheannulla la derivata prima:

Φ′(xc + αrc) = αrTc Arc − rTc rc.

Quindi

Φ′(xc + αrc) = 0 ⇒ α =rTc rcrTc Arc

.

Abbiamo quindi il seguente algoritmo:

k = 0x0 arbitario

r0 = b− Ax0

while rk 6= 0αk = (rTk rk)/(r

TkArk)

xk+1 = xk + αkrkrk = b − Axkk = k + 1

end

Il vettore xk rappresenta un’approssimazione della soluzione. La convergen-za del metodo dipende dallo spettro della matrice A ed in particolare dalrapporto tra gli autovalori estremi di A, cioe

λ1(A)

λn(A).

Quando tale rapporto e grande la convergenza e lenta. Inoltre quanto piuil vettore xk si avvicina alla soluzione tanto piu il residuo rk ha componenti


piccole e quindi la velocita di convergenza tende a diminuire.Una modifica del metodo consiste nel minimizzare la funzione Φ(x) lungouna direzione pk di descrescita della funzione, cioe tale che

pk∇φ(xk) < 0

determinando cosı l’approssimazione

xk+1 = xk + αkpk

tale cheφ(xk+1) = min

α∈R

Φ(xk + αkpk).

Ripetendo i passaggi in modo gia visto si ottiene

Φ′(xk + αpk) = (xk + αpk)TApk − bTpk

da cui, calcolando il punto stazionario, si ottiene

αk =(b − Axk)

TpkpTkApk

=rTkpk

pTkApk. (7.4)

Poiche rTkpk > 0 segue che αk > 0. Calcolando il residuo al passo k + 1 siottiene l’espressione:

b − Axk+1 = b − A(xk + αkpk) = b − Axk + αkApk

e quindirk+1 = rk − αkApk

e, dalla (7.4):

rTk+1pk = (rk − αkApk)Tpk = rTkpk − αkApTkpk = 0,

quindi ad ogni passo il residuo rk+1 e ortogonale al vettore direzione pk delpasso precedente. Il problema e la scelta delle direzioni pk. In particolaresi puo definire un metodo in cui viene scelta ad ogni passo una direzionepk lungo cui muovere le approssimazioni che tenga conto delle direzionip0,p1, . . . ,pk−1 calcolate ai passi precedenti. Un metodo che utilizza tale


strategia e il metodo del gradiente coniugato in cui le direzioni sono definitedai seguenti vettori:

pk =

r0 k = 0

rk + βkpk−1 k ≥ 1,

dove βk ∈ R viene scelto in modo tale che risulti:

pTkApk−1 = 0

cioe pk risulta A−coniugato rispetto al vettore pk−1. Imponendo tale propri-eta risulta:

(rk + βkpk−1)TApk−1 = 0

rTkApk−1 + βkpTk−1Apk−1 = 0

βk = −rTkApk−1

pTk−1Apk−1

, k ≥ 1. (7.5)

La direzione pk e una direzione di decrescita di Φ(x), infatti

−pTk∇Φ(xk) = pTk rk = rTk rk + βkpTk−1rk = rTk rk > 0.

Il valore di αk viene trovato nello stesso modo descritto nel metodo di piuripida discesa.

Φ(xk + αkpk) =1

2(xk + αkpk)

TA(xk + αkpk) − (xk + αkpk)Tb

=1

2(xTkAxk + αkp

TkAxk + αkx

TkApk + α2

kpTkApk) − xTkb − αkp

Tkb

=1

2

[

α2kp

TkApk + 2αk

(

pTkAxk − pTkb)]

+1

2xTkAxk − xTkb

=1

2

(

α2kp

TkApk − 2αkp

Tk rk)

+1

2xTkAxk − xTkb.

Il valore di αk che minimizza la funzione e

αk =pTk rk

pTkApk=

rTk rkpTkApk

.


OrarTk rk−1 = rTkpk−1 − βk−1r

Tkpk−2 = −βk−1r

Tkpk−2

erTkpk−2 = rTk−1pk−2 − αkp

Tk−1Apk−2 = 0

da cui segue cherTk rk−1 = 0, (7.6)

cioe ogni vettore residuo e ortogonale al precedente. Inoltre

pTk rk−1 = rTk rk−1 + βkpTk−1rk−1 = βkpk−1rk−1 = βkr

Tk−1rk−1 (7.7)

epTk rk−1 = pTk rk + αk−1p

TkApk−1 = pTk rk = rTk rk (7.8)

ottenendo una nuova espressione per βk:

βk =rTk rk

rTk−1rk−1

.

Possiamo riassumere il metodo del gradiente coniugato nel seguente algorit-mo:

k = 0x0 arbitario

r0 = b− Ax0

while rk 6= 0βk = (rTk rk)/(r

Tk−1rk−1) (β0 = 0 se k = 0)

pk = rk + βkpk−1 (p0 = r0 se k = 0)αk = (rTk rk)/(p

TkApk)

xk+1 = xk + αkpkrk+1 = rk − αkApkk = k + 1

end

Capitolo 1 Equazioni alle Derivate...

Documents

Transcript of Capitolo 1 Equazioni alle Derivate...