Rapid Evaluation of Nonreflecting Boundary Kernels for Time-Domain Wave Propagation.
We present a systematic approach to the computation of exact nonreflecting boundary conditions for the wave equation. In both two and three dimensions, the critical step in our analysis involves convolution with the inverse Laplace transform of the logarithmic derivative of a Hankel function. The main technical result in this paper is that the logarithmic derivative of the Hankel function $H_\nu^{(1)}(z)$ of real order $\nu$ can be approximated in the upper half $z$-plane with relative error $\varepsilon$ by a rational function of degree $d \sim O (\log|\nu|\log\frac{1}{\varepsilon}+ \log^2 |\nu| + | \nu |^{-1} \log^2\frac{1}{\varepsilon} )$ as $|\nu|\rightarrow\infty$, $\varepsilon\rightarrow 0$, with slightly more complicated bounds for $\nu=0$. If N is the number of points used in the discretization of a cylindrical (circular) boundary in two dimensions, then, assuming that $\varepsilon < 1/N$, $O(N \log N\log\frac{1}{\varepsilon})$ work is required at each time step. This is comparable to the work required for the Fourier transform on the boundary. In three dimensions, the cost is proportional to $N^2 \log^2 N + N^2 \log N\log\frac{1}{\varepsilon}$ for a spherical boundary with N2 points, the first term coming from the calculation of a spherical harmonic transform at each time step. In short, nonreflecting boundary conditions can be imposed to any desired accuracy, at a cost dominated by the interior grid work, which scales like N3 in two dimensions and N2 in three dimensions.
argument, z, satisfying Im(z)  0. The number of poles is bounded by O log ||

log 1 +log2 ||+||1 log2 1 . A similar representation for derived which
is valid for Im(z)  >0 requiring O log 1 log 1 +log 1 log log 1 +log 1 log log 1
poles.
In section 2, we introduce nonreflecting boundary kernels. In section 3 we collect
background material in a form convenient for the subsequent development. Section 4
contains the analytical and approximate treatment of the logarithmic derivative, while
a procedure for computing these representations is presented in Section 5. The results
of our numerical computations are contained in section 6, and we present our
conclusions in section 7.
2. Nonreflecting boundary kernels. Let us first consider the wave equation
in a two-dimensional annular domain 0 <<1. The general solution can be
expressed as
where Kn and In are modified Bessel functions (see, for example, [17, section 9.6]),
the coecients an and bn are arbitrary functions analytic in the right half-plane, L
denotes the Laplace transform
denotes the inverse Laplace transform
ds.
Likewise, for the wave equation in a three-dimensional domain r0 <r<r1, the
general solution can be expressed as

rs/c rs/c
If we imagine that is to be used as a nonreflecting boundary,
then we can assume there are no sources in the exterior region and the coecients
(or bnm(s)) are zero. Let us now denote by un(, t) the function satisfying
Then

un (, s)=an(s)   Kn (s/c)
s Kn (s/c)
c Kn(s/c)
so that

c Kn(s/c)
where  denotes Laplace convolution
convolution kernel in (2.9) is a generalized function. Its singular part is easily
removed, however, by subtracting the first two terms of the asymptotic expansion
c Kn(s/c) c 2
From the assumption un(, t)=0fort  0 and standard properties of the Laplace
transform we obtain the boundary condition
where

which we impose at
Remark. The solution to the wave equation in physical space is recovered on the
nonreflecting boundary from un by Fourier transformation:
assuming N points are used in the discretization.
The analogous boundary condition in three dimensions is expressed in terms of
the functions unm(r, t) satisfying
rs/c
After some algebraic manipulation, assuming unm(r, t)=0fort  0, we have
where

(rs/c)which we impose at
Note that the boundary conditions (2.12) and (2.16) are exact but nonlocal, since
they rely on a Fourier (or spherical harmonic) transformation in space and are history
dependent. The form of the history is simple, however, and expressed, for each separate
mode, in terms of a convolution kernel which is the inverse Laplace transform of
a function defined in terms of the logarithmic derivative of a modified Bessel function
d K (z)
log K(z)= .
dz K(z)
Remark. In three dimensions, the required logarithmic derivative of Kn+ 1 (z)isa ratio of polynomials, so that one can recast the boundary condition in terms of a
dierential operator of order n. The resulting expression would be equivalent to those
derived by Sofronov [7] and Grote and Keller [8].
The remainder of this paper is devoted to the approximation of the logarithmic
derivatives (2.18) as a ratio of polynomials of degree O(log ), from which the convolution
kernels n and n can be expressed as a sum of decaying exponentials. This
representation allows for the recursive evaluation of the integral operators in (2.12)
and (2.16), using only O(log n) work per time step (see [18]). We note that, by Par-
sevals equality, the L2 error resulting from convolution with an approximate kernel
is sharply bounded by the L error in the approximation to the kernels transform.
Precisely, approximating the kernel B(t) by the kernel A(t)wefind

