two-dimensional wave equation system 1
M. Lukacova-Medvid'ova, 2
G. Warnecke 3
and Y. Zahaykah 3
Abstract
ThesubjectofthepaperistheanalysisofstabilityoftheevolutionGalerkin(EG)methods
forthetwo-dimensionalwaveequationsystem. WeapplyvonNeumannanalysisanduse
theFouriertransformationtoestimatethestabilitylimitsofboththerstandthesecond
order EGmethods.
Keywords: hyperbolicsystems, wave equation,evolutionGalerkinschemes, discreteFourier
transformation,amplicationmatrix,CFL condition.
1 Introduction
EvolutionGalerkin methods(EG-methods) were proposedto approximaterst order hyper-
bolicproblems. These schemes were introduced by Morton, see, e.g., [9 ] forscalar problems
and [10 ] for one-dimensional systems. The rst generalization to two-dimensional systems
was made in [11 ] by Ostkamp for the wave equation system as well for theEuler equations
ofgasdynamics. In [4 ]Lukacova, Mortonand Warneckesystematically studiedapproximate
evolution operators and constructed further EG-schemes with better accuracy and stability
properties. FurtherEGschemesaswellastheapproximateevolutionoperatorofthesolution
for thewave equationsystem inthreespace dimensionswere derived in[13 ]. These methods
or their nitevolume versions were applied to the nonlinearEuler equations, see [1 ], [6 ], as
wellasto thelinearizedEuler equationsand Maxwellequations [7 ]. Higherorder nite vol-
ume EG-methodshave beenintroduced andstudied in[3],[5 ],[6] and [8].
The main objective of this paper is the analysis of the stability of the evolution Galerkin
schemes. We derive a necessary and suÆcient stabilitycondition forthe evolution Galerkin
scheme (EG3 scheme) applied to the wave equation system in two space dimensions. The
discreteFouriertransformationisusedto obtaintheamplicationmatricesoftheseschemes.
Using estimates of thespectral norm we nda suÆcient stability condition. We derive am-
plication matrices for the rst and the second order nitevolume schemes (FVEG) based
on the approximate evolution operators. The spectral radius of the amplication matrices
is estimated experimentally by abuilt-in Matlab procedure. Hence the stabilitylimitof the
schemesisestimated also numerically.
The outline ofthispaperis asfollows: inthenext sectionwe surveythegeneral theorythat
1
ThisresearchwassupportedbytheVolkswagenStiftungAgency,bytheDeutscheForschungsgemeinschaft
GrantNo.Wa633/6-2andpartiallybytheEuropeannetworkHYKE,fundedbytheECascontractHPRN-
CT-2002-00282. Authorsgratefully acknowledgethesesupports.
2
TUHamburg-Harburg,ArbeitsbereichMathematik,Schwarzenbergstrasse95,21073Hamburg,Germany,
email: lukacova@tu-harburg.de
3
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
schemes. The exact integral equations as well as the approximate evolution operators for
thetwo-dimensionalwave equation systemare given inSection4. InSection 5we introduce
the discrete Fourier transformation as well as the spectral norm that serve as tools in our
analysis. In Section 6 we present the derivation of the suÆcient and necessary condition
and compare the theoretical limit, that we obtained by means of the Fourier analysis with
the experimental one. In Section 7 we consider the rst order nitevolume schemes based
on theapproximateevolution operatorE const
. We determinethe amplicationmatrices and
estimate theirstability limits. Finally inSection 8 we determine the amplication matrices
ofthesecondordernitevolumeschemesbasedontheapproximateevolutionoperatorE bil in
and estimate thestabilitylimits.
2 General theory
In this section we recall the exact integral equations for a general linear hyperbolicsystem
usingtheconcept ofbicharacteristics. Consider ageneral formof linearhyperbolicsystem
U
t +
d
X
k=1 A
k U
x
k
=0; x=(x
1
;:::;x
d )
T
2R d
; (2.1)
wherethecoeÆcientmatricesA
k
;k=1;:::;d areelementsofR pp
andthedependentvariables
are U=(u
1
;:::;u
p )
T
=U(x;t)2R p
. Let A(n)= P
d
k=1 n
k A
k
bethe pencil matrix,where
n = (n
1
;:::;n
d )
T
is a unit vector in R d
. Since the system (2.1) is hyperbolic the matrix
A(n) has p real eigenvalues
k
, k = 1;:::;p, and p corresponding linearlyindependent right
eigenvectorsr
k
=r
k
(n);k=1;:::;p. LetR=[r
1 jr
2 j:::jr
p
]bethematrixofrighteigenvectors.
WedenethecharacteristicvariableW =W (n)as@W (n)=R 1
@U. Sincethesystem(2.1)
has constant coeÆcient matrices A
k
we have W=R 1
U orU=RW .
Transformingsystem(2.1) bymultiplyingitwith R 1
from theleftweget
R 1
U
t +
d
X
k=1 R
1
A
k RR
1
U
x
k
=0: (2.2)
Let B
k
=R 1
A
k R =(b
k
ij )
p
i;j=1
,where k =1;2;:::;d, then thesystem (2.2) can berewritten
in thefollowingformusingcharacteristic variables
W
t +
d
X
k=1 B
k W
x
k
=0:
NowwedecomposeB
k
into thediagonalpartD
k
and therestpartB 0
k ,i.e.B
k
=D
k +B
0
k . We
obtain
W
t +
d
X
k=1 D
k W
x
k
= d
X
k=1 B
0
k W
x
k
=:S: (2.3)
The i-th bicharacteristiccorresponding to thei-th equation of(2.3) is denedby
dx
i
d
~
t
=b
ii
(n)=(b 1
ii
;b 2
ii
;:::;b d
ii )
T
;
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . .
. . . . . . . . . . .
. . . . . . .
. . . . . . . .
. . . . .
. . . . . .
. . . . .
. . . .
. . . .
. . . .
. . .
. . .
. . .
. .
. .
. .
. .
. .
. .
. .
.
. .
. .
.
. .
..
. .
. ..
. .
. ..
...
. ..
....
....
. ...
. ....
. ...
...
...
. ...
. ...
. ...
... . ...
... ...
... ...
... . ... ... . .. . ... . .. ... . . . ... .. .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. ...
. . . . . . . . . .
. . . . ... ... . .. .... ... . ... ... . .. ... . ...
. . . .. ... . . . .
. . . . . . . . . . . .
.
. .
. .
.
.
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. . .
.
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
.
. .
.
. .
..
. .
..
. .
..
.
. .
. .
..
. .
..
.
. .
..
. .
. .
.
..
. .
..
. .
..
. .
.
. .
. .
. .
..
. .
.
..
. .
..
. .
. .
. .
.
. .
..
. .
..
. .
..
.
. .
. .
. .
.
. .
..
. .
..
. .
..
.
. .
. .
..
. .
. .
. .
.
. .
. .
. .
. .
. .
. .
.
. .
..
. .
..
. .
.
. .
. .
. .
. .
.
. .
. .
. .
..
. .
. .
.
..
. .
. ... ... ...
. .... . .... .... . .... ...
... .... ...
... . .... . .... .... . .... ..
x y t
P 0
Q
i (n)
Figure 1: Bicharacteristics alongthe Mach conethroughP andQ
i
(n), d=2.
where i = 1;:::;p. The diagonal entries b k
ii
of the matrices B
k
, k = 1;:::;d, i = 1;:::;p,
create the so-called ray velocity vector b
ii
. We considerthe bicharacteristics backwards in
time and set the initial conditions x
i
(t+t;n) = x for all n 2 R d
and i = 1;:::;p, i.e.
x
i (
~
t ;n)=x b
ii
(n)(t+t
~
t).
We will integrate the i-th equation of the system (2.3) from the point P down to the point
Q
i
(n),where thebicharacteristic hitstheplane at timet, see Figure1. Note that bicharac-
teristicsarestraight lines because thesystem is linearhavingconstant coeÆcients. Nowthe
i-th equationreads
@w
i
@t +
d
X
k=1 b
k
ii
@w
i
@x
k
= 0
@ d
X
j=1;i6=j
b 1
ij
@w
j
@x
1 +b
2
ij
@w
j
@x
2
+:::+b d
ij
@w
j
@x
d
1
A
=S
i
; (2.4)
whereP (x;t+t)2 R p
R
+
is taken to be a xed point,while Q
i
(n)= (x
i
(n;t);t) =
(x tb
ii
;t). Takinga vector
i
=(b 1
ii
;b 2
ii
;:::;b d
ii
;1), we candene thedirectionalderivative
dw
i
d
i
=
@w
i
@x
1
;
@w
i
@x
2
;:::;
@w
i
@x
d
;
@w
i
@t
i
=
@w
i
@t +b
1
ii
@w
i
@x
1 +b
2
ii
@w
i
@x
2
+:::+b d
ii
@w
i
@x
d :
Hence thei-th equation (2.4) can berewritten asfollows
dw
i
d
i
=S
i
= d
X
j=1;i6=j
b 1
ij
@w
j
@x
1 +b
2
ij
@w
j
@x
2
+:::+b d
ij
@w
j
@x
d
:
Integration from P to Q
i
(n) gives
w
i
(P) w
i (Q
i
(n)) =S 0
i
; (2.5)
where
S 0
i
= Z
t+t
t S
i (x
i (
~
t ;n);
~
t ;n)d
~
t= Z
t
0 S
i (x
i
(;n);t+t ;n)d:
Reverse transformationof (2.5) into thesystem writteninoriginalphysicalvariablesisdone
by multiplication with R from the left and (d 1) dimensional integration of the variable
n over the unit sphereO in R . This leadsto the integral representation of the solutionat
pointat time t+t
U(P)=U(x;t+t)= 1
jOj Z
O R(n)
0
B
B
B
B
B
@ w
1 (Q
1 (n);n)
w
2 (Q
2 (n);n)
w
3 (Q
3 (n);n)
.
.
.
w
p (Q
p (n);n)
1
C
C
C
C
C
A dO+
~
S; (2.6)
where
~
S =(
~
S
1
;
~
S
2
;:::;
~
S
p )
T
= 1
jOj Z
O
R(n)S 0
dO= 1
jOj Z
O Z
t
0
R(n)S(t+t ;n)ddO
and jOjcorresponds to themeasure ofthedomain ofintegration.
3 Exact integral equations and approximate evolution opera-
tors for the wave equation system
Inthis section we illustratethe application of general theory of the bicharacteristics forthe
two-dimensionalsystemofwave equation. We recalltheexact integralequationsand present
their possible approximation, the so-called EG3 approximate evolution operator. Consider
thetwo dimensionalwave equationsystem
t +c(u
x +v
y
)= 0;
u
t +c
x
= 0;
v
t +c
y
= 0;
(3.1)
where c is a given constant representing the speed of sound. We will recall here the exact
integralequationsderivedin[4 ]. LetP =(x;y;t+t), P 0
=(x;y;t),Q=(x+ctcos;y+
ctsin;t)=(x+ctn();t)and theso-called sourcetermbe given as
S = c
u
x sin
2
(u
y +v
x
)sincos+v
y cos
2
; (3.2)
thenexact integralequations forthewave equation system(3.1) are given as
P
= 1
2 Z
2
0 (
Q u
Q
cos v
Q
sin)d+
~
S
1
; (3.3)
u
P
= 1
2 u
P 0
+ 1
2 Z
2
0 (
Q
cos+u
Q cos
2
+v
Q
sin cos)d+
~
S
2
; (3.4)
v
P
= 1
2 v
P 0
+ 1
2 Z
2
0 (
Q
sin+u
Q
cos sin+v
Q sin
2
d+
~
S
3
; (3.5)
~
S
1
= 1
2 Z
2
0 Z
4t
0
S(x+cn();t+t ;)d d;
~
S
2
= 1
2 Z
2
0 Z
4t
0
cosS(x+cn();t+t ;)d d
1
2 Z
2
0 Z
4t
0
c
x
(x;t+t )sin 2
c
y
(x;t+t )sincos
dd;
~
S
3
= 1
2 Z
2
0 Z
4t
0
sinS(x+cn();t+t ;)d d
1
2 Z
2
0 Z
t
0
c
y
(x;t+t )cos 2
c
x
(x;t+t )sincos
d d:
The above integral equations give us an implicit formulation of the solution at the point
P =(x;y;t n+1
). In order to obtainan explicitnumericalscheme itis necessary to use some
numerical quadratures in order to approximate the time integral from 0 to t. Using the
backwardrectangle ruleleadstoan O(t 2
) approximationof thetimeintegralsappearingin
~
S
1 ,
~
S
2 and
~
S
3
. Furtherwe usethefollowingresult [4 ,Lemma 2.1]
t Z
2
0
S(t;)d = Z
2
0
(ucos+vsin)d;
t Z
2
0
S(t;)cosd = Z
2
0
(ucos2+vsin2)d;
t Z
2
0
S(t;)sind = Z
2
0
(usin2+vcos2)d: (3.6)
Note that these formulaeallow to replace thederivatives ofour dependent variablesinS by
the variablesthemselves. Rectangle rule approximation forthe timeintegraland (3.6) yield
the so-called EG3 approximate evolution operator. We refer a reader to [4 ], [13] for other
approximateevolution operatorsEG1-EG4.
Approximate evolution operator for EG3
P
= 1
2 Z
2
0 (
Q 2u
Q
cos 2v
Q
sin)d+O(t 2
) (3.7)
u
P
= 1
2 u
P 0
+ 1
2 Z
2
0
( 2
Q
cos+(3cos 2
1)u
Q cos
2
+3v
Q
sincos)d+O(t 2
)(3.8)
v
P
= 1
2 v
P 0
+ 1
2 Z
2
0
( 2
Q
sin+3u
Q
sincos+(3sin 2
1)v
Q sin
2
)d+O(t 2
):(3.9)
Inthissection we describethe evolution Galerkinschemesinthe nitedierenceframework
as well as the nite volume evolution Galerkin schemes. The main idea of the evolution
Galerkinschemesisthefollowing. Transportedquantitiesareshiftedalongthebicharacteris-
ticsandthenprojectedonto aniteelementspace. Thesemethodsconnect theniteelement
ideaswith the theory of bicharacteristics. In the nite volume framework the approximate
operators are used only in order to compute uxes on cell interfaces. Thus, instead of the
one-dimensionalRiemann solvers,whichwork onlyinthenormal directionsto thecellinter-
faces,wecomputetheapproximatesolutionatcellinterfacesbyamulti-dimensionalevolution
operator. Thiscan be consideredasa predictorstep. In thecorrectorstep thenitevolume
update ismade.
ConsiderameshinR 2
,whichconsists ofthesquare meshcells
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
;
where,h>0isthemeshsizeparameter,k;l2Z. LetusdenotebyH
(R 2
)theSobolevspace
of distributionswith derivatives up to the order in L 2
space, where 2 N. Consider the
generalhyperbolicsystemgivenbytheequation(2.1). LetusdenotebyE(s):(H
(R 2
)) p
!
(H
(R 2
)) p
the exact evolution operatorforthesystem(2.1), i.e.
U(:;t+s)=E(s)U(:;t): (4.1)
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 a 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
consider E
:(L
1
l oc (R
2
)) p
! (H
(R 2
)) p
to be a suitable approximate evolution operator for
E(). In practice we willuse restrictionsofE
to thesubspace S m
h
form 0. Then we can
dene thegeneralclass ofevolutionGalerkinmethods.
Denition 4.2 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
; (4.3)
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: (4.4)
We denote byR
h :S
m
h
!S r
h
a recovery operator, r >m 0 and considerour approximate
evolutionoperatorE
onS
r
h
. Inthispaperwewilllimitourconsiderationsto thecaseswhere
m = 0. In this case the integrals that we obtained from the projectionare evaluated either
exactly using the factthat the approximate values U n
are piecewiseconstants or bymeans
of some numericalquadratures. Usingpiecewiseconstantstheresultingschemes willonlybe
of rst order, even when E
is approximated to a higher order. Higher order accuracy can
be obtainedeitherbytakingm>0, orbyinsertinga recovery stage R
h
beforetheevolution
step inequation (4.3) to give
U n+1
=P
h E
R
h U
n
: (4.5)
This approach seemed to be hardly feasible for eective derivation and implementation of
higher order methods. A simplicationthatwe usedwastheapplication of themultidimen-
sionalevolutiononlyonthecellinterfaces. Thisleadsto thenitevolumeevolutionGalerkin
methods.
Denition 4.6 StartingfromsomeinitialdataU 0
2S m
h
,thenitevolumeevolutionGalerkin
method (FVEG) isrecursively dened by means of
U n+1
=U n
1
h Z
t
0 2
X
j=1 Æ
x
j f
j (
~
U n+
t
)d; (4.7)
where Æ
x
j f
j (
~
U n+
t
) represents an approximation to the edge uxdierenceand Æ
x
is dened
by Æ
x
v(x) = v(x+ h
2
) v(x h
2
). The cell boundary value
~
U n+
t
is evolved using the
approximate 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
; (4.8)
where
kl
is the characteristic functionof @
kl .
Formore detailson the higher order nitevolume evolution Galerkin (FVEG) schemes, see
[1 ], [6], [8], where the error analysis aswell as numerical experiments are presented. Using
the L 2
-projection (4.4), the approximate evolution operator E
, and (4.7), (4.8) the nite
dierence formulaeforboththeEGand theFVEG schemescan bewritten intheform
U n+1
kl
=U n
kl +
1
X
r=1 1
X
s= 1 C
rs U
n
k+rl +s
; (4.9)
where
C
rs
= 0
@
1
rs
1
rs
1
rs
2
rs
2
rs
2
rs
3
rs
3
rs
3
rs 1
A
: (4.10)
Heretheentries m
rs
; m
rs
; m
rs
,m=1;2;3aretakenappropriatelyaccordingtotheapproximate
evolutionoperatorE
used. IntheAppendixthestencilmatrices m
; m
and m
,m=1;2;3
arewrittenforsome EG schemes.
5 Basic tools
Aswe mentioned above ourstability considerations are based on Fourier analysis. We rst
recall some basicconcepts. Letf n
kl g
1
k;l = 1
bea two dimensionalsequence in`
2 .
Denition 5.1 The discrete Fourier transformation of f n
kl g 2 `
2
is the function
^n
2
L
2
h
;
h
2
dened by
^ n
=h 2
1
X
k= 1 1
X
l = 1 n
kl exp
ih(k+l )
:
andParseval'sequality.
Lemma5.2 (Inverse formula) If f n
kl g 2`
2 and
^n
is the discrete Fourier transformation
of f n
kl g, then
n
kl
= 1
4 2
Z
h
h Z
h
h
^ n
exp
ih(k+l )
dd:
Lemma 5.3 (Parseval'sequality)Iff n
kl g2`
2 and
^n
isthediscreteFouriertransformation
of f n
kl g, then
jj
^ n
jj=jj n
kl jj;
wheretherstnormistheL
2
-normon
h
;
h
h
;
h
andthesecondnormisthe`
2 -norm.
Hencewe have thefollowingresult.
Lemma5.4 Thesequencef n
kl
gisboundedin `
2
ifandonlyifthe sequencef
^ n
g isbounded
in L
2
h
;
h
h
;
h
.
Inorderto studythestabilityoflinearnumericalschemestheFouriertransformationisused.
Thisleadstothe estimationof thespectralradiusof theso-calledamplicationmatrix. The
spectralradiusof asquare complexmatrixA witheigenvalues
i
isdened to be
(A)=max
i j
i
j: (5.5)
Thespectral normof thematrixA isdened as
jjAjj=sup
x6=0 jjAxjj
jjxjj
: (5.6)
The norms on the right handside of equation (5.6) are the Euclideannorms of the vectors
Ax and x, respectively. Notethat forthe spectralnorm, asforanymatrixnorm, we always
have jjAjj(A).
6 Estimate of the stability limit
In [4 , Lemma 5.1] Lukacova et al. proved the following stability result forthe EG-schemes.
Thereexists
max
<1suchthattheEGschemesforthetwo-dimensionalwaveequationsystem
(3.1) are stable forany such that 0
max
, where =c t
h
is the CFLnumber. The
goal ofthissectionisto estimate precisely
max
fortheEG3 scheme. Analogouscalculations
can be donealso for other EG-schemes of type EG1-EG4 as wellas forthe FVEGschemes.
We applythediscreteFourier transformationon bothsidesofequation (4.9).
^
U n+1
=
^
U n
+h 2
1
X
k= 1 1
X
l = 1 1
X
r= 1 1
X
s= 1 C
rs U
n
k+rl +s
!
exp
ih(k+l )
: (6.1)
Bymaking thechangeof variablesk =k+r and l =l+swe get
h 2
1
X
k= 1 1
X
l = 1 1
X
r= 1 1
X
s= 1 C
rs U
n
k+rl +s
!
exp
ih(k+l )
= 1
X
r= 1 1
X
s= 1 C
rs exp
ih(r+s)
h 2
1
X
k 0
= 1 1
X
l 0
= 1 U
n
k 0
l 0
exp ih(k
0
+l 0
)
!
= 1
X
r= 1 1
X
s= 1 C
rs exp
ih(r+s)
^
U n
(6.2)
Thususingthisexpressionintheequation (6.1) we get
^
U n+1
= I+ 1
X
r= 1 1
X
s= 1 C
rs exp
ih(r+s)
!
^
U n
; (6.3)
whereI is theidentitymatrix. The coeÆcient of
^
U n
intheequation (6.3),
T(;) =I+ 1
X
r= 1 1
X
s= 1 C
rs exp
ih(r+s)
; (6.4)
is called theamplication matrix of the nite dierencescheme (4.9). Applying recursively
theresultof equation (6.3) n+1times yields
^
U n+1
= I+ 1
X
r= 1 1
X
s= 1 C
rs exp
ih(r+s)
!
n+1
^
U 0
=T n+1
(;)
^
U 0
: (6.5)
Wenote thatifjjT(;)jj1 then jj^
u n+1
jj jj^
u 0
jj,whichmeansthatthe f^
u n
g isL 2
-stable.
ConsidertheEG3 scheme,i.e.thenumericalscheme basedonequations(3.7) -(3.9), cf. also
stencilmatrices in the Appendix. After some calculation we obtain the following entries of
theamplicationmatrixT(;)
T
11
(;) = 1+
2
4
+
2
cos (h)cos(h)+
2
2
(cos (h)+cos (h));
T
12
(;) = i
4 2
3
sin(h)cos(h)+
4
2
3
sin(h)
;
T
13
(;) = i
4 2
3
cos (h)sin(h)+
4
2
3
sin(h)
;
T
22
(;) = 1 2
+
2
2 +
2
2
cos (h)cos(h)+
2
2
2
cos (h)
2
2
cos (h);
T
23
(;) = 3
2
8
sin(h)sin(h);
T
33
(;) = 1 2
+
2
2 +
2
2
cos (h)cos(h)+
2
2
2
cos (h)
2
2
cos (h);
T
21
(;) = T
12
(;); T
31
(;)=T
13
(;); T
32
(;)=T
23 (;):
Using the substitutions S
= sin(h), s
= sin(
2 ), S
= sin(h) and s
= sin(
2 ) the
amplicationmatrixT =T(;) can bewrittenas
T = 0
@ C
11
iC
iC
iC
C
22
2
C
iC
2
C
C
33 1
A
;
where
C
11
= 1 4
(s
2
+s
2
)+
4 2
s
2
s
2
;
C
= S
(1
8
3 s
2
);
C
= S
(1
8
3 s
2
);
C
= 3
8 S
S
;
C
22
= 1 4
s
2
+
2 2
s
2
s
2
;
C
33
= 1 4
s
2
+
2 2
s
2
s
2
:
SetE = 0
@
i 0 0
0 1 0
0 0 1 1
A
; then Q = 0
@ C
11
C
C
C
C
22
2
C
C
2
C
C
33 1
A
= E 1
TE; which means that
T andQ aresimilarmatricesandthustheyhavethesameeigenvalues. Moreover, thematrix
Q can be decomposedas
Q=I ( D+C)+ 2
~
C;
where
D =
0
@
d+f 0 0
0 d 0
0 0 f
1
A
; C= 0
@
0 C
C
C
0 0
C
0 0
1
A
;
~
C= 0
@
0 0 0
0 0 C
0 C
0
1
A
;
d = 4
s
2
2
s
2
s
2
= 2
s
2
(2 s 2
); f = 4
s
2
2
s
2
s
2
= 2
s
2
(2 s 2
):
Since
jjQ (I (D+C))jj 2
jC
j
3 2
8
=O(
2
): (6.6)
we can approximate the eigenvaluesof Q bythe eigenvalues of H =I (D+C) with the
O(
2
) error. Note that for all [;] 2 [
h
;
h ][
h
;
h
] the entries of H are bounded. To
estimate thestability limitof the EG3 scheme we needto estimate spectral radius of Q for
all choices of , and ,0 1. Since 0s 2
1 and 0s 2
1 and 1 thend 0
andf 0. Nowthematrices D,C arereal and C isantisymmetricaswell. HenceD+C has
eitherthree real eigenvaluesorone real eigenvalueand two complexconjugate eigenvalues.
Considerareal eigenvalue,say=
r
. Letv=(v
1
;v
2
;v
3
) bethecorresponding eigenvector,
thenv T
(D+C)v=v T
r
v . Since C isantisymmetricthen v T
Cv=0. Hence we get
(d+f
r )v
2
1
+(d
r )v
2
2
+(f
r )v
2
3
=0: (6.7)
0min(d;f)
r
d+f: (6.8)
Let
r
be a real eigenvalue of H then
r
= 1
r
. Hence j
r
j 1 is equivalent to 1
1
r
1. Usinginequality(6.8) we get
1 4
s
2
+s
2
+ 4
2
s
2
s
2
1
r 1:
To getj
r
j1weneed
1 4
s
2
+s
2
+ 4
2
s
2
s
2
1:
The lastinequalityreads
2
4
s
2
s
2
4
s
2
+s
2
+20: (6.9)
Now, wewant to bound suchthat 2 4
s 2
+s
2
0. Henceinequality(6.9) yields
2
4
s
2
s
2
4
s
2
+s
2
+2 2 4
s
2
+s
2
2
8
0:
The lastinequalitygives
4
0:7854: (6.10)
Now let usassumethat
c
isa complexeigenvalue ofH . Then
c
=1
c
,where
c is a
complexeigenvalue ofthe matrixD+C. This implies
j
c j
2
=1 2R e(
c )+
2
j
c j
2
:
Thusj
c j
2
1 isequivalent to 2
j
c j
2
2R e(
c
)0. Supposethat
r
>0then
2
r j
c j
2
2
r R e(
c
)0: (6.11)
Letb=C
andc=C
. It is wellknown that
det(D+C) = d 2
f+f 2
d+b 2
f+c 2
d=
r j
c j
2
;
Tr(D+C) = 2(d+f)=
r +
c +
c
=
r
+2R e(
c );
Henceinequality(6.11) reads
p(
r )=
2
r
2(d+f)
r +(d
2
f +f 2
d+b 2
f +c 2
d)0: (6.12)
Thediscriminantof p gives
2
=4(d+f) 2
4(d 2
f+f 2
d+b 2
f+c 2
d)=4(d 2
+f 2
)+8fd 4(d 2
f+f 2
d+b 2
f +c 2
d):
8fd 4(d 2
f+f 2
d+b 2
f +c 2
d)0;
whichleadsto 2
4(d 2
+f 2
)0 and we obtain
8fd= 32
2
s 2
s
2
(2 s 2
)(2 s 2
) =
32
2
s 2
s
2
(4 2(s 2
+s
2
)+
2
s 2
s
2
)
32
2
s 2
s
2
(4 2(s 2
+s
2
)):
Hence,
8fd 32
2
s 2
s
2
(4 4)= 128
2
s 2
s
2
(1 ): (6.13)
Furtherwe have
d 2
f = 8
3
s 2
s
2
(2 s 2
)
2
(2 s 2
)
64
3
s 2
s
2
64
3
;
b 2
f = S
2
(1
8 2
3 s
2
)
2 2
s
2
(2 s 2
)
4
S
2
s
2
4
:
Analogously,we obtain
f 2
d 64
3
and c
2
d 4
:
Therefore,
4(d 2
f+f 2
d+b 2
f+c 2
d) 4
128+8 2
3
: (6.14)
Putting inequalities(6.13) and (6.14) together we get
8df 4(d 2
f+f 2
d+b 2
f +c 2
d)
128
2
s 2
s
2
(1 ) 4
128+8 2
3
= 128
2
s 2
s
2
128
2
s 2
s
2
+4
128+8 2
3
0:
Thelastinequalityimplies
128
2
s 2
s
2
4
128+8 2
3
+ 128
2
s 2
s
2
128
2
4
128+8 2
3
+ 128
2
s 2
s
2
:
Since
4
128+8 2
3
+ 128
2
s 2
s
2
4
128+8 2
3
we thenhave
1
4
128+8 2
3
+ 128
2
s 2
s
2
1
4
128+8 2
3
:
Thereforewe get
128
2
4(128+8 2
)
3
=
32
128+8 2
0:4858: (6.15)
Thus we have obtained a suÆcient condition on for 2
0. For 0:4858 we have
2
4(d 2
+f 2
)0.
Now=0correspondstod +f =0,whichimpliesthat
r
=0,cf. (6.8). Sinceweassume
r
>0then p(
r
)has two distinctrealrootsr
1 and r
2
,where
r
1
=(d+f)
2
; r
2
=(d+f)+
2 :
Inequality(6.8) gives
r r
2
. To show that
r r
1
note thatfrom 4(d 2
+f 2
) we have
r
1
(d+f) p
d 2
+f 2
. Furthermore p
d 2
+f 2
max (d;f). Therefore
r
1
(d+f) p
d 2
+f 2
(d+f) max (d;f)=min(d;f)
r :
Hence
r 2[r
1
;r
2
]. Thisimpliesthatp(
r
)0,what wewantedto show, cf. (6.12).
Now considerthe case
r
= 0, then either d =0 or f = 0. Suppose d = 2
s
2
(2 s 2
) = 0
thens
=0 and
D+C= 0
@ 4
s
2
0 S
0 0 0
S
0
4
s
2
1
A
:
The eigenvaluesofD+C are0;
4
s
2
iS
. Now j
c j
2
=j1
c j
2
1 gives
1
4
s
2
+iS
1
4
s
2
iS
=
1 4
s
2
2
+ 2
S 2
!
= 1 8
s
2
+
16 2
2
s 4
+
2
S 2
1:
Thisleadsto
8
s
2
+
16
2
s 4
+S
2
0:
Supposes
6=0,otherwise kH k=1 forany ,thenwe have
8
+
16
2
+
S
s
2
!
0;
16
2
+
S
s
2
!
8
:
Now thelastinequalityyields
8
16
2
+
S
s
2
2
1:5708: (6.16)
Finallyinequalities(6.10), (6.15) and (6.16) implythat if
32
128+8 2
0:4858; (6.17)
thenthespectralradius ofthe matrixHis lessthanorequalto 1.
stability limit of the scheme EG3
0
32
128+8 2
0:4858
isless than or equal to 3
8
32
128+8 2
2
0:0885.
Hence we have proved thefollowingresult.
Lemma 6.18 Consider the evolution Galerkin scheme EG3. Then, this scheme is stable if
0
max
=
32
128+8 2
0:4858. Theerror corresponding to this estimation is less than
or equal to 3
8
32
128+8 2
2
0:0885.
Remark 6.19 In the next table we have estimated the stability limit of the scheme EG3
using the standard MATLAB procedure for the eigenvalues of the matrix T. Note that the
upper boundof ourtheoretical result, i.e.0.4858+0.0885=0.5743, matches very well with the
experimental stability 0.58.
ct
h
;
(T(;))forEG3
0.10 1.00000000000000 0
0.20 1.00000000000000 0
0.30 1.00000000000000 0
0.40 1.00000000000000 0
0.50 1.00000000000000 0
0.58 1.00000000000000 0
0.59 1.00000324446152 1
0.60 1.00011223611144 8
0.70 1.00847404969631 9
Table1: Stabilitylimitusing
;
(T(;)).
In Figure 2 left we plot the eigenvalues of the matrix H as well as the unit circle. Similar
plotwithdierentscaleis showninFigure2 right. Withthissimulationwe illustratethatit
ispossibleto includeall eigenvalues, possiblyexcept the eigenvalue 1, inside theunit circle.
Since the entries of the of the matrix H are bounded, the condition stated in (6.17) is a
necessary and suÆcient stabilitycondition, see Richtmyer and Morton [12]. In Figure 3 we
showa sequence ofplots, withdierent scales,ofthe eigenvaluesof theamplicationmatrix
correspondingto therstorderEG3 scheme. Againusingsimilarargumentweconcludethat
thescheme isstableto CFL=0.58.
7 Approximate evolution operator E
const
for piecewise con-
stant data
In[2 ]Lukacova,MortonandWarneckeproposednewapproximateevolutionoperatorsE const
undE bil in
forthe two-dimensionalwave equation systemand for theEuler equationsof gas
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
Eigenvalues of the matrix H, CFL=0.48
0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 1.02 1.04 1.06
0.2 0.25 0.3 0.35 0.4
Eigenvalues of the matrix H, CFL=0.48
Figure2: Eigenvaluesofthe matrixH ,CFL=0.48.
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
EG3 scheme, CFL=0.55
0.86 0.88 0.9 0.92 0.94 0.96 0.98 1
0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
EG3 scheme, CFL=0.55
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
EG3 scheme, CFL=0.58
0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 1.005
0.1 0.15 0.2 0.25
EG3 scheme, CFL=0.58
Figure 3: Eigenvalues of theamplication matrix of the rst order EG3 scheme, CFL=0.55
(top),CFL=0.58 (bottom).
improvethestabilityoftheFVEG-schemesconsiderably. Wewillshowthatforspecialchoices
of discretization techniques stability limits close to the natural limit of 1 can be achieved.
Numerical experiments, presented in [2], for these FVEG schemes conrm high accuracy
as well as good multidimensional behaviour of new FVEG schemes. The key idea of the
development of these new operators was to exploit the fact that the exact explicitsolution
to the one-dimensional wave equation system is available. Our new approximate operators
are constructedinsucha waythat thisexact solutionis reproducedexactlyfora given one-
dimensional data. Thus, the approximate evolution operator E const
calculates exactly any
one-dimensionalwave,whichisrepresentedbyapiecewiseconstantdataandpropagateseither
in x or y direction. Analogous situation holdsfor theoperator E bil in
and approximated
waves by means of continuous piecewisebilinear data. The approximate evolution operator
E const
for piecewiseconstant datareads,cf. [2 ]
P
= 1
2 Z
2
0 (
Q u
Q
sgncos v
Q
sgnsin)d; (7.1)
u
P
= 1
2 Z
2
0
Q
sgncos+u
Q
1
2 +cos
2
+v
Q
sincos
d; (7.2)
v
P
= 1
2 Z
2
0
Q
sgnsin+u
Q
sincos+v
Q
1
2 +sin
2
d: (7.3)
Integrationsfrom0to2 aroundthesoniccircle in(7.1)- (7.3)areevaluatedexactly. Inthis
way all innitelymanydirections of wave propagation are taken into account explicitly. For
thecellinterfaceintegration along@in(4.8)wehavetwo possibilities. Theseedgeintegrals
can either be computed exactly or numerically. Exact cell interface integration yields, e.g.
fortheverticaledge, thefollowingintermediate values
~
n+
1
2
edge
=
1+
2 Æ
2
y
x
n
1
2 +
4 Æ
2
y
Æ 2
x U
n
x
y Æ
y V
n
;
~
U n+
1
2
edge
=
1
2 +
4 Æ
2
y
Æ
x
n
+
1+ 5
12 Æ
2
y
2
x U
n
+
6 Æ
x
y Æ
y V
n
; (7.4)
where
x f(x)=
1
2
f x+ h
2
+f x h
2
; Æ 2
x
f(x)=f(x+h) 2f(x)+f(x h):
The stencilmatricesofthisFVEGscheme aregiven intheAppendix. Anotherpossibilityto
evaluatethecellinterfaceintegrals isto usesome numerical quadrature. Inthisway,further
simplication in the evaluation of integrals can be made. Instead of the two-dimensional
integrals along thecell interfaces and around the sonic circle, onlythe sonic circle integrals
needto beevaluatedexactly. Inourexperimentswe usedthetrapezoidalruleand Simpson's
ruleforthe cellinterfaceintegration. Thus, we needto determine
~
U n+
1
2
~
n+
1
2
vertex
=
x
y
n 1
2
y Æ
x U
n 1
2
x Æ
y V
n
;
~
U n+
1
2
vertex
= 1
2
y Æ
x
n
+
x
y U
n
+ 1
4 Æ
x Æ
y V
n
;
~
n+
1
2
midpoint
=
x
n 1
2 Æ
x U
n
;
~
U n+
1
2
midpoint
= 1
2 Æ
x
n
+
x U
n
: (7.5)
forthecellinterfaceare given inthe Appendixaswell.
Analogous to theSection 6, we can show that the amplicationmatrix T of the rst order
FVEGscheme withexact edge integrals issimilarto thematrix
Q=I ( D+C)+ 2
~
C;
wherethematrixD isdened asbefore with
d=2s 2
1 2
s
2
; f =2s 2
1 2
s
2
:
Thematrices C and
~
C aregiven as
C= 0
@
0 C
C
C
0 0
C
0 0
1
A
;
~
C= 0
@ 0
1
S
1
S
0 0 C
0 C
0
1
A
;
where
C
=S
1 2
s
2
; C
=S
1 2
s
2
;
C
= 1
S
S
:
Since the matrix C is symmetric it is now not possible to carry out the analysis similarto
theSection 6 inorder to estimatethe stabilitylimit. Instead we use a MATLABprocedure
to estimate the stability limit. The results are given in Table 2. In Column 2 we present
the stabilitylimit of the rst order FVEG scheme with exact edge integrals. The stability
limitof this scheme is improved considerably, thescheme is stable approximatelyup to the
CFL=0.89. Column3demonstratesthattherstorderschemebasedonthetrapezoidalrule
isstableto thenaturalstabilitylimit1. Column4 shows that thestabilityofthe rst order
scheme based on Simpson's ruleis also increased, the scheme is stable approximately up to
theCFL=0.75.
ct
h
Exact Trapezoidal Simpson
0.70 1.0000000000 1.0000000000 1.0000000000
0.75 1.0000000000 1.0000000000 1.0000000000
0.76 1.0000000000 1.0000000000 1.0206666667
0.80 1.0000000000 1.0000000000
0.89 1.0000000000 1.0000000000
0.90 1.0007993640 1.0000000000
1.00 1.0000000000
1.01 1.0200000000
Table 2: Stabilitylimitusing
;
(T(;))
In Figures 4 and 5 we plot, using dierent scales, the eigenvalues of the amplication ma-
trices corresponding to the rst order FVEG schemes based on the operator (7.1) - (7.3).
bottom the trapezoidal rule was used to approximate the interface integrals. Analogous to
theprevioussection,itis possibleto encircleall eigenvalues,with exceptionof 1,intheunit
circle. Since the entries of the amplication matrices are bounded the estimated stability
limitsarenecessaryconditionsand suÆcientaswell. InFigure 5we have plottedeigenvalues
ofthe amplicationmatrix ofthe FVEG3 scheme withSimpson's ruleapproximation of the
cellinterfaceintegrations.
8 Approximate evolution operator E bilin
for piecewise bilinear
data
Inthissectionweinvestigatethestabilityofthesecondordernitevolumeschemesproposed
by Lukacova et al. in [2]. These schemes are based on theapproximate evolution operator
E bil in
,which isgiven as
P
=
1
2
0
P +
1
2 Z
2
0
2
Q
2cosu
Q
2sinv
Q
d+O(t 2
); (8.1)
u
P
=
1
4
u 0
P +
1
2 Z
2
0
2cos
Q +
2 3cos
2
1
u
Q +
3
2
sincosv
Q
d (8.2)
+O(t 2
);
v
P
=
1
4
v 0
P +
1
2 Z
2
0
2sin
Q +
3
2
sincosu
Q +
2 3sin
2
1
v
Q
d (8.3)
+O(t 2
):
Analogous to E const
,thisapproximateevolution operatorisdesigned such thatit computes
any one-dimensional linear plane wave propagating in x or y directionexactly, for more
detailssee[2 ]. Inorder to obtainsecond ordernitevolumeschemeswe carryouta recovery
stage before applyingthe approximate evolution operator, see Denition 4.6. The following
twotypesofbilinearrecoveries havebeenconsideredin[2 ]
R C
h U
k l
=
2
x
2
y +
x x
k
h
x
2
y Æ
x +
y y
l
h
2
x
y Æ
y +
(x x
k )(y y
l )
h
x
y Æ
x Æ
y
U
k l
;(8.4)
R D
h U
k l
=
1+ x x
k
h
x
2
y Æ
x +
y y
l
h
2
x
y Æ
y +
(x x
k )(y y
l )
h
x
y Æ
x Æ
y
U
k l
: (8.5)
Note, that the recovery (8.4) is continuous while the recovery (8.5) is discontinuous and
conservative. We use the midpoint rule to approximate the time integral in the equation
(4.7). Denoting the cellinterfaceintermediate value, that iscomputed in thepredictor step
(4.8),byU n+
1
2
wecan obtainthefollowingschemes
schemeA U
n+
1
2
= E
bil in
R
C
h U
n
+E const
(1
2
x
2
y )U
n
;
scheme B U
n+
1
2
= E
bil in
R
C
h U
n
;
scheme C U
n+
1
2
= E
bil in
R
D
h U
n
:
Each of these schemes hasfurthertwo types according to theevaluation of thecellinterface
integrals. We used thesubscripts 1,2 to distinguishbetween them. Thus, 1 correspondsto
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
Eigenvalues: FVEG scheme based on operator (7.1−7.3), exact interface integrals,CFL=0.89
0.94 0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04
0.1 0.15 0.2 0.25 0.3
Eigenvalues: FVEG scheme based on operator (7.1−7.3), exact interface integrals,CFL=0.89
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
Eigenvalues: FVEG scheme based on operator (7.1−7.2), trapezoidal rule, CFL=0.9
0.88 0.9 0.92 0.94 0.96 0.98 1 1.02 1.04
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Eigenvalues: FVEG scheme based on operator (7.1−7.2), trapezoidal rule, CFL=0.9
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
Eigenvalues: FVEG scheme based on operator (7.1−7.3), trapezoidal rule, CFL=0.95
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
Eigenvalues: FVEG scheme based on operator (7.1−7.3), trapezoidal rule, CFL=1.0
Figure 4: Eigenvalues of the amplication matrices: top: exact interface integration
(CFL=0.89),middleand bottom: interfaceintegralsapproximatedusingthetrapezoidalrule
fortheCFL=0.9,0.95, 1.0.
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
Eigenvalues: FVEG scheme based on operator (7.1−7.3), Simpson’s rule, CFL=0.70
0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1 1.01 1.02
0.1 0.15 0.2 0.25 0.3
Eigenvalues: FVEG scheme based on operator (7.1−7.3), Simpson’s rule, CFL=0.70
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1
1.5 Eigenvalues: FVEG schemebased on operator (7.1−7.3), Simpson’s rule, CFL=0.75
0.92 0.94 0.96 0.98 1 1.02 1.04 1.06
0.05 0.1 0.15 0.2 0.25 0.3
Eigenvalues: FVEG schemebased on operator (7.1−7.3), Simpson’s rule, CFL=0.75
Figure 5: Eigenvalues of the amplication matrix, interface integrals approximated using
Simpson'srule,CFL=0.7 (top), CFL=0.75 (bottom).
2
valuesalong theright cellinterfaceare
~
n+
1
2
=
1+
32 +
16
Æ 2
x
2
y +
32 +
16
Æ 2
y
2
x +
32
8 +
2
32
Æ 2
x Æ
2
y
x
y
n
+
2
+
1
2
8
2
x
2
y +
1
8
16
Æ 2
y
2
x +
1
2 +
8 +
4
2
6
2
x Æ
2
y
Æ
x
y U
n
+
2
+
1
8
16
Æ 2
x
2
y +
1
2
8
2
x
2
y +
1
2 +
8 +
4
2
6
2
y Æ
2
x
Æ
y
x V
n
~
U n+
1
2
=
2
+
1
2
8
2
x
2
y +
1
8
16
Æ 2
y
2
x +
1
2 +
8 +
4
2
6
2
x Æ
2
y
Æ
x
y
n
+
1+
64 +
16
Æ 2
x
2
y
64 Æ
2
y
2
x +
64
16 +
2
64
Æ 2
x Æ
2
y
x
y U
n
+
1
8 +
1
16
8 +
2
64
2
x
2
y
3Æ
x Æ
y V
n
withtheequationfor
~
V n+
1
2
thatisanalogoustothatof
~
U n+
1
2
:Further,wecanexpressanalo-
gously thepredictedvaluesforothercellinterfacesaswellasforotherschemes. Substituting
thepredictedvaluesinthecorrectorstep(4.7)yieldsforallsecondordernitevolumeschemes
FVEG-A, B, C
U n+1
kl
= U
n
kl +
1
X
r= 1 1
X
s= 1 C
rs U
n
k+rl +s +C
x
rs U
n
x
k +rl+s +C
y
rs U
n
y
k +rl+s +C
xy
rs U
n
xy
k +rl+s
;
(8.6)
whereC x
rs ,C
y
rs andC
xy
rs
arethecoeÆcientmatricescorrespondingtotheapproximationofx ,
y ,and xy slopes. Moreover,
U n
x
k +rl+s
=
x
2
y Æ
x U
n
k+rl +s
; U n
y
k +rl+s
= 2
x
y Æ
y U
n
k+rl +s
; U n
xy
k +rl+s
=
x
y Æ
x Æ
y U
n
k+rl +s :
Applyingthevon Neumannanalysisand theFouriertransformationwederivetheamplica-
tionmatrices T:It shouldbepointedoutthattheir structureistoo complicated in orderto
applythesimilarestimatesofthespectralradiusaswedidinSection6fortherstorderEG3
scheme. Anyway,wecan usethestandardMATLABprocedure to determinetheeigenvalues
ofT:This yieldsthestabilitylimitsgiven inTable 3.
InFigures6and7weplot,usingdierentscales,theeigenvaluesoftheamplicationmatrices
corresponding to the second order FVEG schemes; scheme A
i , B
i
and C
i
, where i = 1;2.
Similar to the previous cases, these plots indicate that all eigenvalues, possibly except of
1, can be bounded inside the unit circle. Thus the stability limits which we obtained are
necessaryaswellassuÆcient conditions.
Further,it followsfrom Figure6 thatthesecond order FVEGscheme based onthe operator
(8.1) - (8.3) with the continuous non-conservative recovery (8.4) using Simpson's rule to
approximate the cell interface integrals, i.e. scheme B
1
, is unconditionally unstable. This
scheme A 0.94 0.75
scheme B 0.78 -
scheme C 0.78 0.58
Table 3: Stabilitylimitsofthesecond order FVEGschemes.
fact has also been conrmed by other numerical tests for the wave equation system with
discontinuous solution, see the Problem 3 in [4 ]. We have found for all CFL numbers, no
matterhow smallthey werechosen, instabilitiesinthe solutionforneenough meshes. We
should note that all other CFL limits given in the Table 3 have also been conrmed by
particular numericalexperiments.
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1
1.5 Eigenvalues: Scheme A1, CFL=0.75
0.9 0.95 1 1.05 1.1
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Eigenvalues: Scheme A1, CFL=0.75
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
Eigenvalues: Scheme B1, CFL=0.1
0.992 0.994 0.996 0.998 1 1.002 1.004 1.006
−0.1
−0.05 0 0.05 0.1
Eigenvalues: Scheme B1, CFL=0.1
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1
1.5 Eigenvalues: Scheme C1, CFL=0.58
0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 1.08
0.05 0.1 0.15 0.2 0.25 0.3
Eigenvalues: Scheme C1, CFL=0.58
Figure 6: Eigenvalues corresponding to the amplication matrices of the scheme A
1
(CFL=0.75), scheme B
1
(CFL=0.1), and thescheme C
1
(CFL=0.58).
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1
1.5 Eigenvalues: Scheme A2, CFL=0.94
0.85 0.9 0.95 1 1.05 1.1 1.15 1.2
−0.2
−0.15
−0.1
−0.05 0 0.05 0.1 0.15 0.2
Eigenvalues: Scheme A2, CFL=0.94
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1 1.5
Eigenvalues: Scheme B2, CFL=0.78
0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1
−0.3
−0.2
−0.1 0 0.1 0.2 0.3
Eigenvalues: Scheme B2, CFL=0.78
−1.5 −1 −0.5 0 0.5 1 1.5
−1.5
−1
−0.5 0 0.5 1
1.5 Eigenvalues: Scheme C2, CFL=0.78
0.5 0.6 0.7 0.8 0.9 1
0.05 0.1 0.15 0.2 0.25 0.3 0.35
Eigenvalues: Scheme C2, CFL=0.78
Figure 7: Eigenvalues corresponding to the amplication matrices of the scheme A
2
(CFL=0.94), scheme B
2
(CFL=0.78), scheme C
2
(CFL=0.78).
EG3scheme
1
:=
8
>
>
>
>
>
<
>
>
>
>
>
:
2
4
2
2
2
4
2
2 4
+
2
2
2
2
4
2
2
2
4 9
>
>
>
>
>
=
>
>
>
>
>
;
; 1
:=
8
>
>
>
>
>
<
>
>
>
>
>
:
2
3 0
2
3
2 2
2
3 0
2 +
2 2
3
2
3 0
2
3 9
>
>
>
>
>
=
>
>
>
>
>
;
;
1
:=
8
>
>
>
>
>
<
>
>
>
>
>
:
2
3
2 +
2 2
3
2
3
0 0 0
2
3
2 2
2
3
2
3 9
>
>
>
>
>
=
>
>
>
>
>
;
; 2
:=
8
>
>
>
>
>
<
>
>
>
>
>
:
2
3 0
2
3
2 2
2
3 0
2 +
2 2
3
2
3 0
2
3 9
>
>
>
>
>
=
>
>
>
>
>
;
;
2
:=
8
>
>
>
>
>
<
>
>
>
>
>
:
2
8
2
4
2
8
2
4 2
+
2
2
2
4
2
8
2
4
2
8 9
>
>
>
>
>
=
>
>
>
>
>
;
; 2
:=
8
>
>
>
>
>
<
>
>
>
>
>
: 3
2
32 0
3 2
32
0 0 0
3 2
32 0
3 2
32 9
>
>
>
>
>
=
>
>
>
>
>
;
;
3
:=
8
>
>
>
>
>
<
>
>
>
>
>
:
2
3
2 +
2 2
3
2
3
0 0 0
2
3 +
2 2
2
3
2
3 9
>
>
>
>
>
=
>
>
>
>
>
;
; 3
:=
8
>
>
>
>
>
<
>
>
>
>
>
: 3
2
32 0
3 2
32
0 0 0
3 2
32 0
3 2
32 9
>
>
>
>
>
=
>
>
>
>
>
;
;
3
:=
8
>
>
>
>
>
<
>
>
>
>
>
:
2
8
2
4
2
8
2
4 2
+
2
2
2
4
2
8
2
4
2
8 9
>
>
>
>
>
=
>
>
>
>
>
; :
