Simulations of Acoustic Wave Phenomena Using High-Order Finite Difference Approximations.
Numerical studies of hyperbolic initial boundary value problems (IBVP) in several space dimensions have been performed using high-order finite difference approximations. It is shown that for wave propagation problems, where the wavelengths are small compared to the domain and long time integrations are needed, high-order schemes are superior to low-order ones. In fact, in two dimensions an acoustic lens is simulated, leading to large scale computations where high-order methods and powerful parallel computers are necessary if an accurate solution is to be obtained.
Introduction
In this paper we present results of numerical simulations using high-order finite difference
methods. In particular we are interested in hyperbolic wave propagation
problems where accurate long time integrations are needed. For such problems
high accuracy in both time and space are needed if an accurate solution is to be
obtained. The efficiency of high-order methods has been studied by Kreiss and
Oliger in [4] and by Swartz and Wendroff in [12]. These studies are done on a
scalar hyperbolic model problem with periodic boundary conditions. In [7] and
numerical experiments on hyperbolic IBVP show the efficiency of high-order
methods compared with low-order ones. In this paper we will use two different
types of high-order finite difference approximations for hyperbolic IBVP. The first
method uses difference operators that satisfy a summation by parts rule. The numerical
boundary conditions are built in into the operator. To treat the analytic
boundary conditions a projection operator is used [8], [9]. For the second method
numerical boundary conditions are obtained by extrapolation of the outgoing characteristic
variables and by using the analytic boundary conditions and the PDE
for the ingoing characteristic variables [3], [11].
In section 2 and 3 these difference methods are used to solve hyperbolic IBVP.
In section 2 the efficiency of high-order approximations is compared to second-order
approximations on a scalar model equation. The number of points per
wavelength and points per period needed to obtain a certain accuracy for long
Department of Scientific Computing, Uppsala University, Uppsala, Sweden
time integrations is studied. These results are then compared with corresponding
results by Kreiss and Oliger [4] and Swartz and Wendroff [12] where periodic
boundary conditions are used.
In section 3 we simulate plane waves that are refracted when going through a
convex lens. The refracted waves then intersects at a point called the focus. This
can be simulated by solving the two-dimensional wave equation in a channel with
a plane wave entering the domain. The convex lens is modeled with a variable
velocity of sound in the channel. To be able to distinguish the focus it is important
to have high accuracy in both space and time to make errors from dispersion small.
1.1 Summation by Parts and Projections
We discretize by dividing the x-axis into intervals of length h, and use for
the notation
d is a grid vector function. The basic idea is to approximate @=@x
with a discrete operator Q, satisfying a summation by parts rule
where
N ); are grid vectors and
In reality Q is a banded matrix with centered difference stencils in the interior and
one-sided difference stencils near the boundaries. The scalar product and norm is
defined as
where h is the grid spacing and H is the norm matrix having the form
I
The existence of such operators has been proved by Kreiss and Scherer [6], [5]. In
[10] Strand has used their theory to construct such operators. These operators
are distinguished with respect to the form of the norms for the summation by
parts rule. Here we will use two different kinds of norms, diagonal norms and
restricted full norms. The corresponding difference operators are called diagonal
norm difference operators and restricted full norm difference operators. For the
diagonal norm the matrices H (1) ; H (2) are diagonal and for the restricted full norm
they have the form
are scalars.
To define difference operators in two space dimensions we denote by
functions and by Q 1 and Q 2 the difference operators
approximating @=@x 1 and @=@x 2 . The difference operators are written as
are the grid spacings in the x 1 - and x 2 -direction, and the q's are
defined such that summation by parts and accuracy conditions holds. In [8] it is
showed that the difference operators defined by (2) satisfy a summation by parts
rule in two dimensions if the summation by parts norm is a diagonal norm.
In [8], [9] Olsson has proved stability for high-order approximations of hyperbolic
and parabolic systems by using such difference operators. To treat the
analytic boundary conditions in a correct way they are represented as an orthogonal
projection. In several space dimensions stability can only be proved if the
difference operators satisfy a summation by parts rule in a diagonal norm. For
the restricted full norms such results can not be proved. However, in section 3
numerical computations of the two-dimensional wave equation show that for this
problem the restricted full norm difference operator results in a stable scheme.
In the following we will refer to this type of method as the summation by parts
method.
1.2 Strongly Stable Approximations
For the second type of method used here, the extra numerical boundary conditions
needed to close the high-order finite difference approximations near the
boundaries, are obtained by extrapolation of the outgoing characteristic variables,
and by differentiating the analytic boundary conditions and using the PDE for the
ingoing characteristic variables. This technique was first proposed by Gustafsson,
Kreiss and Oliger [3] who proved strong stability for a fourth-order approximation
of systems of hyperbolic PDEs in one space dimension. This scheme was later
generalized to general order of accuracy 2r by Strand [11]. Strong stability means
that an estimate of the solution at any given time is obtained in terms of the
forcing function, initial data and boundary data. In several space dimensions this
technique is applicable although no stability results have been proved. We will
refer to this type of method as the strongly stable method.
2 Efficiency of High-Order Methods
In this section we will compare the efficiency of high-order finite difference approximations
with second-order ones when long time integrations are needed. This is
done by computing the errors of the numerical approximations using the methods
mentioned in section 1 when applied to a scalar hyperbolic IBVP. It is then shown
that with non-periodic boundary conditions, there is no restriction on the number
of points per period and points per wavelength needed to obtain a certain error,
compared with the case of periodic boundary conditions. In the periodic case we
refer to classical results by Kreiss and Oliger [4], for semi-discrete approximations,
and Swartz and Wendroff [12] for fully discrete approximations.
2.1 Periodic Case
Here we present some of the results in [4] and [12] regarding the periodic problem
We approximately solve the problem by
dt
where
r
is a centered difference operator. The solution to this system is
!h
and has the phase error
Let P denote the number of periods we want to compute in time , the
wavelength and -
h the number of points per wavelength. Then
c and we
define
c
as the phase error per period.
For the fully discrete approximation we proceed as in [12] except that we here
use the standard fourth-order Runge-Kutta method in time. The system (6) can
be written as
dt
j be the approximate solution of the differential equation (3-5) obtained by
using the fourth-order Runge-Kutta method on (8). That is,
where \Deltat is the time step. Then the solution at
where
Denote by
!c\Deltat the number of time intervals per period. With
the exact solution to (3), we have
where we used the stability of the Runge-Kutta method, jRj - 1, in the third step.
By using
Thus from (10) we conclude
As in [12] an error function, valid for small errors, is now defined as