A  u  B

where we assume that A, B, and u are all regular for Re(s) > 0. For finite times we
may let s have a positive real part, :

We therefore concentrate our theoretical developments on L approximations. For
ease of computation, however, we compute our rational representations by least
squares methods. These do generally lead to small relative errors in the maximum
norm, as will be shown.
Since Hankel functions are more commonly used in the special function literature,
we will write the logarithmic derivatives as

d d (1) i/2 H(1) zei/2
log K(z)= log H
dz dz H(1) zei/2
We are, then, interested in approximating logarithmic derivative of the Hankel function
on and above the real axis.
3. Mathematical preliminaries. In this section we collect several well-known
facts concerning the Bessel equation, the logarithmic derivative of the Hankel func-
tion, and pole expansions, in a form that will be useful in the subsequent analytical
development.
3.1. Bessels equation. Bessels dierential equation
dz2 z dz z2
(1) (2)
for   R, has linearly independent solutions H and H , known as Hankels
functions. These can be expressed by the formulae
J(z)  eiJ(z) J(z)  eiJ(z)
where the Bessel function of the first kind is defined by
z  (z2/4)k
The expressions in (3.2) are replaced by their limiting values for integer values of .
(1) (2)
(See, for example, [17, section 9.1].) For general , the functions H and H have
a branch point at and it is customary to place the corresponding branch cut
on the negative real axis and impose the restriction <arg z  . We shall find
it more convenient, however, to place the branch cut on the negative imaginary axis,
with the restriction
Hankels functions have especially simple asymptotic properties. In particular (see,
for example, [19, section 7.4.1]),
z zk
z zk 2z z
as z , with
(42  12)(42  32) (42  (2k  1)2)
and the branch of the square root is determined by
Finally we note the symmetry
We also make use of the modified Bessel functions K(z) and I(z). These are linearly
independent solutions of the equation obtained from (3.1) by the transformation
z  iz. Their Wronskian satisfies
Moreover we have for positive r [20]
Asymptotic expansions of K(r) and I(r) for r small and large are also known [17,
sections 9.6 and 9.7]. For real r and   0wehave
r
log ,=0,

Here  =0.5772 . is the Euler constant.
Finally, we note the uniform expansions of Bessel functions for  given in
[17]. For Hankel function and derivative we have


as  , where we restrict z to |arg(z)|/2 and define
3 z
Here, Ai(t) denotes the Airy function [17, section 10.4]. Note that
Large  approximations of the modified Bessel functions for real arguments, r, are
given by

where
r
3.2. Hankel function logarithmic derivative. We denote the logarithmic
(1)
derivative of H by G,
(1)
(3.21) G(z)= log H(1)(z)= .
dz  H(1)(z)
The following lemma states a few fundamental facts about G that we will use below.
Lemma 3.1. The function G(z), for   R, satisfies the formulae
z is the complex conjugate of z. Asymptotic approximations to
G are
where  is the Euler constant,
zk 2z z zk
(1)
Fig. 3.1. Curve z() defined by (3.18) near which the scaled zeros of H lie (see Lemma 3.2).
(1)
The branch cut of H is chosen (3.4) on the negative imaginary axis.
where Ak() is defined in (3.7), and

