2010.03901.pdf

Two-dimensional Fourier Continuation and applications

Oscar P. Bruno

†

and Jagabandhu Paul

†

Computing and Mathematical Sciences, Caltech, Pasadena, CA 91125

Abstract

This paper presents a “two-dimensional Fourier Continuation” method (2D-FC) for construction of

bi-periodic extensions of smooth non-periodic functions defined over general two-dimensional smooth

domains. The approach can be directly generalized to domains of any given dimensionality, and even

to non-smooth domains, but such generalizations are not considered here. 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 func-

tion 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 two-dimensional Fourier Continuation method are illustrated with appli-

cations 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.

Keywords

: Two-dimensional Fourier Continuation, Poisson Equation, Wave Equation, FC Solver, Fourier

Forwarding, FFT.

1 Introduction

This paper presents a “two-dimensional Fourier Continuation” method (2D-FC) for construction of bi-

periodic extensions of smooth non-periodic functions defined over general two-dimensional smooth domains.

The approach can be directly generalized to domains of any given dimensionality, and even to non-smooth

domains, but such generalizations are not considered here. The usefulness and performance of the proposed

two-dimensional Fourier Continuation 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.

The periodic-extension problem has actively been considered in the recent literature, in view, in par-

ticular, of its applicability to the solution of various types of Partial Differential Equations (PDE) [1, 3,

4, 6, 8, 13, 16, 20, 23]. The contributions [3, 4, 8, 20], in particular, utilize the Fourier Continuation (FC)

method in one dimension in conjunction with dimensional splitting for the treatment of multidimensional

arXiv:2010.03901v2 [math.NA] 15 Oct 2020

PDE problems. The dimensional splitting is also used in [11] to produce Fourier extensions to rectangular

domains in two dimensions, where the Fourier Continuation is effected by separately applying the one-

dimensional FC-Gram method [3, 4, 8] first to the columns and then to the rows of a given data matrix

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 approach to periodic function extension presented in [6, 23] is based on the solution of a high-

order PDE, where the extension shares the values and normal derivatives along the domain boundary.

Reference [13], in turn, presents a function-extension method based on use of radial basis functions. In

that approach, overlapping circular partitions, or patches, are placed along the physical boundary of the