Table

1-3 show the error e defined by (12) for second-, fourth-, and sixth-order
approximation in space combined with the fourth-order Runge-Kutta method in
time. The error is shown for various number of points per wavelength, -, and CFL
h . The number of points per period can then be expressed as
-=CFL. Furthermore we note that the approximation (9) is stable for CFL -
2.85, 2.06, 1.79 for second-, fourth-, and sixth-order approximation in space.
Here we are interested in very long time integrations and table 1-3 show the
theoretical estimates on the error for up to 900 periods. These estimates will be
used later to compare with numerical results when non-periodic boundary conditions
are used.

Table

1: e, for second-order centered difference approximation of u t +cu

Table

2: e, for fourth-order centered difference approximation of u t +cu
Table

3: e, for sixth-order centered difference approximation of u t +cu
2.2 Non-Periodic Boundary Conditions
Consider the problem
where
is the wavelength and 10-. The solution, given
by u(x; is a wave packet, containing 10 wavelengths, traveling with
speed one along the x-axis.
We now want to solve (13) numerically by using the methods presented in
section 1. We will here use the method of lines approach by first defining semi-discrete
approximations and then using the standard fourth-order Runge-Kutta
method for integration in time. We discretize as in the periodic case and turn to
the semi-discrete approximations of (13).
2.2.1 Summation by Parts Method
Denote by vector and use the following discrete form of
the analytic boundary condition
To treat the analytic boundary condition in a correct way we will use a method
by Olsson [8], [9] where the boundary condition is represented as an orthogonal
projection with respect to the scalar product
(\Delta; \Delta) h . In [8] it is shown that the projection becomes
where I is the identity matrix and H the summation by parts norm. A semi-discrete
approximation is then defined as
~
where P is the projection operator and Q a difference operator that satisfy a
summation by parts rule in the scalar product (\Delta; \Delta) h . In [8] it is proved that if
the initial condition satisfies the boundary condition this approximation is stable.
This is obviously true for our problem and we have a stable scheme.
2.2.2 A Strongly Stable Method
For the second method mentioned in section 1 we add extra grid points, x
r. The differential equation (13)
is approximated for by a centered finite difference scheme of order
dt
where Q is the difference operator defined by (7). In [11] it is proved that with
the numerical boundary conditions
at the inflow boundary and
at the outflow boundary, the approximation (19)-(21) is strongly stable and the
error of the solution is of order h 2r . By using these numerical boundary conditions
to modify the difference operator close to the boundary we can write the
approximation (19)-(21) as an ODE system
where ~
Q is the (N \Theta N)-matrix obtained from the difference
operator Q and the boundary conditions (20), (21), and F is a vector containing
g(t) and its derivatives.
2.2.3 Numerical Results
Here we have used the fourth-order Runge-Kutta method to integrate the semi-discrete
systems (18) and (22) in time for 450 and 900 periods. This is done for
different number of points per wavelength, -, and different CFL numbers. We
assume that is the number of grid points in the interval
is the wavelength. For the summation by parts method we have
used a diagonal norm difference operator of order three at the boundary and six
in the interior, and a restricted full norm operator of order three at the boundary
and four in the interior. We will refer to these operators as the fourth-order
diagonal norm operator (D4) and the fourth-order restricted full norm operator
(RF4). Here the order of the difference operator refers to the order of the global
accuracy that the theory by Gustafsson [1], [2] predicts. There it is proved that
boundary conditions of at least order must be imposed to retain pth-order
global accuracy. For the strongly stable method we have used a second-, fourth-,
and sixth-order scheme in space.
Let u denote the exact solution and u h the approximate solution and define by
the error at a given point. Thus, e 1 corresponds to the theoretical
estimate of the error in the periodic case, e.

Table

4-6 displays the error e 1 when 450 and 900 periods in time have been
computed and figure 1-3 shows approximate and exact solution for the strongly
stable schemes. Comparing with the periodic case, table 1-3, we see that e in
general is larger then the computed error e 1 . An explanation for this is that we
maybe made a too rough estimate in the third step in (10) where we replaced
jR(i\Deltatfi)j with one. For large \Deltatfi we will have jR(i\Deltatfi)j ! 1 and the error will
be less than what is predicted by the estimate. Furthermore, the error function
(12) is only valid for small phase errors e f . When we increase the number of grid
points per wavelength the agreement is better and (12) can be used to predict the
number of grid points per wavelength needed, even when the boundary conditions
are non-periodic.
To compare the efficiency of the schemes we define as in [12] the work per
period per wavelength to obtain a certain error by
where d is the number of operations per mesh point per time step. d, of course
depends on the spatial accuracy. In [12] the optimal values of - and M were
obtained by minimizing w for a fixed error. Here we will use table 4-6 to obtain an
optimal combination (-; CFL;w) for a given error. These tables do not in any way
cover all combinations (-; CFL) and we will not claim that we found the optimal
combination over all (-; CFL).
As table 4 shows the errors for the second-order scheme are large even for
140 points per wavelength. If we compute for 450 periods with the strongly stable
schemes and want to have e 1 - 0:85, optimal schemes are for second-order
(140,2,245000), for fourth-order (20,1,13200), and for sixth-order (10,0.5,8200).
Thus, the fourth-order scheme are about times more efficient than the second-order
0one, and the sixth-order scheme about times more efficient. It should
be pointed out that this is a very large error and yet the high-order schemes are
much more efficient than the second-order one. If we compute for 900 periods
the errors for the second-order scheme are larger than one even for 140 points per
wavelength. If we compare the high-order schemes for an error e 1 - 0:65, optimal
schemes are for fourth-order (25,1,20625) and for sixth-order (20,1,16400), i.e. the
sixth-order scheme is 1.2 times more efficient than the fourth-order scheme. Table
5 clearly shows that the efficiency of the sixth-order scheme compared with the
fourth-order scheme increases with decreasing error. Here we have only considered
the strongly stable schemes but the schemes satisfying a summation by parts rule
gives the same errors for the same number of points per wavelength and points
per period.

Table

4: for the second-order strongly stable scheme for CFL=1.0 and
CFL=2.0.

Table

5: strongly stable schemes for CFL=0.5 and
CFL=1.0.
fourth-order sixth-order
Table

summation by parts schemes for CFL=0.5 and
CFL=1.0.
fourth-order fourth-order
diagonal restricted full
-o- approx
exact
-o- approx
exact
-o- approx
exact
-o- approx
exact

Figure

1: Approximate solution after 450 periods obtained with the second-order
-o- approx
exact
-o- approx
exact
-o- approx
exact
-o- approx
exact

Figure

2: Approximate solution after 450 periods obtained with the fourth-order
strongly stable scheme.
-o- approx
exact
-o- approx
exact
-o- approx
exact
-o- approx
exact

Figure

3: Approximate solution after 450 periods obtained with the sixth-order
3 The Acoustic Lens Problem
As an application where high-order accurate approximations are needed we will
simulate plane acoustic waves that are refracted to a focus when they are traveling
through a convex lens. This can be simulated by solving the two-dimensional wave
equation in a channel
where have the components pressure and velocity in the x 1 - and
-direction and
The velocity of sound c(x 1 defined by
simulates a convex lens. At x we have the boundary
cu) is the local inflow component and
ae sin 2 ( -
a
\Gamma3 is the wavelength and a = 10P , where
the period. This boundary condition simulates a incoming plane wave traveling to
the right. At x set the inflow component to zero, i.e.2 (p \Gamma cu)(1; x
At
The problem is symmetric with respect to the x 1 -axis and for the numerical
computations we will use this symmetry and only compute the solution in the
domain 1=15. The symmetry conditions are given for
3.1 Numerical Methods
To solve the problem (23) numerically we first define semi-discrete approximations
by using the summation by parts method and the strongly stable method.
To integrate the semi-discrete approximations in time the standard fourth-order
Runge-Kutta method is used.
We discretize in space and leave time continuous by dividing the domain in
intervals of length h 1 and h 2 in the x 1 - and x 2 -directions. For
we use the notation
is a vector grid function.
3.1.1 Summation by Parts Method
Here we will use the difference operators defined in section 1.1. However, we modify
the difference operator approximating @=@x 2 , at grid points near x
we let the interior scheme and the symmetry conditions (28) define the difference
operator instead of the one sided difference stencils. To define a semi-discrete
approximation of (23) we proceed as in [8]. We express the analytic boundary
conditions as
@\Omega and
ae
and
Define a grid vector by u
The boundary conditions are discretized and written as
ij
with the non-zero element being the ith entry and ~
1. The boundary conditions can then be written as
~
~
In [8] it is shown that when the summation by parts norm is a diagonal norm the
projection operator (17) is independent of the norm, i.e.
A semi-discrete approximation is then defined as in [8]
where A 1 and A 2 here denotes the grid matrix representations of A fl
are the difference operators approximating @=@x 1 and
with
In [8] stability is proved for systems of the type (32) when the difference operator
satisfy a summation by parts rule with respect to a diagonal norm. Here we
will also use restricted full difference operators. In this case the operator does not
satisfy a summation by parts rule in two dimensions and a projection operator can
not be defined. However, in one dimension the projection operator is the same
for diagonal and restricted full norm difference operators and we therefore use
the projection operator defined by (31) also for the restricted full norm difference
operators.
3.1.2 A Strongly Stable Method
Here we need extra grid points outside the domain 0 - x 1 -
we will apply centered difference stencils of order 2r at the first interior points.
Therefore we add grid points, x 1 i
in the x 1 -direction and x 2 j
in the x 2 -direction.
We approximate (23) by a centered difference approximation of order 2 r
ij and
r
To derive extra boundary conditions for the inflow components we regard the
velocity locally constant at the boundaries x
since the derivatives of c are small. At x need 2r conditions on p; u
and at x need 2r conditions on p; v. At x we note that
the boundary conditions (24), (26) do not depend on x 2 , therefore v(j; x
By differentiating the boundary
conditions with respect to time and using the differential equation we have for
@x k1
where g (k) (t) is the k-th derivative of g(t). At x the boundary condition
gives us v The differential equation,
Furthermore, by p x2
But, us
Integration of (36) with respect to t and the initial condition implies
Thus we have arrived at the condition
Repeating this procedure we get at x
@ 2k
We approximate (34) and (35) for even k's, thus for
where D+ fl
are the usual forward and backward differences in the the x
direction. At x
(D+2
For the locally outgoing characteristic variables we use extrapolation to derive
extra conditions. Thus for
Finally, at x we use the symmetry conditions (28), i.e. for
Summing up, to integrate the ODE system (33) in time we proceed as follows.
Given the solution, use (38),
(40) to compute p; u at the boundary and outside the domain at x
(39), (40) to compute p; v at the boundary and outside the domain at x
At the symmetry conditions (41) is used to compute p; v outside the
domain. This will give us p at all points except at the two corner points
, u at all points except at the boundary x
and v at all points except at the boundaries x 1. Thus, all points that need
to be specified for the integration of (33) are well defined.
3.2 Simulations on a Parallel Computer
The solution to our problem is a plane wave packet traveling essentially in the
-direction. The focus will occur at the end of the channel and the waves have to
be propagated through the whole channel. The sharpness of the focus will depend
on the number of wavelengths that can be contained in the x 2 -direction. When
the wavelength goes to zero the focus will be a point with infinite amplitude.
By numerical experiments we have found that at least 30 wavelengths in the x 2 -
direction are needed. In our case we have which means that the size
of the domain is 500 wavelengths in the x 1 -direction, and 33.3 wavelengths in the
-direction. Thus, to propagate the wave to the focus we have to compute up to
500 periods in time.
Since the wavelength is very small compared to the size of the domain and
the number of periods we have to compute is large, we need a powerful computer
with large memory to be able to solve this problem. For example, from the one-dimensional
problem we needed at least 10 points per wavelength for a fourth- or
sixth-order scheme when we computed 450 periods. In the x 2 -direction we do not
need that many points, since in that direction the solution is smoother. Thus,
if we chose 10 points per wavelength in the x 1 -direction and 3 points in the x 2 -
direction as a lower limit on the number of points needed, we will have 1:5
unknowns. Since h 1 - h 2 we define the CFL number as the quotient between
the time step and the space step in the x 1 -direction. Thus, for CFL= 0:5 we
have to do 10000 iterations in time. However, since the solution is a wave packet
traveling to the right we only compute those points in the x 1 -direction where the
solution is non-zero. Numerical results show that about 50-100 wavelengths in the
-direction have to be computed, giving about 0:15 unknowns
in each iteration. This is a lower limit on the number of unknowns and it grows
rapidly if we increase the number of points per wavelength.
The implementation was made on a 96 processors SP2 distributed memory
computer from IBM, at the Center for Parallel Computers at KTH, equipped with
66.7 MHz RS/6000 processors, giving a theoretical peak performance per processor
of 266 MFlop/s and 26 GFlop/s for the whole machine. The program was written
in Fortran 90 with MPI (Message Passing Interface) for communication between
the processors.
The computations are made in a moving window, figure 4, that follows the
wave through the domain. The size of the window is 100-200 wavelengths in the
-direction and the whole domain in the x 2 -direction. In the window only those
points in the x 1 -direction where the absolute value of the solution is larger than
are computed. As shown in figure 4 the window is divided into p strips
oriented in the x 1 -direction and processor q in the parallel computer is assigned
to
subdomain\Omega q . Since we use explicit difference operators in space and explicit
time integration on a uniform grid, the program is trivially parallel. The main
communication between the processors is to exchange p and v, needed to compute
-derivatives, at the interior boundary between two sub-domains.
In figure 5 contour plots of p 2 at different times are displayed, obtained by
integrating (32) in time with the standard fourth-order Runge-Kutta method and
the fourth-order diagonal norm difference operator in space, with
denotes the number of points per wavelength in the x fl -direction.
The wave packet is refracted and the focus occur after about 450 periods, i.e after
0.9 sec. To find the exact location of the focus this scheme and the sixth-order
strongly stable scheme were run with and the location of the
focus was defined as the time, t, where This happened at the
same time, schemes and the solutions, p 2 , is shown in figure 6.

Figure

7-9 show the focus, p 2 at different number of points per
wavelength, obtained with the summation by parts method and the strongly stable
method, combined with the fourth-order Runge-Kutta method in time. All
schemes work very well although stability results in two dimensions exists only
for the summation by parts method combined with a diagonal norm difference
operator. To compare these schemes an error, e, is defined as the number of
wavelengths that the center of the focus for the approximative solutions is behind
the center of the focus for the "exact" solutions obtained above. The center
of the focus is defined as follows. Let - x denote the point on the x 1 -axis where
occured. The center of the focus is then defined as the
on the x 1 -axis, where - x
x is defined as the largest
point where
x is defined as the
smallest point where
In table 7 the error, e,
and the consumed CPU time is compared for the high-order schemes. All computations
in table 7 were made on 10 processors. First of all one can note that
the consumed CPU time for the strongly stable schemes are lower than for the
summation by parts schemes. This is because of different implementations. A
notable observation is that for the same number of points per wavelength, the
sixth-order schemes often need less CPU time to reach the focus compared with
the fourth-order schemes. The explanation for this is that with few grid points per
wavelength the dispersion for sixth-order schemes are lower than for fourth-order
schemes and we do not need to compute as many wavelengths in the x 1 -direction
with the sixth-order schemes as with the fourth-order schemes. By comparing
the CPU time needed to obtain a certain error we see that for this problem the
sixth-order schemes are much more efficient than the fourth-order schemes. In this
application the boundaries have a relatively small influence on the solution, and
the scheme (D4) behaves as a sixth-order scheme, which table 7 shows. This is
because it is the sixth-order stencils in the interior that are significant. In figure
7 the focus is computed with a second-order scheme with
25. The CPU time needed to compute these solutions was in the
former case 5.1 hours on 20 processors, and in the latter case 10.8 hours on processors
.The corresponding errors of the focus was 2.36 and 0.91 wavelengths. By
comparing the latter solution with the strongly stable schemes, with
we see that although they have a smaller error and only 10 processors were
used, they require about 11 and 23 times less CPU time than the second-order
scheme needs on processors. For a smaller error the efficiency of the high-order
schemes compared to the second-order scheme should be much more pronounced.
In fact, for small errors of order 0.1, say, it is very doubtful if one can afford to
use second-order methods even on powerful parallel computers.
Conclusions
Numerical studies on hyperbolic IBVP in one and two space dimensions have been
performed using two different high-order finite difference methods. For the summation
by parts method a fourth-order diagonal norm difference operator, (D4),
and a fourth-order restricted full norm difference operator, (RF4), was used. For
the strongly stable method fourth- and sixth-order difference operators was used.
In one space dimension the number of points per wavelength and points per period
needed to obtain a certain accuracy was studied on a scalar hyperbolic IBVP. It
was shown that with non-periodic boundary conditions there is no restriction on
the number of points per wavelength and points per period needed to obtain a
certain error compared with the periodic case. In two space dimensions we simulated
plane waves that was refracted to a focus when going through a convex lens.
The problem was solved by an implementation, using Fortran 90 and MPI, on a 96
processors SP2 distributed memory computer. All schemes worked well although
stability results only exists for the summation by parts scheme combined with a
diagonal norm difference operator. The efficiency of high-order schemes compared
with second-order schemes was demonstrated by comparing the CPU time needed
to compute the focus. In this application the sixth-order schemes turned out to
be the most efficient ones. The reason for this is that for few grid points per
wavelengths it has lower dispersion than fourth-order schemes which means that
fewer wavelengths in the x 1 -direction need to be computed.

Table

7: The error, e, and consumed CPU hours, T cpu , for high-order schemes and
different computations were made on 10 processors.
D4 RF4 fourth-order sixth-order
strongly stable strongly stable

Acknowledgments

I would like to thank my adviser Professor Bertil Gustafsson for helpful support.
Also, I would like to thank Docent Leif Abrahamsson for many stimulating and
helpful discussions.
computational domain are between
broken lines
propagating
wave
computation in
moving window
-x 0x 20
e
computed points are between
s and x 0
e

Figure

4: The computations are done in a moving window. The window is divided
into p strips oriented in the x 1 -direction and the processor q in the parallel
computer is assigned the
subdomain\Omega q .
x
y
x
y
100 periods0.70.28 0.300.07
x
y
x
y
200 periods1.00.48 0.500.07
x
y
x
y
300 periods2.30.68 0.700.07
x
y
x
y
x
y
periods

Figure

5: Contour plots of p 2 at different times for the fourth-order diagonal norm
scheme with
6th strongly stable
y
D4 summation by parts

Figure

Focus obtained by solving the wave equation with the fourth-order
diagonal norm scheme and the sixth-order strongly stable scheme with
y
100x25 points/wavelength
y

Figure

7: Focus obtained by solving the wave equation with the second-order
diagonal norm scheme with
y
RF4, 10x3 points/wavelength
y
RF4, 20x5 points/wavelength
y
y
D4, 20x5 points/wavelength

Figure

8: Focus obtained by solving the wave equation with high-order summation
by parts schemes for different number of points per wavelength and
4th, 10x3 points/wavelength
y
4th, 20x5 points/wavelength
y
6th, 10x3 points/wavelength
y
6th, 20x5 points/wavelength

Figure

9: Focus obtained by solving the wave equation with high-order strongly
stable schemes for different number of points per wavelength and



--R

The convergence rate for difference approximations to mixed initial boundary value problems.
The convergence rate for difference approximations to general mixed initial boundary value problems.
Time Dependent Problems and Difference Methods.
Comparison of accurate methods for the integration of hyperbolic equations.
Finite element and finite difference methods for hyperbolic partial differential equations.
On the existence of energy estimates for difference approximations for hyperbolic systems.

Summation by parts
Summation by parts
Summation by parts for finite difference approximations for d/dx.

The relative efficiency of finite difference methods.
--TR

--CTR
Bernhard Mller , H. C. Yee, Entropy Splitting for High Order Numerical Simulation of Vortex Sound at Low Mach Numbers, Journal of Scientific Computing, v.17 n.1-4, p.181-190, December 2002
hyperbolic initial boundary value problems;wave propagation;high-order finite difference;acoustics