where  is defined in (3.18). Furthermore, the function u defined by
satisfies the recurrence
Proof. Equations (3.22) and (3.23) and asymptotic expansion (3.24) follow im-
(1)
mediately from the definitions (3.2) through (3.4) of J and H . The asymptotic
expansion (3.25) follows from (3.5) and (3.6), while (3.26) is a consequence of (3.16)
and (3.17). The recurrence (3.28) from standard Bessel recurrences [17, section
9.1.27].
(1)
The zeros of H (z) are well characterized [17, 20]; they lie in the lower half z-plane
near the curve shown in Figure 3.1, obtained by transformation [21] of Bessels
equation. In terms of the asymptotic approximation (3.16), this curve corresponds to
negative, real arguments of the Airy function.
(1)
Lemma 3.2. The zeros h,1,h,2,. of H (z) in the sector /2  arg z  0
are given by the asymptotic expansion
uniformly in n, where n is defined by the equation
z() is obtained from inverting (3.18), and an is the nth negative zero of Airy function
Ai. The zeros in the sector   arg z  3/2 are given by h,1, h,2,.In
particular,
where
3.3. Pole expansions. A set of poles in a finite region defines a function that
is smooth away from the region, with the smoothness increasing as the distance in-
creases. This fact leads to the following approximation related to the fast multipole
method [22, 23].
Lemma 3.3. Suppose that q1,.,qn are complex numbers and z1,.,zn are
complex numbers with |zj|1 for j =1,.,n. The function
z  zj
can be approximated for Re(z)=a>1 by the m pole expansion
z  j
is a root of unity and j is defined by
The error of the approximation is bounded by
where
z  zj
Proof. We use the geometric series summation
z  v zk+1 zm z  v
to obtain
zm z  zj z  j
All m terms of the first summation vanish, due to the combination of (3.34) and the
equality mk0. For the error term we obtain


Re
|z| (a  1)2 1+a2 (a  1)2 |z| 1  zj/z

(a  1)2 z  zj (a  1)2
and

Moreover, repeating the computations of (3.39), we find
(a  1)2
Now the combination of (3.38) through (3.41) and the triangle inequality gives
(3.35).
Inequality (3.35) remains valid if we assume instead that |zj|b and Re(z)=
ab>b, for arbitrary b  R, b>0; this fact leads to the next two results whose proofs
mimic that of Lemma 3.3 and are omitted.
Lemma 3.4. Suppose n, p are positive integers, q1,.,qn are complex numbers,
and z1,.,zn are complex numbers contained in disks D1,.,Dp of radii r1,.,rp,
centered at c1,.,cp, respectively. The function
z  zj
can be approximated for z satisfying Re(z ci)  ari >ri for i =1,.,pby the mp
pole expansion
ri j),
where ij is defined by
ri
with jiDj. The error of the approximation is bounded by
(am  1)(a  1)2
1148 BRADLEY ALPERT, LESLIE GREENGARD, AND THOMAS HAGSTROM
where
z  zj
Lemma 3.5. Suppose that the discrete poles of Lemma 3.4 are replaced with a
density q defined on a curve C with C  specifically

C z
which is finite for z outside Up, and that gm is defined by (3.43) with ij defined by
ri
with Then the bound (3.45) holds as before. Lemma 3.3 enables us
to approximate, with exponential convergence, a function defined as a sum of poles.
The fundamental assumption is that the region of interest be separated from the
pole locations. The notion of separation is eectively relaxed by covering the pole
locations with disks of varying size in an adaptive manner. In Lemmas 3.4 and 3.5,
we use this approach to derive our principal analytical result.
4. Rational approximation of the logarithmic derivative. The Hankel
functions logarithmic derivative G(z) defined in (3.21) approaches a constant as
z and is regular for finite z  C, except at z = 0, which is a branch point, and at
(1)
the zeros of H (z), all of which are simple. We can therefore develop a representation
for G analogous to that of the Mittag-Leer theorem; the only addition is due to
the branch cut on the negative imaginary axis. It will be convenient to work with
u(z), defined in (3.27), for which approximations to be introduced have simple error
bounds.
Theorem 4.1. The function u(z)=zG(z), where G is defined for   R by
(3.21) with the branch cut defined by (3.4), is given by the formula
for z  C not in {0,h,1,h,2,.,h,N } and not on the negative imaginary axis.
(1)
Here h,1,h,2,.,h,N denote the zeros of H (z), which number N.
Proof. The case of the spherical Hankel function, where
is simple and we consider it first. Here u(z) is a ratio of polynomials in iz with real
coecients, which is clear from the observation that u1/2(z)=iz1/2 in combination
with the recurrence (3.28). Hence
z  h,n
where p is a polynomial and ,n is the residue of u at h,n,
zh,n
(1)
Cm
(2)
Cm
Re(z)
Fig. 4.1. Integration contour Cm, with inner circle radius 1/m and outer radius m +1.
by lHopitals rule. We see from (3.25) .Noting that u(iy)  R for y  R, and combining (4.2), (4.3), and (4.5), we obtain
(4.1).
We now consider the case Z, for which the origin is a branch point.

