DECOMPOSITION AND CONFORMAL MAPPING TECHNIQUES FOR THE QUADRATURE OF NEARLY SINGULAR INTEGRALS WILLIAM MITCHELL ABBIE NATKIN PAIGE ROBERTSON MARIKA SULLIVAN

2025-05-06 0 0 3.69MB 24 页 10玖币

侵权投诉

DECOMPOSITION AND CONFORMAL MAPPING TECHNIQUES

FOR THE QUADRATURE OF NEARLY SINGULAR INTEGRALS

WILLIAM MITCHELL, ABBIE NATKIN, PAIGE ROBERTSON, MARIKA SULLIVAN,

XUECHEN YU, AND CHENXIN ZHU

Abstract. Gauss-Legendre quadrature and the trapezoidal rule are power-

ful tools for numerical integration of analytic functions. For nearly singular

problems, however, these standard methods become unacceptably slow. We

discuss and generalize some existing methods for improving on these schemes

when the location of the nearby singularity is known. We conclude with an

application to some nearly singular surface integrals of viscous ﬂow. numerical

quadrature, nearly singular integral, Stokes ﬂow

1. Introduction

It is usually easy to compute deﬁnite integrals when the integrand is analytic.

Gauss-Legendre and Clenshaw-Curtis quadratures generally converge rapidly for

aperiodic problems, while the humble trapezoid rule is equally powerful for mea-

suring the area under one period of a periodic integrand. Unhappily, this is not

the case for the class of problems known as nearly singular integrals, wherein the

integrand fails to be analytic somewhere near but not on the interval of integra-

tion in the complex plane. The trapezoidal and Gauss-Legendre rules still provide

geometric convergence – that is, the logarithm of the error decreases linearly with

the number of quadrature nodes – but with an unacceptably small slope. Stan-

dard theorems relate this slope to the domain of analyticity of the integrand. In

this paper we survey the existing methods of accelerating the convergence of these

standard methods and we introduce several new ones, assuming that the location

of the nearby singularity is known, but without using any additional information

on the nature of the singularity.

While our motivation comes from the quadrature problem in the context of

solving integral equations, some very similar issues also arise in the context of

interpolation from samples of a nearly singular function. This is an important

component of spectral methods for diﬀerential equations when the desired solution

has abrupt fronts or peaks. In fact, some of the same acceleration techniques have

been independently discovered by researchers in the two communities. For example,

Johnston and Elliot’s hyperbolic sine transformation [7] was also discovered by Tee

and Trefethen [14], while Jafari’s transformation [6] is very similar to the repeated

sinh method of [3]. One of our aims is to assemble the results from these two

communities in one place.

We focus on two families of acceleration strategies. The ﬁrst group of methods

split the domain into carefully chosen subintervals and then solve subproblems on

Date: Received: date / Accepted: date.

arXiv:2210.09954v1 [math.NA] 18 Oct 2022

2 MITCHELL ET AL.

each of them. The other strategy that we consider is to change variables with a

complex-analytic transformation so that the singularity lies farther away.

An early example of a splitting method for aperiodic problems was given by

Ma and Kamiya [10], who subdivide at the real part of the singularity, assum-

ing that this lies within the integration interval and the problem is aperiodic (in

fact they subsequently use an exponential change of variables for each of the new

subproblems, thereby combining both of the principal strategies considered here).

Subsequently, Driscoll and Weideman [2] gave a formula for the optimal splitting

location in terms of the location of the singularity, again for the aperiodic case.

Their formula is especially useful in cases where the singularity is near the end-

point of the interval in the complex plane or on the real line outside the interval,

where simply using the real part of the singularity is ineﬀective or impossible. We

develop an analogous method for periodic problems, replacing the trapezoid rule

for a full period with individual Gauss-Legendre integrations on two subintervals

of diﬀerent sizes.

Splitting methods are extremely simple and, as we demonstrate, can yield dra-

