Introduzione ai metodi numerici per problemi differenziali ai limiti o al contorno

67 Pages·2017·0.4 MB·Italian

Introduzione ai metodi numerici per problemi
differenziali ai limiti o al contorno
Alessandra Sestini
20 novembre 2017
1 Introduzione
In queste note si indichera` con Ω un dominio in IRd (aperto, non vuoto) che
assumeremo limitato e connesso, mentre si indichera` con ∂Ω la sua frontiera.
Un generico problema differenziale al contorno pu`o allora essere formalmente
formulato come segue:
data la funzione f : Ω → IR, trovare u : Ω → IR tale che
Lu = f in Ω + boundary conditions su ∂Ω,
dove L indica un operatore differenziale e le condizioni al contorno spesso
sono di Dirichlet (u assegnata su ∂Ω) o di Neumann (derivata normale di u
assegnata su ∂Ω).
Imetodinumericisipongonol’obiettivodiapprossimarelasoluzioneudel
problema al contorno dato (che naturalmente deve essere ben posto) con una
sua approssimazione opportunamente (univocamente) definita e costruibile
con l’ausilio di un elaboratore elettronico. A questo scopo tutti utilizzano
una discretizzazione del problema differenziale che risulta definibile dopo che
`estatafissataunamesh, dove, genericamente, possiamopensareadunamesh
¯
come un ricoprimento di Ω (o anche di una sua porzione o estensione in certi
casi) fatto mediante politopi d–dimensionali (vedi sotto per ,d = 1,2). In
genere ad ogni mesh si associa un parametro di finezza che si indica con
h e che rappresenta il diametro massimo dei politopi che compongono la
mesh. Si noti che una condizione necessaria per l’utilizzo di un qualsiasi
metodo numerico `e che esso risulti convergente, ossia che l’approssimazione
1
da esso definita, se calcolata in aritmetica esatta, tenda alla soluzione esatta
u quando h tende a zero.
Alcuni metodi numerici si pongono l’obiettivo di approssimare u solo sui
vertici della mesh mentre altri (per esempio il metodo degli elementi finiti) ne
costruiscono un’approssimazione u in spazi funzionali a dimensione finita la
h
cui definizione dipende anche dalla mesh scelta. I metodi alle differenze finite
appartengono al primo gruppo e quindi costruiscono solo le approssimazioni
u ,i = 1,...,N (N = #verticidellamesh)diu(V )(doveV indical’i–esimo
i v v i i
vertice della mesh). Piu` precisamente essi definiscono il problema discreto
utilizzandodelleapprossimazionimediantedifferenzefinitedellederivatediu
che compaiono in L e eventualmente nelle condizioni al bordo. Considerando
qui solo problemi al contorno, tutti i valori u ,i = 1,...,N saranno calcolati
i v
simultaneamente risolvendo un sistema lineare o nonlineare (a seconda se il
problema differenziale `e lineare o no).
2 Problemi 1D
Quando Ω `e un intervallo I limitato della retta reale parleremo di problemi
differenziali ai limiti.
2.1 Un problema ai limiti modello
In questo caso senza perdita di generalita` possiamo assumere Ω = I := (0, 1)
e si considera il seguente problema differenziale,
−u(cid:48)(cid:48) + σu = f in I, (1)
combinatoconcondizioniassegnateagliestremi, dovef eσ sonoduefunzioni
assegnate, con σ ≥ 0. In particolare le condizioni agli estremi possono essere
di Dirichlet,
u(0) = g , u(1) = g , (2)
0 1
di Neumann
u(cid:48)(0) = g(cid:48) , u(cid:48)(1) = g(cid:48) . (3)
0 1
o miste,
u(0) = g , u(cid:48)(1) = g(cid:48) . (4)
0 1
2
Si noti che le ipotesi σ ≥ 0 e σ,f ∈ C0[0, 1], garantiscono l’unicit`a della
soluzione purch`e si evitino le condizioni di Neumann quando σ ≡ 0. Facen-
do esplicitamente riferimento al problema di Dirichlet, possiamo infatti far
vedere che, se z `e tale che −z(cid:48)(cid:48) + σz = 0 e z(0) = z(1) = 0, allora risulta
z ≡ 0. Questo si ottiene subito moltiplicando ambo i membri dell’equazione
differenziale per z stesso e integrando in I. Infatti, integrando per parti, si
ottiene che
(cid:90) 1 (cid:90) 1
(z(cid:48))2(x)dx+ σ(x)z2(x)dx = 0,
0 0
la quale, anche se fosse σ ≡ 0, considerando che deve essere z(0) = z(1) = 0,
implica che z ≡ 0. Quindi il problema omogeneo con condizioni di Dirichlet
omogenee ammette una e una sola soluzione, la soluzione nulla. Se ne pu`o
dedurre che anche il problema (1) non omogeneo con condizioni di Dirichlet
generiche ammette una e una sola soluzione quando σ e f sono entrambe
C0[0, 1]. Infatti, se u (x)+c z (x)+c z (x) indica l’integrale generale dell’
0 1 1 2 2
equazione differenziale lineare considerata (con u che indica una soluzione
0
particolareez ,z duesoluzionilinearmenteindipendentidell’equazioneomo-
1 2
genea), per trovare la soluzione cercata occorre determinare c ,c in modo
1 2
che
(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
z (0) z (0) c g −u (0)
1 2 1 = 0 0 .
z (1) z (1) c g −u (1)
1 2 2 1 0
Poich´e in particolare la soluzione del problema omogeneo con condizioni di
Dirichlet omogenee z `e la sola funzione identicamente nulla, necessariamen-
te la matrice dei coefficienti del precedente sistema lineare deve essere non
singolare. Da qui l’esistenza e unicit`a della soluzione.
Osserviamo poi che condizioni agli estremi non omogenee possono essere
sostituite dalle analoghe omogenee senza perdita di generalita`. Infatti basta
considerare il problema differenziale verificato da uˆ := u − R dove R `e
g g
un polinomio di grado 1 (2 per le Neumann pure) che dicesi rilevamento
dei dati sul bordo in quanto verifica le condizioni agli estremi, per esempio
R (x) := g (1−x)+g x nel caso di condizioni di Dirichelet e,
g 0 1
−uˆ(cid:48)(cid:48) +σuˆ = (f −σR ) in I
g
Perquantoriguardalaregolarita`, siosserviche, sef ,σ ∈ Ck[0, 1],k ≥ 0,
allora u ∈ Ck+2([0, 1]).
3
2.2 Metodo alle differenze finite
Definiamo uno schema alle differenze finite innanzi tutto fissando una mesh
su I, ossia un ricoprimento di [0, 1] mediante i segmenti definiti da una
partizione 0 = x < ··· < x = 1. Per semplicita` supponiamo che tale
0 N
partizione sia uniforme, ossia x = jh, j = 0,...,N, dove h = 1/N. Osser-
j
viamo che, supponendo che u ∈ C4[0, 1], dagli sviluppi di Taylor in x di
j
u(x ) = u(x ±h) fino al terzo ordine, si ottiene che
j±1 j
u(x )−2u(x )+u(x ) 1
u(cid:48)(cid:48)(x ) = j+1 j j−1 − u(4)(ξ )h2,
j h2 12 j
doveξ `eunpuntoopportuno in(x , x ).Quindi, considerandoilproble-
j j−1 j+1
ma modello con condizoni di Dirichlet, per j = 1,...,N −1, si pu`o scrivere
che
−u(x )+2u(x )−u(x ) 1
j+1 j j−1 + u(4)(ξ )h2 + σ(x )u(x ) = f(x ).
h2 12 j j j j
Se introduciamo la notazione τ = 1 u(4)(ξ )h2 (errore di troncamento
j 12 j
locale) e poniamo
u = (u(x ),··· ,u(x ))T,
1 N−1
τ = (τ ,··· ,τ )T,
h 1 N−1
b = (f(x )+g /h2,f(x ),··· ,f(x ),f(x )+g /h2)T ,
h 1 0 2 N−2 N−1 1
possiamo allora scrivere in forma matriciale,
A u = b +τ , (5)
h h h
dove A = 1 tridiag(−1,2,−1)+diag(σ ,··· ,σ ) e σ = σ(x ). Il metodo
h h2 1 N−1 j j
alle differenze finite consiste allora nel determinare un’approssimazione u di
h
u andando a risolvere il sistema lineare
A u = b . (6)
h h h
Osserviamo che u risulta ben definito in quanto A `e una matrice non sin-
h h
golare, essendo sdp. Infatti la sua simmetria `e evidente e si pu`o verificare
che, ∀y ∈ IRN−1 non nullo risulta yTA y > 0. Infatti si ha che
h
N−1 N−1
(cid:88) (cid:88)
yTA y = σ y2 +[y2 +y2 + (y −y )2]/h2.
h j j 1 N−1 j j−1
j=1 j=2
4
Risultando che τ tende a zero quando h tende a zero, il metodo dicesi
h
consistente. In particolare nel nostro caso risulta τ = O(h2). Tuttavia la
h
consistenza non assicura da sola la convergenza del metodo. Per studiarne la
convergenza dobbiamo considerare il comportamento dell’errore e = u −u
h h
quando h tende a zero. Dato che risulta A e = τ , e quindi e = A−1τ ,
h h h h h h
possiamo scrivere
(cid:107)e (cid:107) ≤ (cid:107)A−1(cid:107)(cid:107)τ (cid:107).
h h h
Vogliamo allora far vedere che, lavorando in norma infinito, siamo in grado di
trovare una costante che, per ogni h, maggiora (cid:107)A−1(cid:107). A questo scopo osser-
h
viamochesipu`odimostrarechesiaA chelamatriceA = 1 tridiag(−1,2,−1)
h 0h h2
hanno inversa non negativa (vedi Appendice C) e si ha
A−1 −A−1 = A−1(A −A )A−1 ≥ 0,
0h h 0h h 0h h
per le ipotesi su σ. Questo implica naturalmente che (cid:107)A−1(cid:107) ≤ (cid:107)A−1(cid:107) . Ora
h ∞ 0h ∞
osserviamo che(cid:107)A−1(cid:107) = max (A−1e) , dove e indica il vettore di tutti 1.
0h ∞ j 0h j
Osservando che la soluzione esatta del problema −u(cid:48)(cid:48) = 1, u(0) = u(1) = 0,
`e il polinomio di secondo grado φ(x) = x(1 − x)/2, possiamo concludere
che (A−1e) = φ(x ) e quindi che (cid:107)A−1(cid:107) ≤ (cid:107)A−1(cid:107) ≤ max |φ(x)|.
0h j j h ∞ 0h ∞ 0≤x≤1
Questo risultato di stabilit`a ci permette di concludere che l’errore e ha lo
h
stesso ordine dell’errore di troncamento τ e di conseguenza che il metodo `e
h
convergente del secondo ordine.
Sinotichel’uniformelimitatezzadellanormadiA−1 implicacheilmetodo
h
numerico sia stabile come vedremo nel caso generale nonlineare.
2.3 Metodo delle differenze finite per un problema ai
limiti piu` generale
Per problemi ai limiti piu` generali, se si considerano condizioni agli estremi di
Dirichlet, per quanto riguarda l’esistenza e l’unicita` della soluzione si pu`o far
ricorso al teorema sotto riportato che si riferisce alla seguente formulazione
generale di un problema ai limiti del secondo ordine (detto talvolta problema
dei due punti):
(cid:26)
u(cid:48)(cid:48) = F(x,u,u(cid:48)), x ∈ I
(7)
u(0) = g , u(1) = g .
0 1
5
Teorema 1. Sia F(x,y,z) : D → IR con D := I×(−∞,+∞)2 una funzione
continua con derivate parziali F e F continue in D e tale che esistono finiti
y z
max |F | e max |F |.
y z
(x,y,z)∈D (x,y,z)∈D
Se risulta
F ≥ q > 0 ∀(x,y,z) ∈ D,
y
allora il problema (7) ammette una e una sola soluzione.
Consideriamo per esempio il seguente problema lineare, detto problema di
diffusione–trasporto–reazione che per brevita` nel seguito indicheremo come
problema ADR (acronimo dell’inglese Advection–Diffusion–Reaction),
−u(cid:48)(cid:48) + γu(cid:48) + σu = f in I, (8)
dove si assume sempre che i dati γ,σ e f siano funzioni continue in [0, 1] e
dove al solito chiediamo σ ≥ 0 in [0, 1]. Riscritto come in (7), si ha che
F(x,y,z) = σ(x)y +γ(x)z −f(x).
Se quindi l’equazione differenziale `e combinata con condizioni ai limiti di
Dirichlet e si assume σ > 0 in [0, 1], siamo nelle ipotesi del precedente
teorema e quindi `e garantita l’esistenza e unicita` della soluzione (soluzione
forte). Vedremopiu` avantiche, passando allaformulazionedeboledelproble-
ma differenziale, sar`a possibile tornare solo a chiedere σ ≥ 0 per dimostrare
l’esistenza e unicit`a della soluzione in senso debole. Osserviamo inoltre che
se γ `e una costante, procedendo come fatto per il problema modello, si puo`
dimostrare l’esistenza e unicit`a della soluzione forte anche se si assume solo
σ ≥ 0.
Vediamo quindi come si puo` risolvere il problema ADR numericamente,
avvalendoci del metodo delle differenze finite. Approssimando la derivata
prima in x col la differenza finita centrata (molecola a tre punti), se u si
j
suppone di classe C3 si ha che
u(x )−u(x ) 1
u(cid:48)(x ) = j+1 j−1 + u(3)(η )h2.
j j
2h 6
Di conseguenza si ottiene che vale sempre il sistema in (5), dove ora A `e
h
sempre tridiagonale ma tale che, ∀1 ≤ i ≤ N −1 risulta
1 hγ 1 1 hγ
A (i,i−1) = − (1+ i), A (i,i) = (2+h2σ ), A (i,i+1) = − (1− i),
h h2 2 h h2 i h h2 2
6
e si ha,
(cid:18) 1 hγ 1 hγ (cid:19)T
1 N−1
b = f + (1+ )g , f ,··· ,f ,f + (1− )g .
h 1 h2 2 0 2 N−2 N−1 h2 2 1
Per quanto riguarda l’errore di troncamento locale si ha che ora
1 1
τ = u(4)(ξ )h2 + γ u(3)(η )h2
j j j j
12 6
che quindi implica che τ `e sempre del secondo ordine.
h
Al solito u sara` determinata andando a risolvere (6).
h
Teorema 2. Se σ,γ,f ∈ C0[0, 1], con σ ≥ σ > 0, allora ∀h > 0la
min
matrice tridiagonale A risulta a predominanza diagonale in senso forte (e
h
quindi nonsingolare) purch´e risulti
hγ (h)
max
≤ 1, (9)
2
dove γ (h) := max |γ |.
max i
i=1,...,N−1
Dimostrazione : Osserviamo che nell’ipotesi fatta risulta 1−1hγ (h) ≥ 0
2 max
e quindi risulta anche che ∀j si ha 1 ± 1hγ ≥ 0. Ne consegue che tutti gli
2 j
elementi extradiagonali di A sono non positivi e su tutte le righe interne
h
la somma dei loro valori assoluti `e pari a 2 e quindi strettamente minore
h2
del corrispondente elemento diagonale, essendo per ipotesi σ ≥ σ > 0.
min
Questo a maggior ragione sar`a vero sulla prima e ultima riga. Quindi A `e a
h
predominanza diagonale in senso forte.
Il seguente teorema ci assicura la stabilit`a del metodo e quindi la sua
convergenza del secondo ordine se la soluzione u `e di classe C4.
Teorema 3. Se σ,γ e f sono continue in [0, 1] e risulta
• σ ≥ σ > 0,
min
• vale la disuguaglianza in (9)
si ha anche che
1
max |u (j)−u(x )| ≤ (cid:107)τ (cid:107) .
h j h ∞
1≤j≤N−1 σmin
7
Dimostrazione : Poniamo e := u (j)−u(x ) e supponiamo che sia |e | ≥
j h j k
|e |,j = 1,...,N −1. Supponiamo inoltre che sia 1 < k < N −1 (altrimenti
j
si pu`o adattare il ragionamento). Poich´e risulta che A e = τ , tenendo
h h h
presente i segni degli elementi non nulli di A , si ha che
h
(2+σ h2)|e | ≤ (1+ 1hγ )|e | + (1− 1hγ )||e |+h2|τ |
k k 2 k k−1 2 k k+1 k
≤ 2|e |+h2|τ |
k k
Questo quindi implica che
σ h2|e | ≤ h2|τ |
k k k
e quindi la tesi.
Osservazione 1. Dal precedente teorema possiamo dedurre che il metodo
delle differenze finite (con derivata prima approssimata con formula alle dif-
ferenze finite centrate) risulta solo condizionatamente stabile per il problema
ADR.
Osserviamoche, nell’ipotesiγ ≥ 0inI,ilvincolodistabilit`apuo`essereri-
mosso se nel metodo alle differenze finite si utilizza un rapporto incrementale
all’indietro (schema upwind) invece che centrato per approssimare la derivata
prima, al prezzo di avere un errore locale di troncamento τ dell’ordine di h
h
invece che di h2. In tal caso la matrice A risulta essere
h
1 1 1
A (i,i−1) = − (1+hγ ), A (i,i) = (2+hγ +h2σ ), A (i,i+1) = − ,
h h2 i h h2 i i h h2
che `e a diagonale dominante in senso forte. Il fatto che in questo caso il
metodo risulti incondizionatamente stabile possiamo dimostrarlo osservando
che si pu`o scrivere
u −u u −u −u +2u −u
i i−1 i+1 i−1 i−1 i i+1
= +µ ,
h 2h h h2
con
h
µ := , viscosita` artificiale.
h
2
Allora la i–esima istanza dell’equazione alle differenze si pu`o riscrivere come
segue,
−u +2u −u u −u
i−1 i i+1 i+1 i−1
(1+γ µ ) +γ +σ u = f .
i h h2 i 2h i i i
8
Posto
ˆ
γˆ := γ /(1+γ µ ), σˆ := σ /(1+γ µ ), f := f /(1+γ µ )
i i i h i i i h i i i h
dal teorema precedentemente dimostrato si ha che la condizione di stabilita``e
hγˆ (h) ≤ 2, dove γˆ (h) = maxγˆ . Ricordando che µ = h/2, si verifica
max max i h
che tale condizione `e verificata per ogni passo h.
2.4 Il problema di diffusione–trasporto a coefficienti
costanti
Consideriamo il seguente problema di diffusione–trasporto a coefficienti co-
stanti:
−µu(cid:48)(cid:48) +γu(cid:48) = 0, x ∈ I, u(0) = 0, u(1) = 1,
con µ e γ costanti positive rispettivamente chiamate coefficiente di viscosit`a
e di trasporto. Si osservi che in questo caso `e possibile calcolare facilmente
la soluzione esatta del problema considerato. Infatti si ha che l’integrale
generale dell’equazione differenziale lineare omogenea a coefficienti costanti
diordine2consideratarisultaesserec +c exp(γ x) equindilasoluzioneu(x)
1 2 µ
si trova determinando le costanti c e c in modo che valgano le condizioni
1 2
ai limiti, ottenendo nel nostro caso la seguente finzione crescente di tipo
esponenziale,
γ γ
u(x) = (exp( x)−1)/(exp( )−1).
µ µ
Se γ `e molto grande tale soluzione presenta un boundary layer di ampiezza
µ
O(µ) in corrispondenza di x = 1 e in tale regione la derivata assume valori
γ
prossimi a γ. Se si usa il metodo DFC (differenze finite centrate) per appros-
µ
simare la soluzione di questo problema si ottiene la seguente equazione alle
differenze lineare:
hγ hγ
−(µ+ )u + 2µu − (µ− )u = 0, i = 1,...,N −1.
i−1 i i+1
2 2
Dividendo per µ, posto
hγ
Pe := , numero di Peclet locale (10)
2µ
si pu`o piu` brevemente riscrivere
−(1+Pe)u +2u −(1−Pe)u = 0, i = 1,...,N −1.
i−1 i i+1
9
In questo caso siamo quindi in grado di scrivere a priori l’espressione ana-
litica della soluzione del problema discreto. Infatti la soluzione generale
dell’equazione alle differenze considerata `e la seguente,
u = a λi +a λi ,
i 1 1 2 2
dove λ e λ sono le radici del polinomio (Pe−1)λ2+2λ−(Pe+1) e quindi
1 2
(cid:112)
−1± 1+(Pe2 −1)
λ =
1,2
Pe−1
e quindi λ = 1, λ = 1+Pe . La soluzione del nostro problema discreto
1 2 1−Pe
si trova determinando a e a in modo che valgano lo condizioni al limiti
1 2
assegnate, ottenendo le equazioni
a +a = 0, a +a λN = 1
1 2 1 2 2
e quindi ponendo −a = a = 1/(λN −1). In definitiva si ha
1 2 2
λi −1
u = 2 , i = 0,...,N .
i λN −1
2
La soluzione discreta trovata `e crescente in i solo se Pe < 1 in quanto solo
sotto questa condizione risulta λ > 1. Se invece si ha Pe > 1 risulta λ < −1
2 2
e quindi la soluzione discreta presenter`a delle oscillazioni intorno allo zero
tanto piu` marcate quanto piu` ci si avvicina al boundary layer. Osserviamo
che la condizione Pe < 1, `e proprio una istanza della condizione di stabilita`
introdotta per il problema ADR generale. Notiamo infine che comunque la
convergenza per h → 0 `e assicurata (tenendo presente che N = 1/h e come `e
definito Pe, verificare che lim λN = exp(γ)). Tuttavia l’andamento della
h→0 2 µ
soluzione numerica non `e affatto conforme a quello della soluzione continua
se Pe > 1.
2.5 Metodo delle differenze finite per problemi ai limiti
non lineari
Per il problema (7) il metodo DFC porta alla risoluzione di un sistema non
lineare di N −1 equazioni in altrettante incognite la cui i–esima equazione `e
−u +2u −u u −u
i−1 i i+1 i+1 i−1
+ F(x , u , ) = 0.
h2 i i 2h
10