For m =1, 2,., we define Cm to be the simple closed curve, shown in Figure 4.1,
(1)
which proceeds counterclockwise along the circle Cm of radius centered at the
origin from arg z = /2to3/2, to the vertical segment
(2)
to the circle Cm of radius 1/m centered at the origin from arg z =3/2to/2, to
the vertical segment z = rei/2, back to the first circle. Since none of the zeros of
(1)
H lies on the imaginary axis, Cm encloses them all if m is suciently large. For
(1)
such m, and z  C inside Cm with H the residue theorem gives

(1) (2)
We now consider the separate pieces of the contour Cm. For the circles Cm and Cm ,
we use the asymptotic expansion (4.4) about infinity and (3.24) about the origin to
obtain


Fig. 4.2. Plot of Re u(rei/2) , containing the zero crossing, and Im u(rei/2) , for
Now exploiting the symmetry u(re3i/2)=u(rei/2) from (3.23) for the vertical
segments, we obtain
which, when combined with (4.6), yields (4.1) and the theorem.
The primary aim of this paper is to reduce the summation and integral of (4.1)
to a similar summation involving dramatically fewer terms. To do so, we restrict z
to the upper half-plane and settle for an approximation. Such a representation is
(1)
possible, for the poles of u (zeros of H ) lie entirely in the lower half-plane and do
not cluster near the real axis. We first examine the behavior of u on the negative
imaginary axis.
The qualitative behavior of u on the branch cut is illustrated by the case of
shown in Figure 4.2. The plot changes little with changing , except for the

sign of Im u(z) and the sharpness of its extremum.
Lemma 4.2. For   R, Z, the function u(rei/2) is infinitely

dierentiable on r  (0, ) and has imaginary part satisfying the following formulae:
cos()
r2||,=0,

where  is defined in (3.20).
(1)
Proof. Infinite dierentiability of u(z) follows from the observation that H (z)
on the negative imaginary axis. To derive (4.9) we recall (3.11) to obtain

then apply (3.10). The remaining formulas follow from the asymptotic forms of K(r)
and I(r) for small and large r, and the uniform large  expansions given in (3.12)
through (3.15) and (3.19). Here we use the symmetry u. Note that (4.10) is
valid for r/||0. The approximation (4.12) is nonuniform for   2k  1/2 and
Lemma 4.3. Given 0 > 0 there exist constants c0 and c1 such that for all   R,
Z, and all z satisfying Im(z)  0, the function

satisfies the bounds
Moreover, there exists >0 such that for all   R, ||0, and  with 0 < <1/2,
f(z) admits an approximation g(z) that is a sum of d   1+||1 log(1/) log(1/)
poles, with
Proof. We assume and begin by changing variables,

so that
From the nonvanishing of z and its asymptotic behavior in w, it is clear that (4.15)
holds for ||(0,1) and any fixed 1 >0. Using (4.12) for || large but bounded
away from 2k 1/2 for integral k, an application of Watsons lemma to (4.14) focuses
on the unique positive zero, w,of defined in (3.20). As the derivative of this
function is positive, we conclude
cos()
where  is a function of w, so that (4.15) clearly holds. However, as   2k  1/2,
the denominator on the right-hand side of (4.12) may nearly vanish at w and the
expansion loses its uniformity. Setting cos()= in these cases, we see that the
denominator has a minimum which is bounded below by O(2). Hence in an O(||1)
neighborhood of the minimum which includes w,wehave
which by the change of variables is seen to satisfy the upper bound in
uniformly in . As the rest of the integral is small, the upper bound holds.
We now move on to the approximation. For a positive integer m and a positive

number w0, we define intervals
and
where f0,f1, and f2 are defined by the formulae
We will now choose w0 and m so that f0 and f2 can be ignored and then use Lemma
3.5 to approximate f1. Using (4.10) and (4.12) and taking w0 suciently small we
have, for some constant c2 independent of ,
(4.22) |f0(z)| w2||1dw  w0 .
Hence, a choice of
suces to guarantee
|f0(z)| |f(z)|in the closed upper half-plane. Now using (4.11) and (4.12) and assuming m su-
ciently large we have, for some constant c3 independent of ,