matic improvements in the convergence rate. However, we can ﬁnd even better

convergence rates using methods that avoid dividing the available quadrature nodes

among two or more subproblems. There are a vast array of possibilities for con-

formal mappings that distort the neighborhood of a real interval while ﬁxing the

endpoints and remaining real-valued and monotone within the interval. The Jacobi

elliptic functions are a powerful tool for designing transformations that use all of

the analyticity we have assumed; we list existing methods for periodic and aperiodic

problems and give a new version for the aperiodic case when the singularity lies on

the real line outside the integration interval. This is a small extension of results

by Trefethen, Tee and Hale [13, 14, 5, 4]. However, the methods employing elliptic

functions are less robust than some elementary alternatives when the integrand has

other challenging features like additional distant singularities or rapid growth away

from the real line. We therefore list or introduce some good elementary alterna-

tives such as the sinh transformation for aperiodic problems [7, 14], our iterated

sine map for periodic problems, and our quadratic transformation for aperiodic

problems with real singularity.

The organization of the paper is as follows. We state and comment on the stan-

dard theorems describing the convergence of the trapezoidal and Gauss-Legendre

rules in Sec. 2. We consider periodic problems in Sec. 3. For aperiodic problems

we ﬁrst consider singularities oﬀ the real line in Sec. 4 and then singularities that

occur on the real line (but outside the interval of integration) in Sec. 5. We test

all of the methods on a suite of integrands with various properties in Sec. 6. As

an application, we then compute some nearly singular surface integrals arising in

Stokes ﬂow in Sec. 7, followed by concluding remarks.

2. Discussion of the standard theorems

For a periodic integrand, the convergence rate of the trapezoid rule depends on

the distance from the singularity to the real line.

Theorem 1 (Geometric convergence of trapezoid rule [18]).Suppose that a 2π-

periodic function fis analytic and satisﬁes |f(z)|< M within the strip |=(z)|< λ.

Then the diﬀerence between the integral I=R2π

0f(x)dx and its n-point trapezoidal

DECOMOPOSITION AND CONFORMAL MAPPING TECHNIQUES 3

rule approximation In= (2π/n)Pn

j=1 f(2πj/n)satisﬁes

(1) |I−In| ≤ 4πM

exp(λn)−1.

A similar theorem holds for Gauss-Legendre quadrature, with an ellipse replacing

the inﬁnite strip.

Theorem 2 (Geometric convergence of Gauss-Legendre quadrature [17]).Suppose

that fis analytic with |f(z)|< M on the interior of the Bernstein ellipse Eρ, whose

foci are ±1and whose semimajor and semiminor axis lengths sum to ρ > 1. Let C

be the constant C= 64ρ2/(15(1 −ρ−2)), let I=R1

−1f(x)dx be the exact value of

the integral, and let Inbe its n-point Gauss-Legendre rule approximation. Then

(2) |I−In| ≤ CM

ρ2n.

To employ this result in practice, we will often need to ﬁnd the value of the

ellipse parameter ρso that the ellipse passes through a given point z∈C\[−1,1].

We quote af Klinteberg and Barnett’s useful formula [8],

(3) ρ=ρ(z) = |z±pz2−1|,with sign chosen so that ρ > 1.

When we apply the preceding theorems, the twin goals of maximizing λor ρwhile

minimizing Mare in tension. If fis entire or has only distant singularities and nis

ﬁxed, this leads to an optimization problem for the value of λor ρthat will produce

the strongest statement from the theorem. Of course, the relative importance of M

decreases as ngrows. In the nearly singular case, where λ≈0 or ρ≈1, our priority

is to improve the convergence rate even if this results in a large increase in M. In

this paper we are assuming knowledge about the locations of the singularities of

the integrand, so λand ρare known. However, we make no assumptions about the

nature of those singularities or about the growth of |f(z)|away from the integration

interval, so Mis unavailable. Therefore we present a range of options instead of

searching for an optimal strategy. The relatively cautious methods we consider do

