Capitolo 4 Approssimazione di PDE con il metodo della

Capitolo 4
Approssimazione di PDE con il metodo
della decomposizione di domini (DD)
Prof. Alfio Quarteroni
4.1 Introduzione
Sia Ω un dominio di dimensione d, per d = 2, 3, con frontiera Lipschitziana ∂Ω; si consideri su
Ω il seguente problema differenziale:

Lu = f in Ω,




u=ϕ
su ΓD ,
(4.1)


∂u


= ψ su ΓN ,
∂n
dove L è un operatore ellittico su Ω (ad esempio, il Laplaciano o l’operatore di diffusione e trasporto), mentre ϕ, ψ sono due funzioni assegnate su ΓD e ΓN , rispettivamente, con ΓD ∪ ΓN =
∂Ω.
Vogliamo calcolare una soluzione per il problema modello (4.1), utilizzando uno schema numerico efficiente. A tal fine, come primo passo, introduciamone un’approssimazione alla Galerkin.
Utilizzeremo quindi le tecniche della decomposizione di domini (DD) per calcolare la soluzione
di tale problema approssimato.
L’idea alla base dei metodi DD è quella di suddividere il dominio globale Ω in due o più sottodomini su cui risolvere dei problemi di dimensione minore rispetto a quello iniziale, utilizzando
possibilmente degli algoritmi paralleli. In particolare, esistono due modi differenti con cui fare
una decomposizione del dominio Ω: si possono cioè considerare dei sotto-domini con o senza
sovrapposizione (si veda la Figura 4.1). Tale scelta individuerà metodi differenti per la risoluzione del problema assegnato.
95
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
96
Ω1
Γ
Ω
Ω2
PSfrag replacements
Ω1
Γ
2
Γ
Γ12
Ω
Ω2
1
Figura 4.1: Esempi di partizione del dominio Ω senza e con sovrapposizione.
Come lettura di riferimento per le tecniche di decomposizione di domini, rimandiamo, ad esempio, a [30, 41, 36, 25].
4.2 Tre classici metodi iterativi basati su DD
Introduciamo in questa sezione tre diversi schemi iterativi partendo dal problema modello (4.1)
fatte le scelte ϕ = 0 e ΓN = ∅: trovare u : Ω → R tale che
(
Lu = f in Ω,
u=0
su ∂Ω,
(4.2)
con L operatore ellittico generico del second’ordine.
4.2.1 Il metodo di Schwarz
Consideriamo una decomposizione del dominio Ω in due sotto-domini Ω 1 e Ω2 tali che Ω =
Ω1 ∪ Ω2 , Ω1 ∩ Ω2 = Γ12 6= ∅ (si veda la Figura 4.1) e sia Γi = ∂Ωi \ (∂Ω ∩ ∂Ωi ).
(0)
Consideriamo il seguente metodo iterativo: dato u2 su Γ1 , si risolvano i seguenti problemi per
k ≥ 1:

(k)
Lu1 = f
in Ω1 ,



(k)
(k−1)
(4.3)
u1 = u 2
su Γ1 ,


 (k)
u1 = 0
su ∂Ω1 \ Γ1 ,
PSfrag replacements
4.2. TRE CLASSICI METODI ITERATIVI BASATI SU DD
97
Ω2
γ2
a
γ1
b
Ω1
Figura 4.2: Esempio di decomposizione con sovrapposizione in dimensione 1.

(k)
Lu2 = f



(k)



(k)
u2 = u 1
(k)
u2 = 0
in Ω2 ,
(4.4)
su Γ2 ,
su ∂Ω2 \ Γ2 .
Abbiamo così due problemi ellittici con condizioni al bordo di Dirichlet per i due sotto-domini
(k)
(k)
Ω1 e Ω2 , e vogliamo che le due successioni {u1 } e {u2 } tendano alle rispettive restrizioni della
soluzione u del problema (4.2) quando k tende a infinito:
(k)
lim u1 = u|Ω1 e
k→∞
(k)
lim u2 = u|Ω2 .
k→∞
Si può dimostrare che il metodo di Schwarz applicato al problema (4.2) converge sempre alla
soluzione del problema di partenza, con una velocità che aumenta all’aumentare della misura
|Γ12 | di Γ12 . Mostriamo questo risultato in un caso semplice monodimensionale.
Esempio 4.1. Sia Ω = (a, b) e siano γ1 , γ2 ∈ (a, b) tali che a < γ2 < γ1 < b (si veda la Figura
4.2). I due problemi (4.3) e (4.4) diventano

(k)
Lu1 = f
a < x < γ1 ,



(k)
(k−1)
(4.5)
u1 = u 2
x = γ1 ,


 (k)
u1 = 0
x = a,

(k)
Lu2 = f γ2 < x < b,



(k)
(k)
u2 = u 1 x = γ 2 ,


 (k)
u2 = 0
x = b.
Per dimostrare la convergenza di tale schema, consideriamo il problema semplificato
(
−u00 (x) = 0
a < x < b,
u(a) = u(b) = 0,
(4.6)
(4.7)
ovvero il problema modello (4.2) con L = −d2 /dx2 e f = 0, la cui soluzione è evidentemente
u = 0 in (a, b).
(1)
(1)
Sia k = 1; dal momento che (u1 )00 = 0, u1 (x) è una funzione lineare e, in particolare, coincide
(0)
con la retta che assume il valore zero in x = a e u2 in x = γ1 . Conoscendo dunque il valore
PSfrag replacements
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
98
(0)
u2
(1)
u1
(2)
u1
a
γ2
(1)
u2
(2)
u2
γ1
b
Figura 4.3: Alcune iterazioni del metodo di Schwarz per il problema (4.7).
(1)
di u1 in γ2 si può risolvere il problema (4.6) che, a sua volta, è caratterizzato da una soluzione
lineare. Si procede poi in modo analogo. Mostriamo in Figura 4.3 alcune iterazioni: si vede
chiaramente che il metodo converge e che la velocità di convergenza si riduce al ridursi della
lunghezza dell’intervallo (γ2 , γ1 ).
•
Osserviamo che il metodo iterativo di Schwarz (4.3)-(4.4) richiede, ad ogni passo, la risoluzione
di due sotto-problemi con condizioni al bordo dello stesso tipo di quelle del problema di partenza: infatti si parte da un problema con condizioni di Dirichlet omogenee in Ω e si risolvono due
problemi con condizioni al bordo ancora di Dirichlet su Ω1 e Ω2 .
Se il problema differenziale (4.2) fosse stato completato invece da una condizione di tipo Neumann su tutto il bordo ∂Ω, ci si sarebbe ricondotti alla risoluzione di un problema misto (di tipo
Dirichlet-Neumann) su ciascuno dei sotto-domini Ω1 e Ω2 .
4.2.2 Il metodo di Dirichlet-Neumann
Decomponiamo ora il dominio Ω in due sotto-domini senza sovrapposizione (si veda la Figura
4.1): siano dunque Ω1 , Ω2 ⊂ Ω, tali che Ω1 ∪ Ω2 = Ω, Ω1 ∩ Ω2 = Γ e Ω1 ∩ Ω2 = ∅. Nel seguito
indicheremo con ni la normale esterna al dominio Ωi , utilizzando la seguente convenzione n =
n1 = −n2 .
Si può dimostrare il seguente risultato (si veda [30]):
Teorema 4.1 (d’equivalenza). La soluzione u del problema (4.2) è tale che u |Ωi = ui per i =
1, 2, dove ui è la soluzione del problema
(
Lui = f in Ωi ,
(4.8)
ui = 0
su ∂Ωi \ Γ,
con condizioni di interfaccia
e
u1 = u 2
(4.9)
∂u2
∂u1
=
∂nL
∂nL
(4.10)
su Γ, avendo indicato con ∂/∂nL la derivata conormale.
4.2. TRE CLASSICI METODI ITERATIVI BASATI SU DD
99
Grazie a tale risultato di equivalenza, si può decomporre il problema (4.2) utilizzando le condizioni di accoppiamento (4.9)-(4.10) come condizioni al bordo sull’interfaccia Γ. In particolare si
(0)
può costruire il seguente metodo iterativo: assegnato u2 su Γ, si risolvano i seguenti problemi
per k ≥ 1:

(k)
Lu1 = f
in Ω1 ,



(k)
(k−1)
(4.11)
u1 = u 2
su Γ,


 (k)
u1 = 0
su ∂Ω1 \ Γ,

(k)
Lu2 = f
in Ω2 ,




(k)
(k)
∂u1
∂u2
=
su Γ,


∂n
∂n

 (k)
u2 = 0
su ∂Ω2 \ Γ.
(4.12)
Si è utilizzata la (4.9) come condizione di Dirichlet su Γ per il sotto-problema associato a Ω 1 e
la (4.10) come condizione di Neumann su Γ per il problema assegnato su Ω 2 .
Si nota dunque che, a differenza del metodo di Schwarz, il metodo di Dirichlet-Neumann introduce un problema di Neumann sul secondo sotto-dominio Ω 2 . La soluzione del problema di
partenza è dunque ottenuta risolvendo, successivamente, un problema di Dirichlet e un problema
misto sui due sotto-domini. Inoltre il Teorema di equivalenza 4.1 garantisce che, se le successioni
(k)
(k)
{u1 } e {u2 } convergono, allora convergono sempre alla soluzione esatta del problema (4.2).
Il metodo di Dirichlet-Neumann è dunque consistente. Tuttavia la convergenza di tale metodo
non è sempre garantita. Andiamo a verificarlo con l’aiuto di un semplice esempio.
Esempio 4.2. Sia Ω = (a, b), γ ∈ (a, b), L = −d2 /dx2 e f = 0. Si hanno dunque i due seguenti
sotto-problemi:

(k)
−(u1 )00 = 0 a < x < γ,



(k)
(k−1)
(4.13)
u1 = u 2
x = γ,


 (k)
u1 = 0
x = a,

(k)
−(u2 )00 = 0
γ < x < b,



(k)
(k)
(4.14)
(u2 )0 = (u1 )0 x = γ,


 (k)
u2 = 0
x = b.
Procedendo come nell’Esempio 4.1, si può dimostrare che le successioni ottenute convergono
solamente se γ > (a + b)/2 come mostrato nelle Figure 4.4 e 4.5.
•
In generale, per i problemi di dimensione generica d > 1, si deve avere che la misura del sottodominio Ω1 sia più grande di quella del dominio Ω2 al fine di garantire la convergenza del metodo
(4.11)-(4.12). Tuttavia questo rappresenta un vincolo molto forte e difficile da soddisfare, soprattutto nel momento in cui si debbano utilizzare parecchi sotto-domini.
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
100
(0)
u2
PSfrag replacements
(2)
u2
a+b
2
a
γ
b
(1)
u2
Figura 4.4: Esempio di iterazione convergente per il metodo di Dirichlet-Neumann in 1D.
(0)
u2
PSfrag replacements
a
γ
a+b
2
b
(1)
u2
Figura 4.5: Esempio di iterazione non convergente per il metodo di Dirichlet-Neumann in 1D.
4.2. TRE CLASSICI METODI ITERATIVI BASATI SU DD
101
Tale limitazione viene superata introducendo una variante del metodo iterativo Dirichlet-Neumann,
rimpiazzando la condizione di Dirichlet (4.11)2 nel primo sotto-problema con la seguente
(k)
(k−1)
u1 = θu2
(k−1)
+ (1 − θ)u1
su Γ,
(4.15)
introducendo cioè un rilassamento che dipende da un parametro positivo θ. Si osservi che in
questo modo è sempre possibile ridurre l’errore da un’iterata alla successiva.
Nel caso rappresentato in Figura 4.5 si può facilmente dimostrare che, se si sceglie
(k−1)
θopt = −
u1
(k−1)
u2
(k−1)
− u1
,
il metodo converge in una sola iterata.
Più in generale si può dimostrare che, in dimensione 1 e 2, esiste un intervallo in cui scegliere il
parametro θ in modo da garantire la convergenza del metodo di Dirichlet-Neumann e, in particolare, è possibile determinare un valore massimo θmax < 1 per il parametro di rilassamento tale
che l’intervallo di cui sopra sia (0, θmax ).
4.2.3 Il metodo di Neumann-Neumann
Si consideri ancora una partizione del dominio Ω senza sovrapposizione e sia λ il valore della
soluzione u in corrispondenza dell’interfaccia Γ. Introduciamo lo schema iterativo seguente:
dato λ(0) su Γ, si risolvano i seguenti problemi per k ≥ 0:

(k+1)
= f in Ωi ,

 −4ui

(k+1)
(4.16)
ui
= λ(k)
su Γ,


 (k+1)
ui
=0
su ∂Ωi \ Γ,

(k+1)
=0
in Ωi ,
−4ψi




(k+1)
(k+1)
(k+1)
∂u1
∂u2
∂ψi
(4.17)
=
−
su Γ,


∂n
∂n
∂n

 (k+1)
ψi
=0
su ∂Ωi \ Γ,
con
(k+1)
(k+1)
,
(4.18)
λ(k+1) = λ(k) − θ σ1 ψ1|Γ − σ2 ψ2|Γ
essendo θ un parametro d’accelerazione positivo, σ1 e σ2 due coefficienti positivi.
Esempio 4.3. Fatte per i dati le stesse scelte dell’Esempio 4.2, si risolvono ora i seguenti quattro
sotto-problemi:

(k+1) 00
−(u1
) = 0 a < x < γ,



(k+1)
(4.19)
u1
= λ(k)
x = γ,


 (k+1)
u1
=0
x = a,
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
102

(k+1) 00
−(u2
) = 0 γ < x < b,



(k+1)
u2
= λ(k)
x = γ,


 (k+1)
u2
=0
x = b,
con

(k+1) 00
−(ψ1
) =0



(k+1) 0
(k+1) 0
(k+1) 0
(ψ1
) = (u1
) − (u2
)


 (k+1)
ψ1
=0

(k+1) 00
−(ψ2
) =0



(k+1) 0
(k+1) 0
(k+1) 0
(ψ2
) = (u1
) − (u2
)


 (k+1)
ψ2
=0
a < x < γ,
x = γ,
(4.21)
x = a,
γ < x < b,
x = γ,
(4.22)
x = b,
(k+1)
(k+1)
λ(k+1) = λ(k) − θ σ1 ψ1
− σ 2 ψ2
θ, σ1 e σ2 definiti come sopra.
(4.20)
x = γ,
(4.23)
•
4.3 Formulazione multi-dominio del problema di Poisson ed
equazioni di interfaccia
In questa sezione scegliamo L = −4 e consideriamo il problema di Poisson con condizioni al
bordo di Dirichlet omogenee:
(
−4u = f in Ω,
(4.24)
u=0
su ∂Ω.
Il Teorema di equivalenza 4.1 ci permette di riscrivere questo problema nella formulazione multidominio seguente:

−4u1 = f in Ω1 ,






u1 = 0
su ∂Ω1 \ Γ,





 −4u2 = f in Ω2 ,
(4.25)
u2 = 0
su ∂Ω2 \ Γ,





u1 = u 2
su Γ,






 ∂u1 = ∂u2 su Γ.
∂n
∂n
4.3.1 L’operatore di Steklov-Poincaré
Indichiamo ora con λ il valore incognito della soluzione u del problema (4.24) sull’interfaccia Γ,
cioè λ = u|Γ . Se si conoscesse a priori il valore di λ su Γ, si potrebbero risolvere i due problemi
4.3. FORMULAZIONE MULTI-DOMINIO
103
seguenti con condizioni al bordo di Dirichlet, in corrispondenza di Γ, pari a λ:

−4wi = f in Ωi
(i = 1, 2),


wi = 0
su ∂Ωi \ Γ,


wi = λ
su Γ.
(4.26)
Si noti che tali problemi possono essere risolti indipendentemente uno dall’altro.
Il problema è dunque come ottenere il valore di λ su Γ. Possiamo scomporre intanto le funzioni
wi , soluzioni di (4.26), nel seguente modo:
wi = wi∗ + u0i ,
dove wi∗ e u0i sono le soluzioni dei due problemi

−4wi∗ = f in Ωi
(i = 1, 2),


wi∗ = 0
su ∂Ωi ∩ ∂Ω,

 ∗
wi = 0
su Γ,
e
(4.27)

−4u0i = 0 in Ωi
(i = 1, 2),


su ∂Ωi ∩ ∂Ω,
u0i = 0

 0
ui = λ
su Γ,
(4.28)
rispettivamente. Osserviamo che le funzioni wi∗ (i = 1, 2) dipendono solamente dal dato f ,
mentre u0i (i = 1, 2) dipende soltanto dal valore λ su Γ. Inoltre, entrambe queste dipendenze
sono lineari, e quindi possiamo scrivere wi∗ come wi∗ = Gi f , con Gi operatore lineare, e u0i come
u0i = Hi λ, essendo Hi l’operatore d’estensione armonico di λ sul dominio Ωi , per i = 1, 2.
Ora, confrontando formalmente il problema (4.25) con (4.26), si osserva che l’uguaglianza
ui = wi∗ + u0i
(i = 1, 2)
vale se e solamente se le funzioni wi soddisfano la condizione (4.25)6 sulle derivate normali su
Γ, ovvero se e solamente se
∂w1
∂w2
=
su Γ.
∂n
∂n
Utilizzando le notazioni introdotte sopra, possiamo riscrivere quest’ultima condizione come
∂
∂
(w1∗ + u01 ) =
(w ∗ + u02 ),
∂n
∂n 2
ovvero
e quindi
∂
∂
(G1 f + H1 λ) =
(G2 f + H2 λ),
∂n
∂n
∂H1 ∂H2
−
∂n
∂n
λ=
∂G2 ∂G1
−
∂n
∂n
f
su Γ.
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
104
λ
H1 λ
H2 λ
l1
a
Ω1
γ
l2
Ω2
b
Figura 4.6: Estensioni armoniche in una dimensione.
Abbiamo ottenuto in tal modo un’equazione sull’interfaccia Γ per l’incognita λ, nota come
equazione d’interfaccia di Steklov-Poincaré, espressa, in forma più compatta, come
Sλ = χ
su Γ.
(4.29)
S è l’operatore pseudo-differenziale di Steklov-Poincaré definito formalmente come
2
X ∂
∂
∂
Hi µ,
H1 µ −
H2 µ =
Sµ :=
∂n
∂n
∂n
i
i=1
(4.30)
mentre χ è un funzionale lineare dipendente dai dati del problema:
2
χ :=
X ∂
∂
∂
G2 f −
G1 f = −
Gi f.
∂n
∂n
∂n
i
i=1
(4.31)
Ricordiamo che ni è la normale esterna al sotto-dominio Ωi , per i = 1, 2.
∂
Si noti che S = S1 + S2 , dove Si =
Hi , per i = 1, 2.
∂ni
Esempio 4.4. Consideriamo un caso semplice monodimensionale per fornire un esempio di operatore S. Sia Ω = (a, b) ⊂ R come illustrato in Figura 4.6 e Lu = −u00 . Se suddividiamo Ω in
due sotto-domini senza sovrapposione, l’interfaccia Γ si riduce ad un punto solo γ ∈ (a, b), e
l’operatore di Steklov-Poincaré S diventa
dH1 dH2
1
1
Sλ =
λ=
+
λ,
−
dx
dx
l1 l2
con l1 = γ − a e l2 = b − γ.
•
Usando la formula di Green si può fornire la seguente caratterizzazione variazionale dell’operatore di Steklov-Poincaré:
< Si λ, µ >= ai (Hi λ, Hi µ),
1/2
(4.32)
∀λ, µ ∈ Λ.
(4.33)
∀λ, µ ∈ Λ = H00 (Γ)
per i = 1, 2. Pertanto
< Sλ, µ >=
2
X
i=1
ai (Hi λ, Hi µ),
4.3. FORMULAZIONE MULTI-DOMINIO
105
4.3.2 Proprietà dell’operatore di Steklov-Poincaré
Sia Ωi ⊂ Rd un dominio limitato e Lipschitziano, con i = 1, 2. Definiamo lo spazio
Vi = {v ∈ H 1 (Ωi )| v = 0 su ∂Ωi \ Γ}.
Teorema 4.2 (Disuguaglianza di Poincaré).
∀v ∈ Vi
kvk2L2 (Ωi ) ≤ CΩi k∇vk2L2 (Ωi ) .
(4.34)
Corollario 4.1.
(
p
1 + CΩi )−1 kvkH 1 (Ωi ) ≤ k∇vkL2 (Ωi ) ≤ kvkH 1 (Ωi ) ,
∀v ∈ Vi .
(4.35)
Per ogni funzione v ∈ H 1 (Ωi ) è possibile definire l’operatore di traccia
γ : Vi → L2 (Γ) t.c.
γw = w|Γ
∀w ∈ C 0 (Ω) ∩ Vi .
È noto che, se d > 1, H 1 (Ω) * C 0 (Ω) ma tale definizione può essere estesa per densità, ovvero
esiste un operatore γ continuo tale che
kv|Γ kL2 (Γ) = kγvkL2 (Γ) ≤ CkvkH 1 (Ωi ) .
Inoltre, poiché γ(Vi ) ( L2 (Γ), si definisce l’immagine di γ in L2 (Γ) introducendo lo spazio di
traccia Λ = Im(γ).
Si può dimostrare che Λ ⊂ L2 (Γ) e che è possibile definire una norma sullo spazio Λ
kηkΛ := inf kvkH 1 (Ωi )
v∈Vi
v|Γ =η
rispetto alla quale (Λ, k · kΛ ) sia uno spazio di Banach. L’operatore γ : Vi → Λ è continuo
rispetto a tale norma, cioè si ha che
∃C > 0 :
kv|Γ kΛ ≤ CkvkH 1 (Ωi ) ,
∀v ∈ H 1 (Ωi ).
Consideriamo ora alcune proprietà del rilevamento armonico Hi η.
Si consideri il problema armonico seguente

−4(Hi η) = 0 in Ωi ,


Hi η = η
su Γ,


Hi η = 0
su ∂Ωi \ Γ.
(4.36)
La forma debole associata è data da: trovare Hi η ∈ Vi , con Hi η = η su Γ ed η ∈ Λ, tale che
Z
∇(Hi η) · ∇v dΩi = 0
∀v ∈ Vi0 ,
(4.37)
Ωi
dove Vi0 = {v ∈ H 1 (Ωi )|v = 0 su ∂Ωi }.
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
106
Teorema 4.3 (Teorema di estensione). Esistono due costanti C1 , C2 > 0 tali che
C1 kηkΛ ≤ kHi ηkH 1 (Ωi ) ≤ C2 kηkΛ
i = 1, 2,
∀η ∈ Λ.
(4.38)
Dimostrazione. Grazie alla (4.36)2 , si ha che γ(Hi η) = Hi η |Γ = η e dunque la disuguaglianza a
sinistra segue dalla continuità di γ:
kηkΛ ≤ CkHi ηkH 1 (Ωi ) .
D’altro canto, poiché η ∈ Λ = Im(γ), esiste una funzione u∗ ∈ H 1 (Ωi ) tale che γu∗ = η.
Possiamo allora riscrivere Hi η come la somma Hi η = u∗ + u0 , essendo u0 una funzione di
H01 (Ωi ) tale che γu0 = 0. La forma debole del problema (4.36) può essere riformulata come:
trovare u0 ∈ H01 (Ωi ) tale che
Z
Z
0
∇u · ∇v dΩi = −
∇u∗ · ∇v dΩi
∀v ∈ H01 (Ωi ).
Ωi
Ωi
Sia v = u0 ; dunque
k∇u0 k2L2 (Ωi )
=−
Z
∇u∗ · ∇u0 dΩi .
Ωi
Grazie alla disuguaglianza di Poincaré (4.34), si ha
1
ku0 k2H 1 (Ωi ) ≤ k∇u∗ kL2 (Ωi ) k∇u0 kL2 (Ωi ) ≤ ku∗ kH 1 (Ωi ) ku0 kH 1 (Ωi ) ,
1 + C Ωi
e dunque
ku0 kH 1 (Ωi ) ≤ (1 + CΩi )ku∗ kH 1 (Ωi ) .
Ma
kHi ηkH 1 (Ωi ) = ku0 + u∗ kH 1 (Ωi ) ≤ ku0 kH 1 (Ωi ) + ku∗ kH 1 (Ωi ) ≤ (2 + CΩi )ku∗ kH 1 (Ωi )
e tale disuguaglianza è vera per ogni u∗ ∈ Vi tale che u∗ |Γ = η. Da ciò segue
kHi ηkH 1 (Ωi ) ≤ (2 + CΩi ) inf
ku∗ kH 1 (Ωi ) = (2 + CΩi )kηkΛ .
∗
u ∈Vi
u∗ |Γ =η
Come abbiamo visto, l’operatore di Steklov-Poincaré S può essere riscritto come la somma di
due operatori S1 e S2 tali che
Si : η ∈ Λ →
∂
(Hi η) ∈ Λ0
∂ni
Definiamo la forma bilineare seguente su Λ:
Z
∂
< Si η, µ >:= µ
(Hi η) dγ
∂ni
Γ
i = 1, 2.
∀η, µ ∈ Λ.
4.3. FORMULAZIONE MULTI-DOMINIO
107
Scriviamo la formulazione debole del problema (4.36) scegliendo una funzione test v ∈ V i :
trovare Hi η ∈ Vi tale che
Z
∇(Hi η) · ∇v dΩi =
Ωi
Z
Γ
∂
(Hi η) v dγ.
∂ni
(4.39)
Per ogni µ ∈ Λ, esiste Ri µ ∈ Vi tale che Ri µ|Γ = µ. Scegliendo allora v = Ri µ e grazie alla
(4.39), si ha
Z
< Si η, µ >=
∇(Hi η) · ∇(Ri µ) dΩi ,
(4.40)
Ωi
e dunque
< Sη, µ >=
2
X
< Si η, µ >=
i=1
2 Z
X
i=1
∇(Hi η) · ∇(Ri µ) dΩi .
(4.41)
Ωi
Infine, essendo Ri un operatore d’estensione qualsiasi, possiamo fare la scelta particolare R i =
Hi , da cui:
Z
< Si η, µ >=
∇(Hi η) · ∇(Hi µ) dΩi .
(4.42)
Ωi
Proprietà 4.1. La forma bilineare (4.42) è simmetrica, coerciva e continua.
Dimostrazione. La simmetria è evidente. Per quanto riguarda la coercività, grazie alla (4.35) ed
alla (4.38), si ha
< Si η, η > = k∇(Hi η)k2L2 (Ωi ) ≥
1
2
kHi ηk2H 1 (Ωi ) ≥
kηk2Λ ,
1 + C Ωi
1 + C Ωi
ovvero
< Si η, η > ≥ αkηk2Λ
∀η ∈ Λ,
dove α è la costante positiva
α :=
2
.
1 + C Ωi
Infine, grazie ancora a (4.38),
< Si η, µ >=
Z
∇(Hi η) · ∇(Hi µ) dΩi ≤ k∇(Hi η)kL2 (Ωi ) k∇(Hi µ)kL2 (Ωi )
Ωi
≤ C22 kηkΛ kµkΛ
e dunque la forma è continua con costante di continuità β = C22 > 0.
∀η, µ ∈ Λ,
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
108
4.3.3 Equivalenza tra il metodo di Dirichlet-Neumann e il metodo di Richardson
Il metodo di Dirichlet-Neumann può essere reinterpretato come un metodo iterativo di Richardson precondizionato per risolvere l’equazione d’interfaccia di Steklov-Poincaré.
Consideriamo ancora un dominio Ω suddiviso in due sotto-domini senza sovrapposizione Ω 1 e
Ω2 con interfaccia Γ.
Nella Sezione 4.2.2 abbiamo considerato il metodo di Dirichlet-Neumann, e dopo aver introdotto
l’equazione d’interfaccia di Steklov-Poincaré abbiamo anticipato che tale metodo è equivalente
allo schema di Richardson applicato alla risoluzione dell’equazione di Steklov-Poincaré su Γ.
Vogliamo dimostrare ora tale risultato di equivalenza sul problema di Laplace.
Ricordiamo innanzitutto il metodo di Dirichlet-Neumann: dato λ0 , si risolva, per k ≥ 1,

(k)
−4u1 = f1 in Ω1 ,



(k)
(4.43)
u1 = λ(k−1) su Γ,


 (k)
u1 = 0
su ∂Ω1 \ Γ,

(k)

−4u2 = f2 in Ω2 ,




(k)
(k)
∂u2
∂u1
=
su Γ,

∂n2
∂n2



 (k)
u2 = 0
su ∂Ω2 \ Γ,
(k)
λ(k) = θu2
|Γ
+ (1 − θ)λ(k−1) .
(4.44)
(4.45)
Ricordiamo anche il metodo di Richardson precondizionato per l’equazione d’interfaccia di
Steklov-Poincaré (4.29):
P (λ(k) − λ(k−1) ) = θ(χ − Sλ(k−1) ),
(4.46)
essendo P un precondizionatore opportuno.
Si ha il risultato seguente:
Teorema 4.4. Il metodo di Dirichlet-Neumann (4.43)-(4.45) è equivalente al metodo di Richardson precondizionato (4.46) con precondizionatore P = S 2 = ∂(H2 µ)/∂n2 . In versione
algebrica tale metodo è equivalente allo schema di Richardson precondizionato
Ph (λ(k) − λ(k−1) ) = θ(χΓ − Σλ(k−1) ),
con precondizionatore Ph = Σ2 .
Dimostrazione. Consideriamo innanzitutto il caso differenziale. Abbiamo già visto che la solu(k)
zione u1 di (4.43) può essere scritta come
(k)
u1 = H1 λ(k−1) + G1 f1 .
(4.47)
4.3. FORMULAZIONE MULTI-DOMINIO
109
Inoltre G2 f2 soddisfa il problema differenziale
(
−4(G2 f2 ) = f2 in Ω2 ,
G2 f2 = 0
su ∂Ω2 .
(k)
(k)
Dunque, dalla (4.44) si ha che u2 − G2 f2 soddisfa l’equazione −4(u2 − G2 f2 ) = 0 sul
(k)
dominio Ω2 , assieme alla condizione u2 − G2 f2 = 0 su ∂Ω2 \ Γ, e
(k)
∂
∂u
∂
(k)
(u2 − G2 f2 ) = 2 −
(G2 f2 )
∂n2
∂n2
∂n2
(k)
(k)
=
∂
∂u
∂
∂u1
−
(G2 f2 ) = − 1 +
(G2 f2 )
∂n2
∂n2
∂n
∂n
su Γ,
(k)
con n = n1 = −n2 , essendo ni la normale esterna al dominio Ωi . Si ha allora che u2 − G2 f2 è
la soluzione del problema differenziale

(k)

−4(u2 − G2 f2 ) = 0
in Ω2 ,





(k)
∂
∂u1
∂
(k)
(4.48)
(u
−
G
f
)
=
−
+
(G2 f2 ) su Γ,
2 2

2

∂n
∂n
∂n

2


 (k)
u2 − G 2 f 2 = 0
su ∂Ω2 \ Γ.
(k)
(k)
In particolare u2
|Γ
= (u2 − G2 f2 )|Γ . Ora ricordiamo che
S i : µ |Γ →
∂
(Hi µ)|Γ
∂ni
(i = 1, 2),
cioè che l’operatore di Steklov-Poincaré ci permette di passare da un dato di Dirichlet ad un dato
di Neumann su Γ. Il suo inverso Si−1 , partendo da un dato di Neumann, ci fornisce invece un
dato di Dirichlet su Γ.
Si può esprimere ciò dicendo anche che S2−1 η = µ se esiste w tale che w|Γ = µ e w è la soluzione
del seguente problema:

 −4w = 0 in Ω2 ,


∂w
(4.49)
= η su Γ,

∂n

2

w=0
su ∂Ω2 \ Γ.
Allora, se si pone
(k)
∂u1
∂
+
(G2 f2 ),
∂n
∂n
e se si confronta la (4.48) con la (4.49), si può concludere che
η=−
(k)
u2
(k)
|Γ
= (u2 − G2 f2 )|Γ = S2−1
!
(k)
∂
∂u1
+
(G2 f2 ) .
−
∂n
∂n
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
110
Ma, grazie alla (4.47), si ottiene
∂
∂
∂
(k)
−1
(k−1)
u2 |Γ = S 2
− (H1 λ
)−
(G1 f1 ) +
(G2 f2 ) = S2−1 (−S1 λ(k−1) + χ),
∂n
∂n
∂n
ricordata la definizione (4.31) di χ. Utilizzando la (4.45), possiamo dunque scrivere che
λ(k) = θ S2−1 (−S1 λ(k−1) + χ) + (1 − θ)λ(k−1)
λ(k) − λ(k−1) = θ S2−1 (−S1 λ(k−1) + χ) − λ(k−1) .
Ma −S1 = S2 − S, dunque
λ(k) − λ(k−1) = θ S2−1 ((S2 − S)λ(k−1) + χ) − λ(k−1)
= θS2−1 (χ − Sλ(k−1) ),
ovvero
S2 (λ(k) − λ(k−1) ) = θ(χ − Sλ(k−1) ).
4.4 Approssimazione con elementi finiti del problema di Poisson e formulazione per sotto-domini
In questa sezione introduciamo l’approssimazione discreta del problema (4.24) basata sullo
schema degli elementi finiti, e studiamo la forma algebrica della corrispondente formulazione
multi-dominio (4.25) e dell’equazione d’interfaccia di Steklov-Poincaré (4.29).
Innanzitutto identifichiamo V con lo spazio H01 (Ω) definito in (1.50). Dal Capitolo 1, ne segue che, indicata con Th una triangolazione del dominio Ω, la formulazione alla Galerkin del
problema (4.24) può essere scritta come
trovare uh ∈ Vh : a(uh , vh ) = F (vh )
con
a(w1 , w2 ) =
Z
∀vh ∈ Vh ,
∇w1 · ∇w2 dΩ ∀w1 , w2 ∈ V,
(4.50)
(4.51)
Ω
F (v) =
Z
f v dΩ ∀v ∈ V,
(4.52)
Ω
k
h
= vh ∈ Xhk : vh |∂Ω = 0 spazio degli elementi finiti di grado k con base {ϕj }N
Vh = Xh,∂Ω
j=1 .
(1)
Introduciamo la partizione seguente dei nodi del dominio Ω: siano {x j , 1 ≤ j ≤ N1 } i nodi
(2)
(Γ)
ubicati dentro il sotto-dominio Ω1 , {xj , 1 ≤ j ≤ N2 } quelli in Ω2 e, infine, {xj , 1 ≤ j ≤
NΓ } quelli posizionati sull’interfaccia Γ. Ripartiamo in modo analogo anche le funzioni di base:
4.4. APPROSSIMAZIONE CON ELEMENTI FINITI
111
(1)
(1)
(2)
indicheremo dunque con ϕj le funzioni di base associate ai nodi xj , con ϕj quelle associate
(2)
(Γ)
(Γ)
ai nodi xj , e con ϕj quelle relative ai nodi xj di interfaccia. Ciò significa che
δij 1 ≤ i, j ≤ Nα se α = β,
(α) (β)
ϕj (xj ) =
0
se α 6= β,
con α, β = 1, 2, Γ, e δij il simbolo di Kronecker.
Scelto ora in (4.50) vh coincidente con una funzione test, possiamo dare la seguente formulazione
equivalente della (4.50): trovare uh ∈ Vh tale che

(1)
(1)
a(uh , ϕi ) = F (ϕi ) ∀i = 1, . . . , N1



(2)
(2)
(4.53)
a(uh , ϕj ) = F (ϕj ) ∀j = 1, . . . , N2



(Γ)
(Γ)
a(uh , ϕk ) = F (ϕk ) ∀k = 1, . . . , NΓ .
Sia ora
ai (v, w) :=
Z
∇v · ∇w dΩ
∀v, w ∈ V, i = 1, 2
Ωi
la restrizione della forma bilineare a(., .) al sotto-dominio Ωi e Vi,h = {v ∈ H 1 (Ωi ) | v =
0 su ∂Ωi \ Γ} (i = 1, 2). Analogamente sia Fi la restrizione del funzionale lineare F al sotto(i)
dominio Ωi , per i = 1, 2. Indichiamo infine con uh = uh|Ωi la restrizione di uh al sotto-dominio
Ωi , con i = 1, 2.
(2)
(1)
Il problema (4.53) può così essere riscritto nella forma equivalente: trovare u h ∈ V1,h , uh ∈
V2,h tali che

(1)
(1)
(1)

∀i = 1, . . . , N1 ,
a1 (uh , ϕi ) = F1 (ϕi )





(2)
(2)
(2)
∀j = 1, . . . , N2
a2 (uh , ϕj ) = F2 (ϕj )
(4.54)

(1)
(Γ)
(2)
(Γ)

a1 (uh , ϕk |Ω1 ) + a2 (uh , ϕk |Ω2 )



 = F (ϕ(Γ) | ) + F (ϕ(Γ) | )
∀k = 1, . . . , NΓ .
1
2
k Ω1
k Ω2
Si osservi che la condizione (4.25)5 di continuità della soluzione all’interfaccia è automatica(i)
mente verificata grazie alla definizione delle funzioni uh . Osserviamo inoltre che le equazioni
(4.54)1 -(4.54)3 corrispondono alla discretizzazione ad elementi finiti delle (4.25) 1 -(4.25)6 , rispettivamente.
Possiamo ora esprimere la funzione uh come una somma rispetto alla base dello spazio Vh :
uh (x) =
N1
X
(1)
(1)
uh (xj )ϕj (x)
j=1
+
NΓ
X
+
N2
X
(2)
(2)
uh (xj )ϕj (x)
j=1
(Γ)
(Γ)
uh (xj )ϕj (x),
(4.55)
j=1
(α)
dove uh (xj ), con α = 1, 2, Γ, sono i coefficienti della combinazione lineare, indicati d’ora in
(α)
(α)
avanti come uj = uh (xj ), per j = 1, . . . , Nα e α = 1, 2, Γ.
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
112
Utilizzando la relazione (4.55), possiamo riscrivere il problema (4.54) nel seguente modo:

NΓ
N1

X
X

(Γ)
(Γ)
(1)
(1)
(1)
(1)
(1)


uj a1 (ϕj , ϕi ) = F1 (ϕi ) ∀i = 1, . . . , N1 ,
uj a1 (ϕj , ϕi ) +


 j=1

j=1


NΓ
N2

X
X


(Γ)
(Γ)
(2)
(2)
(2)
(2)
(2)


uj a2 (ϕj , ϕi ) = F2 (ϕi ) ∀i = 1, . . . , N2 ,
uj a2 (ϕj , ϕi ) +



j=1
 j=1
N
Γ
h
i
X (Γ)
(4.56)
(Γ)
(Γ)
(Γ)
(Γ)

u
a
(ϕ
,
ϕ
)
+
a
(ϕ
,
ϕ
)

1
2
j
j
i
j
i



j=1



N1
N2

X
X

(2)
(2)
(Γ)
(1)
(1)
(Γ)

 +
uj a1 (ϕj , ϕi ) +
uj a2 (ϕj , ϕi )




j=1
j=1


 = F (ϕ(Γ) | ) + F (ϕ(Γ) | )
∀i = 1, . . . , NΓ .
1
Ω1
2
Ω2
i
i
Introduciamo ora le seguenti matrici e vettori:
(1)
(1)
(A1Γ )ij = a1 (ϕj , ϕi ),
(2)
(2)
(A2Γ )ij = a2 (ϕj , ϕi ),
(A11 )ij = a1 (ϕj , ϕi ),
(A22 )ij = a2 (ϕj , ϕi ),
(Γ)
(Γ)
(1)
(Γ)
(Γ)
(1)
(Γ)
(2)
(Γ)
(Γ)
(2)
(Γ)
(A1ΓΓ )ij = a1 (ϕj , ϕi ), (A2ΓΓ )ij = a2 (ϕj , ϕi ),
(AΓ1 )ij = a1 (ϕj , ϕi ),
(1)
e sia
(f1 )i = F1 (ϕi ),
(Γ)
f1Γ i = F1 (ϕi ),
(AΓ2 )ij = a2 (ϕj , ϕi ),
(2)
(f2 )i = F2 (ϕi ),
(Γ)
(1)
f2Γ i = F2 (ϕi , ϕi ),
u = (u1 , u2 , λ)T ,
(1)
(2)
(Γ)
con u1 = uj , u2 = uj e λ = uj .
Il problema (4.56) può dunque essere scritto in forma algebrica nel seguente modo:

A11 u1 + A1Γ λ = f1 ,



A22 u2 + A2Γ λ = f2 ,


 A u + A u + A(1) + A(2) λ = f Γ + f Γ ,
Γ1 1
Γ2 2
1
2
ΓΓ
ΓΓ
o, in forma matriciale, come Au = f , ovvero


 

u1
f1
A11 0 A1Γ
 0 A22 A2Γ   u2  =  f2  ,
λ
fΓ
AΓ1 AΓ2 AΓΓ
(1)
(2)
avendo indicato con AΓΓ = AΓΓ + AΓΓ e con fΓ = f1Γ + f2Γ .
(4.57)
(4.58)
4.4. APPROSSIMAZIONE CON ELEMENTI FINITI
113
4.4.1 Il complemento di Schur
Consideriamo ora l’equazione all’interfaccia di Steklov-Poincaré (4.29). Vogliamo approssimarla con gli elementi finiti. Poiché λ rappresenta il valore di u su Γ, essa corrisponde, in dimensione
finita, al vettore λ dei valori di uh sull’interfaccia.
Applicando il metodo di eliminazione Gaussiana al sistema (4.58), possiamo scrivere un problema algebrico per l’incognita vettoriale λ.
Le matrici A11 e A22 sono invertibili in quanto associate a due problemi di Dirichlet omogenei
per l’equazione di Laplace, e quindi si può facilmente ricavare
u1 = A−1
11 (f1 − A1Γ λ) ,
(4.59)
u2 = A−1
22 (f2 − A2Γ λ) .
(4.60)
e
Si può allora rimpiazzare (4.59) e (4.60) in (4.57), ottenendo così:
(1)
(2)
−1
(f
−
A
λ)
+
A
+
A
−
A
λ)
+
A
A
AΓ1 A−1
(f
2
2Γ
1
1Γ
Γ2
22
11
ΓΓ
ΓΓ λ = fΓ ,
ovvero
h
i
(1)
(2)
−1
AΓΓ − AΓ1 A−1
A
+
A
−
A
A
A
λ=
Γ2 22 2Γ
11 1Γ
ΓΓ
fΓ −
AΓ1 A−1
11 f1
−
(4.61)
AΓ2 A−1
22 f2 .
Introduciamo ora le seguenti definizioni:
(i)
Σi := AΓΓ − AΓi A−1
ii AiΓ
i = 1, 2,
(4.62)
Σ := Σ1 + Σ2 ,
(4.63)
−1
χΓ := fΓ − AΓ1 A−1
11 f1 − AΓ2 A22 f2 .
(4.64)
e
Dalla (4.61) possiamo allora derivare l’approssimazione ad elementi finiti desiderata per l’equazione di Steklov-Poincaré (4.29):
Σλ = χΓ ,
(4.65)
rappresentando Σ e χΓ le approssimazioni alla Galerkin di S e χ, rispettivamente.
La matrice Σ rappresenta il cosiddetto complemento di Schur della matrice A rispetto alle incognite vettoriali u1 e u2 , mentre le matrici Σi sono i complementi di Schur relativi ai sotto-domini
Ωi (i = 1, 2).
Il problema (4.65) è un sistema lineare nella sola incognita λ. Risolto tale sistema rispetto
all’incognita λ, in virtù delle (4.59) e (4.60), possiamo calcolare u 1 e u2 . Tale calcolo equivale a risolvere due problemi di Poisson in corrispondenza dei due sotto-domini Ω 1 e Ω2 , con
(i)
condizione al bordo di Dirichlet uh |Γ = λh , (i = 1, 2), sull’interfaccia Γ.
Per quanto concerne le proprietà del complemento di Schur Σ rispetto alla matrice A, si può
dimostrare il risultato seguente:
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
114
Lemma 4.1. La matrice Σ soddisfa le proprietà seguenti:
1. se A è una matrice singolare, allora Σ è singolare;
2. se A (rispettivamente, Aii ) è simmetrica, allora Σ (rispettivamente, Σi ) è simmetrica;
3. se A è definita positiva, allora Σ è definita positiva.
Per quanto riguarda il numero di condizionamento delle matrici A e Σ, si ricorda che A è una
matrice mal condizionata, il suo numero di condizionamento essendo dato da
K(A) '
C
.
h2
Quanto alla matrice Σ, si può invece dimostrare che
K(Σ) '
C
.
h
(4.66)
Osserviamo infine che, nel caso specifico preso in considerazione, la matrice A (e dunque la
matrice Σ, grazie al Lemma 4.1) è una matrice simmetrica definita positiva. È possibile dunque
utilizzare il metodo del Gradiente Coniugato per risolvere il sistema (4.65). Ad ogni passo,
bisognerà risolvere due problemi di Dirichlet indipendenti sui sotto-domini Ω i .
4.4.2 L’operatore di Steklov-Poincaré discreto
Ci proponiamo di trovare l’operatore discreto associato al complemento di Schur. A tal fine de0
finiamo, oltre allo spazio Vi,h introdotto all’inizio di questa sezione, lo spazio Vi,h
generato dalle
(i)
funzioni {ϕj } associate ai soli nodi interni al dominio Ωi , e lo spazio Λh generato dall’insieme
(Γ)
{ϕj | }. La formulazione discreta del problema (4.37) può così essere espressa come: trovare
Γ
Hi,h ηh ∈ Vi,h , con Hi,h ηh = ηh su Γ, tale che
Z
0
∇(Hi,h ηh ) · ∇vh dΩi = 0
∀vh ∈ Vi,h
.
(4.67)
Ωi
Riscriviamo la soluzione Hi,h ηh come combinazione lineare delle funzioni di base
Hi,h ηh =
Ni
X
j=1
(i) (i)
uj ϕ j
+
NΓ
X
(Γ)
η k ϕk
k=1
|Ω i ,
arrivando così a riscrivere la (4.67) sotto forma matriciale come
Aii u(i) = −AiΓ η.
(4.68)
Il teorema d’estensione uniforme enunciato nel caso continuo vale analogamente nel caso discreto:
4.4. APPROSSIMAZIONE CON ELEMENTI FINITI
115
Teorema 4.5. Esistono due costanti Ĉ1 , Ĉ2 > 0 che dipendono dalle dimensioni relative di Ω1 e
Ω2 ma independenti da h, tali che
Ĉ1 kηh kΛ ≤ kHi,h ηh kH 1 (Ωi ) ≤ Ĉ2 kηh kΛ
i = 1, 2,
∀ηh ∈ Λh .
(4.69)
Dimostrazione. Poiché ηh è la traccia sull’interfaccia Γ sia di H1,h ηh che di H2,h ηh , esistono,
grazie alla disuguaglianza di traccia, due costanti C1∗ e C2∗ positive, tali che
kηh kΛ ≤ Ci∗ kHi,h ηh kH 1 (Ωi )
∀ηh ∈ Λh ,
i = 1, 2.
La prima disuguaglianza della (4.69) è dunque immediatamente verificata fatta la scelta
Ĉ1 = min
1 1 ,
.
C1∗ C2∗
D’altro canto, abbiamo che
kHi,h ηh kH 1 (Ωi ) ≤ kHi,h ηh − Hi ηh kH 1 (Ωi ) + kHi ηh kH 1 (Ωi ) .
(4.70)
Ora, grazie al Teorema 4.3, si ha che
kHi ηh kH 1 (Ωi ) ≤ C2 kηh kΛ
i = 1, 2.
Poiché ηh è una funzione polinomiale a tratti e continua sull’interfaccia Γ, la soluzione H i ηh
∗
appartiene ad un opportuno spazio di Sobolev H 1+s (Ωi ), con s∗ > 1/2. Vale dunque la seguente
stima di regolarità:
kHi ηh kH 1+s∗ (Ωi ) ≤ Ckηh kH 1/2+s∗ (Γ) .
Sul primo termine della (4.70) utilizziamo invece la stima per l’errore di discretizzazione degli
elementi finiti:
∗
kHi,h ηh − Hi ηh kH 1 (Ωi ) ≤ Chs kHi ηh kH 1+s∗ (Ωi ) .
Infine, grazie alla disuguaglianza inversa standard, si ha che
∗
hs kηh kH 1/2+s∗ (Γ) ≤ Ckηh kΛ
∀ηh ∈ Λh ,
con C costante indipendente da h. Possiamo quindi concludere che, per un’opportuna costante
Ĉ2 indipendente da h,
kHi,h ηh kH 1 (Ωi ) ≤ Ĉ2 kηh kΛ
ovvero dimostrare la seconda disuguaglianza della (4.69).
Partendo da tale teorema, è facile verificare il seguente risultato
Corollario 4.2. Gli operatori discreti di Steklov-Poincaré S1,h , S2,h e Sh = S1,h + S2,h sono
continui e coercivi in Λh , uniformemente rispetto ad h. Inoltre sono equivalenti spettralmente in
modo uniforme (ovvero le costanti dell’equivalenza spettrale non dipendono da h).
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
116
Corollario 4.3. Esistono due costanti K1 , K2 > 0, independenti da h, tali che
K1 kH1,h ηh kH 1 (Ω1 ) ≤ kH2,h ηh kH 1 (Ω2 ) ≤ K2 kH1,h ηh kH 1 (Ω1 )
∀ηh ∈ Λh .
(4.71)
Possiamo a questo punto fornire la seguente
Definizione 4.1. Definiamo l’operatore Si,h : Λh → Λ0h dato da
Z
< Si,h ηh , µh >=
∇(Hi,h ηh ) · ∇(Hi,h µh ), ∀ηh , µh ∈ Λh .
(4.72)
Ωi
Lemma 4.2. L’operatore di Steklov-Poincaré discreto può essere espresso in funzione del complemento di Schur come
< Si,h ηh , µh >= µT Σi η
dove
ηh =
NΓ
X
(Γ)
η k ϕ k |Γ ,
µh =
(Γ)
µk ϕ k
|Γ
k=1
k=1
e
NΓ
X
η = (η1 , . . . , ηNΓ )T ,
µ = (µ1 , . . . , µNΓ )T .
Pertanto, Sh = S1,h + S2,h verifica la relazione
< Sh ηh , µh >= µT Σ η.
(4.73)
Dimostrazione. Si ha:
< Si,h ηh , µh > = ai (Hi,h ηh , Hi,h µh )
= ai
NΓ
X
(i)
uj ϕ j
+
j=1
=
NΓ
X
NΓ
X
k=1
(Γ)
η k ϕ k |Ω ,
i
(i)
(i)
wl ai (ϕj , ϕl )uj
+
NΓ
X
(i)
w l ϕl
+
NΓ
X
(Γ)
wl ai (ϕk
k,l=1
(i)
|Ω i
, ϕl )ηk +
µm ϕ(Γ)
m |Ω
m=1
l=1
(i)
i
NΓ
X
(Γ)
µm ai (ϕk
k,m=1
(Γ)
|Ωi , ϕm |Ωi )ηk
(i)
= wT Aii u + µT AΓi u + wT AiΓ η + µT AΓΓ η.
Ma, usando la (4.68), si ottiene:
(i)
T
T
< Si,h ηh , µh > = −wT AiΓ η − µT AΓi A−1
ii AiΓ η + w AiΓ η + µ AΓΓ η
(i)
T
−1
= µ AΓΓ − AΓi Aii AiΓ η
= µT Σi η.
i
µm ai (ϕj , ϕ(Γ)
m |Ω )uj
j,m=1
j,l=1
+
NΓ
X
NΓ
X
4.4. APPROSSIMAZIONE CON ELEMENTI FINITI
117
Grazie al Teorema 4.5 d’estensione uniforme per il caso discreto ed alla (4.72), si può dimostrare
l’esistenza di due costanti K̂1 , K̂2 > 0, indipendenti da h, tali che
K̂1 < S1,h µh , µh > ≤ < S2,h µh , µh > ≤ K̂2 < S1,h µh , µh >
∀µh ∈ Λh .
(4.74)
Infine, utilizzando il Lemma 4.2 e la (4.74), si può affermare che esistono due costanti K̃1 , K̃2 >
0, indipendenti da h, tali che
∀µ ∈ RNΓ .
(4.75)
K̃1 µT Σ1h µ ≤ µT Σ2h µ ≤ K̃2 µT Σ1h µ
4.4.3 Equivalenza tra il metodo di Dirichlet-Neumann e il metodo di Richardson precondizonato: il caso algebrico
Dimostriamo ora il risultato di equivalenza del Teorema 4.4 nel caso algebrico. L’approssimazione con elementi finiti dei problemi (4.43)-(4.44) si scrive in forma matriciale come
(k)
A11 u1 = f1 − A1Γ λ(k−1) ,
A22 A2Γ
(2)
AΓ2 AΓΓ
"
(k)
u2
λ(k−1/2)
#
=
(4.76)
f2
(k)
(1)
fΓ − AΓ1 u1 − AΓΓ λ(k−1)
,
(4.77)
mentre la (4.45) diventa
λ(k) = θλ(k−1/2) + (1 − θ)λ(k−1) .
(4.78)
(k)
Eliminiamo u2 dalla (4.77) usando il metodo di Gauss:
(2)
(k)
(1)
−1
AΓΓ − AΓ2 A22 A2Γ λ(k−1/2) = fΓ − AΓ1 u1 − AΓΓ λ(k−1) − AΓ2 A−1
22 f2 .
Grazie alla definizione (4.62) di Σ2 ed alla (4.76), si ha
Σ2 λ
(k−1/2)
= fΓ −
AΓ1 A−1
11 f1
−
AΓ2 A−1
22 f2
−
(1)
AΓΓ
−
AΓ1 A−1
11 A1Γ
ovvero, usando la definizione (4.62) di Σ1 e la (4.64),
(k−1)
λ(k−1/2) = Σ−1
χ
−
Σ
λ
.
1
Γ
2
Ora, in virtù della (4.78), ricaviamo
(k−1)
λ(k) = θΣ−1
χ
−
Σ
λ
+ (1 − θ)λ(k−1) ,
1
Γ
2
cioè, poiché −Σ1 = −Σ + Σ2 :
(k−1)
(k−1)
λ(k) = θΣ−1
χ
−
Σλ
+
Σ
λ
+ (1 − θ)λ(k−1) ;
2
Γ
2
λ(k−1) ,
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
118
Dunque abbiamo
(k−1)
λ(k) − λ(k−1) = θΣ−1
χ
−
Σλ
.
Γ
2
Σ2 (λ(k) − λ(k−1) ) = θ(χΓ − Σλ(k−1) ),
ovvero il metodo di Richardson per il sistema Σλ = χΓ con precondizionatore Σ2 .
Osservazione 4.1. Il precondizionatore del metodo di Dirichlet-Neumann è sempre l’operatore
di Steklov-Poincaré associato al sotto-dominio su cui si risolve il problema di Neumann. Dunque,
se si risolve il problema di Dirichlet su Ω2 e quello di Neumann su Ω1 , il precondizionatore per
il metodo di Richardson sarà S1 anziché S2 .
Per quanto riguarda il condizionamento del sistema Σλ = χΓ , dalla (4.66) si ha che la matrice
Σ è mal condizionata. L’uso di un precondizionatore per la risoluzione di tale sistema è dunque
fortemente raccomandato.
Nella Sezione 4.6 considereremo il problema della costruzione di un opportuno precondizionatore per l’equazione di Steklov-Poincaré. Possiamo comunque proporre già fin d’ora i precondizionatori Σi , i = 1, 2, per i quali si può dimostrare il seguente risultato:
Teorema 4.6. Se si sceglie come precondizionatore Ph = Σ1 o Ph = Σ2 , allora il numero di
condizionamento spettrale K(Ph−1 Σ) è indipendente da h, ovvero
K(Σ−1
i Σ) = O(1) ,
i = 1, 2.
Corollario 4.4. Il numero d’iterazioni del metodo di Dirichlet-Neumann è indipendente da h.
Osservazione 4.2. Per quel che concerne il metodo di Neumann-Neumann introdotto nella Sezione 4.2.3, si può dimostrare un risultato dello stesso tipo di quello enunciato nel Teorema 4.4,
il che ci garantisce che questo metodo è equivalente al metodo di Richardson per l’equazione di
Steklov-Poincaré con precondizionatore
P −1 = σ1 S1−1 + σ2 S2−1 .
Nel caso algebrico
−1
Ph−1 = σ1 Σ−1
1 + σ 2 Σ2
con σ1 e σ2 coefficienti positivi opportuni. Si può inoltre dimostrare che
−1
K((σ1 Σ−1
1 + σ2 Σ2 )Σ) = O(1)
ovvero che, anche per lo schema di Neumann-Neumann, il numero di iterazioni è indipendente
dal parametro h.
Ricordiamo la definizione seguente:
4.5. GENERALIZZAZIONE AL CASO DI PIÙ SOTTO-DOMINI
119
Definizione 4.2. Sia Ax = b un sistema lineare con A ∈ RN ×N . Un precondizionatore P per
A è detto ottimale se il numero di condizionamento di P −1 A è uniformemente limitato rispetto
alla dimensione N della matrice A.
Possiamo allora concludere che, per la risoluzione del sistema Σλ = χ Γ , si può ricorrere ai
seguenti precondizionatori ottimali:
Ph−1
 −1
Σ
per il metodo di Dirichlet-Neumann,

 2
Σ−1
per il metodo di Neumann-Dirichlet,
=
1


−1
σ1 Σ−1
per il metodo di Neumann-Neumann.
1 + σ 2 Σ2
(4.79)
Anticipiamo il fatto che nel caso di molti sotto-domini tale risultato non risulta più essere valido
nel senso che, se si considera, ad esempio, il metodo di Dirichlet-Neumann, si dovranno risolvere
M/2 problemi di Dirichlet in Ω1 e M/2 problemi di Neumann in Ω2 , ma il precondizionatore
Σ2 non è più ottimale per questo algoritmo. In questo caso il numero di condizionamento della
matrice Σ è così maggiorabile
H
K(Σ) ≤ C
,
(4.80)
2
hHmin
essendo H il diametro massimo dei sotto-domini e Hmin quello minimo.
Ricordiamo che il metodo di Dirichlet-Neumann è equivalente ad uno schema di Richardson con
precondizionatore Σ2 :
Σ2 (λ(k+1) − λ(k) ) = θ(χΓ − Σλ(k) ),
mentre il metodo di Neumann-Neumann è equivalente allo schema di Richardson con precondi−1
−1
−1
zionatore Ph,N
N = (Σ1 + Σ2 )/2 (si rimanda alla (4.79) per σ1 = σ2 = 1/2).
Ora, dalla teoria della convergenza del metodo iterativo di Richardson, sappiamo che la velocità
di convergenza ottimale è data da
ρ=
K(Ph−1 Σ) − 1
,
K(Ph−1 Σ) + 1
da cui, ricordando la (4.75), segue il risultato seguente:
Teorema 4.7. Il numero di condizionamento di Ph−1 Σ è indipendente da h quando Ph = Σ2
ovvero Ph = Ph,N N . Conseguentemente, la velocità di convergenza non dipende da h, cioè
questi due precondizionatori sono ottimali.
4.5 Generalizzazione al caso di più sotto-domini
Vogliamo ora generalizzare i risultati ottenuti nelle sezioni precedenti al caso in cui il dominio
Ω sia suddiviso in un numero M > 2 arbitrario di sotto-domini (vedremo svariati esempi nel
seguito).
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
120
Indichiamo con Ωi , con i S
= 1, . . . , M , gli M sotto-domini senza sovrapposizione tali che
Ω, Γi = ∂Ωi \ ∂Ω e Γ = Γi .
S
Ωi =
A livello differenziale si può generalizzare la formulazione (4.25) del problema di Poisson (4.24)
affermando che è equivalente alla seguente formulazione

−4ui = f in Ωi , i = 1, . . . , M,





 ui = u k
su Γik ,
(4.81)
∂ui
∂uk


=
su
Γ
,

ik

∂ni

 ∂ni
ui = 0
su ∂Ωi ∩ ∂Ω,
con Γik = ∂Ωi ∩ ∂Ωk 6= ∅ e dove ni indica la normale esterna al sotto-dominio Ωi .
A livello discreto, seguendo le idee presentate nella Sezione 4.4 ed indicando con u = (u I , uΓ )T
il vettore delle incognite scomposto in quelle relative ai nodi interni (u I ) e a quelli sull’interfaccia
Γ (uΓ ), si può pervenire alla formulazione algebrica seguente
fI
uI
AII AIΓ
,
(4.82)
=
fΓ
uΓ
AΓI AΓΓ
essendo AΓI = ATIΓ . Osserviamo che ogni sotto-dominio interno non comunica con gli altri, e
dunque AII coincide con la matrice diagonale a blocchi


A11 0 . . .
0
.. 
..

.
. 
 0
AII =  .
.
..
 ..
.
0 
0 . . . 0 AM M
Al contrario AIΓ è una matrice a banda, poiché ci sono delle intersezioni solamente con le
interfacce locali Γi , ma non con le interfacce locali non adiacenti al sotto-dominio considerato.
Su ciascun sotto-dominio ci troveremo dunque a risolvere un problema di Neumann che è rappresentato dalla matrice locale
Aii AiΓ
Ai =
(i)
AΓi AΓΓ ,
con
(Aii )lj = ai (ϕj , ϕl ) 1 ≤ l, j ≤ Ni ,
(i)
(AΓΓ )sr = ai (ψr , ψs ) 1 ≤ r, s ≤ NΓi ,
(AiΓ )lr = ai (ψr , ϕl )
1 ≤ r ≤ N Γi
1 ≤ l ≤ Ni ,
essendo Ni il numero di nodi in Ωi , NΓi quello dei nodi sull’interfaccia Γi , ϕj e ψr le funzioni di
base associate ai nodi interni e di interfaccia, rispettivamente.
La matrice AII è non singolare, e dunque dalla (4.82) possiamo ricavare
uI = A−1
II (fI − AIΓ uΓ ).
(4.83)
4.5. GENERALIZZAZIONE AL CASO DI PIÙ SOTTO-DOMINI
121
Eliminando l’incognita uI del sistema (4.82), si ha
AΓΓ uΓ = fΓ − AΓI A−1
II (fI − AIΓ uΓ ),
ovvero
Ora ponendo
−1
AΓΓ − AΓI A−1
II AIΓ uΓ = fΓ − AΓI AII fI .
(4.84)
Σ := AΓΓ − AΓI A−1
II AIΓ
e
χΓ := fΓ − AΓI A−1
II fI ,
ed introducendo λ = uΓ , la (4.84) diviene
Σλ = χΓ .
(4.85)
L’equazione (4.85) è dunque l’approssimazione agli elementi finiti del problema all’interfaccia
di Steklov-Poincaré nel caso di M sotto-domini.
Un algoritmo generale per risolvere il problema di Poisson su Ω potrà dunque essere così formulato:
1. calcolare la soluzione di (4.85) per ottenere il valore di λ sull’interfaccia Γ;
2. risolvere (4.83) e, poiché AII è una matrice diagonale a blocchi, ciò implica la risoluzione parallela di M problemi indipendenti di dimensione ridotta su ciascun sotto-dominio,
ovvero Aii uiI = gi , i = 1, . . . , M .
4.5.1 Alcuni risultati numerici
Consideriamo il problema modello
(
−∆u = f
u =0
in Ω = (0, 1)2 ,
su ∂Ω,
(4.86)
la cui formulazione ad elementi finiti coincide con la (4.50).
Decomponiamo il dominio Ω in M regioni quadrate Ωi , di dimensione caratteristica H, tali che
∪M
i=1 Ωi = Ω. Un esempio di tale decomposizione basata su 4 sotto-domini è riportato in Figura
4.7 (sulla sinistra).
In Tabella 4.1 riportiamo i valori numerici di K(Σ) relativi al problema in esame. Possiamo
osservare che K(Σ) cresce linearmente con 1/h e con 1/H, come indicato dalla formula (4.80)
per Hmin = H. In Figura 4.7 (sulla destra) riportiamo il pattern di sparsità della matrice Σ
associata alle scelte h = 1/8 e H = 1/2. La matrice ha una struttura a blocchi, che tiene
conto delle interfacce Γ1 , Γ2 , Γ3 e Γ4 , più il contributo dovuto al punto di intersezione Γc delle
quattro interfacce. Si osservi che Σ è densa. Per questo motivo, quando si utilizzano metodi
iterativi per risolvere il sistema (4.85), non è conveniente, dal punto di vista dell’occupazione
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
122
K(Σ)
h=1/8
h=1/16
h=1/32
h=1/64
H = 1/2 H = 1/4 H = 1/8
9.77
14.83
25.27
21.49
35.25
58.60
44.09
75.10
137.73
91.98
155.19
290.43
Tabella 4.1: Numero di condizionamento del complemento di Schur Σ.
0
Γ1
5
PSfrag replacements
Γ4
PSfrag replacements
Γ3
Γ1
Γc
Γ2
Γ3
Γ1
Γ4
Γc
Γ2
10
Γ2
Γ
15
c
Γ
3
20
Γ4
25
30
0
5
10
15
nz = 597
20
25
30
Figura 4.7: Esempio di decomposizione del dominio Ω = (0, 1) 2 in quattro sotto-domini quadrati (a
sinistra). Pattern di sparsità del complemento di Schur Σ (sulla destra) associato alla decomposizione di
domini riportata sulla sinistra.
di memoria, calcolare in modo esplicito gli elementi di Σ. Al contrario, usando l’Algoritmo 4.1
è possibile calcolare il prodotto ΣxΓ , per ogni vettore xΓ . In tale algoritmo abbiamo indicato
con RΓi : Γ → Γi = ∂Ωi \ ∂Ω un opportuno operatore di restrizione, mentre x ← y indica
l’operazione x = x + y.
Algoritmo 4.1 (applicazione del complemento di Schur). Dato xΓ , calcolare yΓ =
ΣxΓ nel modo seguente:
a. Porre yΓ = 0
b. For i = 1, . . . , M Do in parallelo:
c.
x i = R Γ i xΓ
d.
zi = AiΓi xi
e.
zi ← A−1
ii zi
f.
sommare nel vettore globale yΓi ← AΓi Γi xi − AΓi i zi
g.
sommare nel vettore globale yΓ ← RΓTi yΓi
h. EndFor
Si osservi che l’Algoritmo 4.1 è completamente parallelo, nel senso che non è richiesta alcuna
comunicazione tra i vari sotto-domini.
Prima di utilizzare per la prima volta il complemento di Schur, è necessario avviare una fase di
startup descritta nell’Algoritmo 4.2. Questo è un calcolo off-line.
4.6. PRECONDIZIONATORI NEL CASO DI PIÙ SOTTO-DOMINI
123
Algoritmo 4.2 (fase di startup per la risoluzione del sistema associato al complemento di
Schur). Dato xΓ , calcolare yΓ = ΣxΓ nel modo seguente:
a. For i = 1, . . . , M Do in parallelo:
b.
Costruire la matrice Ai
c.
Riordinare Ai come
Aii AiΓi
Ai =
A Γi i A Γi Γi
d.
e.
ed estrarre le sotto-matrici Aii , AiΓi , AΓi i e AΓi Γi
Calcolare la fattorizzazione LU di Aii
EndFor
4.6 Precondizionatori nel caso di più sotto-domini
Vediamo come si possono ottenere dei precondizionatori scalabili per il caso in cui il dominio Ω
sia partizionato in più sotto-domini. Forniamo dapprima la seguente
Definizione 4.3. Si definisce scalabile un precondizionatore che permette di ottenere una matrice
precondizionata Ph−1 Σ con un numero di condizionamento indipendente dal numero di sottodomini.
I metodi iterativi con precondizionatori scalabili consentono velocità di convergenza indipendenti
dal numero dei sotto-domini. Tale proprietà è senza dubbio auspicabile nel momento in cui si
ricorra ad un calcolo parallelo su più processori.
Sappiamo che ad un problema globale su Ω si associa una matrice globale A e che a ciascun
problema su un sotto-dominio Ωi con i = 1, . . . , M , corrisponde una matrice locale
Aii AiΓ
Ai =
.
(i)
AΓi AΓΓ
Definiamo un operatore di restrizione Ri che associa ad un vettore vh del dominio globale Ω la
sua restrizione al sotto-dominio Ωi :
i
Ri : vh|Ω → vh|
.
Ω ∪Γ
i
Sia inoltre RiT ,
i
i
RiT : vh|
→ v h |Ω
Ω ∪Γ
i
i
l’operatore d’estensione (o prolungamento) a zero di vhi . In forma algebrica Ri può essere
rappresentato da una matrice che coincide con la matrice identità in corrispondenza del sottodominio Ωi a cui essa è associata:


0 ... 0
0 ... 0 1

.. . . ..  .
..
Ri =  ... . . . ...
.
. . 
.
1 0 ... 0
0 ... 0
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
124
Ω
e3
e2
e1
v1
e6
e8
v3
Ωi
e4
v2
e5
e7
e9
v4 e10
e12
e11
Figura 4.8: Esempio di decomposizione in pù sotto-domini.
|
{z
Ωi
Analogamente possiamo definire il precondizionatore
Ph =
M
X
}
RΓTi Σi RΓi ,
i=1
che agisce sui vettori i cui valori sono associati ai soli nodi della griglia che vivono sull’unione
delle interfacce.
L’idea seguita per costruire un precondizionatore per il complemento di Schur Σ della matrice
globale è di combinare opportunamente i contributi dovuti ai precondizionatori locali, ovvero che
si possono costruire sui singoli sotto-domini, in modo da fornire un contributo globale ottenuto
a partire da una griglia più grossolana (coarse) di Ω, ad esempio quella, di ampiezza H, i cui
elementi sono i sotto-domini stessi di Ω. Possiamo esprimere questa idea attraverso la seguente
relazione
M
X
−1
RΓi + RΓT PH−1 RΓ .
(Ph )−1 =
RΓTi Pi,h
i=1
Abbiamo indicato con H la massima ampiezza dei diametri Hi dei sotto-domini Ωi . Esistono
diverse scelte possibili per il precondizionatore Pi,h che daranno luogo ad andamenti differenti
per quel che concerne la convergenza dell’algoritmo adottato per risolvere il problema.
4.6.1 Il precondizionatore di Jacobi
Sia {e1 , . . . , em } l’insieme dei lati e {v1 , . . . , vn } l’insieme dei vertici di una partizione del dominio Ω (si veda la Figura 4.8).
4.6. PRECONDIZIONATORI NEL CASO DI PIÙ SOTTO-DOMINI
125
Per la matrice Σ possiamo allora fornire una rappresentazione a blocchi della forma seguente:


Σee Σev


,
Σ=


ΣTev Σvv
con

Σe1 e1 . . . Σe1 em


..
..
Σee =  ...
,
.
.
Σem e1 . . . Σem em

Σev
e

Σe1 v1 . . . Σe1 vn


..
..
=  ...

.
.
Σem v1 . . . Σem vn


Σv1 v1


Σvv = 

0
..
.
0
0 ...
0
..
..
.
.
..
.
0
. . . 0 Σ vn vn



.

Il precondizionatore di Jacobi del complemento di Schur Σ è dunque definito come
Σ̂ee 0
J
Ph =
0 Σvv
dove Σ̂ee coincide con la stessa Σee o con una sua opportuna approssimazione. Tale precondizionatore non tiene conto in alcun modo delle interazioni tra i lati e i vertici della griglia essendo
caratterizzato da una matrice che è diagonale. Inoltre anche la matrice Σ̂ee è diagonale essendo
data da


Σ̂e1 e1 0 . . .
0


..
..

 0
.
.
,

Σ̂ee =  .

..
.
.
 .
0 
0
. . . 0 Σ̂em em
con Σ̂ek ek coincidente con Σek ek o con una sua opportuna approssimazione.
Si può dimostrare il risultato seguente:
Teorema 4.8.
K
con H = maxi (diam(Ωi )) .
(PhJ )−1 Σ
≤ CH
−2
H
1 + log
h
2
,
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
126
Il numero delle iterazioni necessarie per la convergenza del metodo del Gradiente Coniugato è dunque proporzionale ad H −1 , e la presenza di H fa perdere la scalabilità desiderata per
l’algoritmo.
Osserviamo inoltre che la presenza del termine log(H/h) introduce un legame tra la dimensione
di ciascun sotto-dominio e quella degli elementi della griglia computazionale T h . Ciò genera
una propagazione dell’informazione tra i vari sotto-domini caratterizzata da una velocità finita
anziché infinita.
Il precondizionatore PhJ può essere anche espresso in funzione degli operatori di restrizione e di
prolungamento nel seguente modo:
PhJ
−1
=
m
X
T −1
ReTk Σ̂−1
ek ek Rek + Rv Σvv Rv ,
(4.87)
k=1
dove Rek e Rv sono gli operatori di restrizione sui lati e sui vertici, rispettivamente.
4.6.2 Il precondizionatore di Bramble-Pasciak-Schatz
Per accelerare la velocità di propagazione delle informazioni all’interno del dominio Ω, bisogna
ricorrere ad un accoppiamento globale. Dopo aver decomposto Ω in sotto-domini, si può considerare questa stessa decomposizione come una griglia grossolana T H del dominio. In Figura 4.9,
per esempio, abbiamo una griglia coarse composta da 9 elementi e 4 nodi interni. Si può dunque
associare a tale griglia una matrice di rigidezza AH (nel caso in Figura 4.9 sarà una matrice di
dimensione 4 × 4).
Possiamo introdurre inoltre un operatore di restrizione RH : Γh → ΓH , definito sui nodi delle
T
interfacce Γh a valori sui nodi interni della griglia grossolana. Il suo trasposto R H
è sempre l’operatore di estensione. La matrice AH è quella che permette un accoppiamento globale all’interno
del dominio Ω.
Possiamo allora definire il precondizionatore seguente, detto di Bramble-Pasciak-Schatz:
(PhBP S )−1
=
m
X
T −1
ReTk Σ̂−1
ek ek R ek + R H AH R H .
(4.88)
k=1
A differenza del precondizionatore di Jacobi (4.87), nel secondo addendo non compare più una
matrice legata ai vertici (che non fornisce alcun accoppiamento interno, dal momento che ciascun
vertice comunica soltanto con se stesso rendendo Σvv diagonale), ma la matrice globale AH .
Per il precondizionatore di Bramble-Pasciak-Schatz valgono i seguenti risultati:
i) nel caso 2D
K
ii) nel caso 3D
(PhBP S )−1 Σ
H
≤ C 1 + log
h
H
K (PhBP S )−1 Σ ≤ C .
h
2
;
4.6. PRECONDIZIONATORI NEL CASO DI PIÙ SOTTO-DOMINI
127
PSfrag replacements
Figura 4.9: Esempio di griglia grossolana su Ω.
Come si può osservare è stata eliminata la dipendenza della velocità di convergenza dal parametro H. Il numero di iterazioni del metodo del Gradiente Coniugato è dunque proporzionale a
log(H/h) in 2D.
4.6.3 Il precondizionatore di Neumann-Neumann
Il precondizionatore di Bramble-Pasciak-Schatz introduce un miglioramento sul numero di condizionamento della matrice associata al problema di Schur precondizionato, tuttavia nel caso 3D
si ha ancora una dipendenza da H e da h.
Si consideri allora il precondizionatore seguente, detto di Neumann-Neumann:
(PhN N )−1
=
M
X
RΓTi Di Σ−1
i Di R Γ i ,
(4.89)
i=1
essendo RΓi ancora l’operatore di restrizione di Γ ai valori sulle interfacce locali Γi . Sia invece
Di


d1


..
Di = 
,
.
dn
una matrice diagonale di pesi positivi, con dj > 0, per i = 1, . . . , n, essendo n il numero di nodi
su Γi ; per la precisione, dj coincide con l’inverso del numero di sotto-domini che condividono lo
stesso nodo j-esimo.
Esempio 4.5. Se si considerano i quattro nodi interni in Figura 4.9, si avrà d j = 1/4.
•
Teorema 4.9. Per il precondizionatore (4.89) si ha
K
(PhN N )−1 Σ
≤ CH
−2
H
1 + log
h
2
.
L’introduzione di Di e RΓi in (4.89) implica solamente dei prodotti tra matrici. Al contrario si
vorrebbe poter esprimere Σ−1
i in funzione delle matrici Ai . È possibile far ciò nel modo seguente:
sia q il vettore delle incognite su ciascuna interfaccia locale Γi ; allora
T
−1
Σ−1
i q = [0, I]Ai [0, I] q.
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
128
In particolare, [0, I]T q = [0, q]T , e





|
nodi
interni
bordo
{z
A−1
i


0
  .. 
 . 
 
 0 
q
}
corrisponde alla risoluzione sul sotto-dominio Ωi del problema di Neumann

 −4wi = 0 in Ωi ,
 ∂wi = q
∂n
sur Γi .
(4.90)
Ora, poiché Σ−1
associa un dato di Dirichlet sull’interfaccia ad un dato di Neumann, bisogna
i
considerare soltanto i valori della soluzione wi del problema (4.90) sull’interfaccia. Si moltiplica
dunque tale vettore per [0, I], la qual cosa annulla le sue componenti associate all’interno di Ω i ,
conservando soltanto quelle relative all’interfaccia locale Γi .
Si osservi anche che se Ωi è un sotto-dominio interno ad Ω, allora Γi = ∂Ωi . In tal caso il
problema di Neumann (4.90) non risulta ben posto. Bisognerà rimpiazzarlo, per esempio, con il
problema −(4 + αI)wi = 0 in Ωi , con α coefficiente piccolo opportunamente scelto.
Algoritmo 4.3 (precondizionatore Neumann-Neumann). Dato il vettore rΓ , calcolare zΓ = (PhN N )−1 rΓ nel modo seguente:
a. Porre zΓ = 0
b. For i = 1, . . . , M Do in parallelo:
c.
restringere il residuo su Ωi : ri = RΓi rΓ
T
d.
zi = [0, I]A−1
i [0, ri ]
e.
Aggiungere nel residuo globale: zΓ ← RΓTi zi
f. EndFor
Anche in questo caso è richiesta una fase di startup: questa consiste semplicemente nel preparare
il codice per la risoluzione del sistema lineare con matrice dei coefficienti A i . Si osservi che, nel
caso del nostro problema modello (4.86), la matrice Ai non è definita se ∂Ωi \ ∂Ω = ∅. Una delle
seguenti strategie deve essere dunque adottata:
1. calcolare la fattorizzazione LU di Ai +I, con > 0 opportunamente piccolo ed assegnato;
1
2. calcolare la fattorizzazione LU di Ai + 2 Mi , dove Mi è la matrice di massa le cui entrare
H
sono date da
Z
Mi (k, j) =
ϕk ϕj dΩi ;
Ωi
4.6. PRECONDIZIONATORI NEL CASO DI PIÙ SOTTO-DOMINI
129
Convergence history using PNN
0
10
−1
10
−2
scaled residual
10
−3
10
H=1/8
−4
10
H=1/4
−5
10
H=1/2
−6
10
−7
10
0
2
4
6
8
10
Conjugate gradient iterations
12
14
16
Figura 4.10: Storia di convergenza per il metodo del gradiente coniugato precondizionato con la matrice
PhN N e per h = 1/32.
3. calcolare la decomposizione ai valori singolari di Ai ;
4. utilizzare un solutore iterativo per calcolare A−1
i ri , per un certo ri .
Per i nostri risultati numerici abbiamo adottato il terzo approccio. La storia di convergenza per
K((PhN N )−1 Σ)
h = 1/16
h = 1/32
h = 1/64
h = 1/128
H = 1/2 H = 1/4 H = 1/8 H = 1/16
2.55
15.20
47.60
3.45
20.67
76.46
194.65
4.53
26.25
105.38
316.54
5.79
31.95
134.02
438.02
Tabella 4.2: Numero di condizionamento per la matrice precondizionata (P hN N )−1 Σ.
il metodo del gradiente coniugato precondizionato con la matrice P hN N in corrispondenza della
scelta h = 1/32 è riportata in Figura 4.10. Osserviamo che il precondizionatore NeumannNeumann non ha quasi nessun effetto nel caso in cui il numero M dei sotto-domini sia grande. In Tabella 4.2 riportiamo infine i valori del numero di condizionamento per la matrice
precondizionata (PhN N )−1 Σ in corrispondenza di scelte differenti per H.
4.6.4 Il precondizionatore di Neumann-Neumann con bilanciamento
Osserviamo che nella definizione del precondizionatore di Neumann-Neumann non compare
alcun termine globale. Introducendo un termine di tale natura, si ottiene il cosiddetto precondizionatore di Neumann-Neumann con bilanciamento:
−1
N N −1
T −1
T −1
PhBN N
= [(I −RΓT A−1
) +(RΓT A−1
H RΓ Σ)(Ph
H RΓ )](I −ΣRΓ AH RΓ )+RΓ AH RΓ . (4.91)
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
130
Per passare da λ(k) a λ(k+1) utilizzando tale precondizionatore, si devono seguire i passi seguenti:
i)
(k)
λ(k+1/3) = λ(k) + RΓT A−1
H RΓ (χΓ − Σλ )
ii)
λ(k+2/3) = λ(k+1/3) + (PhN N )−1 (χΓ − Σλ(k+1/3) )
(k+2/3)
iii) λ(k+1) = λ(k+2/3) + RΓT A−1
)
H RΓ (χΓ − Σλ
Teorema 4.10. Sia in 2D che in 3D si ha
K
(PhBN N )−1 Σ
H
≤ C 1 + log
h
2
.
Si è ottenuto un miglioramento per il numero di condizionamento rispetto al precondizionatore
di Bramble-Pasciak-Schatz per il quale, nel caso 3D, si aveva ancora una dipendenza da H.
Tuttavia anche il precondizionatore di Neumann-Neumann con bilanciamento non garantisce la
scalabilità ottimale propria dell’algoritmo di Dirichlet-Neumann, ma soltanto una dipendenza
logaritmica da H e h.
Il precondizionatore di Neumann-Neumann con bilanciamento aggiunge una matrice coarse A 0
alle correzioni locali introdotte dal precondizionatore di Neumann-Neumann. La matrice A 0
viene costruita utilizzando l’Algoritmo 4.4.
Algoritmo 4.4 (costruzione della matrice coarse per PhBN N ). Dato il vettore rΓ ,
calcolare zΓ = (PhN N )−1 rΓ nel modo seguente:
a. Costruire l’operatore di restrizione R̄0 che restituisce, per
ogni sotto-dominio, la somma pesata dei valori di tutti i
nodi sul bordo di quel sotto-dominio.
I pesi sono determinati dall’inverso del numero di sotto-domini che contengono ogni nodo i
b. Costruire la matrice A0 = R̄0 S R̄0T
Il passo a. dell’Algoritmo 4.4 è computazionalmente molto economico. Al contrario il passo
b. richiede molti (ad esempio, `) prodotti matrice-vettore con il complemento di Schur Σ. Ciò
significa, nel caso in cui Σ non sia costruita esplicitamente, dover risolvere ` × M problemi di
Dirichlet per generare la matrice A0 . Osserviamo inoltre che l’operatore di restrizione introdotto
al passo a. definisce uno spazio coarse le cui funzioni di base sono costanti a tratti su ciascun
Γi . Per questo motivo, il precondizionatore di Neumann-Neumann con bilanciamento è spesso
una scelta obbligata quando o la griglia o la decomposizione in sotto-domini (o entrambe) è non
strutturata (si veda, ad esempio, la Figura 4.11).
Confrontando i risultati numerici ottenuti per il complemento di Schur e il precondizionatore di
Neumann-Neumann, con e senza bilanciamento, possiamo trarre le seguenti conclusioni:
• anche se meglio condizionata rispetto ad A, Σ è mal-condizionata, e perciò si deve ricorrere
ad un opportuno precondizionatore;
4.7. I METODI ITERATIVI DI SCHWARZ
131
Plotmesh
1
0.8
0.6
0.4
Y − Axis
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
X − Axis
0.5
1
Figura 4.11: Esempio di decomposizione non strutturata in 8 sotto-domini.
K((PhBN N )−1 Σ)
h = 1/16
h = 1/32
h = 1/64
h = 1/128
H = 1/2 H = 1/4 H = 1/8 H = 1/16
1.67
1.48
1.27
2.17
2.03
1.47
1.29
2.78
2.76
2.08
1.55
3.51
3.67
2.81
2.07
Tabella 4.3: Numero di condizionamento di (P hBN N )−1 Σ al variare di H.
• il precondizionatore di Neumann-Neumann può essere applicato in modo soddisfacente utilizzando un numero relativamente piccolo di sotto-domini (ad esempio, fino a 16),
mentre, per grandi valori di M , K((PhN N )−1 Σ) > K(Σ);
• il precondizionatore di Neumann-Neumann con bilanciamento è quasi ottimale e scalabile,
e dunque dovrebbe esser una scelta obbligata nel caso in cui il numero dei sotto-domini sia
grande.
4.7 I metodi iterativi di Schwarz
Il metodo di Schwarz è stato in origine proposto come schema iterativo per dimostrare l’esistenza
della soluzione delle equazioni ellittiche definite su domini a forma di L (si veda la Figura 4.12).
Sebbene venga ancora utilizzato come metodo di soluzione per tali equazioni su domini di forma
generica, è oggigiorno più comunemente impiegato come precondizionatore di tipo domain decomposition per lo schema del gradiente coniugato (o per i metodi di Krylov) nella risoluzione
dei sistemi algebrici derivanti dalla discretizzazione di problemi ai valori al bordo.
Riprendiamo l’idea generale già esposta nella Sezione 4.2.1. Una caratteristica essenziale del
metodo di Schwarz è che si basa su di una suddivisione del dominio computazionale in sottodomini con sovrapposizione. Tuttavia, per semplificare le notazioni, indicheremo nel seguito
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
132
Ω1
Γ
12
Ω2
Figura 4.12: Esempio di dominio a forma di L.
ancora con Ωm il sotto-dominio “esteso”. Inoltre, per semplicità di esposizione, consideriamo
come problema di riferimento il problema di Poisson (4.24).
Per incominciare, supponiamo che Ω sia un dominio bidimensionale e siano Ω 1 ed Ω2 due sottodomini con sovrapposizione la cui unione coincide con Ω stesso. Denotiamo allora con Γ 1 =
∂Ω1 \ ∂Ω e con Γ2 = ∂Ω2 \ ∂Ω.
Il metodo classico di Schwarz prevede la generazione di una successione di funzioni convergente
(0)
alla soluzione del problema (4.24) ed è definito nel seguente modo: data û 2 , per k ≥ 0, risolvere

(k+1)
−4û1
= f
in Ω1 ,



(k+1)
(k)
û1
= û2 su Γ1 ,



(k+1)
û1
= 0
su ∂Ω1 \ Γ1 ,

(k+1)

−4û2
= f



( (k+1)


û1
(k+1)
û
=
2

(k)

û1



(k+1)

û2
= 0
(4.92)
in Ω2 ,
su Γ2 ,
(4.93)
su ∂Ω2 \ Γ2 .
Le funzioni ûi sono definite solamente su Ωi . Come si osserva, in (4.93) vengono contemplate
(k+1)
(k+1)
due alternative: la scelta û2
= û1
porta ad una versione di tipo sequenziale, detta molti(k+1)
(k)
plicativa, mentre ponendo û2
= û1 individuiamo una versione parallela detta additiva. Il
motivo dietro tale terminologia è fornito nell’Osservazione 4.3.
4.7.1 Interpretazione variazionale del metodo di Schwarz
Al metodo (4.92)-(4.93) può essere data un’interpretazione variazionale introducendo opportuni
operatori di proiezione. Tale interpretazione, oltre ad essere essenziale per effettuare un’analisi
di convergenza del metodo, suggerirà il modo per generalizzare il metodo di Schwarz a for-
4.7. I METODI ITERATIVI DI SCHWARZ
133
me più interessanti. Ricordiamo che il problema in esame (4.24) ammette come formulazione
variazionale la (4.50), rimpiazzati uh e vh con u e v.
Facciamo notare inltre che il teorema che stiamo per dimostrare è molto tecnico e non essenziale
per la comprensione degli aspetti algoritmici che vedremo nel seguito.
Teorema 4.11. Un’iterazione moltiplicativa di Schwarz può essere interpretata come un’iterazione del metodo di Richardson per il problema
Qm u = g m .
(4.94)
Analogamente, un’iterazione additiva di Schwarz è equivalente ad un’iterazione del metodo di
Richardson per il problema
Qa u = g a ,
(4.95)
con
Q m = P 1 + P 2 − P 2 P1 ,
Qa = P 1 + P 2 ,
gm e ga termini noti dipendenti da f , e P1 e P2 operatori di proiezione opportuni.
Dimostrazione. Siano
V = H01 (Ω),
Vi0 = H01 (Ωi ) e
Vi∗ = {v ∈ V | v = 0 in Ω \ Ωi }.
Definiamo dunque, per il caso moltiplicativo,
( (k+1)
û1
in Ω1 ,
u(k+1/2) =
(k)
û2 |Ω2 \Ω1 in Ω \ Ω1 ,
u
(k+1)
=
(
(k+1)
û2
(k+1)
û1
in Ω2 ,
|Ω1 \Ω2 in Ω \ Ω2 .
(k)
Ne segue che u(k+1/2) = u(k) + w̃1 , essendo
( (k)
w1 in Ω1 ,
(k)
w̃1 =
0
in Ω \ Ω1 .
(k)
La funzione w1 ∈ V10 soddisfa su Ω1 il seguente problema di Dirichlet:
(k)
aΩ1 (w1 , v1 ) = (f, v1 )Ω1 − aΩ1 (u(k) , v1 )
(k)
Inoltre, u(k+1) = u(k+1/2) + w̃2 con
(k)
w̃2
=
(
(k)
w2
in Ω2 ,
0
in Ω \ Ω2 ,
∀v1 ∈ V10 .
(4.96)
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
134
(k)
e w2 ∈ V20 soluzione su Ω2 del problema di Dirichlet:
(k)
aΩ2 (w2 , v2 ) = (f, v2 )Ω2 − aΩ2 (u(k+1/2) , v2 )
∀v2 ∈ V20 .
(4.97)
Dalla (4.96), per ogni v ∈ V1∗ , otteniamo
(k)
(k)
(k)
a(u(k+1/2) − u(k) , v) = a(w̃1 , v) = aΩ1 (w̃1 , v) = aΩ1 (w1 , v|Ω1 )
= (f, v|Ω1 )Ω1 − aΩ1 (u(k) , v|Ω1 ) = (f, v)Ω − a(u(k) , v)
= a(u, v) − a(u(k) , v) = a(u − u(k) , v).
Analogamente, partendo dalla (4.97), si ha
a(u(k+1) − u(k+1/2) , v) = a(u − u(k+1/2) , v).
Per i = 1, 2, sia ora Pi∗ : V → Vi∗ l’operatore di proiezione ortogonale rispetto al prodotto
scalare a(., .), cioè
∀w ∈ V,
Pi∗ w ∈ Vi∗ :
a(Pi∗ w − w, v) = 0 ∀v ∈ Vi∗ .
(4.98)
Abbiamo allora che
u(k+1/2) − u(k) = P1∗ (u − u(k) )
(4.99)
u(k+1) − u(k+1/2) = P2∗ (u − u(k+1/2) ).
(4.100)
e
Sia ora Ji : Vi∗ → V un operatore di iniezione e defiamo Pi come Pi = Ji Pi∗ : V → V . Sia
inoltre G : H −1 (Ω) → V definito nel modo seguente: Gg = ψ se e solo se
(
−4ψ = g in Ω
ψ = 0 su ∂Ω.
Quindi u = Gf . Dalla (4.99) e dalla (4.100), ricaviamo
u(k+1/2) = u(k) + P1∗ (u − u(k) ) = (I − P1 )u(k) + P1 Gf
= u(k) + P1 G(f − Lu(k) ),
u(k+1) = (I − P2 )u(k+1/2) + P2 Gf = u(k+1/2) + P2 G(f − Lu(k+1/2) ),
(4.101)
(4.102)
avendo posto L = −4. Dunque
u(k+1) = (I − P2 )[(I − P1 )u(k) + P1 Gf ] + P2 Gf
= (I − P2 )(I − P1 )u(k) + (I − P2 )P1 Gf + P2 Gf,
ed otteniamo
u(k+1) = u(k) + Qm (Gf − u(k) ).
(4.103)
4.7. I METODI ITERATIVI DI SCHWARZ
Posto gm = Qm Gf , si ha
135
u(k+1) − u(k) = gm − Qm u(k) .
(4.104)
Si osservi che il termine di destra di tale equazione è il residuo associato al problema Q m u = gm .
Ciò prova la (4.94).
In modo del tutto analogo, possiamo dimostrare che un’iterazione additiva di Schwarz può essere
scritta come
u(k+1) − u(k) = ga − Qa u(k)
(4.105)
con ga = Qa Gf , il che prova la (4.95).
Per definizione
ga = Qa Gf = P1 Gf + P2 Gf.
Posto gi = Pi Gf , per i = 1, 2 e sfruttata la definizione di Pi , otteniamo
∀vi ∈ Vi∗
a(gi , vi ) = a(Gf, vi ) = a(u, vi ) = (f, vi )
(i = 1, 2).
Il calcolo di ga richiede dunque la soluzione di due problemi di Dirichlet indipendenti, su Ω 1 ed
Ω2 , rispettivamente.
Osservazione 4.3. Gli aggettivi moltiplicativa ed additiva nascono, rispettivamente, dalla struttura moltiplicativa di Qm ed additiva di Qa .
Per studiare la convergenza delle iterazioni di Schwarz, introduciamo l’errore e k = u − uk
associato alla k-esima iterata. Utilizzando la (4.104), troviamo la seguente formula ricorsiva
ek+1 = (I − P2 )(I − P1 )ek ,
k≥0
per la versione moltiplicativa, mentre dalla (4.105) si ha che
ek+1 = (I − P1 − P2 )ek ,
k≥0
per il caso additivo.
Utilizzando queste formule ricorsive e il principio del massimo per gli operatori di proiezione, si
può verificare che, in entrambi i casi,
kek kL∞ (Ω) ≤ C(δ)2k ke0 kL∞ (Ω)
per un’opportuna costante 0 < C(δ) < 1, che diventa tanto più prossima ad 1 quanto più la
misura δ della zona di sovrapposizione si riduce a zero. Nel caso moltiplicativo, ad esempio, se
Ω coincide con un intervallo monodimensionale, abbiamo
C(δ) =
1 − δ/2
,
1 + δ/2
espressione che conferma appunto il vantaggio nell’aumentare la sovrapposizione tra i sottodomini.
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
136
4.7.2 Forma algebrica dei metodi di Schwarz per una discretizzazione ad
elementi finiti
I metodi di Schwarz possono essere adattati alla soluzione dei sistemi algebrici che derivano
dall’approssimazione numerica di problemi ai valori al bordo. Consideriamo il sistema Au = f
che supponiamo generato dalla discretizzazione ad elementi finiti del problema di Poisson.
Supponiamo ancora che il dominio Ω sia decomposto in due sotto-domini, Ω 1 ed Ω2 , con sovrapposizione come mostrato in Figura 4.1 (in basso).
Indichiamo con Nh il numero totale dei nodi interni di Ω, e con n1 e n2 , rispettivamente, i
nodi interni di Ω1 e Ω2 . Osserviamo che Nh ≤ n1 + n2 e che l’uguaglianza vale soltanto se
la sovrapposizione si riduce ad una singola striscia di elementi. Sia I = {1, . . . , N h } l’insieme
degli indici dei nodi di Ω, mentre denotiamo con I1 e I2 quelli associati, rispettivamente, ad Ω1 ed
Ω2 . Abbiamo dunque che I = I1 ∪ I2 , mentre I1 ∩ I2 = ∅ solo nel caso in cui la sovrapposizione
si riduca ad una singola fascia di elementi.
Ordiniamo i nodi nei tre blocchi in modo tale che il primo blocco corrisponda ai nodi di Ω 1 \ Ω2 ,
il secondo a Ω1 ∩ Ω2 , e il terzo a Ω2 \ Ω1 . La matrice di stiffness A contiene due sotto-matrici, A1
ed A2 , che corrispondono, rispettivamente, alle matrici di stiffness locali associate ai problemi
di Dirichlet in Ω1 e Ω2 (si veda la Figura 4.13). Possono essere ottenute nel modo seguente:
n1
PSfrag replacements
A1
Nh
A
A2
n2
Figura 4.13: Le sotto-matrici A1 ed A2 della matrice di stiffness A.
A1 = R1 AR1T ∈ Rn1 ×n1 e A2 = R2 AR2T ∈ Rn2 ×n2 , dove Ri ed RiT , con i = 1, 2, sono gli
operatori di restrizione e di prolungamento la cui rappresentazione matriciale è data da




1 ... 0


 .. . . .. 
0


 .
. . 








Nh ×n1
T
T
0
.
.
.
1
,
R2 =  1 . . . 0  ∈ RNh ×n2 .
R1 = 
∈R




 .. . . .. 


 .


. . 
0
0 ... 1
Se v è un vettore di RNh , allora R1 v è un vettore di Rn1 le cui componenti coincidono con le
prime n1 componenti di v. Se v è invece un vettore di Rn1 , allora R1T v è un vettore di dimensione
Nh le cui ultime Nh − n1 componenti sono nulle.
4.7. I METODI ITERATIVI DI SCHWARZ
137
Sfruttando tali definizioni, possiamo rappresentare un’iterata del metodo di Schwarz, nella versione moltiplicativa, applicata al sistema Au = f nel seguente modo:
k
uk+1/2 = uk + R1T A−1
1 R1 (f − Au ),
(4.106)
k+1/2
uk+1 = uk+1/2 + R2T A−1
),
2 R2 (f − Au
(4.107)
in analogia alle (4.101) e (4.102). Equivalentemente, posto
Pi = RiT A−1
i Ri A,
con i = 1, 2, abbiamo, in analogia con la (4.103),
uk+1/2 = (I − P1 )uk + P1 u,
uk+1 = (I − P2 )uk+1/2 + P2 u = (I − P2 )(I − P1 )uk + (P1 + P2 − P2 P1 )u.
In modo analogo, un’iterata del metodo di Schwarz, nella versione additiva, diventa
T −1
(k)
u(k+1) = u(k) + (R1T A−1
1 R1 + R2 A2 R2 )(f − Au ),
(4.108)
u(k+1) = (I − P1 − P2 )u(k) + (P1 + P2 )u.
(4.109)
ovvero
Se ora introduciamo le matrici
−1
Qi := RiT A−1
i Ri = P i A ,
per i = 1, 2, allora dalle (4.106) e (4.107) deriviamo che
uk+1 = uk + Q1 (f − Auk ) + Q2 [f − A(uk + Q1 (f − Auk ))]
= uk + (Q1 + Q2 − Q2 Q1 )(f − Auk )
per il metodo di Schwarz moltiplicativo, mentre, per il caso additivo, ricaviamo dalla (4.108) che
uk+1 = uk + (Q1 + Q2 )(f − Auk ).
(4.110)
Quest’ultima formula può essere facilmente estesa al caso di una decomposizione di Ω in M ≥ 2
sotto-domini {Ωi } con sovrapposizione (si veda la Figura 4.14), avendosi in tal caso
u
k+1
k
=u +
M
X
i=1
Qi (f − Auk ).
(4.111)
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
138
4.7.3 Il metodo di Schwarz come precondizionatore
Indicato con
Pas =
M
X
Qi
i=1
−1
,
(4.112)
segue immediatamente dalla (4.111) che un’iterata del metodo di Schwarz additivo corrisponde
ad un’iterata dello schema di Richardson precondizionato con la matrice P as applicato alla risoluzione del sistema lineare Au = f . Per questo motivo, Pas è noto come precondizionatore di
Schwarz additivo.
Analogamente, se introduciamo la matrice precondizionata
Qa =
−1
Pas
A
=
M
X
Pi ,
i=1
un’iterata del metodo di Schwarz additivo corrisponde ad un’iterata dello schema di Richardson
−1
applicato al sistema Qa u = ga , con ga = Pas
f.
Lemma 4.3. La matrice Pi è simmetrica e non negativa rispetto al prodotto scalare indotto da
A
(w, v)A := (Aw, v)
∀w, v ∈ RNh .
Dimostrazione. Abbiamo, per i = 1, 2,
T −1
(Pi w, v)A = (APi w, v) = (RiT A−1
i Ri Aw, Av) = (Aw, Ri Ai Ri Av)
= (w, Pi v)A ,
∀v, w ∈ RNh .
Inoltre,
Nh
−1
.
(Pi v, v)A = (APi v, v) = (RiT A−1
i Ri Av, Av) = (Ai Ri Av, Ri Av) ≥ 0. ∀v ∈ R
Lemma 4.4. La matrice precondizionata Qa del metodo di Schwarz additivo è simmetrica e
definita positiva rispetto al prodotto scalare indotto da A.
Dimostrazione. Dimostriamo, per primo, la simmetria: per ogni u, v ∈ R Nh , in virtù della
simmetria di A e Pi , otteniamo
X
(Qa u, v)A = (AQa u, v) = (Qa u, Av) =
(Pi u, Av)
i
=
X
(Pi u, v)A =
i
X
(u, Pi v)A = (u, Qa v)A .
i
Per quel che riguarda la positività, scegliendo nelle precedenti identità u = v, otteniamo:
X
X
X
(Qa v, v)A =
(Pi v, v)A =
(RiT A−1
(A−1
i Ri Av, Av) =
i qi , qi ) ≥ 0
i
i
i
4.7. I METODI ITERATIVI DI SCHWARZ
139
con qi = Ri Av. Segue che (Qa v, v)A = 0 se e solo se qi = 0, per ogni i, ovvero se e solo se
Av = 0; poiché A è definita positiva, ciò si verifica se e solo se v = 0.
Si può dunque ottenere un metodo iterativo più efficiente per risolvere il sistema lineare Au = f
utilizzando lo schema del gradiente coniugato precondizionato con il precondizionatore di Schwarz additivo Pas . Tale precondizionatore non è tuttavia ottimale poiché il numero di condizionamento della matrice precondizionata Qa cresce al ridursi della misura dei sotto-domini. Infatti
si ha:
1
−1
,
K(Pas
A) ≤ C
δH
con C costante indipendente da h, H e δ, ed essendo H = max i=1,...,M {diam(Ωi )}. Ciò è
dovuto al fatto che lo scambio di informazione avviene soltanto tra sotto-domini vicini, in quanto
l’applicazione di (Pas )−1 coinvolge soltanto risolutori locali. Tale limite può essere superato
introducendo un problema globale coarse definito sull’intero dominio Ω e in grado di garantire
una comunicazione globale tra tutti i sotto-domini. Questa correzione ci porterà a considerare
una strategia generale di tipo multi-livello, come avremo modo di vedere nella Sezione 4.7.4.
Per applicare i metodi di decomposizione di domini con sovrapposizione, dividiamo Ω in M
M
sotto-domini {Ωi }M
i=1 tali che ∪i=1 Ωi = Ω. Gli Ωi condividono una sovrapposizione almeno
+
pari a δ = ξh, con ξ ∈ N . In particolare, ξ = 1 corrisponde alla sovrapposizione minimale,
cioè al caso in cui la sovrapposizione si riduce ad una sola fila di elementi. Da un punto di vista
algoritmico, si procede nel seguente modo:
Algoritmo 4.5 (definizione di sotto-domini con sovrapposizione).
a. Costruire una triangolazione Th del dominio computazionale Ω
b. Suddividere la griglia Th in M sotto-domini {Ω̂i }M
i=1 senza sovrapposizione tali che ∪M
i=1 Ω̂i = Ω
c. Estendere ogni sotto-dominio Ω̂i aggiungendo tutte le strisce
di elementi di Ω entro una certa distanza δ da Ω̂i . Sono così
individuati i domini Ωi
Riferiamo alla Figura 4.14 per una rappresentazione di una regione bidimensionale rettangolare
suddivisa in 9 regioni Ω̂i senza sovrapposizione (sulla sinistra) ed un esempio di sotto-dominio
esteso (sulla destra).
Per applicare il precondizionatore di Schwarz (4.112), procediamo come indicato nell’Algoritmo
4.6. Ricordiamo che ni indica il numero dei nodi contenuti in Ωi , RiT ∈ Rn×ni è una matrice
di prolungamento la cui azione è quella di estendere a zero un vettore definito su Ω i , Ri è la sua
trasposta e Ai = Ri ARiT . Riportiamo in Figura 4.15 un esempio di pattern di sparsità per la
matrice Ri .
Algoritmo 4.6 (fase di startup per applicare Pas ).
a. Su ogni sotto-dominio Ωi , costruire Ri ed RiT
b. Costruire la matrice A corrispondente alla discretizzazione
ad elementi finiti sulla griglia Th
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
PSfrag
140
replacements
Ω̂1
Ω̂4
Ω̂7
Ω̂2
Ω̂5
Ω̂8
Ω̂3
Ω̂6
Ω̂9
Ω5
Figura 4.14: Partizione di una regione rettangolare 2D Ω suddivisa in 9 regioni senza sovrapposizione
Ω̂i (sulla sinistra), e un esempio di sotto-dominio esteso Ω 5 (sulla destra).
0
5
10
15
20
25
0
10
20
30
40
nz = 27
50
60
70
80
Figura 4.15: Esempio di pattern di sparsità della matrice R i per una partizione del dominio in 4 sottodomini.
c. Su ogni Ωi , costruire le sotto-matrici locali Ai = Ri ARiT
d. Su ogni Ωi , predisporre il codice alla risoluzione di un
sistema lineare con matrice associata Ai . Per esempio,
si può calcolare la fattorizzazione LU di Ai
Alcuni commenti generali sugli Algoritmi 4.5 e 4.6 sono d’obbligo:
• i passi a. e b. dell’Algoritmo 4.5 possono essere effettuati in ordine inverso, ovvero si
può prima suddividere il dominio computazionale (utilizzando, ad esempio, informazioni
tratte dal problema fisico in esame), e poi triangolare;
• a seconda di quella che è la struttura generale del codice, si potrebbe pensare di unire i
passi b. e c. dell’Algoritmo 4.6 al fine di ottimizzare sia le risorse di memoria che i
tempi di CPU;
• nel passo d. dell’Algoritmo 4.6, la fattorizzazione esatta potrebbe essere rimpiazzata da
una fattorizzazione incompleta.
L’implementazione finale del precondizionatore Pas può essere effettuata nel seguente modo:
Algoritmo 4.7 (precondizionatore di Schwarz ad un livello). Dato un vettore r,
−1
calcolare z = Pas
r nel modo seguente:
4.7. I METODI ITERATIVI DI SCHWARZ
141
a. Porre z = 0
b. For i = 1, . . . , M Do in parallelo:
c.
restringere il residuo su Ωi : ri = Ri r
d.
calcolare zi = A−1
i ri
e.
sommare nel residuo globale: z ← RiT zi
f. EndFor
Osserviamo che nel passo b. si considerano soltanto i nodi interni di Ω i dal momento che stiamo
imponendo condizioni al bordo omogenee su ∂Ωi \∂Ω. Inoltre, nel passo d., il calcolo di zi viene
spesso effettuato in modo approssimato. Per esempio, si può ricorrere ad una fattorizzazione
incompleta di Ai . In alternativa, si potrebbe utilizzare un risolutore iterativo annidato (munito
di un precondizionatore locale), arrestato in corrispondenza di una tolleranza relativamente alta.
Quest’ultima soluzione richiede tuttavia cautela dal momento che il precondizionatore risultante
non è costante ad ogni applicazione, ed il risolutore iterativo che ne deriva potrebbe dunque
divergere.
In Tabella 4.4, analizziamo il caso di una sovrapposizione minimale (δ = h) considerando diversi
valori per il numero M di sotto-domini. Per quel che riguarda la decomposizione di domini,
consideriamo dei quadrati Ωi , con sovrapposizione, ciascuno di area H 2 . Come previsto dalla
−1
teoria, si osserva che K(Pas
A) ≤ C/(Hδ).
−1
K(Pas
A)
h = 1/16
h = 1/32
h = 1/64
h = 1/128
H = 1/2 H = 1/4 H = 1/8 H = 1/16
15.95
27.09
52.08
31.69
54.52
104.85
207.67
63.98
109.22
210.07
416.09
127.99
218.48
420.04
832.57
−1 A al variare di h e H.
Tabella 4.4: Numero di condizionamento di Pas
4.7.4 Metodi di Schwarz a due livelli
Come anticipato in Sezione 4.7.3, il limite principale del metodo di Schwarz sta nel fatto che
l’informazione si propaga soltanto tra i sotto-domini adiacenti. Per ovviarvi, come già precedentemente fatto per il metodo di Neumann-Neumann, si può introdurre un termine di scambio di
informazioni globale all’interno del dominio Ω. L’idea è sempre quella di considerare i sottodomini Ωi come dei macro-elementi su Ω costituenti una griglia coarse TH a cui si può associare
la matrice AH . Si può dunque introdurre l’operatore QH associato alla soluzione ad elementi
finiti globale:
T −1
QH = R H
AH RH ,
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
142
dove RH denota l’operatore di restrizione sulla griglia coarse. Posto Q0 := QH , possiamo allora
definire una nuova matrice di precondizionamento
−1
Pcas
=
M
X
Qi ,
(4.113)
i=0
per la quale è possibile dimostrare il risultato seguente:
Teorema 4.12. Esiste una costante C > 0, indipendente da h ed H, tale che
−1
K(Pcas
A) ≤ C
H
.
δ
−1
Si osserva dunque che, utilizzando il precondizionatore Pcas
, si perviene ad un metodo scalabile
e quindi parallelizzabile.
L’applicazione di Pcas richiede le stesse operazioni necessarie per l’utilizzo di Pas , più le seguenti
altre:
Algoritmo 4.8 (fase di startup per applicare Pcas ).
a. Eseguire l’Algoritmo 4.6
b. Definire una triangolazione grossolana T0 i cui elementi sono
dell’ordine di H, e porre n0 = dim(V0 ). Supponiamo che Th sia
annidata in T0 (si veda, per esempio, la Figura 4.16 o 4.17)
c. Costruire l’operatore di restrizine R0 ∈ Rn0 ×n .
questa matrice sono date da
Le entrate di
R0 (i, j) = Φi (xj ),
dove Φi è la funzione di base associata al nodo i della griglia
grossolana mentre con xj indichiamo le coordinate del nodo j
sulla griglia fine
d. Costruire la matrice coarse A0 . Ciò può esser fatto
discretizzando il problema variazionale originale su T0 ,
ovvero calcolando A0 come
Z X
d
∂Φi ∂Φj
A0 (i, j) = a(Φj , Φi ) =
,
Ω `=1 ∂x` ∂x`
oppure come
A0 = R0 AR0T .
4.7. I METODI ITERATIVI DI SCHWARZ
143
Figura 4.16: Esempio di griglia coarse per un dominio 2D basata su una mesh strutturata. I triangoli
della griglia fine sono in linea sottile mentre le linee spesse identificano i triangoli associati alla mesh
coarse.
Figura 4.17: Esempio di griglia coarse per un dominio 2D basata su una mesh non strutturata. I triangoli
della griglia fine sono in linea sottile mentre le linee rosse identificano i triangoli associati alla mesh
coarse.
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
144
Se il dominio computazionale ha una forma “semplice” (come nel caso del problema che stiamo
considerando), viene solitamente costruita prima la griglia grossolana T 0 e poi, raffinandola un
certo numero di volte, si genera la griglia fine Th . In presenza invece di domini caratterizzati da
una forma complessa o di griglie non uniformi, la generazione della griglia coarse costituisce l’aspetto più problematico degli algoritmi a due livelli. In tal caso, si può costruire prima la griglia
fine Th e poi, deraffinandola opportunamente, si genera la griglia grossolana T 0 . Operativamente
ciò significa che i nodi di T0 sono un sotto-insieme dei nodi di Th . Gli elementi costituenti T0
verranno quindi costruiti a partire da tale insieme di nodi, operazione questa abbastanza semplice
nel caso bi-dimensionale, meno banale nel caso 3D.
Come alternativa si potrebbero generare le due griglia Th e T0 indipendentemente una dall’altra
e poi costruire gli operatori R0 ed R0T (si osservi che, in tal caso, le due mesh potrebbero non
essere annidate). È importante inoltre ricordare che il problema sulla griglia coarse dovrà essere
completato da condizioni al bordo opportune, che tengano conto di quanto assegnato sulla griglia
fine.
L’implementazione finale di Pcas può dunque essere realizzata nel seguente modo:
−1
Algoritmo 4.9 (applicazione di Pcas ). Dato un vettore r, calcolare z = Pcas
r
nel modo seguente:
a.
b.
c.
d.
e.
f.
g.
h.
Porre z = 0
For i = 1, . . . , M Do in parallelo:
restringere il residuo su Ωi : ri = Ri r
calcolare zi = A−1
i ri
sommare nel residuo globale: z ← RiT zi
EndFor
Risolvere z0 = A−1
0 (R0 r)
sommare nel residuo globale: z ← R0T z0
I passi a.-f. dell’Algoritmo 4.9 sono esattamente gli stessi dell’Algoritmo 4.7. Perciò, una
volta che R0 ed A0 sono stati generati, è estremamente facile modificare una function che
implementa Pas includendo la correzione basata sulla griglia coarse.
−1
In Tabella 4.5 riportiamo il numero di condizionamento di Pcas
A. La sovrapposizione è identificata da δ = h. Osserviamo che il numero di condizionamento è pressoché costante per valori
fissati del rapporto H/δ.
−1
K(Pcas
A)
h = 1/32
h = 1/64
h = 1/128
h = 1/256
H = 1/4 H = 1/8 H = 1/16 H = 1/32
7.03
4.94
12.73
7.59
4.98
23.62
13.17
7.66
4.99
45.33
24.34
13.28
-
−1 A al variare di h e H.
Tabella 4.5: Numero di condizionamento di Pcas
4.7. I METODI ITERATIVI DI SCHWARZ
145
È possibile definire lo spazio coarse in modo alternativo a quanto fatto fino ad ora. Supponiamo
che gli elementi identificanti l’operatore R0 di restrizione dalla griglia fine a quella grossolana
siano dati da
(
1 se il nodo j ∈ Ωi ,
R̂0 (i, j) =
0 altrimenti,
e che Â0 = R̂0 AR̂0T . Questa procedura è detta aggregazione (o agglomerato), poiché gli elementi
di Â0 sono costruiti semplicemente sommando gli elementi di A. Tale procedimento non richiede
l’introduzione di una griglia coarse. Il precondizionatore risultante è dato da
−1
Paggre
= R̂0T Â−1
0 R̂0 + Pas .
Si può dimostrare che
−1
K(Paggre
A)
H
≤C 1+
δ
.
Riportiamo in Tabella 4.6 alcuni risultati numerici associati a tale precondizionatore.
−1
Paggre
A
h = 1/16
h = 1/32
h = 1/64
h = 1/128
H = 1/4 H = 1/8 H = 1/16
13.37
8.87
26.93
17.71
9.82
54.33
35.21
19.70
109.39
70.22
39.07
−1 A al variare di h e H.
Tabella 4.6: Numero di condizionamento di Paggre
Occupiamoci ora dei tempi di CPU richiesti dai metodi a due livelli su architetture parallele.
Utilizziamo a tal fine il seguente precondizionatore di Schwarz a due livelli:
−1
−1
Phybrid2 = Pas
+ B0 − Pas
AB0 ,
con B0 = R0T A−1
0 R0 . Questo precondizionatore è non simmetrico; si dovrà perciò ricorrere, in questo caso, ad opportuni risolutori alla Krylov, quali il GMRES, il TFQMR, il CGS, o
il Bi-CGSTAB. Quest’ultimo schema è quello da noi adottato nei test numerici a cui faremo
riferimento.
Il problema di riferimento è
(
−∆u = xey
in Ω = (0, 1)2 ,
u = −xey
su ∂Ω.
Come criterio globale di convergenza per lo schema iterativo abbiamo adottato il seguente:
kAuk − f k2
≤ tol,
kuk k2
(4.114)
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
146
essendo tol ≤ 10−6 la tolleranza richiesta. Sui vari sotto-domini abbiamo risolto il problema
utilizzando il metodo CG precondizionato con RILU. Il criterio locale di convergenza è stato
ancora identificato con la (4.114), scelto però tol = 10−4 . Per la discretizzazione del problema è
stata utilizzata una griglia fine costituita da 480×480 elementi Q1 , ed una griglia coarse costituita
da 60 × 60 elementi Q1 .
Riportiamo:
• il numero di iterazioni richiesto per la convergenza:
Pas
Phybrid2
M =4 M =9
91
112
11
11
M = 16 M = 25 M = 36
132
145
157
11
12
11
• il tempo di CPU in ambiente parallelo:
Pas
Phybrid2
M =4 M =9
558.25 174.89
169.32 40.83
• il corrispondente speed-up1 , definito come
speedup ideale
Pas
Phybrid2
M =4
4
4
4
M = 16 M = 25 M = 36
79.07
45.53
31.96
18.25
11.16
5.87
T4procs ∗ 4
:
Tnprocs
M = 9 M = 16 M = 25 M = 36
9
16
25
36
12.76
28.24
49.05
69.86
16.48
37.05
61.45
115.37
• il numero di iterazioni richiesto per la convergenza fissata una sovrapposizione pari ad h:
Pas
Phybrid2
1
M =4 M =9
91
112
11
11
M = 16 M = 25 M = 36
132
145
157
11
12
11
Lo speedup di un algoritmo è definito come
sp =
Tp
,
T1
dove Tp è il tempo richiesto per risolvere il problema considerato servendosi di p processori, mentre T 1 è l’analogo
tempo servendosi però di un solo processore. A volte viene considerato invece lo speedup relativo, definito come:
s0p =
Tp
,
T` × `
essendo ` il numero minimo di processori utilizzabile per la risoluzione del problema in esame.
4.8. ALCUNE ESTENSIONI
147
• il numero di iterazioni richiesto per la convergenza fissata una sovrapposizione pari ad
Pas
Phybrid2
M =4 M =9
35
51
4
4
1
:
20
M = 16 M = 25 M = 36
67
81
96
5
6
6
Osserviamo che, come anticipato dalla teoria, il numero di iterazioni richiesto per la convergenza
aumenta all’aumentare del numero di sotto-domini se si adotta il precondizionatore P as . Per i
metodi a due livelli abbiamo invece un numero di condizionamento limitato da un fattore costante
indipendente da H ed h, poiché H/δ =costante. In questo caso perciò il precondizionatore a due
livelli è ottimale (nel senso che il numero di condizionamento non dipende da h) e scalabile (nel
senso che il numero di condizionamento non dipende da H).
In base ai risultati presentati per i precondizionatori di Schwarz a uno e a due livelli, possiamo
concludere che:
- i metodi a due livelli sono sempre caratterizzati da performance migliori di quelli ad un
solo livello, anche se, per h piccoli e per un ridotto numero di sotto-domini, P as fornisce
risultati comunque soddisfacenti;
- per valori di M grandi è fondamentale utilizzare metodi a due livelli, di tipo additivo o
ibrido;
- possono essere adottate tecniche di aggregazione per la definizione della griglia coarse nel
caso in cui tale costruzione risulti difficile o computazionalmente onerosa.
4.8 Alcune estensioni
Le tecniche di decomposizione in sotto-domini illustrate nelle Sezioni precedenti per il problema
di Poisson, possono essere applicate, in modo analogo, a problemi qualsiasi, anche evolutivi,
della forma
(
Lu = f
in Ω,
(4.115)
condizioni al bordo su ∂Ω,
o

∂u


+ Lu = f
in Ω × (0, T ),

 ∂t
condizioni al bordo
su ∂Ω × (0, T ),




condizione iniziale
in Ω, con t = 0,
(4.116)
con L operatore differenziale qualsiasi, lineare o non lineare. Identifichiamo, ad esempio, l’operatore L in (4.115) con
Lu = −∇ · (a∇u) + b · ∇u + cu,
(4.117)
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
148
dove a, b e c sono funzioni assegnate opportune. Considerando una decomposizione senza
sovrapposizione, avremo le seguenti condizioni d’interfaccia su Γ:
u1 = u 2 ,
(a∇u1 ) · n = (a∇u2 ) · n.
Nel caso evolutivo, identificato ancora L in (4.116) con
seguenti problemi, per ogni t > 0:

∂u1


+ Lu1 = f



∂t



 ∂u2
+ Lu2 = f
∂t




u1 = u 2




 (a∇u ) · n = (a∇u ) · n
1
2
(4.117), bisognerà risolvere invece i
in Ω1 ,
in Ω2 ,
su Γ,
su Γ.
Bisognerà in questo caso discretizzare la derivata temporale utilizzando un metodo implicito
piuttosto che esplicito. Nel caso in cui si utilizzi, ad esempio, il metodo di Eulero implicito si
avrà
un+1 − un
+ Lun+1 = f
∆t
e poi andranno applicate le condizioni di interfaccia.
4.8.1 Condizioni all’interfaccia per problemi differenti
Introduciamo in questa sezione le condizioni di interfaccia associate alla risoluzione con il
metodo della decomposizione di domini di alcuni problemi interessanti per le applicazioni.
La Tabella 4.7 riassume le condizioni da imporsi sull’interfaccia Γ per alcuni tra i più diffusi
problemi differenziali.
Forniamo qualche spiegazione riguardo ai problemi di trasporto e di Stokes.
Problemi di trasporto. Consideriamo il problema
Lu = ∇ · (bu) + a0 u = f
in Ω,
(4.118)
completato su ∂Ω da condizioni al bordo opportune. È un problema del prim’ordine e dunque
all’interfaccia si avrà una sola condizione. Introduciamo una partizione dell’interfaccia Γ (si
veda la Figura 4.18): Γ = Γin ∪ Γout , dove
Γin = {x ∈ Γ| b(x) · n(x) > 0 }
e Γout = Γ \ Γin . Si può utilizzare allora il seguente metodo di tipo Dirichlet-Neumann (con
rilassamento):
(
Luk+1
=f
in Ω1 ,
1
(b · n)uk+1
= (b · n)uk2 su Γout ,
1
4.8. ALCUNE ESTENSIONI
149
PSfrag replacements
Ω1
Ω
Γin
n
Γ
Γout
Ω2
Figura 4.18: Decomposizione del dominio per il problema di trasporto (4.118).
(
Luk+1
=f
2
in Ω2 ,
(b · n)uk+1
= (b · n)uk1 su Γin .
2
Il problema di Stokes. Osserviamo che per il problema di Stokes ci sono due incognite: la
velocità e la pressione del fluido. Tuttavia all’interfaccia si impone solamente la continuità della
velocità e non quella della pressione. Ciò è dovuto al fatto che p ∈ L 2 (Ω) e che, se una funzione
ϕ appartiene ad L2 (Ω), allora le sue restrizioni ai sotto-domini Ω1 ed Ω2 appartengono ancora
ad L2 (Ω1 ) e L2 (Ω2 ), e viceversa. Al contrario se f ∈ H 1 (Ωi ), per i = 1, 2, allora in generale
f∈
/ H 1 (Ω): questo giustifica l’imposizione della continuità su Γ per la velocità.
Il problema che si considera è il seguente:

−ν4u + ∇p = f in Ω,


∇·u=0
in Ω,


u=0
su ∂Ω.
Quello che si fa è imporre sull’interfaccia Γ la continuità della velocità u e del tensore di Cauchy
∂u
− pn.
∂n
Sia S l’operatore di Stokes. Si può allora applicare uno schema di tipo Dirichlet-Neumann e, in
tal caso, si è ricondotti alla risoluzione dei due seguenti problemi:

k+1

S(uk+1
in Ω2 ,

2 , p2 ) = f



∂uk1
∂uk+1
(4.119)
=
ν
− pk1 su Γ,
ν 2 − pk+1
2

∂n
∂n



 uk+1 = 0
su ∂Ω \ Γ,
ν
2
2

k+1

S(uk+1
in Ω1 ,
1 , p1 ) = f


k+1
k+1
u1 = θu2 + (1 − θ)uk1 su Γ,


 uk+1 = 0
su ∂Ω1 \ Γ.
1
(4.120)
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
150
Osserviamo che, se si considera come condizione al bordo per la velocità u = 0, la pressione
p non è univocamente determinata, ovvero è definita a meno
R di una costante additiva. Bisogna
dunque fissare tale costante, imponendo, ad esempio, che Ω p = 0.
Ora la questione è: come si può garantire che tale vincolo venga rispettato quando k tende
all’infinito in (4.119) e (4.120)?
Quando si risolve il problema di Neumann (4.119) sul sotto-dominio Ω 2 , la velocità uk+1
e
2
k+1
la pressione p2 sono determinate univocamente. Successivamente si risolve il problema di
Dirichlet (4.120) su Ω1 , ma, poiché la pressione è definita a meno di una costante additiva, si
deve aggiungere la condizione
Z
Z
k+1
p1 dΩ1 = −
pk+1
dΩ2 .
2
Ω1
Ω2
Ora è evidente che quando le successioni {uk1 }, {uk2 }, pk1 e pk2 convergono, la condizione sulla
media nulla per la pressione è automaticamente verificata.
Consideriamo ora il metodo di Schwarz per il problema di Stokes. In tal caso si impone soltanto
la continuità della velocità u:

k+1

S(uk+1
1 , p1 ) = f in Ω1 ,


uk+1
= uk2
su Γ1 ,
(4.121)
1


 uk+1 = 0
su ∂Ω \ Γ ,
1
1
1

k+1

S(uk+1
2 , p2 ) = f in Ω2 ,


uk+1
= uk+1
su Γ2 ,
2
1


 uk+1 = 0
su ∂Ω2 \ Γ2 .
2
(4.122)
Come si vede non vi è alcuna condizione sulla pressione ora.
La velocità del fluido deve avere divergenza nulla in Ω. Ora, quando si risolve (4.121), si ha
∇ · u1 = 0 in Ω1 e, grazie alla formula di Green, ne segue che
Z
u1 · n = 0.
∂Ω1
Dunque uk2 in (4.121)2 non può assumere un valore prefissato, ma deve soddisfare piuttosto la
condizione seguente:
Z
Z
Z
Z
Z
0=
u1 · n dγ =
u1 · n dγ +
u1 · n dγ =
u1 · n dγ =
uk2 · n dγ. (4.123)
∂Ω1
∂Ω1 \Γ1
Γ1
Γ1
Γ1
Alla prima iterazione si può scegliere uk2 in modo tale che la (4.123) sia soddisfatta, ma successivamente si perde del tutto il controllo su tale quantità. Inoltre, per lo stesso motivo, quando si
risolve (4.122), si deve avere che uk+1
soddisfa la condizione
1
Z
uk+1
· n = 0.
(4.124)
1
Γ
4.8. ALCUNE ESTENSIONI
151
Bisogna dunque garantire a priori che tale condizione sia sempre soddisfatta. Fortunatamente, il
metodo di Schwarz garantisce automaticamente tale condizione. Sia infatti Γ 12 = Ω1 ∩ Ω2 ; ora
in Γ12 si ha ∇ · u2 = 0, ma su Γ12 \ (Γ1 ∪ Γ2 ) uk2 è nulla in virtù delle condizioni al bordo. Allora
ne segue che
Z
Z
Z
0=
∂Γ12
uk+1
· n dγ =
1
Γ1
· n dγ +
uk+1
1
Γ2
uk+1
· n dγ,
1
essendo nullo il primo integrale a causa della (4.123), e dunque la (4.124) è verificata.
152
Operatore
CAPITOLO 4. DECOMPOSIZIONE DI DOMINI
Tabella 4.7: Condizioni di accoppiamento su Γ
Problema
Continuità (D)
Laplace
−4u = f
Elasticità
−∇ · (σ(u)) = f
con
σkj = µ̂(Dk uj + Dj uk ) + λ̂divuδkj ,
u spostamento della membrana nel piano
P
− kj Dk (Akj Dj u) + div(bu) + a0 u = f
Diffusione-trasporto
Trasporto
Stokes (fluidi viscosi
incomprimibili)
Stokes (fluidi viscosi
comprimibili)
Stokes (fluidi non viscosi
comprimibili)
Maxwell (regime armonico)
u
−α2 εE + iασE = f
σ(u) · n
u
u
X
∂u
akj Dj u · nk
=
∂nL
k
b ·nu
div(bu) + a0 u = f
−divT(u, p) + (u∗ · ∇)u = f
divu = 0
con
Tkj =
ν(Dk uj + Dj uk ) − pδkj ,
(Stokes)
 0
u∞ (Oseen)
u∗ =

u
(Navier-Stokes)
αu − divT̂(u, σ) = f
ασ + divu = g
con
T̂kj = ν(Dk uj +Dj uk ),
−βσδkj + g − 2ν
divuδkj ,
d
ρ = densità del fluido = log σ
αu + β∇σ = f
ασ +
divu = 0
1
rotE
rot
µ
Continuità (N)
∂u
∂n
u
T(u, p) · n
u
T̂(u, σ) · n
u·n
σ
n×E
n×
1
rotE
µ