From (4.23), choosing
for appropriate m0 and m1 independent of  and  leads to
|f(z)|.Finally, we apply Lemma 3.5 to the approximation of f1. The error involves the
function but we note that Using p poles for
each j we produce a p  m-pole approximation g(z) with an error estimate, again for
Im(z)  0, given by(4.28) |f1(z)  g(z)| |f1(z)|.
A choice of
enforces
combining (4.24), (4.27), (4.30), and the triangle inequality, we obtain (4.16) with
the number of poles, satisfying the stated bound.
The case treatment. First, the direct application of the
preceding arguments leads to a significantly larger upper bound on the number of
poles. Second, we note that so that relative error bounds near z =0
require a vanishing absolute error. Finally, the lack of regularity of u0(z)atz =0
precludes uniform rational approximation, as discussed in [10]. Therefore, we relax the
condition Im(z)  0toIm(z)  >0. By (2.20) this will lead to good approximate
convolutions for times T  1.
Lemma 4.4. There exists >0 such that for all , 0 <<1/2 and , 0 <<
1/2, the function f(z)=u0(z)iz +1/2 admits an approximation g(z) that is a sum
of d    log(1/) log log(1/)  log(1/) poles, with
Proof. Note that since u0(z) has no poles, f(z) is given by (4.14) and satisfies
(4.15). Define intervals

Now
where f1 and f2 are defined by the formulae
We will now choose m so that f2 can be ignored and then use Lemma 3.5 to approximate
f1. Using (4.11) and assuming m suciently large we have, for some constant c,
Hence, choosing

log log(1/)
for appropriate m0 independent of  and  leads to
Finally, we apply Lemma 3.5 to the approximation of f1. Using p poles for each j we
produce a p  m-pole approximation g(z) with an error estimate for Im(z)   given
A choice of
enforces
(4.39), and the triangle inequality, (4.31) is achieved with the number of
poles, satisfying the stated bound.
We now consider the contribution of the poles.
Lemma 4.5. There exist constants C0, C1, >0 such that for all ,  R with
2 || and 0 <<1/2 the function
z  h,n
(1)
where h,1,.,h,N are the roots of H , satisfies the inequalities
and admits an approximation g(z) that is a sum of d    log ||log(1/) poles, with
Proof. The curve C defined in Lemma 3.2, near which h,1/||,.,h,N /|| lie,
is contained in disks separated from the real axis. If we denote the disk of radius r
centered at c by D(r, c), then the disks

for example, contain C\{+1, 1}. From (3.31), the root h,1 closest to the real axis
satisfies

hence it is contained in a disk of (4.43) with n  log2 24/331/2(a1)1||2/3 , and
all of the roots are contained in O(log ||)of the disks. Now applying Lemma 3.4 we
obtain (4.42) with |h| replaced by To obtain the upper
bound in (4.41) for both h and H we note first that it is trivial except for |z/|1.
A detailed analysis of the roots as described by Lemma 3.2 shows that

Hence, for |z/|1,
z  h,j
The lower bound in (4.41) is again obvious except for |z/|1. Then, however, we
note that
Since, from (3.26), by (4.15) the
right-hand side is dominated by iz and
The combination of Theorem 4.1 and Lemmas 4.3 and 4.5 suces to prove our
principal analytical result.
Theorem 4.6. Given 0 > 0 there exists >0 such that for all   R, ||0,
and 0 <<1/2 there exists d with

and complex numbers 1,.,d and 1,.,d, depending on  and , such that the
function
approximates u(z) with the bound

provided that Im(z)  0. Furthermore
Proof. We first note the lower bound
c||
(4.52) u(z)  iz +1/2  .
For >0 the function is nonvanishing and has the correct asymptotic behavior, so
we need only consider the case of || large. The result then follows from (3.26). This
proves (4.51) and (4.50) with u replaced by u  iz +1/2 on the right-hand side.
From (3.26) we have

so that the final result follows from the scaling  ||1/3.
The number of poles in (4.48) required to approximate u(z) to a tolerance
depends on both  and . The asymptotic dependence on  is proportional to
||1 log2(1/). We will see in the numerical examples, however, that this term is important
only for small ||; otherwise the dominant term is the first, for an asymptotic

dependence of O log ||log(1/) . As we generally have  ||1 in practice, the
term log2 || is of less importance.
Similarly, Lemma 4.4 leads to the following theorem for  =0.
Theorem 4.7. There exists >0 such that for all , 0 <<1/2 and ,

0 <<1/2 there exists d    log(1/)  log(1/) log log(1/)
and complex numbers 1,.,d and 1,.,d, depending on  and , such that the
function
approximates u0(z) with the bound

provided that Im(z)  . Furthermore
Proof. Again we already have (4.55) with u0(z)iz +1/2 on the right-hand side.
By (3.24) we find

