Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
SIAM J. S
CI
. C
OMPUT
.
©
2022 Society for Industrial and Applied Mathematics
Vol. 44, No. 2, pp. A964–A992
TWO-DIMENSIONAL FOURIER CONTINUATION AND
APPLICATIONS
∗
OSCAR P. BRUNO
†
AND
JAGABANDHU PAUL
†
Abstract.
This paper presents a fast “two-dimensional Fourier continuation” (2D-FC) method
for construction of biperiodic extensions of smooth nonperiodic functions defined over general two-
dimensional smooth domains. The approach, which runs at a cost of
O
(
N
log
N
) operations for an
N
-
point discretization grid, can be directly generalized to domains of any given dimensionality, but such
generalizations are not considered in this contribution. The 2D-FC extensions are produced in a two-
step procedure. In the first step the
one-dimensional
Fourier continuation method is applied along a
discrete set of outward boundary-normal directions to produce, along such directions, continuations
that vanish outside a narrow interval beyond the boundary. Thus, the first step of the algorithm
produces “blending-to-zero along normals” for the given function values. In the second step, the
extended function values are evaluated on an underlying Cartesian grid by means of an efficient,
high-order boundary-normal interpolation scheme. A Fourier continuation expansion of the given
function can then be obtained by a direct application of the two-dimensional fast Fourier transform
(FFT). Algorithms of arbitrarily high order of accuracy can be obtained by this method. The
usefulness and performance of the proposed 2D-FC method are illustrated with applications to the
Poisson equation and the time-domain wave equation within a bounded domain. As part of these
examples the novel “Fourier forwarding” solver is introduced which,
propagating plane waves as they
would in free space
and relying on certain boundary corrections, can solve the time-domain wave
equation and other hyperbolic partial differential equations
within general domains
at computing
costs that grow
sublinearly
with the size of the spatial discretization.
Key words.
two-dimensional Fourier continuation, Poisson equation, wave equation, FC solver,
Fourier forwarding, FFT
AMS subject classifications.
35J05, 78A45, 35L05, 42A10, 42A15
DOI.
10.1137/20M1373189
1. Introduction.
This paper presents a fast “two-dimensional Fourier Contin-
uation” (2D-FC) method for construction of biperiodic extensions of smooth non-
periodic functions defined over general two-dimensional (2D) smooth domains. The
algorithm, which runs at a cost of
O
(
N
log
N
) operations for an
N
-point discretization
grid, converges with a user-prescribed
d
th order of convergence, and can be directly
generalized to domains of any given dimensionality, but such generalizations are not
considered here. The usefulness and performance of the proposed 2D-FC method
are illustrated with applications to the Poisson equation and the time-domain wave
equation within a bounded domain. As part of these examples the novel “Fourier
forwarding” solver is introduced which,
propagating plane waves as they would in free
space
and relying on certain boundary corrections, can solve the time-domain wave
equation and other constant-coefficient hyperbolic partial differential equations
within
general domains
at computing costs that grow
sublinearly
with the size of the spatial
discretization.
∗
Submitted to the journal’s Methods and Algorithms for Scientific Computing section October 13,
2020; accepted for publication (in revised form) November 18, 2021; published electronically April
25, 2022.
https://doi.org/10.1137/20M1373189
Funding:
This work was supported by the NSF through grants DMS-1714169, DMS-2109831,
by DARPA through grant HR00111720035, by AFOSR through grant FA9550-21-1-0373, and by the
Vannewar Bush Foundation through grant N00014-16-1-2808.
†
Computing and Mathematical Sciences, Caltech, Pasadena, CA 91125 USA (obruno@
caltech.edu, jpaul@caltech.edu).
A964
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
TWO-DIMENSIONAL FOURIER CONTINUATION
A965
The periodic-extension problem has actively been considered in the recent lit-
erature, in view, in particular, of its applicability to the solution of various types
of partial differential equations (PDEs) [1, 3, 4, 6, 7, 8, 10, 11, 14, 17, 20, 21,
24]. The contributions [3, 4, 11, 20], in particular, utilize the Fourier continua-
tion (FC) method in one dimension in conjunction with dimensional splitting for
the treatment of multidimensional PDE problems. The dimensional splitting is also
used in [12] to produce Fourier extensions to rectangular domains in two dimen-
sions, where the FC is effected by separately applying the one-dimensional FC-Gram
method [3, 4, 11] first to the columns and then to the rows of a given data ma-
trix of function values. The method does assume that the given smooth function is
known on a rectangular region containing the domain for which the continuation is
sought. The 2D-FC algorithm proposed in this paper does not rely on such a stringent
requirement.
The approach to periodic function extension presented in [6, 24] is based on the
solution of a high-order PDE, where the extension shares the values and normal deriv-
atives along the domain boundary. [14], in turn, presents a function-extension method
based on the use of radial basis functions (RBFs). In that approach, overlapping cir-
cular partitions, or patches, are placed along the physical boundary of the domain,
and a local extension is defined on each patch by means of RBFs. A second layer of
patches is placed outside the first, on which the local values are set to vanish. The
zero patches are used in conjunction with a partition of unity function to blend the
local extensions into a global counterpart. The choice of functions used to build the
partition of unity determines the regularity of the extended function. [21] proposes
an accurate and superalgebraically convergent function approximation method using
Fourier extension frames in general 2D domains at a cost of
O
(
N
2
log
2
N
) operations:
a significantly higher cost than the
O
(
N
log
N
) cost required by the 2D-FC method
introduced in the present contribution. In fact, as indicated in that reference, the
O
(
N
2
log
2
N
) cost reported only improves upon the
O
(
N
3
) cost required by a cor-
responding full SVD-based algorithm for discretizations containing more than 8,100
points. As an example, cost and accuracy comparisons presented in Example 3.3
below, which concerns approximation of a function exhibiting 36 oscillations on a cir-
cular domain, and for both high (7
·
10
−
8
) and low (1
.
7
·
10
−
3
) accuracy, the 2D-FC
method is approximately 10,000 times faster and requires over 20 times less memory
than the approach in [21], with significant increases in improvement factors expected
as discretizations are enlarged further.
The 2D-FC extensions proposed in this paper are produced in a two-step proce-
dure. In the first step the
one-dimensional
FC (1D-FC) method [4] is applied along a
discrete set of outward boundary-normal directions to produce, along such directions,
continuations that vanish outside a narrow interval beyond the boundary. Thus, the
first step of the algorithm produces “blending-to-zero along normals” for the given
function values. In the second step, the extended function values are evaluated on
an underlying Cartesian grid by means of an efficient, high-order boundary-normal
interpolation scheme. Since the continuation-along-normals procedure is a fixed cost
operation for each normal line, the cost of this procedure grows only linearly with
the size of the boundary discretization. An FC expansion of the given function can
then be obtained by a direct application of the 2D FFT algorithm. Algorithms of
arbitrarily high order of accuracy can be obtained by this method. In view of its
construction on the basis of combined 1D-FC approximation and regular polynomial
interpolation, the convergence theory for the 2D-FC method follows directly from the
theoretical results presented in [20].
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
A966
OSCAR P. BRUNO AND JAGABANDHU PAUL
As mentioned above, this paper demonstrates the usefulness of the proposed
general-domain 2D-FC technique via applications to both the Poisson problem for
the Laplace equation and the time-domain wave equation. In the Poisson case the
2D-FC method is utilized to obtain a
particular solution
for a given right-hand side;
the boundary conditions are then made to match the prescribed boundary data by
adding a solution of the Laplace equation which is produced by means of boundary-
integral methods. The Fourier forwarding approach, in turn, uses the 2D-FC method
to solve the spatio-temporal PDE in the interior of the domain and it then corrects
the solution values near the boundary by means of a classical time-stepping solver.
The overall procedure, which utilizes large time steps for the interior solver and small
CFL-constrained time steps for the near-boundary solver, runs in computing times per
small time step that grow sublinearly with the size of the spatial discretization mesh.
It is interesting to note that the primary continuation device in the 2D-FC
method, namely, continuation along normals to the domain boundary, is a
one-
dimensional procedure
. This one-dimensional continuation procedure can be utilized
in a generalization of the method to
n
-dimensional domains with
n >
2. This is
in contrast to other extension methods mentioned above. For example, the RBF-
based extension method [14] requires solution of boundary problems of increasing
dimensionality as the spatial dimension grows, which, given the method’s reliance
on dense-matrix linear algebra for the local-extension process, could have a signifi-
cant impact on computing costs. Similar comments apply to PDE-based extension
methods such as [24].
The proposed 2D-FC algorithm performs favorably in the context of existing
related approaches. Specific comparisons with results presented in [14] are provided
in subsection 4.1.1 for a Poisson problem considered in that reference. The recent
contribution [13], in turn, presents an FFT-based high-order solver for the Poisson
problem for
rectangular
domains, namely, Cartesian products of one-dimensional (1D)
intervals in either 2D or three-dimensional space. The present 2D-FC based Poisson
solver achieves,
for general domains
, a similar performance (similar accuracy and
computing time) to that demonstrated in [13, Tables 3 and 4] under the Cartesian-
domain assumption.
As discussed in [3, 4], and particularly e.g. in [4, Figure 12], in view of the spec-
tral character of the 1D-FC approach, which incurs errors in finite-degree polynomial
approximation near boundaries, but which propagates the polynomial boundary ac-
curacy to the domain interior via Fourier series, enjoys excellent dispersion charac-
teristics as well: using FC derivatives as part of a PDE solver results in errors that,
asymptotically, do not grow as the number of oscillations in the PDE domain and the
discretization size are simultaneously and proportionally increased. This dispersion
characteristic is also enjoyed by the 2D-FC algorithm, as illustrated in the top graph
of Figure 11.
This paper is organized as follows. After a brief review of the 1D-FC method
presented in section 2, the proposed 2D-FC method is introduced in section 3. The two
main applications considered, namely, solution of the Poisson and Fourier-forwarding
for the wave equation, are presented in subsections 4.1 and 4.2. Finally our conclusions
are presented in section 5.
2. Background: 1D “blending-to-zero” FC algorithm.
A brief review of
the 1D-FC method [3, 4] is presented in this section, with an emphasis on one of
its key components, the
blending-to-zero
procedure—which is employed in the normal
direction continuation portion of the proposed 2D-FC approach presented in section 3.
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
TWO-DIMENSIONAL FOURIER CONTINUATION
A967
2.1. 1D-FC algorithm: Outline.
Given the vector
φ
D
= (
φ
0
,...,φ
N
−
1
)
t
of
values of a smooth function
φ
: [0
,
1]
→
C
on the equispaced grid
D
=
{
x
j
=
jk
: 0
≤
j
≤
N
−
1
}
of step size
k
= 1
/
(
N
−
1), the 1D-FC method [3, 4] of order
d
(with, e.g., 4
≤
d
≤
12)
produces, at first, an (
N
+
C
)-dimensional vector
φ
c
of discrete continuation function
values (including the
N
given function values) over an extended interval [0
,b
]
,b >
1.
To do this, the algorithm utilizes the
d
-dimensional vectors
φ
ℓ
= (
φ
0
,...,φ
d
−
1
)
t
and
φ
r
= (
φ
N
−
d
,...,φ
N
−
1
)
t
of values of the function
φ
on the left and right “matching-
point” sets
D
ℓ
=
{
x
0
,...,x
d
−
1
}
and
D
r
=
{
x
N
−
d
,...,x
N
−
1
}
, respectively, each one
of which is contained in a small subinterval of length (
d
−
1)
k
near the corresponding
endpoint of the containing interval [0
,
1]. In order to obtain the
C
necessary continua-
tion values, the 1D-FC method blends
φ
ℓ
and
φ
r
to zero (see subsection 2.2), towards
the left and right, respectively, resulting in two zero-blending vectors of length
C
.
The sum of these two vectors is then utilized as a rightward discrete continuation to
the set
D
c
=
{
x
j
= 1 +
jk
: 1
≤
j
≤
C
}
of points in the interval (1
,b
]—as described
in subsection 2.3. As indicated in that section, the overall 1D-FC procedure is then
completed via an application of the FFT algorithm to the (
N
+
C
)-dimensional vec-
tor
φ
c
(cf. (2.7) below) of “smoothly periodic” discrete continued function values.
The following two subsections describe the blending-to-zero and 1D-FC approaches,
respectively.
2.2. Blending-to-zero algorithm.
In our description of the order-
d
blending-
to-zero algorithm [4] we only present details for the
rightward
blending-to-zero tech-
nique, since the leftward blending-to-zero problem can easily be reduced to the right-
ward problem. Thus, given the column vector
F
D
= (
F
0
,...,F
d
−
1
)
t
of values of a
complex-valued function
F
on the set
D
=
{
x
0
,x
1
,...,x
d
−
1
}
, the rightward blending-
to-zero approach starts by producing a polynomial interpolant for
F
over the interval
[
x
0
,x
d
−
1
] relying on the Gram polynomial basis
G
d
=
{
g
0
(
x
)
,g
1
(
x
)
,...,g
d
−
1
(
x
)
}
(2.1)
for this interval. The functions
g
j
(
x
) (
j
= 0
,...,d
−
1) are the polynomials with real
coefficients that are obtained as the Gram–Schmidt orthogonalization procedure is
applied, in order of increasing degree, to the polynomials in the set
{
1
,x,x
2
,...,x
d
−
1
}
,
with respect to the discrete scalar product
(
g,h
) =
d
−
1
X
j
=0
g
(
x
j
)
h
(
x
j
)
.
Discrete values of the Gram polynomials on the set
D
can be computed on the basis
of the
QR
factorization [15]
P
=
QR
of the matrix
P
= (
x
i
j
−
1
)
0
≤
i,j
≤
d
−
1
.
(2.2)
(Note that the
j
th column of
Q
contains the values of the
j
th Gram polynomial on
the set
D
.) Following [4] we obtain the necessary
QR
factorization by applying the
stabilized Gram–Schmidt orthogonalization method to the matrix
P
.
In order to closely approximate each one of the Gram polynomials in
G
d
, through-
out the continuous interval [0
,
(
d
−
1)
k
] containing
D
, by corresponding trigonometric
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
A968
OSCAR P. BRUNO AND JAGABANDHU PAUL
polynomials, as described below, we use a certain “oversampled matching” method.
According to this method the polynomials in
G
d
are oversampled to an equispaced
set of discretization points with step size
k/n
os
(containing
n
os
(
d
−
1) + 1 points),
where
n
os
denotes the oversampling factor, and where the oversampled values are
used as part of a certain low-dimensional SVD matching procedure described in
what follows. Note that the aforementioned oversampled values on the refined grid
D
os
:
=
{
̃
x
j
=
jk/n
os
: 0
≤
j
≤
n
os
(
d
−
1)
}
coincide with the columns of the matrix
Q
os
=
P
os
R
−
1
,
(2.3)
where
P
os
is the Vandermonde matrix of size (
n
os
(
d
−
1) + 1)
×
d
corresponding to the
oversampled discretization
D
os
, and where
R
is the upper triangular matrix in (2.2).
The aforementioned SVD matching procedure, which is one of the crucial steps
in the FC approach [4], produces a band-limited Fourier series of the form
g
c
j
(
x
) =
J
X
m
=
−
J
a
j
m
e
2
πimx
(
d
+2
C
+
Z
−
1)
k
(2.4)
for each polynomial
g
j
∈
G
d
(0
≤
j
≤
d
−
1), where
C
is the number of blending-
to-zero values to be produced, and
Z
being the number of “zero-matching” points.
The Fourier coefficients are selected so as to match, in the least-squares sense, both
the oversampled polynomial values over the interval [0
,
(
d
−
1)
k
], and identically zero
values on an equally fine discretization of the “zero matching” interval
[(
d
+
C
)
k,
(
d
+
C
+
Z
−
1)
k
]
of length (
Z
−
1)
k
. The coefficients in (2.4) are taken to equal the solution
a
of the
minimization problem
min
a
=(
a
−
J
,...,a
J
)
T
B
os
a
−
q
j
os
0
2
,
(2.5)
where
B
os
is a matrix whose entries are values of (2.4) at all points in the set
D
os
as well as all the set of
k/n
os
-spaced points in the zero matching interval mentioned
above (which, in particular, contains the endpoints (
d
+
C
)
k
and (
d
+
C
+
Z
−
1)
k
). The
minimizing Fourier coefficients
a
are then found via an SVD-based [15] least-squares
approach. Once the coefficients
a
have been obtained, the resulting Fourier expan-
sions (2.4) are used to produce a certain “continuation matrix”
A
∈
C
C
×
d
, whose
columns equal the values of the expression (2.4) at the
C
(unrefined)
k
-discretization
points in the interval [
dk,
(
d
+
C
−
1)
k
] (cf. Remark 2.2 below). The desired vector
F
r
of rightward blending-to-zero function values at the
C
continuation points in the
interval [
dk,
(
d
+
C
−
1)
k
] is then given by the expression
F
r
=
AQ
T
F
.
(2.6)
2.3. 1D-FC algorithm.
As outlined in subsection 2.1, the 1D-FC algorithm re-
quires use of a certain rightward (resp., leftward) blending-to-zero vector
φ
r
r
(resp.,
φ
ℓ
ℓ
)
for a given matching-point vector
φ
r
(resp.,
φ
ℓ
). In view of (2.6) we define
φ
r
r
=
AQ
T
φ
r
. To obtain the leftward extension
φ
ℓ
ℓ
, in turn, we first introduce the “order
reversion” matrix
R
e
∈
C
e
×
e
(
e
∈
N
) by
R
e
(
g
0
,g
1
...,g
e
−
2
,g
e
−
1
)
t
= (
g
e
−
1
,g
e
−
2
,...,g
1
,g
0
)
t
,
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
TWO-DIMENSIONAL FOURIER CONTINUATION
A969
and we then define
φ
ℓ
ℓ
=
R
C
AQ
T
R
d
φ
ℓ
. A vector
φ
c
containing both the
N
given
values in the vector
φ
= (
φ
0
,φ
1
,...,φ
N
−
1
)
t
as well as the
C
“continuation” function
values is constructed by appending the sum
φ
ℓ
ℓ
+
φ
r
r
at the end of the vector
φ
, so
that we obtain
φ
c
j
=
φ
j
for 0
≤
j
≤
N
−
1
,
φ
ℓ
ℓ
+
φ
r
r
(
j
−
N
)
for
N
≤
j
≤
N
+
C
−
1
.
(2.7)
Following the various stages of the construction of the vector
φ
c
it is easy to
check that, up to numerical error, this vector contains point values of a smoothly
periodic function defined over the interval [0
,b
]. An application of the FFT algorithm
to this vector therefore provides the desired continuation function in the form of a
trigonometric polynomial,
φ
c
(
x
) =
(
N
+
C
)
/
2
X
ℓ
=
−
(
N
+
C
)
/
2
ˆ
φ
c
ℓ
e
2
πiℓx
b
,
(2.8)
which closely approximates
φ
in the interval [0
,
1]. In fact, as demonstrated in previous
publications (including [3, 4]), for sufficiently smooth functions
φ
, the corresponding
1D-FC approximants converge to
φ
with order
O
(
k
d
)—so that, as expected, the num-
ber
d
of points used in the blending-to-zero procedures determines the convergence
rate of the algorithm, while, as mentioned in section 1, giving rise to essentially
dispersionless numerics in the context of the PDE solution. The 2D-FC algorithm
introduced in the following section also relies on the 1D blending-to-zero procedure
described in subsection 2.2, and its convergence in that case is once again of the order
O
(
k
d
), and, when applied to solution of PDE problems, essentially dispersionless.
It is important to note that, for a given order
d
, the matrices
A
and
Q
can be
computed once and permanently stored on disc for use whenever an application of the
blending-to-zero algorithm is required—as these matrices do not depend on the point
spacing
k
. A graphical demonstration of various elements of the 1D-FC procedure is
presented in Figure 1.
Remark
2.1 (extra vanishing values). The 1D-FC implementations [3, 4] allow
for an additional number
E
≥
0 of identically zero “extra” function values to be
added on an (unrefined)
k
-discretization of the interval [(
d
+
C
)
k,
(
d
+
C
+
E
−
1)
k
],
as illustrated in Figure 1(c), to obtain a desired overall number of discrete function
values (including the given function values and the continuation values produced) such
as, e.g., a power of two or a product of powers of small prime numbers, for which the
subsequent application of the FFT is particularly well suited. The corresponding use
of extra vanishing values for the 2D continuation problem is mentioned in Remark 3.2.
Remark
2.2 (blending-to-zero on a refined grid). As indicated above in the present
section, the 2D-FC procedure introduced in section 3 utilizes the 1D blending-to-
zero strategy described above in this section to extend a function given on a 2D
domain Ω along the normal direction to Γ; the continuation values obtained at all
normals are then utilized to obtain the continuation function on the Cartesian grid
by interpolation. As detailed in subsection 3.4, in order to prevent accuracy loss
in the 2D interpolation step we have found it necessary to use 1D normal-direction
grids finer than the grids inherent in the blending-to-zero process itself. To easily
provide the necessary fine-grid values, a modified fine-grid continuation matrix
A
r
∈
C
C
r
×
d
is constructed, where
C
r
> C
denotes the number of fine-grid points utilized.
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
A970
OSCAR P. BRUNO AND JAGABANDHU PAUL
(a) Demonstration of the 1D-FC method. The continuation values are
computed as the sum of the blended-to-zero rightward and leftward ex-
tensions.
(b) Inset of Figure 1a.
(c) Numbers
E
of leftward and right-
ward extra zeroes.
Fig. 1
.
Illustration of the
1
D-FC procedure. Figure
1
(a) depicts the FC of the nonperiodic
function
φ
: [0
,
1]
→
R
given by
φ
(
x
) = exp(sin(5
.
4
πx
−
2
.
7
π
)
−
cos(2
πx
))
−
sin(2
.
5
πx
) + 1
. Fig-
ure
1
(b) presents a close up of the right continuation region
[1
−
(
d
−
1)
k,b
]
. Subsequently Figure
1
(c)
illustrates the use of a number
E
of extra zeroes in the blending-to-zero process, to yield a continu-
ation mesh containing FFT-friendly numbers (products of powers of small prime numbers) of point
values.
The modified continuation matrix
A
r
can be built on the basis of the minimizing
coefficients
a
in (2.5): the corresponding columns of the fine-grid continuation matrix
A
r
are obtained by evaluating (2.4) on the given fine-grid points in the interval ((
d
−
1)
k,
(
d
+
C
−
1)
k
]. The necessary blending-to-zero function values at the
C
r
fine-grid
points are given by
A
r
Q
T
φ
D
.
3. 2D-FC method.
This section presents the proposed volumetric FC method
on 2D domains Ω
⊂
R
2
with a smooth boundary Γ =
Ω
\
Ω, some elements of which
are illustrated in Figure 2. Let a smooth function
f
:
Ω
→
C
be given; we assume
that values of
f
are known on a certain uniform Cartesian grid within
Ω as well as a
grid of points on the boundary Γ. The 2D-FC algorithm first produces 1D blending-
to-zero values for the function
f
along directions normal to the boundary
Γ, yielding
continuation values on a certain 2D tangential-normal curvilinear grid around Γ, as
detailed in sections 3.1 and 3.2 and illustrated in Figure 2. These continuation values,
which are produced on the basis of the corresponding blending-to-zero procedure
presented in subsection 2.2 in the context of the 1D-FC method, are then interpolated
onto a Cartesian grid around the domain boundary, to produce a 2D blending-to-zero
continuation of the function
f
. The necessary interpolation from the curvilinear grid
to the Cartesian grid is accomplished by first efficiently obtaining the foot of the
normal that passes through a given Cartesian grid point
r
= (
x,y
) exterior to Ω and
near the boundary Γ (section 3.3), and then using a local 2D interpolation procedure
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
TWO-DIMENSIONAL FOURIER CONTINUATION
A971
Fig. 2
.
Geometrical constructions underlying the
2
D-FC procedure, with reference to the various
regions defined in subsection
3.1
.
to produce the corresponding continuation value at the point
r
(section 3.4). Once
the interpolated values have been obtained throughout the Cartesian mesh around Γ,
the desired 2D-FC function
f
c
(
x,y
) =
N
x
/
2
X
ℓ
=
−
N
x
/
2+1
N
y
/
2
X
m
=
−
N
y
/
2+1
ˆ
f
c
ℓ,m
e
2
πi
ℓx
L
x
+
my
L
y
(3.1)
(where
L
x
and
L
y
denote the period in the
x
and
y
directions, respectively) is ob-
tained by means of a 2D FFT. Following the algorithmic prescriptions presented in
sections 3.1 through 3.4, a summary of the overall 2D-FC approach is presented in
section 3.5.
3.1. 2D tangential-normal curvilinear grid.
The necessary curvilinear grids
around Γ can be produced on the basis of a (smooth) parametrization
r
=
q
(
θ
) = (
x
(
θ
)
,y
(
θ
))
,
0
≤
θ
≤
2
π,
(3.2)
of the boundary Γ. We assume the boundary to be of class
C
d
+1
to ensure the inter-
polation processes described in subsection 3.4 enjoy the necessary
O
(
h
d
) convergence
order. However, per Example 3.2, use of interpolation of order
M
=
d
+ 2 in that
context yields somewhat improved accuracy. In order to achieve such improvements
it is necessary to assume the boundary Γ is a
C
d
+3
curve.
In view of their intended application (blending-to-zero along the normal direction
in accordance with subsection 2.2), the curvilinear grids are introduced within interior
and exterior strips
V
−
and
V
+
(illustrated in Figure 2) given by
(
V
−
=
{
q
(
θ
)
−
n
(
θ
)
γ
: 0
≤
θ
≤
2
π
and 0
≤
γ
≤
(
d
−
1)
k
1
}
,
V
+
=
{
q
(
θ
) +
n
(
θ
)
γ
: 0
≤
θ
≤
2
π
and 0
≤
γ
≤
Ck
1
}
,
(3.3)
where
n
(
θ
) = (
n
x
(
θ
)
,n
y
(
θ
)) denotes the unit normal to the boundary Γ at
q
(
θ
), and
where
d
,
C
, and
k
1
denote, respectively, the number of matching points, the number
of blending-to-zero points, and the step size used, in the present application of the
1D blending-to-zero procedure described in subsection 2.2. Using, in addition, the
uniform discretization
I
B
=
{
θ
p
=
pk
2
: 0
≤
p < B
}
, k
2
=
2
π
B
,
(3.4)
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
A972
OSCAR P. BRUNO AND JAGABANDHU PAUL
of the interval [0
,
2
π
], we then construct a curvilinear 2D discretization
V
−
B,d
=
{
r
p,q
:
r
p,q
=
q
(
θ
p
) +
n
(
θ
p
)(
q
−
d
+ 1)
k
1
; 0
≤
p < B
and 0
≤
q
≤
d
−
1
}
,
(3.5)
V
+
B,C
r
=
{
s
p,q
:
s
p,q
=
q
(
θ
p
) +
n
(
θ
p
)
qk
1
/n
r
; 0
≤
p < B
and 0
≤
q
≤
C
r
}
,
(3.6)
within
V
−
and
V
+
respectively, for the given step size
k
1
, where
C
r
=
Cn
r
for certain
integer (refinement factor)
n
r
; note that the points in
V
−
B,d
for
q
=
d
−
1 and the points
in
V
+
B,C
r
for
q
= 0 coincide and that they lie on Γ. Here the constants
d
,
C
, and
n
r
are independent of
B
. The continuation function is constructed so as to vanish at all
points
s
p,q
∈
V
+
B,C
r
with
q
=
C
r
. Let now
R
= [
a
0
,a
1
]
×
[
b
0
,b
1
] denote the smallest
closed rectangle containing Ω
∪
V
+
, and consider the equispaced Cartesian grid of
step size
h
,
H
=
{
z
i,j
= (
x
i
,y
j
) :
x
i
=
a
0
+
ih
;
y
j
=
b
0
+
jh
: 0
≤
i < N
x
,
0
≤
j < N
y
}
(3.7)
on
R
, where the 2D continuation function values are to be computed. We note that
the size of the rectangle
R
along with the strips
V
−
and
V
+
decrease as the step size
k
2
is decreased.
3.2. Computation of FC values on
V
+
B,C
r
.
A continuation of the function
f
to the exterior of Ω is obtained via application of the blending-to-zero procedure
presented in subsection 2.2 (cf. Remark 2.2) along each one of the normal directions
inherent in the definition of the set
V
+
B,C
r
. For given
p
, the
d
equidistant points
s
p,q
∈
V
−
B,d
(0
≤
q
≤
d
−
1), which are indicated by the solid circles in Figure 3,
constitute a set
D
p
of matching points that are used to effect the blending-to-zero
procedure per the prescriptions presented in subsection 2.2. To obtain the desired
continuation function values it is necessary to first obtain the vector
f
D
p
of the values
of the function
f
(or suitable approximations thereof) on the set
D
p
. In the proposed
method, the needed function values
f
D
p
are computed on the basis of a two-step
polynomial interpolation scheme, using polynomials of a certain degree (
M
−
1), as
briefly described in what follows.
With reference to the right image of Figure 3, and considering first the case
|
n
x
(
θ
p
)
| ≥ |
n
y
(
θ
p
)
|
, the algorithm initially interpolates vertically the function values
at
M
open-circle Cartesian points selected as indicated in the figure, onto the points
of intersection, shown as red stars, of the normal line and the vertical Cartesian
grid lines. For intersection (red-star) points close enough to the boundary, boundary
function values at boundary points shown as squares in the figure, are utilized in the
1D interpolation process as well. Once the red-star function values are acquired, the
function value at the matching solid-black point is effected by interpolation from the
M
red-star point values previously obtained, on the basis of a polynomial of degree
(
M
−
1).
The case
|
n
x
(
θ
p
)
|
<
|
n
y
(
θ
p
)
|
is treated similarly, substituting the initial interpola-
tion along vertical Cartesian lines, by interpolation along horizontal Cartesian lines;
the algorithm then proceeds in an entirely analogous fashion.
3.3. Proximity map.
As described below in section 3.4, the 2D-FC algorithm
interpolates the FC values on
V
+
B,C
r
onto the Cartesian mesh
H
. The interpolation al-
gorithm used in that section relies on a certain “proximity map”
P
:
H
∩
V
+
→
V
+
B,C
r
which associates a curvilinear grid point
s
p,q
=
P
(
z
i,j
)
∈
V
+
B,C
r
in the “proximity”
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
TWO-DIMENSIONAL FOURIER CONTINUATION
A973
Fig. 3
.
Interpolation scheme for evaluation of
f
D
p
in the case
|
n
x
(
θ
p
)
| ≥ |
n
y
(
θ
p
)
|
. Left:
Black solid circles indicate the matching points that define the set
D
p
. Right: Known function
values at the Cartesian points (denoted by the empty circles) and, in some cases, at the point of
intersection (represented by an empty square) of the vertical grid lines with
Γ
, are used to interpolate
the matching function values at the red-star intersection points of the normal with the vertical grid
lines. The function values at the red points are then used to obtain, by interpolation, the function
values at the matching points.
of each given Cartesian grid point
z
i,j
. The algorithm described in the next sec-
tion primarily uses the associated point
s
p,q
along with the corresponding discretized
boundary parameter value
θ
p
to interpolate the foot of the normal for the Cartesian
point
z
i,j
. The proximity function we use is obtained by first associating with each
curvilinear discretization point
s
p,q
the nearest Cartesian point, a procedure that re-
sults in a set
P
0
⊆
(
H
∩
V
+
)
×
V
+
B,C
r
of pairs of points, one in the Cartesian grid and
the other in the curvilinear grid. (The initial set
P
0
can easily be obtained by using
the “integer part”
floor
(
⌊
.
⌋
) and the
ceil
(
⌈
.
⌉
) operators.) The set
P
0
is then modi-
fied by removing multiple associations for a given Cartesian point, and, if necessary,
by adding a “next-nearest” curvilinear neighbor to Cartesian points that previously
remained unassociated. The next-nearest curvilinear neighbor can be obtained simply
by subtracting one from the floor value and adding one to the ceil value and, taking
all the combinations of the results for each point in
V
+
B,C
r
. The resulting set
P
defines
the desired function.
3.4. FC values on the Cartesian grid.
Once FC values on
V
+
B,C
r
have been
obtained, per the procedure presented in subsection 3.2, the 2D-FC scheme can be
completed by (a) interpolation onto the set
H
∩
V
+
of outer Cartesian grid points;
and (b) subsequent evaluation of the corresponding FCs in (3.1) by means of an FFT.
(Note that since the continuation function
f
c
is a smooth function which vanishes
outside
Ω
∪
V
+
, this function can be viewed as the restriction of a smooth and biperi-
odic function with periodicity rectangle
R
—whose Fourier series approximates
f
c
and, therefore
f
, accurately.) The efficiency of the interpolation scheme is of the ut-
most importance in this context—since interpolation to a relatively large set
H
∩
V
+
of Cartesian points is necessary. An accurate and efficient interpolation strategy is
obtained by combining two 1D (local) interpolation procedures based on nearby nor-
mal directions. The first interpolation procedure produces the parameter value
θ
of the foot of the normal to Γ passing through a given Cartesian point; the second
procedure then approximates the continuation function value utilizing the parameter
value
θ
just mentioned and the continuation function values at the points in
V
+
B,C
r
around the given Cartesian point. A detailed description of the combined interpola-
tion methodology is presented in what follows. Specifically, we describe the strategy
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
A974
OSCAR P. BRUNO AND JAGABANDHU PAUL
we use to interpolate the continuation function onto each point
Q
=
z
i,j
∈
H
∩
V
+
and, to do this, we first obtain the foot
F
(
Q
) of this point. Using the proximity map
P
described in subsection 3.3, the algorithm utilizes the curvilinear discretization point
s
p,q
=
P
(
z
i,j
)
∈
V
+
B,C
r
as well as the corresponding boundary discretization parame-
ter value
θ
p
; according to (3.6), the point
q
(
θ
p
) equals the foot of the normal passing
through
s
p,q
:
q
(
θ
p
) =
F
(
s
p,q
). The algorithm then seeks approximation of the foot
F
(
Q
) and the corresponding parameter value
θ
=
θ
Q
via a preliminary interpolation
step, as indicated in what follows. The foot
F
(
Q
) and the corresponding parameter
value
θ
=
θ
Q
are then used by the algorithm to produce the desired interpolated
continuation value at
Q
.
In order to obtain
F
(
Q
) the algorithm uses the
M
boundary parameter values in
the set
S
θ
p
=
{
θ
p
−
K
ℓ
,θ
p
−
K
ℓ
+1
,...,θ
p
,θ
p
+1
,...,θ
p
+
K
r
}⊂
I
B
(3.8)
around
θ
p
(where
K
r
+
K
ℓ
+ 1 =
M
, and where
K
r
=
K
ℓ
if
M
is odd and
K
r
=
K
ℓ
+ 1
if
M
is even. Parameter values
θ
k
with negative values of
k
, which may arise in (3.8),
are interpreted by periodicity:
θ
k
=
θ
B
+
k
).
The algorithm then utilizes the line
L
⊥
Q
passing through
Q
that is orthogonal to
the normal vector
n
(
θ
p
) (see left image in Figure 4), together with the parametrization
ℓ
tr
Q
(
τ
) of
L
⊥
Q
, where the parameter
τ
represents the signed distance of the points on
L
⊥
Q
from the point
Q
. Clearly, then,
ℓ
tr
Q
(0) =
Q
. Each point of intersection of
L
⊥
Q
with
the normals
n
(
θ
j
) (
θ
j
∈
S
θ
p
), on the other hand, equals
ℓ
tr
Q
(
τ
j
), where
τ
j
denotes the
signed distance between
Q
and the corresponding intersection point. Thus, defining
the function
θ
=
T
(
τ
), where
T
(
τ
) gives the parameter value of the foot of the normal
through the point
ℓ
tr
Q
(
τ
), we clearly have
θ
j
=
T
(
τ
j
)
, p
−
K
ℓ
≤
j
≤
p
+
K
r
.
(3.9)
It follows that a 1D interpolation procedure on the function
T
(
τ
) can be used to obtain
the desired approximation of the value
θ
Q
=
T
(0) of the parameter corresponding to
the foot of the point
Q
=
z
i,j
:
F
(
Q
) =
q
(
θ
Q
) =
q
(
T
(0)).
Fig. 4
.
Interpolation schemes utilized to obtain the continuation function values on
H
∩
V
+
on
the basis of the continuation values on
V
+
B,C
r
. Left: Evaluation of the boundary parameter value for
the foot of the normal line passing through
Q
, depicted as a finely dotted line passing through that
point. The left image also displays the set of red-star interpolation points along the dashed-orange
line
L
⊥
Q
. Right: Interpolation of continuation values from the curvilinear mesh to a point
Q
on the
Cartesian mesh.
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy
Copyright
©
by SIAM. Unauthorized reproduction
of this article is prohibited.
TWO-DIMENSIONAL FOURIER CONTINUATION
A975
Once we have the corresponding foot parameter value
θ
Q
for the given Cartesian
point
Q
, the distance
η
Q
of the point
Q
to the boundary Γ is easily computed. Let
S
θ
Q
⊂
I
B
be a set of
M
boundary parameter values, similar to the set
S
θ
p
defined
in (3.8), but with
S
θ
Q
“recentered” around
θ
Q
. In order to obtain the continuation
function value at the point
Q
using the continuation function values on
V
+
B,C
r
, we
employ a local 2D polynomial interpolation scheme based on the set
S
θ
Q
and the
distance
η
Q
, as indicated in what follows and illustrated in the right image of Figure 4.
First, the continuation values are obtained at each point (marked by blue asterisks in
the figure) at a distance
η
Q
from the boundary Γ along the normal grid lines in
V
+
B,C
r
that correspond to boundary parameter values in the set
S
θ
Q
. Each one of these values
is obtained via 1D interpolation of the continuation function values on
V
+
B,C
r
along
the corresponding normal grid line in
V
+
B,C
r
. The desired continuation value at the
point
Q
, then, is obtained via a final 1D degree-(
M
−
1) interpolation step based on
the parameter set
S
θ
Q
and values at the “blue-asterisk” points just obtained. Finally,
by applying the 2D FFT to the continuation function values computed above we then
obtain the desired Fourier series expression in (3.1) for the continuation function.
Remark
3.1 (function values on Γ). I. For definiteness, in this paper we have
assumed that the boundary data are provided in the form of values of the given
function—which correspond to the the Dirichlet boundary data in the PDE context.
But the approach is also applicable in cases for which the boundary data are given as
the normal derivative of the function, (Neumann boundary data), or even a combi-
nation of function and normal derivative values (Robin data) by relying on a slightly
modified blending-to-zero procedure of the type presented in [4]. II. If no boundary
data are available, the 2D-FC method can still be utilized on the basis of interior data
only, albeit with a certain reduction in accuracy near the boundary.
Remark
3.2 (extra vanishing values in 2 dimensions). As in the 1D case, prior
to the FFT procedure the grid
H
can be enlarged, with vanishing function values
assigned to the added discretization points to obtain a discretization containing a
number of discretization points equal to a power of two (or a product of powers of
small prime numbers) along each Cartesian direction, which leads to especially fast
evaluations by means of the FFT.
3.5. Summary of the 2D-FC procedure.
This section presents a summary
of the 2D-FC procedure described in sections 3.1 through 3.4 for a function
f
given
on a uniform Cartesian grid
H
∩
Ω within the domain of definition Ω, where
H
is
a Cartesian mesh over the rectangle
R
containing both Ω and the near-boundary
outer region
V
+
; see subsection 3.1 and, in particular, Figure 2. The construction
of the continuation function
f
c
for the given function
f
relies on the use of three
main parameters associated with the 1D blending-to-zero approach presented in sec-
tion 2, namely,
d
(number of points in the set
D
),
C
(number of unrefined discrete
continuation points), and
n
os
(oversampling factor for the 1D blending-to-zero FC
procedure), together with the parameters
n
r
(refinement factor for the discrete con-
tinuation points; Remark 2.2 and subsection 3.1), and
M
−
1 (degree of the inter-
polating polynomials; subsections 3.2 and 3.4). Additionally, the 2D-FC procedure
utilizes the precomputed matrices
A
r
and
Q
, which, with reference to subsection 2.2,
are obtained as per the description provided in Remark 2.2.
Using the aforementioned parameters and matrices, the algorithm proceeds in two
main steps, namely, step (a) A
geometrical setup
precomputation procedure (com-
prising points 1 through 4 below); and step (b) A
Cartesian evaluation
procedure
Downloaded 06/02/22 to 131.215.70.177 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy