CUBO A Mathematical Journal

Vol.11, No¯ 04, (29–48). September 2009

Towards accurate artificial boundary conditions for nonlinear
PDEs through examples

Xavier Antoine,

Institut Elie Cartan Nancy,
Nancy-Université,

CNRS, INRIA Corida Team,
Boulevard des Aiguillettes B.P. 239 F-54506 Vandoeuvre-lès-Nancy, France.

email: Xavier.Antoine@iecn.u-nancy.fr
Christophe Besse

Projet Simpaf-Inria Futurs,
Laboratoire Paul Painlevé, Unité Mixte de Recherche CNRS (UMR 8524),

UFR de Mathématiques Pures et Appliquées,
Université des Sciences et Technologies de Lille,

Cité Scientifique, 59655 Villeneuve d’Ascq Cedex, France.
email: Christophe.Besse@math.univ-lille1.fr

and

Jérémie Szeftel 1

Department of Mathematics,
Princeton University, Fine Hall,

Washington Road, Princeton NJ 08544-1000, USA. and
CNRS, Mathématiques Appliquées de Bordeaux,

Université Bordeaux 1,
351 cours de la Libération, 33405 Talence cedex France.

email: jszeftel@math.princeton.edu

ABSTRACT

The aim of this paper is to give a comprehensive review of current developments re-

lated to the derivation of artificial boundary conditions for nonlinear partial differential

equations. The essential tools to build such boundary conditions are based on pseudod-

ifferential and paradifferential calculus. We present various derivations and compare

1partially supported by NSF Grant DMS-0504720.



30 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

them. Some numerical results illustrate their respective accuracy and analyze the po-

tential of each technique.

RESUMEN

La meta de este artículo es entregar una revisión comprensiva de los desarrollos ac-

tuales relacionados con la derivación de condiciones de borde artificiales para ecua-

ciones diferenciales parciales nolineales. Las herramientas esenciales para construir

tales condiciones de borde se basan en el cálculo pseudodiferencial y paradiferencial.

Presentamos varias derivaciones y las comparamos. Algunos resultados numéricos ilus-

tran su precisión respectiva y se analiza el potencial de cada técnica.

Key words and phrases: Nonlinear PDEs, wave equation, Schrödinger equation, artificial bound-
ary conditions for nonlinear PDEs, numerical schemes.

Math. Subj. Class.: 35A21, 35A27, 35L05, 35Q55, 35S50, 65M99.

1 Introduction

The subject of designing artificial boundary conditions for linear scalar and systems of Partial Dif-
ferential Equations (PDEs) has been studied since more than thirty years now. Essentially, the goal
of these boundary conditions is to truncate an infinite domain into a finite one for computational

considerations. To this end, a fictitious boundary Γ is introduced delimiting therefore a finite do-

main Ω. All the difficulty is then to build an accurate, robust and easy-to-implement approximate

boundary condition on this fictitious boundary. These boundary conditions can be found in the

literature under different names (which in fact have different subtle meanings) like artificial bound-

ary conditions, absorbing boundary conditions, non-reflecting or transparent boundary conditions

and sometimes Dirichlet-to-Neumann operators. Among the major contributions written on the

topic and without being exhaustive, let us quote e.g. the papers by Engquist and Majda [9, 10],

Bayliss and Turkel [4], Mur [15] or also Bérenger [5]. A few review papers have also been published

to establish the current state-of-the-art on the subject (see e.g. [2, 11, 12, 22]).

While many improvements have been achieved over the past years, most of them are devel-

oped for linear equations. Practically, most available methods for linear equations do not apply to

nonlinear equations since they often rely on the explicit computation of the transparent boundary

condition by using the Fourier or Laplace transforms (note however that in the particular case

of integrable equations, one may use the inverse scattering transform to explicitly compute the

transparent boundary condition, see e.g. [23]). Now, Nonlinear problems have recently received
some increasing special care because of their importance in applications like e.g. in wave propa-

gation, quantum mechanics, fluid mechanics,... The aim of this paper is to give a comprehensive



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 31

introduction to possible solutions to handle such nonlinear situations. They are mainly based on

pseudodifferential [21] and paradifferential operators techniques [6].

The paper is organized as follows. In Section 2, we analyze in detail the way of construct-

ing approximate artificial boundary conditions for a general one-dimensional wave equation with

variable coefficients using pseudodifferential calculus. Then, we test numerically these absorbing

boundary conditions on a model problem showing that they yield accurate computations at least

for small times. In a third Section, we consider a general one-dimensional nonlinear Schrödinger

equation. We present several ways to extend the method of Section 2 to this nonlinear equation

depending on the kind of nonlinearity involved in the equation. The various types of absorbing

boundary conditions are obtained using either the pseudodifferential or the paradifferential calcu-

lus. In a fourth Section, some numerical comparisons are developed to test the accuracy of the

various absorbing boundary conditions. The last Section draws a conclusion and suggests some

future directions of research.

2 Artificial boundary conditions for linear variable coeffi-

cients equations: the case of the wave equation

2.1 The case of the constant coefficients wave equation

Before directly going to the case of the wave equation with variable coefficients, let us first consider

the simple wave equation

(∂
2

t − ∂2x)u = 0, (2.1)

with initial data u(0,x) = u0(x) and ∂tu(0,x) = u1(x), where the field u = u(t,x) is defined on the