log(1/)u0(z).
The theorem follows from the scaling   log1(1/).
As we must take we see that the number of poles required may grow
like log(1/)  log T log T  log log T. However, this is only for the mode in the
two-dimensionsal case. In short, the T dependence is insignificant in practice.
5. Computation of the rational representations. Analytical error bound
estimates developed in the previous sections are based on maximum norm errors
as in (2.19) and (2.20). In numerical computation it is often convenient, however, to
obtain least squares solutions. Our method of computing a rational function U, that
satisfies (4.50) is to enforce (4.51). An alternative approach would be to use rational
Chebyshev approximation as developed by Trefethen and Gutknecht [24, 25, 26].
In the numerical computations, we work with
and its sum-of-poles approximation U,(z)=U,(z)  iz +1/2. In particular, we
have the nonlinear least squares problem
for P,Q polynomials with deg(P) Problem (5.2) is not only
nonlinear, but also very poorly conditioned when P,Q are represented in terms of
their monomial coecients. We apply two tactics for coping with these diculties:
linearization and orthogonalization.
We linearize the problem by starting with a good estimate of Q and updating
P,Q iteratively. In particular, we solve the linear least squares problem
where the integral is replaced by a quadrature. The initial values P(0), Q(0) are
obtained by exploiting the asymptotic expansion (3.25) and the recurrence (3.28).
We find that two to three iterations of (5.3) generally suce.
The quadrature for (5.3) is derived by first changing variables,
where 1,.,m and w1 .,wm denote appropriate quadrature nodes and weights.
The transformed integrand is periodic on the interval [/2,/2], so the trapezoidal
rule (or midpoint rule) is an obvious candidate. The integrand is infinitely continously
dierentiable, except at its regularity is of order 2||.For|| > 8 (say),
the trapezoidal rule delivers at least 16th-order convergence and is very eective.
For small ||, however, a quadrature that adjusts for the complicated singularity at
needed. Here we can successively subdivide the interval near the singularity,
applying high-order quadratures on each subinterval (see, for example, [27]).
The quadrature discretization of (5.3) cannot be solved as a least squares problem
by standard techniques, due to its extremely poor conditioning. We avoid forming the
corresponding matrix; rather we solve the least squares problem by Gram-Schmidt
orthogonalization. The 2d
are orthogonalized under the real inner product
Ref(x)g(x)
to obtain the orthogonal functions
u(x),n=1,
where
Now

Table
Number d of poles to represent the Laplace transform of nonreflecting boundary kernels n and
n, for various values of .
=106
52- 86 11
52- 86 11
729-1024
15-
19- 22
26 17
ndnd46- 54 21
119- 145 26
178- 216 28
266- 324
487- 595 33
596- 728 34
19-
26-
37- 44 20
45- 53 21
119- 144 26
177- 216 28
265- 324
486- 594 33
595- 727 34
891-1024 36
=108
21- 28 11
59- 84 14
419- 638 19
20- 28 11
58-
419- 637 19
972-1024 21
ndnd
ndnd
so P(i+1) and Q(i+1) are computed from the recurrence coecients cnj by splitting
into even- and odd-numbered parts.
For some applications, including nonreflecting boundary kernels, it is convenient
to represent P/Q as a sum of poles,
P(z) d n
Q(z) z  n
We compute 1,.,d (zeros of Q) by Newton iteration with zero suppression (see,

Table
Laplace transform of cylinder kernel n defined in (2.13), approximated as a sum of d poles,
for n =1,.,4 and  =106.
nd
Pole Coecient
Re Im
Pole Location
Re Im 0.426478E  02
for example, [28]) by the formula

where 1,.,n1 are the previously computed zeros of Q. Then 1,.,d are
computed by the formula P(n)/Q(n). The derivative Q(z) is obtained by
dierentiating the recurrence (5.7).
6. Numerical results. We have implemented the algorithm described in section
5 to compute the representations of n and n through their Laplace transforms.
Recall that for the cylinder kernels, n, we have while for the sphere kernels, n,
we have presents the sizes of the representations for  =106,
108, and 1015 in (4.51). For the cylinder kernels, which are aected by the branch
cut, the number of poles for small n is higher than for the sphere kernels. This dis-
crepancy, however, rapidly vanishes as n increases and the asymptotic performance
ensues. The log(1/) dependence of the number of poles for
For  =108 we have also computed the maximum norm relative errors which
appear in (2.19) by sampling on a fine mesh. For the cylinder kernel with n =0,
we expect an O(1) error in a small interval about the origin due to (4.10). However,
errors of less than  are achieved for |s| > 5  107. This implies a similar accuracy
in the approximation of the convolution for times of order 106. For all other cases the
maximum norm relative errors are of order .
Finally, Table 2 presents poles and coecients for the cylinder kernels for
1,.,4 and  =106 to allow comparison by a reader interested in repeating our cal-
culations. Note that the pole locations are written in terms of Extensive tables
will be made available on the Web at http://math.nist.gov/mcsd/Sta/BAlpert.
Remark. Our approximate representation of the nonreflecting boundary kernel
could be used to reduce the cost of the method introduced by Grote and Keller [8]. The
dierential operators of degree n obtained in their derivation need only be replaced by
the corresponding dierential operators of degree log n for any specified accuracy. It
is interesting to note that in the two-dimensional case, where the approach of [8] does
not apply, the analysis described above can be used to derive an integrodierential
formulation in the same spirit.
7.

