A New Theoretical Approach to Absorbing Layers.
We consider the construction of absorbing boundary layers using asymptotic expansions, the small parameter being the width of the layer. This allows us to define the order of the layer. We compute layers of various orders and analyze their stability properties. Finally, we perform some numerical tests to evaluate the practical utility of our layers.
Introduction
Classically when one wants to simulate numerically the propagation of waves in
an unbounded domain in space, there are (at least) two methods for limiting the
computational domain. One is to use absorbing boundary conditions which consist in
prescribing boundary conditions specifically derived in order to minimize the artificial
reflections inside the computational domain. Such conditions are now very well known
and have been studied intensively from the physical, mathematical and numerical
point of view in the last 15 years; see for instance [1], [3], [4], [7], [15] and references
therein. A technically simpler and older approach is to use absorbing layers. This
is equivalent to surrounding the domain of interest by some geometrical layer where
one adds physical attenuation terms to the propagation equations in order to absorb
the waves penetrating this region. Of course to get good absorption the thickness of
the absorbing layer must be sufficiently large and to avoid increasing too much the
size of the computational domain and therefore the computational cost it must be
sufficiently small. Moreover, the transition between the propagation medium and the
absorbing region should not be too abrupt in order to allow the waves to penetrate
the absorbing layer.
The design of such absorbing layers being essentially based on physical or engineering
considerations, it is not surprising that most of the articles concerning this
subject can be found in the physical literature, see for instance [2], [8], [9], [10] and
[13]. On the other hand, this domain has not retained too much the attention of the
applied mathematicians (perhaps because of the development of the more attractive
1991 Mathematics Subject Classification. primary 65M99, 35C20 secondary 65M06, 35L05.
Key words and phrases. absorbing layers, absorbing boundary conditions, wave equation.
absorbing boundary conditions). As far as we know, the only paper of a mathematical
nature treating the problem is by Israeli and Orszag, [8], where one can find
a complete presentation of absorbing layers in the one dimensional case. However,
their article does not really contain a theoretical analysis and in particular it does
not address the following question: is it possible to establish a mathematical link
between absorbing layers and absorbing boundary conditions? In this paper we intend
to develop a tentative answer (or at least an embryo of an answer) for the one
dimensional case to begin with. Our objective is to interpret the absorbing boundary
condition as a limit case of some absorbing layer when the thickness " of this layer
tends to zero. In fact, this approach leads to mathematical criteria for the design
of the absorbing layer, especially when one wants to improve the accuracy of the
approximation with respect to the parameter ".
There exists a nice formal solution to this that we shall present in section 2. It is
based on an asymptotic analysis, the starting point of which is the generalized image
principle for absorbing boundary conditions recently introduced by Ha Duong and
Joly [5]. However, theoretical problems occur when dealing with the justification of
the asymptotic expansions we use, and more precisely with the well-posedness and
the stability of the approximate problems. We shall evoke these problems in section
3 and we shall only give partial and not completely satisfactory solutions to them.
In section 4 we will prove the existence of absorbing layers of any order and compute
some examples. Finally in section 5 we will present various numerical simulations
illustrating on one hand the pertinence of our approach and on the other hand the
difficulties of its implementation.
Although we think that the interest of this paper is a priori essentially theoretical,
we also think that the underlying ideas should be useful for designing numerical
schemes for absorbing boundary conditions when one uses high order finite difference
schemes for the interior propagation equation. This was in any case the motivation
for the introduction of the generalized image principle. Moreover, as this principle
can be generalized to two and three dimensional cases, our absorbing layers could
also be useful in these situations.
To conclude this introduction we would like to emphasize a great difference between
what we are doing here and what is usually done with absorbing layers. When
one wants to improve the absorption properties of such layers it is natural to increase
their size (and consequently the computational cost is increased). The situation here
is exactly the opposite since we approach the exact transparent boundary condition
when the thickness of the absorbing layer tends to zero. From the computational
point of view this is of course very interesting and may appear miraculous. How-
ever, one should not forget that our criteria for designing absorbing layers are purely
mathematical while previous approaches were based on physical considerations. On
the practical side let us finally note that it is the number of points in the absorbing
layer which should be small and not its thickness as such; unfortunately we shall see
that this number must be quite high to obtain a good absorption and so these layers
are not so interesting in practice as we had hoped.