not make use of all of the analyticity we have assumed for f(z), and are therefore

less vulnerable to the danger of fast growth in |f(z)|away from the integration

interval. In contrast, the more aggressive methods make use of all of the analyticity

we have assumed (these generally involve the Jacobi elliptic functions, as we will

see below). The aggressive methods will boast better theoretical convergence rates,

although the errors may reach machine precision before they actually decrease at

the advertised rate. In contrast, the more cautious methods have (slightly) smaller

convergence rates but are more likely to actually achieve these rates for challenging

integrands.

3. Methods for periodic problems

We begin with the case where the integrand is 2π-periodic. We suppose that

fis analytic except at the points x= 2πk ±Bi for k∈Z, or along branch cuts

extending vertically from these points to inﬁnity. The ordinary trapezoidal rule will

be an eﬀective choice if Bis large, since Theorem 1 predicts errors of size e−Bn with

nfunction evaluations. We therefore assume that Bis small but nonzero, so the

integral is nearly singular, and we consider several methods to use our knowledge

of the domain of analyticity of fin order to improve on the performance of the

4 MITCHELL ET AL.

trapezoid rule. We begin with a decomposition method and then continue with

several conformal mappings. Numerical examples demonstrating these techniques

appear in Sec. 6.

3.1. Subdivision. One way to integrate over a full period of fis to choose an

appropriately small δ > 0 and then split the integral into subproblems on [−δ, δ]

and [δ, 2π−δ]. The subproblems are not periodic, so we apply Gauss-Legendre

integration for each of them. Rescaling the subintegrals linearly to [−1,1], we have

Zπ

−π

f(x)dx =δZ1

−1

f(δt)dt + (π−δ)Z1

−1

f(π+ (π−δ)w)dw.(4)

The two subintegrals are singular at t=Bi/δ and at w=Bi±π

π−δ, respectively. If

we want the overall integration procedure to converge as quickly as possible and

we assume nis large, we should choose δso that both subproblems have the same

asymptotic convergence rate. Equivalently, both Bi/δ and Bi±π

π−δshould lie on the

same ellipse with foci ±1 in the complex plane. Letting mdenote the semiminor

axis length, we have the equations

1 + m2+(B/δ)2

m2= 1,(π/(π−δ))2

1 + m2+(B/(π−δ))2

m2= 1.

After some simpliﬁcation we arrive at the cubic equation 0 = 2δ3+2B2δ−B2π. This

has a unique real solution which we write in terms of hyperbolic sines as follows:

(5) δ=2B

√3sinh 1

3arcsinh 3π√3

4B!!.

The error of the trapezoid rule on the original integral decays like e−Bn for large

n, where nis the number of quadrature nodes. With subdivision as in (4) and

Gauss-Legendre quadrature with n/2 nodes on each subinterval,1the total error

will decay like ρ−2(n/2) =ρ−n, where ρ= (B+√B2+δ2)/δ. Therefore, the

splitting procedure is advantageous when log ρ>B. We found numerically that

this holds for B < 0.95, although in practice, for ﬁnite values of n, it may be

advisable to use the trapezoid rule for slightly smaller values of B.

3.2. Conformal maps for periodic problems. A more powerful method for

periodic problems is to compute the integral using the n-point trapezoidal rule

following a transformation x(t) that preserves the interval [−π, π]:

(6) Zπ

−π

f(x)dx =Zπ

−π

f(x(t))x0(t)dt ≈2π

j=1

f(x(tj)) x0(tj)

where tj=−π+ 2πj/n. In view of Theorem 1, the error will decay like e−λn as

long as the new integrand f(x(t))x0(t) is periodic and analytic on the strip Sλ=

{t∈C:|=(t)|< λ}. We therefore seek a transformation x(t) which carries a wide

strip surrounding the real line in the complex t-plane into the domain of analyticity

of fin the complex x-plane, avoiding the branch cuts; another consideration is that