Summary

. In this paper we have introduced new representations for the
logarithmic derivative of a Hankel function of real order, that scale in size as the
logarithm of the order. An algorithm to compute the representations was presented
and our numerical results demonstrate that the new representations are modest in
size for orders and accuracies likely to be of practical interest.
The present motivation for this work is the numerical modeling of nonreflect-
ing boundaries for the wave equation, discussed briefly here and in more detail in
[18]. Maxwells equations are also susceptible to similar treatment as outlined in [29].
The new representations enable the application of the exact nonreflecting boundary
conditions, which are global in space and time, to be computationally eective.
8.

Appendix

Stability of exact and approximate conditions. In this ap-
pendix, we consider the stability of our approach to the design of nonreflecting boundary
conditions. Given that we are approximating the exact conditions uniformly, it
is natural to expect that our approximations possess similar stability characteristics.
This is, indeed, the case. Oddly enough, however, the exact boundary conditions
themselves do not satisfy the uniform Kreiss-Lopatinski conditions which are necessary
and sucient for strong well-posedness in the usual sense [30]. This may seem
paradoxical since the unbounded domain problem itself is strongly well-posed. The
diculty is that the exact reduction of an unbounded domain problem to a bounded
domain problem gives rise to forcings (inhomogeneous boundary terms) which live in
a restricted subspace. The Kreiss-Lopatinski conditions, on the other hand, require
bounds for arbitrary forcings. In that setting, our best estimates result in the loss of
1/3 of a derivative in terms of Sobolev norms. In practice we doubt that this fact is
of any significance, and have certainly encountered no stability problems in our long
time numerical simulations.
To fill in some of the details, consider a spherical domain  of radius one, within
which the homogeneous wave equation with homogeneous initial data is satisfied. At
the boundary we have
unm skn
r
for the exact condition and is uniformly small when we use our approx-
imations. Here gnm is the spherical harmonic transform of an arbitrary forcing g.
Following Sakamoto, we seek to estimate
where

while 0, denotes the usual L2 norm. On the boundary, , we will make use
of fractional Sobolev norms, most easily defined in terms of the spherical harmonic
coecients:

Strong well-posedness would follow from showing that
g(,t)20,dt.Instead, we can show that
dt.
1/3,To prove this, let .Bounded solutions within the sphere are given by
Precisely, setting
we find
where
(1)
We now estimate norms of the solution. First note that the products in the definition

(1) (1)
of n, J(z)H (z), zJ(z)H (z), are uniformly bounded for Im(z)  0. (See the
limits z  0, z , and  .) Therefore, as mentioned above, the error term,
so long as its small, has no eect on the estimates we derive, and we simply ignore
it. That is, we set n =1.
We concentrate on the boundary terms in H, as they are both the most straight-forward
to compute and the most ill behaved. In transform space we have

(Here and throughout, c will denote a positive constant independent of all variables.)
We first note that as the only singularities of Bessel functions occur at zero and
infinity, we need only consider the limits z  0, z , and  . The first two
are straightforward:


For large  we use the uniform asymptotic expansions of Bessel functions due to Olver
which yield
z
From Parsevals relation, we conclude that
g(,t)1/3,dt.
The estimation of the spatial integrals is more involved, as for r<1 the solution has
two transition zones, z   and rz  , and there are a number of cases to consider.
However, the estimates follow along the same lines and lead to the same result.
It is interesting to note that the loss-of-derivative phenomenon is suppressed when
one looks at the error due to the approximation of the boundary condition. In that
case the transform of the exact solution near the boundary is
(1)
hn (rz)
(1)
hn (z)
so that the error, e, satisfies the problem above with gnm given by

(1)
zhn (z)
(1)
hn (z)
Now the best estimate of n takes the form
which, in combination with (8.6), would lead to an estimate of the 1-norms of the
error in terms of the 4/3-norms of the solution. However, using again the large
asymptotics, a direct calculation shows
Thus n is smaller than its maximum by O(1/3) in the transition region where
O(1/3). Hence we find for the error
In other words, the 1-norms of the error are controlled by the 1-norms of the solution.
We have, of course, ignored discretization error, which could conceivably cause
diculties in association with the lack of strong well-posedness. To rule them out
would require a more detailed analysis. In practice we have encountered no diculties,
even for very long time simulations. We should also note that strong well-posedness
could be artificially recovered by perturbing the approximate conditions for large n.
allowing high accuracy to be maintained for smooth solutions. Finally, we note that
a similar analysis can be carried out in two dimensions.


--R


On the accurate long-time solution of the wave equation in exterior domains: Asymptotic expansions and corrected boundary conditions
On high-order radiation boundary conditions
Artificial boundary conditions of absolute transparency for two- and three-dimensional external time-dependent scattering problems
Fast discrete polynomial transforms with applications to data analysis for distance transitive graphs
A fast transform for spherical harmonics

Modulus and phase of the reduced logarithmic derivative of the Hankel function
Handbook of Mathematical Functions
Nonreflecting Boundary Conditions for the Time-Dependent Wave Equation
Asymptotics and Special Functions
The asymptotic expansion of bessel functions of large order
The asymptotic solution of linear di
THOMAS HAGSTROM
A fast algorithm for particle simulations
An implemenation of the fast multipole method without multipoles
Rational Chebyshev approximation on the unit disk
Real and complex Chebyshev approximation on the unit disk and interval
Rational Carath
On the Numerical Solution of One-Dimensional Integral and Dierential Equa- tions
Introduction to Numerical Analysis
Accurate boundary treatments for Maxwell
Hyperbolic Boundary Value Problems
--TR

--CTR
Laurence Halpern , Olivier Lafitte, Dirichlet to Neumann map for domains with corners and approximate boundary conditions, Journal of Computational and Applied Mathematics, v.204 n.2, p.505-514, July, 2007
Marcus J. Grote , Christoph Kirsch, Nonreflecting boundary condition for time-dependent multiple scattering, Journal of Computational Physics, v.221 n.1, p.41-62, January, 2007
Andreas Atle , Bjorn Engquist, On surface radiation conditions for high-frequency wave scattering, Journal of Computational and Applied Mathematics, v.204 n.2, p.306-316, July, 2007
Bjrn Sjgreen , N. Anders Petersson, Perfectly matched layers for Maxwell's equations in second order formulation, Journal of Computational Physics, v.209 n.1, p.19-46, 10 October 2005
John Visher , Stephen Wandzura , Amanda White, Stable, high-order discretization for evolution of the wave equation in 1 Journal of Computational Physics, v.194 n.2, p.395-408, March 2004
Bradley Alpert , Leslie Greengard , Thomas Hagstrom, Nonreflecting boundary conditions for the time-dependent wave equation, Journal of Computational Physics, v.180 n.1, p.270-296, July 20, 2002
Isaas Alonso-Mallo , Nuria Reguera, Discrete absorbing boundary conditions for Schrdinger-type equations: practical implementation, Mathematics of Computation, v.73 n.245, p.127-142, January 2004
Jing-Rebecca Li , Leslie Greengard, High order marching schemes for the wave equation in complex geometry, Journal of Computational Physics, v.198 n.1, p.295-309, 20 July 2004
S. V. Tsynkov, Artificial boundary conditions for the numerical simulation of unsteady acoustic waves, Journal of Computational Physics, v.189 August
Stephen R. Lau, Rapid evaluation of radiation boundary kernels for time-domain wave propagation on blackholes: theory and numerical methods, Journal of Computational Physics, v.199 n.1, p.376-422, 1 September 2004
S. I. Hariharan , Scott Sawyer , D. Dane Quinn, A Laplace transorm/potential-theoretic method for acoustic propagation in subsonic flows, Journal of Computational Physics, v.185 n.1, p.252-270, 10 February
Dan Givoli , Beny Neta, High-order non-reflecting boundary scheme for time-dependent waves, Journal of Computational Physics, v.186 n.1, p.24-46, 20 March
high-order convergence;wave equation;approximation;radiation boundary condition;nonreflecting boundary condition;absorbing boundary condition;maxwell's equations;bessel function