Acknowledgements

We would like to thank Olivier Vacus (INRIA, Rocquen-
court) for his help in carrying out the numerical experiments in the final section.
2. Transparent boundary condition as the limit case of absorbing
layers: an asymptotic approach
2.1. Setting of the problem. We shall consider the very simple one dimensional
model problem which consists of reducing to the half line the
computation of the solution u of the Cauchy problem for the wave equation with
constant velocity equal to one. We also assume that Cauchy data have compact
support in R \Gamma . 8
It is well known that the restriction of u to R \Gamma is characterized as being the unique
solution of the following initial boundary value problem
The boundary condition is the transparent boundary condition. The
alternative solution is to replace u by the solution u a of the new problem
u a
some condition at
In the preceding problem L denotes the thickness of the absorbing layer (which is
here the interval I = [0; L]) and the functions oe and - are positive functions which
model the physical attenuation introduced in the absorbing layer. These functions
are supposed to satisfy: 8 !
Note that except for the fact that we have written the third order term in the conservative
form, the damped wave equation appearing in (2.2) is exactly the same
that one finds in the article by Israeli and Orszag [8]. The damping terms oeu t and
\Gamma(-u xt ) x can be interpreted physically as fluid and viscous damping terms.
The dissipation effect of these terms is easily seen in the following energy identity
(assuming for instance that u x (L;
d
dt
(ju a
Z L
(oe ju a
which proves that the energy
E(u a
Z L
(ju a
x
is a decreasing function of time. The identity (2.3) should be compared to the one
associated to the 'exact' solution u of problem (2.1)
d
dt
We are now naturally confronted with the following questions
(i) How to choose functions oe and -?
(ii) How to choose the thickness L?
(iii) Which boundary condition should be taken at
(iv) In what sense can u a be considered as an approximation to u?
There is a priori no clear answer to these questions and usually the solutions proposed
are based more on engineering know-how than on scientific argument. Recently, Ha
Duong and Joly ([5]) introduced the so-called image principle, with help of which
the solution u of (2.1) can also be considered as the restriction to R \Gamma of the unique
solution ~ u of the following problem
~
2ffi(x)\Omega ~
where ffi denotes the Dirac distribution at
denote the
extensions to R of u 0 and u 1 by evenness, i.e.
~
~
~
appears then as the solution of a wave equation with a punctual (and thus singular)
damping term, namely
ffi(x)\Omega ~
. It can be shown that the solution ~
u has never better
regularity than H 1 in the neighborhood of origin. It is then a natural idea to regularize
~
u by regularizing ffi, which is equivalent to approximating ~
u by the solution ~
" of the
following problem
where the function i " (x) is supposed to have the following properties.
If these conditions are satisfied then it is well known that in the sense of distributions
one has
It is then natural to approximate u, as it is the restriction of ~
u to R \Gamma , by u " defined
to be the restriction to R \Gamma of ~
we have chosen i " to be even and because ~
and ~ u 1 are also even, it is clear that ~
can be characterized as the unique solution of
~
~
~
Note that under the additional hypothesis i - 0 (which we will abandon in the
sequel), the energy identity associated with (2.4) proves that
which gives a priori estimates on u " which are independent of " and allow us to prove
that at least weakly
"&0
In fact it is not difficult to prove that if u 0 and u 1 have compact support we get
strong convergence in the space L 1 (0; T ; L 2 ). Now we realize that solving (2.4) is
equivalent simply to solving the wave equation in (\Gamma1; \Gamma"] with an absorbing layer
0]. Up to the change of variable x coincides with (2.2) if we set
6 PATRICK JOLY AND JUKKA TUOMELA
We have thus given a possible answer to questions (ii) and (iii) (since we take the
Neumann boundary condition) and (iv). However, the answer to question (i) still
remains unknown. The purpose of the next section is to clarify this point.
2.2. Asymptotic analysis of the approximation procedure. Before entering
the details of the analysis we shall restrict ourselves to a specific, but classical, class
of functions i " which are obtained by scaling. First we rewrite the problem (2.4)
after having applied the change of variable x our domain of interest
is (\Gamma1; 0], independent of ", while the absorbing layer is [0; "]. We thus consider
which is the unique solution of the following problem
The function i is supposed to satisfy
Note that we relax the positivity condition for reasons which will become clear later.
We now wish to study the asymptotic behavior of u " when " tends to zero. In order
to do that we first use a change of variable in the strip [0; "] to get equations defined
in a domain independent of ". So let us introduce a new function v " such that
The problem (2.5) is then equivalent to the following transmission problem
find
x (0; t)="
transmission conditions
Note that the transmission conditions express nothing but the continuity of u " and u "
x
at denotes the solution of the original problem (2.5). Because of the terms
with 1=" and 1=" 2 occurring in the equation for v " and in the second transmission
condition, we have to deal with a singular perturbation problem. To analyze it we
make the ansatz that both functions u " and v " admit an expansion in powers of ".
We assume these to be asymptotic expansions in the sense that
Assuming further that this is valid in a sufficiently strong sense the calculations
which follow are rigorous, otherwise they must be regarded as formal. Of course we
are essentially interested in the behavior of u " , that of v " being purely anecdotal.
Indeed, we want to find conditions for i which guarantee that
u 0 is a solution of (2.1)
(2.
where n is as large as possible. According to our asymptotic expansion this is the
same as requiring u us introduce
Definition 2.1. If the expansions (2.7) are such that they have the property (2.8),
then we say that the boundary layer i is of order n.
To have explicit conditions on i we substitute (2.7) into (2.6) and compare coeffi-
cients; this is the purpose of the following section.
2.2.1. Characterization of the zeroth and first order terms. Looking at the O(" \Gamma2 )
term in the equation satisfied by v " we obtain v 0
Taking into account the
boundary condition v 0
using the fact that v 0 (0;
that
Recall that we want u 0 to be a solution of (2.1) with absorbing boundary condition
we have to calculate v 1 . The O(" \Gamma1 ) term of the equation
where we have used the identity v 0
t (0; t). The transmission conditions are
x (0; t) and they can be interpreted as initial
conditions for the above ODE. Now to have a reasonable absorbing layer we of course
require that u 1 is identically zero so that the first condition is simply v 1 (0;
Then direct integration gives
We have not yet used the boundary condition at taking the derivative gives
We have proved the following expected result
Lemma 2.1. The boundary layer is of first order if and only if
This result expresses the fact that
Let us then consider the higher order terms. For v 2 we get the following differential
equation
xx
Recall that v 0 (x;
Introducing the function D 1
s)i(s)ds we obtain
Using the initial conditions v 2 (0;
immediately
ds
Then the boundary condition gives
Then we see that
Lemma 2.2. i is of order two if and only if
This result calls for two important comments.
Remark 2.1. Condition (ii) of the previous lemma forces us to abandon the assumption
of the positivity of the function i because of the following result.
Lemma 2.3. There exists no positive function i satisfying simultaneously conditions
(i) and (ii) of the lemma 2.2.
Proof. Writing down the second condition in detail gives
Z 1i(x)
Now if i is positive then the first term is positive, and if the condition (i) is satisfied
the second term is less than one in absolute value, so the sum of the three terms is
positive.
Remark 2.2. We notice that the condition (ii) is nonlinear with respect to i. Generally
conditions for such functions to be good approximations to ffi are linear and
concern the successive moments of i, that is,
. This is what
occurs in the analysis of the approximation of linear advection equations by particle
methods, see [12]. This is due to the fact that the solution in this case depends linearly
on the Dirac distributions which are regularized and which are in this case the
initial data for the advection equation. Here the situation is different since ffi appears
as a coefficient in the damped wave equation and it is well known that the solution
of a partial differential equation depends nonlinearly on its coefficients, even if the
equation itself is linear.
We now proceed to the generalization to higher order terms.
2.2.2. Characterization of higher order terms. The next result will be established by
induction. Let us make the assumption that
where D k denote some functions to be determined. By previous calculations we
know that this is true when the assumptions (i) and (ii) of lemma 2.2
are satisfied and if D 1 is defined as above. It is clear that u n (like u 1 ) satisfies the
homogeneous wave equation with zero initial data and zero boundary condition so
that it is identically zero. So we just need to compute v n . As before we get the
differential equation
Now, using the induction hypothesis, we get
i(x)D
Comparing with (2.10), we add the convention that D Combined with the
transmission conditions this leads to the following differential equation for D n
with initial conditions D n
\Gamma1. This guarantees that u n is identically zero. Then the boundary condition can
be expressed simply as D 0
computations show that this is
equivalent to Z 1i(x)D
We have thus proved the following result
Theorem 2.1. Suppose that (2.7) holds in a sufficiently strong sense. Define the
sequence of functions D k by the induction
i(s)D
ds
Then only if D 0
or more explicitly if
and only if
R 1'
i(x)D
Moreover the function v " satisfies
Remark 2.3. As the sequence D k depends on i we should have denoted it D k (i; x).
Moreover, the dependence of D k with respect to i is nonlinear and in some sense the
nonlinearity increases with k. Indeed, if we fix i, the function - ! D k (-i; x) is a
polynomial of degree k.
2.3. Generalization to the case - 6= 0. Let us write down the general problem.
find
xt
x
x (0; t)="
transmission conditions
We substitute the same asymptotic expansions to the equations as before. Starting
with the zeroth order term we get the following differential equation.
We easily deduce that v 0 (x; Now proceeding to higher order terms as
before we require that the u n are identically zero when n - 1. For the first order
term we get the following equation.
The term containing - vanishes because v 0
xt is zero. The solution of the above problem
is then
Now imposing the boundary condition v 1
equivalently requiring that
Let us then define the functions d
R x
s)i(s)ds and
make the induction hypothesis v k (x;
Writing
the equation for higher order terms (n - 2) we have
xt
This has the solution
xt
ds
ds
Z x-(s)v
ds
xt (0; t)
i(s)d
ds
Now using the induction hypothesis and the facts that d 0
when
i(s)d
ds
(s)ds
ds
Consequently we define
i(s)d
ds
(s)ds
ds
3. Note that if there would be no need to distinguish the two
cases.
Of course d defined in the previous section and other d n 's
are characterized also by the following differential equations
with initial conditions d n
Then remebering the boundary condition at we arrive at
Theorem 2.2. Suppose that (2.7) holds in a sufficiently strong sense. Define the
sequence of functions d k as above. Then only if the n conditions
are satisfied, or more explicitly if and only if
R 1'
i(x)d
Moreover the function v " satisfies
Corollary 2.1. There are no positive functions i and - such that the conditions
are satisfied.
Proof. Note that d 1 does not depend on -.
14 PATRICK JOLY AND JUKKA TUOMELA
Remark 2.4. Note that if (which is rather natural) the first two conditions
do not depend on -.
3. Some considerations about the stability of the absorbing layers
In this section we shall consider the problem of the mathematical justification of
section 2. Of course, as we said before, the existence and uniqueness of u " for a given
value of " is not a problem and the real question is the justification of the expansion
(2.7). As we deal with a linear problem, this can be reduced to the derivation
of appropriate estimates of u " which are independent of " in appropriate norms.
This is classical in the analysis of numerical schemes: section 2 would correspond to
the consistency and accuracy analysis while this section corresponds to the stability
analysis. Then it suffices to refer to the the well known Lax equivalence theorem:
consistency stability implies convergence. The only case where we know how to
get such estimates easily is the case i(x) - 0 (cf. section 2), which is not really
interesting since it does not allow us to get a better result than
appears that in the case where i has a variable sign, the possibility of getting uniform
estimates of u " is an open and non trivial question and from numerical experiments
one concludes that the answer is positive for some i's and negative for some others
even when one respects the condition
The calculations we are going to develop now are mostly formal and still need
to be justified. Nevertheless we think that they can be considered as a first step
towards the understanding of this complicated stability problem. As our calculations
are formal, we shall not make precise any functional framework.
Let us introduce the operator T " acting on functions OE(t) of the only variable t ? 0
and defined by :
denotes the unique solution of :
By construction, u " is solution of :
In fact we are essentially interested in obtaining uniform estimates for u " and what
happens to v " does not really interest us. Moreover it is clear that the properties of
are linked to those of T " . Indeed we have the energy identity :
d
dt
Therefore it is clear that a sufficient condition for obtaining of uniform estimates for
" is that T " has the following positivity property :
Now we first remark that, since the coefficients in (3.2) do not depend on time, T " is
a causal convolution operator so that we can write, at least formally :
Moreover using the time scaling t ! t=", one easily sees that :
corresponds to " = 1. Now assume that G 1 (t) admits a Fourier transform
here by ! the dual variable of t), then it is clear that -
is solution of :
Remark 3.1. (i) We have chosen as the Fourier transform :
Z
f(t)e \Gammai!t dt
when f is an L 1 function. With this convention the Fourier transform of the
Dirac distribution at is equal to 1.
(ii) We have supposed implicitly that (3.9) is uniquely solvable when ! is real.
This is not so obvious for a general function i(x).
Let us now go back to the operator T " . If we assume that T " maps L 2 (R) into itself
then Plancherel's theorem says that :
Z
Let us apply this result to the function OE T defined by :
Z
Therefore using the causality of the operator
Z
Therefore we get the positivity property for T " if we assume that:
Remark 3.2. (i) Condition (3.14) appears, assuming the validity of our compu-
tations, as a sufficient condition for obtaining of uniform L 2 estimates of u " .
Note that we should have denoted -
it depends on i
in a very non linear way through the solution of problem (3.9).
(ii) It is not difficult to derive the identity :
(just multiply (3.9) by -
V and integrate by parts). This proves in particular
that (3.14) is satisfied if i(x) - 0, which is not astonishing. However the
positivity of i(x) is not necessary for (3.14) to hold.
According to comment (i), our condition (3.14) is not really explicit in terms of
the function i and in practice it can only be verified by numerical calculations. So in
order to get more explicit information, we have chosen to embed our original problem
into a more general class of problems depending on one real parameter -. The idea
is to come back to our original problem where we replace ffi(x) by -ffi(x) for - ? 0
and thus to consider the equation :
~
~
~
This problem is obviously well-posed for any - 0. It is now natural to approximate
its solution u by the solution u
~
~
~
If we wish to analyze the convergence of u " to u along the same lines as in the case
1, we are naturally lead to introduce the time convolution operator T -
" whose
" (t) is given by:
The Fourier transform of G -(t) is given by:
x (0; !)
is the solution of:
x
For given - we obtain a sufficient condition for the stability of the approximation
process by expressing the positivity of the operator T -
" , which gives the condition:
where we also have, according to (3.15):
It is natural to ask that condition (3.19) be satisfied for any - ? 0. Up to now we
have not progressed very far in terms of simplicity except that we can now say that
if (3.19) holds for any -, it is satisfied in particular for small values of -.It is not
difficult to check that -
is the solution of :
x
and thus given by :
Now, from (3.20), we deduce that, when - goes to 0:
so that we deduce from
I
Condition (3.21) appears as a necessary condition for our more general condition
(3.19) to be satisfied. It is a priori not sufficient since it means only that (3.19)
will be satisfied only for small -. However, it has the advantage of being linear and
rather explicit in terms of i and we shall see below that it is sufficient to explain
some instability phenomena we have observed numerically. Before presenting these
examples let us remark :
(i) One checks immediately that (3.21) is automatically satisfied as soon as i is
positive.
(ii) The inequality I i (!) - 0 will be automatically satisfied for large ! as soon as
with respect to ! and:
lim
I
Let us apply this to two particular examples of second order accurate layers in the
sense of lemma 2.2. Looking for a linear function and plugging this
expression in the two equations (i) and (ii) of lemma 2.2, we get a quadratic system
of equations in a and b which can be solved easily. The two solutions and thus the
two distinct layers are given by:
ABSORBING LAYERS 19
In particular we note that none of the layers is positive. Now, if
computes that:
I
from which we are going to show that i 1
d
does satisfy the condition (3.21) while i 1
du
does not. Indeed we note that :
I
du
which proves that i 1
du violates condition (3.21). On the other hand, setting
sin !=!, we can express I i (!) as a polynomial of degree 2 with respect to X
I i (!) =2
We easily check that the discriminant is negative for
d which proves that I i 1
d
is positive for any !. In the final section we will show that in actual numerical
experiments
d gives a stable method while using i 1
du leads to instability.
4. Analysis of the existence of the absorbing layers
4.1. Argument from control theory. We analyse the existence of the absorbing
layers by using nonlinear control theory, see [6], [14]. First of all let us introduce
some notations. Let y k
(x). Then the
differential equations (2.11) can be written as
where f and g are vector fields defined by
y 2n1
Initial conditions are then y n+1
conditions are given by the final condtions y n+k From the
point of view of control theory y is a state and i is control. The final state is not
completely given because y k (1) are arbitrary for In any case we shall
prove that system satisfies everywhere a local controllability criterium.
We will need first some definitions.
Definition 4.1. Given two (smooth) vector fields f 1 and f 2 , their bracket [f
the vector field J(f 2 )f is the Jacobian of f i . A vector space
(of vector fields) with the bracket as a multiplication operator is called Lie algebra. If
f k is a collection of vector fields then the Lie algebra generated by f k is the smallest
Lie algebra which contains f k .
These expressions usually are given on manifolds, but since our state space is just
R 2n that formalism is not needed and the above formula gives globally a vector field
in R 2n . Note that a Lie algebra is non-associative and non-commutative.
Definition 4.2. The system (4.1) is said to have accessibility property at y 2 R 2n
if the set of states which can reached from y in finite time contains an open set. The
system (4.1) is said to satisfy the Lie algebra rank condition at y if the Lie algebra
generated by f and g has rank 2n at y.
This formal property is interesting because of the following result (see [6], [14]).
Proposition 4.1. If a control system satisfies the Lie algebra rank condition at y
then it has accessibility property at y.
Then we can prove
Lemma 4.1. System (4.1) satisfies the Lie algebra rank condition for all y 2 R 2n
and for all n.
Proof. Starting with vector fields f and g we generate the following fields:
ABSORBING LAYERS 21
Inductive argument shows that
f
\Gammay 1
\Gammay n\Gammak0
In f 2k\Gamma1 the element \Gamma1 is in k'th position and y n+1 is
in f 2k we have 1 as 1'th element. Then we just evaluate the determinant
\Gamman . This proves that (4.1)
satisfies the Lie algebra rank condition for any n. Moreover the value of the determinant
does not depend on the point y so the condition is satisfied everywhere.
Corollary 4.1. The system (2.14) formulated as a control problem also satisfies the
Lie algebra rank condition.
Proof. Introducing another control function - means that we have more vector fields,
so the rank cannot diminish.
These results suggest that the systems are controllable, which would imply the existence
of absorbing layers for any n. However, we have not been able to pass from
these local results to global ones.
4.2. Computation of some absorbing layers. To effectively calculate the layers
we proceed as follows (we consider first the case us look for i in the
following
are some known functions and c i are
constants to be determined.
Lemma 4.2. The conditions (2.15) can be written as a polynomial system F
is a polynomial
of degree k. Consequently the equation F (c) = 0 has at most n! (complex) solutions.
Proof. By induction f k is a polynomial of degree k. Then the classical theorem of
B'ezout says that there cannot be more than n! solutions.
22 PATRICK JOLY AND JUKKA TUOMELA
Remark 4.1. The structure of the equations appears to be such that there are exactly
n! solutions, see [11].
Of course in the present case we are only interested in real solutions. It seems to be
difficult to prove directly that the system has some real solutions and we have not
been able to find the proof in the general case. However, numerically the solutions can
be calculated, although the system tends to become more and more ill-conditioned
when n grows. Before giving the actual layers, let us formulate some conjectures
concerning the form of the absorbing layers.
Conjecture 4.1. If
is of order n then c i c
2.
Conjecture 4.2. If
is of order n then it has in (0; 1).
We have analysed the following types of layers. First we have two piecewise constant
layers whose orders are 3 and 4, and the graphs are given in figure 4.1. Note
that even discontinuous layers satisfy the conjeture 4.2; the definition of the zero of
the discontinuous functions being obvious from the picture!
0:5461829583590986 if 2=5
Then we have two layers which have the discontinuity at the origin (d for discontin-
uous). The order is the number of coefficients and the curves are in figure 4.2.
d
105)x
d
Finally there are six layers with
-248

Figure

4.1. Layers
pc and i 2
pc .

Figure

4.2. Layers i 1
d and i 2
d .
s
To obtain the equations for the coefficients we used Mathematica [16]. The equations
get quickly very cumbersome so they cannot be given here. Then we solved the
equations with Newton's method. For higher order layers the condition number of the
Jacobian becomes quite large so for the sixth order layer we represented the layer using
Legendre polynomials instead of the power form. This improved the conditioning,
which might also be guessed from the fact the coefficients of the Legendre polynomials
were all less than ten in absolute value.
The initial guess for the Newton's method is not a priori easy, but in the present
case the successive layers behaved very regularly which allowed to predict a good
starting point. In the pictures 4.3 and 4.4 there are the layers i 2
s to i 6
s and examining
numerically their zeros and the maximum value at allows one to guess the
location of zeros of the next layer and the maximum. With this trick Newton's
method converges without any problem.
-5
Figure

4.3. Layers i 2
s to i 6
s .

Figure

4.4. Layers i 2
s
to i 6
s
26 PATRICK JOLY AND JUKKA TUOMELA
Figure

4.5. Stability curves for
d and i 3
s .
Next to test the layers with - we have calculated the following pairs of functions.
Note that when the first two conditions do not depend on - and so we use
the same i 2
s as above. Then raising the order of the layer we get the - i given above.
The computation of the coefficients proceeds similarly as in the case where
4.3. Stability considerations. Let us give some examples of the stability curves
using the necessary stability condition (3.21). In figure 4.5 we have plotted examples
of typical curves in piecewise constant, discontinuous and smooth case. Clearly the
'stability margin' improves with smoothness, which is in some sense 'natural'. Finally
we would like to formulate some conjectures. All the computed examples satisfy them
and in some sense they seem to be rather natural.
Conjecture 4.3. For all n and l there exists a layer of the form
which is of order n and which satisfies the necessary stability condition.
Conjecture 4.4. Let
be of order n; then it satisfies the necessary
stability condition if and only if c
5. Numerical results
We will present some numerical simulations with the boundary layers derived
above. Let us first describe the general setting. We consider the wave equation
in the interval I = [0; 10] and add the absorbing layer at the right end (the boundary
condition at the other end is irrelevant). Consequently the computational domain
will be I is the length of the layer. As initial conditions we
used
~
We tried also other signals with different smoothness properties, but the results were
so similar that we shall not present them in detail. Of course the reflections were
bigger for irregular signals, but the effect was same for all layers. Also varying the
pulsewidth had similar expected effects.
The initial signal then separates into right and left going components and consequently
the amplitude of the signal which hits the absorbing layer is half of the
amplitude of the initial signal. We discretized the equation with the follwing simple
second order scheme
ffit=h, the time-step divided by the space-step. Evidently
where i is one of the layers calculated above for the unit interval. The stability
condition for the scheme without the absorbing layer is ff - 1, and as is well known,
the most accurate solution is obtained at the stability limit. Of course the layer
changes a little this stability limit; we have not tried to determine the exact condition
in each case, but instead have used the value throughout to be able to compare
the results easily. In this first series of tests we adjust the space discretization step
such that for all " the number of points in the layer is about 50. We measure the
reflexion in four different norms: L 2 , L 1 , H 1 and the energy norm. The latter is
defined by kuk E
and it is in some sense the most natural with
wave propagation problems.
28 PATRICK JOLY AND JUKKA TUOMELA

Figure

5.1. Solution at
du .
Before actual comparison let us first demonstrate that if the necessary stability
condition (3.21) is violated we really get an unstable layer. Let us take the linear
du defined in (3.23). In figure 5.1 we have plotted the solution at
even though the whole signal has not yet entered the layer the unstable reflection is
clearly visible already. After that the instability grows rapidly to the overflow.
5.1. Comparison of different layers. Let us start with the piecewise continuous
layers
pc and i 2
pc given in the previous sections. Recall that the number of points in
the layer is about 50 for all ". The norms of the reflected signal were computed at
when the whole signal had already left the absorbing layer. The results
of simulations are presented in table 1. The norms of the initial signal u i were as
11:37. We see
that for small " the schemes are unstable. We do not know if they are unstable for
all ff ? 0, but at least they were unstable even with quite small ff so that in any case
they would not be practical. Recall that in figure 4.5, it was seen that the necessary
stability condition was satisfied with very small margin, so that in some sense one
'expected' them to be unstable. Note, however, that fourth order method is clearly
better than the third order method for bigger ".
So we have to abandon the piecewise constant layers. Then we demonstrate that
the discontinuity at the origin, although creating no problems of stability, is bad for
the accuracy. Consider the results given in table 2. The orders are the same, that
is
d and i 2
s are of the second order and i 2
d and i 3
s are of the third order. It is clear
that the continuous layers are better.
Now let us give the corresponding results for layers i 4
s and i 6
s . These are
presented in table 3. We see that the absorption improves with increasing order.
However, higher order layers tend to be more unstable. Recall that we used the
0:1 91:753 7:3
0:2 0:0358 0:4498 0:0263 0:3190
pc 0:3 0:1080 0:8975 0:0719 0:6387
0:4 0:2213 1:6186 0:1449 1:1537
0:5 0:3543 2:5733 0:2426 1:8357
0:2 0:0296 2:0367 0:0085 1:5075
pc 0:3 0:0425 0:4542 0:0236 0:3220
0:4 0:0800 0:7561 0:0539 0:5373
0:5 0:1390 1:2481 0:1024 0:8884

Table

1. Norms of the reflected signal at
pc and i 2
pc
same step size throughout; now when " is small then the absorption term begins to
dominate in the stability condition and this phenomenon is more pronounced with
higher order layers. Consequently, for given " there is 'optimal' order of the layer
which should be used.
Let us show some pictures of the reflected signal. We just give the representative
cases of the phenomena. First in figure 5.2 we have the reflections using i 3
s with
different value of ". We observe that the reflected signal becomes more and more
irregular when " tends to zero. Then in figure 5.3 we keep " constant, but use
different i. This time the reflected signal becomes more and more irregular when the
order grows, and funnily the number of the extremum points of the reflected signal
'tends' to be the order of the layer plus one.
Then turning to the case - 6= 0, we observe an interesting phenomenon: the
absorption does not improve with increasing order! At first, looking at the numerical
data only we thought that there is no reasonable effect at all. However, looking
at the figure of the reflected signal we noticed that there is indeed an effect but
unfortunately it does not improve the accuracy. Consider figure 5.4. In the upper
left corner there is the case and then the cases
(other -'s gave similar results). We see that - makes the first 'hill' disappear, so
it really does something it should do, but the other 'hills' do not diminish, so that
measuring the reflexion in some norm shows no improvement.
All in all, in the case absorption properties as predicted
by formal theory. However, the orders of convergence measured with the norms we
have used and the formal orders seem to have no simple relation. This suggests that
the asymptotic series we have used are divergent or that the convergence takes place
in such an exotic space that it is irrelevant in practice.
5.2. Practical tests. In the tests above we have kept the number of points in the
layer independent of ". Of course in practice one proceeds differently and one chooses
h to be sufficiently small so that the physical phenomenon is well represented. Then
one expects that by adding few points (say 1 % of total number of points) in the layer
using the same step size as in the initial domain one would get a good absorption.
Here the results are rather disappointing: if we take only a reasonable number of
points in the layer the absorption is not very good.
Consider table 4. On the left we have taken the number of points per pulsewidth
the number of points in the layer. In this case then the number
of points in the initial domain is On the right the corresponding values are
In both cases the simplest layer is the best, the advantage of
higher order layers appearing only with high values of N - . In each case the reflections
are quite big, and they would be unacceptable in practice.
Finally let us remark that because the number of zeros of the layer increases with
the order, they become more and more oscillatory and hence more difficult to dis-
cretize. This explains at least partly why one needs so many points in the layer.



--R

Radiation boundary conditions for wave-like equations

Absorbing boundary conditions for the numerical simulation of waves
Radiation boundary conditions for acoustic and elastic wave calcu- lations
A generalized image principle for the wave equation with absorbing boundary conditions and applications to fourth order schemes
Nonlinear control systems
Numerical absorbing boundary conditions for the wave equation
Approximations of radiation boundary conditions

A general non-reflecting boundary condition for discrete wave-propagation algorithms
Nonlinear homotopies for solving deficient polynomial systems with param- eters
An analysis of particle methods
Etudes de sch'emas num'eriques pour des mod'eles de propagation d'ondes en milieux h'et'erog'enes
A general theorem on local controllability

Mathematica: a system for doing mathematics by computer
--TR
wave equation;absorbing layers;absorbing boundary conditions