Post on 18-Jan-2019
POLITECNICO DI MILANO
Facoltà di Ingegneria dei SistemiCorso di Laurea Specialistica in Ingegneria Matematica
SPECTRAL COLLOCATION METHODS ONSEMI–INFINITE DOMAINS AND APPLICATION
TO OPEN BOUNDARY CONDITIONS
Relatore CandidatoDott. Luca Bonaventura Tommaso Benacchio
Matr. 711145
Anno Accademico 2009–2010
Se chiudete la porta a tutti gli errori,anche la verità ne resterà fuori.
R. Tagore
Ai miei genitori e a mia nonna
Abstract. We introduce a spectral approach to the problem of numerical open boundary conditions with
applications in environmental modelling. When differential problems are defined on unbounded domains,
an artificial boundary is introduced to encompass a finite region of interest, and the equations are solved
in this region by standard techniques. Suitable conditions have to be enforced on the boundary that
avoid the occurrence of spurious or unphysical phenomena; however, conventional rigid lids and constant
pressure surfaces entail energy reflection, while the classical approaches of radiative boundary conditions
and absorbing/sponge layers are either difficult to derive or computationally expensive and inefficient.
We tackle the problem in a more direct way, approximating with spectral methods the equations in a
semi–infinite domain attached to the finite region of interest. We review the main results on polynomial
approximation and interpolation on the real positive half line with Laguerre polynomials and Laguerre
functions, notably taking into account a special class of scaled Laguerre functions which enable to represent
different spatial scales in the same framework. Next, we apply these results to the approximation of open
channel flow equations with spectral methods on semi–infinite domains. The linear and nonlinear shallow
water equations are discretized on a half line with collocation methods in space and semi–Lagrangian,
semi–implicit time discretization techniques. The resulting stable, spectrally accurate method yields a
tool which can be interfaced with a numerical discretization of the equations in a finite domain and, using
low–order quadrature rules, is expected to outperform the cost of traditional absorbing layers. To give an
example, we consider a finite volume discretization in a finite domain and attach to the right endpoint
a semi–infinite discretization. We show that a reasonable number of spectral base functions is sufficient
to reach the accuracy necessary for practical implementation of open boundary conditions, along with an
accurate treatment of the solution behaviour at infinity.
Sommario. Il lavoro si propone di studiare un approccio con metodi spettrali per il trattamento di
condizioni al contorno aperte applicate alla modellistica ambientale. Nella risoluzione di problemi diffe-
renziali su domini illimitati viene introdotta una frontiera artificiale intorno alla regione finita di interesse
in cui le equazioni del modello in esame vengono risolte con metodi di discretizzazione standard. Su tale
frontiera vanno imposte opportune condizioni al contorno tali da evitare l’insorgere di fenomeni spuri,
tuttavia scelte tradizionali come pareti rigide o superfici a pressione costante causano riflessioni di ener-
gia, mentre approcci classici quali condizioni di tipo radiazione o strati assorbenti/spugne sono difficili da
ricavare o computazionalmente inefficienti. Nel presente lavoro il problema viene affrontato direttamente,
approssimando le equazioni con un metodo spettrale in un dominio semi–infinito collegato alla regione
finita. Si studiano i principali risultati di approssimazione e interpolazione polinomiale sulla semiretta
reale positiva con polinomi e funzioni di Laguerre, considerando anche una particolare classe di funzioni di
Laguerre riscalate che permettono di rappresentare diverse scale spaziali. Si applicano poi questi risultati
per l’approssimazione delle equazioni per fluidi in canale aperto su domini semi–infiniti. Si discretizzano le
equazioni delle acque basse lineari e non lineari su una semiretta con un metodo di collocazione spettrale
in spazio e tecniche di tipo semi–implicito semi–Lagrangiano in tempo. Il metodo che ne risulta è stabile
e accurato e fornisce uno strumento che, collegato a una discretizzazione finita, grazie all’uso di formule
di quadratura di basso ordine, ha un minor costo computazionale rispetto agli strati assorbenti standard.
Come esempio si considera una approssimazione ai volumi finiti in un dominio finito e si collega nell’e-
stremo destro la discretizzazione semi–infinita. Si mostra che è sufficiente un numero relativamente basso
di funzioni di base spettrali per raggiungere l’accuratezza necessaria all’implementazione di condizioni al
contorno aperte, con un trattamento accurato della soluzione all’infinito.
Preface
The aim of the work is to study numerical methods based on Laguerre functions expansions
for the approximation of wave propagation problems on semi–infinite domains. The purpose
of this development is twofold.
On one hand, in many environmental applications, such as for example climate or ocean mod-
elling, it is of interest to consider computational domains that span over very large length scales.
In climate modelling, increasing attention has been devoted to stratospheric phenomena (see e.g.
[HPS+03, GMR+06]), which requires to extend the typical computational domain of standard
climate model by a substantial amount. It is therefore of interest to assess whether spectral or
pseudo–spectral approaches can be used to derive accurate and efficient discretizations for do-
mains with arbitrarily large length scales. This goal will be pursued both by direct application
of Laguerre spectral methods and by their coupling to standard finite volume discretizations on
domains of finite size by a simple domain decomposition approach. As a result, a semi–implicit and
semi–Lagrangian discretization has been developed for the one dimensional shallow water model
equations on semi–infinite domains, along with a coupling approach based on imposing continuity
of the mass fluxes at the interface between finite size and unbounded domain, that can be in
principle extended to couple the Laguerre spectral method to an arbitrary discretization on the
finite size domain.
The semi–implicit, semi–Lagrangian approach (see e.g. [Cas90, SC91, RCB05]) has been chosen
because of its well known efficiency and accuracy features and because a large number of envi-
ronmental models use these time discretization techniques. In this first attempt, an explicit finite
volume discretization has been considered for the finite size domain, both for the sake of simplicity
and in order to show that such a coupling is feasible also if two essentially different discretiza-
tion approaches are employed on either domain. Numerical tests show that spurious reflections
due to the proposed coupling approach have very small amplitude and allow to use the proposed
method, at least for environmental applications, with no substantial loss of accuracy with respect
to discretizations using a single computational domain.
viii Preface
On the other hand, for computational reasons, it is common to restrict the simulation to a
bounded domain of interest and to apply some kind of open boundary condition in order to allow
for wave propagation out of the computational domain. While conventional boundary conditions
such as rigid lids or constant pressure surfaces entail total reflection of wave energy, two classi-
cal approaches to this problem, radiative boundary conditions and absorbing/sponge layers (see
e.g. [Giv08]), have been extensively applied in the last three decades to atmospheric models for
numerical weather prediction and climate. A review of most commonly used approaches can be
found in [Ras86]. The other main goal of this work is to use Laguerre spectral methods in order to
implement absorbing layer open boundary conditions. By coupling discretizations on domains of
finite size to Laguerre spectral methods employing a relatively small number of nodes, an efficient
approach to absorbing boundary conditions is achieved, that allows to account for large absorbing
regions at a relatively low computational cost.
In chapter 1 we first review the main results concerning polynomial approximation and inter-
polation on the real positive half line by means of Laguerre polynomials and Laguerre functions,
with specific reference to quadrature formulas and spectral derivatives [She01, SW09]. In partic-
ular, we consider a special class of generalized Laguerre functions dependent on a parameter β,
tuning which one is able to represent functions on different spatial scales in a unified framework
[WGW09].
In chapter 2 the shallow water equations are introduced and their derivation from the Navier–
Stokes equations for free–surface fluids is briefly discussed. Furthermore, different numerical ap-
proaches to the approximation of open boundary conditions are reviewed and discussed, with
special attention devoted to techniques employed in environmental modelling and, in particular,
to absorbing boundary conditions.
In chapter 3 spectral and pseudo–spectral numerical methods based on the scaled Laguerre
functions are introduced, in order to show that the unsteady advection–diffusion–reaction equation
and the shallow water equations can be approximated accurately and efficiently on semi–infinite
domains by extensions of semi–implicit, semi–Lagrangian methods that are customary in many
environmental applications.
In chapter 4 the coupling strategy is described, that was proposed to connect finite volume
discretizations of the shallow water equations on a finite size domain to the previously introduced
pseudo–spectral, semi–implicit semi–Lagrangian method.
In chapter 5 we present the results of a number of numerical tests. A first goal of these tests
is to validate the pseudo–spectral, semi–implicit and semi–Lagrangian discretization approach
on semi–infinite domains as an independent tool. The coupling approach is then validated, by
comparison of the coupled finite volume/pseudo–spectral model results to those of reference runs
of the finite volume model on a large domain. Finally, absorbing boundary conditions have been
Preface ix
validated, as implemented in the context of the coupled finite volume/pseudo–spectral model.
Numerical simulations show that a reasonable number of spectral base functions is sufficient to
reach the accuracy necessary for practical implementations of open boundary conditions, along
with an accurate treatment of the solution behaviour at infinity.
Some conclusions on the present work are drawn in a final chapter, where the perspectives for
further development of this approach are also discussed.
Contents
1 Polynomial approximation on semi–infinite domains . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Orthogonal polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Approximation with orthogonal polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Gaussian quadrature formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Zeros of orthogonal polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.5 Discrete polynomial transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.6 Spectral differentiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.7 The case of semi–infinite domains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.8 Rescaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2 Free–surface flow in shallow water assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1 Geometrical setup and notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Model derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3 Advective formulation and linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.4 The problem of open boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3 Spectral discretization of boundary–value problems on semi–infinite domains 39
3.1 The spectral G–NI and collocation approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Advection–diffusion–reaction equation: G–NI approximation . . . . . . . . . . . . . . . . . . . 41
3.3 Linear shallow water model: collocation approximation . . . . . . . . . . . . . . . . . . . . . . . 45
3.4 Nonlinear shallow water model: semi–Lagrangian approximation . . . . . . . . . . . . . . . 50
4 Finite–domain approximation and coupling with the semi–infinite part . . . . . . 55
4.1 Finite–volume approximation of the shallow water equations . . . . . . . . . . . . . . . . . . 55
4.2 Coupling the finite and semi–infinite discretizations . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.1 Scaling properties of scaled Laguerre functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
xii Contents
5.2 Unsteady advection–diffusion–reaction equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.3 Validation of the spectral collocation method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.4 Validation of the coupling approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.5 Validation of the efficiency for absorbing boundary conditions . . . . . . . . . . . . . . . . . 84
Conclusions and future perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
A Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
1
Polynomial approximation on semi–infinite domains
The aim of this first chapter is to present a review of the theory of approximation with
orthogonal polynomials and Gaussian numerical quadrature formulas. We present first clas-
sical results in a general setting. Our focus will be then directed to approximation and numerical
quadrature in semi–infinite domains. For the first part we refer to [QSS07, ST06, San04], for the
second to [She01, GS00, CHQZ06].
1.1 Orthogonal polynomials
In this section we analyze the properties of orthogonal polynomials. Let I := (a, b) ⊂ R ∪
−∞,+∞ be an interval and ω(x) a weight function, that is ω(x) > 0 in (a, b), ω ∈ L1(a, b). Let
L2ω(I) be the space of functions u such that
‖u‖ω :=
(∫ b
a
|u(x)|2ω(x)dx
)1/2
is finite. L2ω(I) is a Hilbert space, equipped with the inner product
(u, v)ω :=∫ b
a
u(x)v(x)ω(x)dx.
Furthermore, we denote by PN (I) the space of polynomials of degree at most N in I, that is
PN (I) =
u | u(x) =
N∑n=0
cnxn, cn ∈ R
.
Definition 1 (Orthogonality). Two functions f and g are orthogonal in L2ω(I) if
(f, g)ω = 0.
A sequence of polynomials pnn∈N with deg(pn) = n:
pn(x) = cnxn + cn−1x
n−1 + . . .+ c1x+ c0, cn 6= 0 (leading coefficient)
2 1 Polynomial approximation on semi–infinite domains
is said to be orthogonal in L2ω(I) if
(pn, pm)ω =∫ b
a
pn(x)pm(x)ω(x)dx = γnδmn,
where δmn is the Kronecker delta, and γn = ‖pn‖2ω 6= 0.
A first important result (whose proof can be found in [ST06]) is that orthogonal polynomials can
be constructed via a three–term recurrence relation:
Theorem 1. (Gram-Schmidt) Let pnn∈N be a sequence of orthogonal polynomials in L2ω(I) with
deg(pn) = n, and let cn be the leading coefficient of pn. Then we have:
pn+1 = (anx− bn) pn − knpn−1, n ≥ 0, (1.1)
with p−1 := 0, p0 := 1 and
an = cn+1
cn, bn = cn+1
cn
(xpn, pn)ω‖pn‖2
ω
, kn = cn+1cn−1
c2n
‖pn‖2ω
‖pn−1‖2ω
.
The converse is also true, that is: if a sequence of polynomials pnn∈N satisfies the three–term
recurrence relation (1.1), then pnn∈N is orthogonal in L2ω.
1.2 Approximation with orthogonal polynomials
By the Stone–Weierstrass theorem, any continuous function in a bounded interval is uniformly
approximated by polynomials. A similar result holds for unbounded intervals, too (see [ST06]):
Theorem 2. If u ∈ C ([0,∞)) and for certain δ > 0 u satisfies
u(x)e−δx → 0, as x→∞,
then, for any ε > 0, there exist n ∈ N, pn ∈ PN , x0 ∈ R such that, for any x > x0
|u(x)− pn(x)| e−δx ≤ ε.
In general, though, the derivation of the best uniform approximation polynomial (i.e. p∗N ∈ PNs.t. ‖u− p∗N‖L∞(I) = inf
pN∈PN‖u− pN‖L∞(I)) is not easy. However, we can easily construct the best
L2ω–approximation polynomial:
Theorem 3. Let I = (a, b) be a bounded or unbounded interval. Then, for u ∈ L2ω(I) and n ∈ N
there exists a unique q∗N ∈ PN such that
‖u− q∗N‖ω = infqN∈PN
‖u− qN‖ω, (1.2)
where
1.3 Gaussian quadrature formulas 3
q∗N (x) =N∑k=0
ukpk(x),
and
uk = 1‖pk‖2
ω
∫ b
a
u(x)pk(x)ω(x)dx, (1.3)
pnNn=0 being an L2ω–orthogonal basis of PN .
(1.3) is the polynomial transform of u. We denote the best approximation polynomial q∗N by πNu,
the L2ω–orthogonal projection of u on PN , defined by
(πNu, ϕ)ω = (u, ϕ)ω ∀ϕ ∈ PN ,
and, as seen above, we have:
πNu =N∑k=0
ukpk,
the N–th order truncated series of u, whereas the formal series of u ∈ L2ω in terms of pkk∈N is:
Su =∞∑k=0
ukpk, (1.4)
where uk is as in (1.3). The (1.2), seen in terms of πNu, reads:
‖u− πNu‖L2ω(I) = inf
qN∈PN‖u− qN‖L2
ω(I),
and is the expression of the projection theorem. Furthermore, we have that:
limN→∞
‖u− πNu‖ = 0,
and (Parseval’s identity):
‖u‖2ω =
∞∑k=0
(uk)2 ‖pk‖2ω.
1.3 Gaussian quadrature formulas
We review here Gaussian quadrature formulas. The basic idea is to choose the quadrature nodes in
order to reach the best numerical approximation of the integral of a smooth function. In general, a
numerical quadrature rule to approximate the integral of a function f over an interval (a, b) takes
the form: ∫ b
a
f(x)ω(x)dx =N∑j=0
f(xj)ωj + EN [f ],
where xjNj=0, ωjNj=0 are the quadrature nodes and weights, respectively, and EN [f ] is the
approximation error. If EN [f ] ≡ 0, the quadrature rule is exact for f .
Given N+1 distinct nodes xjNj=0 ⊆ [a, b], polynomial quadrature rules are exact by construction
for polynomials up to degree N . Indeed, if we set:
4 1 Polynomial approximation on semi–infinite domains
ωj =∫ b
a
hj(x)ω(x)dx, (1.5)
with:
hj(x) =
∏i 6=j
(x− xi)∏i 6=j
(xj − xi),
being the Lagrange interpolation polynomial associated to the nodes xjNj=0, we have:∫ b
a
p(x)ω(x)dx =N∑j=0
p(xj)ωj ∀p ∈ PN .
Picking as quadrature nodes the zeros of orthogonal polynomials in (a, b) enables to reach a higher
degree of exactness, as stated in this classical result (proof can be found, for instance, in [ST06]):
Theorem 4 (Gaussian quadrature). Let xjNj=0 be the set of zeros of the orthogonal polynomial
pN+1. Then there exists a unique set of quadrature weights ωjNj=0 such that:∫ b
a
p(x)ω(x)dx =N∑j=0
p(xj)ωj ∀p ∈ P2N+1. (1.6)
In particular, ωj > 0 ∀j = 0, . . . , N and
ωj = cN+1
cN
‖pN‖2ω
pN (xj)p′N+1(xj)0 ≤ j ≤ N,
cj being the leading coefficient of the polynomial pj.
Gaussian quadrature is optimal because it is exact for all polynomials up to degree 2N + 1, and
it is impossible to find a set of nodes xjNj=0 and weights ωjNj=0 such that (1.6) holds for all
p ∈ P2N+2 (see also [QSS07]).
In view of the use of quadrature formulas in the context of boundary–value problems on a semi–
infinite interval we analyze also the Gauss–Radau quadrature rules, whose nodes include the left
endpoint x = a. To include both endpoints, one should use the Gauss–Lobatto nodes and weights,
which we do not consider here (see for instance [CHQZ06, QV08]).
First, we define
qN (x) := pN+1(x) + αNpN (x)x− a
,
with
αN = −pN+1(a)pN (a)
.
We have qN ∈ PN , and for any rN−1 ∈ PN−1 orthogonality yields:∫ b
a
qN (x)rN−1(x)ω(x)(x− a)dx =∫ b
a
(pN+1(x) + αNpN (x))rN−1(x)ω(x)dx = 0.
Hence, qNN∈N is a sequence of orthogonal polynomials in L2ω(a, b), where ω(x) := ω(x)(x− a).
Note also that the leading coefficient of qN is cN+1. We can now state the following result (see
[ST06]):
1.4 Zeros of orthogonal polynomials 5
Theorem 5 (Gauss–Radau quadrature). Let x0 = a and xjNj=0 be the zeros of qN . Then there
exists a unique set of quadrature weights ωjNj=0, defined by (1.5) such that:∫ b
a
p(x)ω(x)dx =N∑j=0
p(xj)ωj ∀p ∈ P2N . (1.7)
Moreover, ωj > 0 ∀j = 0, . . . , N and:
ω0 = 1qN (a)
∫ b
a
qN (x)ω(x)dx, (1.8)
ωj = 1xj − a
cN+1
cN
‖qN−1‖2ω
qN−1(xj)q′N (xj)1 ≤ j ≤ N. (1.9)
1.4 Zeros of orthogonal polynomials
We review here two ways of computing the zeros of orthogonal polynomials. We follow the pre-
sentation in [ST06].
Eigenvalue method
Theorem 6. The zeros xjnj=0 of the orthogonal polynomial pn+1 are the eigenvalues of the
following symmetric tridiagonal matrix:
An+1 =
α0√β1
√β1 α1
√β2
. . . . . . . . .√βn−1 αn−1
√βn
√βn αn
,
where:
αj = bjaj, for j ≥ 0; βj = kj
aj−1aj, for j ≥ 1.
with ak, bk, kk being the coefficients of the three–term recurrence relation given in Theorem 1.
Efficient eigensolvers can be used to compute Gaussian quadrature weights (see e.g. [QSS07])
Iterative root–finding method
With this approach (for which we refer to [ST06]), given an initial guess for the zero location, we
want to improve it by an iterative procedure. If x0j is an initial approximation to the j–th zero of
pn+1, the iterative scheme takes the form:xk+1j = xkj +D(xkj ), k ≥ 0
given x0j , and a terminal rule
(1.10)
6 1 Polynomial approximation on semi–infinite domains
The deviation term D(·) classifies different types of iterative schemes (see [ST06]). For instance,
the Newton’s method is a second–order scheme with:
D(x) = −pn+1(x)p′n+1(x)
,
whereas the Laguerre method is a third–order scheme with:
D(x) = −pn+1(x)p′n+1(x)
−pn+1(x)p′′n+1(x)
2(p′n+1(x))2 .
In practice, few iterations are often sufficient to reach a good estimate of the zeros. As usual, one
should pay attention to the fact that, if the initial estimation of the zeros is not sufficiently close,
the algorithm may converge to unwanted zeros.
1.5 Discrete polynomial transforms
In this section we give the definition of interpolating polynomial and the representation of a
function in terms of its values at the quadrature nodes. Hereafter, xjNj=0 and ωjNj=0 will be
Gauss or Gauss–Radau quadrature nodes and weights, though all the statements also hold for
Gauss–Lobatto quadrature rule. The corresponding discrete inner product and norm are:
(u, v)N,ω :=N∑j=0
u (xj) v (xj)ωj , ‖u‖N,ω := (u, v)2N,ω .
In particular, for the Gauss–Radau quadrature we have, thanks to (1.7),
(u, v)N,ω = (u, v)ω ∀u, v s.t. uv ∈ P2N .
Definition 2. Let Λ := [a, b). The operator IN : C(Λ)→ PN such that:
INu(xj) = u(xj) ∀j = 0, . . . N (1.11)
holds for all u ∈ C(Λ) is the interpolation operator.
(1.11) implies that INp = p ∀p ∈ PN , as well as (INu, v)N,ω = (u, v)N,ω ∀v ∈ C(a, b); that is,
the interpolant INu is the orthogonal projection of u upon PN with respect to (·, ·)N,ω.
Moreover, since INu ∈ PN , we have:
(INu)(x) =N∑k=0
ukpk(x),
which is called the discrete polynomial series. It follows that:
u(xj) = (INu)(xj) =N∑k=0
ukpk(xj), (1.12)
1.6 Spectral differentiation 7
(the inverse discrete polynomial transform), while the uk are the discrete polynomial coefficients
and have the expression:
uk = 1γk
N∑j=0
u(xj)pk(xj)ωj , 0 ≤ k ≤ N, (1.13)
where:
γk = ‖pk‖2ω, 0 ≤ k ≤ N.
Thanks to (1.12) and (1.13), it is easy to switch from a representation in physical space, given by
u(xj), to a representation in frequency space, given by uk. Moreover, the coefficients ukNk=0
and ukNk=0 (see (1.4)) are related by:
uk = uk + 1γk
∑l>N
(pl, pk)N,ωul k = 0, . . . , N.
1.6 Spectral differentiation
The definition of INu given above enables also to approximate the derivatives of u. The derivatives
of the function are approximated by the analytical derivatives of the interpolating polynomial (see
[ST06] and [QV08]). This fact turns out to be very useful when we deal with spectral discretization
of boundary–value problems.
The pseudo–spectral derivative DN (u) of a continuous function u is defined by:
DN (u) = D [INu] ,
that is, DN (u) is the derivative of the interpolating polynomial of u. Moreover, DN can be ex-
pressed in terms of a matrix, the pseudo–spectral derivation matrix D:
D = [dij ]i,j=0,...,N ∈MN+1(R). (1.14)
Indeed, given the nodes xjNj=0, an approximation u ∈ PN of an unknown function and hjNj=0,
the Lagrange interpolation polynomials associated to the points xj , differentiating m times the
expression
u(x) =N∑j=0
u(xj)hj(x)
yields:
u(m)(xk) =N∑j=0
h(m)j (xk)u(xj) 0 ≤ k ≤ N.
If we define:
u(m) :=(u(m)(x0), u(m)(x1), . . . , u(m)(xN )
)T, u := u(0),
D(m) :=(d
(m)ij := h
(m)j (xi)
)0≤i,j≤N
,
d(m)ij := h
(m)j (xi),
8 1 Polynomial approximation on semi–infinite domains
then:
D := D(1), dij = d(1)ij .
We now state two important results. The first ensures that it is sufficient to compute the first
order differentiation matrix, the second gives the general expression of its entries.
Theorem 7.
D(m) = D ·D · · ·D = Dm, m ≥ 1.
Theorem 8. The entries of D are determined by
dij :=
Q′(xi)Q′(xj)
1xi − xj
i 6= j
Q′′(xi)2Q′(xi)
i = j
where Q(x) = pN+1(x) or Q(x) = (x − a)qN (x) are the quadrature polynomials of Gauss and
Gauss–Radau rules, respectively.
Proof. Lagrange polynomials can be expressed as:
hj(x) = Q(x)Q′(xj)(x− xj)
, 0 ≤ j ≤ N.
Differentiating this expression, and using Q(xj) = 0 leads to:
dij = h′j(xi) = Q′(xi)Q′(xj)
1xi − xj
∀i 6= j.
Furthermore, we have:
dii = limx→xi
h′i(x) = 1Q′(xi)
limx→xi
Q′(x)(x− xi)−Q(xi)(x− xi)2 = Q′′(xi)
2Q′(xi).
ut
The former result implies that:
u(m) = Dmu, m ≥ 1,
the latter that differentiation in the physical space can be performed through matrix–matrix and
matrix–vector products.
1.7 The case of semi–infinite domains
In this section we reformulate the previous results in the case of the semi–infinite interval I =
(0,+∞). We review the basic properties of classical orthogonal systems in this interval (Laguerre
polynomials and functions) and the corresponding results on interpolation, discrete transforms,
quadrature rules and spectral differentiation. We refer to [She01, GS00, ST06, CHQZ06].
1.7 The case of semi–infinite domains 9
1.7.1 Laguerre polynomials
Laguerre polynomials Lnn∈N are a system of orthonormal polynomials in L2ω(R+), where ω(x) =
e−x and R+ := (0,+∞), that is:
(Ln(x),Lm(x))ω = δnm.
In this case, the three–term recurrence relation reads:
(n+ 1)Ln+1(x) = (2n+ 1− x)Ln(x)− nLn−1(x), (1.15a)
L0 = 1, L1 = 1− x.
The leading coefficient of Ln(x) is:
cn = (−1)n
n!,
and Ln(0) = 1. An explicit expression for Ln is the following(hereafter ∂x := d
dx
):
Ln(x) =n∑k=0
(−1)k
k!
(n
k
)xk, (1.16)
or
Ln(x) = 1n!
ex∂nx(xne−x
).
Laguerre polynomials satisfy the following recurrence relations:
∂xLn(x) = −n−1∑k=0
Lk(x),
Ln(x) = ∂xLn(x)− ∂xLn+1(x), (1.17a)
x∂xLn(x) = n [Ln(x)−Ln−1(x)] . (1.17b)
Furthermore, we have the following property (orthogonality of the derivatives):∫ +∞
0∂xLn(x)∂xLm(x)xω(x)dx = nδnm.
Since Lnn∈N is an orthonormal basis of L2ω(R+), any function u ∈ L2
ω(R+) can be expanded in
series of Laguerre polynomials:
Su =∞∑k=0
ukLk, SNu =N∑k=0
ukLk,
uk = (u,Lk)ω .
For large x, Ln(x) grows like its leading term (−1)n
n!xn. Moreover (see [ST06]):
|Ln(x)| ≈ x−1/4ex/2 for n 1, x > 0, (1.18)
10 1 Polynomial approximation on semi–infinite domains
that is, the expansion polynomials are unbounded for x → +∞. This fact, along with the obser-
vation that ω(x) → 0 exponentially for x → +∞, makes the approximation of u ∈ L2ω(R+) with
its truncated Laguerre series SNu less and less accurate as x → +∞. In other words, although
Laguerre polynomials can in principle approximate functions with growth at infinity, this approx-
imation often turns out to be inaccurate.
On the other hand, we are mainly interested in the approximation of functions which enjoy some
decay properties at infinity. In this setting it is more appropriate to expand the function in series
of Laguerre functions.
1.7.2 Laguerre functions
Laguerre functions Lnn∈N are defined by:
Ln(x) := e−x/2Ln(x) x ∈ R+, n ≥ 0.
They are orthogonal in L2(R+). Introducing the derivative operator:
∂x = ∂x + 12,
we have that
∂xLn(x) = ex/2∂xLn(x),
and Laguerre functions’ properties read as “hatted” Laguerre polynomials’ ones:
• Three–term relation:
(n+ 1)Ln+1(x) = (2n+ 1− x)Ln(x)− nLn−1(x),
L0 = e−x/2, L1 = (1− x)e−x/2.
• Recurrence formulas:
∂xLn(x) = −n−1∑k=0
Lk(x),
Ln(x) = ∂xLn(x)− ∂xLn+1(x),
x∂xLn(x) = n[Ln(x)− Ln−1(x)
].
• Orthogonality of the derivatives:∫ +∞
0∂xLn(x)∂xLm(x)xdx = nδnm.
• Asymptotic behaviour:
|Ln(x)| → 0 as x→ +∞, |Ln(x)| ≤ 1 ∀x ∈ [0,+∞).
1.7 The case of semi–infinite domains 11
Since Lnn∈N is an orthonormal basis of L2(R+), any function u ∈ L2(R+) can be expanded in
series of Laguerre functions, namely:
Su =∞∑k=0
ukLk, SNu =N∑k=0
ukLk,
uk =(u, Lk
).
As usual, SNu is the N–th order truncated of Su, while uk is the k–th Laguerre coefficient.
Graphs of the first six Laguerre polynomials and functions are reported in Figure 1.1.
1 2 3 4 5 6
-6
-4
-2
2
4
6
8
L5
L4
L3
L2
L1
L0
5 10 15 20
-0.4
-0.2
0.2
0.4
0.6
0.8
1.0
L`
5
L`
4
L`
3
L`
2
L`
1
L`
0
Fig. 1.1: The first six Laguerre polynomials (left) and functions (right).
1.7.3 Approximation with Laguerre polynomials and functions
We review here some results on the approximation of functions in R+ with Laguerre polynomials
and functions. We begin with Laguerre polynomials; our discussion follows the presentation in
[GS00].
First, let us define the Sobolev spaces used in the approximation results:
Hmω (R+) :=
u | ∂kxu ∈ L2
ω(R+), 0 ≤ k ≤ m,
with the corresponding semi–norm and norm:
|u|m,ω = ‖∂mx u‖ω, ‖u‖m,ω =
(m∑k=0|u|k,ω
)1/2
,
(‖ · ‖ω is the L2ω–norm). Hereafter, we will restrict to the case m ∈ N. In view of boundary
conditions on the left endpoint, we define:
12 1 Polynomial approximation on semi–infinite domains
H10,ω(R+) :=
u | u ∈ H1
ω(R+), u(0) = 0.
As above, let PN be the space of polynomials of degree at mostN , and let P0N = PN∩H1
0,ω(R+). Let
πN : L2ω(R+)→ PN be the L2
ω–orthogonal projection operator, i.e. such that for any u ∈ L2ω(R+):
(u− πNu, ϕ)ω = 0 ∀ϕ ∈ PN
holds.
Next, for ξ ∈ N we define:
Hmω,ξ(R+) :=
u ∈ Hm
ω (R+) | xξ/2u ∈ Hmω (R+)
,
with the associated norm:
‖u‖m,ω,ξ :=∥∥∥u(1 + x)ξ/2
∥∥∥m,ω
.
The subscript ω will be omitted when (e.g. Laguerre functions) ω = 1. Furthermore, all the
constants (C,C1, C2 and so on) on the right hand side of the inequalities are generic positive
constants independent of N.
Lemma 1 (Lemma 2.1, [GS00]). Let m ≥ µ ≥ 0. Then
‖u− πNu‖m,ω ≤ CNµ−m/2‖u‖m,ω,m ∀u ∈ Hmω,m(R+),
with C ∈ R+.
Following [GS00], if we define the H1ω–projector π1
N : H1ω → PN by:
(∂x(u− π1
Nu), ∂xϕ)ω
+(u− π1
Nu, ϕ)ω
= 0 ∀ϕ ∈ PN , u ∈ H1ω(R+),
and the H10,ω–projector π
1,0N u : H1
0,ω(R+)→ P0N by:(
∂x(u− π1,0N u), ∂xϕ
)ω
= 0 ∀ϕ ∈ H10,ω(R+) ∩Hm
ω,m−1(R+),
the following result holds:
Lemma 2. Let m ≥ 1. Then:
1)∥∥u− π1
Nu∥∥
1,ω ≤ C1N(1−m)/2‖u‖m,ω,m−1 ∀u ∈ Hm
ω,m−1(R+),
2)∥∥∥u− π1,0
N u∥∥∥
1,ω≤ C2N
(1−m)/2‖u‖m,ω,m−1 ∀u ∈ H10,ω(R+) ∩Hm
ω,m−1(R+),
with C1, C2 ∈ R+.
As for Laguerre functions (orthonormal in L2(R+)), we make use of non–weighted Sobolev spaces.
Following [She01], we set:
PN =u | u = ve−x/2, v ∈ PN
, P0
N =u | u = ve−x/2, v ∈ P0
N
,
1.7 The case of semi–infinite domains 13
and
πNu = e−x/2πN
(uex/2
), π1,0
N u = e−x/2π1,0N
(uex/2
).
The last two operators are characterized by the following result:
Lemma 3.
(u− πNu, vN ) = 0 ∀vN ∈ PN , u ∈ L2(R+),(∂x
(u− π1,0
N u), ∂xvN
)+ 1
4
(u− π1,0
N u, vN
)= 0 ∀vN ∈ P0
N , u ∈ H10 (R+).
Next, we have the following result (cfr. [She01]):
Lemma 4. Let m ≥ 1. Then:
1) ‖u− πNu‖µ ≤ C1Nµ−m/2‖u‖m,m ∀u ∈ Hm
m (R+),
2)∥∥∥u− π1,0
N u∥∥∥
1≤ C2N
(1−m)/2∥∥∥∂xu∥∥∥
m−1,m−1∀u ∈ H1
0 (R+) s.t. ∂xu ∈ Hm−1m−1 (R+).
where 0 ≤ µ ≤ m, C1, C2 ∈ R+, ∂x(·) = ∂x(·)+ 12(·) and the subscript ω has been omitted in Hm
ω,m
since ω ≡ 1 for Laguerre functions.
1.7.4 Gauss–Laguerre quadrature formulas
We specify here the results for Gaussian quadrature to the case of the integration over [0,+∞) with
weight function e−x. In this setting, we talk about Gauss–Laguerre (GL) and Gauss–Laguerre–
Radau (GLR) nodes and weights.
The following result holds (see [She01, SW09]):
Theorem 9 (Gauss–Laguerre quadrature). Let xjNj=0, ωjNj=0 be the set of GL or GLR quadra-
ture nodes and weights. Then,∫ +∞
0p(x)e−xdx =
N∑j=0
p(xj)ωj ∀p ∈ P2N+δ,
where δ = 1, 0 for GL and GLR quadrature, respectively.
• For GL quadrature:
xjNj=0 are the zeros of LN+1(x);
ωj = xj(N + 1)2[LN (xj)]2
, 0 ≤ j ≤ N.
• For GLR quadrature:
x0 = 0, xjNj=1 are the zeros of ∂xLN+1(x);
14 1 Polynomial approximation on semi–infinite domains
ωj = 1(N + 1)[LN (xj)]2
, 0 ≤ j ≤ N.
We see that the weights are exponentially small for large xj . As for GL and GLR quadrature
rules based on Laguerre functions (Gaussian quadrature in [0,+∞) with ω(x) = 1), we have the
following result (see [ST06]):
Theorem 10 (Modified Gauss–Laguerre quadrature). Let xjNj=0,ωjNj=0 be the nodes and
weights given in Theorem 9. Denote:
ωj := exjωj , 0 ≤ j ≤ N. (1.19)
Then, the following quadrature formula holds:∫ +∞
0p(x)dx =
N∑j=0
p(xj)ωj ∀p ∈ P2N+δ,
where δ = 1, 0 for the modified GL and GLR quadrature, respectively.
The nodes are the same as in the polynomial case. Moreover, it turns out that to obtain the
weights ωj it is sufficient to replace ∂x with ∂x in Theorem 9. For instance, the modified GLR
weights are:
ωj = 1(N + 1)[LN (xj)]2
, 0 ≤ j ≤ N,
which, unlike ωj, for N 1 behave like (cfr. Eq. (1.18), see also [ST06])
ωj =√xj
N + 1.
1.7.5 Computation of nodes and weights
The general setting seen in Section 1.4 can be adapted to the semi–infinite case as follows. The
zeros xjNj=0 of LN+1(x) are the eigenvalues of the symmetric tridiagonal matrix:
AN+1 =
α0√β1
√β1 α1
√β2
. . . . . . . . .√βN−1 αN−1
√βN
√βN αN
,
with:
αj = 2j + 1 0 ≤ j ≤ N,
βj = j2 1 ≤ j ≤ N.
1.7 The case of semi–infinite domains 15
The weights ωjNj=0 can be computed as the first component of the orthonormal eigenvectors of
AN+1. However, for the sake of numerical stability, it is preferable to compute the ωj ’s from the
ωj ’s by means of the inverse of (1.19), whose computation is always stable, that is:
ωj = e−xj ωj , 0 ≤ j ≤ N.
1.7.6 Interpolation and discrete transforms
Following [ST06], we specify the results seen above on polynomial interpolation and discrete trans-
forms to the case of the unbounded interval [0,+∞). We restrict ourselves to Laguerre functions,
though the discussion holds also for Laguerre polynomials.
Let xjNj=0,ωjNj=0 be the nodes and weights as in Theorem 10. The corresponding discrete inner
product and norm are:
(u, v)N,ω =N∑j=0
u(xj)v(xj)ωj , ‖u‖N,ω = (u, u)1/2N,ω
It follows from Theorem 10 that (ω ≡ 1) :
(p, q)N,ω = (p, q)ω ∀p, q s.t. pq ∈ P2N+δ,
where, as above, δ = 1, 0 for modified GL and GLR rules, respectively.
The interpolation operator IN : C ([0,+∞))→ PN is defined by:
INu(x) =N∑n=0
unLn(x) ∈ PN ,
the coefficients unNn=0 being determined by the forward discrete transform:
un =N∑j=0
u(xj)Ln(xj)ωj 0 ≤ j ≤ N,
while the backward discrete transform enables to pass from the frequency space representation
un to the physical values u(xj).
In order to study the interpolation error, we introduce the following weighted space (see [SW09]):
Bm(R+) :=u | xk/2∂kxu ∈ L2
ω(R+), 0 ≤ k ≤ m,
where ω(x) = e−x. Bm is equipped with the following norm and semi–norm:
‖u‖Bm :=
(m∑k=0
∥∥∥xk/2∂kxu∥∥∥2
ω
)1/2
, |u|Bm :=∥∥∥xm/2∂mx u
∥∥∥ω.
In a similar way, we define:
Bm(R+) :=u | xk/2∂kxu ∈ L2(R+), 0 ≤ k ≤ m
,
16 1 Polynomial approximation on semi–infinite domains(∂x = ∂x + 1
2
), and the corresponding norm and semi–norm:
‖u‖Bm
:=
(m∑k=0
∥∥∥xk/2∂kxu∥∥∥2)1/2
, |u|Bm
:=∥∥∥xm/2∂mx u
∥∥∥ .Note that, unlike the space Hm
ω , the space Bm has different weight functions for different derivative
orders. Next, for the interpolation estimates we also need the two following spaces:
−Bm(R+) :=u | x(k−1)/2∂kxu ∈ L2
ω(R+), 0 ≤ k ≤ m,
−Bm(R+) :=u | x(k−1)/2∂kxu ∈ L2(R+), 0 ≤ k ≤ m
,
with norms and semi–norms:
‖u‖−Bm :=
(m∑k=0
∥∥∥x(k−1)/2∂kxu∥∥∥2
ω
)1/2
, |u|−Bm :=∥∥∥x(m−1)/2∂mx u
∥∥∥ω,
and:
‖u‖−Bm :=
(m∑k=0
∥∥∥x(k−1)/2∂kxu∥∥∥2)1/2
, |u|−Bm := ‖x(m−1)/2∂mx u‖.
The interpolation operator IN : C ([0,+∞))→ PN is defined in terms of IN by the relation:
IN
(uex/2
)= ex/2
(INu
)(x).
The following two results hold (see [GWW06, WGW09]):
Theorem 11. If u ∈ −Bm(R+) ∩Bm(R+), then for m ≥ 1,m ∈ N:
‖u− INu‖ω ≤ CN(1−m)/2
[|u|−Bm + 2
√lnN |u|Bm
].
Moreover:
|u− INu|1,ω ≤ CN(3−m)/2
[(1 +N−1/2)|u|−Bm + 2
√lnN |u|Bm
]In the case of the GLR nodes, the latter estimate holds if u ∈ −Bm−1(R+).
Theorem 12. If u ∈ −Bm(R+) ∩ Bm(R+), then for m ≥ 1,m ∈ N:∥∥∥u− INu∥∥∥ ≤ CN (1−m)/2[|u|−Bm + 2
√lnN |u|
Bm
].
Moreover, if ∂xu ∈ Bm−1(R+), then:∣∣∣u− INu∣∣∣1≤ CN (3−m)/2
[ (1 +N−1)N−1/2 |∂xu|Bm−1 + |u|−Bm + 2
√lnN |u|
Bm
].
The proof of the latter follows directly from the former’s.
1.7 The case of semi–infinite domains 17
1.7.7 Laguerre spectral differentiation
We specify to the Laguerre polynomials and functions case the general setting seen above for
spectral differentiation (see [ST06]).
As for Laguerre polynomials, let xjNj=0 be the GL or GLR nodes and u ∈ PN . Let hj(x)Nj=0
be the Lagrange interpolation polynomials relative to xjNj=0. From Theorem 7 we have:
u(m) = Dmu, m ≥ 1,
where:
D = (dij)i,j=0,...,N =(h′j(xi)
)i,j=0,...,N ,
u(m) =(u(m)(x0), u(m)(x1), . . . , u(m)(xN )
)T, u = u(0).
In particular, we can state the following result:
Theorem 13. The entries of the differentiation matrix D associated to the GL and GLR points
xjNj=0 have the following form:
• GL points:
dij =
LN (xi)(xi − xj)LN (xj)
if i 6= j
xi −N − 22xi
if i = j
• GLR points:
dij =
LN+1(xi)(xi − xj)LN+1(xj)
if i 6= j
12
if i = j 6= 0
−N2
if i = j = 0
Proof. Following Theorem 8, for GL points (zeros of LN+1(x)) we have Q(x) = LN+1(x). There-
fore, if i 6= j we have:
dij = Q′(xi)Q′(xj)
1xi − xj
= LN (xi)LN (xj)
1xi − xj
,
where we used the fact that LN+1(xi) = 0 and the recurrence relation (1.17b).
Next, if i = j:
18 1 Polynomial approximation on semi–infinite domains
dii = Q′′(xi)2Q′(xi)
= (∂xLN+1(x))′
2∂xLN+1(x)
∣∣∣x=xi
(1.17b)=−x(N+1
x LN (x))′
2x∂xLN+1(x)
∣∣∣x=xi
=x(− 1
x2 LN (x) + 1x∂xLN (x))
2LN (xi)
∣∣∣x=xi
(1.17b)= −LN (x) +NLN (x)−NLN−1(x)2xLN (xi)
∣∣∣x=xi
(1.15a)= xi −N − 22xi
.
As for GLR points (0 and the zeros of ∂xLN+1(x)), we have Q(x) = LN+1(x)−LN (x). Therefore,
if i 6= j we have:
dij = Q′(xi)Q′(xj)
1xi − xj
= ∂xLN (xi)∂xLN (xj)
1xi − xj
= LN+1(xi)LN+1(xj)
1xi − xj
,
where we used the fact that ∂xLN+1(xi) = 0, and relations (1.17a), (1.17b).
Next, if i 6= j we have:
dii = Q′′(xi)2Q′(xi)
= (∂xLN+1(x)− ∂xLN (x))′
2(∂xLN+1(x)− ∂xLN (x))
∣∣∣x=xi
(1.17a)= ∂xLN (xi)2LN (xi)
.
If i = j = 0 then ∂xLN (x0) = ∂xLN (0) = −N , LN (0) = 1 so that dii = −N2 ; on the other hand,
if i = j 6= 0 from (1.17a) we have ∂xLN (xi) = LN (xi) so that dii = 12 , which completes the
proof. ut
As for Laguerre functions, given u ∈ PN we have:
u(x) =N∑j=0
u(xj)hj(x),
hj
Nj=0
being the Lagrange interpolation functions defined by:
hj ∈ PN , hj(xi) = δij 0 ≤ i, j ≤ N.
In particular, if hj is as above, we have that:
hj(x) = e−x/2
e−xj/2hj(x),
and therefore:
h′j(xi) = e−xi/2
e−xj/2
(h′j(xi)−
12δij
).
In view of this fact, we can give the following result, whose proof follows straightforwardly from
Theorem 13:
1.8 Rescaling 19
Theorem 14. The entries of the differentiation matrix D associated to GL and GLR points
xjNj=0 have the following form:
• GL points:
dij =
LN (xi)(xi − xj)LN (xj)
if i 6= j
−N + 22xi
if i = j
• GLR points:
dij =
LN+1(xi)(xi − xj)LN+1(xj)
if i 6= j
0 if i = j 6= 0
−N + 12
if i = j = 0
1.8 Rescaling
In this section we summarize the main results concerning a more general type of orthogonal
polynomials (and functions) in R+. In particular, it is shown (see e.g. [WGW09]) that using the
weight function ωβ = e−βx, β > 0 (shortcut for ωα,β = xαe−βx in the case of α = 0) enables to
approximate solutions of differential equations with different asymptotic behaviours. In our case,
we will set β = 1L , where L stands for a typical length scale for environmental applications, so
that our interest will be for β 1.
We review the definition of scaled Laguerre polynomials and functions, their properties and the
main results on approximation, interpolation and quadrature rules, referring to [GW07, WGW09].
1.8.1 Scaled Laguerre polynomials
Let ωβ(x) = e−βx, β > 0. Scaled Laguerre polynomials (SLPs)
L(β)n
n∈N
are a system of
orthogonal polynomials in L2ωβ
(R+), that is:∫ +∞
0L (β)n (x)L (β)
m (x)ωβ(x)dx = 1βδnm.
Referring to section 1.7.1, we have that:
L (β)n (x) = Ln(βx) n ∈ N. (1.20)
SLPs’ properties follow from those of Lnn∈N. In particular, the three–term recurrence relation
reads:
20 1 Polynomial approximation on semi–infinite domains
(n+ 1)L (β)n+1(x) = (2n+ 1− βx)L (β)
n (x)− nL(β)n−1(x), (1.21a)
L(β)0 = 1, L
(β)1 = 1− βx.
An explicit expression for L(β)n is therefore:
L (β)n (x) =
n∑k=0
(−β)k
k!
(n
k
)xk, (1.22)
or
L (β)n (x) = 1
n!eβx∂nx (xne−βx) n ∈ N,
with:
L (β)n (0) = 1 ∀n ∈ N.
Moreover, we have the following recurrence relations (for n ≥ 1):
∂xL(β)n (x) = −β
n−1∑k=0
L(β)k (x),
L (β)n (x) = 1
β
(∂xL
(β)n (x)− ∂xL (β)
n+1(x)), (1.23a)
x∂xL(β)n (x) = n
[L (β)n (x)−L
(β)n−1(x)
]. (1.23b)
Since L (β)n n∈N is an orthonormal basis of L2
ωβ(R+), any function u ∈ L2
ωβ(R+) can be expanded
in series of SLPs:
Su =∞∑k=0
ukL(β)k , SNu =
N∑k=0
ukL(β)k ,
uk = β(u,L
(β)k
)ωβ.
1.8.2 Scaled Laguerre functions
Scaled Laguerre functions (SLFs) L (β)n n∈N are defined by:
L (β)n (x) := e−βx/2L (β)
n (x) n ∈ N. (1.24)
They are orthogonal in L2(R+), that is:∫ +∞
0L (β)n (x)L (β)
m (x)dx = 1βδnm.
SLFs’ properties follow from SLPs’:
• Three–term relation:
(n+ 1)L (β)n+1 (x) = (2n+ 1− βx)L (β)
n (x)− nL(β)n−1 (x), n ≥ 1
L(β)
0 = e−βx/2, L(β)
1 = (1− βx)e−βx/2.
1.8 Rescaling 21
• Recurrence formulas:
∂xL(β)n (x) + 1
2βL
(β)n (x) = −β
∑n−1k=0 L
(β)k (x),
12β[L
(β)n (x) + L
(β)n+1 (x)
]= ∂xL
(β)n (x)− ∂xL (β)
n+1 (x),
x∂xL(β)n (x) = 1
2
[(n+ 1)L (β)
n+1 (x)− L(β)n (x)− nL
(β)n−1 (x)
].
Since L (β)n n∈N is an orthogonal basis of L2(R+), any function can be expanded in series of
SLFs, namely:
Su =∞∑k=0
ukL(β)k , SNu =
N∑k=0
ukL(β)k ,
uk = β(u, L
(β)k
),
where SNu and uk are the N–th order truncated of Su and the k–th Fourier–Laguerre coefficient,
and (·, ·) is the dot product in L2(R+).
1.8.3 Approximation with SLPs and SLFs
Following [WGW09], we recall some results on the approximation with SLPs and then proceed
with the SLFs. We remark that the estimates will contain negative powers of βN = N/L, where
L is a characteristic length scale. Therefore, since N/L ≈ 1/∆x, the estimates will make sense for
a sufficiently fine grid spacing, given by the combination of the scaling factor β and the number of
nodes N (see below, section 5.1, for a discussion on the scaling properties of SLFs). Moreover, in
approximating large domains, i.e. using small β, one has to be sure to employ a sufficient number
of nodes in order to maintain the approximation properties of scaled base functions.
Let π(β)N : L2
ωβ(R+)→ PN be the L2
ωβ–orthogonal projector, i.e. such that for any u ∈ L2
ωβ(R+):(
u− π(β)N u, ϕ
)ωβ
= 0 ∀ϕ ∈ PN .
We introduce here the non–uniformly weighted Sobolev space:
Bmβ (R+) :=u | xk/2∂kxu ∈ L2
ωβ(R+), 0 ≤ k ≤ m
.
Note that for β = 1 we have Bmβ = Bm (see subsection 1.7.6). Bmβ is equipped with the following
norm and semi–norm:
‖u‖Bmβ
:=
(m∑k=0‖xk/2∂kxu‖2
ωβ
)1/2
, |u|Bmβ
:= ‖xm/2∂mx u‖2ωβ.
We have the following result:
22 1 Polynomial approximation on semi–infinite domains
Lemma 5. For any u ∈ Bmβ (R+), m ∈ N, 0 ≤ µ ≤ m:∥∥∥u− π(β)N u
∥∥∥Bµβ
≤ C(βN)(µ−m)/2 |u|Bµβ.
Above and hereafter C denotes a positive constant independent of N , β and any function.
As for SLFs (orthogonal in L2(R+)), we set (see [WGW09]):
P(β)N u :=
u | u = ve−βx/2, v ∈ PN
,
and
π(β)N u = e−βx/2π
(β)N
(ueβx/2
)∈ P(β)
N .
Since: (u− π(β)
N u, vN
)=(π
(β)N
(ueβx/2
)− ueβx/2, vNeβx/2
)ωβ
= 0 ∀vN ∈ P(β)N ,
π(β)N u is the L2(R+) orthogonal projector, π(β)
N : L2(R+)→ P(β)N .
Next, we define the derivative operator:
∂x,β = ∂x + β
2,
and introduce the Sobolev space:
Bmβ (R+) :=u | xk/2∂kx,βu ∈ L2(R+), 0 ≤ k ≤ m
.
Note that for β = 1 we have Bmβ = Bm (see subsection 1.7.6). Bmβ is equipped with the following
norm and semi-norm:
‖u‖Bmβ
:=
(m∑k=0
∥∥∥xk/2∂kx,βu∥∥∥2)1/2
, |u|Bmβ
:=∥∥∥xm/2∂mx,βu
∥∥∥2.
A first result is the following:
Lemma 6. For any u ∈ Bmβ (R+), µ,m ∈ N, 0 ≤ µ ≤ m∥∥∥u− π(β)N u
∥∥∥Bµβ
≤ C(βN)(µ−m)/2 |u|Bmβ
.
However, in view of the approximation of solutions of boundary–value problems, we need an
estimate in Hm(R+), which we give after the following inverse inequality between the L2 and Hm
norms in P(β)N , referring for the proofs to [WGW09]:
Lemma 7. For any ϕ ∈ P(β)N and m ≥ 0:
‖ϕ‖m ≤ C(βN)m‖ϕ‖.
Lemma 8. If ∂kxu ∈ Bm−kβ (R+) for µ,m ∈ N, k = 0, 1, . . . , µ, then∣∣∣u− π(β)N u
∣∣∣µ≤ C
(β−1/2 (βN)−µ/2
)(βN)µ−m/2
µ∑k=1
∣∣∂kxu∣∣Bm−kβ
,
where | · |m denotes the semi–norm in Hm(R+).
1.8 Rescaling 23
1.8.4 SGL quadrature formulas
If we use ωβ(x) instead of ω(x) as a weight function, we have the SGL (Scaled Gauss–Laguerre) and
SGLR (Scaled Gauss–Laguerre–Radau) nodes and weights. Theorems 9 and 10 must be modified
as follows:
Theorem 15. Letx
(β)j
Nj=0
,ω
(β)j
Nj=0
be the SGL or SGLR quadrature nodes and weights.
Then: ∫ +∞
0p(x)e−βxdx =
N∑j=0
p(x
(β)j
)ω
(β)j ∀p ∈ P2N+δ,
where δ = 1, 0 for SGL and SGLR quadrature, respectively.
• For SGL quadrature:x
(β)j
Nj=0
are the zeros of L(β)N+1(x);
ω(β)j =
x(β)j
(N + 1)2[L
(β)N (x(β)
j )]2 , 0 ≤ j ≤ N .
• For SGLR quadrature:
x(β)0 = 0,
x
(β)j
Nj=1
are the zeros of ∂xL (β)N+1(x);
ω(β)j = 1
β(N + 1)[L
(β)N (x(β)
j )]2 , 0 ≤ j ≤ N .
Next, as for SGL and SGLR quadrature rules based on SLFs, we have:
Theorem 16. Letx
(β)j
Nj=0
,ω
(β)j
Nj=0
be the nodes and weights given in Theorem 15. Denote:
ω(β)j = eβx
(β)j ω
(β)j , 0 ≤ j ≤ N.
Then: ∫ +∞
0p(x)dx =
N∑j=0
p(x
(β)j
)ω
(β)j ∀p ∈ P(β)
2N+δ,
where δ = 1, 0 for the modified SGL and SGLR quadrature, respectively.
As for the computation of nodes and weights, Theorem 6 has to be modified as follows: the zeros
x(β)j Nj=0 of L
(β)N+1(x) are the eigenvalues of the symmetric tridiagonal matrix:
AN+1 =
α0√β1
√β1 α1
√β2
. . . . . . . . .√βN−1 αN−1
√βN
√βN αN
,
24 1 Polynomial approximation on semi–infinite domains
with:
αj = 2j + 1β
0 ≤ j ≤ N,
βj = j2
β1 ≤ j ≤ N.
The weights ω(β)j Nj=0 can be computed as the first component of the orthonormal eigenvectors
of AN+1. However, for the sake of numerical stability, it is preferable to compute the ω(β)j ’s from
the ω(β)j ’s by means of the inverse of (1.19), whose computation is always stable, that is:
ω(β)j = e−βx
(β)j ω
(β)j 0 ≤ j ≤ N.
1.8.5 Interpolation
We specify here the results given in subsection 1.7.6 to the scaled case. As above, we restrict
ourselves to SLFs, though the results also hold for SLPs.
Letx
(β)j
Nj=0
,ω
(β)j
Nj=0
be the nodes and weights as in Theorem 16. The corresponding discrete
inner product and norm are:
(u, v)(β)N =
N∑j=0
u(x
(β)j
)v(x
(β)j
)ω
(β)j , ‖u‖(β)
N =[(u, u)(β)
N
]1/2(1.25)
It follows from Theorem 16 that:
(p, q)(β)N = (p, q) ∀p, q s.t. pq ∈ P(β)
2N+δ,
where δ = 1, 0 for modified SGL and SGLR rules, respectively.
In particular,
‖u‖(β)N = ‖u‖ ∀u ∈ P(β)
N .
The interpolation operator I(β)N : C([0,+∞))→ P(β)
N is defined by:
I(β)N u(x) =
N∑n=0
u(β)n L (β)
n (x) ∈ P(β)N ,
the coefficients u(β)n Nn=0 being determined by the forward discrete transform:
u(β)n =
N∑j=0
u(x
(β)j
)L (β)n
(x
(β)j
)ω
(β)j 0 ≤ j ≤ N,
while the backward discrete transform enables to pass from the frequency space u(β)n to the
physical values u(x(β)j ).
We define here the two spaces:
−Bmβ (R+) :=u | x(k−1)/2∂kxu ∈ L2
ωβ(R+), 0 ≤ k ≤ m
,
−Bmβ (R+) :=u | x(k−1)/2∂kx,βu ∈ L2(R+), 0 ≤ k ≤ m
,
1.8 Rescaling 25
with norms and semi–norms:
‖u‖−Bmβ
:=
(m∑k=0
∥∥∥x(k−1)/2∂kxu∥∥∥2
ωβ
)1/2
, |u|−Bmβ
:=∥∥∥x(m−1)/2∂mx u
∥∥∥ωβ,
and:
‖u‖−Bmβ
:=
(m∑k=0
∥∥∥x(k−1)/2∂kx,βu∥∥∥2)1/2
, |u|−Bmβ
:=∥∥∥x(m−1)/2∂mx,βu
∥∥∥ .The interpolation operator I(β)
N : C ([0,+∞))→ P(β)N is defined in terms of I(β)
N by:
I(β)N
(ueβx/2
)= eβx/2
(I
(β)N u
)(x).
The following two results hold for SGL and SGLR–based interpolations (see [GWW06, WGW09]):
Theorem 17. [SLP interpolation]
If u ∈ −Bmβ (R+) ∩Bmβ (R+) then for m ≥ 1,m ∈ N:∥∥∥u− I(β)N u
∥∥∥ωβ≤ C (βN)(1−m)/2
[β−1 |u|−Bm
β+(1 + β−1/2
)√lnN |u|Bm
β
].
Moreover:∣∣∣u− I(β)N u
∣∣∣1,ωβ≤ C(βN)(3−m)/2
[β−1
(1 +N−1/2
)|u|−Bm
β+(1 + β−1/2
)√lnN |u|Bm
β
].
In the case of the SGLR nodes, the latter estimate holds if u ∈ −Bm−1β (R+).
Theorem 18. [SLF interpolation]
If u ∈ −Bmβ (R+) ∩ Bmβ (R+) then for m ≥ 1,m ∈ N:∥∥∥u− I(β)N u
∥∥∥ ≤ C (βN)(1−m)/2[β−1|u|−Bm
β
+(1 + β−1/2
)√lnN |u|
Bmβ
].
Moreover, if ∂xu ∈ Bm−1β (R+), then:
∣∣∣u− I(β)N u
∣∣∣1≤ C (βN)(3−m)/2 [ (
β−1 + (βN)−1)N−1/2 |∂xu|Bm−1β
+ β−1 |u|−Bmβ
+(1 + β−1/2
)√lnN |u|
Bmβ
].
1.8.6 Scaled Laguerre spectral differentiation
We generalize to the case β 6= 1 the results on spectral differentiation seen in subsection 1.7.7,
giving the expression of the entries of the differentiation matrices relative to SLPs and SLFs in
the following two results, which, for β = 1, are equivalent to the ones seen above. We begin with
the SLPs.
Letx
(β)j
Nj=0
be the SGL or SGLR nodes and u ∈ P(β)N . Let (hβ)j(x)Nj=0 be the Lagrange
interpolation polynomials relative tox
(β)j
Nj=0
. From Theorem 7 we have:
26 1 Polynomial approximation on semi–infinite domains
u(m)β = Dm
β uβ , m ≥ 1,
where:
Dβ =[(dβ)ij
]i,j=0,...,N
=[(hβ)′j
(x
(β)i
)]i,j=0,...,N
,
u(m)β =
(u
(m)β (x0), u(m)
β (x1), . . . , u(m)β (xN )
)T, uβ = u
(0)β .
Next, we have:
Theorem 19. The entries of the differentiation matrix Dβ associated to the SGL and SGLR
pointsx
(β)j
Nj=0
have the following form:
• SGL points:
(dβ)ij =
L(β)N
(x
(β)i
)(x
(β)i − x(β)
j
)L
(β)N
(x
(β)j
) if i 6= j
βx(β)i −N − 2
2x(β)i
if i = j
• SGLR points:
(dβ)ij =
L(β)N+1
(x
(β)i
)(x
(β)i − x(β)
j
)L
(β)N+1
(x
(β)j
) if i 6= j
β
2if i = j 6= 0
−βN2
if i = j = 0
The proof is similar to the one of Theorem 13 and uses recurrence relations (1.21a), (1.23a) and
(1.23b).
We now turn to the SLFs case. Given u ∈ P(β)N we have:
u(x) =N∑j=0
u(x
(β)j
)(hβ
)j(x),
(hβ
)j
Nj=0
being the Lagrange interpolation functions defined by:
(hβ
)j∈ P(β)
N ,(hβ
)j
(x
(β)i
)= δij 0 ≤ i, j ≤ N.
In particular, if (hβ)j is as above, we have that:(hβ
)j(x) = e−βx/2
e−βx(β)j/2
(hβ)j (x),
1.8 Rescaling 27
and therefore: (hβ
)′j
(x
(β)i
)= e−βx
(β)i/2
e−βx(β)j/2
[(hβ
)′j
(x
(β)i
)− β
2δij
].
Theorem 19 has the following generalization:
Theorem 20. The entries of the differentiation matrix Dβ associated to SGL and SGLR pointsx
(β)j
Nj=0
have the following form:
• SGL points:
(dβ
)ij
=
L(β)N
(x
(β)i
)(x
(β)i − x(β)
j
)L
(β)N
(x
(β)j
) if i 6= j
−N + 22x(β)
i
if i = j
• SGLR points:
(dβ
)ij
=
L(β)N+1
(x
(β)i
)(x
(β)i − x(β)
j
)L
(β)N+1
(x
(β)j
) if i 6= j
0 if i = j 6= 0
−βN + 12
if i = j = 0
Again, the proof is similar to the one of Theorem 14. In the previous sections we have reviewed
the main results of polynomial approximation and Gaussian integration in R+. These results will
be useful to approximate differential boundary–value problems.
2
Free–surface flow in shallow water assumptions
In this chapter we define the model equations we will consider for discretization on a semi–
infinite domain. The shallow water equations describe waves arising from perturbations of the
free surface of a shallow layer of fluid within a constant gravitational field.
Classical derivations can be found in [Sto92, Whi99], where the equations are obtained with an
approximation to the lowest order of a development of the variables in series of the small param-
eter ε, the ratio of the water depth to the longitudinal length scale.
A recent, in–depth analysis can be found in [DBMS09], where a second order approximation
of section–averaged shallow water equations is derived from Reynolds–Averaged Navier–Stokes
model, taking into account eddy viscosity, bottom friction and variable channel geometry.
As for our work, we follow the presentation in [Tor01], where the shallow water model is de-
rived from Navier–Stokes conservation–of–mass and conservation–of–momentum equations for an
inviscid, incompressible flow.
2.1 Geometrical setup and notation
We consider a Cartesian frame of reference (x, y, z) in R3, with the time variable denoted by t. For
a general function dependent on space and time f = f(x, y, z, t), we denote its time derivative by:
ft = ∂f
∂t= ∂tf,
and its derivative with respect to x by:
fx = ∂f
∂x= ∂xf.
(analogous symbols for y and z). Next, given two vectors u = (u1, u2, u3), v = (v1, v2, v3) in R3,
their scalar product is denoted by:
u · v =3∑i=1
uivi.
The gradient of a scalar f = f(x, y, z) ∈ R is the vector field defined by:
30 2 Free–surface flow in shallow water assumptions
gradf = ∇f := (fx, fy, fz) ,
whereas the divergence of a vector field V ∈ R3 is the scalar defined by:
divV = ∇ · V := ∂v1
∂x+ ∂v2
∂y+ ∂v3
∂z.
Since our interest is directed to fluid motion, we also define the Lagrangian or material derivative
for a quantity f = f(x, y, z, t) by:DfDt
:= ∂f
∂t+ V · ∇f, (2.1)
where V is the fluid velocity. (2.1) represents the rate of change of f along the paths of fluid
particles passively advected in the flow field with velocity V .
2.2 Model derivation
The conservation of mass and conservation of momentum equations for a compressible flow can
be written as: ρt +∇ · (ρv) = 0
∂
∂t(ρv) +∇ · [ρv ⊗ v +Σ] = ρf
(2.2)
Here ρ is the fluid density, v = (u, v, w) = (v1, v2, v3) the fluid velocity, ρv the momentum, Σ the
stress tensor, f = (f1, f2, f3) represents the external forces acting on the fluid, and:
v ⊗ v ∈M3(R), (v ⊗ v)ij := vivj , ∀i, j = 1, 2, 3.
The (2.2) are four equations in the unknowns ρ, u, v, w. If we neglect viscous stresses, a reasonable
assumption for the phenomena we take into account, the stress tensor reduces to the pressure
acting on the fluid, namely:
∇ ·Σ = ∇ · (pI) = ∇p,
where I is the identity tensor, Iij := δij . With this hypothesis, equations (2.2) read component-
wise as follows:
ρt + uρx + vρy + wρz + ρ (ux + vy + wz) = 0
ut + uux + vuy + wuz = −1ρpx + f1
vt + uvx + vvy + wvz = −1ρpy + f2
wt + uwx + vwy + wwz = −1ρpz + f3
(2.3)
2.2 Model derivation 31
We consider the motion of a fluid with a free surface. As said above, the equations will be written in
a Cartesian right–handed three–dimensional frame of reference, in which the x and y coordinates
determine the horizontal plane, z the vertical direction. In an atmospheric setting, referring e.g.
to a plane tangent to the Earth surface in a point of the Northern Hemisphere, without loss of
generality we can set the x coordinate to point eastwards, the y northwards, and the z away from
Earth surface perpendicularly to the x–y plane.
The bottom boundary (the terrain) is defined by the function:
z = b(x, y),
the free surface by:
z = η(x, y, t) = b(x, y) + h(x, y, t),
where h(x, y, t) represents the thickness of fluid layer, i.e. the vertical distance between the terrain
and the free surface position. Neglecting the effects of friction and turbulent viscosity and assuming
constant density in a constant gravitational field, Eqs. (2.3) read:
ux + vy + wz = 0 (2.4a)
ut + uux + vuy + wuz = −1ρpx (2.4b)
vt + uvx + vvy + wvz = −1ρpy (2.4c)
wt + uwx + vwy + wwz = −1ρpz − g (2.4d)
where an external forces vector f = (0, 0,−g) has been assumed, thus neglecting Coriolis and
other contributions (g = 9.81 m/s2 is the constant acceleration of gravity).
In many instances of geophysical motions, the typical wavelength of the disturbances can be con-
sidered large compared with the average fluid layer thickness: this is the shallow water assumption,
which in a hydraulic context corresponds to consider a layer of water whose depth is small with
respect to the perturbation wavelength. We remark that here we are not making any assumption
about the perturbation amplitude.
Next, we impose two conditions on the free surface:
1. Kinematic condition:DDt
[z − η(x, y, t)] = 0, (2.5)
which forces a particle on the free surface to remain on it throughout the motion. Assuming
bottom impermeability, a similar condition must hold for the terrain, that is:
DDt
[z − b(x, y, t)] = 0. (2.6)
2. Dynamical condition:
32 2 Free–surface flow in shallow water assumptions
p(x, y, z, t)|z=η(x,y) = 0, (2.7)
that is, we neglect the effects of surface tension and assume a constant (without loss of gener-
ality, null) pressure at the free surface.
As for velocity variables, the model equations will be written for vertical means of u and v, defined
by:
U = 1h
∫ η
b
u(x, y, z, t)dz
V = 1h
∫ η
b
v(x, y, z, t)dz
As a consequence of the shallow water assumption, the vertical acceleration is negligible with
respect to gravity acceleration, DwDt
= 0. Plugging this condition into (2.4d) gives:
pz = −ρg (2.7)=⇒ p = ρg(η − z) (2.8)
that is, the pressure is hydrostatic (pressure acting on any point in the fluid is equal to the weight
of the fluid above it). Derivation of (2.8) with respect to x and y yields:
px = ρgηx, py = ρgηy. (2.9)
As U and V have been defined as vertically–averaged quantities, their derivatives with respect to
z are null, and therefore Equations (2.4b–2.4c) become:ut + uux + vuy = −gηx
vt + uvx + vvy = −gηy
where for simplicity we go on denoting by u (resp. v) the vertically averaged velocity U (resp. V ).
Next, we consider Equation (2.4b); with a vertical integration between the bottom b and the free
surface η, we obtain:
w|z=η − w|z=b +∫ η
b
uxdz +∫ η
b
vydz = 0. (2.10)
Kinematic conditions on free and bottom surface read:
(ηt + uηx + vηy − w) |z=η = 0,
(ubx + vby − w) |z=b = 0.
Substitution of w|z=η and w|z=b into (2.10) gives:
(ηt + uηx + vηy) |z=η − (ubx + vby) |z=b +∫ η
b
uxdz +∫ η
b
vydz = 0. (2.11)
2.2 Model derivation 33
Leibniz formula1 implies that:∫ η
b
uxdz = ∂
∂x
∫ η
b
udz − u|z=ηηx + u|z=bbx,
a similar relation holding also for the last integral in (2.10). Finally, we obtain:
ηt + ∂
∂x
∫ η
b
udz + ∂
∂y
∫ η
b
vdz = 0,
which, in view of the fact that the velocities are vertically averaged and that bt = 0 (fixed terrain),
η = b+ h, reduces to:
ht + (hu)x + (hv)y = 0. (2.13)
Moreover, multiplying by h Equations (2.4b–2.4c) and adding (2.13) multiplied by u and v, re-
spectively, we have:
(hu)t +(hu2 + 1
2gh2)x
+ (huv)y = −ghbx (2.14)
(hv)t + (huv)x +(hv2 + 1
2gh2)y
= −ghby (2.15)
Equations (2.13–2.14–2.15), the Shallow Water Equations, are a system of nonlinear conservation
laws, which can be written in the following, more compact form:
U t + F (U)x +G (U)y = S (U) (2.16)
where U = (h, hu, hv)′ is the vector of conserved variables (fluid depth and discharge), F (U) =(hu, hu2 + 1
2gh2, huv
)′and G (U) =
(hv, huv, hv2 + 1
2gh2)′
are the flux vectors and S (U) =
(−ghbx,−ghby)′ is the source term vector.
Furthermore, assuming that U is sufficiently smooth, system (2.16) can be written in quasi–linear
form, that is:
U + F ′ (U)Ux +G′ (U)Uy = S (U) , (2.17)
where F ′ := dFdU
and G′ := dGdU
are the flux Jacobian matrices (see also [LeV02]):
F ′(u) =
0 1 0
−u2 + gh 2u 0
−uv v u
G′(u) =
0 0 1
−uv 2v u
−v2 + gh 0 2v
(2.18)
Classical solution strategies for the shallow water equations begin with the study of the associated1
d
dα
∫ ξ2(α)
ξ1(α)f(ξ, α)dξ =
∫ ξ2(α)
ξ1(α)
∂f
∂αdξ + f(ξ2, α)∂ξ2
∂α− f(ξ1, α)∂ξ1
∂α(2.12)
34 2 Free–surface flow in shallow water assumptions
Riemann Problem, given by equations (2.16) or (2.17) with piecewise constant initial data. We
remark that, in contrast to the linear advection case, solutions to nonlinear hyperbolic equations
can develop finite–time discontinuities also starting from smooth initial data (see e.g. [LeV92,
LeV02]). However, the presence of frictional and forcing terms (e.g. Coriolis) usually does not allow
discontinuities to develop in environmental flows, often characterized by low Froude numbers.
2.3 Advective formulation and linearization
We review here a different form of the shallow water equations which highlights the coupling
between the free surface elevation and the fluid velocity. Although it adopts a different convention
for the depth–related variables, this formulation is widely used in the literature on numerical
hydraulics and fluid dynamics of the atmosphere (see e.g. [Cas90]).
A derivation similar to the one seen above leads to the following equations:
∂η
∂t+H
∂u
∂x+ u
∂H
∂x+H
∂v
∂y+ v
∂H
∂y= 0
∂u
∂t+ u
∂u
∂x+ v
∂u
∂y= −g ∂η
∂x
∂v
∂t+ u
∂v
∂x+ v
∂v
∂y= −g ∂η
∂y
(2.19)
Here, differently from above, η(x, y, t) is the water surface elevation measured from the undisturbed
water surface z = 0. H(x, y, t) = h(x, y) + η(x, y, t) is the total water depth, h(x, y) the water
depth measured from the undisturbed water surface. With these conventions the Shallow Water
assumption reads H λ, where λ is the wavelength of the free surface perturbation. Figure 2.1
gives a representation of the notation for the conservative and advective formulation of free–surface
flows.
If we consider free surface perturbations having also small amplitude with respect to their wave-
length, we obtain a system of linear equations. To this end, for the sake of simplicity we consider
the monodimensional version of (2.19):
∂η
∂t+ u
∂H
∂x+H
∂u
∂x= 0
∂u
∂t+ u
∂u
∂x= −g ∂η
∂x
(2.20)
To linearise these equations, with a standard procedure (see e.g. [Whi99, BK00]) we suppose that
2.4 The problem of open boundary conditions 35
x
z
b(x, y)
h(x, y, t)
z = η(x, y, t)
z = b(x, y)
0
x
z
H(x, y, t) = η + h
z = η(x, y, t)
z = −h(x, y, t)
0
Fig. 2.1: Vertical section of free–surface flow, conservative (left) and advective (right) formulation nota-
tion.
all the quantities involved are composed by a constant average part plus a small perturbation,
that is:
η(x, t) = η + η′(x, t), u(x, t) = u+ u′(x, t), H(x, t) = H + η′(x, t),
where∣∣∣∣η′η∣∣∣∣ = O(ε),
∣∣∣∣u′u∣∣∣∣ = O(ε), where ε 1 denotes a small parameter.
Plugging these expressions into (2.20) and taking into account that the average parts are constant,
we obtain:
∂η′
∂t+ H
∂u′
∂x+ u
∂η′
∂x= 0
∂u′
∂x+ (u+ u′) ∂u
′
∂x= −g ∂η
′
∂x
(2.21)
Finally, neglecting O(ε2) terms and dropping the primes for simplicity, we find:
∂η
∂t+ H
∂u
∂x+ u
∂η
∂x= 0
∂u
∂t+ u
∂u
∂x= −g ∂η
∂x
(2.22)
One can easily check (see e.g. [Sto92, Whi99, LeV02]) that the above system is equivalent to
second–order linear wave equations for both the elevation and velocity variables, whose solutions
can be expressed in the form of plane waves (combination of sinusoidal functions) propagating in
both directions at velocities ±c = u±√gH relatively to the fluid.
2.4 The problem of open boundary conditions
Having stated the mathematical formulation of free–surface flows, we review the problem of open
boundary conditions, which arises every time that wave propagation is modeled in an unbounded
36 2 Free–surface flow in shallow water assumptions
domain. Several examples include phenomena in acoustics, geophysics, oceanography and meteo-
rology.
If there is no physical boundary on which boundary conditions can be imposed, to encompass the
finite region of interest one must employ an artificial boundary B, which separates this region from
the external space. However, as the presence of B has no physical counterpart, things must be set
up so that, on one hand, B is transparent to the waves propagating in either direction, on the
other hand, that no spurious waves are generated at the boundary and no reflections of outgoing
waves occur into the finite computational domain of interest. To this end, an accurate analytical
and numerical design of B has to be carried out.
This issue has been receiving the attention of applied mathematicians’ community for the last four
decades, though, as pointed out in [MDD10], it has not yet received a general answer. Reviews on
different strategies for the solution of the problem can be found in [Ras86] and [Giv08], who classes
the difficulties encountered in constructing a consistent open boundary condition as follows:
I Well–posedness: the differential problem set in the finite domain with the open boundary
condition on B is mathematically well–posed;
I Accuracy on the continuous level: problem in the finite domain with the boundary conditions
on B is a good approximation of the original problem in the infinite domain;
I Scheme compatibility: the open boundary condition on B is compatible with the numerical
scheme used in the finite domain;
I Stability: numerical method employed including the open boundary condition on B must result
in a stable numerical scheme;
I Accuracy on the discrete level: the error generated by the numerical scheme due to the use of
the discrete open boundary conditions on B is small.
I Efficiency: the use of the open boundary conditions on B does not involve a large computational
effort;
I Ease of implementation;
I Generality.
Referring for example to atmospheric modelling, outgoing waves reflecting from the boundaries can
spread to the region of interest, corrupting the forecast (see e.g. [McD03]). Early suggestions for
boundary conditions at the artificial “ceiling” in atmospheric models include setting the vertical
velocity to zero at some finite height or, in a pressure coordinate, at p = 0 (infinite height).
However, both solutions entail total energy reflection or lead to ill–posed formulations (see [Ras86]).
Thus, it is in this context that the two most popular approaches to open boundary conditions have
been developed: radiation boundary conditions and absorbing/sponge layers. The basic idea of the
radiation or characteristic boundary conditions (see [EM77, IO81]) is that disturbances propagate
2.4 The problem of open boundary conditions 37
across B as waves that can be described by a simplification of the full model dynamics (see
also [KD83]). On the other hand, relaxation methods extend the computational domain to allow
disturbances to leave and penetrate into a limited region (the absorbing layer) where they can be
damped by viscous or reaction terms towards the external solution (see e.g. [LT08, KL78, DK83]).
If the external solution is zero, this approach amounts to introduce a simple damping term in the
layer (see also [KDH08]).
The choice of the relaxation coefficient is a trade–off between the aim to damp efficiently outgoing
perturbances on one side and limited spatial resolution, which does not allow high spatial variations
of wave amplitude, on the other. Various functional forms for the parameters have been tested,
but, as pointed out in [MDD10], the choice of this coefficients and the thickness of the layer
has been often directed by empirical criteria. Moreover, a major drawback of this approach is its
high computational cost, as considerable effort must be put into simulation of modified equations
outside the area of interest (see [Giv08]).
Following [LT08] and [MDD10], we report here a sponge layer formulation for the shallow water
equations in advective form, which read:
∂η
∂t+ ∂ (Hu)
∂x+ ∂ (Hv)
∂y= −γ(η − ηe)
∂u
∂t+ u
∂u
∂x+ v
∂u
∂y+ g
∂η
∂x= −γ(u− ue)
∂v
∂t+ u
∂v
∂x+ v
∂v
∂y+ g
∂η
∂y= −γ(v − ve)
where ηe, ue, ve are the external elevation and velocity fields, γ is the absorption coefficient, and
for simplicity we have neglected Coriolis terms. Advanced formulations (see e.g. [NNH04]) have
been developed from this simple sponge, leading to perfectly matched layers, which split variables
into contributions parallel and perpendicular to the boundary, introducing parallel and perpen-
dicular velocity sponge coefficients γx and γy. In the classical test of the evolution of an initially
motionless Gaussian elevation perturbation, the boundary treatment should allow the train of
outgoing waves to leave the domain undisturbed with no reflections.
The approach we introduce in this work, which will be explained in the sequel, adopts a different
perspective: using spectral collocation methods with scaled Laguerre functions, the equations will
be directly approximated in the semi–infinite domain out of the artificial boundary. On one hand,
this method enjoys spectral accuracy which ensures numerical consistency, stability and conver-
gence, on the other hand, the use of a small number of base functions (i.e. collocation nodes) would
significantly reduce the whole computational cost and, ultimately, yield a more efficient tool for
the treatment of open boundary conditions.
38 2 Free–surface flow in shallow water assumptions
Finally, it is worth mentioning here another approach to the approximation of problems in un-
bounded domains and to open boundary conditions, which is the infinite element method. Similar
to finite element, infinite element formulations are local, but their shape functions are chosen to
mimic the asymptotic behaviour of the solution at infinity. A clear advantage of such discretiza-
tions is the sparseness of the coefficient matrices, while stability and ill–conditioning are seen as
critical factors in determining their utility, as pointed out in the extensive review in [Ast00]. In
particular, the accuracy of an infinite element relies on the choice of the shape functions towards
infinity and the order of approximation. On the other hand, infinite elements can be used to
discretize the whole exterior domain without truncating it, and from the point of view of open
boundaries, they can be seen as a local absorbing boundary condition.
Infinite elements have been widely employed in numerical models of wave propagation phenomena
from scattering or radiating objects, see e.g. [Ger00]. In this context, usually there is also an inner
conventional finite element mesh close to the object, the choice of the location of the interface be-
tween the two meshes being critical to the accuracy properties of a coupled finite/infinite elements
discretization. We will face the same problem when we will deal with our coupled finite/semi–
infinite discretization.
3
Spectral discretization of boundary–value problems on
semi–infinite domains
In this chapter we see how to apply the theoretical results for polynomial approximation on
semi–infinite domains to discretize initial boundary value problems on the real half line with
spectral methods. First, we review the spectral G–NI and collocation method in a general setting
and give a simple discretization example. Then, to validate the method, we give a discretization
algorithm for parabolic advection–diffusion–reaction equation, based on the weak form of the
equation (G–NI method in space, θ–method in time). Next, we proceed to the discretization of
linear and nonlinear shallow water equations with a collocation scheme applied directly to the
strong form. Moreover, we also give a formulation of the equations with a damping term.
3.1 The spectral G–NI and collocation approximations
In this general presentation we restrict to steady problems, and consider the linear, monodimen-
sional case, referring to [CHQZ06, QV08]. The non–stationary case can be derived with little effort
adopting implicit time discretizations. Thus, let us consider a domain (open connected set) I ⊂ R,
and the boundary–value problem: Lu = f in I (3.1a)
Bu = 0 in ∂I
where L is a linear differential operator defined in a suitable Hilbert space V on I, and B is a set
of linear boundary differential operators on ∂I.
With a standard procedure, multiplying (3.1a) by a test function v ∈ V , we obtain:
(Lu, v) = (f, v) ∀v ∈ V,
where (·, ·) denotes the scalar product in V . Integration by parts leads to the following weak form
of the problem, as opposed to (3.1a), the strong form:
40 3 Spectral discretization of boundary–value problems on semi–infinite domains
Problem 1. Find u ∈ V s.t. Bu = 0 on ∂I and:
a(u, v) = F (v) ∀v ∈ V
where a is a bilinear form on V × V and F is a linear functional on V .
A spectral Galerkin method for the problem (3.1a) aims to find an approximation uN of u in a
finite dimensional subspace VN ⊂ V , that is:
Problem 2. Find uN ∈ VN such that:
a(uN , vN ) = F (vN ) ∀vN ∈ VN (3.2)
where the test functions are chosen in VN , too. In our case, the space VN is the space of polynomials
of degree less or equal to N in each variable on I.
Next, theGalerkin method with Numerical Integration (G–NI) is obtained replacing all the integrals
in (3.2) with Gaussian quadrature formulas, which yields:
Problem 3. Find uN ∈ VN such that BuN = 0 and:
aN (uN , vN ) = FN (vN ) ∀vN ∈ VN . (3.3)
Again, the essential boundary conditions are satisfied exactly by the numerical solution uN .
Then, let xiNi=0 ∈ I be the nodes of a Gaussian integration formula; depending on the interval,
one would consider Gauss, Gauss–Lobatto or Gauss–Radau points (see above, section 1.3 and
[CHQZ06]). Back to the strong form (3.1a), the spectral collocation solution is a function uN ∈ VNsatisfying:
LNuN = f at xi ∈ I \ ∂I (3.4a)
BNuN = 0 at xi ∈ ∂I (3.4b)
Above, LN and BN are obtained from L and B replacing every exact derivative with the pseudo–
spectral derivative, that is the derivative of the interpolating polynomial in the nodes xi (see
section 1.6): for that reason, the collocation method is also called the pseudo–spectral method.
The equation is enforced at the internal nodes, and the boundary conditions are satisfied at the
nodes lying on ∂I. Indeed, one can integrate by parts (3.4a) and obtain a weak form that, in the
case of fully Dirichlet boundary conditions, is equivalent to the G–NI formulation (3.3).
As for the algebraic formulation, let hj be the Lagrange interpolation polynomials relative to the
points xiNi=0, that is:
hj ∈ VN : hj(xi) = δij ∀i, j = 0, . . . , N.
We recall that as the hj ’s are a basis of VN , for any uN ∈ VN we can write:
3.2 Advection–diffusion–reaction equation: G–NI approximation 41
uN (x) =N∑j=0
uN (xj)hj(x) x ∈ I.
The collocation method then reads:
N∑j=0
(LNhj)(xi)uN (xj) = f(xi) ∀xi ∈ I \ ∂I
N∑j=0
(BNhj)(xi)uN (xj) = 0 ∀xi ∈ ∂I
which is a linear system for the unknowns uN (xj), j = 0, . . . , N . The matrix of this system can be
expressed as a linear combination of the pseudo–spectral derivative matrices, and is equivalent to
the G–NI matrix in the case of Dirichlet boundary conditions, obtained taking the interpolation
polynomials hi as test functions in the weak formulation (3.3). For a stability and convergence
analysis relative to the standard collocation method we refer to [GO77, CHQZ06, QV08].
3.2 Advection–diffusion–reaction equation: G–NI approximation
Next, we discretize the following unsteady advection–diffusion–reaction equation:
∂u
∂t(x, t)− ∂
∂x
(µ(x, t)∂u
∂x(x, t)
)+ b(x, t)∂u
∂x(x, t) + γ(x, t)u(x, t) = f(x, t)
x ∈ R+, t ∈ (0, T )
u(0, t) = 0, limx→+∞
u(x, t) = 0
u(x, 0) = u0(x) x ∈ R+
(3.5)
µ(x, t), b(x, t), γ(x, t) are the diffusion, advection and reaction coefficients, respectively, and for
the sake of simplicity we have assumed Dirichlet boundary conditions, homogeneous in the left
endpoint and asymptotic homogeneous at infinity.
For the space variable we use a standard G–NI method, whereas the time variable is discretized
with a finite difference approach, a brief review of which is given below. We refer to [CHQZ06,
QSS07, QV08]. In the sequel, the dependence of u upon x will be assumed. For the sake of simplicity
we will also assume the coefficients to be independent of time. Moreover, let us define the objects
we use in the sequel:
V := H10 (R+), a : V × V → R, F : V → R.
V is a Hilbert space, a a bilinear form and F a linear functional. The variational formulation for
(3.5) is the following:
42 3 Spectral discretization of boundary–value problems on semi–infinite domains
Problem 4. Find u(x, t) ∈ V such that:
(∂u
∂t(x, t), v
)+ a(u(x, t), v) = F (v) ∀v ∈ V, t ∈ (0, T )
u(0) = u0
(3.6)
where:
a(u(x, t), v) =∫ +∞
0
[µ(x)∂u
∂x(x, t)∂v
∂x+ b(x)∂u
∂x(x, t)v + γ(x)u(x, t)v
]dx
F (v) =∫ +∞
0f(x, t)vdx,
and (·, ·) is the inner product in L2(R+).
Let u0N be a suitable approximation of u0 and:
VN = H10 (R+) ∩ PN ⊂ V.
The Spectral Galerkin formulation for (3.5) is:
Problem 5. Find uN ∈ VN such that:
(∂uN∂t
(t), vN)
+ a (uN (t), vN ) = F (vN ) ∀vN ∈ VN , t ∈ (0, T )
uN (0) = u0N
(3.7)
Note that (3.7) is a semi–discretization of (3.5). To obtain a complete discretization, let us consider
a uniform mesh for the time variable t defining (see e.g. [QV08]):
tn := n∆t n = 0, 1, . . . ,N ,
where ∆t > 0 is the time step, N := [N/∆t] the integer part of N/∆t.
In a finite difference method, one constructs a general sequence unN (x) ≈ u(tn, x) replacing the
time derivatives with suitable quotients. Let us recall the fundamentals of this method considering
the following general Ordinary Differential Equation (see [QV08]):dydt
(t) = Ψ (t, y(t)) t ∈ (0, T )
y(0) = y0
(3.8)
The θ–method discretizes (3.8) as follows:yn+1 − yn
∆t= θΨ(tn+1, yn+1) + (1− θ)Ψ(tn, yn)
y0 = y0
(3.9)
3.2 Advection–diffusion–reaction equation: G–NI approximation 43
We recall that for θ = 0 we obtain the so–called Forward Euler method, while the choices θ = 1
and θ = 12are relative to the Backward Euler and Crank–Nicolson, respectively.
We can now proceed to the time discretization of (3.5). Using the Backward Euler method and
replacing the continuous inner product with the discrete ones we obtain the following discrete
G–NI formulation:
Problem 6. For any n = 0, 1, . . . ,N − 1, find unN ∈ VN such that:
(un+1N − unN∆t
, vN
)N
+ aN (un+1N , vN ) = FN (vN ) ∀vN ∈ VN
u0N = u0,N ∈ VN
(3.10)
where:
aN (u, v) : = (µ(x)u′, v′)N + (b(x)u′, v)N + (γ(x)u, v)N
=N∑i=0
µ(xi)u′(xi)v′(xi)ωi +N∑i=1
b(xi)u′(xi)v(xi)ωi
+N∑i=1
γ(xi)u(xi)v(xi)ωi,
and:
FN (v) =(fn+1, v
)N
=N∑i=1
fn+1(xi)ωi,
where we have denoted with xiNi=1 the internal GLR nodes and with u′ the derivative ∂u∂x
.
As hjNj=1 are a base of VN , the solution unN at the time step tn can be written as:
unN (x) =N∑j=1
unN (xj)hj(x).
With the choice vN = hi, i = 1, . . . , N , we obtain:(un+1N − unN∆t
, hi
)N
+ aN(un+1N , hi
)= FN (hi) ,
that is:N∑j=1
un+1N (xj)− unN (xj)
∆t(hj , hi)N +
N∑j=1
un+1N (xj)aN (hj , hi) = FN (hi) i = 1, . . . , N.
The latter expression can be written as:
Wun+1N − unN∆t
+Aun+1N = W fn+1, (3.11)
with:
44 3 Spectral discretization of boundary–value problems on semi–infinite domains
unN = (unN (xi))i=1,...,N , A = (aij)i,j=1,...,N ,
aij = aN (hj , hi)
=N∑k=0
µ(xk)h′j(xk)h′i(xk)ωk +N∑k=1
b(xk)h′j(xk)hi(xk)ωk
+N∑k=1
γ(xk)hj(xk)hi(xk)ωk
=N∑k=0
µ(xk)h′j(xk)h′i(xk)ωk + b(xi)h′j(xi)ωi + γ(xi)δijωi,
and:
f = (fi)i=1,...,N , fi = f(xi),
W = diag (ω1, . . . , ωN ) .
Note that we have h′j(xk) = djk. Referring to [QV08], for any n the expression (3.11) brings to a
linear system of the form:
Kun+1 = W(un +∆tfn+1) , (3.12)
with:
K := W +∆tA.
This system can be solved applying the LU decomposition (see e.g. [QSS07]), that allows to write
K = LU , where U (resp. L) is an upper (resp. lower) triangular matrix . Therefore the solution of
(3.12) at t = tn is equivalent to the solution of the two following triangular systems at each time
step: Ly = W
(un +∆tfn+1)
Uun+1 = y
3.2.1 Stability analysis of discrete formulation
Following the discussion in [CHQZ06], we now prove that the discrete approximation seen above
is unconditionally stable for 12≤ θ ≤ 1. We also refer to [Sal08]. For the sake of simplicity, we
restrict our analysis to the homogeneous case with constant coefficients.
As seen above, the discrete approximation (θ–method in time, G–NI in space) is:(un+1N − unN , vN
)N
+∆taN (θun+1N + (1− θ)unN , vN ) = 0 ∀vN ∈ VN
u0N = u0,N ∈ VN
(3.13)
3.3 Linear shallow water model: collocation approximation 45
where u0N = INu0, the interpolant of u0 in the GLR nodes, and VN = H1
0 (R+) ∩ PN . As GLR
quadrature formulas have degree of exactness 2N , we have that ‖u‖2N = ‖u‖2 if u ∈ PN and then,
choosing vN = θun+1N + (1− θ)unN in (3.13) and thanks to the coercivity of aN , we have:
θ∥∥un+1
N
∥∥2N
+ (1− 2θ)(un+1N , unN )N − (1− θ) ‖unN‖
2N + α∆t‖θ(un+1
N )x
+ (1− θ)(unN )x‖2 ≤ 0, (3.14)
where α is the coercivity constant of aN . From the fact that (1− 2θ) ≤ 0 and from the Cauchy–
Schwarz and Young inequalities it follows that:
(1− 2θ)(un+1N , unN )N ≥
(12− θ)(∥∥un+1
N
∥∥2N
+ ‖unN‖2N
),
and then, substituting in (3.14):
∥∥un+1N
∥∥2N
+ α∆t
2∥∥θ(un+1
N )x + (1− θ)(unN )x∥∥2 ≤ ‖unN‖
2N .
Finally, for any n ≥ 0,
‖unN‖2N + α
∆t
2
n−1∑j=0
∥∥θ(un+1N )x + (1− θ)(unN )x
∥∥2 ≤∥∥u0
N
∥∥2N
= ‖u0‖2N ,
which proves that the discrete formulation is unconditionally stable for θ ∈ [ 12 , 1].
As for the convergence estimates, we refer to [CHQZ06, QV08], just observing that there will be a
contribution of the spatial discretization (dependent on N) and another due to the temporal one
(dependent on ∆t). As we will see in the numerical tests, the typical spectral methods geometric
decay for the space error will be in general bounded by the temporal approximation error, which
tends to decay more slowly. In particular, it is well known that the discrete formulation with
the θ–method is first–order accurate in time for any θ 6= 12, second–order accurate for θ = 1
2.
Numerical results obtained with the implementation of this formulation can be found in section
5.2.
3.3 Linear shallow water model: collocation approximation
We now apply the theory reviewed above to the discretization of the shallow water equations on
the spatial and temporal domain R+× [0,+∞), starting from the linear approximation, see (2.22),
which, together with boundary and initial conditions, reads as follows:
46 3 Spectral discretization of boundary–value problems on semi–infinite domains
∂η
∂t+ H
∂u
∂x+ u
∂η
∂x+ γη = 0 x ∈ R+, t > 0
∂u
∂t+ u
∂u
∂x+ g
∂η
∂x+ γu = 0 x ∈ R+, t > 0
η(0, t) = ϕ(t), u(0, t) = ψ(t)
limx→+∞
η(x, t) = limx→+∞
u(x, t) = 0
η(x, 0) = η0(x)
u(x, 0) = u0(x)
(3.15)
In the formulation above, we have already taken into account a damping term in both equa-
tions and Dirichlet–type boundary conditions are considered: inhomogeneous at the left endpoint,
asymptotically homogeneous for both variables at infinity.
Letx
(β)0 = 0, x(β)
1 , . . . , x(β)N
be the the SGLR quadrature nodes, defined in subsection 1.8.4. The
pseudo–spectral approximation of (3.15) is:
Problem 7. ∀t > 0, find ηN (t), uN (t) ∈ P(β)N such that:
(∂ηN∂t
+ H∂uN∂x
+ u∂ηN∂x
+ γηN
)(x
(β)j , t
)= 0, ∀j = 1, . . . , N(
∂uN∂t
+ u∂uN∂x
+ g∂ηN∂x
+ γuN
)(x
(β)j , t
)= 0 ∀j = 1 . . . , N
ηN (x(β)0 , t) = ϕ(t)
uN (x(β)0 , t) = ψ(t)
ηN
(x
(β)j , 0
)= η0
(x
(β)j
)∀j = 0, . . . , N
uN
(x
(β)j , 0
)= u0
(x
(β)j
)∀j = 0, . . . , N
(3.16)
(3.16) is a semi–discretization of (3.15). To obtain a fully discrete approximation, we consider a
uniform grid for the time variable defining:
tn := n∆t, n = 0, 1, . . . ,N ,
∆t > 0 is the time step, N :=[N
∆t
].
In a finite difference approximation, we construct two sequences ηnN (x) ≈ ηN (tn, x), unN (x) ≈
uN (tn, x), and time derivatives are approximated with suitable quotients. In particular, here we
use the Implicit Euler method, which gives an unconditionally stable, first–order accurate approx-
imation (see e.g. [QV08]). Therefore, the fully discrete approximation reads:
Problem 8. ∀n = 0, 1, . . . ,N − 1, find (ηn+1N , un+1
N ) ∈ P(β)N × P(β)
N such that:
3.3 Linear shallow water model: collocation approximation 47
(ηn+1N − ηnN∆t
+ H∂un+1
N
∂x+ u
∂ηn+1N
∂x+ γηn+1
N
)(x
(β)j
)= 0 ∀j = 1, . . . , N (3.17a)(
un+1N − unN∆t
+ u∂un+1
N
∂x+ g
∂ηn+1N
∂x+ γun+1
N
)(x
(β)j
)= 0 ∀j = 1, . . . , N
ηn+1N
(x
(β)0
)= ϕn+1 (3.17b)
un+1N
(x
(β)0
)= ψn+1 (3.17c)
η0N = η0,N ∈ P(β)
N
u0N = u0,N ∈ P(β)
N
where η0,N ≈ η0, u0,N ≈ u0 are suitable approximations to the initial data, and ϕn ≈ ϕ(tn).
As seen above (see subsection 1.8.6), if a function u(x) is identified by the vector containing its
values in the SGLR nodes, i.e. u =(u(x
(β)j
))Nj=0
, the m–th order derivative of u(x) with respect
to x is approximated by the vector obtained multiplying u by the pseudo–spectral differentiation
matrix, that is:
u(m) = Dmβ u, m ≥ 1,
where we recall that:
(Dβ)ij =(hβ
)′j
(x
(β)i
)i, j = 0, . . . , N.
In particular, for the first derivative with respect to x, we have:
u′(xi) =N∑j=1
(hβ)′j(xi)u (xj) i = 1, . . . , N =⇒ u′ = (Dβ)N u,
where we split u and the matrix Dβ in the following way:
u =
u0
u
, Dβ =
(Dβ)00 (Dβ)01 · · · (Dβ)0N
(Dβ)10... (Dβ)N
(Dβ)N0
=
(Dβ)00 dβ
dβ (Dβ)N
.
In terms of its rows and columns, respectively, Dβ reads:
Dβ =
(dβ)0
(dβ)1
...
(dβ)N
=[(dβ)0 (dβ)1 (dβ)N
]
We also define the matrix D1β in the following way:
D1β =
dβ (Dβ)N
48 3 Spectral discretization of boundary–value problems on semi–infinite domains
With these conventions, we remark that
Dβ ∈MN+1, (Dβ)N ∈MN , dβ ∈ R1×N , (dβ)i ∈ R1×(N+1), dβ ∈ RN×1, (dβ)i ∈ R(N+1)×1,
i = 0, . . . , N.
As for Dirichlet boundary conditions, we strongly impose the elevation and velocity boundary
values in the left endpoint according to (3.17b–3.17c). Then the equations for the internal nodes
read:
un+1N − unN∆t
+ uD1βu
n+1N + gD1
βηn+1N + Γ un+1
N = 0
ηn+1N − ηnN∆t
+ HD1βu
n+1N + uD1
βηn+1N + Γ ηn+1
N = 0
(3.18)
having defined the diagonal matrix:
Γ ∈MN , Γ :=
γ1
γ2. . .
γN
,
where:
γi = γ(x(β)i ) ∀i = 1, . . . , N.
Highlighting the contributions of the (known) first component of the vectors we obtain:IΓ u
n+1N = unN − u∆t(Dβ)N un+1
N − u∆tdβ(un+1N )0 − g∆t(Dβ)N ηn+1
N − g∆tdβ(ηn+1N )0
IΓ ηn+1N = ηnN −∆tH(Dβ)N un+1
N −∆tHdβ(un+1N )0 − u∆t(Dβ)N ηn+1
N − u∆tdβ(ηn+1N )0
(3.19)
where (ηn+1N )0 = ηn+1
N
(x
(β)0
)and IΓ = IN +∆tΓ . Rearranging the terms, we have:
BΓ un+1N + g∆t(Dβ)N ηn+1
N = fnu
∆tH(Dβ)N un+1N +BΓ η
n+1N = fnη
(3.20)
where:
B = IΓ + u∆t(Dβ)N ,
fnu = unN −∆t[u(un+1
N )0 + g(ηn+1N )0
]dβ , (3.21)
fnη = ηnN −∆t[H(un+1
N )0 + u(ηn+1N )0
]dβ . (3.22)
Finally, we have:
3.3 Linear shallow water model: collocation approximation 49BΓ u
n+1N = −g∆t(Dβ)N ηn+1
N + fnu
FΓ ηn+1N = fnη −∆tH(Dβ)NB−1fnu
(3.23)
where:
FΓ = BΓ − g∆t2H(Dβ)NB−1(Dβ)N . (3.24)
Solution of (3.23) is found solving an N ×N linear system for ηN at each time step and replacing
the obtained solution into the equation for uN , which gives rise to a N ×N system to be solved
at each time step, too.
Next, we consider transmissive boundary conditions at the left endpoint:
∂η
∂x(x0) = ∂u
∂x(x0) = 0 (3.25)
that in a collocation sense read:
(ηn+1N )0 = − 1
(Dβ)00dβη
n+1N , (un+1
N )0 = − 1(Dβ)00
dβun+1N .
Plugging these expressions into (3.19) and defining the matrix:
D1 = (Dβ)N −1
(Dβ)00dβdβ
we obtain: B1Γ u
n+1N + g∆tD1η
n+1N = unN
∆tHD1un+1N +BΓ η
n+1N = ηnN
(3.26)
where here:
B1Γ = IΓ + u∆tD1,
and therefore: B1Γ u
n+1N = unN − g∆tD1η
n+1N
F 1Γ η
n+1N = ηnN −∆tHD1
(B1Γ
)−1unN
(3.27)
where:
F 1Γ = B1
Γ − gH∆t2D1(B1Γ
)−1D1η
n+1N .
As (3.23), (3.27) is solved with a N × N linear system for ηN at each time step, replacing the
obtained solution in the equation for uN , and solving another N ×N linear system.
Finally, let us now consider the model given by (3.15), letting u = 0 for the sake of simplicity, and
impose Neumann boundary conditions on the fluid discharge in the left endpoint, that is:
H∂u
∂x(0, t) = f(t). (3.28)
50 3 Spectral discretization of boundary–value problems on semi–infinite domains
This condition at discrete level (Problems 7 and 8) reads:
H∂uN∂x
(x(β)0 , t) = f(t) =⇒ H
∂unN∂x
(x(β)0 ) = fn ∀n = 0, . . . ,N
that is, in a collocation sense:
H (dβ)0un+1N = H
[(Dβ)00
(un+1N
)0 + dβun+1
N
]= fn+1. (3.29)
With this boundary condition we obtain the value of the free surface elevation in the left endpoint
by substitution in (3.17a), that is:
(1 + γ0∆t)(ηn+1N
)0 = (ηnN )0 −∆tf
n+1, (3.30)
where γ0 = γ(x(β)0 ). Indeed, we are left with a Dirichlet condition on the left endpoint elevation.
Then, we can obtain the left endpoint velocity by:
(1 + γ0∆t)(un+1N
)0 = (unN )0 − g∆t (Dβ)00
(ηn+1N
)0 − g∆tdβη
n+1N (3.31)
and replace it in the equations for the internal nodes. Neglecting advection, Eqs. (3.19) read:IΓ u
n+1N + g∆t(Dβ)N ηn+1
N = unN − g∆tdβ(ηn+1N
)0 (3.32a)
IΓ ηn+1N = ηnN −∆tH (Dβ)N u
n+1N −∆tHdβ
(un+1N
)0
Replacing un+1N in the second equation we obtain:
ηn+1N = ηnN −∆tHI−1
Γ (Dβ)N I−1Γ
[unN − g∆t (Dβ)N η
n+1N − g∆tdβ
(ηn+1N
)0
]+
− ∆tH
1 + γ0∆tdβ
[(unN )0 − g∆t (Dβ)00 (ηnN )0 + g∆t2 (Dβ)00 f
n+1 − g∆tdβηn+1N
]. (3.33)
and rearranging the terms:
ηn+1N
[IΓ − gH∆t2 (Dβ)N I
−1Γ (Dβ)N −
gH
1 + γ0∆t∆t2dβdβ
]=
= ηnN −∆tH (Dβ)N I−1Γ unN −
∆tH
1 + γ0∆tdβ (unN )0
+ gH∆t2(ηn+1N
)0
[(Dβ)N I
−1Γ dβ +
(Dβ)00 dβ
1 + γ0∆t
]. (3.34)
Again, the solution of this linear system at each time step and the replacement of ηn+1N in (3.32a)
brings to the solution of the problem.
3.4 Nonlinear shallow water model: semi–Lagrangian approximation
We now proceed to the discretization of the nonlinear problem, that is, we consider the following
model:
3.4 Nonlinear shallow water model: semi–Lagrangian approximation 51
∂η
∂t+ ∂(Hu)
∂x+ γη = 0 x ∈ R+, t > 0 (3.35a)
∂u
∂t+ u
∂u
∂x+ g
∂η
∂x+ γu = 0 x ∈ R+, t > 0 (3.35b)
∂(Hu)∂x
(0, t) = f(t) (3.35c)
η(x, 0) = η0(x) (3.35d)
u(x, 0) = u0(x) (3.35e)
where we impose Neumann boundary conditions on the fluid discharge at the left endpoint, al-
ready taking into account the damping term in both equations.
As for the continuity equation (3.35a), we discretize it with a θ–method in time and spectral
collocation in space.
On the other hand, we get rid of the advective nonlinearity in the momentum equation (3.35b) by
adopting a semi–Lagrangian approach (for a review see e.g. [SC91]). Before applying the method
to the nonlinear shallow water equations, we review it in the following simple linear advection case:
∂u
∂t+ a
∂u
∂x= 0 x ∈ R, t > 0
u(x, 0) = u0(x) x ∈ R
(3.36)
with a ∈ R being the constant. The solution to (3.36) is:
u(x, t) = u0(x− at), (3.37)
which is constant on the characteristic lines x = at. As before, the numerical solution is denoted
by uni ≈ u(x
(β)i , tn
), where x(β)
i are the SGLR nodes, tn = n∆t are the time levels.
The semi–Lagrangian method exploits the representation (3.37) of the exact solution in terms of
the initial datum. If a > 0 without loss of generality, we have:
u(x
(β)i , tn+1
)= u0
(x
(β)i − a(n+ 1)∆t
)(3.38)
= u0(x(β)i − a∆t− an∆t) = u
(x
(β)i − a∆t, tn
). (3.39)
The value of the solution at time n + 1 in the mesh points is expressed in terms of the solution
values at time n in the points of a mesh that within a single time step will be transported by
the flow onto the computational mesh. Therefore, the definition of the semi–Lagrangian method
is obtained by the following expression:
∂u
∂t+ a
∂u
∂x= Du
Dt≈ un+1 − un∗
∆t,
where we have defined:
52 3 Spectral discretization of boundary–value problems on semi–infinite domains
un∗ (x) := un (x− a∆t) , (3.40)
and used the Lagrangian derivative notation (see (2.1)). The nodal values un∗ in (3.40) are computed
by interpolation. Specifically, for any i = 1, . . . , N , we set:
x∗ := x(β)i − a∆t, ui = u
(x
(β)i
), i = 0, . . . , N (3.41)
α :=
x(β)i − x∗
x(β)i − x(β)
i−1
if uni ≥ 0
x(β)i+1 − x∗
x(β)i+1 − x
(β)i
if uni < 0
Then, the velocity values uni in the point x∗, the foot of the characteristic, can be computed for
all i = 1, . . . , N by linear interpolation as follows:
(un∗ )i =
αuni−1 + (1− α)uni if uni ≥ 0
αuni+1 + (1− α)uni if uni < 0(3.42)
A more accurate value of the velocity (un∗ )i at the foot of the characteristic can be obtained with
cubic interpolation(see e.g. [QSS07]
):
(un∗ )i ≈
i+1∑k=i−2
u(x
(β)k
) i+1∏j=i−2j 6=k
x∗ − x(β)j
x(β)k − x(β)
j
if uni ≥ 0
i+2∑k=i−1
u(x
(β)k
) i+2∏j=i−1j 6=k
x∗ − x(β)j
x(β)k − x(β)
j
if uni < 0
(3.43)
Formulas (3.43) hold for the nodes for which the notations is well–defined, that is, for all i =
4, . . . , N−1 whereas for the most external nodes the linear interpolation strategy (3.42) is adopted.
As for the nonlinear momentum equation (3.35b), an extension of what seen above gives, in the
spectral collocation case:∂u
∂t+ u
∂u
∂x:= Du
Dt≈ un+1 − un∗
∆t,
where now:
un∗ = un (x∗) , (3.44)
x∗ = x− un(x)∆t.
As above, the components of the vector (unN )∗ are computed by linear or cubic interpolation.
The total depth variable H in the mass flux expression (3.35a) is treated explicitly, which amounts
3.4 Nonlinear shallow water model: semi–Lagrangian approximation 53
to a linearization in time. With a semi–Lagrangian discretization of the material derivative in the
momentum equation, the θ–method in time, spectral collocation in space discretization of the
system (3.35a–3.35e) is:
(ηn+1N − ηnN∆t
+ θ∂
∂x
(HnNu
n+1N
)+ (1− θ) ∂
∂x(Hn
NunN ) + γηn+1
N
)(x
(β)j
)= 0 ∀j = 1, . . . , N(
un+1N − (unN )∗
∆t+ gθ
∂ηn+1N
∂x+ g(1− θ)∂η
nN
∂x+ γun+1
N
)(x
(β)j
)= 0 ∀j = 1, . . . , N
∂HnNu
n+1N
∂x(x(β)
0 ) = fn+1
η0N = η0,N ∈ P(β)
N
u0N = u0,N ∈ P(β)
N
In particular, the boundary condition reads collocationwise as:
(dβ)0HnNu
n+1N =
[(Dβ)00
(HnNu
n+1N
)0 + dβ
(H
n
N un+1N
)]= fn+1, (3.45)
where we have used the following notation:
MN 3 Hn
N =
(Hn
N )1
(HnN )2
. . .
(HnN )N
, MN+1 3HnN =
(Hn
N )0 0 · · · 0
0... H
n
N
0
.
These discrete boundary conditions imply as above that:
(1 + γ0∆t)(ηn+1N
)0 = (ηnN )0 −∆t
[θfn+1 − (1− θ)fn
](1 + γ0∆t)
(un+1N
)0 = (unN )0 − g∆tθ
[(Dβ)00
(ηn+1N
)0 + dβηn+1
N
]−g∆t(1− θ)
[(Dβ)00 (ηnN )0 + dβηnN
](3.46)
Therefore, the discrete equations for the internal nodes read:
IΓ un+1N = (unN )∗ − g∆tθ
[(Dβ)N η
n+1N + dβ
(ηn+1N
)0
]− g∆t (1− θ)
[(Dβ)N η
nN + dβ (ηnN )0
]
IΓ ηn+1N = ηnN − θ∆t
[(Dβ)N H
nN u
n+1N + dβ
(HnNu
n+1N
)0
]− (1− θ)∆t
[(Dβ)N H
nN u
nN
+dβ (HnNu
nN )0
](3.47)
54 3 Spectral discretization of boundary–value problems on semi–infinite domains
where (see (3.44)):
(unN )∗ =[(unN )∗]i
i=1,...,N .
Following what done in the linear case, we obtain, for the elevation at the internal nodes:
IΓ ηn+1N = ηnN − θ∆t
(Dβ)N H
nNI−1Γ
(unN )∗ − g∆tθ (Dη)n+1
N − g∆t (1− θ) (Dη)nN
+ dβ(Hn
N )01 + γ0∆t
(unN )0 − g∆tθ (Dη)n+1
0 − g∆t(1− θ) (Dη)n0
− (1− θ)∆t[(Dβ)N H
nN u
nN + dβ (Hn
NunN )0
], (3.48)
where:
(Dη)nN := (Dβ)N ηnN + dβ (ηnN )0 ,
(Dη)n0 := (Dβ)00 (ηnN )0 + dβηnN .
After further reductions, the linear system for the elevation internal nodes reads:
IΓ − θ2g∆t2G
ηn+1N =
IN + θ(1− θ)g∆t2G
ηnN
− θ∆t (Dβ)N HnNI−1Γ (unN )∗ + θ2g∆t2
(ηn+1N
)0 G0 − (1− θ)∆t (Dβ)N H
nN u
nN
−∆t (HnN )0 (unN )0 dβ
(1 + (1− θ)γ0∆t
1 + γ0∆t
)+ θ(1− θ)g∆t2 (ηnN )0 G0, (3.49)
where:
G := (Dβ)N HnNI−1Γ (Dβ)N +
(HnN )0
1 + γ0∆tdβdβ ,
G0 := (Dβ)N HnNI−1Γ dβ +
(HnN )0
1 + γ0∆t(Dβ)00 dβ .
Again, replacing this expression in (3.47) gives the solution of the problem at time n+1. Assuming
a constant average depth H (i.e. assuming a flat bottom profile), the total water depth is updated
according to:
Hn+1N = H + ηnN .
4
Finite–domain approximation and coupling with the
semi–infinite part
The objective of this chapter is to describe a classical finite–volume approach to the numeri-
cal approximation of the shallow water equations in a finite domain and to show how such a
discretization can be coupled to the semi–infinite discretization described in chapter 3, with the ul-
timate goal of using the latter to extend the finite volume discretization to an unbounded domain.
The finite volume method has been chosen as it is a standard approximation technique for the
shallow water and Euler equations(see e.g. [LeV92, Tor01, LG08]
). First, we review the classical
conservative finite–volume numerical methods for the approximation of shallow water equations,
notably Lax–Friedrichs’ method, Godunov’s method and its spatial higher–order extensions. Next,
we detail how to attach to this discretization a spectral discretization in a half line, in order to
employ the obtained setting as an absorbing boundary condition. To this end, the finite volume
discretization described in the first part of the chapter is used as an example, and we remark that
different methods can be employed in the finite domain.
4.1 Finite–volume approximation of the shallow water equations
Following [CDMW98] and [LG08], the monodimensional version of (2.16) reads:
∂U
∂t+ ∂F (U)
∂x= S (U) , (4.1)
where:
U =
hhu
, F (U) =
hu
hu2 + 12gh2
, S (U) =
0
−ghbx
. (4.2)
We recall that U is the vector of the conserved quantities (depth and discharge), F (U) is the flux
function, and S (U) is the source term. The homogeneous version of (4.1)(with S (U) ≡ 0
)is
obtained, under suitable regularity assumptions on the solution, from the following integral form(see [LG08]
):
56 4 Finite–domain approximation and coupling with the semi–infinite part
d
dt
∫ x2
x1
U(x, t)dx = F (U (x1, t))− F (U (x2, t)) , ∀x1, x2 (4.3)
which states that in any interval [x1, x2] the total mass of U(x, t) is modified only by fluxes passing
through the edges.
If the solution is sufficiently smooth, we can consider the quasi–linear form of (4.1):
∂U
∂t+A∂U
∂x= S (U) , (4.4)
where A = F ′ (U) is the Jacobian of the flux. In particular, we have:
A =
0 1
−u2 + gh 2u
. (4.5)
The eigenvalues and eigenvectors of A are given by:
λ1 = u− c, r1 =
1
u− c
,λ2 = u+ c, r2 =
1
u+ c
,where c =
√gh is the celerity. Unless h = 0
(the case of a dry bed, which we do not address
here, see e.g. [Amb95, Tor01]), λ1 and λ2 are real and distinct so that the system (4.4) is strictly
hyperbolic. We also define the Froude number :
Fr := |u|c,
which represents the magnitude of u relative to c.
We discretize the shallow water equations using conservative finite volume methods, starting with
the homogeneous problem. We want to solve (4.3) in the spatial and temporal domain [x1, x2] ×
[0, T ]. We divide the interval [x1, x2] in Nx cells of uniform length ∆x, [xi− 12, xi+ 1
2], i = 1, . . . , Nx,
the centers of the cells being the points xi, i = 1, ..., Nx.
The numerical solution is approximated by cell averages of the solution over grid cells. These
averages are then updated at each time step through the fluxes crossing the cell boundary points
xi− 12, xi+ 1
2, the interfaces.
As for the time discretization, as above we use a uniform mesh tn = n∆t, n = 0, . . . , T/∆t and a
fully explicit method. We define Uni , the numerical approximation to the cell average at time tn
over the i–th cell, by:
4.1 Finite–volume approximation of the shallow water equations 57
Uni ≈
1∆x
∫ xi+ 1
2
xi− 1
2
U (x, tn) dx, (4.6)
and the approximation to the flux at the interface xi− 12by:
F ni− 12≈ 1∆t
∫ tn+1
tnF(U(xi− 1
2, t))
dt. (4.7)
F ni− 1
2is the numerical flux. A finite volume method in conservation form reads
(see e.g. [LeV02]
):
Un+1i = Un
i −∆t
∆x
(F ni+ 1
2− F ni− 1
2
). (4.8)
An example is the finite volume method obtained by Lax Friedrichs numerical flux:
F ni− 12
= 12(F(Uni−1)
+ F (Uni ))− ∆x
2∆t(Uni −Un
i−1). (4.9)
Next, we give a brief review of Godunov’s method. A Riemann problem for the shallow water
equations consists of the hyperbolic system (4.1) together with piecewise constant initial data
with a single jump discontinuity located in x = xd:
U0(x) =
U l if x < xd
U r if x ≥ xd
Since the cell averages are discontinuous over the interfaces xi− 12, i = 2, . . . , Nx and they form a
piecewise constant function in [x1, x2], Godunov’s method can be written in the form (4.8) with
the numerical flux given by:
F ni− 12
= F(U i− 1
2(0)),
where U i− 12(0) denotes the exact solution of the Riemann problem in xi− 1
2with initial data given
by U l = Uni−1, U r = Un
i . Defining the numerical flux in this way, discontinuities propagate
with the correct velocities and no oscillations, and we obtain a stable and monotonic scheme, as
pointed out in [LeV02, CDMW98]. However, solving exactly the Riemann problem at each cell
interface requires the solution of a set of nonlinear equations. Therefore, approximate Riemann
solvers have been constructed, and Roe solver [Roe81] is one of them. Instead of (4.1), we solve,
in the homogeneous case∂U
∂t+ Ai− 1
2
∂U
∂x= 0, (4.10)
where U denotes an approximation to the exact solution of the Riemann problem and Ai− 12an
approximation to A as in (4.4) in a neighbourhood of U i,U i−1, which must have the following
properties:
1. Ai− 12is diagonalizable with real eigenvalues (hyperbolicity);
2. Ai− 12−→ F ′(U) as U i−1,U i −→ U (consistency);
3. Ai− 12(U i −U i−1) = F (U i)− F (U i−1).
58 4 Finite–domain approximation and coupling with the semi–infinite part
The last relation is indeed consistent with the Rankine–Hugoniot condition for the speed s of a
shock wave connecting U i and U i−1, which is:
F (U i)− F (U i−1) = s(U i −U i−1). (4.11)
In particular, Roe solver sets (see [CDMW98, Tor01]):
U =h− 1
2i−1U i−1 + h
− 12
i U i
h− 1
2i−1 + h
− 12
i
,
so that:
h = hi−1 + hi2
,
u =√hi−1ui−1 +
√hiui√
hi−1 +√hi
.
Roe matrix in our case reads (see also [LeV02]):
Ai− 12
=
0 1
−u2 + gh 2u
, (4.12)
and satisfies the properties seen above. The Roe matrix is simply the matrixA as in (4.5) computed
in the state U . The eigenvalues and eigenvectors are:
λ1 = u− c, r1 =
1
u− c
,λ2 = u+ c, r2 =
1
u+ c
,where c =
√gh. Across the discontinuity at the interface we have:
U i −U i−1 =2∑p=1
αpi− 1
2rp,
which gives, setting δ ≡ U i −U i−1([LeV02]
),
α1i− 1
2= (u+ c) δ1 − δ2
2c, α2
i− 12
= − (u− c) δ1 + δ2
2c.
Then, Roe numerical flux is:
F ni− 12
= F(Uni−1)
+2∑p=1
λ−p αp
i− 12rp, (4.13)
where λ−p = min(λp, 0
). Since the time advancing is performed with an explicit procedure, a CFL
condition must hold:
C = ∆tλmax∆x
≤ 1,
where λmax is the maximum eigenvalue for the cell (for a definition of the Courant number see
below, (5.3)). Godunov’s method with Roe’s approximate Riemann solver is first–order accurate
in space and time, conservative and monotonic.
4.1 Finite–volume approximation of the shallow water equations 59
Entropy fix
Roe linearization gives the correct speeds of propagation for shock waves. As for rarefaction waves,
they are replaced by Roe solver with discontinuities. Let λp,i, λp,i−1(p = 1, 2) be the values of the p–
th eigenvalue computed to the right and to the the left of the interface placed in xi− 12, respectively;
an entropy fix (i.e. a modification of the speed made to obtain an entropic solution) is needed only
in the transonic case, that is when λp,i−1 < 0, λp,i > 0. Omitting the details(see [LeV92]
), the
modified numerical flux reads:
F ni− 12
= F(Uni−1)
+2∑p=1
λ∗pαp
i− 12rp,
where:
λ∗p := λ−p,i−1
(λ+p,i − λp
λ+p,i − λ
−p,i
),
and the subscript i means that the eigenvalue is computed from the averages in the i–th cell.
Higher–order extension
We extend Godunov’s method to second order spatial accuracy using slope limiters. In this ap-
proach(see also [LeV92, CDMW98, LeV02]
), instead of piecewise constant cell averages, we con-
struct a piecewise linear representation of U = (h, hu)T :
U (x, tn) = Uni + σni (x− xi) , x ∈ [xi− 1
2, xi+ 1
2]
and then use this reconstruction in Godunov’s scheme. In general, though, unless the slopes σniare properly limited, this procedure results in a growth of total variation of the numerical solution
over the subsequent time steps; therefore, in order to have that:
TV(U (·, tn)
)≤ TV (Un
i ) ,
so that the scheme is globally TVD, one of the following options has been adopted:
1. MINMOD: σni = 1∆x
minmod(Uni −Un
i−1,Uni+1 −Un
i
), where the minmod function is:
minmod(a, b) := 12
(sgn(a) + sgn(b)) min (|a|, |b|) ;
and has been applied componentwise to U .
2. SUPERBEE: σni = 1∆x
maxmod(σ
(1)i , σ
(2)i
), where:
σ(1)i = minmod
(Uni+1 −Un
i , 2(Uni −Un
i−1))
σ(2)i = minmod
(2(Uni+1 −Un
i
),Un
i −Uni−1).
60 4 Finite–domain approximation and coupling with the semi–infinite part
The resulting numerical scheme reads(for details see [LeV02]
):
Un+1i = Un
i −∆t
∆x
(F ni+ 1
2− F ni− 1
2
)− ∆t
∆x
(F i+ 1
2− F i− 1
2
)(4.14)
where(see also [LeV92]
):
F i− 12
= 12
2∑p=1
∣∣∣λpi− 12
∣∣∣ (1− ∆t
∆x
∣∣∣λpi− 12
∣∣∣)(αpi− 12rpi− 1
2
)LTD
, (4.15)
where the last factor represents the limited version of the p–th wave arising in the (Roe approxi-
mated) solution to the Riemann problem at the interface. This limited version is obtained applying
one of the slope limiters to the jump αi− 12rpi− 1
2and the jump αI− 1
2rpI− 1
2in the upwind direction:
I =
i− 1 if λpi− 1
2> 0
i+ 1 if λpi− 1
2< 0
Note that the speeds λpi− 1
2are the ones corrected with the entropy fix. The resulting scheme is
first–order accurate near discontinuities, but is second–order accurate in regions of smooth flow.
4.2 Coupling the finite and semi–infinite discretizations
Having shown a standard discretization for the shallow water equations in a finite domain, in
order to use the spectral semi–infinite discretization described in chapter 3 as an absorbing layer
to extend the finite domain discretization, we need a strategy to couple the two domains, that is,
to connect the right endpoint of the finite domain to the left endpoint of the semi–infinite one, the
contact point being the interface. A representation of the coupling setting is given in Figure 4.1.
x 32
x 12
xNx+ 12
x2x1 xNx
x(β)0 x
(β)1 x
(β)2 x
(β)3 x
(β)4
Fig. 4.1: Nodes (circles) and cell interfaces (ticks) of the finite interval discretization, along with the
nodes (squares) of the semi–infinite part.
Multidomain discretizations have been extensively studied in the context of Partial Differential
Equations, and in the hyperbolic case the constraints at the interface involve suitable conditions
4.2 Coupling the finite and semi–infinite discretizations 61
on the fluxes or the characteristic variables. For a detailed discussion on domain decomposition
methods in the context of spectral approximations to hyperbolic problems we refer to [Qua90,
CHQZ07, YCC10, QV99].
In our case, we found that a simple condition on continuity of the mass flux is sufficient for the
interface to be transparent as waves pass through it. On one hand, this proves to be an effective
way to attach our semi–infinite discretization with a flux–based finite volume scheme; on the
other hand, this approach enables to use our semi–implicit method in the semi–infinite domain
and does not force us to use a characteristic–based numerical scheme, as we should have done
using a classical characteristic coupling. Specifically, in sections 3.3–3.4 the spectral collocation
formulation of the equations in the semi–infinite part with Neumann condition on the discharge
at the left endpoint has been considered. Therefore, in the coupling process one has simply to put
into the expression of the boundary condition (3.45):
(dβ)0HnNu
n+1N =
[(Dβ)00
(HnNu
n+1N
)0 + dβ
(H
n
N un+1N
)]= fn+1
as fn+1 the datum coming from the left domain, which is an approximation of the derivative of the
discharge in the point xNx+1/2. For simplicity, to obtain this approximation we adopt a first–order
one–sided finite difference approach, i.e. we impose:
fn+1 = (hu) (xNx)− (hu) (xNx−1)∆x
, (4.16)
where ∆x is the finite domain uniform grid spacing. The quantities on the right hand side of (4.16)
result from time advancing the finite volume scheme according to (4.8).
Dually, the finite volume scheme’s numerical flux at the last interface of the left domain xNx+1/2
is computed taking as right states the depth and discharge values at the left endpoint x(β)0 of the
the semi–infinite domain.
Therefore, the solution of the coupled scheme entails the execution of the following algorithm:
Algorithm 1. For each tn, n = 0, . . . , T/∆t− 1:
1. Advance the finite volume scheme by (4.8), having computed FnNx+ 1
2assuming a right state:
UnNx+1 =
(h(x
(β)0
), h(x
(β)0
)u(x
(β)0
)).
2. Compute:
fn+1 = (hu) (xNx)− (hu) (xNx−1)∆x
,
and replace it in:
(dβ)0HnNu
n+1N =
[(Dβ)00
(HnNu
n+1N
)0 + dβ
(H
n
N un+1N
)]= fn+1.
3. Advance the spectral collocation scheme by solving system (3.47).
5
Numerical Results
In this chapter we present the results obtained with the implementation of the numerical
schemes previously described. Algorithms have been coded in MATLAB R©1 and different tests
have been simulated.2 First, we deal with the spectral discretization, and demonstrate the reliabil-
ity of Scaled Laguerre functions by showing that different choices of the number of nodes and the
scaling parameter β enable to represent accurately different spatial scales and to obtain different
grid spacing. Then, we test the accuracy of a stand–alone spectral collocation approximation of
shallow water equations described in sections 3.3–3.4. Next, we validate the coupling approach
described in section 4.2 by matching a semi–infinite spectral discretization with a standard finite
volume approximation of the shallow water equations and compare the results with those of a
full finite volume approximation, showing that the interface is sufficiently transparent to the the
propagation of waves. These results will give motivation for the use of the spectral semi–infinite
discretization as an absorbing layer for the approximation of open boundary conditions for waves
departing from the finite domain. In particular, we show that a reasonable number of spectral base
functions to the right of the artificial interface is sufficient to damp the waves without reflections
either at the interface or at infinity.
We give here some definitions that will be used throughout the chapter. Given two functions uNand u defined on Ω, we define the absolute and relative Lp error of uN with respect to u as:
Eabsp (uN ) := ‖u− uN‖p ∀ 1 ≤ p ≤ ∞ (5.1)
Erelp (uN ) :=‖u− uN‖p‖u‖p
∀ 1 ≤ p ≤ ∞ (5.2)
where we recall the definition of the Lp norms:
‖u‖p :=(∫
Ω
|u|p dx)1/p
∀ 1 ≤ p <∞
‖u‖∞ := supx∈Ω |u(x)|
1 Vers. 7.9.0.529 (R2009b)2 The simulations have been performed on a personal computer with processor IntelR©CoreTM2 Duo CPUP9500 @ 2.53 GHz–2.53 GHz, 3072 MB RAM
64 5 Numerical Results
Dealing with discrete functions, we also recall the p–norm definition for a vector u = ujNj=1:
‖u‖p :=
∑j
|uj |p1/p
∀ 1 ≤ p <∞
‖u‖∞ := maxj (|uj |)
Given a matrix A ∈Mn(R), we define its condition number in the 2–norm as (see [QSS07]):
K2(A) = ‖A‖2∥∥A−1∥∥
2 ,
where:
‖A‖2 = supx6=0
‖Ax‖2
‖x‖2.
Next, we define the local Courant number as:
Ci = |uni |∆t∆x
,
where ∆t is the time step, ∆xi is the distance between the (i − 1)–th and the i–th spatial grid
point, uni is the total velocity of the numerical scheme within the cell. Hence, referring to the
advective formulation (2.19) of the shallow water equations, we define the global Courant number
for a discretization Ij of a spatial interval Ω as:
C := ∆tmaxj
[∣∣unj + cnj∣∣
∆xj
], (5.3)
where cnj :=√gHn
j .
5.1 Scaling properties of scaled Laguerre functions
Introduced in chapter 1 (see (1.24)), Scaled Laguerre functions (SLFs) are an orthonormal basis
of L2(R+), and can therefore be used to approximate functions and integrals on the semi–infinite
real half line. Table 5.1 shows the position of the first two internal nodes for unscaled (β = 1)
GLR quadrature formulas, as well as the position of the last node, i.e. the actual domain length,
for different polynomial orders. In particular, the size of the domain increases with the order of
the Laguerre function. Taking 1 m as the unity, for N = 4 the last node is located 10.95 m away
from the origin, while with N = 125 the domain length approaches 500 m.
5.1 Scaling properties of scaled Laguerre functions 65
N 4 6 8 10 12 15
0.75, 2.57 0.53, 1.8 0.41, 1.39 0.33, 1.13 0.28, 0.95 0.23, 0.77
10.95 17.65 24.59 31.68 38.89 49.84
N 20 30 50 75 100 125
0.17, 0.59 0.12, 0.4 0.07, 0.24 0.05, 0.16 0.036, 0.122 0.029, 0.098
68.38 106.05 182.62 279.45 376.94 474.81
Table 5.1: Position of the first two internal nodes x(β)1 –x(β)
2 and the last node x(β)N , for unscaled GLR
quadrature formulas.
The position of the nodes for SGLR quadrature formulas (i.e. for β 6= 1) is obtained multiplying
the values in the previous table by 1β
= L, which can be seen as a measure of the domain length.
Indeed, see expressions (1.20–1.21–1.22–1.23a), the value of the n–th scaled Laguerre polynomial
L(β)n or function L
(β)n evaluated in the point x is equal to the value of its unscaled version Ln
evaluated in the point βx. For example, in table 5.2 we see the values for β = 0.0025, i.e. L = 400.
N 4 6 8 10 12 15
2.97E02,1.03E03 2.11E02,7.19E02 1.64E02,5.54E02 1.34E02,4.51E02 1.13E02,3.81E02 9.18E01,3.09E02
4.38E03 7.06E03 9.83E03 1.27E04 1.56E04 1.99E04
N 20 30 50 75 100 125
7E01,2.35E02 4.74E01,1.59E02 2.88E01,9.65E01 1.93E01,6.48E01 1.45E01,4.87E01 1.17E01,3.91E01
2.74E04 4.24E04 7.3E04 1.12E05 1.51E05 1.9E05
Table 5.2: Position of the first two internal nodes x(β)1 –x(β)
2 and the last node x(β)N , for scaled GLR
quadrature formulas, β = 0.0025.
For this value of β, taking again 1 m as the unity, we see that for instance in the case N = 100
the first node is 15 m, the second one 49 m away from the origin, whereas the last node is located
151 km away. Keeping β fixed and taking an increasing number of nodes refines the spectral
grid spacing, extending the range of the computational domain as the last node is located further
and further away from the origin. Moreover, keeping N fixed and reducing the value of β (i.e.
increasing L) the range of the computational domain is still extended, but this time the spectral
grid is coarsened, as pointed out by figure 5.1, which shows the Scaled Laguerre function L(β)
5
for four different values of β.
66 5 Numerical Results
Fig. 5.1: The function L (β)5 for
four different values of β.
50 100 150
-0.4
-0.2
0.2
0.4
0.6
0.8
1.0
Β 0.001
Β 0.01
Β 0.1
Β 1
5.2 Unsteady advection–diffusion–reaction equation
We now present some numerical tests which validate the θ–method/G–NI approximation to
parabolic advection–diffusion–reaction equation described in 3.2.
Non–homogeneous problem
We present the numerical results obtained with the implementation of θ–method, starting from
the non–homogeneous problem. In this case the forcing term is chosen in such a way that the exact
solution is u(x) = xe−x sin(x− t)2, in order to have unitary wave velocity and convergence to zero
at +∞ as required by a G–NI spatial approximation with Laguerre functions. Equation (3.5) is
solved in the spatial–temporal domain (0,+∞)× (0, T ). In figure 5.2 we see the plot of the exact
solution at the initial (left) and final (right) time, along with the computed solution obtained with
the θ–method at the final time T = 10 s. We observe a good agreement with few time steps.
Figure 5.3 shows the relative L2(R+) error decay (computed at time T = 10 s) as a function of the
time steps (for the definition see below, (5.2)). As expected, the Crank–Nicolson method produces
a solution much more accurate than the one computed by Backward Euler method.
Figure 5.4 shows the relative L2(R+) error decay (computed at time T = 10 s) as a function of the
number of nodes used for the numerical integration in G–NI approximation. The level attained
by the error corresponds to the temporal part in the error estimate, which bounds the accuracy
of the scheme. In order to improve our results, one could employ higher order methods for the
time discretizations, which we do not consider here. As spectral formulations are characterized
by full matrices, a major concern about their stability is the matrices’ condition numbers (see
e.g. [CHQZ06, QSS07, QV08]). Table 5.3 reports the condition numbers of the matrix K of G–NI
5.2 Unsteady advection–diffusion–reaction equation 67
0 2 4 6 8 100
0.05
0.1
0.15
0.2
0.25
0.3
0.35
x
Solution at initial time
0 2 4 6 8 100
0.05
0.1
0.15
0.2
0.25
0.3
0.35
x
Solution at time T
Fig. 5.2: Exact solution (solid line) at time t0 = 0 and T = 10. Solution computed with the θ–method
(right, dashed, with N = 64, ∆t = 0.25, θ = 1).
101 102 103 10410−8
10−7
10−6
10−5
10−4
10−3
10−2
10−1Temporal accuracy
1/∆t
Rel
ativ
e L2 E
rror
θ = 1θ = 0.5
Fig. 5.3: θ–method’s relative L2 error
as a function of the num-
ber T/∆t of time steps with
N = 128 for the parabolic
problem (3.5) with µ = 10,
b = 1, γ = 4.
formulation (3.12) for T/∆t = 1000 and different choice of the number of nodes. The values are
all quite low, and the observed growth is sublinear.
Homogeneous problem
Next we consider the case f = 0 and test that the numerical scheme we implemented gives the
awaited results. First of all, we verify that the scheme simulates the correct wave velocity.
68 5 Numerical Results
Fig. 5.4: θ–method’s relative L2 er-
ror as a function of N with
T = 10 s, ∆t = 0.01 s for
the parabolic problems (3.5)
with µ = 10, b = 1, γ = 4. 0 20 40 60 80 100 120 140 160 18010−6
10−5
10−4
10−3
10−2
10−1
100
101 Spatial accuracy
N
Rel
ativ
e L2 E
rror
θ=1θ=0.5
N 4 8 16 32 64 128 180
θ1 4.56 8.43 14.99 27.17 49.06 88.25 117.66
0.5 4.94 10.00 17.78 32.09 57.99 104.37 139.17
Table 5.3: Condition numbers of the matrix K for non–homogeneous advection–diffusion–reaction prob-
lem.
Results in figure 5.5 are good in this direction: using the meters as spatial, seconds as temporal
units of measurement, letting the scheme advance for T = 120 s with a coefficient b = 0.1 (0.1 m/s)
we find that the initial impulse (solid line), given as above by the function u0(x) = xe−x sin2(x),
has covered 12 meters and has been damped by diffusion µ = 0.01 (the reaction coefficient γ is
also equal to 0.01). Numerical solution at time T is the dashed–dotted line. Condition numbers of
the matrix K of G–NI formulation (3.12) are reported in table 5.4 for different number of spectral
nodes. In both cases, a slightly superlinear growth is observed. However, the values do not reach
worrying levels.
N 4 8 16 32 64 128 180
θ1 5.29 12.1 24.96 44.21 79.00 142.30 189.86
0.5 5.36 12.67 28.39 54.30 95.58 172.11 229.71
Table 5.4: Condition numbers of the matrix K for undamped homogeneous advection–diffusion–reaction
problem.
5.2 Unsteady advection–diffusion–reaction equation 69
0 5 10 15 20 25 30 35 40 45−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
x
Fig. 5.5: Evolution of the numerical
solution for the parabolic
problem (3.5) with f =
0, µ = 0.01, b = 0.1,
γ = 0.01 computed with θ–
method with T = 120 s,
∆t = 1.2 s, N = 64.
Finally, we take into account the homogeneous case with a variable reaction coefficient, acting as
a wave damper; in particular the function we use to model the damping is the following:
γ(x) = γD(∆γ,α,w, L0, x) := ∆γ
1 + exp(αL0 − x
σ
) . (5.4)
The above expression is a sigmoid connecting the level γ = 0 with γ = ∆γ, where ∆γ is the value
of the absorption that we want to reach at the end of the domain. The sigmoid is centered in the
point x = αL0 and the transition area between the two values is w = σ log(99) wide. In general
L0 is the width of the semi–infinite part, and thus depends on L = 1/β, the scale coefficient of
the SGLR quadrature formulas. Besides, adjusting the parameters α ∈ [0, 1], the width w and the
strength ∆γ we are able to model different damping terms, see figure 5.6, where, taking 1m as
the unit, the x coordinate spans a interval of 40km. These lengths are suitable for a geophysical
interpretation. In this test we consider the following coefficient:
0 0.5 1 1.5 2 2.5 3 3.5 4
x 104
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
x
γ(x)
0 0.5 1 1.5 2 2.5 3 3.5 4
x 104
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
x
γ(x)
Fig. 5.6: The functions γD(2, 0.5, 5000, 40000, x) and γD(2, 0.8, 2000, 40000, x).
70 5 Numerical Results
0 5 10 15 20 25 30 35 40 45−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35T=100
x0 5 10 15 20 25 30 35 40 45
−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35T=200
x
0 5 10 15 20 25 30 35 40 45−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35T=250
x0 5 10 15 20 25 30 35 40 45
−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35T=300
x
Fig. 5.7: Evolution of the initial value u(x) = xe−x sin(x)2 with the θ–method. In solid line the initial
value, in dashed line the solution with γ(x) = 0, in dashed–dot line the solution with γ(x) =
γD(1, 0.5, 15, 60, x).
γ(x) = γD(1, 0.5, 15, 60, x).
and show the results of the evolution of the computed solution. In figure 5.7 we see the result of
the time evolution of the initial value u0(x) = xe−x sin(x)2 with the θ–method (θ = 1, ∆t = 0.01,
N = 64, µ(x) = 0.01, b(x) = 0.1), for different values of T , with reaction coefficient fixed to the
value 0 (dashed line) or equal to a sigmoid centered in x = 30 m. We observe that far from the
middle of the sigmoid the two computed solutions are practically undistinguishable, but as time
goes on, the scheme with the reaction coefficient damps more and more the initial value, and in
x = 30 m (i.e. T ≈ 300 s) it is virtually null, whereas the damping in the scheme without reaction
is only due to the diffusion coefficient.
5.3 Validation of the spectral collocation method 71
5.3 Validation of the spectral collocation method
The objective of these first set of tests is to validate the use of Laguerre–type base functions, and
specifically scaled Laguerre functions, to approximate hyperbolic initial boundary–value problems.
To this end, as a first step towards the fully nonlinear discretizations, we show the results obtained
with the spectral collocation discretization of the linear shallow water equations in a semi–infinite
domain, referring for the theoretical discussion to section 3.3. In particular, an implicit Euler time
discretization will be adopted.
As for the domain data, setting β = 10−3 and N = 100, we have x(β)N = 3.77 · 105 m so that,
taking 1 m as the spatial unity, the domain is 377 km wide. We consider two boundary conditions
on the left endpoint x = 0:
1. Homogeneous Neumann boundary condition on the discharge, i.e. expression (3.28) with f = 0.
2. Transmissive boundary conditions (3.25)
As initial condition on the free surface elevation, we consider a bell–shaped perturbation:
η0(x) = η exp
[−(x− x0
σ
)2]
(5.5)
that is, a Gaussian hump of amplitude η centered in x = x0 with a standard deviation σ. The initial
velocity will be set to zero. Specifically, we center the hump at x0 = 80000 m, with σ = 20000.
We solve linear shallow water equations with the following data:
T = 15000 s, ∆t = 6 s, H = 10 m, u = 0 m/s, η = 0.3 m
In particular, as η H, the use of the linear model (see also (2.22)) is justified. The corresponding
global Courant number is 1.635. Figure 5.8 shows the initial depth, figure 5.9 its evolution with
the first boundary conditions. The hump splits into two crests propagating at the correct velocity
c =√gH. When the left crest reaches the left endpoint, it is reflected and begins to travel
rightwards.
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.9
9.95
10
10.05
10.1
10.15
10.2
10.25
10.3
10.35
10.4
x
H(x
)
Initial depth
Fig. 5.8: Initial depth for validation of
spectral–collocation scheme.
72 5 Numerical Results
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x)
t = 1500
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 3000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 3750
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 6000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 7500
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = T
Fig. 5.9: Evolution of an initial Gaussian depth perturbation with implicit Euler–spectral collocation
scheme, C ≈ 1.64, homogeneous Neumann boundary conditions.
Figure 5.10 shows the evolution with the transmissive boundary conditions. Now the left crest
leaves the domain as if the left boundary were not there, while the right crest keeps propagating
rightwards at the correct velocity. Table 5.5 shows the condition numbers of the matrix F as of
(3.24) (with Γ ≡ 0) for different choices of the number of nodes. A sublinear growth is observed
both in the case of Neumann and in the case of non–reflecting boundary conditions and the values
are acceptable.
5.3 Validation of the spectral collocation method 73
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4t = 1500
x
H(x
)
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 3000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 3750
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 6000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x)
t = 7500
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x)
t = T
Fig. 5.10: Evolution of an initial Gaussian depth perturbation with implicit Euler–spectral collocation
scheme, C ≈ 1.64, transmissive boundary conditions.
N 4 8 16 32 64 128 180
BC1 1.07 1.25 1.92 4.50 14.61 80.99 106.67
2 1.005 1.017 1.063 1.26 2.35 11.79 34.4
Table 5.5: Condition numbers of the matrix F for spectral–collocation approximation of linear shallow
water equations, homogeneous Neumann (1) and non–reflecting (2) boundary conditions.
74 5 Numerical Results
Fig. 5.11: Damping coefficient for spectral–
collocation scheme validation. 0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5x 10−3
x
γ D(x
)
Damping coefficient
Finally, in view of the use of the semi–infinite discretization as an absorbing layer, we consider
the linear shallow water equations with the damping coefficient (see (3.15)), in particular with
homogeneous Neumann boundary conditions and damping given by the function (see above, (5.4)
in 5.2):
γ(x) = γD(∆γ,α,w, L0, x).
Here we take:
L0 = 4 · 105 m, α = 0.7, ∆γ = 0.005, w = L0/2,
so that the damper is centered in x = 280000 m (see figure 5.11). Moreover, now we set T =
40000 s, ∆t = 16 s, so that C ≈ 4.88. Figure 5.12 shows the results of the evolution of the initial
hump; as in figure 5.9, the left crest is reflected at the left endpoint, but now the right one begins
to be damped as expected until all the perturbations have been erased by the damper, as can be
seen by the plot of final time depth.
5.3 Validation of the spectral collocation method 75
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
x
H(x
)
t = 4000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 5000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x)
t = 10000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 16000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = 20000
0 0.5 1 1.5 2 2.5 3 3.5 4
x 105
9.5
9.6
9.7
9.8
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
t = T
Fig. 5.12: Evolution of an initial Gaussian depth perturbation with implicit Euler–spectral collocation
scheme with damping, C ≈ 4.88, homogeneous Neumann boundary conditions.
Therefore, in contrast with [Fra72], these tests show that orthogonal functions in the real positive
half line can be used to represent solutions to hyperbolic problems in environmental modelling.
76 5 Numerical Results
5.4 Validation of the coupling approach
To test the accuracy of the coupling strategy described in the previous chapter, we compare the
coupled finite/semi–infinite fully nonlinear discretization with a stand–alone finite volume approx-
imation. To this end, we consider again the initial condition (5.5) on free surface elevation for the
shallow water model, the initial velocity being again set to zero. In particular, we consider two
configurations: a quasi–linear case, in which the initial perturbation is such that η H, and a
nonlinear case. For both of them, we compare the results of a full finite volume approximation
with the ones of the split finite volume/spectral collocation approximation, considering the full
finite volume scheme as the reference solution. We recall that the finite volume approximation is
performed on the conservative formulation of the shallow water equations (see above, Eqs. (4.1)
and (4.8)) with Lax–Friedrichs, Godunov or minmod slope–limited spatial second order numerical
fluxes, see (4.9–4.13–4.15). As for the spectral part, the spectral collocation scheme described in
chapter 3 is used for space discretization. We recall that the nonlinearity in the continuity equation
is handled by linearizing in time the total depth variable in the mass flux, whereas the advective
nonlinearity is handled with a semi–Lagrangian approach with a cubic interpolation technique (see
section 3.4) for the reconstruction of the velocity at the foot of the characteristic lines.
In this setting, the aim is to show that no spurious phenomena occur as a wave travelling in either
direction crosses the interface. Indeed, we check that in the finite domain the relative error of the
coupled discretization with respect to the full finite volume one is small enough to validate the
coupling approach. The lack of spurious phenomena arising from the interface will give a motiva-
tion for the use of the spectral semi–infinite discretization as a possible alternative to standard
absorbing layers.
5.4.1 Linear wave
First we consider the case of an initial perturbation of the form (5.5), with:
η = 0.3 m, H0N ≡ 10 m
so that we are in a case of linear wave motion. We consider the coupled approximation described
in section 4.2, with the following data:
x 12
= 0, xNx+ 12
= x(β)0 = 10000 m, Nx = 1250, N = 180, x(β)
N = 2.86 · 105,
where we recall that x 12is the left endpoint of the finite domain, xNx+ 1
2= x
(β)0 is the interface
between the right endpoint of the finite domain and the left endpoint of the semi–infinite one,
Nx is the number of cells of the finite volume discretization of the left domain, N is the number
5.4 Validation of the coupling approach 77
of nodes of the spectral semi–infinite collocation approximation, x(β)N is the last spectral node. In
this and the following tests, the parameters have to be set up so that the semi–infinite domain
grid spacing is of the same order of the finite volume grid spacing at the interface, that is:
xNx+ 12− xNx− 1
2≈ x(β)
1 − x(β)0 .
Indeed, with several tests we have checked that if the above condition is violated, spurious phe-
nomena are more likely to occur at the interface, corrupting the solution in the finite domain.
As for the domain data, taking 1 m as the unity, the finite domain is 10 km wide, whereas the
semi–infinite part spans 276 km, for a total length domain of 286 km. A high number of nodes is
employed in order to validate the coupling strategy by sending a wave through the interface.
Test 1
We center the initial perturbation at x0 = 12000 m, to the right of the interface. Figure 5.13 shows
the initial depth for two different values of σ.
0 0.5 1 1.5 2
x 104
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
Initial total depth
σ=500σ=1000
0 0.5 1 1.5 2
x 104
9.9
10
10.1
10.2
10.3
10.4
x
H(x
)
Initial total depth
σ=500σ=1000
Fig. 5.13: Initial depth perturbations for linear wave test 1 (left) and 2 (right).
78 5 Numerical Results
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 104
9.95
10
10.05
10.1
10.15
Depth at time T
x
H(x
)
FV−SC
FV
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 104
−0.15
−0.1
−0.05
0
0.05
0.1
0.15
Velocity at time T
x
u(x)
FV−SC
FV
Fig. 5.14: Comparison of coupled (solid black) and full finite volume (dashed red) schemes’ computed
depth (left) and velocity (right), linear wave.
We apply algorithm 1 of section 4.2 until the final time T = 400 s. The initial hump splits into
two waves travelling in opposite directions. The left crest crosses the interface and, propagating
at the correct velocity√gH ≈ 10 m/s, reaches at t = T the point x ≈ 8000 m, whereas the right
crest symmetrically travels the same length in the positive direction.
The coupled configuration is compared with a full finite volume approximation in a domain 20 km
wide with the same initial data. Grid spacing has been suitably chosen so that the first 10 km are
covered by the same number of nodes of the finite part of the coupled scheme, i.e. Nx = 1250 in the
first 10 km. In view of the use of our coupled discretization to simulate open boundary conditions,
we want reflections into the finite domain generated at the interface to be minimized. Table 5.6
shows, for σ = 500 and C ≈ 0.6 (T/∆t = 825) , the L1, L2 and L∞ errors on free surface elevation
and velocity in the first 10 km computed by the coupled model relative to the same quantities
computed by the full finite volume discretization, taken as reference. We consider different choices
(Lax–Friedrichs, Godunov, minmod slope limited) of the numerical flux for the finite volume
approximations, with fully implicit (θ = 1) and semi–implicit (θ = 0.6) time discretizations of the
coupled scheme.
Results for a longer wave (σ = 1000) are reported in table 5.7, while figure 5.14 shows comparison
of computed final depth and velocity at the final time by the coupled and full finite volume schemes.
For this example, we considered C ≈ 0.9, θ = 1, σ = 1000 and Lax–Friedrichs flux. The solution
computed by the coupled formulation is practically indistinguishable from the one computed by
the reference method. Results for C ≈ 0.9 (T/∆t = 550),and C ≈ 0.3 (T/∆t = 1650) are reported
in the appendix (tables A.1–A.2 for σ = 500, tables A.3–A.4 for σ = 1000). Relative errors are
5.4 Validation of the coupling approach 79
higher for the shorter wave and for the velocity variable. Nonetheless, on the whole they are all
acceptable in an environmental modelling context (cfr. e.g. [RCB05]).
Table 5.6: Test 1 results, C ≈ 0.6, σ = 500, linear wave.
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 3.38E-05 2.34E-05 5.2E-05 6.8E-05 3.01E-05 2.2E-05
u 0.025 0.018 0.039 0.051 0.023 0.017
Erel2η 7.5E-05 5.62E-05 1.24E-04 1.54E-04 6.7E-05 5.18E-05
u 0.021 0.015 0.033 0.043 0.018 0.014
Erel∞η 3.19E-04 2.2E-04 5.62E-04 6.67E-04 2.77E-04 1.88E-04
u 0.024 0.016 0.038 0.051 0.02 0.013
Table 5.7: Test 1 results, C ≈ 0.6, σ = 1000, linear wave.
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 1.86E-05 1.55E-05 2.33E-05 3.47E-05 1.41E-05 8.1E-06
u 0.007 0.006 0.009 0.013 0.005 0.003
Erel2η 3.46E-05 2.77E-05 4.72E-05 6.58E-04 2.69E-05 1.46E-05
u 0.007 0.005 0.009 0.013 0.005 0.003
Erel∞η 1.07E-04 6.59E-04 1.43E-04 2.07E-04 7.97E-05 3.89E-05
u 0.008 0.005 0.01 0.015 0.006 0.003
Test 2
In this second test we center the initial perturbation η0, with η = 0.3 m, at x0 = 5000 m, i.e.
in the finite domain, see figure 5.13. We advance the coupled and the full finite volume scheme
taking as final time T = 1000 s. At this time, as c =√gH ≈ 10 m/s all the waves are out of the
finite domain (the right crest reaching x ≈ 15000 m, the left one x = −5000 m) and the elevation
and velocity variables are expected to be null within the finite domain. The other parameters are
the same as in the previous case. This test aims to check that no spurious waves are generated at
the interface and bounce back into the finite domain. Absolute L1, L2 and L∞ errors for C ≈ 0.3
80 5 Numerical Results
(T/∆t = 4100) are reported in table 5.8 for σ = 500, in table 5.9 for σ = 1000. L1 and L2 errors
have been rescaled by the domain length and its square root, respectively, to obtain average values.
Results for C ≈ 0.9 (T/∆t = 1375) and C ≈ 0.6 (T/∆t = 2050) are reported in the appendix
(tables A.5–A.6 for σ = 500, A.7–A.8 for σ = 1000). In this case, too, the errors are all relatively
low. Furthermore, we report in table 5.10 the condition numbers of the matrix F = IΓ − θ2g∆t2G
(see (3.49)) of the semi–infinite spectral collocation system computed at final time for θ = 1,
σ = 1000. The obtained values are very low.
Table 5.8: Test 2 results, C ≈ 0.3, σ = 500, linear wave.
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 1.7E-06 1.41E-05 1.95E-05 2.2E-05 2E-05 2.8E-05
u 1.7E-06 1.4E-05 1.9E-05 3.7E-06 2E-05 2.8E-05
Eabs2η 1.3E-05 1E-05 1.2E-05 1.8E-05 1.4E-05 1.73E-05
u 1.3E-05 1E-05 1.2E-05 1.4E-05 9E-06 1.7E-05
Eabs∞η 2.25E-04 0.0015 0.0015 3.21E-04 0.002 0.0021
u 2.23E-04 0.0014 0.0015 3.18E-04 0.002 0.0021
Table 5.9: Test 2 results, C ≈ 0.3, σ = 1000, linear wave.
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 6.1E-07 4.6E-07 5.2E-07 6.1E-07 5.6E-07 6.4E-07
u 6.1E-07 4.6E-07 5.1E-07 6E-07 5.6E-07 6.4E-07
Eabs2η 2.43E-06 2.93E-06 3.21E-06 2.52E-06 3.86E-06 4.27E-06
u 2.4E-06 2.9E-06 3.18E-06 2.49E-06 3.83E-06 4.23E-06
Eabs∞η 2.65E-05 5.68E-05 5.69E-05 3.24E-05 7.58E-05 7.61E-05
u 2.62E-05 5.63E-05 5.63E-05 3.21E-05 7.51E-05 7.54E-05
N 4 8 16 32 64 128 180
1.0008 1.0026 1.0095 1.036 1.14 1.553 2.09
Table 5.10: Condition numbers of the matrix F of (3.49) for spectral–collocation approximation of non-
linear shallow water equations, test 1/linear waves.
5.4 Validation of the coupling approach 81
5.4.2 Nonlinear wave
Next, we repeat the same tests 1 and 2 of the above subsection for a larger perturbation, giving
rise to nonlinear waves. Specifically we set η = 1 m with H0N ≡ 10 m. Table 5.11 shows the relative
errors of test 1 at the final time T = 400 s for C ≈ 0.6 (T/∆t = 840) and σ = 500. Relative errors
for the case σ = 1000 are reported in table 5.12. As for test 1, initial depth is shown in figure
5.15 and a comparison between depth and velocity at final time computed by the coupled and full
finite volume scheme is shown in figure 5.16 in the case C = 0.9, θ = 1, σ = 1000, Lax–Friedrichs
flux. Results for C ≈ 0.9 (T/∆t = 560) and C ≈ 0.3 (T/∆t = 1680) are reported in the appendix
(tables A.9–A.10 for σ = 500 A.11–A.12 for σ = 1000). The reported relative errors are higher
that in the linear case, as the higher is the amplitude of the wave, the higher is the amplitude of
the phenomena occurring at the interface, which however do not affect the correct propagation
of the wave through the semi–infinite domain in the coupled scheme. We remark that as above
results are better in the case of the longer wave.
0 0.5 1 1.5 2
x 104
9.6
9.8
10
10.2
10.4
10.6
10.8
11
11.2
x
H(x
)
Initial total depth
σ=500σ=1000
0 0.5 1 1.5 2
x 104
9.6
9.8
10
10.2
10.4
10.6
10.8
11
11.2
x
H(x
)
Initial total depth
σ=500σ=1000
Fig. 5.15: Initial depth perturbations for nonlinear wave test 1 (left) and 2 (right).
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 104
10
10.1
10.2
10.3
10.4
10.5
10.6Depth at time T
x
H(x)
FV−SC
FV
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 104
−0.5
−0.4
−0.3
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5Velocity at time T
x
u(x)
FV−SC
FV
Fig. 5.16: Comparison of coupled (solid black) and full finite volume (dashed red) schemes’ computed
depth (left) and velocity (right), nonlinear wave.
82 5 Numerical Results
Table 5.11: Test 1 results, C ≈ 0.6, σ = 500, nonlinear wave.
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 1.49E-04 1.33E-04 2.23E-04 2.62E-04 1.92E-04 1.93E-04
u 0.034 0.03 0.05 0.059 0.044 0.044
Erel2η 3.19E-04 3.15E-04 4.99E-04 5.74E-04 4.22E-04 4.36E-04
u 0.027 0.026 0.04 0.0037 0.0027 0.0028
Erel∞η 0.0012 0.0012 0.0024 0.0021 0.0018 0.0018
u 0.029 0.027 0.051 3.93E-04 3.47E-04 3.41E-04
Table 5.12: Test 1 results, C ≈ 0.6, σ = 1000, nonlinear wave.
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 5.9E-05 5.1E-05 6.92E-05 1.09E-04 4.6E-05 2.1E-05
u 0.0067 0.0058 0.008 0.012 0.005 0.0024
Erel2η 1.07E-04 9.22E-05 1.41E-04 2.02E-04 8.71E-05 4.06E-05
u 0.0062 0.0053 0.008 6.47E-04 2.79E-04 1.3E-04
Erel∞η 2.91E-04 2.42E-04 4.54E-04 5.93E-04 2.28E-04 1.26E-04
u 0.0063 0.0051 0.01 5.55E-05 2.16E-05 1.18E-05
As before, in the second test we center the initial perturbation η0 at x0 = 5000 m, this time with
η = 1 m, see figure 5.15, and iterate the coupled and the full finite volume scheme taking as final
time T = 1000 s. The other parameters are the same as in the previous case. Absolute L1, L2
and L∞ errors for C ≈ 0.3 (T/∆t = 4200) are reported in table 5.13 for σ = 500. Results for
σ = 1000 are reported in table 5.14. Results for C ≈ 0.9 (T/∆t = 1400) C ≈ 0.6 (T/∆t = 2100)
are reported in the appendix (tables A.13–A.14 for σ = 500, A.15–A.16 for σ = 1000). The results
are consistent with what observed above, noting that the errors are higher than in the linear wave
case, though remaining low. Furthermore, we report in table 5.15 the condition numbers of the
matrix F = IΓ − θ2g∆t2G (see (3.49)) of the semi–infinite spectral collocation system computed
at final time with θ = 1, σ = 1000. As above, the obtained values are very low.
5.4 Validation of the coupling approach 83
Table 5.13: Test 2 results, C ≈ 0.3, σ = 500, nonlinear wave.
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 7.5E-04 9.4E-04 1.41E-03 8.95E-04 1.25E-03 1.8E-03
u 7.4E-04 9.3E-04 1.4E-03 8.8E-04 1.24E-03 1.7E-03
Eabs2η 3.1E-03 3.9E-03 5.2 3.7E-03 5.7E-03 6.5E-03
u 3E-03 3.8E-03 5.1E-03 3.7E-03 5.6E-03 6.5E-03
Eabs∞η 0.028 0.032 0.034 0.037 0.256 0.045
u 0.028 0.032 0.033 0.037 0.257 0.045
Table 5.14: Test 2 results, C ≈ 0.3, σ = 1000, nonlinear wave.
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 1.7E-05 4.7E-05 6.2E-05 2E-05 6.1E-05 8.2E-05
u 1.6E-05 4.6E-05 6.2E-05 2E-05 6.1E-05 8.1E-05
Eabs2η 1.2E-04 3.3E-04 3.9E-04 1.5E-04 4.4E-04 5.2E-04
u 1.2E-04 3.3E-04 3.9E-04 1.5E-04 4.3E-04 5.2E-04
Eabs∞η 0.0025 0.0053 0.0053 0.0033 0.007 0.0071
u 0.0023 0.0052 0.0053 0.0031 0.0069 0.007
N 4 8 16 32 64 128 180
1.0007 1.0025 1.009 1.035 1.135 1.533 2.05
Table 5.15: Condition numbers of the matrix F of (3.49) for spectral–collocation approximation of non-
linear shallow water equations, test 1/nonlinear waves.
84 5 Numerical Results
5.5 Validation of the efficiency for absorbing boundary conditions
In the previous section we have validated our coupling strategy by showing that the amplitude of
reflected waves at the interface is very small, carrying out our analysis in the idealized situation
of a high number of nodes. As the ultimate purpose of our coupled discretization is to use it to
implement absorbing layer–type open boundary conditions, in this section we present the numerical
results obtained in this direction. In our case, the absorbing layer will be simulated by a damping
reaction coefficient placed in the semi–infinite part. Specifically, we want to show that on one
hand, reducing the number of spectral collocation nodes the coupling strategy does not give rise
to numerical instability at the interface, on the other hand that waves can be effectively damped
in the semi–infinite domain in a finite time interval towards the external solution, that we will
assume to be zero both for the elevation and the velocity. To this end, we consider again a finite
domain of length 10 km and evolve through the interface an initially motionless Gaussian elevation
perturbation whose profile is given by the function (5.5). We expect both elevation and velocity to
be erased in the damping process which, however, should not produce spurious waves neither in the
finite nor in the semi–infinite part. The performance of this new absorbing layer open boundary
condition will be assessed in terms of maximum values of residual elevation and velocity after
the damping. As before, we consider a linear and a nonlinear wave, corresponding to the choices
η = 0.3 m and η = 1 m, with the initial average depth given by H0N ≡ 10 m. Moreover, for
the sake of simplicity we will consider only the case of the implicit Euler method (θ = 1), Lax–
Friedrichs flux for the finite volume approximation and linear interpolation for the computation of
the characteristics in semi–Lagrangian method, see section 3.4. For a decreasing number of nodes,
we report the maximum norm of the elevation and velocity vector for the two wavelengths σ = 500
and σ = 1000 and two central values x01 = 4000 m, x02 = 7000 m for the initial perturbation.
Furthermore, to compare results we will take T = 5000 s as a final time. In each test, the grid
spacing in the semi–infinite part will be set by tuning the value of β: as above, the resulting grid
spacing must be of the same order of the finite volume grid spacing at the interface, that is:
xNx+ 12− xNx− 1
2≈ x(β)
1 − x(β)0 .
We recall that the damping is modeled by the coefficient:
γ(x) = γD(∆γ,α,w, L0, x),
a sigmoid of width w centered in x = αL0 connecting the values 0 and ∆γ. In our case we will
set L0 ≈ x(β)N , and w = L0/2. The value of α will be set according to the domain length, in order
to let the damping process begin around the point x0 = 30000 m, which, for the present test,
corresponds approximately to t = T/2.
5.5 Validation of the efficiency for absorbing boundary conditions 85
N = 40
Setting L = 800 (β = 0.00125), we impose L0 ≈ x(β)N ≈ 1.15 · 105, α = 0.6, Nx = 100 and
T/∆t = 2000, so that C ≈ 0.3. Absolute elevation and velocity values are in table 5.16 for linear
wave (η = 0.3 m), in table 5.17 for nonlinear wave (η = 1 m). Errors are low in both cases, with
better results for the shorter wave. Moreover, as expected, a stronger damping entails smaller final
values, whereas the initial position of the wave seems not to have a clear influence on the residual
elevation and velocity after the damping. Indeed, one would expect to have smaller errors for the
wave starting from x02, as it spends more time in the damped region. However, especially for the
longer wave, as the propagation begins near the interface, it may occur that the right crest is not
yet well–formed when it crosses the interface, which may lead to slightly larger errors.
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 1.23E-03 1.74E-03 6.63E-04 8.46E-04
1000 9.14E-04 2.63E-03 9.86E-04 1.35E-03
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 1.7E-05 2.39E-03 4.96E-04 1.12E-03
1000 1.64E-03 3.77E-03 7.48E-04 1.69E-03
Table 5.16: Maximum elevation (left) and velocity (right) norm, N = 40, linear wave
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 4.6E-03 5.32E-03 2.43E-03 2.73E-03
1000 7.72E-03 1.10E-02 3.51E-03 5.05E-03
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 3.39E-03 8.85E-03 1.45E-03 4.23E-03
1000 8.52E-03 1.33E-02 4.06E-03 6.15E-03
Table 5.17: Maximum elevation (left) and velocity (right) norm, N = 40, nonlinear wave
N = 30
Taking again L = 800 (β = 0.00125), we have L0 ≈ x(β)N ≈ 8.48 ·104 m and again we take α = 0.65,
T/∆t = 2000, Nx = 100, so that C ≈ 0.25. Absolute elevation and velocity values are in table
5.18 for linear wave (η = 0.3 m), in table 5.19 for nonlinear wave (η = 1 m). As before, results are
better for the shorter wave and for a stronger damper.
86 5 Numerical Results
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
6σ
500 1.9E-03 3.69E-03 1.19E-03 2.47E-03
1000 3.18E-03 5.76E-03 2.13E-03 3.93E-03
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 1.67E-03 1.76E-03 1.08E-03 9.15E-04
1000 2.67E-03 3.46E-03 1.6E-03 1.9E-03
Table 5.18: Maximum elevation (left) and velocity (right) norm, N = 30, linear wave
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 6.02E-03 1.24E-02 4.09E-03 8.58E-03
1000 9.82E-03 1.78E-02 7.33E-03 1.3E-02
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 6.25E-03 6.56E-03 3.01E-03 3.35E-03
1000 1.04E-02 1.31E-02 4.94E-03 9.86E-03
Table 5.19: Maximum elevation (left) and velocity (right) norm, N = 30, nonlinear wave
N = 20
Taking again L = 800 (β = 0.00125), we have L0 ≈ x(β)N ≈ 5.47 · 104 m, α = 0.7, and take
T/∆t = 1200, Nx = 55, so that C ≈ 0.25. Absolute elevation and velocity values are in table 5.20
for linear wave (η = 0.3 m), in table 5.21 for nonlinear wave (η = 1 m). As before, results are better
for the shorter wave and for a stronger damper, though the damping starts to be difficult with
the coarse resolution in the last part of the semi–infinite domain. A fine tuning of the coefficient γ
with different parameters can improve results, here we choose to be consistent with the previous
tests.
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 1.04E-03 8.8E-04 4.47E-04 5.38E-04
1000 2.27E-03 1.96E-03 9.19E-04 1.01E-03
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 2.01E-03 3.97E-03 1.5E-03 1.93E-03
1000 3.52E-03 6.92E-03 2.79E-03 3.39E-03
Table 5.20: Maximum elevation (left) and velocity (right) norm, N = 20, linear wave
5.5 Validation of the efficiency for absorbing boundary conditions 87
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 3.87E-03 3.46E-03 1.38E-03 1.77E-03
1000 8.85E-03 8.08E-03 2.75E-03 4.05E-03
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 6.26E-03 1.31E-02 4.98E-03 6.66E-03
1000 9.92E-03 2.19E-02 9.02E-03 1.16E-02
Table 5.21: Maximum elevation (left) and velocity (right) norm, N = 20, nonlinear wave
N = 10
In this last test, the parameters will be slightly modified because of the very coarse grid in the
semi–infinite part. We take L = 700 (β = 0.00143), so that L0 ≈ x(β)N ≈ 2.22 · 104 m and α = 0.7,
T/∆t = 2500, Nx = 60, so that C ≈ 0.1. Absolute elevation and velocity values are in table
5.22 for linear wave (η = 0.3 m), in table 5.23 for nonlinear wave (η = 1 m). As before, results
are better for the shorter wave and for a stronger damper. These results on the use of our semi–
infinite discretization as an absorbing layer, and especially this last case of 10 nodes, will make
possible a comparison with currently used absorbing layers (see e.g. [LT08, MDD10]). However,
we remark that standard layers are very thin with respect to the finite domain length, whereas
our scaled approach would enable to represent a larger semi–infinite part approximately at the
same computational cost (employing the same number of nodes).
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 9E-04 1.18E-03 1.07E-03 7.36E-04
1000 2.33E-03 1.99E-03 2.67E-03 1.4E-03
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 1.42E-03 4.99E-04 1.11E-03 6.93E-04
1000 3.56E-03 1.09E-03 3.25E-03 1.22E-03
Table 5.22: Maximum elevation (left) and velocity (right) norm, N = 10, linear wave
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 2.85E-03 3.81E-03 3.48E-03 2.45E-03
1000 7.25E-03 6.16E-03 8.64E-03 4.58E-03
∆γ = 1 ∆γ = 2
x01 x02 x01 x02
σ500 4.7E-03 1.89E-03 3.59E-03 2.16E-03
1000 1.18E-02 4.31E-03 1.05E-03 3.54E-03
Table 5.23: Maximum elevation (left) and velocity (right) norm, N = 10, nonlinear wave
Conclusions and future perspectives
In this thesis we have studied a numerical method based on Laguerre functions expansions
for the approximation of wave propagation problems on semi–infinite domains. In particular,
a semi–implicit and semi–Lagrangian discretization has been developed for the one dimensional
shallow water model equations on semi–infinite domains, along with a coupling approach based
on imposing continuity of the mass fluxes at the interface between finite size and unbounded do-
main, that can be in principle extended to couple the Laguerre spectral method to an arbitrary
discretization on the finite size domain. The semi–implicit, semi–Lagrangian approach has been
chosen because of its well known efficiency and accuracy features and because a large number of
environmental models use these time discretization techniques. On the other hand, in this first at-
tempt, an explicit finite volume discretization has been considered for the finite size domain, both
for the sake of simplicity and in order to show that such a coupling is feasible also if two essentially
different discretization approaches are employed on either domain. The proposed pseudo–spectral
method has also been employed to implement absorbing layer open boundary conditions. By cou-
pling discretizations on domains of finite size to Laguerre spectral methods employing a relatively
small number of nodes, an efficient approach to absorbing boundary conditions is achieved, that
allows to account for large absorbing regions at a relatively low computational cost.
A number of numerical tests has been carried out. Firstly, the pseudo–spectral, semi–implicit and
semi–Lagrangian discretization approach on semi–infinite domains has been validated as an inde-
pendent tool, considering various kinds of boundary conditions at the origin of the half line. Then,
the coupling approach has been validated, by comparison of the coupled finite volume/pseudo–
spectral model results to those of reference runs of the finite volume model on a large domain.
Numerical tests have shown that spurious reflections due to the proposed coupling approach have
very small amplitude and allow to use the proposed method, at least for environmental applica-
tions, with no substantial loss of accuracy with respect to discretizations using a single computa-
tional domain. Finally, absorbing boundary conditions have been validated, as implemented in the
context of the coupled finite volume/pseudo–spectral model. Numerical simulations have shown
90 Conclusions and future perspectives
that a reasonable number of spectral base functions is sufficient to reach the accuracy necessary
for practical implementations of open boundary conditions, along with an accurate treatment of
the solution behaviour at infinity.
All the results are quite promising and appear to confirm that the proposed approach to absorbing
layer boundary conditions is advantageous for a number of reasons. It is therefore of interest to
consider various natural extensions of the present thesis work and future developments that could
lead to exploit the full benefit of the proposed technique in practical applications. In particular,
the following points should be investigated:
• The proposed semi–implicit, semi–Lagrangian pseudo–spectral method on semi–infinite do-
mains can be coupled to different discretizations on the finite size domain; in particular, ef-
ficient and commonly used semi–implicit, semi–Lagrangian techniques could be considered in
the finite domain, based on either a finite difference or finite element spatial discretization.
• Two dimensional discretizations of the shallow water equations can be developed, in which a
semi–infinite strip is discretized using Laguerre basis functions in the x direction and some
other discretization approach in the y direction. For example, tensor product–based multi–
dimensional polynomial spaces as in [ZSX10] could be employed, or, more simply, finite differ-
ence or finite volume discretizations could be used. This would naturally lead to comparison
with the results in [LT08, MDD10] on two dimensional test cases.
• Two dimensional, vertical slice atmospheric models could be considered (see e.g. [Bon00, GR08,
GR10]), to which an extension to unbounded domains could be coupled along the lines described
in this work. The main target of this further development would be the comparison of the
standard absorbing layer approaches (see e.g. [KD83, KDH08]) to those built following the
ideas proposed in this work. For these numerical models, the width of the necessary upper
absorbing layer is known to increase with the increase in horizontal resolution, thus leading
to an extra computational cost that can be especially high for models employing height based
coordinates, such as many high resolution numerical weather prediction models. Thus, it would
be interesting to assess whether a substantial decrease of the computational cost for these
models can be achieved by extending the techniques proposed in this thesis.
A
Appendix
In this appendix we report the extended numerical results for the validation of coupling ap-
proach of sections 4.2 and 5.4. Results for the linear wave are reported in tables A.1–A.2,
A.3–A.4 for test 1 (C ≈ 0.9, C ≈ 0.3, respectively), in tables A.5–A.6, A.7–A.8 for test 2 (C ≈ 0.9,
C ≈ 0.6, respectively). For the same Courant numbers, results for the nonlinear wave are reported
in tables A.9–A.10,A.11–A.12 for test 1 (C ≈ 0.9, C ≈ 0.3, respectively) in tables A.13–A.14,
A.15–A.16 for test 2 (C ≈ 0.9, C ≈ 0.6, respectively). We recall that test 2 absolute L1 and L2
errors are rescaled by the domain length.
Table A.1: Test 1 results, C ≈ 0.9, σ = 500, linear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 5.78E-05 6.59E-05 7.34E-05 2.36E-05 2.42E-05 2.58E-05
u 0.043 0.05 0.056 0.018 0.018 0.019
Erel2η 1.41E-04 1.61E-04 1.79E-04 4.83E-05 5.33E-05 6.18E-05
u 0.038 0.043 0.048 0.013 0.014 0.017
Erel∞η 6.16E-04 7.05E-04 7.97E-04 1.64E-04 1.88E-04 2.43E-04
u 0.043 0.048 0.054 0.012 0.0129 0.017
92 A Appendix
Table A.2: Test 1 results, C ≈ 0.3, σ = 500, linear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 1.42E-04 3.45E-05 3.21E-05 1.56E-04 5.28E-05 2.01E-05
u 0.107 0.026 0.024 0.117 0.04 0.015
Erel2η 3.2E-04 7.6E-05 7.14E-05 3.53E-04 1.17E-04 4.63E-05
u 0.098 0.021 0.019 0.108 0.033 0.012
Erel∞η 0.001 3.23E-04 5.62E-04 0.001 5.17E-04 1.82E-04
u 0.114 0.024 0.022 0.126 0.038 0.013
Table A.3: Test 1 results, C ≈ 0.9, σ = 1000, linear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 3.01E-05 3.31E-05 3.55E-05 1.18E-05 1.3E-05 1.36E-05
u 0.011 0.013 0.013 0.004 0.005 0.005
Erel2η 5.97E-04 6.65E-05 7.18E-04 2.14E-05 2.32E-05 2.43E-05
u 0.011 0.013 0.014 0.004 0.004 0.005
Erel∞η 1.68E-04 1.88E-04 7.97E-04 5.27E-04 5.61E-05 5.95E-05
u 0.011 0.013 0.014 0.004 0.004 0.004
Table A.4: Test 1 results, C ≈ 0.3, σ = 1000, linear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 8.94E-05 1.73E-05 1.12E-05 9.78E-05 2.55E-05 3.024E-06
u 0.034 0.007 0.004 0.037 0.01 0.001
Erel2η 1.72E-04 3.25E-05 2.29E-05 1.88E-04 4.81E-05 6.15E-06
u 0.034 0.006 0.004 0.037 0.009 0.001
Erel∞η 5.17E-04 1.01E-04 7.63E-05 5.65E-04 1.98E-04 2.39E-05
u 0.038 0.007 0.005 0.041 0.014 0.002
A Appendix 93
Table A.5: Test 2 results, C ≈ 0.9, σ = 500, linear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 1.18E-05 1.48E-05 1.59E-05 3.45E-05 4.34E-05 4.68E-05
u 1.17E-05 1.47E-05 1.57E-05 3.42E-05 4.31E-05 4.64E-05
Eabs2η 7E-05 9E-05 9E-05 2.2E-04 2.6E-04 2.7E-04
u 7E-05 9E-05 9E-05 2.1E-04 2.6E-04 2.7E-04
Eabs∞η 8.95E-04 0.001 0.001 0.0026 0.003 0.003
u 8.87E-04 0.001 0.001 0.0026 0.003 0.003
Table A.6: Test 2 results, C ≈ 0.6, σ = 500, linear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 5.9E-06 1.4E-05 1.75E-05 1.17E-05 2.85E-05 3.58E-05
u 5.9E-06 1.38E-05 1.73E-05 1.16E-05 2.83E-05 3.54E-05
Eabs2η 5E-05 9E-05 1.1E-04 9E-05 1.9E-04 2.2E-04
u 4E-05 9E-05 1E-04 9E-05 1.9E-04 2.1E-04
Eabs∞η 6.85E-04 0.0012 0.0012 0.0014 0.0025 0.0026
u 6.78E-04 0.0012 0.0012 0.0014 0.0025 0.0025
Table A.7: Test 2 results, C ≈ 0.9, σ = 1000, linear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 3.2E-07 3.4E-07 3.4E-07 6.1E-07 6.8E-07 7E-07
u 3.2E-07 3.3E-07 3,4E-07 6E-07 6.8E-07 7E-07
Eabs2η 1.83E-06 2.02E-06 2.05E-06 4.24E-06 4.79E-06 4.88E-06
u 1.82E-06 2E-06 2.03E-06 4.2E-06 4.74E-06 4.83E-06
Eabs∞η 3.35E-05 3.62E-05 3.62E-05 7.66E-05 8.46E-05 8.47E-05
u 3.31E-05 3.58E-05 3.58E-05 7.58E-05 8.38E-05 8.39E-05
94 A Appendix
Table A.8: Test 2 results, C ≈ 0.6, σ = 1000, linear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 4E-07 3.9E-07 4.2E-07 5E-07 6E-07 6.7E-07
u 4E-07 3,9E-07 4.2E-07 4.9E-07 6E-07 6.6E-07
Eabs2η 1.98E-06 2.42E-06 2.56E-06 2.97E-06 4.26E-06 4.55E-06
u 1.96E-06 2.4E-06 2.54E-06 2.94E-06 4.22E-06 4.5E-06
Eabs∞η 3.59E-05 4.53E-05 4.54E-05 6.08E-05 8.03E-05 8.05E-05
u 3.5E-05 4.49E-05 4.49E-05 5.93E-05 7.96E-05 7.97E-05
Table A.9: Test 1 results, C ≈ 0.9, σ = 500, nonlinear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 2.24E-04 2.5E-04 2.74E-04 1.76E-04 1.83E-04 1.93E-04
u 0.051 0.057 0.062 0.04 0.042 0.044
Erel2η 5.27E-04 5.84E-04 6.39E-04 3.99E-04 4.18E-04 4.41E-04
u 0.043 0.047 0.051 0.0026 0.0027 0.0028
Erel∞η 0.003 0.0029 0.0032 0.0017 0.0017 0.0017
u 0.056 0.062 0.068 3.33E-04 3.32E-04 3.31E-04
Table A.10: Test 1 results, C ≈ 0.3, σ = 500, nonlinear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 4.72E-04 1.74E-04 1.93E-04 5.23E-04 2.33E-04 1.97E-04
u 0.106 0.039 0.044 0.118 0.053 0.045
Erel2η 0.001 3.8E-04 4.22E-04 0.0011 5.13E-04 4.43E-04
u 0.094 0.032 0.034 0.0072 0.0033 0.0029
Erel∞η 0.004 0.0015 0.0016 0.0043 0.0019 0.0019
u 0.108 0.036 0.035 7.98E-04 3.64E-04 3.53E-04
A Appendix 95
Table A.11: Test 1 results, C ≈ 0.9, σ = 1000, nonlinear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 9.79E-05 1.06E-04 1.11E-04 4.11E-05 4.4E-05 4.16E-05
u 0.011 0.012 0.013 0.0047 0.005 0.0047
Erel2η 1.95E-04 2.13E-04 2.23E-04 7.45E-05 7.89E-05 7.68E-05
u 0.011 0.012 0.013 2.39E-04 2.54E-04 2.47E-04
Erel∞η 5.86E-04 6.31E-04 6.67E-04 1.91E-04 2.11E-04 2.21E-04
u 0.012 0.013 0.014 1.81E-05 1.99E-05 2.08E-05
Table A.12: Test 1 results, C ≈ 0.3, σ = 1000, nonlinear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Erel1η 2.81E-04 5.51E-05 3.3E-05 3.07E-04 8.16E-05 1.88E-05
u 0.032 0.0063 0.0037 0.035 0.009 0.0021
Erel2η 5.3E-04 1.03E-04 7.08E-05 5.81E-04 1.6E-04 3.58E-05
u 0.032 0.006 0.004 0.0019 5.15E-04 1.15E-04
Erel∞η 0.0015 2.87E-04 2.48E-04 0.017 0.0015 1.05E-04
u 0.035 0.006 0.0052 1.59E-04 1.41E-04 9.91E-06
Table A.13: Test 2 results, C ≈ 0.9, σ = 500, nonlinear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 8.24E-04 9.1E-04 9.84E-04 1.36E-03 1.49E-03 1.62E-03
u 8.14E-04 8.98E-04 9.7E-04 1.35E-03 1.48E-03 1.6E-03
Eabs2η 3.1E-03 3.46E-03 3.91E-03 5.24E-03 5.64E-03 6.2E-03
u 3E-03 3.4E-03 3.84E-03 5.2E-03 5.6E-03 6.1E-03
Eabs∞η 0.044 0.07 0.154 0.037 0.041 0.137
u 0.043 0.069 0.15 0.036 0.04 0.134
96 A Appendix
Table A.14: Test 2 results, C ≈ 0.6, σ = 500, nonlinear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 6.2E-04 8.4E-04 1.15E-03 1E-03 1.3E-03 1.7E-03
u 6.1E-04 8.3E-04 1.13E-03 1E-03 1.3E-03 1.7E-03
Eabs2η 2.6E-03 3.2E-03 4.1E-03 4.2E-03 5.2E-03 6.3E-03
u 2.5E-03 3.2E-03 4.1E-03 4.2E-03 5.1E-03 6.2E-03
Eabs∞η 0.022 0.023 0.059 0.038 0.04 0.043
u 0.022 0.023 0.058 0.038 0.04 0.042
Table A.15: Test 2 results, C ≈ 0.9, σ = 1000, nonlinear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 3.4E-05 3.9E-05 4.2E-05 7.8E-05 9E-05 9.6E-05
u 3.4E-05 3.9E-05 4.1E-05 7.7E-05 8.9E-05 9.5E-05
Eabs2η 2.2E-04 2.5E-04 2.6E-04 5.1E-04 5.7E-04 5.9E-04
u 2.2E-04 2.5E-04 2.6E-04 5.1E-04 5.7E-04 5.9E-04
Eabs∞η 0.003 0.0034 0.0034 0.0068 0.0076 0.0076
u 0.003 0.0033 0.0033 0.0068 0.0075 0.0075
Table A.16: Test 2 results, C ≈ 0.6, σ = 1000, nonlinear wave
θ = 1 θ = 0.6
LF GD SL–MM LF GD SL–MM
Eabs1η 2.64E-05 4.2E-05 5.1E-05 4.4E-05 7.3E-05 8.9E-05
u 2.6E-05 4.2E-05 5E-05 4.4E-05 7.2E-05 8.8E-05
Eabs2η 1.9E-04 2.9E-04 3.2E-04 3.3E-04 5E-04 5.6E-04
u 1.9E-04 2.8E-04 3.2E-04 3.3E-04 4.9E-04 5.5E-04
Eabs∞η 0.0033 0.0042 0.0042 0.0057 0.0073 0.0074
u 0.0032 0.0042 0.0042 0.0055 0.0072 0.0073
References
[Amb95] D. Ambrosi. Approximation of shallow water equations by Roe’s Riemann solver.
International Journal for Numerical Methods in Fluids, 20(2):157–168, 1995.
[Ast00] R.J. Astley. Infinite elements for wave problems: a review of current formulations
and an assessment of accuracy. International Journal for Numerical Methods in En-
gineering, 49(7):951–976, 2000.
[BK00] J. Billingham and A.C. King. Wave motion, volume 24 of Cambridge texts in applied
mathematics. Cambridge University Press, 2000.
[Bon00] L. Bonaventura. A Semi–implicit Semi–Lagrangian Scheme Using the Height Coor-
dinate for a Nonhydrostatic and Fully Elastic Model of Atmospheric Flows. Journal
of Computational Physics, 158(2):186–213, 2000.
[Cas90] V. Casulli. Semi–implicit finite difference methods for the two–dimensional shallow
water equations. Journal of Computational Physics, 86(1):56–74, 1990.
[CDMW98] S. Chippada, C.N. Dawson, M.L. Martinez, and M.F. Wheeler. A Godunov–type
finite volume method for the system of shallow water equations. Computer Methods
in Applied Mechanics and Engineering, 151(1–2):105–129, 1998.
[CHQZ06] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods. Fun-
damentals in Single Domains. Springer–Verlag, Berlin Heidelberg, 2006.
[CHQZ07] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral methods. Evo-
lution to Complex Geometries and Application to Fluid Dynamics. Springer–Verlag,
Berlin Heidelberg, 2007.
[DBMS09] A. Decoene, L. Bonaventura, E. Miglio, and F. Saleri. Asymptotic derivation of the
section–averaged shallow water equations for natural river hydraulics. Mathematical
Models and Methods in Applied Sciences, 19(3):387–417, 2009.
[DK83] D.R. Durran and J.B. Klemp. A compressible model for the simulation of moist
mountain waves. Monthly Weather Review, 111(12):2341–2361, 1983.
98 References
[EM77] B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simu-
lation of waves. Mathematics of Computation, 31(139):629–651, 1977.
[Fra72] P.E. Francis. The possible use of Laguerre polynomials for representing the vertical
structure of numerical models of the atmosphere. Quarterly Journal of the Royal
Meteorological Society, 98(417):662–667, 1972.
[Ger00] K. Gerdes. A review of infinite element methods for exterior Helmholtz problems.
Journal of Computational Acoustics, 8(1):43–62, 2000.
[Giv08] D. Givoli. Computational Absorbing Boundaries. In S. Marburg and B. Nolte, edi-
tors, Computational Acoustics of Noise Propagation in Fluids – Finite and Boundary
Element Methods. Springer Berlin Heidelberg, 2008.
[GMR+06] M.A. Giorgetta, E. Manzini, E. Roeckner, M. Esch, and L. Bengtsson. Cilmatology
and forcing of the quasi–biennal oscillation in the MAECHAM5 model. Journal of
Climate, 19(16):3882–3901, 2006.
[GO77] D. Gottlieb and S.A. Orszag. Numerical analysis of spectral methods: theory and
applications. SIAM–CBMS, Philadelphia, 1977.
[GR08] F.X. Giraldo and M. Restelli. A study of spectral element and discontinuous
Galerkin methods for the Navier–Stokes equations in nonhydrostatic mesoscale at-
mospheric modeling: Equation sets and test cases. Journal of Computational Physics,
227(8):3849–3877, 2008.
[GR10] F.X. Giraldo and M. Restelli. High–order semi–implicit time–integrators of a trian-
gular discontinuous Galerkin oceanic shallow water model. International Journal for
Numerical Methods in Fluids, 63(9):1077–1102, 2010.
[GS00] B.-Y. Guo and J. Shen. Laguerre–Galerkin method for nonlinear partial differential
equations on a semi–infinite interval. Numerische Mathematik, 86(4):635–654, 2000.
[GW07] B.-Y. Guo and Z.-Q. Wang. Numerical integration based on Laguerre–Gauss interpo-
lation. Computer Methods in Applied Mechanics and Engineering, 196(37–40):3726–
3741, 2007.
[GWW06] B.-Y. Guo, L.-L. Wang, and Z.-Q. Wang. Generalized Laguerre interpolation and
pseudospectral method for unbounded domains. SIAM Journal on Numerical Anal-
ysis, 43(6):2567–2589, 2006.
[HPS+03] T. Horinouchi, S. Pawson, K. Shibata, U. Langematz, E. Manzini, M.A. Giorgetta,
F. Sassi, R.J. Wilson, K. Hamilton, J. De Grandpré, and A.A. Scaife. Tropical cumulus
convection and upward–propagating waves in middle–atmospheric GCMs. Journal of
the Atmospheric Sciences, 60(22):2765–2782, 2003.
[IO81] M. Israeli and S.A. Orszag. Approximation of radiation boundary conditions. Journal
of Computational Physics, 41(1):115–135, 1981.
References 99
[KD83] J.B. Klemp and D.R. Durran. An upper boundary condition permitting internal
gravity wave radiation in numerical mesoscale models. Monthly Weather Review,
111(3):430–444, 1983.
[KDH08] J.B. Klemp, J. Dudhia, and A.D. Hassiotis. An upper gravity–wave absorbing layer
for NWP applications. Monthly Weather Review, 136(10):3987–4004, 2008.
[KL78] J.B. Klemp and D.K. Lilly. Numerical simulation of hydrostatic mountain waves.
Journal of the Atmospheric Sciences, 35(1):78–107, 1978.
[LeV92] R.J. LeVeque. Numerical methods for conservation laws. Lectures in Mathematics
ETH Zürich. Birkhäuser, 1992.
[LeV02] R.J. LeVeque. Finite volume methods for hyperbolic problems, volume 31 of Cambridge
texts in applied mathematics. Cambridge University Press, Cambridge, 2002.
[LG08] R.J. LeVeque and D.L. George. High–resolution finite volume methods for the shallow
water equations with bathymetry and dry states. In P. L.-F. Liu, H.H. Yeh, and
C. Synolakis, editors, Advanced numerical models for simulating tsunami waves and
runup, volume 10 of Advances in coastal and ocean engineering. World Scientific, 2008.
[LT08] J.W. Lavelle and W.C. Thacker. A pretty good sponge: Dealing with open boundaries
in limited–area ocean models. Ocean Modelling, 20(3):270–292, 2008.
[McD03] A. McDonald. Transparent boundary conditions for the shallow–water equations:
Testing in a nested environment. Monthly Weather Review, 131(4):698–705, 2003.
[MDD10] A. Modave, É. Deleersnijder, and É.J.M. Delhez. On the parameters of absorbing
layers for shallow water models. Ocean Dynamics, 60(1):65–79, 2010.
[NNH04] I.M. Navon, B. Neta, and M. Y. Hussaini. A perfectly matched layer approach to the
linearized shallow water equations models. Monthly Weather Review, 132(6):1369–
1378, 2004.
[QSS07] A. Quarteroni, R. Sacco, and F. Saleri. Numerical mathematics, volume 37 of Texts
in applied mathematics. Springer, Berlin Heidelberg, Second edition, 2007.
[Qua90] A. Quarteroni. Domain decomposition methods for systems of conservation laws:
Spectral collocation approximations. SIAM Journal on Scientific and Statistical Com-
puting, 11(6):1029–1052, 1990.
[QV99] A. Quarteroni and A. Valli. Domain decomposition methods for partial differential
equations. Oxford University Press, 1999.
[QV08] A. Quarteroni and A. Valli. Numerical Approximation of Partial Differential Equa-
tions, volume 23 of Springer series in computational mathematics. Springer–Verlag,
Berlin Heidelberg, 2008.
100 References
[Ras86] P.J. Rasch. Toward atmospheres without tops: Absorbing upper boundary condi-
tions for numerical models. Quarterly Journal of the Royal Meteorological Society,
112(474):1195–1218, 1986.
[RCB05] G. Rosatti, D. Cesari, and L. Bonaventura. Semi–implicit, semi–Lagrangian modelling
for environmental problems on staggered Cartesian grids with cut cells. Journal of
Computational Physics, 204(1):353–377, 2005.
[Roe81] P.L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes.
Journal of Computational Physics, 43(2):357–372, 1981.
[Sal08] S. Salsa. Partial differential equations in action: from modelling to theory. Springer–
Verlag Italia, Milan, 2008.
[San04] G. Sansone. Orthogonal Functions. Courier Dover Publications Inc., New York, 2004.
[SC91] A. Staniforth and J. Côté. Semi–Lagrangian integration schemes for atmospheric
models–a review. Monthly Weather Review, 119(9):2206–2231, 1991.
[She01] J. Shen. Stable and efficient spectral methods in unbounded domains using Laguerre
functions. SIAM Journal on Numerical Analysis, 38(4):1113–1133, 2001.
[ST06] J. Shen and T. Tang. Spectral and high order methods with applications. Science
Press, Beijing, 2006.
[Sto92] J.J. Stoker. Water waves: the mathematical theory with applications, volume 36 of
Wiley classics library. John Wiley and Sons, 1992.
[SW09] J. Shen and L.-L. Wang. Some recent advances on spectral methods for unbounded
domains. Communications in Computational Physics, 5(2–4):195–241, 2009.
[Tor01] E.F. Toro. Shock–capturing methods for free–surface shallow flows. John Wiley, New
York, 2001.
[WGW09] Z.-Q. Wang, B.-Y. Guo, and Y.-N. Wu. Pseudospectral method using generalized
Laguerre functions for singular problems on unbounded domains. Discrete and Con-
tinuous Dynamical Systems Series B, 11(4):1019–1038, 2009.
[Whi99] G. B. Whitham. Linear and nonlinear waves. Wiley, 1999.
[YCC10] C. Yang, J. Cao, and X.-C. Cai. A fully implicit domain decomposition algorithm for
shallow water equations on the cubed–sphere. SIAM Journal on Scientific Computing,
32(1):418–438, 2010.
[ZSX10] Q. Zhuang, J. Shen, and C. Xu. A coupled Legendre–Laguerre spectral–element
method for the Navier–Stokes equations in unbounded domains. Journal of Scientific
Computing, 42(1):1–22, 2010.