two-dimensional wave equation system
M. Lukacova-Medvid'ova 1 2
, G.Warnecke 1
and Y. Zahaykah 1
Abstract
Thesubjectofthepaperisthestudyofsomenonreectingandreectingboundarycondi-
tionsfortheevolutionGalerkinmethods(EG)whichareappliedforthetwo-dimensional
wave equation system. Dierent known tools are used to achieve this aim. Namely,
the method of characteristics,the method of extrapolation, theLaplace transformation
method, and the perfectly matched layer(PML) method. We show that the absorbing
boundaryconditionswhicharebasedontheuseoftheLaplacetransformationleadtothe
Engquist-Majdarst andsecond order absorbingboundaryconditions, see[3]. Further,
following Berenger[1] we consider the PML method. We discretize the wave equation
system withthe leap-frogschemeinside the PML whilethe evolutionGalerkin schemes
areusedinsidethecomputationaldomain. Numericaltestsdemonstratethatthismethod
producesmuch lessunphysical reectedwavesaswellasthe best resultsin comparison
withothertechniquesstudiedin thepaper.
Keywords: hyperbolicsystems,waveequation,evolutionGalerkinschemes,absorbingbound-
aryconditions,reectingboundaryconditions,perfectlymatchedlayer, Laplace transforma-
tion.
1 Introduction
It is well known that to solve a dierential equation numerically in an unbounded domain,
it is necessary to restrict the computations to a bounded domain . Therefore we need to
introduceboundaryconditionson articialboundaries. Asaconsequencethequestionarises
\is there a boundary condition such that the solution of the problem in together with
thisboundaryconditioncoincides exactlywiththerestrictionto of thesolutionintheun-
boundeddomain". Suchboundaryconditionsarecallednonreecting boundaryconditions
orabsorbing boundaryconditions.
On the other hand in the case of a physical boundary it is necessary to reect the solu-
tion(wave) fromtheboundarycompletely. Theseboundaryconditionsare calledreecting
boundaryconditions.
Regarding the above question, many techniques were developed either for the scalar wave
equationorforlinearsystemsofhyperbolicequations,see, e.g.GroteandKeller[5 ,6,7],En-
gquistandMajda[3 ],Higdon[9 ],Thompson[19 ,20 ]andBradly,GreengardandHagstrom[2].
Moreoversuch absorbingboundaryconditionsappearofteninelectromagneticcomputations
where many such conditions were developed, see e.g. the matched layer [10 ] which is based
1
Institut fur Analysis und Numerik, Otto-von-Guericke-Universitat Magdeburg, Univer-
sitatsplatz 2, 39106 Magdeburg, Germany, emails: Gerald.Warnecke@mathematik.uni-magdeburg.de,
Yousef.Zahaykah@mathematik.uni-magdeburg.de
2
TU Hamburg-Harburg, Arbeitsbereich Mathematik, Schwarzenbergstrasse 95, 21073 Hamburg
email:lukacova@tu-harburg.de
impedance matches that of the free space. Others are the Mur conditions [16 ], a perfectly
matchedlayerfortheabsorptionof electromagneticwaves [1 ]and manyothers.
In this paper known nonreecting and reecting boundary conditions for the two dimen-
sional wave equation system are considered and applied for the evolution Galekin methods
(EG). The EG methods were studied by Lukacova, Morton and Warnecke in [12 ], [13 ] for
the multidimensional hyperbolic systems. These methods belong to the class of genuinely
multidimensionalschemes. TypicallythedimensionalsplittingFVmethodsapproximateso-
lutiononlyinnormaldirectionstocellinterfaces. OntheotherhandtheEGmethodsusethe
approximateevolutionGalerkinoperators,whicharebasedontheuseofthebicharacteristics
of multi-dimensionalhyperbolicsystems,suchthatalloftheinnitelymanydirectionsofthe
wave propagationare taken into account, see Section 2. Until now, EG methods had been
applied to easy testcases with periodicorother simpleboundarydata. Forrealy interesting
physicalapplicationsitisnecessarytodemonstratethattheEGschemesarecompatiblewith
more sophisticatednumerical boundarytreatments.
The outline of thepaperwillbe asfollows: inthe nextsection wegive a briefdescriptionof
EG methods. For moredetailson these schemesthereaderisreferred to [12 ],[13 ], [14 ],[15 ],
[17 ]. Inthe secondsection we followThompson[19 ], [20]and usethemethod of characteris-
ticsto derivenonreectingboundaryconditionsaswellasreectingboundaryconditionsfor
thewave equationsystem. We alsoapplydirectlyextrapolationtothewave equationsystem
to producean absorbingboundarycondition. Thisdoesnotalways producesatisfactory nu-
merical results, seee.g. [14]. Thereforewe considerother numerical techniques. Moreoverin
Section 3we linktheLaplacetransformation andFourierseriesto thetwo-dimensionalwave
equation system andderive Engquist-Majda,[3 ], absorbingboundaryconditions. Inthelast
part ofSection 3webriey considerthe absorbingboundaryconditionsbased ona perfectly
matchedlayer(PML)whichdevelpedbyBerenger [1 ]. Inthethirdsectionwepresentresults
of numerical experiments. The above boundary conditions will be combined with the EG
methods inside the computational domain. We use particularly the EG3, EG4 rst order
schemesand theFVEG3second order scheme,see Section2. Finally,we emphasizethatthe
best resultswereobtained whenthe PMLis used.
2 Evolution Galerkin methods
EvolutionGalerkinmethods,EG-methods, wereproposedto approximatehyperbolicconser-
vation laws. An important featureof such methodsisthat they take into account better the
innitelymanydirectionsofpropagationofwavesinmultidimensionalcases. Itiswell-known,
see [12 ], [13 ], [14 ], [15], [17 ], that a basic tool to derive these schemes is the general theory
of bicharacteristics of linear hyperbolic systems. This theory is used to derive the system
of integral equations which is equivalent to the concerned rst order system, e.g. the wave
equation system. Using certain types of quadratures, these integral equations lead to the
approximateevolution operatorthatbuildup theevolution Galerkinscheme.
Let thegeneralform of linearhyperbolicsystembe givenas
U
t +
d
X
j=1 A
j U
x
j
=0; x=(x
1
;:::;x
d )
T
2R d
(2.1)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . .
. . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . .
. . . . . .
. . . . .
. . . .
. . . .
. . . .
. . .
. . .
. . .
. .
. .
. .
. .
. .
. .
. .
.
. .
. .
.
. .
..
. .
. ..
. .
. ..
...
. ..
....
....
. ...
. ....
. ...
...
...
. ...
. ...
. ...
... . ...
... ...
... ...
... . ... ... . .. . ... . .. ... . . . ... .. .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. ...
. . . . . . . . . .
. . . . ... ... . .. .... ... . ... ... . .. ... . ...
. . . .. ... . . . .
. . . . . . . . . . . .
.
. .
. .
.
.
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. .
.
. .
..
. .
..
. .
..
.
. .
. .
..
. .
..
.
. .
..
. .
. .
.
..
. .
..
. .
..
. .
.
. .
. .
. .
..
. .
.
..
. .
..
. .
. .
. .
.
. .
..
. .
..
. .
..
.
. .
. .
. .
.
. .
..
. .
..
. .
..
.
. .
. .
..
. .
. .
. .
.
. .
. .
. .
. .
. .
. .
.
. .
..
. .
..
. .
.
. .
. .
. .
. .
.
. .
. .
. .
..
. .
. .
.
..
. .
. ... ... ...
. .... . .... .... . .... ...
... .... ...
... . .... . .... .... . .... ..
x y t
P 0
Q
i (n)
Figure 1: Bicharacteristicsalong theMach conethroughP and Q
i (n)
wherethecoeÆcientmatricesA
j
;j=1;:::;dareelementsofR pp
andthedependentvariables
are U = (u
1
;:::;u
p )
T
2 R p
. Let A(n) = P
d
j=1 n
j A
j
be the pencil matrix with n =
(n
1
;:::;n
d )
T
being a directional vector in R d
. Then using the eigenvectors of A(n) system
(2.1) can be rewritten in a characteristic form via the substitution W = R 1
U, where the
columns of the matrix R are the linearly independent right eigenvectors of A(n). Since
the coeÆcients of the original system are constants, the bicharacteristics of the resulting
characteristicsystemarestraightlinesPQ
i
andPP 0
generatingthemantleoftheMachcone,
seeFigure1. Diagonalizingthissystemandintegratingalongthebicharacteristicsleadtothe
followingsystemof integralequations
U(P)= 1
jOj Z
O R(n)
2
6
4 W
1 (Q
1 (n);n)
.
.
.
W
p (Q
p (n);n)
3
7
5 dO+
1
jOj Z
O Z
t
0
R(n)S(t+;n)ddO:
Where O is the unit sphere inR d
, jOj its surface measure and S is a nontrivialterm which
we callthesourceterm. Formore detailssee[13].
Nextwegivetheapproximateevolutionoperatorsforthetwo-dimensionalwaveequation
system (3.1) oftwoEG schemes, namely theEG3and theEG4 schemes.
Approximate evolution operators:
EG3 scheme:
P
= 1
2 Z
2
0 (
Q 2u
Q
cos 2v
Q
sin)d+O(t 2
);
u
P
= 1
2 u
P 0+
1
2 Z
2
0
( 2
Q
cos+u
Q (3cos
2
1)+3v
Q
sincos)d
+O(t 2
);
v
P
= 1
2 v
P 0
+ 1
2 Z
2
0
( 2
Q
sin+3u
Q
sincos+v
Q (3sin
2
1))d
+O(t 2
): (2.2)
P
= 1
2 Z
2
0 (
Q 2u
Q
cos 2v
Q
sin)d+O(t 2
);
u
P
= 1
2 Z
2
0
( 2
Q
cos+2u
Q cos
2
+2v
Q
sincos)d+O(t 2
);
v
P
= 1
2 Z
2
0
( 2
Q
sin+2u
Q
sincos+2v
Q sin
2
)d+O(t 2
):
(2.3)
Evoution Galerkin schemes:
Forsimplicityletustaked=2. Considerh>0to bethemeshsizeparameter. We construct
a meshforR 2
,which consistsof thesquaremeshcells
kl
=
(k 1
2
)h;(k+ 1
2 )h
(l 1
2
)h;(l+ 1
2 )h
=
x
k h
2
;x
k +
h
2
y
l h
2
;y
l +
h
2
;
wherek;l2Z. Letusdenote byH
(R 2
) the Sobolevspace of distributionswith derivatives
up to order in L 2
space, where 2 N. Consider the general hyperbolic system given by
the equation (2.1). Let us denote by E(s) : (H
(R 2
)) p
! (H
(R 2
)) p
the exact evolution
operator forthesystem(2.1), i.e.
U(:;t+s)=E(s)U(:;t): (2.4)
We suppose that S m
h
is a nite element space consisting of piecewise polynomials of order
m 0 with respect to the square mesh. Assume constant time step, i.e. t
n
= nt. Let
U n
be an approximation in thespace S m
h
to the exact solutionU(:;t
n
) at time t
n
0. We
considerE
:L
1
l oc (R
2
)!(H
(R 2
)) p
tobeasuitableapproximateevolutionoperatorforE().
In practice we willuse restrictionsofE
to thesubspace S m
h
form0. Then we can dene
the generalclassof evolution Galerkinmethods.
Denition 2.5 StartingfromsomeinitialdataU 0
2S m
h
attimet=0,anevolutionGalerkin
method (EG-method) is recursively dened by means of
U n+1
=P
h E
U
n
; (2.6)
where P
h
is the L 2
projection givenby the integral averages in the following way
P
h U
n
j
k l
= 1
j
kl j
Z
k l
U(x;y;t
n )d xd y:
We denote byR
h :S
m
h
!S r
h
a recovery operator, r m 0 and considerour approximate
evolutionoperatorE
onS
r
h
. Wewilllimitourfurtherconsiderationstothecasewherem=0
and r=1. Takingpiecewiseconstants theresulting schemeswillonlybe of rst order,even
when E
isapproximatedtoa higherorder. Higherorderaccuracy canbeobtainedeitherby
taking m>0,orbyinsertingarecovery stage R
h
beforethe evolutionstep inequation(2.6)
to give
U n+1
=P
h E
R
h U
n
: (2.7)
This approach involves the computation of multiple integrals and becomes quite complex
for higher order recoveries. To avoid this we will consider higher order evolution Galerkin
schemesbasedon thenitevolume formulation instead.
Denition2.8 StartingfromsomeinitialdataU 2S
h
,thenitevolumeevolutionGalerkin
method (FVEG) isrecursively dened by means of
U n+1
=U n
1
h Z
t
0 2
X
j=1 Æ
xj f
j (
~
U n+
t
)d; (2.9)
where Æ
x
j f
j (
~
U n+
t
) represents an approximation to the edge uxdierenceand Æ
x
is dened
by Æ
x
=v(x+ h
2
) v(x h
2
). Thecellboundaryvalue
~
U n+
t
isevolved using theapproximate
evolution operator E
tot
n
+ and averaged a long the cell boundary, i.e.
~
U n+
t
= X
k;l 2Z
1
j@
kl j
Z
@
k l E
R
h U
n
dS
kl
; (2.10)
where
kl
is the characteristic functionof @
kl .
InthisformulationarstorderapproximationE
totheexactoperatorE()yieldsanoverall
higherorderupdatefromU n
toU n+1
. Toobtainthisapproximationinthediscreteschemeit
isonlynecessarytocarryoutarecoverystageateachleveltogenerateapiecewisepolynomial
approximation
~
U n
= R
h U
n
2 S r
h
from the piecewise constant U n
2 S 0
h
, to feed into the
calculation oftheuxes. ToconstructthesecondorderFVEGschemes, forexample, wetake
the rst order accurate approximate evolution operator and dene a bilinear reconstruction
R
h
. Amongmanypossiblerecoveryschemes,whichcanbeused,wewillchooseadiscontinuous
bilinear recovery usingfourpoint averagesat each vertex. It isgiven as
R
h Uj
k l
= U
kl +
(x x
k )
4h (
0x U
kl +1 +2
0x U
kl +
0x U
kl 1 )
+
(y y
l )
4h (
0y U
k+1l
+2
0y U
kl +
0y U
k 1l )
+
(x x
k )(y y
l )
h 2
0y
0x U
kl
;
where
0z
v(z) = 1
2
(v( z+h) v( z h) ): Note that in the updating step (2.9) some nu-
merical quadratures are used instead of the exact time integration. Similarly, to evaluate
the intermediate value
~
U n+
t
in (2.10) either the two dimensional integrals along the cell-
interfaceand aroundtheMach coneareevaluatedexactlyorbymeansof suitablenumerical
quadratures. Since theapproximateevolution operator E
contains integral a longthe base
cone(circlefrom 0to 2),seeFigure1 andequations(2.2) and (2.3),itisdiÆculttoinvolve
boundary conditions in Denitions 2.5 and 2.8. In the next section we will use well known
techniquesto derive suitableboundaryconditions to be usedwith EG-methods.
3 Boundary conditions
Consider theCauchy problemforc>0 and (x;y)2R 2
@
@t
cr
u
v
= 0;
@
@t
u
v
cr = 0;
(3.1)
0
u(x;y;0) = u
0 (x;y);
v(x;y;0) = v
0 (x;y):
(3.2)
Letustake thecomputational domainto be =[0;][0;]R 2
. Moreover, assume that
thereis no physicalboundary,i.e. isobtained onlybytruncationof thewholedomain R 2
.
Aswementioned,to derive appropriateboundaryconditionsforthissystemwewillconsider
knowntechniquesasfollows:
3.1 Boundary conditions based on characteristic methods
In this subsection we use the method of characteristics, Thompson [19 ], [20 ], and derive
articial boundary conditions at the articial interfaces x =0; x =; y =0; and y = of
thecomputationaldomain. Now system(3.1) can berewritteninthe matrixform
U
t +A
1 U
x +A
2 U
y
=0; (3.3)
where
U:=
0
@
u
v 1
A
; A
1 :=
0
@
0 c 0
c 0 0
0 0 0
1
A
; A
2 :=
0
@
0 0 c
0 0 0
c 0 0
1
A
:
PutC=A
2 U
y
and denote thematrix ofthe right eigenvectorsof the coeÆcient matrix A
1
byS=[r
1 jr
2 jr
3
]. It holds
S 1
A
1 S=
1
; where
1ij
=
0 if i6=j;
i
if i=j; i;j=1;2;3:
Withthese notations system(3.3) can bewritten as
S 1
@U
@t +
1 S
1
@U
@x +S
1
C=0
orwe can writeit intheform
l T
i
@U
@t +
i l
T
i
@U
@x +l
T
i
C=0; (3.4)
where l T
i
is aleft eigenvector of A
1 and
i
isthecorrespondingeigenvalue. DeneL
i to be
L
i :=
i l
T
i
@U
@x
; (3.5)
then equation (3.4) hastheform
l T
i
@U
@t +L
i +l
T
i
C=0;
orequivalently
@U
@t
+S L+C=0; (3.6)
where L := (L
1
;L
2
;L
3 )
T
. An easy calculation shows that the eigenvalues of the matrix A
1
are c;0; c,andthe corresponding left andright eigenvectorsare
l T
1
=
1
p
2
; 1
p
2
;0
; l T
2
=(0;0;1); l T
3
=
1
p
2
; 1
p
2
;0
;
r
1
= B
@ 1
p
2
1
p
2
0 C
A
; r
2
= 0
@ 0
0
1 1
A
; r
3
= B
@ 1
p
2
1
p
2
0 C
A
;
respectively. UsingthedenitionofL
i
,equation(3.5),andtheabove lefteigenvectorsweget
L
1
:=
1 l
T
1
@U
@x
= c
p
2
@
@x +
@u
@x
;
L
2
:=
2 l
T
2
@U
@x
=0;
L
3
:=
3 l
T
3
@U
@x
= c
p
2
@
@x +
@u
@x
: (3.7)
Moreoverwehave
S L= 0
B
@ 1
p
2 0
1
p
2
1
p
2 0
1
p
2
0 1 0
1
C
A 0
@ L
1
L
2
L
3 1
A
= 0
B
@ 1
p
2 (L
1 L
3 )
1
p
2 (L
1 +L
3 )
L
2
1
C
A
: (3.8)
Now substitutingfrom (3.8) into (3.6) leadsto thesystem
@
@t +
1
p
2 (L
1 L
3 ) c
@v
@y
= 0; (3.9)
@u
@t +
1
p
2 (L
1 +L
3
) = 0; (3.10)
@v
@t +L
2 c
@
@y
= 0: (3.11)
This system is used forthe updateof the solutionon the interfaces x =0 and x =, while
in the interior domain f(x;y) : (x;y) 2]0;[]0;[g we solve system (3.1) together with
the initial data (3.2). Now the transverse derivatives , i.e. the y derivatives, in (3.9) and
(3.11) are evaluated as usual. The values of the L
i
are determined depending on the sign
of thecorresponding eigenvalues
i
at thecorresponding interface. Namely, at the interface
x=0theeigenvaluesare c; 0; c,whichfrom(3.7) immediatelygivesthatL
2
mustbezero.
The wave corresponding to the eigenvalue
1
= c < 0 is an outgoing one that leaves the
computational domain. Thus L
1
is evaluated directly from its denition (3.7). Since l
3 is
matched with thepositive eigenvalue
3
=c the corresponding wave is an incoming one. It
is this lack of informationthat prevents us from usingthe denitionto evaluate L
3
. Let us
dropthetransversederivativefromequation (3.9),thensubtracting(3.9) from(3.10) we end
withtherelation
@( +u)
@t +
2
p
2 L
3
=0:
Now +uistheone-dimensionalcharacteristicwavethatentersthecomputationaldomain.
Preventingthiswavefromenteringthecomputationaldomainisequivalenttotakingitstime
derivative equal to zero, see [8 ]. This forces L
3
to be equal to zero. A similar argument
forthe interfacex = impliesthat L
1
and L
2
must be zero while L
3
is evaluated from its
denition(3.7).
Forthecasey=0andy=weuseasimilarargumentasabove. Atcornerswecombineboth
onlyabsorbing boundaryconditions butalso reecting ones. To showthis, assume that the
interfacex=0is a reector. Todene areected boundaryconditionat x=0 thefunction
mustvanishforalltimet>0atx=0. Thenequation(3.9)impliesthatL
3
=L
1 p
2c
@v
@y .
Remark 3.12 One of the simplest absorbing boundary conditions is to extrapolate the data
from the interior domain to the cells that lie on the boundary. Mathematically, for the wave
equation system (3.1), this can be done by taking
@'
@s
= 0, where ' 2 f;u;vg and the
derivativewithrespecttosstands forthedirectionalderivativeinsomedirection parametrized
by the parameter s. We will show however in our numerical experiments that this approach
can lead to considerable unphysical reections in the solution.
3.2 Absorbing boundary conditions based on the Laplace transformation
In thissubsection we recall the use of Fourier seriesand theLaplace transformation to con-
structanabsorbingboundaryconditiontobeused,say,attheboundaryf(;y) :0yg.
We will expand the solution of the wave equation system in a complex Fourier series. Sub-
stituting into the wave equation system and applying the Laplace transformation allows a
derivationof Engquist-Majda [3] absorbingboundaryconditions.
Considerthewave equationsystem(3.1)together withtheinitialcondition(3.2) inR 2
. Sup-
posethattheinitialfunctions
0
( x;y),u
0
( x;y)andv
0
(x;y)areallsmoothandvanishoutside
thedomain=[0;][0;] . Welookforaboundaryconditionthatholdsattheboundaryof
thedomain@suchthatthesolutionofsystem(3.1) accompaniedwiththeinitialdata(3.2)
on the domain together withthisconditioncoincides withthesolutionon theunbounded
domain. The Laplace transformation of a function f(t) will be denoted by
^
f(s) = L(f(t)),
and its inversebyf(t)=L 1
(
^
f(s)). Theyare denedasfollows
^
f(s):=L(f(t))= Z
1
0
f(t)exp( st) dt;
f(t):=L 1
(
^
f(s))= 1
2i Z
a+i1
a i1
^
f(s)exp(st) dt:
To produceanexact absorbingboundaryconditionusingtheLaplace transformationweuse
the Fourier series and expand the solution U(x;y;t) at the points to the right of the line
x= as
U(x;y;t)= 1
X
k= 1 U
k
(x;t)expi(k
y) (3.12)
Substituting(3.12) into thewave equation system(3.3) yields
U k
t +A
1 U
k
x +ik
A
2 U
k
=0: (3.13)
ApplyingtheLaplacetransformationto(3.13)andusingthepropertiesL(f 0
(t))=sL(f(t))
f 0
(0) and L
@f
@x
=
@
@x
( L(f)) as well as the assumption that the initial data vanishin the
exteriorof ,we get
s
^
U k
( x;s)+A
1
^
U k
x
(x;s)+ik
A
2
^
U k
(x;s)=0: (3.14)
Expressing u^ and ^v in terms of
^
we end with the following second order dierential
equation
^
k
xx
"
s 2
c 2
ik
2
#
^
k
=0: (3.15)
The boundednessof thesolutionforx!1 impliesthat thesolutionof(3.15) is
^
k
(x;s)=A( s;y)exp 0
@ s
c s
1
ick
s
2
x 1
A
: (3.16)
The valuesofu^ k
and ^v k
read
^ u k
(x;s) = s
1
ick
s
2
A(s;y)exp 0
@ s
c s
1
ick
s
2
x 1
A
; (3.17)
^ v k
(x;s) = ick
s
A(s;y)exp 0
@ s
c s
1
ick
s
2
x 1
A
: (3.18)
Ifweapplythe operator
L= 0
@
@
@x +
s
c s
1
ick
s
2 1
A
(3.19)
to thefunctions
^
k
; u^ k
and ^v k
we getL
^
k
=0,L^u k
=0 and L^v k
=0. Hence,at x=,we
have thefollowing boundaryconditionforeach component u k
(x;t)
L 1
n
L
^
k
o
=0; L 1
n
L^u k
o
=0; L 1
n
L^v k
o
=0: (3.20)
Let g(t) = L 1
s q
1 a
2
s 2
, where a = ick
, then at x = condition (3.20) can be
rewrittenas
k
x
(x;t)+ 1
c Z
t
0
g(t ) k
(x;)d =0;
u k
x
(x;t)+ 1
c Z
t
0
g(t )u k
(x;)d =0;
v k
x
(x;t)+ 1
c Z
t
0
g(t )v k
(x;)d =0: (3.21)
It is clear from equations (3.21) that this boundary condition is local in position but un-
fortunately not local in time. Using q
1 a
2
x 2
= 1+O 1
x 2
we approximate the operator
L as L =
@
@x +
s
c
: Hence (3.20) gives the following local in timeEngquist-Majda rst order
absorbingboundaryconditionforthefunctionat theinterfacex=.
@
@x +
1
c
@
@t
=0: (3.22)
Using 1 a
2
x 2
=1 a
2
2x 2
+O 1
x 4
we approximatethe operatorL asL=
@
@x +
s
c
1 a
2
2s 2
:
If we substitute sL into (3.20) then we obtain the following Engquist-Majda second order
absorbingboundaryconditionforat x=.
@ 2
@x@t +
1
c
@ 2
@t 2
c
2
@ 2
@y 2
=0: (3.23)
Note thattheabsorbingboundaryconditionsforother interfacescan bederivedinananalo-
gous way. Forthe descritizationof(3.22) and(3.23) we refereto [3].
3.3 Absorbing boundary condition based on a perfectly matched layer
In the previous two subsections absorbing boundary conditions were dened using certain
types of operators. Although such boundary conditions absorb the incoming waves, the
absorption is not complete for some cases, as we will see in the numerical examples in the
next section. By construction, such boundaryconditionsare derived to absorbthe reected
wavesprovidedthattheincidentplanewavesarenormaltotheinterfaces,formoredetailssee
Grote [4 ]. Inthissubsectionwe willapplyboundaryconditionsthat arevalidindependently
of the directionof the incidentwave orits frequency. Thistype ofboundaryconditions was
introduced by Berenger [1 ] for electromagnetic problems. Afterwards many authors have
beenusingsuchtechnique,forexample, Fang [11]and Quanand Geers[18]. Themainpoint
is to surround the computational domain by a layer that has the same impedance, roughly
speakingresistanceto thewave propagation,asthefreedomain, thedomainwhere thewave
equation system(3.1) valid. Thusforplanewaves, atleasttheoretically,Berengerhasshown,
fortheelectromagneticcase, thatthere isno reection. Formore detailsonPMLmethodsee
[1 ], [11 ], [18 ]. To denePML layerswe splitthe functioninto two components x
and y
,
i.e. we write= x
+ y
. Then we dene
Denition 3.24 A layer(R 2
) in which thesystem,damped waveequation system,
@
@t
x
y
+
?
x 0
0
?
y
x
y
c
@u
@x
@v
@y
=0;
@
@t
u
v
+
x 0
0
y
u
v
cr(
x
+ y
)=0: (3.25)
issatisedwhere
x
;
?
x
;
y
;and
?
y
arefunctionsof xandy,dampingparameters,iscalled
a perfectly matchedlayer.
LetPML1,characterisedby
x
1
;
?
x
1
;
y
1
;
?
y
1
andPML2,characterisedby
x
2
;
?
x
2
;
y
2
;
?
y
2
,
be two perfectlymatchedlayers. Supposethatx=0 istheinterfacebetweenthetwo layers.
If =
?
and if
y1
=
y2
then Berenger has shown that the two layers PML1 and PML2
producenoreectionfromtheoutgoingwavesattheinterfacex=0. Thedampingparameter
isdenedas()=
m
Æ
2
where
m
isdeterminedbythetheoretical reectionatnormal
incidence, isthedistancefromtheinterfaceandÆ isthethicknessofthelayer, seeBerenger
[1 ]. We implementthe perfectly matchedlayermethod asfollows: intheinteriordomainwe
usedthe evolution Galerkin scheme, i.e. by means of (2.6) in thecase of EG-schemes or by
(2.9), (2.10) in the case of nite volume EG-schemes. In the PML layer we used leap-frog
scheme,seeBerenger[1 ] orQuan and Geers[18 ] forthe discretisationofthe PMLlayer.
Inthe followingnumerical experimentswe test and compare theabsorbing boundarycondi-
tionsusingthemethodsofcharacteristics, extrapolation,theLaplacetransformationandthe
perfectlymatchedlayer(PML). The initialdata aregiven eitherasplane waves oraspulses
centered at certain pointsinthe calculation domain. Thesenumericaltests indicatethat all
thedevelopedabsorbingboundaryconditionsworkedcorrectlyinthecaseofonedimensional
problems. However it is well known that there are big problems with waves having an inci-
dent angle notnormal to the boundary. Therefore, for multidimensionalproblemsthe PML
methodgivesthebestresults.
Example4.1
Totestthereectingandtheabsorbingboundaryconditionsthatarebasedonthemethodof
characteristicsweconsiderthewaveequation system(3.1) andtakex=y=1,t=0:5,
c=1 and=[0;100][0;100]. InFigure2theinitialincidenthalfsinewavewhichistaken
from the wave
0
=sin(
x
24 ); u
0
=sin(
x
24
) and v
0
=0 is moved to the left and is reected
at the boundaryx =0. In Figure 3 theGaussian wave
0
=0:5exp ( ln(2)(
x 10
3 )
2
); u
0
=
0:5exp ( ln(2)(
x 10
3 )
2
) and v
0
=0 is moved to theright and isabsorbed at the boundary
x = 100. We use the rst order EG4 and the second order FVEG3 schemes, respectively.
In Figure 2 we see that the one-dimensional half sine wave propagates to the left. Upon
arrival at the left wall boundary it is reected back. This happensafter T=75:0 (150 time
steps). The phase of the reected wave is oppositeto that of theincident wave, asrequired
bythereectedboundarycondition. In Figure3theinitialGaussianwave propagates tothe
right. After the arrivalat the right articial boundaryit leaves the domain as we see after
T=62:5 (125timesteps). We have madeseveral other 1Dexamples fordierentinitial data.
All absorbing boundary conditions described above yield analogous results. These results
indicatethat such absorbing boundary conditions are doing wellif the initial data are one-
dimensionaldata. Unfortunately,thisisnotthecase iftwo-dimensionalinitialdataareused
aswe willseeintheExample4.3.
Example4.2
In this test we apply the PML method. We consider the two dimensional wave equation
systeminsidethedomain =[ 1;1][ 1;1]. The initialdataaretaken to be
0
(x;y)= exp( 30((x+0:85) 2
+y 2
)); u
0
(x;y)=0=v
0 (x;y):
Inthesimulationweuse a100100 mesh,the thickness ofthematchedlayeris taken to be
10cellsand
m
is chosen suchthatthetheoretical reectionat thenormalincidenceis10 5
.
Inside the computational domain we update using the rst order EG3 scheme. Inside the
PMLlayerweuseleapfrogscheme. We taketheCourant-Friedrichs-LewynumberCFL=0.4,
we deneit asCFL= ct
h
; whereh=x=y. We considera timesequence, 0.2,0.4, 0.6,
1.0, of absolute times (T). The isolines of the solution are shown in Figure 4. Moreover,
Figure5 shows theresult of thesame problemusingextrapolation forabsolute timesT=0.2
and T=1.0. Figure 5 clearly shows that, in general, using extrapolation as an absorbing
boundaryconditionmayproduced unphysicalreectedwaves. It isclear fromFigure 4 that
thePMLmethodproducesonlyamarginalamountofreectionfromthearticialboundaries.
Thus inthe case of the wave equation system this method is the preferableone in order to
implement absorbing boundaryconditions.
0 10 20 30 40 50 60 70 80 90 100
−1 0 1
Initial Data (T=0.0)
0 10 20 30 40 50 60 70 80 90 100
−1 0 1
T=12.5
0 10 20 30 40 50 60 70 80 90 100
−1 0 1
T=25.0
0 10 20 30 40 50 60 70 80 90 100
−1 0
1 T=50.0
0 10 20 30 40 50 60 70 80 90 100
−1 0
1 T=75.0
0 10 20 30 40 50 60 70 80 90 100
−1 0
1 T=87.5
Figure2: The numerical solutionof thecomponent usingtheEG4scheme together witha
reectorat x=0. Halfsinewave initialdata.
0 10 20 30 40 50 60 70 80 90 100
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 0.2 0.3 0.4
0.5 T=2.5
0 10 20 30 40 50 60 70 80 90 100
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 0.2 0.3 0.4 0.5
T=12.5
0 10 20 30 40 50 60 70 80 90 100
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 0.2 0.3 0.4
0.5 T=25.0
0 10 20 30 40 50 60 70 80 90 100
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 0.2 0.3 0.4
0.5 T=37.5
0 10 20 30 40 50 60 70 80 90 100
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 0.2 0.3 0.4
0.5 T=50.0
0 10 20 30 40 50 60 70 80 90 100
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 0.2 0.3 0.4
0.5 T=62.5
Figure 3: The initial data(- -), theexact solution(-) and the approximate solution(-.-.) of
thecomponent usingthesecond orderFVEG3scheme. Absorberat x=100. Gaussianpuls
initial data.
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
PML Method; T=0.2; first order EG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
PML Method; T=0.4; first order EG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
PML Method; T=0.6; first order EG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
PML Method; T=1.0; first order EG3 scheme
Figure4: Isolinesof thecomponent of theapproximatesolutionusingthePML absorbing
boundarycondition.
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Extrapolation; T=0.2; first order EG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Extrapolation; T=1.0; first order EG3 scheme
Figure5: Isolinesof the component of theapproximatesolution usingextrapolation asan
absorbingboundarycondition.
In Example 4.1 we used absorbing boundary conditions that based on the methods of ex-
trapolation and characteristics. The initial data were one-dimensional and such boundary
conditions give good results with respect to preventing the outgoing wave from reentering
the computational domain. However if the initial data are not one-dimensional then these
boundaryconditionsproducereectedwavesaswehave alreadyseeninExample4.2. Inthis
examplewe showthatabsorbingboundaryconditionsbasedonthemethodofcharacteristics
and the Laplace transformation can also produce reected waves ingeneral. Again we con-
siderthetwo-dimensionalwaveequation systeminsidethedomain=[ 1;1][ 1;1]. The
initialdataaretaken to be
0
(x;y)= exp( 15((x+0:5) 2
+y 2
)); u
0
(x;y)=0=v
0 (x;y):
WeuserstandsecondorderEG3andFVEG3methodsandcombineitwithdierentbound-
ary conditions. We take the CFL=0.55, the absolute time T=0.825 and consider100100
mesh. Figures 6 top-left and top-right show the isolines of the solution for the absorbing
boundarycondition usingPML and extrapolation methods,respectively. Figures6 bottom-
left and bottom-right show theisolines usingcharacteristics and Engquist-Majda rst order
condition methods, respectively. In Figure 7 we show the isolines of the solution for the
reectingboundaryconditionbasedonthemethodofcharacteristicsatx= 1. Theseresults
demonstrate that theabsorbing boundaryconditionsbased on themethod ofcharacteristics
aswellastheonebasedontheLaplacetransformationcanproducemuchstrongerunphysical
oscillationsandbackwardreections. PicturesatthebottomofFigures6aswellaspictureat
top-right are infactcomparableto thepictureinFigure 7,where theresultof thereecting
boundaryconditionis shown. In Figures8 weusethesecond orderFVEG3 methodtogether
withPML(top-left), extrapolation(top-right)and Engquist-Majdaabsorbingbounderycon-
dition ofsecond order (bottom). The eect ofthe reectionis clearintop-rightaswellasin
bottom pictures. FinallyinFigure9wedrawtheisolinesofforabsolute timesT=0.55 and
T=0.99 when second order FVEG3 method is used together with extrapolation, PML and
Engquist-Majdasecond order conditionmethods, respectively. We see thatless eect of the
boundaryinthecase of PMLmethod.
Conclusions
In this paper we derive some absorbing and reecting boundary conditions for the two-
dimensional wave equation system. These conditions are based on the known methods of
characteristics,extrapolation,Laplacetransformation(Engquist-Majdaconditions)andPML.
We use them with EG-schemes. Inside the computational domain we always used the EG-
schemeswhicharegiven inSection2. The numericalexperimentswhichwe presentedinthis
paper showed clearly thatfor one-dimensionalcases all methodsgive good results. Further,
for two-dimensional problems we noticed that the PML method produced less nonphysical
oscillation and reected waves when used together with EG-methods. In order to compare
computational costs of dierent ABC techniques applying to EG-methods, we have taken
one particular EG-method, e.g. thesecond order FVEG3. We carry out the calculationsin
Example 8 with meshes 100, 200 and 400, respectively. The results are tabulated below in
Table1. Weseethatthecheapestmethodisthemethodofextrapolationwhiletheexpensive
one isclearlythePML method.
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
PML; ABC; T=0.825; first order EG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Extrapolation; ABC; T=0.825; first order EG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Engquist−Majda first order condition; ABC; T=0.825; first order EG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Method of Characteristics; ABC; T=0.825; first order EG3 scheme
Figure6: Isolinesofthecomponentoftheapproximatesolution,absorberbasedatx= 1:
top-left: based on the PML, top-right: based on extrapolation, bottom-left: based on the
methodof characteristics and bottom-right: based onLaplace transform.
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Method of Characteristics; Reflecting BC; T=0.825; first order EG3 scheme
Figure 7: Isolines of the component of the approximate solution, reector based on the
methodof characteristics at x= 1.
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
PML; ABC; T=0.825; second order FVEG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Extrapolation; ABC; T=0.825; second order FVEG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Engquist−Majda second order condition; ABC; T=0.825; second order FVEG3 scheme
Figure 8: Isolines of the component of the approximate solution, absorber based on the
PML (top-left), on extrapolation (top-right) and on Engquist-Majda second order condition
(bottom).
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Extrapolation; ABC; T=0.55; second order FVEG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Extrapolation; ABC; T=0.99; second order FVEG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
PML; ABC; T=0.55; second order FVEG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
PML; ABC; T=0.99; second order FVEG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Engquist−Majda second order codition; ABC; T=0.55; second order FVEG3 scheme
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8
Engquist−Majda second order condition; ABC; T=0.99; second order FVEG3 scheme
Figure9: Isolinesofthecomponent oftheapproximatesolutionusingextrapolation,PML
andEngquist-Majda second order conditionmethods.
100100 cost 0.92s 0.01s 0.43s
200200 cost 6.00s 0.06s 3.26s
400400 cost 44.84s 0.31s 25.82s
Table1: Computational timeresults.
With thispaperwe havedemonstrated that itis feasableto implementboundaryconditions
withinthe framework of evolution Galerkinschemesbased onbicharacteristics.
Acknowledgements
This research wassupportedunder theDFG grant No. Wa 633/6-2 of Deutsche Forschungs-
gemeinschaft, by the grants GACR 201/00/0557 and CZ 39001/2201 of the University of
Technology Brno aswellas by theVolkswagenStiftungand DAAD Agencies. We would like
to thank the referees for a large number of helpful comments which enabled us to improve
the paper.
References
[1] J.-P. Berenger. A perfectly matched layer for the absorbtionof electromagneticwaves.
J. Comput. Phys.,114:185{200, 1994.
[2] A. Bradly, L.Greengard, and T. Hagstrom. Nonreecting boundary conditions for the
time-dependentwave equation. J.Comput. Phys.,180(1):270{296, 2000.
[3] B.EngquistandA.Majda. Absorbingboundaryconditionsforthenumericalsimulation
of waves. Math.Comput., 31(139):629{651, 1977.
[4] M. J. Grote. Nonreecting boundary conditions for time dependent wave propagation.
SeminarfurAngewandteMathematik(SAM),ResearchreportNo.2000-04,Zurich,2000.
[5] M. J. Grote and J. B. Keller. Exact nonreecting boundary conditions for the time
dependentwave equation. SIAM Journal of Appl. Math., 55(2):280{297, 1995.
[6] M. J. Grote and J. B. Keller. Nonreecting boundary conditions for time-dependent
scattering. J.Comput. Phys., 127:52{65, 1996.
[7] M.J.GroteandJ.B.Keller.NonreectingboundaryconditionsforMaxwell'sequations.
J. Comput. Phys.,139:327{342, 1998.
[8] G.W. Hedstrom. Nonreectingboundaryconditionsfornonlinearhyperbolicsytems. J.
Comput. Phys.,30:222{237, 1979.
[9] R. L. Higdon. Numerical absorbing boundaryconditions for the wave equation. Math.
Comput., 49(179):60{90, 1987.
[10] R.Holland and J.W. Williams. Total-eld versusscattered-eld nite-dierencecodes:
A comparativeassessment. IEEE Transaction on Nuclear Science,NS-30(6):4583{4588 ,
1983.
perfectlymatchedlayer. J.Comput. Phys.,129:201{219, 1996.
[12] M. Lukacova,K.W.Morton,andG. Warnecke. FinitevolumeevolutionGalerkinmeth-
ods for Euler equations of gas dynamics. Int. J. Num. Methods in Fluids, accepted,
2001.
[13] M. Lukacova, K.W.Morton, and G. Warnecke. Evolution Galerkinmethods forhyper-
bolicsystemsintwospace dimensions. Math. Comput., 69(232):1355{13 84 , 2000.
[14] M. Lukacova, J. Saibertova, G. Warnecke, and Y. Zahaykah. On evolution Galerkin
methodsfortheMaxwelland thelinearizedEulerequations. Acceptedto (Appl.Math.),
2003.
[15] M. Lukacova, G. Warnecke, and Y. Zahaykah. Third order nite volume evolution
Galerkin (FVEG) methods for two-dimensional wave equation system. submitted to
East-West J. Numer. Math.,2002.
[16] G. Mur. Absorbingboundaryconditionsfornite-dierenceapproximationof thetime-
domainelectromagneticeldequations. IEEETransactiononElectromagneticCompata-
bility,EMC-23(4):377{382, 1981.
[17] S. Ostkamp. Multidimensional characteristic Galerkinschemes and evolution operators
forhyperbolicsystems. Math. Meth. Appl. Sci., 20:1111{1125, 1997.
[18] Quan Qi and T. L. Geers. Evolution of the perfectlymatched layer for computational
acoustics. J.Comput. Phys.,139:166{183, 1998.
[19] K. W. Thompson. Time dependent boundary conditions for hyperbolic systems. J.
Comput. Phys.,68:1{24, 1987.
[20] K. W. Thompson. Time dependent boundaryconditions for hyperbolic systems, II. J.
Comput. Phys.,89:439{461, 1990.
