POLITECNICODIMILANO FacoltàdiIngegneriadeiSistemibonavent/Homepage_Luca_file/benacchio.pdf ·...

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:


ω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(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(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)


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),


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)


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,+∞).


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:


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

,


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);


ω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.


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.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:


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:


(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:


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

,


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:


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:


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)


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


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.


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:


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:


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:


∂η

∂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


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)


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:


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)


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).


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).


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 β.


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.


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).


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.


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.


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


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.


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.


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


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.


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).


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


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


θ = 1 θ = 0.6


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


(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.


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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.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.


Table 5.11: Test 1 results, C ≈ 0.6, σ = 500, nonlinear wave.

θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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.



θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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.


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.


∆γ = 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


∆γ = 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


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


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


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


∆γ = 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


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


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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


θ = 1 θ = 0.6


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, É. Deleersnĳder, 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, Beĳing, 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.

POLITECNICODIMILANO FacoltàdiIngegneriadeiSistemibonavent/Homepage_Luca_file/benacchio.pdf ·...

Documents

Transcript of POLITECNICODIMILANO FacoltàdiIngegneriadeiSistemibonavent/Homepage_Luca_file/benacchio.pdf ·...