whole space (t,x) ∈ [0; +∞[×R. For simulation purposes, it is standard to introduce a bounded
spatial computational domain setting now (t,x) ∈ [0; +∞[×[xℓ; xr], where xℓ (respectively xr) is
a left (respectively right) fictitious endpoint introduced to get a bounded domain Ω = [xℓ; xr]. Let

us assume that the initial data of our problem, that is u0 and u1, are compactly supported in Ω.

Then, we can define the extension of u (which is still denoted by u) for negative times by fixing its

value to zero so that u is a solution to (2.1) for all times t as long as x ∈ Ωc, where Ωc = R−Ω. Let
us denote by ût(τ,x) = Ft(u)(τ,x) the partial time-Fourier transform of u, where τ ∈ R. Applying
Ft to (2.1) for (t,x) ∈ R × Ωc leads to the Helmholtz-type constant coefficients equation

(∂
2

x + τ
2
)ût(τ,x) = 0, (2.2)

where the wavenumber is τ. The solution of this equation can be written as the superposition of

two waves

ût(τ,x) = A
+
e
iτx

+ A
−
e
−iτx

, (2.3)

where A± are two smooth functions depending on τ. Computing the derivative ∂xût, we obtain

(∂x − iτ)ût = −2iτe−iτxA−, (2.4)



32 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

and

(∂x + iτ)ût = 2iτe
iτx
A

+
, (2.5)

and we obviously check that: (∂x + iτ)(∂x − iτ)ût = 0. We also have the following operator
factorization

∂
2

t − ∂2x = −(∂x + ∂t)(∂x − ∂t). (2.6)

Equation (2.4) (respectively (2.5)) gives a characterization of the right (respectively left) traveling

solution to (2.1) by setting: A− = 0 (respectively A+ = 0). Therefore, the following boundary

condition (2.4)

(∂
n
− iτ)ût = 0, at Γ, (2.7)

acts as a filter in the time-Fourier domain and translates the property that there is no reflection back

into the computational domain Ω, where n is the unit normal vector to Γ = {xℓ; xr}, outwardly
directed to Ω. This is a constraint which forces the wave to be outgoing to Ω. In the time-space

domain, the corresponding boundary condition writes down

(∂
n
− ∂t)u = 0, at [0; +∞[×Γ. (2.8)

Since there is no reflection, this boundary condition is usually called Transparent or Non-Reflecting

Boundary Condition (TBC). Let us remark at this point that another interpretation of writing a

transparent boundary condition is that we require u ∈ Ker(∂
n
− b), setting b(x,t,∂t) = ∂t.

2.2 The case of the variable coefficients wave equation

Let us consider that α, β and γ are three C∞ functions. Writing a TBC for a variable coefficient
model wave equation

(∂
2

t + β(t,x)∂t − ∂2x + γ(t,x)∂x + α(t,x))u = 0, (2.9)

is much more complicate than in the constant coefficients case. Indeed, in such a situation

i) directly applying a time-Fourier transform to the equation (2.9) leads to a convolution equa-

tion which is extremely difficult to write down explicitly,

ii) and even if it is possible to write an inhomogeneous Helmholtz-type equation, solving this

equation for general functions associated with α, β and γ cannot be expected.

Building an accurate boundary condition which approximates the TBC can however be expected

since the pioneering work of Engquist & Majda [9, 10] in the middle of the seventies using the

theory of reflection of singularities at the boundary [14] and pseudodifferential calculus (see for

example [21]).

Let us develop the main ideas. Like in the previous situation with constant coefficients,

we assume that u as been extended by zero for negative times t and that the initial data u0



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 33

and u1 are compactly supported in Ω. Then, Engquist and Majda prove that there exist two

classical pseudodifferential operator a and b of OPS1 such that we get the following Nirenberg-like

factorization [16]

−∂2x + ∂2t + β(t,x)∂t + γ(t,x)∂x + α(t,x) = −(∂x − a(x,t,Dt))(∂x − b(x,t,Dt)) + R, (2.10)

where R is a smoothing operator of OPS−∞. This approximate factorization can be considered

as the extension to the variable coefficients case of the exact form (2.6). Actually, the smoothing

operator R accounts for the fact that the factorization is now true only at high frequencies. The two

pseudodifferential operators a(x,t,Dt) and b(x,t,Dt), with Dt = −i∂t, have respective associated
symbols a(x,t,τ) and b(x,t,τ) of S1 admitting the following asymptotic expansions in homogeneous

symbols

a(x,t,τ) ∼
∑

j≥0

a1−j(x,t,τ) and b(x,t,τ) ∼
∑

j≥0

b1−j(x,t,τ), (2.11)

with classical homogeneous symbols a1−j and b1−j of order 1 − j. This means that we have e.g.
a1−j(x,t,λτ) = λ

1−j
a1−j(x,t,τ), ∀λ > 0. Developing the factorization (2.10), we get

−∂2x + γ(t,x)∂x + ∂2t + β(t,x)∂t + α(t,x) = −∂2x + (a + b)∂x − ab + Op(∂xb) + R (2.12)

since ∂x(bu) = Op(∂xb)u + b∂xu. In the above equation, we designate by Op(σ) the pseudodiffer-

ential operator with symbol σ. If it is possible to compute a and b then, it can be proved that the

TBC for equation (2.9) is given by

(∂
n
− b(x,t,Dt))u = 0, at [0; +∞[×Γ. (2.13)

Indeed, the results in [14] imply that (2.13) annihilates the wave reflected back in the computational

domain. Generally speaking, this TBC, which extends (2.8), cannot be directly implemented since

b is given by an infinite expansion, but it can however be approximated by a k-th order artificial

boundary condition by truncating the series (2.11) to the first k symbolic terms and considering

(∂
n
−
k−1∑

j=0

b1−j(x,t,Dt))u = 0, at [0; +∞[×Γ. (2.14)

The computation of the terms {b1−j}k−1j=0 is therefore needed. To this end, we identify the operators
on both sides of equality (2.12) and we obtain at the symbolic level the system

{
a + b = γ

−a#b + ∂xb = −τ2 + iβτ + α,
(2.15)

which can also be rewritten as

b#b − γb + ∂xb = −τ2 + iβτ + α, (2.16)

by eliminating a. In the above notations, a#b designates the symbol of the composition operator

ab which admits the following expansion (see for example [21]):

b#b ∼
∞∑

m=0

(−i)m
m!

∂
m
τ b ∂

m
t b. (2.17)



34 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

Since b is developable in terms of homogeneous symbols using relation (2.11), we can extract each

term b1−j from (2.16) by identifying the decaying order terms starting from 2 to 2 −j using (2.17).
Beginning with b1, we get that

b1(x,t,τ) = iτ, (2.18)

fixing also the uniqueness of the expansion (2.11). Computing the two next terms gives





b0 =
β + γ

2
,

b−1 =
(γ

2 − β2) + 4α − 2(∂x + ∂t)(γ + β)
8iτ

,

(2.19)

at the right point xr. It directly gives the following proposition.

Proposition 2.1. Let u be the solution to the generalized wave equation (2.9) with C∞ variable
coefficients α, β and γ. Then, the artificial boundary conditions of order k, for k = 1, 2, 3, are
respectively given at the right endpoint xr by

∂xu − ∂tu = 0, at [0; +∞[×{xr},
∂xu − ∂tu −

β + γ

2
u = 0, at [0; +∞[×{xr},

∂xu − ∂tu −
β + γ

2
u − (γ

2 − β2) + 4α − 2(∂x + ∂t)(γ + β)
8

Itu = 0, at [0; +∞[×{xr},

(2.20)

where It is defined by Itu(t) =
∫ t

0

u(s)ds. Similar formulas can be derived at xℓ.

2.3 Short-time vs long-time behavior

Let us now consider that we wish to compute a numerical solution to the problem





(∂
2

t + β(t,x)∂t − ∂2x + γ(t,x)∂x + α(t,x))u = 0, (t,x) ∈]0; T [×Ω,
u(0,x) = u0(x), x ∈ Ω,
∂tu(0,x) = u1(x), x ∈ Ω,

(∂
n
−
k−1∑

j=0

b
p
1−j(x,t,Dt))u = 0, at [0; T ] × {xp},

(2.21)

for a maximal time of computation T and where the k-th artificial boundary condition is defined

by operators {bp
1−j}

k−1
j=0 , for p = ℓ,r, at the left or right fictitious point xp (see e.g. Proposition

2.1). Introducing N intervals of discretization in time, we denote by ∆t the time step defined

by ∆t = T/N. We next seek to compute an approximate solution un(x) ≈ u(tn,x) to system
(2.21), with tn = n∆t, for n ∈ {1, ...,N}. We have seen before that the derivation of the artificial
boundary conditions at the continuous level is made under the high frequency assumption |τ| ≫
1. Since we consider discrete times tn, for n = 1, ...,N, discrete time frequencies τn = π/tn are

then associated and lie in the interval [π/T ; π/∆t]. Hence, the artificial boundary conditions which

work well for high frequencies will be accurate if T ≪ 1. This means that an artificial boundary



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 35

condition is accurate as long as it is used for small times of computation. They may fail for large

computational times which is a known problem of artificial boundary conditions techniques (see

e.g. [8] [13]).

As an illustration, we compare in Figure 1 the performances of the artificial boundary con-

ditions of order 1, 2 and 3 in the case of the wave equation (∂2t + ∂t − ∂2x)u = 0. We give the
relative error in the L2(Ω)-norm for times between 0 and 10. As predicted by the theory, we notice

indeed an improvement for small times by increasing the order. The second order condition is

more efficient than the first order condition for all computed times but the third order condition

is more efficient than the second order condition only for t ≤ 5.2.

0 1 2 3 4 5 6 7 8 9 10
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Figure 1: (∂2t + ∂t − ∂2x)u = 0. Relative error in L2 norm in function of time. abc −− order 1, −−
order 2 and order 3 · − ·.

3 Different approaches for nonlinear equations: the case of

nonlinear Schrödinger equations

3.1 Nonlinear and linear Schrödinger equations

We have seen in Section 2 that it is possible to build accurate artificial boundary conditions using

techniques based on pseudodifferential calculus in the model case of the linear wave equation



36 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

with variable coefficients. The aim of the present Section is to develop some applications of the

pseudodifferential calculus method and its nonlinear version, the paradifferential calculus strategy,

to obtain accurate artificial boundary conditions for nonlinear equations. As a model equation, we

consider the time dependent nonlinear Schrödinger equation

{
i∂tu + ∂

2

xu + α(u)∂xu + β(u)u = 0, (t,x) ∈ [0; +∞[×R,
u(0,x) = u0(x), x ∈ R,

(3.1)

where we assume again that the initial condition u0 has compact support in Ω and that α and β

are two C∞ functions.
Let us now consider the following associated variable coefficients linear Schrödinger equation

{
i∂tu + ∂

2

xu + A(t,x)∂xu + B(t,x)u = 0, (t,x) ∈ [0; +∞[×R,
u(0,x) = u0(x), x ∈ R.

(3.2)

Extending the previous strategy presented for the variable coefficients wave equation in section 2.2

to equation (3.2), one can prove that there exist two pseudodifferential operators a(x,t,Dt) and

b(x,t,Dt) such that we have

∂
2

x + i∂t + A∂x + B = (∂x − a(x,t,Dt))(∂x − b(x,t,Dt)) + R, (3.3)

where again R ∈ OPS−∞. The operators a and b are elements of OPS1/2 admitting the following
expansion in homogeneous symbols

a(x,t,τ) ∼
∞∑

j=0

a(1−j)/2(x,t,τ) and b(x,t,τ) ∼
∞∑

j=0

b(1−j)/2(x,t,τ), (3.4)

where a(1−j)/2 and b(1−j)/2 are homogeneous symbols of order (1 − j)/2. This means that ∀λ > 0,
we have:

a(1−j)/2(x,t,λτ) = λ
(1−j)/2

a(1−j)/2(x,t,τ), b(1−j)/2(x,t,λτ) = λ
(1−j)/2

b(1−j)/2(x,t,τ). (3.5)

If one considers e.g. the right fictitious point xr, then by fixing b1/2(x,t,τ) = −
√
τ, where

√
τ is

defined by:

√
τ =






√
τ if τ ≥ 0,

−i
√
−τ if τ < 0,

(3.6)

it can be shown (see [17, 20]) that the TBC is given by

(∂
n
− b(x,t,Dt))u = 0, at [0; +∞[×{xr}, (3.7)

and that an approximate artificial boundary condition of order k is

(∂
n
−
k−1∑

j=0

b(1−j)/2(x,t,Dt))u = 0, at [0; +∞[×{xr}, (3.8)



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 37

with the convention that the artificial boundary condition of order zero corresponds to the Neumann

boundary condition. The required inhomogeneous symbols can be obtained by adapting relation

(2.16)

b#b + Ab + ∂xb = τ − B, (3.9)

using the suitable substitutions. Then, using the Leibniz symbolic rule we get the four first symbols

b1/2 = −
√
τ,b0 = −

A

2
,b−1/2 = −

1

2
√
τ

(
A

2

4
− B + ∂xA

2
),

b−1 = −
1

8τ
(A∂xA + ∂

2

xA − 2∂xB + i∂tA).
(3.10)

To explain the different strategies which can be considered, we propose now to investigate first

the case α = 0 and next to detail the situation when β = 0, where α and β are the functions in

(3.1).

3.2 Case I: α = 0

3.2.1 Potential strategy

The point of view adopted in this strategy considers that α(u) and β(u) act as potential functions

independent of u. More specifically, they have respectively corresponding functions A and B in

Equation (3.2). If we assume that A = 0, then, the symbols in (3.10) simplify as

b1/2 = −
√
τ,b0 = 0,b−1/2 =

B

2
√
τ
,b−1 =

∂xB

4τ
. (3.11)

Using the definition (3.8), we obtain the following artificial boundary conditions of order k at

[0; +∞[×{xr} 



∂
n
u + e

−i π
4 ∂

1/2
t u = 0, for k = 1, 2,

∂
n
u + e

−i π
4 ∂

1/2
t u − ei

π
4
B

2
I
1/2
t u = 0, for k = 3,

(3.12)

and finally

∂
n
u + e

−i π
4 ∂

1/2
t u − ei

π
4
B

2
I
1/2
t u − i

∂
n
B

4
Itu = 0, for k = 4. (3.13)

In the above equations, the fractional half-order derivative operator ∂
1/2
t , with symbol

√
−iτ, is

given by

∂
1/2
t f(t) =

1√
π
∂t

∫ t

0

f(s)√
t − s

ds, (3.14)

and the half-order integration operator I
1/2
t (with symbol (−iτ)−1/2) is defined by

I
1/2
t f(t) =

1√
π

∫ t

0

f(s)√
t − s

ds. (3.15)



38 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

Following our strategy, we replace B by the nonlinearity β(u) which gives the three following

artificial boundary conditions of order k





∂
n
u + e

−i π
4 ∂

1/2
t u = 0, for k = 1, 2,

∂
n
u + e

−i π
4 ∂

1/2
t u − ei

π
4
β(u)

2
I
1/2
t u = 0, for k = 3,

∂
n
u + e

−i π
4 ∂

1/2
t u − ei

π
4
β(u)

2
I
1/2
t u − i

∂
n
β(u)

4
Itu = 0, for k = 4.

(3.16)

These conditions will be denoted by ABC
β
1,k in the sequel of the paper.

3.2.2 Gauge change strategy

Let us remark that the artificial boundary condition (3.13) is not a transparent boundary condition

even when B is a constant potential. Now, in the case of a time-dependent potential B(x,t) = B(t),

one can get the transparent boundary condition by using the gauge change

v(x,t) = e
−iB(t)

u(x,t), (3.17)

where B(t) = ItB(t), and noticing that v is now solution to the free Schrödinger equation

i∂tv + ∂
2

xv = 0. (3.18)

Then, the transparent boundary condition

∂
n
v + e

−i π
4 ∂

1/2
t v = 0 (3.19)

holds for v and, coming back to the initial unknown u, we obtain the transparent boundary

condition for u

∂
n
u + e

−i π
4 e
iB(t)

∂
1/2
t (e

−iB(t)
u(x,t)) = 0. (3.20)

This boundary condition is clearly not exact if B also depends on x. Nevertheless, we can use a

similar change of gauge, that is

v(x,t) = e
−iB(t,x)

u(x,t), (3.21)

with B(t,x) = ItB(t,x). Then, v is sought as the solution to the variable coefficients Schrödinger
equation

i∂tv + ∂
2

xv + (2i∂xB)∂xv + (i∂2xB − (∂xB)2)v = 0, (3.22)
which is of the general form (3.2) with initial condition v(x, 0) = u0(x). We can therefore apply

the previous general derivation of artificial boundary conditions of Section 3.1 to this equation

of unknown v for suitably defined variable coefficients. As a consequence, if
{
b(1−j)/2

}
j≥0

desig-

nates the symbolic asymptotic expansion of the transparent boundary condition associated with v

solution to (3.22), then an artificial boundary condition of order k is given for u as

∂
n
u −

k−1∑

j=0

e
iB
b(1−j)/2(x,t,Dt)(e

−iB
u) − i∂

n
Bu = 0, at [0; +∞[×{xr}. (3.23)



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 39

More precisely, using (3.10), the computation of the first four symbols gives

b1 = −
√
τ, b0 = −i∂nB, b−1/2 = 0, b−1 =

∂
n
B

4τ
. (3.24)

Finally, we obtain [3] the second- and fourth-order artificial boundary conditions given respectively

by

∂
n
u + e

−i π
4 e
iB
∂

1/2
t (e

−iB
u) = 0, at [0; +∞[×{xr}. (3.25)

and

∂
n
u + e

−i π
4 e
iB
∂

1/2
t (e

−iB
u) − i∂nB

4
e
iB
It(e

−iB
u) = 0, at [0; +∞[×{xr}. (3.26)

Replacing B by β(u), the associated nonlinear artificial boundary conditions are then given by

∂
n
u + e

−i π
4 e
iB(u)

∂
1/2
t (e

−iB(u)
u) = 0, at [0; +∞[×{xr}. (3.27)

and

∂
n
u + e

−i π
4 e
iB(u)

∂
1/2
t (e

−iB(u)
u) − i∂nβ(u)

4
e
iB(u)

It(e
−iB(u)

u) = 0, at [0; +∞[×{xr}, (3.28)

setting B(u)(x,t) = It(β(u))(x,t). These conditions will be referred to as ABC
β
2,j in the sequel,

for j = 2, 4.

Let us develop the connection existing between the artificial boundary conditions ABC
β
m,j, for

m = 1, 2 and j = 2, 4. To this end, let us recall the following Leibniz formula for computing the

fractional derivative of the product of two functions

∂
p
t (fg) =

∞∑

k=0

Γs(p + 1)

k!Γs(p − k + 1)
∂
k
t f∂

p−k
t g, (3.29)

for p > 0. The real-valued function f is supposed to be C∞ and g is a continuous function. The
notation Γs designates the Gamma special function. For p = 1/2, we obtain

∂
1/2
t (fg) = f∂

1/2
t g +

1

2
∂tfI

1/2
t g + R, (3.30)

where R is an error operator in OPS−3/2. Using a similar formula for the integral operator gives

It(fg) = fItg + S (with S ∈ OPS−2). Using these two relations to approximate the half-order
operator in (3.26) by setting f = e−iB and g = u, we see that (3.26) exactly corresponds to (3.13)

up to an operator in OPS−3/2. This error may be not negligible since is involves time derivatives

of the potential, and, in the nonlinear case, of β(u). This difference can therefore be significant

between the two kinds of artificial boundary conditions. This will be more deeply investigated

during the numerical simulations.



40 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

3.3 Case II: β = 0

3.3.1 Potential strategy

Let us now consider the second case where β = 0. Then, the potential strategy consists in replacing

B = 0 in (3.10) leading to

b1/2 = −
√
τ,b0 = −

A

2
,b−1/2 = −

1

2
√
τ

(
A

2

4
+
∂xA

2
),

b−1 = −
1

8τ
(A∂xA + ∂

2

xA + i∂tA).

(3.31)

This gives the following artificial boundary conditions





∂xu + e
−i π

4 ∂
1/2
t u = 0,

∂xu + e
−i π

4 ∂
1/2
t u +

A

2
u = 0,

∂xu + e
−i π

4 ∂
1/2
t u +

A

2
u +

e
i π
4

2
(
A

2

4
+
∂xA

2
)I

1/2
t u = 0,

(3.32)

and

∂xu + e
−i π

4 ∂
1/2
t u +

A

2
u +

e
i π
4

2
(
A

2

4
+
∂xA

2
)I

1/2
t u +

i

8
(A∂xA + ∂

2

xA + i∂tA)Itu = 0, (3.33)

at [0; +∞[×{xr}. Again, replacing A by α(u) yields the nonlinear artificial boundary conditions
of order k






∂xu + e
−i π

4 ∂
1/2
t u = 0, for k = 1,

∂xu + e
−i π

4 ∂
1/2
t u +

α(u)

2
u = 0, for k = 2,

∂xu + e
−i π

4 ∂
1/2
t u +

α(u)

2
u +

e
i π
4

2
(
α(u)

2

4
+
∂xα(u)

2
)I

1/2
t u = 0, for k = 3,

∂xu + e
−i π

4 ∂
1/2
t u +

α(u)

2
u +

e
i π
4

2
(
α(u)

2

4
+
∂xα(u)

2
)I

1/2
t u

+
i

8
(α(u)∂xα(u) + ∂

2

xα(u) + i∂tα(u))Itu = 0, for k = 4.

(3.34)

The set of above j-th order artificial boundary conditions will be called ABCα
1,j is the sequel, for

j = 1, ..., 4.

3.3.2 Linearization strategy

We linearize equation (3.1) with β = 0 around a mean state u and call v its linearization. We

obtain

i∂tv + ∂
2

xv + α(u)∂xv + α
′
(u)∂xuv = 0, (3.35)

which is of the form (3.2) with A = α(u) and B = α′(u)∂xu. Equations (3.10) give

b1/2 = −
√
τ,b0 = −

α(u)

2
,b−1/2 = −

1

2
√
τ

(
α(u)

2

4
− α

′
(u)∂xu

2

)
,

b−1 = −
1

8τ
(α(u)∂xα(u) − ∂2xα(u) + i∂tα(u)).

(3.36)



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 41

This yields absorbing boundary conditions for the linearized problem. To go back to the original

problem, we now need to "unlinearize" these boundary conditions.

The first order absorbing boundary condition for v does not involve u and we immediately

obtain for u

∂xu + e
−i π

4 ∂
1/2
t u = 0, for k = 1. (3.37)

The second order absorbing boundary condition for v reads

∂xv + e
−i π

4 ∂
1/2
t v +

1

2
α(u)v = 0, (3.38)

where α(u)v is the linearization of γ(u), where γ is the primitive of α vanishing at 0. Thus, the

unlinearization of (3.38) is:

∂xu + e
−i π

4 ∂
1/2
t u +

1

2
γ(u) = 0, for k = 2. (3.39)

The unlinearization of the third and fourth order absorbing boundary conditions of v are far more

challenging. We have to unlinearize:

∂xv + e
−i π

4 ∂
1/2
t v +

1

2
α(u)v +

e
i π
4 α(u)

2

8
I
1/2
t (v) −

e
i π
4 α

′
(u)∂xu

4
I
1/2
t (v) = 0, (3.40)

and

∂xv + e
−i π

4 ∂
1/2
t v +

1

2
α(u)v +

e
i π
4 α(u)

2

8
I
1/2
t (v) −

e
i π
4 α

′
(u)∂xu

4
I
1/2
t (v)

+
i(α(u)∂xα(u) − ∂2xα(u) + i∂tα(u))

8
It(v).

(3.41)

To achieve this goal, we rely on the paradifferential operators of J. M. Bony [6] which are general-

ization of pseudodifferential operators well-suited to nonlinear problems. We refer to [17] [20] for

details about these operators and about the rigorous unlinearization of (3.40) and (3.41). We finally

obtain the following nonlinear artificial boundary conditions of order k for u at [0; +∞[×{xr}:





∂xu + e
−i π

4 ∂
1/2
t u = 0, for k = 1,

∂xu + e
−i π

4 ∂
1/2
t u +

γ(u)

2
= 0, for k = 2,

∂xu + e
−i π

4 ∂
1/2
t u +

γ(u)

4
− e

i π
4

4
I
1/2
t (α(u)∂xu) = 0, for k = 3,

∂xu + e
−i π

4 ∂
1/2
t u −

e
i π
4

2
I
1/2
t (α(u)∂xu)

+
i

8
It

(
−α′(u)(∂xu)2 + α(u)2∂xu

)
= 0, for k = 4,

(3.42)

where γ is the primitive of α vanishing at 0. From now, these j-th order artificial boundary

conditions will be referred to as ABCα
2,j, for j = 1, ..., 4.

Remark 1. The unlinearization step of the linearization strategy has been successful not only in
the case of (3.1) with β = 0, but also in the case of the semilinear wave equation with various
nonlinearities (see [19]). However, the unlinearization step of the linearization strategy is not
always successful, as shown by the case of the cubic nonlinear Schrödinger equation (see [20]).



42 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

4 Numerical examples

We consider the model Schrödinger equation (3.1) for the two specific cases I (α = 0) and II (β = 0)

of section 3. More specifically, we choose to present results in Case I for the nonlinearity β(u) = |u|2,
which corresponds to the well-known one-dimensional cubic nonlinear Schrödinger equation, and

in Case II for α(u) = u. We therefore will focus on systems (4.1) and (4.2) respectively given by

{
i∂tu + ∂

2

xu + |u|2u = 0, (t,x) ∈ [0; +∞[×R,
u(0,x) = u0(x), x ∈ R,

(4.1)

and {
i∂tu + ∂

2

xu + u∂xu = 0, (t,x) ∈ [0; +∞[×R,
u(0,x) = u0(x), x ∈ R.

(4.2)

In both cases, we have to our disposal explicit solutions. Concerning system (4.1), we consider the

so-called soliton solution computed by using the inverse scattering theory and given by

uex,α=0(x,t) =
√

2a sech(
√
a(x − ct)) exp(ic

2
(x − ct)) exp(i(a + c

2

4
)t). (4.3)

Concerning system (4.2), adapting the Cole-Hopf transform [20] , we have the explicit solution

uex,β=0(x,t) =

∫
2

0

exp

(
i
(x − y)2

4t

)
u0(y) exp(

∫ y

0

u0(s)

2
ds)dy

×
(√

iπt −
∫ x

0

exp(i
y
2

4t
)dy +

∫
2

0

exp(i
(x − y)2

4t
) exp(

∫ y

0

u0(s)

2
ds)dy

+ exp(

∫
2

0

u0(s)

2
ds)(

√
iπt −

∫
2−x

0

exp(i
y
2

4t
)dy)

)−1
,

(4.4)

which has a compact support in [0, 2] at time t = 0.

In the two situations, we have to solve a nonlinear equation coupled with nonlinear boundary

conditions. The Schrödinger equations are discretized at time tn+1/2 = (tn+1 + tn)/2 by a second-

order approximation. In the sequel, if δt designates the time step, then tn = nδt stands for the

n-th time step, where n ∈ N. The Crank-Nicolson schemes are adapted from the one proposed by
Durán and Sanz-Serna [7] and are given by

i
u
n+1 − un
δt

+ ∂
2

x

u
n+1

+ u
n

2
+

∣∣∣∣
u
n+1

+ u
n

2

∣∣∣∣
2
u
n+1

+ u
n

2
= 0 (4.5)

and

i
u
n+1 − un
δt

+ ∂
2

x

u
n+1

+ u
n

2
+
u
n+1

+ u
n

2
∂x

(
u
n+1

+ u
n

2

)
= 0 (4.6)

respectively for Cases I and II. We denote here by un the approximate value of u at time tn. In

order to reduce the computational time, we set 2vn+1 = un+1 +un, and the schemes read for n ≥ 0

2i
v
n+1 − un
δt

+ ∂
2

xv
n+1

+ |vn+1|2vn+1 = 0, (4.7)



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 43

(a) Case I, a = 2, c = 15, θ = 0 (b) Case II

Figure 2: Exact solutions representations.

and

2i
v
n+1 − un
δt

+ ∂
2

xv
n+1

+ v
n+1

∂xv
n+1

= 0, (4.8)

imposing in both cases that v0 = 0. Clearly, a general form of the previous schemes is

2i
v
n+1 − un
δt

+ ∂
2

xv
n+1

+ v
n+1

f(v
n+1

) = 0, (4.9)

where f designates the map |·|2 or ∂x according to the equation. The Crank-Nicolson approximation
must be coupled to the boundary conditions (3.16), (3.27), (3.28), (3.34) and (3.42). Since the

Jacobian of the maps associated to these nonlinear problems is difficult to obtain, we choose to use

a classical fixed-point method based on the semi-discrete Crank-Nicolson schemes. The choice of a

variational approximation method is thus obvious. Here, we specifically use a P1 linear Lagrange

finite element approximation. The bounded computational domain is the open set Ω =]xl,xr[. The

fictitious boundary is limited to the two endpoints Γ = {xl,xr}. At this point, let us note that the
boundary conditions for the case β = 0 have been given explicitly only at the right endpoint. The

left boundary conditions differ. Concerning the potential strategy, the system of equations (3.34)

is transformed on [0; +∞[×{xl} as:





∂xu − e−i
π
4 ∂

1/2
t u = 0, for k = 1,

∂xu − e−i
π
4 ∂

1/2
t u +

α(u)

2
u = 0, for k = 2,

∂xu − e−i
π
4 ∂

1/2
t u +

α(u)

2
u − e

i π
4

2
(
α(u)

2

4
+
∂xα(u)

2
)I

1/2
t u = 0, for k = 3,

∂xu − e−i
π
4 ∂

1/2
t u +

α(u)

2
u − e

i π
4

2
(
α(u)

2

4
+
∂xα(u)

2
)I

1/2
t u

+
i

8
(α(u)∂xα(u) + ∂

2

xα(u) + i∂tα(u))Itu = 0, for k = 4.

(4.10)



44 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

Concerning the linearization strategy, we get on [0; +∞[×{xl}:





∂xu − e−i
π
4 ∂

1/2
t u = 0, for k = 1,

∂xu − e−i
π
4 ∂

1/2
t u +

γ(u)

2
= 0, for k = 2,

∂xu − e−i
π
4 ∂

1/2
t u +

γ(u)

4
+
e
i π
4

4
I
1/2
t (α(u)∂xu) = 0, for k = 3,

∂xu − e−i
π
4 ∂

1/2
t u +

e
i π
4

2
I
1/2
t (α(u)∂xu)

+
i

8
It

(
−α′(u)(∂xu)2 + α(u)2∂xu

)
= 0, for k = 4,

(4.11)

where γ is the primitive of α vanishing at 0.

The boundary conditions are of memory-type and involve half-order fractional derivatives

and integrals. To preserve the second-order approximation and the unconditional stability of the

Crank-Nicolson schemes, the operators ∂
1/2
t and I

1/2
t are approximated through quadrature rules

which are well suited to the Crank-Nicolson schemes. Namely, we choose the quadrature formulas

derived in [1], which read for the sequence of complex values {fn}n∈N approximating {f(tn)}n∈N,

I
1/2
t f(tn) ≈

√
2δt

2

n∑

k=0

αkf
n−k and ∂

1/2
t f(tn) ≈

2√
2δt

n∑

k=0

βkf
n−k

, (4.12)

where (αk)k∈N and (βk)k∈N designate the sequences defined by





(α0,α1,α2,α3,α4,α5, · · · ) =

(
1, 1,

1

2
,
1

2
,
1 · 3
2 · 4,

1 · 3
2 · 4, · · ·

)
,

βk = (−1)kαk, ∀k ≥ 0.

The composition of the approximation of I
1/2
t with itself gives the approximation of It by the trape-

zoidal rule, which is coherent with the underlying Crank Nicolson scheme. Using these quadratures

formulas, the numerical versions of the boundary conditions (3.16), (3.27), (3.28), (3.34) and (3.42)

are discrete convolutions which may be represented by the following formulation

∂
n
v
n+1

+ e
−i π

4

√
2

δt
v
n+1

+ g(v
n+1

,v
n
,v
n−1

, · · · ,v0) = 0,

where g is a function giving information on the construction of the approximation. For example,

the approximation of ABC
β
1,3 is given by

∂
n
v
n+1

+ e
−i π

4

√
2

δt

n+1∑

k=0

βkv
n+1−k − e

i π
4

2
|vn+1|

√
2δt

2

n+1∑

k=0

αk|vn+1−k|vn+1−k = 0.

The other approximations of ABC
α,β
j,k , j = 1, 2, k = 1, 2, 3 can be found in [3] and [20]. The

complete Crank Nicolson scheme with fixed point procedure therefore takes the form given in

Table 1.



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 45

let w0 = un

s = 0

while ‖ws+1 − ws‖L2(Ωi) > ε do
solve the linear boundary-value problem
∫

Ω

2i

δt
w
s+1

ψdx −
∫

Ω

∂xw
s+1

∂xψdx −
∫

Γ

e
−i π

4

√
2

δt
w
s+1

ψdΓ =

−
∫

Ω

f(w
s
)w

s
ψdx +

∫

Ω

2i

δt
u
n
ψdx +

∫

Γ

g
s
ψdΓ,

setting

g
s

= g(w
s
,v
n
,v
n−1

, · · · ,v0)
and ψ is one of the basis functions of the P1 finite element set.

end while

v
n+1

= w
s+1

u
n+1

= 2v
n+1 − un

Table 1: Fixed-point algorithm for solving the nonlinear Schrödinger equation with nonlinear ABC.

We present below some numerical experiments to show the effectiveness of the different bound-

ary conditions. Since we have an exact solution in Cases I and II, we choose to evaluate the schemes

on the solutions with initial data (4.3) and (4.4) respectively. Concerning Case I, the computational

domain is limited to the open set (−10, 10) discretized with 4000 points. The time step is fixed to
δt = 10

−3. In Case II, the finite domain is (0, 2) discretized with 2000 points. The time step δt is

equal to 2 · 10−3. To analyze the accuracy of the different boundary conditions, we compute the
relative error for the L2(Ω)-norm

‖uex − unum‖0,Ω(t)
‖uex‖0,Ω(t = 0)

,

where unum denotes the numerical solution.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

0.016

0.018

0.02

time

L
2

e
r
r
o
r

n
o
r
m

ABC
β
1,1

ABC
β
1,3

ABC
β
1,4

ABC
β
2,2

ABC
β
2,4

Figure 3: Evolution of the relative error for Case I, a = 2, c = 15, θ = 0.



46 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

For case I, we compare in Figure 3 the ABCs obtained with the potential strategy and the gauge

change strategy for various orders. We can see that increasing the order improves the accuracy.

Moreover, the gauge change strategy (ABC
β
2,k) provides better accuracy compared to the potential

strategy (ABC
β
1,k). Finally, let us also notice that the long-time behaviour of the various ABCs

seems correct. To analyze the accuracy behaviour of the computed solution with respect to the

velocity parameter, we plot on Fig. 4 the evolution of the relative error with respect to c for the

most accurate ABC: ABC
β
2,2. We see that the relative error increases with lower velocities. This is

in agreement with the theory developed in the paper since the boundary conditions are constructed

under a high-frequency assumption.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0

0.005

0.01

0.015

0.02

0.025

0.03

0.035
c=5
c=8
c=12
c=15

time

L
2

e
r
r
o
r

n
o
r
m

Figure 4: Evolution of the relative error for the simulation of the soliton solution with respect to
the velocity c for ABC

β
2,2.

In our last experiment, we focus on case II. We compare in Figure 5 the ABCs obtained with

the potential strategy and the linearization strategy for various orders. Generally, the linearization

strategy leads to the most accurate solutions. Increasing the order of the ABC improves the

accuracy for small times (t ≤ 2.5 in the experiments). After this time, the relative errors cross
and the best results are obtained for ABCα

2,2. The potential strategy is accurate and competitive

for ABCα
1,3 but only for sufficiently large computational times. This shows that each strategy has

his own strengths and weaknesses. Finally, let us mention that ABCα
2,2 has been shown to give

optimal results within a large class of ABCs (see [18]).

5 Conclusion

We presented an analysis of the construction and some numerical validations of ABCs for nonlinear

PDEs considering the example of the nonlinear Schrödinger equation. The methods are mainly

based on pseudo- and paradifferential operator techniques. We show that each strategy can lead

to powerful solutions. However, much work remains to be done. In particular, developing rigorous

extensions for higher dimensions, coupled problems and systems is a complete open problem.



CUBO
11, 4 (2009)

Towards accurate artificial boundary conditions for ... 47

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0

0.005

0.01

0.015

0.02

0.025

time

L
2

e
r
r
o
r

n
o
r
m

ABC
α
1,1

ABC
α
1,2

ABC
α
1,3

ABC
α
2,2

ABC
α
2,3

ABC
α
2,4

Figure 5: Relative error for Case II.

Received: April 2008. Revised: July 2008.

References

[1] Antoine, X. and Besse, C., Unconditionally stable discretization schemes of non-reflecting

boundary conditions for the one-dimensional Schrödinger equation, J. Comput. Phys. 181 (1)

(2003), pp. 157-175.

[2] Antoine, X. Arnold, A. Besse, C. Ehrhardt, M. and Schädle, A., A review of trans-

parent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equa-

tions, Communications in Computational Physics 4 (4) (2008), pp. 729-796 (Open access online

article).

[3] Antoine, X. Besse, C. and Descombes, S., Artificial boundary conditions for one-

dimensional cubic nonlinear Schrödinger equations, SIAM J. Numer. Anal. 43 (6), (2006),

pp.2272-2293.

[4] Bayliss, A. and Turkel, E., Radiation boundary-conditions for wave-like equations, Comm.

Pure Appl. Math. 33 (6), (1980), pp.707-725.

[5] Bérenger, J.P., A perfectly matched layer for the absorption of electromagnetic-waves, J.

Comp. Phys. 114 (2), (1994), pp. 185-200.

[6] Bony, J. M., Calcul symbolique et propagation des singularités pour les équations aux dérivées

partielles non linéaires, Ann. Sci. Ec. Norm. Sup (4 ème série). 14, (1981), pp.209-246.

[7] Durán, A. and Sanz-Serna, J. M., The numerical integration of relative equilibrium solu-

tions. The nonlinear Schrödinger equation, IMA J. Numer. Anal. 20 (2) (2000), pp. 235-261.



48 Xavier Antoine, Christophe Besse and Jérémie Szeftel CUBO
11, 4 (2009)

[8] Engquist, B. and Halpern, L., Long-time behavior of absorbing boundary conditions,

Math. Meth. Appl. Sci. 13, (1990), pp.189-203.

[9] Engquist, B. and Majda, A., Absorbing boundary conditions for the numerical simulation

of waves, Math. Comp. 31, (1977), pp.629-651.

[10] Engquist, B. and Majda, A., Radiation boundary conditions for acoustic and elastic wave

calculations, Comm. Pure Appl. Math. 32, (1979), pp.313-357.

[11] Givoli, D., Nonreflecting boundary conditions, J. Comp. Phys. 94 (1), (1991), pp.1-29.

[12] Hagstrom, T., Radiation boundary conditions for the numerical simulation of waves, Acta

Numerica, (1999), pp. 47-106.

[13] Hagstrom, T. Hariharan, S.I. and MacCamy, R.C., On the accurate long-time solu-

tion of the wave equation in exterior domains: asymptotic expansions and corrected boundary

conditions, Math. Comp. 63, (1994), pp.507-539.

[14] Majda, A. and Osher, S., Reflection of singularities at the boundary, Comm. Pure Appl.

Math. 28, (1975), pp.479-499.

[15] Mur, G., Absorbing boundary-conditions for the finite-difference approximation of the time-

domain electromagnetic-field equations, IEEE Trans. Electromagnetic Compatibility 23 (4),

(1981), pp. 377-382.

[16] Nirenberg, L., Pseudodifferential operators and some applications, CBMS Regional Conf.

Ser. in Math. AMS 17 (1973), Lectures on Linear Partial Differential Equations, pp.19-58.

[17] Szeftel, J., Calcul pseudodifférentiel et paradifférentiel pour l’étude de conditions aux limites

absorbantes et de propriétés qualitatives d’équations aux dérivées partielles non linéaires, Ph.D.

thesis, Université Paris 13, 2004.

[18] Szeftel, J., Absorbing boundary conditions for nonlinear scalar partial differential equations,

Comput. Methods Appl. Mech. Engrg. 195, (2006), pp.3760-3775.

[19] Szeftel, J., A nonlinear approach to absorbing boundary conditions for the semilinear wave

equation, Math. Comp. 75, (2006), pp.565-594.

[20] Szeftel, J., Absorbing boundary conditions for nonlinear Schrödinger equations, Numer.

Math. 104, (2006), pp.103-127.

[21] Taylor, M. E., Pseudodifferential Operators, Princeton University Press, NJ, 1981.

[22] Tsynkov, S.V., Numerical solution of problems on unbounded domains. A review, Appl.

Numer. Math. 27 (4), (1998), pp. 465-532.

[23] Zheng, C., Exact nonreflecting boundary conditions for one-dimensional cubic nonlinear

Schrödinger equations, J. Comput. Phys., 215, (2006), pp.552-565.


	Articulo 3