domain, and a local extension is defined on each patch by means of Radial Basis Functions (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-up the partition of unity determines the regularity of

the extended function.

The 2D-FC extensions proposed in this paper are produced in a two-step procedure. In the first step

the

one-dimensional

Fourier Continuation 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. A Fourier

Continuation expansion of the given function can then be obtained by a direct application of the two-

dimensional FFT algorithm. Algorithms of arbitrarily high order of accuracy can be obtained by this

method. Since the continuation-along-normals procedure is a fixed cost operation, the cost of the method

grows only linearly with the size of the boundary discretization.

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 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

-dimensional domains with

n >

. This

is in contrast with other extension methods mentioned above. For example, the RBF-based extension

method [13] 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 significant impact on computing costs. Similar comments apply to PDE-based extension

methods such as [23].

The proposed 2D-FC algorithm performs favorably in the context of existing related approaches. Spe-

cific comparisons with results presented in [13] are provided in Section 4.1.1 for a Poisson problem con-

sidered in that reference. The recent contribution [12], in turn, presents an FFT-based high-order solver

for the Poisson problem for

rectangular

domains, namely, Cartesian products of one-dimensional intervals

in either two- 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 [12, Tables

3 and 4] under the Cartesian-domain assumption.

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 Sections 4.1 and 4.2.

Finally our conclusions are presented in Section 5.

(a) Demonstration of the 1D-FC method. The continuation values are computed as the sum of the

blended-to-zero rightward and leftward extensions.

(b) Inset of Figure 1a.

of leftward and rightward extra ze-

roes.

Figure 1:

Illustration of the 1D-FC procedure. Figure 1a depicts the Fourier Continuation of the non-periodic function

: [0

→

given by

(

) = exp(sin(5

πx

−

)

−

cos(2

πx

))

−

sin(2

πx

) + 1

. Figure 1b presents a close up of the right

continuation region

−

(

−

k,b

]

. Subsequently Figure 1c illustrates the use of a number

of extra zeroes in the blending

to zero process, to yield a continuation mesh containing FFT-friendly numbers (products of powers small prime numbers) of

point values.

2 Background: 1D “blending-to-zero” FC algorithm

This section presents a brief review of the one-dimensional Fourier Continuation (1D-FC) method [3, 4],

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 two-dimensional Fourier Continuation (2D-FC)

approach presented in Section 3.

2.1 1D-FC algorithm: Outline

Given the vector

= (

,...,φ

−

)

of values of a smooth function

: [0

→

on the equispaced

grid

{

: 0

≤

−

}

of step-size

= 1

(

−

, the 1D-FC method [3,4] of order

(with,

e.g.,

≤

) produces, at first, an (

)-dimensional vector

of discrete continuation function

values (including the

given function values) over an extended interval

]

,b >

. To do this, the

algorithm utilizes the

-dimensional vectors

= (

,...,φ

−

)

and

= (

−

,...,φ

−

)

of values of

the function

on the left and right “matching-point” sets

{

,...,x

−

}

and

{

−

,...,x

−

}

respectively, each one of which is contained in a small subinterval of length

(

−

near the corresponding

endpoint of the containing interval

. In order to obtain the

necessary continuation values, the 1D-

FC method blends

and

to zero (see Section 2.2), towards the left and right, respectively, resulting

in two zero-blending vectors of length

. The sum of these two vectors is then utilized as a rightward

discrete continuation to the set

{

= 1 +

: 1

≤

}}

of points in the interval

]

—as

described in Section 2.3. As indicated in that section, the overall 1D-FC procedure is then completed

via an application of the FFT algorithm to the (

)-dimensional vector

(cf. equation (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-

blending-to-zero algorithm [4] we only present details for the

rightward

blending-to-zero technique, since the leftward blending-to-zero problem can easily be reduced to the right-

ward problem. Thus, given the column vector

= (

,...,F

−

)

of values of a complex-valued function

on the set

{

,...,x

−

}

, the rightward blending-to-zero approach starts by producing a poly-

nomial interpolant for

over the interval

[

−

]

relying on the Gram polynomial basis

{

(

)

(

)

,...,g

−

(

)

}

(2.1)

for this interval. The functions

(

)

(

= 0

,...,d

−

) 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

{

,x,x

,...,x

−

}

, with respect to the discrete scalar product

(

g,h

) =

−

∑

(

)

(

)

Discrete values of the Gram polynomials on the set

can be computed on the basis of the

factoriza-

tion [14]

of the matrix

= (

−

)

≤

i,j

≤

−

(2.2)

(Note that the

column of

contains the values of the

Gram polynomial on the set

.) Following [4]

we obtain the necessary

factorization by applying the stabilized Gram-Schmidt orthogonalization

method to the matrix

In order to closely approximate each one of the Gram polynomials in

, throughout the continuous

interval

(

−

]

containing

, by corresponding trigonometric polynomials, as described below, we use a

certain “oversampled matching” method. According to this method the polynomials in

are oversampled

to an equispaced set of discretization points with stepsize

k/n

(containing

(

−

1) + 1

points) where

denotes the oversampling factor, and where the oversampled values are used as part of a certain Singular

Value Decomposition (SVD) matching procedure described in what follows. Note that the aforementioned

oversampled values on the refined grid

{

jk/n

: 0

≤

(

−

}

coincide with the columns

of the matrix

−

(2.3)

where

is the Vandermonde matrix of size

(

−

1) + 1)

corresponding to the oversampled

discretization

, and where

is the upper triangular matrix in equation (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

(

) =

∑

−

πimx

(

−

(2.4)

for each polynomial

∈

(

≤

−

), where

is the number of blending-to-zero values to be

produced, and

being the number of “zero-matching” points. The Fourier coefficients are selected so as to

match, in the least square sense, both the oversampled polynomial values over the interval

(

−

]

, and

identically zero values on an equally fine discretization of the “zero matching” interval

[(

)

(

−

]

of length

(

−

. The coefficients in (2.4) are taken to equal the solution

of the minimization

problem

min

−

,...,a

)

∥

−

(

0

)

∥

(2.5)

where

is a matrix whose entries are values of (2.4) at all points in the set

as well as all the set of

k/n

-spaced points in the “zero matching” interval mentioned above (which, in particular, contains the

endpoints

(

)

and

(

−

)

. The minimizing Fourier coefficients

are then found via an

SVD-based [14] least-squares approach. Once the coefficients

have been obtained, the resulting Fourier

expansions (2.4) are used to produce a certain “continuation matrix”

∈

, whose columns equal the

values of the expression (2.4) at the

(unrefined)

-discretization points in the interval

[

dk,

(

−

]

(cf. Remark 2.2 below). The desired vector

of rightward blending-to-zero function values at the

continuation points in the interval

[

dk,

(

−

]

is then given by the expression

(2.6)

2.3 1D-FC Algorithm

As outlined in section 2.1, the 1D-FC algorithm requires use of a certain rightward (resp. leftward) blending-

to-zero vector

(resp.

) for a given matching-point vector

(resp.

). In view of equation (2.6) we

define

. To obtain the leftward extension

, in turn, we first introduce the “order reversion”

matrix

∈

(

∈

) by

(

...,g

−

)

= (

−

,...,g

)

and we then define

. A vector

containing both the

given values in the vector

= (

,φ

,...,φ

−

)

as well as the

“continuation” function values is constructed by appending the

sum

at the end of the vector

, so that we obtain







for

≤

−

(

)

(

−

)

for

≤

−

(2.7)

Following the various stages of the construction of the vector

it is easy to check that, up to numerical

error, this vector contains point values of a smoothly periodic function defined over the interval

]

. An

application of the FFT algorithm to this vector therefore provides the desired continuation function in the

form of a trigonometric polynomial,

(

) =

(

)

∑

−

(

)

πi`x

(2.8)

which closely approximates

in the interval

. In fact, as demonstrated in previous publications

(including [3,4]), for sufficiently smooth functions

, the corresponding 1D-FC approximants converge to

with order

(

)

—so that, as expected, the number

of points used in the blending-to-zero procedures

determines the convergence rate of the algorithm. (As discussed in these publications, further, in view of

its spectral character the 1D-FC approach enjoys excellent dispersion characteristics as well.) The 2D-FC

algorithm introduced in the following section also relies on the one-dimensional blending-to-zero procedure

described in section 2.2, and its convergence in that case is once again of the order

(

)

It is important to note that, for a given order

, the matrices

and

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

. 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

≥

of identically zero “Extra” function values to be added on a (unrefined)

-discretization of

the interval

[(

)

(

−

]

, as illustrated in Figure 1c, 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 fast Fourier transform 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

two-dimensional Fourier continuation procedure introduced in Section 3 utilizes the 1D blending-to-zero

strategy described above in this section to extend a function given on a two-dimensional 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 Section 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

∈

is constructed, where

> C

denotes the

number of fine-grid points utilized. The modified continuation matrix

can be built on the basis of the

minimizing coefficients

in (2.5): the corresponding columns of the fine-grid continuation matrix

are

obtained by evaluating (2.4) on the given fine-grid points in the interval

((

−

(

−

]

. The

necessary blending-to-zero function values at the

fine-grid points are given by

3 Two-dimensional Fourier Continuation Method

This section presents the proposed volumetric Fourier continuation method on two-dimensional domains

Ω

⊂

with a smooth boundary

Γ =

Ω

, some elements of which are illustrated in Figure 2. Let a smooth

function

Ω

→

be given; we assume that values of

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 one-dimensional

“blending-to-zero” values for the function

along directions normal to the boundary

, yielding continuation

values on a certain two-dimensional 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 section 2.2 in the context of the 1D-FC method,

are then interpolated onto a Cartesian grid around the domain boundary, to produce a two-dimensional

blending-to-zero continuation of the function

. 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

= (

x,y

)

exterior to

Ω

and near the boundary

(Section 3.3), and then

using a local two-dimensional interpolation procedure to produce the corresponding continuation value at

the point

(Section 3.4). Once the interpolated values have been obtained throughout the Cartesian mesh

around

, the desired two-dimensional Fourier continuation function

(

x,y

) =

∑

−

2+1

∑

−

2+1

`,m

πi

(

)

(3.1)

(where

and

denotes the period in the

and

directions, respectively) is obtained 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 Two-dimensional tangential-normal curvilinear grid

The necessary curvilinear grids around

can be produced on the basis of a (smooth) parametrization

(

) = (

(

)

(

))

≤

π,

(3.2)

of the boundary

. In view of their intended application (blending to zero along the normal direction in

accordance with section 2.2), the curvilinear grids are introduced within interior and exterior strips

−

and

(illustrated in Figure 2) given by

{

−

{

(

)

−

(

)

: 0

≤

and

≤

(

−

}

{

(

) +

(

)

: 0

≤

and

≤

}

(3.3)

where

(

) = (

(

)

(

))

denotes the unit normal to the boundary

(

)

, and where

and

denote, respectively, the number of matching points, the number of blending-to-zero points, and the

stepsize used, in the present application of the 1D blending-to-zero procedure described in Section 2.2.

Using, in addition, the uniform discretization

{

: 0

≤

p < B

}

, k

(3.4)

of the interval

]

, we then construct a curvilinear two-dimensional discretization

−

B,d

{

p,q

(

) +

(

)(

−

+ 1)

; 0

≤

p < B

and

≤

−

}

(3.5)

B,C

{

p,q

(

) +

(

)

; 0

≤

p < B

and

≤

}

(3.6)

within

−

and

respectively, for the given stepsize

where

for certain integer (refinement

factor)

; note that the points in

−

B,d

for

−

and the points in

B,C

for

= 0

coincide and

that they lie on

. Here the constants

and

are independent of

. The continuation function is

constructed so as to vanish at all points

p,q

∈

B,C

with

. Let now

= [

]

[

]

denote

the smallest closed rectangle containing

Ω

∪

, and consider the equispaced Cartesian grid of stepsize

{

i,j

= (

) :

;

: 0

≤

i < N

≤

j < N

}

(3.7)

, where the two-dimensional continuation function values are to be computed. We note that the size

of the rectangle

along with the strips

−

and

decrease as the stepsize

is decreased.

3.2 Computation of FC values on

V

B,C

A continuation of the function

to the exterior of

Ω

is obtained via application of the blending-to-zero

procedure presented in Section 2.2 (cf. Remark 2.2) along each one of the normal directions inherent in

the definition of the set

B,C

. For given

, the

equidistant points

p,q

∈

−

B,d

≤

−

, which

Figure 2:

Geometrical constructions underlying the 2D-FC procedure, with reference to the various regions defined in Sec-

tion 3.1.

are indicated by the solid circles in Figure 3, constitute a set

of matching points that are used to

effect the blending-to-zero procedure per the prescriptions presented in Section 2.2. To obtain the desired

continuation function values it is necessary to first obtain the vector

of the values of the function

(or

suitable approximations thereof) on the set

. In the proposed method, the needed function values

are computed on the basis of a two-step polynomial interpolation scheme, using polynomials of a certain

degree (

−

), as briefly described in what follows.

With reference to the right image of Figure 3, and considering first the case

(

)

| ≥ |

(

)

, the

algorithm initially interpolates vertically the function values at

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 one-dimensional 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

red-star point values previously obtained, on the basis

of a polynomial of degree (

−

The case

(

)

(

)

is treated similarly, substituting the initial interpolation along vertical

Cartesian lines, by interpolation along horizontal Cartesian lines; the algorithm then proceeds in an entirely

analogous fashion.

Figure 3:

Interpolation scheme for evaluation of

in the case

(

)

| ≥ |

(

)

. Left: Black solid circles indicate the

matching points that define the set

. 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.

3.3 Proximity map

As described below in Section 3.4, the 2D FC algorithm interpolates the Fourier continuation values on

B,C

onto the Cartesian mesh

. The interpolation algorithm used in that section relies on a certain

“proximity map”

∩

→

B,C

which associates a curvilinear grid point

p,q

(

i,j

)

∈

B,C

in the “proximity” of each given Cartesian grid point

i,j

. The proximity function we use is obtained by

first associating to each curvilinear discretization point

p,q

the nearest Cartesian point, a procedure that

results in a set

⊆

(

∩

)

B,C

of pairs of points, one in the Cartesian grid and the other in the

curvilinear grid. (The initial set

can easily be obtained by using the “integer part”

floor

(

)

and the

ceil

(

)

operators.) The set

is then modified 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 un-associated. The resulting set

defines the desired function.

3.4 FC Values on the Cartesian Grid

Once Fourier continuation values on

B,C

have been obtained, per the procedure presented in Section 3.2,

the two-dimensional FC scheme can be completed by (a) Interpolation onto the set

∩

of outer

Cartesian grid points; and, (b) Subsequent evaluation of the corresponding Fourier coefficients in equa-

tion Equation (3.1) by means of an FFT. (Note that since the continuation function

is a smooth function

which vanishes outside

Ω

∪

, this function can be viewed as the restriction of a smooth and bi-periodic

function with periodicity rectangle

—whose Fourier series approximates

, and therefore

, accurately.)

The efficiency of the interpolation scheme is of the utmost importance in this context—since interpolation

to a relatively large set

∩

of Cartesian points is necessary. An accurate and efficient interpolation

strategy is obtained by combining two one-dimensional (local) interpolation procedures based on nearby

normal 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 contin-

uation function value utilizing the parameter value

just mentioned and the continuation function values

at the points in

B,C

around the given Cartesian point. A detailed description of the combined interpola-

tion methodology is presented in what follows. Specifically, we describe the strategy we use to interpolate

the continuation function onto each point

i,j

∈

∩

and, to do this, we first obtain the foot

(

)

of this point. Using the proximity map

described in Section 3.3, the algorithm utilizes the curvilinear

discretization point

p,q

(

i,j

)

∈

B,C

as well as the corresponding boundary discretization parameter

value

; according to equation (3.6), the point

(

)

equals the foot of the normal passing through

p,q

(

) =

(

p,q

)

. The algorithm then seeks approximation of the foot

(

)

and the corresponding parame-

ter value

via a preliminary interpolation step, as indicated in what follows. The foot

(

)

and the

corresponding parameter value

are then used by the algorithm to produce the desired interpolated

continuation value at

In order to obtain

(

)

the algorithm uses the

boundary parameter values in the set

{

−

,θ

−

,...,θ

,θ

,...,θ

}⊂

(3.8)

around

(where

+ 1 =

, and where

is odd and

+ 1

is even.

Parameter values

with negative values of

, which may arise in Equation (3.8), are interpreted by

periodicity:

The algorithm then utilizes the line

⊥

passing through

that is orthogonal to the normal vector

(

)

(see left image in Figure 4), together with the parametrization

(

)

⊥

, where the parameter

represents the signed distance of the points on

⊥

from the point

. Clearly, then,

(0) =

. Each

point of intersection of

⊥

with the normals

(

)

(

∈

), on the other hand, equals

(

)

where

denotes the signed distance between

and the corresponding intersection point. Thus, defining the

function

(

)

, where

(

)

gives the parameter value of the foot of the normal through the point

(

)

, we clearly have

(

);

−

≤

(3.9)

It follows that a 1D interpolation procedure on the function

(

)

can be used to obtain the desired

approximation of the value

(0)

of the parameter corresponding to the foot of the point

i,j

(

) =

(

) =

(

(0))

Once we have the corresponding foot parameter value

for the given Cartesian point

, the distance

of the point

to the boundary

is easily computed. Let

⊂

be a set of

boundary parameter

values, similar to the set

defined in ( 3.8), but with

“re-centered” around

. In order to obtain

the continuation function value at the point

using the continuation function values on

B,C

, we employ

a local two-dimensional polynomial interpolation scheme based on the set

and the distance

, 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

from the boundary

along

the normal grid lines in

B,C

that correspond to boundary parameter values in the set

. Each one

of these values is obtained via one-dimensional interpolation of the continuation function values on

B,C

along the corresponding normal grid line in

B,C

. The desired continuation value at the point

, then,

is obtained via a final one-dimensional degree-(

−

) interpolation step based on the parameter set

and values at the “blue-asterisk” points just obtained. Finally, by applying the two-dimensional FFT to

Figure 4:

Interpolation schemes utilized to obtain the continuation function values on

∩

on the basis of the continuation

values on

B,C

. Left: Evaluation of the boundary parameter value for the foot of the normal line passing through

, 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

⊥

. Right: Interpolation of continuation values from the curvilinear mesh to a point

on the Cartesian

mesh.

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 bound-

ary data is provided in the form of values of the given function—which corresponds to the the Dirichlet

boundary data in the PDE context. But the approach is also applicable in cases for which the boundary

data is given as the normal derivative of the function, (Neumann boundary data), or even a combination

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 available, the two-dimensional Fourier

continuation 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 2D).

As in the 1D case, prior to the FFT procedure the grid

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 specially fast evaluations by means

of the fast Fourier transform.

Figure 5:

Numerical errors in

log

scale (left graph) and computing times required (right graph) in the 2D-FC approxima-

tion of the function

in (3.10) in the setting of Example 3.1. The interpolating polynomial degree

+ 3

was used in all

cases. The integer

Ω

equals the number of spatial grid points used over the diameter of the disc. Times reported correspond

to averages over

runs.

3.5 Summary of the 2D-FC procedure

This section presents a summary of the 2D-FC procedure described in the Sections 3.1 through 3.4 for a

function

given on a uniform Cartesian grid

∩

Ω

within the domain of definition

Ω

, where

is a Cartesian

mesh over the rectangle

containing both

Ω

and the near-boundary outer region

; see Section 3.1 and,

in particular, Figure 2. The construction of the continuation function

for the given function

relies

on use of three main parameters associated with the 1D blending-to-zero approach presented in Section 2,

namely

(number of points in the set

(number of unrefined discrete continuation points), and

(oversampling factor for the 1D blending-to-zero FC procedure), together with the parameters

(refinement factor for the discrete continuation points, Remark 2.2 and Section 3.1), and,

−

(degree

of the interpolating polynomials, Sections 3.2 and 3.4). Additionally, the 2D-FC procedure utilizes the

precomputed matrices

and

, which, with reference to Section 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 (comprising points 1. through 4. below); and step

(b) A

“Cartesian Evaluation”

procedure (comprising points 5. through 7. below). Part (a) only concerns

geometry, and, for a given domain and configuration, could be produced once and used for evaluation

of Fourier continuations for many different functions, as is often necessary in practice. The full 2D-FC

algorithm thus proceeds according to the following seven-steps:

1. Discretize the boundary

using a smooth parametrization

{

(

) = (

(

)

(

)) : 0

≤

}

and the uniform discretization

{

0 =

< θ

,...,< θ

−

}

of the interval

]

(Section 3.1).

2. Using the discretization

, construct two curvilinear meshes

−

B,d

and

B,C

in the near-boundary

regions

−

(within

Ω

) and

(outside

Ω

), respectively (equations (3.5) and (3.6)). Note that the

discrete boundary points

(

)

with

∈

are common to both

−

B,d

and

B,C

3. Determine the set

∩

of Cartesian grid points and construct the proximity map

∩

→

B,C

(Section 3.3).

4. For all

∈

∩

obtain the the parameter value

∈

of the foot

(

)

of the normal through

(Section 3.4 and the left image in Figure 4).

5. For each normal grid line (inherent in

−

B,d

and

B,C

) given by the discretization

, compute the

blending-to-zero function values along that normal (Section 3.2).

6. For all the points

∈

∩

, obtain the continuation function value at

by local 2D interpolation

(Section 3.4 and the right image in Figure 4).

7. Apply the two-dimensional FFT once to the continuation function values to obtain the desired Fourier

series in (3.1).

3.6 2D-FC approximation: Numerical Results

This section demonstrates the accuracy and efficiency of the proposed 2D-FC method. Use of the 2D-

FC method requires selection of specific values for each one of the following parameters (all of which are

introduced in Sections 2 and 3):

•

: number of points in the boundary section (Section 2.2).

•

: number of continuation points (Section 2.2).

•

: number of zero matching points (Section 2.2).

•

: oversampling factor used in the oversampled matching procedure (Section 2.2).

•

: refinement factor along the normal directions in

B,C

(Remark 2.2).

•

: the smallest rectangle containing

Ω

∪

(Section 3).

•

: number of points in the uniform spatial grid

(Section 3).

•

: number of points in the boundary discretization (Section 3).

•

−

: interpolating polynomial degree (Sections 3.2 and 3.4).

All the errors reported in this section were computed on a Cartesian grid of step size

within

Ω

. In

all of the numerical examples considered in this article the parameter selections were made in accordance

with Remark 3.3. The computer system used, in turn, is described in Remark 3.4.

Remark 3.3 (Parameter selections).

For a given step-size

in the two-dimensional Cartesian grid

the normal and the boundary step-sizes

and

(Section 3) were taken to coincide with

The parameter values

= 27

= 20

= 12

and

= 6

were used in the evaluation of the matrices

and

(see Section 2.2 and Remark 2.2). And, finally, with exception of the interpolation-degree

experiments presented in Example 3.2 (Table 2) and Example 4.2 (Table 3), the interpolating-polynomial

degree

(

−

1) = (

+ 2)

was used for the various matching-point numbers

considered.

Remark 3.4 (Computer system).

All of the numerical results reported in this paper were run on a

single core of a 3.40 GHz Intel Core i7-6700 processor with 15.4Gb of 2133 MHz memory.

Example 3.1 (Performance and efficiency of the 2D-FC method).

In our first example we consider

a problem of FC approximation of the function

: Ω

→

given by

(

x,y

) =

−

sin(5

πx

) sin(5

πy

)

(3.10)