Approssimazione di equazioni iperboliche in una dimensione spaziale

Capitolo 4
Problemi iperbolici
L’approssimazione numerica delle equazioni iperboliche
1
è un’area di ricerca molto attiva. La
caratteristica principale di questo tipo di problemi è quella che le perturbazioni si propagano
con velocità finita. Un’altra caratteristica è che, per quanto riguarda il bordo, in accordo
con il segno dei coefficienti dell’equazione, le regioni del bordo di "inflow" e di "outflow"
determinano,dacasoacaso,doveoccorredisporrelecondizionialbordo. Senonèimplementata
convenientemente, la realizzazione numerica delle condizioni al bordo è una potenziale sorgente
di instabilità.
4.1 Problemi con dati iniziali
Negli esempi che seguiranno
2
ci occuperemo di un tipo molto significativo di equazioni di
tipo iperbolico: le equazioni lineari di convezione, dette anche equazioni lineari di trasporto.
Affronteremo dunque questi esempi in modo qualitativo, tralasciando di specificare l’esatta
appartenenza delle funzioni incontrate nel prosieguo ad opportuni spazi funzionali. L’esempio
più semplice di problema iperbolico è dato quindi dalla seguente equazione:
@u
@t
(t;x) +a
@u
@x
(t;x) = 0; t> 0;x2R (4.1)
dove il coefficiente a2Rnf0g.
LasoluzionedelproblemadiCauchydefinitodall’equazione(4.1), aggiungendolacondizione
iniziale:
u(x; 0) =u
0
(x); x2R; (4.2)
è semplicemente l’onda viaggiante con velocità a:
u(t;x) =u
0
(x at); t 0:
La soluzione è costante lungo le curve caratteristiche, ovvero le curve x(t) del piano (t;x)
soddisfano le equazioni differenziali ordinarie:
 x
0
(t) =a;
x(0) =x
0
:
t> 0
Nel caso più generale, l’equazione lineare di convezione è data aggiungendo altri elementi
all’equazione (4.1):
1
Per una teoria più approfondita sulle equazioni iperboliche si veda G.B. WHITHAM, Linear and Nonlinear
Waves, John Wiley and Sons, New York (1974).
2
Si farà riferimento a J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations,
Wadsworth and Brooks/Cole, Pacific Grove (1989).
25
26 CAPITOLO 4. PROBLEMI IPERBOLICI
@u
@t
+a
@u
@x
+a
0
u =f; t> 0;x2R (4.3)
dove a;a
0
;f sono funzioni di (t;x) opportunamente date. Definendo ancora le curve carat-
teristiche come la soluzione del sistema di equazioni:
 x
0
(t) =a(t;x(t)) t> 0;
x(0) =x
0
otteniamo che la soluzione dell’equazione (4.1) soddisfa l’equazione differenziale ordinaria:
d
dt
u(t;x(t)) =f(t;x(t)) a
0
(t;x(t))u(t;x(t))
Dunque la soluzione lungo le caratteristiche non è più costante come nel caso precedente,
comunque può essere ottenuta risolvendo due insiemi di equazioni differenziali ordinarie.
Esaminiamo ora un lemma, che si rivelerà molto utile nel seguito di questo capitolo e anche
nei seguenti.
Lemma 4.1 (di Gronwall). Sia f 2 (t
0
;T ) una funzione non negativa, e siano g e ’ due
funzioni continue sull’intervallo chiuso [t
0
;T ]. Se ’ soddisfa:
(i) ’(t) g(t) +
R
t
t
0
f( )’( )d; 8t2 [t
0
;T ],
allora vale la seguente diseguaglianza:
(ii) ’(t) g(t) +
R
t
t
0
f(s)g(s) exp
 R
t
s
f( )d
 ds; 8t2 [t
0
;T ].
Inoltre, se g è anche una funzione non-decrescente, allora abbiamo che vale:
(iii) ’(t) g(t) exp
 R
t
t
0
f( )d
 8t2 [t
0
;T ]
Dimostrazione. Definiamo R(t) =
R
t
t
0
f( )’( )d; 8t2 [t
0
;T ]. Dalla condizione (i) segue che:
d
dt
R(t) =f(t)’(t) f(t)[g(t) +R(t)]
dunque
d
dt
 R(t) exp
  Z
t
t
0
f( )d
 =
 d
dt
R(t) R(t)f(t)
 exp
  Z
t
t
0
f( )d
   f(t)g(t) exp
  Z
t
t
0
f( )d
 ;8t2 [t
0
;T ]
Integrando l’ultima diseguaglianza nell’intervallo , troviamo che:
R(t) R(t) exp
  Z
t
t
0
f( )d
  Z
t
t
0
f(s)g(s) exp
  Z
s
t
0
f( )d
 ds; 8t2 [t
0
;T ] ;
da cui, ricorrendo nuovamente alla condizione (i), segue la validità della diseguaglianza (ii).
Se g è una funzione non-decrescente, allora da (ii) troviamo la seguente diseguaglianza:
’(t) g(t)
 1 +
Z
t
t
0
f(s) exp
 Z
t
s
f( )d
 ds
 ; 8t2 [t
0
;T ];
4.2. PROBLEMI AL CONTORNO 27
ma
d
ds
exp
 Z
t
s
f( )d
 = f(s) exp
 Z
t
s
f( )d
 quindi
Z
t
t
0
f(s) exp
 Z
t
s
f( )d
 ds = exp
 Z
t
s
f( )d
     t
s=t
0
= 1 + exp
 Z
t
t
0
f( )d
 da cui segue la validità della diseguaglianza (iii):
’(t) g(t) exp
 Z
t
t
0
f( )d
 ; 8t2 [t
0
;T ]:
4.2 Problemi al contorno
Vediamo ora di introdurre i problemi con condizioni che vengono assegnate sia al tempo iniziale
che sul bordo del dominio della funzione incognita. A questo scopo esaminiamo prima un
esempio e poi affrontiamo il problema più generale.
4.2.1 Equazione di convezione
Consideriamooraproblemiincuivengonoassegnatiivalorialtempoinizialeequellisulbordo, e
in particolare prendiamo in considerazione quelle equazioni iperboliche definite su un intervallo
spazialmente limitato che supponiamo sia I = (0; 1). Abbiamo ora bisogno di descrivere i dati
al bordo in modo tale che essi siano compatibili con la natura del problema esaminato. Per
l’equazione di convezione (4.1) le condizioni al bordo prendono la forma:
 u(t; 0) =’
 (t); se a> 0
u(t; 1) =’
+
(t); se a< 0
(4.4)
dove ’
 (t) sono funzioni date. Se a > 0, il punto x = 0 sul bordof0; 1g é il punto di inflow,
mentre x = 1 é il punto di outflow. La regola é ribaltata nel caso di a < 0. Similmente, per
l’equazione di convezione più generale (4.3), la soluzione u ha bisogno di essere descritta in
ogni punto del bordo inflow. In accordo con il segno di a(t;x), il numero di punti di inflow può
essere uno, due o nessuno. La soluzione del corrispondente problema con condizioni iniziali e al
bordo, ottenuto combinando l’equazione (4.3) con la condizione iniziale (4.2) e con le condizioni
al bordo (4.4), soddisfa una limitazione a-priori nella norma di L
2
(I). Invero, moltiplicando
(4.3) per u, integrando in x sull’intervallo I = (0; 1), e usando l’integrazione per parti sul
termine di convezione otteniamo per tutti i t> 0:
Z
1
0
@u
@t
(t;x)u(t;x)dx +
Z
1
0
a(t;x)
@u
@x
(t;x)u(t;x)dx +
Z
1
0
a
0
(t;x)u
2
(t;x)dx =
=
Z
1
0
f(t;x)u(t;x)dx
che è equivalente a :
1
2
Z
1
0
@u
2
@t
(t;x)dx +
1
2
Z
1
0
a(t;x)
@u
2
@x
(t;x)dx +
Z
1
0
a
0
(t;x)u
2
(t;x)dx =
Z
1
0
f(t;x)u(t;x)dx
28 CAPITOLO 4. PROBLEMI IPERBOLICI
oppure a :
1
2
d
dt
Z
1
0
u
2
(t;x)dx +
1
2
 a(t;x)u
2
(t;x)
  1
x=0
 Z
1
0
@a
@x
(t;x)u
2
(t;x)dx
 +
Z
1
0
a
0
(t;x)u
2
(t;x)dx =
=
Z
1
0
f(t;x)u(t;x)dx;
semplificando otteniamo:
1
2
d
dt
ku(t;x)k
2
0;I
+
Z
1
0
 a
0
(t;x) 1
2
@a
@x
(t;x)
 u
2
(t;x)dx =
=
Z
1
0
f(t;x)u(t;x)dx +
1
2
  a(t; 1)u
2
(t; 1) +a(t; 0)u
2
(t; 0)
 Poniamo ora  (t;x) :=a
0
(t;x) 1
2
@a
@x
(t;x), e definiamo la funzione ( t) come:
( t) :=  +
(t)a(t; 1)’
2
+
(t) +  (t)a(t; 0)’
2
 (t); dove:
 +
(t) :=
(
1 se a(t; 1)< 0
0 altrimenti,
  (t) :=
(
1 se a(t; 0)> 0
0 altrimenti
da cui segue subito che ( t) 0, ed inoltre:
 a(t; 1)u
2
(t; 1) +a(t; 0)u
2
(t; 0) ( t); 8t> 0
abbiamo che l’equazione diventa:
1
2
d
dt
ku(t;x)k
2
0;I
+
Z
1
0
 (t;x)u
2
(t;x)dx Z
1
0
f(t;x)u(t;x)dx +
1
2
( t)
integrando ora in s su (0;t) otteniamo:
1
2
ku(t;x)k
2
0;I
+
Z
t
0
  (s;x);u
2
(s;x)
 0;I
ds Z
t
0
(f(s;x);u(s;x))
0;I
ds+
1
2
Z
t
0
( s)ds+
1
2
ku(0;x)k
2
0;I
;
che è equivalente a :
ku(t;x)k
2
0;I
 2
Z
t
0
(f(s;x);u(s;x))
0;I
ds+
Z
t
0
( s)ds+ku(0;x)k
2
0;I
 2
Z
t
0
  (s;x);u
2
(s;x)
 0;I
ds;
applicando ora il valore assoluto ad entrambi i membri della disequazione e utilizzando la
diseguaglianza di Hölder otteniamo:
ku(t;x)k
2
0;I
 2
Z
t
0
kf(s;x)k
0;I
ku(s;x)k
0;I
ds +
Z
t
0
( s)ds +ku(0;x)k
2
0;I
+
+
    2
Z
t
0
 j (s;x)j;u
2
(s;x)
 0;I
ds
    :
Sia ora :   (t) := max
x2I
(j (t;x)j): allora, maggiorando la diseguaglianza precedente e tenendo
conto del dato iniziale, otteniamo :
ku(t;x)k
2
0;I
 Z
t
0
2kf(s;x)k
0;I
ku(s;x)k
0;I
ds +
Z
t
0
( s)ds +ku
0
k
2
0;I
+ 2
Z
t
0
  (s)ku(s;x)k
2
0;I
ds;
4.2. PROBLEMI AL CONTORNO 29
usando la diseguaglianza : 2ab a
2
+b
2
;8a;b2R, si ha :
ku(t;x)k
2
0;I
 Z
t
0
 kf(s;x)k
2
0;I
+ku(s;x)k
2
0;I
 ds+
Z
t
0
( s)ds+ku
0
k
2
0;I
+2
Z
t
0
  (s)ku(s;x)k
2
0;I
ds;
raccogliendo ora il termineku(s;x)k
2
0;l
, abbiamo :
ku(t;x)k
2
0;I
k u
0
k
2
0;I
+
Z
t
0
kf(s;x)k
2
0;I
ds +
Z
t
0
( s)ds +
Z
t
0
[1 + 2  (s)]ku(s;x)k
2
0;I
ds;
definiamo ora:
’(t) :=ku(t;x)k
2
0;I
g(t) :=ku
0
k
2
0;I
+
Z
t
0
kf(s;x)k
2
0;I
ds +
Z
t
0
( s)ds; non decrescente;
f(s) := 1 + 2  (s) 0:
’(t);g(t);f(s) verificano tutte le condizioni del lemma di Gronwall, dunque:
’(t) g(t) +
Z
t
0
f(s)’(s)ds ) ’(t) g(t) exp
 Z
t
0
f(s)ds
 :
Il risultato finale è dato da:
ku(t;x)k
2
0;I
  ku
0
k
2
0;I
+
Z
t
0
kf(s;x)k
2
0;I
ds +
Z
t
0
( s)ds
 exp
 Z
t
0
[1 + 2  (s)]ds
 ; 8t2 (0; 1):
4.2.2 Equazioni non lineari
Consideriamo l’equazione scalare:
@u
@t
+
@
@x
F (u) = 0; t> 0;x2R; (4.5)
dove F (u) è una funzione non-lineare di u, e u(0;x) =u
0
(x):
3
Molto spesso, F è una funzione strettamente convessa, cioè F
00
( )> 0;8 2R. Un classico
esempio di questa classe di problemi ci viene fornito dall’equazione di Burgers che ricaviamo
dall’equazione (4.5) assumendo F (u) =
u
2
2
, cioè:
@u
@t
+u
@u
@x
= 0; t> 0;x2R (4.6)
Questa è un’equazione di trasporto non-lineare la cui velocità di convezione è data dal-
la stessa soluzione u. Le curve caratteristiche soddisfano la seguente equazione differenziale
ordinaria:
x
0
(t) =u(t;x(t)); (4.7)
che utilizziamo per calcolare la derivata totale di u lungo ciascuna caratteristica, cioè:
d
dt
u[u(t;x(t))] =
@u
@t
(t;x(t)) +x
0
(t)
@u
@x
(t;x(t)) =
 @u
@t
+u
@u
@x
 (t;x(t)) = 0
3
Una teoria più completa sulle equazioni non lineari può essere trovata in A.QUARTERONI, A.VALLI,
Numerical approximation of partial differential equations, Springer-Verlag, Berlin (1994).
30 CAPITOLO 4. PROBLEMI IPERBOLICI
Quindi, u è costante su ciascuna caratteristica. Ovvero, a causa di (4.7) la pendenza x
0
(t) è
costante, e dunque le caratteristiche sono linee rette la cui pendenza è determinata dal dato
iniziale. Quando quest’ultimo è regolare, questa proprietà comporta che la soluzione u(t;x)
assume la forma:
u(t;x) =u
0
(s)
dove s =s(t;x) soddisfa:
x(t) =s +u
0
(s) t (4.8)
Quest’ultima equazione ha un’unica soluzione per t abbastanza piccoli (in modo che le
caratteristiche non si incrocino. Ogni volta che u
0
0
è negativo in ciascun punto, detto s
0
,
esiste un tempo critico t
c
in cui le prime caratteristiche si incrociano, la funzione u (t
c
;x) ha
una pendenza infinita nel punto x
0
e una forma discontinua. Per l’equazione di Burgers x
0
è
ottenuto da (4.8) prendendo s =s
0
: Si può vedere che:
t
c
= min
 1
u
0
0
(x)
(4.9)
Infatti se u
0
(s
1
)>u
0
(s
2
) per s
1
<s
2
, allora le due caratteristiche:
x
1
(t) =s
1
+u
0
(s
1
) t
x
2
(t) =s
2
+u
0
(s
2
) t
si incontreranno ad un certo tempo
^
t. In tal caso abbiamo che:
^
t =
s
1
 s
2
u
0
(s
2
) u
0
(s
1
)
ma:
lim
s
2
!s
1
s
1
 s
2
u
0
(s
2
) u
0
(s
1
)
=
 1
u
0
0
(s
1
)
Il tempo critico è dunque dato dall’inf dei tempi di incrocio.
4.2.3 Il problema astratto e la sua riformulazione variazionale
Il problema più generale con valori iniziali e al bordo viene dedotto dal problema con valori al
bordo (dunque stazionario) astratto descritto nel secondo capitolo e si presenta nel seguente
modo:
8
<
:
@u
@t
+ Lu =f in Q
T
:= (0;T )  Bu = 0 su   T
:= (0;T ) @  u =u
0
su  ; per t = 0
(4.10)
dove T > 0 è un livello di tempo fissato,
@
@t
denota la differenziazione rispetto al tempo,
u (la funzione incognita) e f sono funzioni di t2 (0;T );x2  e u
0
= u
0
(x) è il dato iniziale
assegnato. L’operatore differenziale L e l’operatore al bordo B possono ora dipendere anche da
t.
Supponiamo ora che esistano tre spazi di Hilbert V;W;H tali che V e W sono contenuti in
H con inclusione continua e densa (solitamente H = L
2
() ). Denotiamo il prodotto scalare di
H con (; )
H
. Assumiamo inoltre che u
0
2 H e f2 L
2
(0;T ; H): Supponiamo poi che esista una
forma bilineare A(; ) continua su W V: La formulazione debole del problema (4.10) si legge
nel seguente modo:
si trovi u2L
2
(0;T ;W )\C
0
([0;T ];H) tale che:
d
dt
(u(t);v)
H
+ A(u(t);v) = F(t;v); 8v2V (4.11)
4.3. DISCRETIZZAZIONE DEL PROBLEMA 31
e u =u
0
a t = 0:
Il membro destro dell’equazione è una notazione compatta per (f(t);v)
H
più altri possibili
termini dipendenti da condizioni al bordo non omogenee. Le equazioni (4.11) sono da intendersi
nel senso delle distribuzioni su (0;T ): Sotto le precedenti assunzioni sui dati tutte le operazioni
indicate in (4.11) hanno senso.
4.3 Discretizzazione del problema
Occorre discretizzare il problema (4.11) rispetto ad entrambe le variabili tempo e spazio.
4.3.1 Approssimazione semi-discreta
Per quanto riguarda la discretizzazione spaziale, essa può essere realizzata allo stesso modo dei
problemi stazionari. Per esempio l’approssimazione di Galerkin viene ottenuta generalizzando
l’approccio (2.15). Il problema risultante si legge nella seguente maniera:
8t2 (0;T ) si trovi u
h
(t; )2V
h
tale che:
d
dt
(u
h
(t);v
h
)
H
+ A (u
h
(t);v
h
) = F (t;v
h
); 8v
h
2V
h
;t2 (0;T ) (4.12)
con u
h
(0) = u
0;h
, il quale è un elemento opportuno di V
h
che approssima il valore iniziale
u
0
:
Diamo ora le seguenti definizioni di convergenza, di stabilità, e di consistenza, dedotte dal
secondo capitolo, con u(t) soluzione esatta del problema (4.11).
Diremo che l’approssimazione (4.12) è:
- convergente se:
ku(t) u
h
(t)k! 0; quando h! 0; per t2 (0;T );
- stabile se:
ku
h
(t)k C
t
; per ogni h> 0; per t2 (0;T ); dove C
t
è una costante opportuna;
- consistente se:
sup
v
h
2V
h
v
h
6=0
  d
dt
(u(t);v
h
)
H
+ A (u(t);v
h
) F (t;v
h
)
  kv
h
k
! 0; quando h! 0; per t2 (0;T );
- fortemente consistente se:
d
dt
(u(t);v
h
)
H
+ A (u(t);v
h
) F (t;v
h
) = 0; 8v
h
2V
h
;t2 (0;T ):
Il problema (4.12) è spesso chiamato una "approssimazione semi-discreta" (od anche "con-
tinua nel tempo") del problema (4.11). Ponendo:
u
h
(t;x) =
N
h
X
j=1
 j
(t)’
j
(x); pert 0; (4.13)
u
0;h
(x) =
N
h
X
j=1
 0;j
’
j
(x);
32 CAPITOLO 4. PROBLEMI IPERBOLICI
dovef’
j
g ancora denota una base di V
h
0 si vede facilmente che (4:12) produce un sistema di
equazioni differenziali ordinarie per la funzione incognita  (t) := ( j
(t)) che si legge così:
M
d (t)
dt
+ A (t) =F(t); 8t 0; (4.14)
 (0) = 0
;
dove:
M
ij
:= (’
j
;’
i
) (4.15)
A
ij
:= A (’
j
;’
i
)
F(t) := (F (t;’
i
))
coni;j = 1;:::;N
h
.
La matrice di stiffness A è la stessa utilizzata nel caso stazionario. E’ indipendente da t
se anche l’operatore L ne è indipendente. La matrice M, simmetrica e definita positiva, che è
sempre indipendente da t, è chiamata invece matrice di massa.
4.3.2 Approssimazione totalmente-discreta
La discretizzazione nel tempo di (4.12) può essere compiuta in diversi modi.
Innanzitutto suddividiamo l’intervallo di tempo [0;T ] in N sotto-intervalli [t
n
;t
n+l
] di lun-
ghezza  t = T=N, con t
0
= 0 e t
N
= T. Quindi cerchiamo una approssimazione di u(t) a
ciascun livello di tempo t
n
, e denotiamo con u
n
h
la corrispondente funzione dimensionalmente
finita. In particolare, u
0
h
è una conveniente approssimazione di u
0
. Un modo classico e molto
diffuso di realizzare una "approssimazione totalmente discreta" del problema (4.12) è quello di
ricorrere ad una discretizzazione delle derivate temporali con uno schema a differenze finite.
Nel prossimo capitolo vedremo di analizzare la stabilità e la convergenza del metodo di
approssimazione di Galerkin, riferito all’equazione lineare di convezione espresso sia in una
forma forte che in una forma debole. Andremo poi a proporre altri metodi di approssimazione,
varianti di quello di Galerkin, che risulteranno migliori e di cui analizzeremo le proprietà di
stabilità e di convergenza.
Nel capitolo successivo affronteremo invece tutto un altro tipo di approccio del problema di
approssimazione; quello dei metodi a differenze finite dei quali vedremo in maniera analoga di
analizzare le proprietà di stabilità e di convergenza.