the derivative x0(t) should not have any singularities for t∈Sλ.

1We experimented with unequal distributions of nodes between the two subintervals but did

not see more than modest improvements. For odd n, one should assign the extra node to the

larger subinterval since errors there are multiplied by (π−δ) instead of δin (4).

DECOMOPOSITION AND CONFORMAL MAPPING TECHNIQUES 5

Two recent works have presented transformations with the general properties we

seek: Tee constructed one using the Jacobi elliptic functions [13], while Berrut and

Elefante gave an elementary formula based on a M¨obius transformation [1]. After

discussing these two methods, we propose a third, which we call the iterated sine

map. We illustrate the conformal maps in Figure 1 with B= 0.3 along with the

resulting predicted convergence rates. We then give numerical examples comparing

these three mappings, as well as the decomposition method from Section 3.1 and

the ordinary trapezoid rule, in Fig. 4.

3.2.1. The Jacobi amplitude map. Tee and Trefethen designed a conformal map

that carries a strip onto the doubly slit region of analyticity of f[13]. Here we

present a new derivation of the same formula using diﬀerential equations rather

than complex analysis. Our approach is based on the intuition that the factor

x0(t) should be small when xis close to the singularity. To respect periodicity,

we wrap the real line onto the unit circle and imagine the singularity as a nearby

point (1,0, B). We then assume that x0(t) is proportional to the distance from

(cos(x),sin(x),0) to the singularity, leading to the boundary value problem

(7) dx

dt =kq(1 −cos(x))2+ sin2(x) + B2, x(−π) = −π, x(π) = π

where the value of kis determined as part of the solution. The exact solution can

be obtained from the Jacobi amplitude function am(t, m), using the relation

(8) d

dt am(t, m) = q1−msin2am(t, m)= dn(t, m).

The solution of the BVP (7) is explicitly

(9) x(t) = −π+ 2 am π+t

πK4

4 + B2,4

4 + B2

where the real quarter-period K(m) is deﬁned for m < 1 by

(10) K(m) = Zπ/2

dθ

p1−msin2θ.

In practice, we suggest computing the derivative dx/dt using the Jacobi elliptic

function dn(t, m) rather than the diﬀerential equation (7):

(11) x0(t) = 2

πK4

4 + p2

ydn K4

4 + p2

y(t−1),4

4 + p2

y.

This leads to an improved convergence rate of

(12) λ=πKB2/(4 + B2)

K(4/(4 + B2)) .

The transformation (9) is the sum of the identity map and a 2π-periodic function,

and it carries the rectangle {x+iy :|x| ≤ π, |y|< λ}onto the doubly slit region

{x+iy :|x| ≤ πand |y|< pyif x= 0}, thereby using all of the analyticity that we

have assumed for f. This is the most aggressive option for the periodic problem.

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

DECOMPOSITIONANDCONFORMALMAPPINGTECHNIQUESFORTHEQUADRATUREOFNEARLYSINGULARINTEGRALSWILLIAMMITCHELL,ABBIENATKIN,PAIGEROBERTSON,MARIKASULLIVAN,XUECHENYU,ANDCHENXINZHUAbstract.Gauss-Legendrequadratureandthetrapezoidalrulearepower-fultoolsfornumericalintegrationofanalyticfunctions.Fornearlysingularprobl...

展开>> 收起<<

DECOMPOSITION AND CONFORMAL MAPPING TECHNIQUES FOR THE QUADRATURE OF NEARLY SINGULAR INTEGRALS WILLIAM MITCHELL ABBIE NATKIN PAIGE ROBERTSON MARIKA SULLIVAN.pdf

共24页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

DECOMPOSITION AND CONFORMAL MAPPING TECHNIQUES FOR THE QUADRATURE OF NEARLY SINGULAR INTEGRALS WILLIAM MITCHELL ABBIE NATKIN PAIGE ROBERTSON MARIKA SULLIVAN

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: