SIAM
J.
NUMER.
ANAL.
Vol.
33,
No.
5,
pp.
1797-1843,
October
1996
()
1996
Society
for
Industrial
and
Applied
Mathematics
006
CONVERGENCE
OF
A
BOUNDARY
INTEGRAL
METHOD
FOR
WATER
WAVES*
J.
THOMAS
BEALEf,
THOMAS
Y.
HOU
t,
AND
JOHN
LOWENGRUB
Abstract.
We
prove
nonlinear
stability
and
convergence
of
certain
boundary
integral
methods
for
time-dependent
water
waves
in
a
two-dimensional,
inviscid,
irrotational,
incompressible
fluid,
with
or
without
surface
tension.
The
methods
are
convergent
as
long
as
the
underlying
solution
remains
fairly
regular
(and
a
sign
condition
holds
in
the
case
without
surface
tension).
Thus,
numerical
instabilities
are
ruled
out
even
in
a
fully
nonlinear
regime.
The
analysis
is
based
on
delicate
energy
estimates,
following
a
framework
previously
developed
in
the
continuous
case
[Beale,
Hou,
and
Lowengrub,
Comm.
Pure
Appl.
Math.,
46
(1993),
pp.
1269-1301].
No
analyticity
assumption
is
made
for
the
physical
solution.
Our
study
indicates
that
the
numerical
methods
must
satisfy
certain
compatibility
conditions
in
order
to
be
stable.
Violation
of
these
conditions
will
lead
to
numerical
instabilities.
A
breaking
wave
is
calculated
as
an
illustration.
Key
words,
water
waves,
boundary
integral
methods,
surface
tension
AMS
subject
classifications.
Primary,
65M12;
Secondary,
76B
15
1.
Introduction.
In
this
paper
we
prove
convergence
of
certain
boundary
integral
meth-
ods
for
time-dependent,
two-dimensional
water
waves,
that
is,
waves
on
the
surface
of
an
inviscid,
incompressible
fluid
in
irrotational
motion
under
the
influence
of
gravity
and
pos-
sibly
surface
tension.
Boundary
integral
methods
have
been
used
extensively
for
calculating
the
motion
of
fluid
interfaces.
A
partial
list
of
works,
discussed
further
below,
includes
[And,
BMO,
BN,
BS,
CS,
Do,
DY,
GAS,
Kerr,
Kral,
Kra2,
LC,
Moor,
NMP,
Pull,
RS,
Rob,
Shel,
SS,
Tryg,
VB,
WBJ,
Ye].
In
boundary
integral
methods,
the
interface
is
tracked
by
markers
which
move
with
the
fluid.
While
there
are
many
variants,
a
key
feature
is
that
only
quantities
on
the
surface
need
to
be
computed
because
of
the
irrotationality.
For
water
waves,
the
velocity
field
is
expressed
as
a
singular
integral,
determined
by
the
position
of
the
interface
and
one
more
scalar
quantity,
which
could
be
the
velocity
potential,
dipole
sheet
strength,
or
vortex
sheet
strength.
The
formulation
used
here
is
closely
related
to
that
of
[BMO,
Bak].
Numerical
instabilities
have
been
encountered
and
dealt
with
in
computations,
and
linear
analysis
at
equi-
librium
has
identified
sources
of
instability
[Rob,
Do,
BN],
but
there
has
been
no
complete
error
analysis
of
boundary
integral
methods.
(An
exception
is
the
ill-posed
case
of
vortex
sheets
for
which
the
solutions
must
be
assumed
analytic
[CL,
HLK].)
The
rigorous
proof
of
convergence
of
the
particular
numerical
methods
for
water
waves
presented
here
means
that
numerical
instabilities
are
ruled
out
even
in
the
fully
nonlinear
regime,
as
long
as
the
solution
remains
fairly
regular.
A
filtering
in
high
wave
numbers
is
used
in
evaluating
the
singular
integrals
occurring
in
the
time
step
in
order
to
maintain
the
numerical
stability;
the
methods
are
designed
so
that
accuracy
is
preserved
even
with
this
filtering.
A
calculation
of
a
breaking
wave
without
surface
tension
is
presented
as
an
example
of
the
capability
of
these
methods
(see
Figures
6-12).
Error
analysis
for
boundary
integral
methods
is
difficult,
especially
far
away
from
equilib-
rium,
because
of
the
nonstandard
nature
of
the
equations.
Our
convergence
proof
uses
energy
*Received
by
the
editors
November
11,
1993;
accepted
for
publication
(in
revised
form)
November
7,
1994.
tMathematics
Department,
Box
90320,
Duke
University,
Durham,
NC
27708
(beale@math.duke.edu).
The
research
of
this
author
was
supported
by
NSF
grant
DMS-9102782.
Research
at
the
M.S.R.I.
was
supported
in
part
by
NSF
grant
DMS-8505550.
tDepartment
of
Applied
Mathematics,
California
Institute
of
Technology,
Pasadena,
CA
91125
(hou@
ama.caltech.edu).
The
research
of
this
author
was
supported
by
a
Sloan
Foundation
Research
Fellowship,
NSF
grant
DMS-9003202,
and
in
part
by
AFOSR
grant
AFOSR-90-0090
and
NSF
grant
DMS-9100383.
Department
of
Mathematics,
University
of
Minnesota,
Minneapolis,
MN
55455
(lowengrb@math.umn.edu).
The
research
of
this
author
was
supported
by
an
NSF
Postdoctoral
Fellowship,
the
University
of
Minnesota
Army
High
Performance
Computing
Research
Center,
and
the
Minnesota
Supercomputer
Center.
1797
1798
J.T.
BEALE,
T.
Y.
HOU,
AND
J.
LOWENGRUB
estimates
in
discrete
Sobolev
spaces
with
norms
measuring
several
derivatives.
In
order
to
maintain
numerical
stability,
balances
among
terms
with
singular
integrals
and
derivatives
have
to
be
taken
into
account.
Examples
indicate
that
poorly
chosen
discretizations
can
indeed
lead
to
numerical
instabilities.
Our
analysis
follows
a
framework
developed
in
the
continuous
case
for
linearized
motion
perturbed
about
an
arbitrary
exact
solution
[BHL1].
A
balance
of
important
terms
was
observed
in
the
continuous
equations.
However,
new
difficulties
arise
at
the
discrete
level.
The
discretization
must
not
introduce
new
instabilities
in
the
high
modes,
and
this
requires
that
certain
compatibility
conditions
be
imposed
in
the
way
various
terms
in
the
equations
are
discretized.
For
this
reason,
the
choice
of
rules
for
integration
and
differen-
tiation,
as
well
as
the
placement
of
the
filtering,
is
interdependent
and
affects
whether
or
not
the
numerical
method
is
stable.
We
now
state
the
water
wave
problem
in
a
form
which
leads
naturally
to
the
numerical
method
of
interest.
This
formulation
is
very
close
to
that
of
[BMO]
for
the
more
general
case
of
an
interface
separating
two
fluids.
For
simplicity
we
assume
that
the
fluid
has
infinite
depth.
The
fluid
interface
is
parametrized
by
a
complex
variable
z(ot,
t)
at
the
time
t.
The
parameter
ot
is
the
Lagrangian
variable.
Further,
we
denote
by
4
(ot,
t)
the
velocity
potential
on
the
interface,
and
the
real
part
and
the
imaginary
parts
of
z
as
x
and
y,
i.e.,
z(ot,
t)
x(ot,
t)
-t-
iy(ot,
t).
To
obtain
a
system
of
evolution
equations,
we
need
to
express
the
velocity
potential
in
the
fluid
domain
in
terms
of
these
two
variables.
Following
[BMO]
and
[BHL1],
we
express
the
complex
potential
by
a
double-layer
representation.
Denote
by/x(ot,
t)
the
dipole
strength
to
be
determined
from
4.
We
can
write
the
complex
potential
in
the
fluid
domain
in
terms
of
/z
as
the
principal
value
integral
2rri
z
z(ot’,
t)
I(ot’,
t)
dz(ot’).
The
real
velocity
potential
in
the
interior
is
0
G(z
z(ot’,
t))lz(ot’,
t)
ds(ot’),
b
ReO
On(ot’)
where
G(z)
(27r)
-1
log
Izl,
The
normal
derivative
is
taken
outward
from
the
fluid
domain.
The
nonnormalized
vortex
sheet
strength
t’
is
given
as
the
Lagrangian
derivative
of
the
dipole
strength,
i.e.,
9/
=/z
0/x/0ot.
For
simplicity,
we
often
drop
the
time
variable
from
now
on,
but
all
the
quantifies
z,/z,
q,
and
7’
will
be
time-dependent.
It
follows
from
the
Plemelj
formula
of
complex
variables,
or
the
properties
of
the
double-layer
potential,
that
the
value
of
4
on
the
interface
is
given
by
1
1
/’
1
b(ot)
/z(ot)
+
Re-/J
z(ot)
z(ot’)
’t(ot’)
dz(ot’)
or
f
(ot)
-lZ(ot)
+
on(oti
G(z(ot)
z(ot’))/z(ot’)
ds(ot’).
Differentiating
both
sides
of
the
4
equation
with
respect
to
ot
and
integrating
by
parts,
we
obtain
q(ot)
t’
-t-
Re
2yri
z(ot)
z(ot’)
The
complex
velocity
w
u
iv
can
be
obtained
by
differentiating
the
complex
potential
with
respect
to
z.
We
get
dO
f
1
w
dz
27ri
z
z(ot’)
?/(ot’)dot’.
A
CONVERGENT
NUMERICAL
METHOD
FOR
WATER
WAVES
1799
Using
the
Plemelj
formula
one
more
time,
we
obtain
the
velocity
on
the
interface
f
,
(or’)
dot’
-
V
(or)
w(ot)
z(ot)
z(ot’)
2z(ot------’
where
we
have
taken
the
limit
of
z
approaching
the
free
surface
z(ot)
from
the
fluid
region.
Since
ot
is
a
Lagrangian
coordinate,
the
velocity
of
the
interface
is
that
of
the
fluid
below,
and
we
obtain
an
evolution
equation
for
the
interface
__Oz
(or,
t)
w*
(or,
t),
Ot
where
the
asterisk
denotes
the
complex
conjugate.
For
the
evolution
of
(ot,
t),
we
use
Bernoulli’s
equation.
If
we
neglect
surface
tension,
the
pressure
is
zero
at
the
interface.
(The
case
with
surface
tension
included
will
be
treated
later.)
Thus,
Bernoulli’s
equation
in
the
Lagrangian
frame
is
1
Ct
lwl
2
+
gY
O.
From
now
on,
with
z(ot,
t)
ot
+
s(ot,
t),
we
assume
that
s(ot,
t)
and
(ot,
t)
are
periodic
in
ot
with
period
27r.
This
implies
that
the
flow
is
at
rest
at
infinity.
The
singular
kernel
1
cot(Z
in
the
velocity
integrand
now
becomes
g
).
To
summarize,
we
obtain
a
system
of
time
evolution
equations
for
z
and
as
follows:
1
f_
(z(ot)-z(c’))dot’_
?’(or)
=w(ot,
t)
(1)
z
47ri
y(ot’)
cot
2
2z(ot)
(2)
(3)
Ct
(u
2+v
2)-g.y,
[z
f_
(z()-z(’))d’]
Y
(od)
cot
+Re
r
2
Equations
(1)-(3)
completely
determine
the
motion
of
the
system.
A
unique
solution
is
specified
by
giving
initial
conditions
for
the
interface
position
z
and
the
velocity
potential
.
It
can
be
shown
that
the
integral
equation
for
the
vortex
sheet
strength
9/
can
be
solved
in
terms
of
[BMO].
This
is
done
at
each
time;
we
then
use
the
interface
equation
(1)
and
the
Bernoulli
equation
(2)
to
update
z
and
.
In
order
to
use
the
system
(1)-(3)
for
a
numerical
algorithm
we
need
to
make
choices
for
a
discrete
derivative
operator
and
a
quadrature
rule.
In
addition,
we
use
a
filtering
of
the
interface
location
z.
These
choices
must
be
made
in
conjunction.
We
find
it
natural
to
use
a
filtering
related
to
the
Fourier
symbol
of
the
derivative
operator.
To
describe
these
choices
further,
we
recall
that
the
discrete
Fourier
transform
is
defined
by
N/2
27r
tk---’
u(vtj)e
-ikaJ
otj--jh,
h=.
j=-N/2+I
N
The
inversion
formula
reads
u
(oj)
N/2
k=-N/2+l
800
J.T.
BEALE,
T.
Y.
HOU,
AND
J.
LOWENGRUB
Here
we
assume
N
is
even.
We
note
that
tk
is
periodic
in
k
with
period
N.
A
discrete
derivative
operator
may
be
expressed
in
the
Fourier
transform
as
(4)
(D
f)k
ikp(kh)
f
k,
k
N2
+1
--’N
We
will
always
assume
p
(zr)
0;
this
is
important
in
some
of
the
arguments
to
follow.
The
choice
of
p
(x)
depends
on
the
derivative
operator.
For
the
second-order
centered
difference
operator,
we
have
sin(kh)
P2(kh)
>
O;
kh
for
the
fourth-order
centered
differencing,
8
sin(kh)
sin(2kh)
p4(kh)
>_
0;
6kh
and
for
the
cubic
spline
approximation,
[sin(kh)
3
1
pc(kh)
kh
2+cos(kh)
>_
0.
It
is
clear
that
for
the
rth-order
difference
operator
p
(kh)
satisfies
I1-
p(kh)l
<
C(kh)
r.
Alternatively,
if
Dh
is
a
spectral
derivative,
applied
directly
in
the
Fourier
transform,
we
can
choose
p
(kh)
to
be
of
infinite
order.
In
this
case
we
will
assume
that
p
satisfies
the
following
properties:
(i)
p(-x)
p(x)
and
p(x)
>
0;
(ii)
p(.)
6
C
2
and
p(zr)
0;
and
(iii)
p(x)
for
Ix
_<
)zr,
where
0
<
2
<
1.
Property
(iii)
guarantees
that
Dh
is
spectrally
accurate.
We
denote
by
Sh
the
pseudospectral
derivative
operator
without
smoothing.
This
corresponds
to
the
case
of
p
(kh)
-=
in
the
definition
of
Dh.
The
filtering
or
smoothing
in
the
interface
variable
will
be
done
by
multiplying
by
p
(kh)
in
the
Fourier
transform.
When
zj
is
an
approximation
to
z
(cj),
sj
zj
cj
will
be
periodic.
We
define
z
as
otj
+
s,
where
^P
(kh)
(5)
s
p
It
is
clear
that
z
p
is
an
rth-order
approximation
to
z
if
p
corresponds
to
the
rth-order
derivative
operator.
Similarly,
if
Dh
is
a
spectral
derivative,
we
take
DhZj
to
mean
DhZj
1
+
Dhsj
1
+
Dh(Zj
--otj).
To
approximate
the
velocity
integral,
given
discrete
functions
zj
,
z(cj)
and
),j
,
(otj),
we
use
the
alternating
trapezoidal
rule,
with
smoothing
in
zj
as
above:
(6)
f_r
(Z(Oli)
Z(Olt)
)
,
(or’)
cot
dot’
.
2
?,j
cot
2h.
j=-N/2+I
2
(j-/)odd
Now
we
can
present
our
numerical
algorithm
for
the
water
wave
equations
(1)-(3)
as
follows:
(7)
dt
4rci
Y
yj
cot
2h
-Jr-
bti
ivi,
=-N/2+1
2
2(DhZ)i
(j-i)odd
A
CONVERGENT
NUMERICAL
METHOD
FOR
WATER
WAVES
1801
(8)
(9)
dqbi
-l
(u/2
+
v2i)
gYi
dt
2
9/j
cot
2h
(j-i)odd
These
equations
can
be
solved
once
initial
conditions
are
specified
for
zi
and
ti
and
a
time
discretization
is
chosen.
The
integral
equation
(9)
must
be
solved
for
9/i
at
each
time
step.
Its
solvability
is
proved
in
Lemma
5
below.
In
practice
it
is
solved
iteratively.
The
version
of
this
algorithm
with
surface
tension
will
be
discussed
later.
Error
estimates
will
be
given
in
terms
of
the
discrete
L2-norm,
given
by
N/2
(10)
Ilzll/2=
Izjlh.
j=-N/2+I
We
now
state
the
convergence
theorem
for
the
numerical
method
without
surface
tension,
followed
by
further
discussion.
THEOREM
1.
Assume
that
an
exact
solution
of
the
water
wave
equations
is
regular
enough
so
that
z(.,
t),
dp(.,
t)
cm+2[O,
2rr]
and
9/(.,
t)
cm+l[o,
27r]
for
m
>_
3,
and
Iz(ot,
t)
z(fl,
t)l
>_
clot
-/31
for
0
<
<
T
and
some
c
>
O.
Furthermore,
assume
the
condition
(11)
(ut,
vt)
n
(0,-g).
n
>_
co
>
0
holds
at
each
point
on
the
interface.
Here
(u,
v)
is
the
Lagrangian
velocity,
n
is
the
normal
vector
to
the
interface
(pointing
out
of
the
fluid
region),
and
co
is
some
constant.
Suppose
the
numerical
solution
z(t),
(t),
9/(t)
of
the
initial
value
problem
is
computed
using
algorithm
(7)-(9).
Then
if
Dh
is
an
rth-order
derivative
approximation
with
r
>
4,
we
have
for
h
<_
ho(T)
IIz(t)
z(.,
t)llt
_<
C(T)h
r,
114(t)
4(’,
t)llt
_<
C(T)h
r,
119/(t)
9/(’,
t)llt2
_<
C(T)h
r-1.
If
Dh
is
a
spectral
approximation
as
above,
we
have
the
same
convergence
result
with
h
replaced
by
h
m
in
the
right-hand
sides.
Condition
(11)
simply
means
that
the
interface
is
not
accelerating
downward,
normal
to
itself,
as
rapidly
as
the
normal
acceleration
of
gravity.
It
can
be
viewed
as
a
generalization
of
the
criterion
of
Taylor
[Tay]
for
horizontal
interfaces
to
rule
out
Rayleigh-Taylor
instabilities.
It
appears
naturally
as
a
sign
condition
in
the
argument
below,
as
well
as
in
the
analysis
of
[BHL1
].
Of
course
the
exact
solution
may
become
singular
at
a
later
time,
and
the
theorem
asserts
convergence
only
as
long
as
the
solution
is
regular.
Existence
results
for
time-dependent
water
waves
with
a
finite
degree
of
smoothness
are
rather
limited.
They
began
with
the
work
of
Nalimov;
see
[Craig]
and
the
references
cited
therein.
The
result
proved
here
could
be
extended
in
several
ways.
In
the
case
of
finite-order
derivative
operators,
we
can
improve
the
results
by
using
asymptotic
error
expansions
in
the
spirit
of
Strang
[Str].
Then
we
can
improve
the
convergence
rate
for
9/
to
the
optimal
order,
i.e.,
h
r.
Also,
Strang’s
argument
would
enable
us
to
prove
convergence
of
the
scheme
corresponding
to
the
second-order
centered
difference
approximation.
While
our
analysis
shows
that
one
set
of
choices
leads
to
a
fully
convergent
method,
it
is
of
course
possible
that