# Arbitrarily High-order Trapezoidal Rules for Functions with Fractional
Singularities in Two Dimensions

Senbao Jiang  Department of Applied Mathematics, Illinois Institute of
Technology  Xiaofan Li  Department of Applied Mathematics, Illinois Institute
of Technology

###### Abstract

In this paper, we introduce and analyze arbitrarily high-order quadrature
rules for evaluating the two-dimensional singular integrals of the forms

$\displaystyle
I_{i,j}=\int_{\mathbb{R}^{2}}\phi(x)\frac{x_{i}x_{j}}{|x|^{2+\alpha}}\operatorname{d}\\!{x},\quad
0<\alpha<2$

where $i,j\in\\{1,2\\}$ and $\phi\in C_{c}^{N}$ for $N\geq 2$. This type of
singular integrals and its quadrature rule appear in the numerical
discretization of fractional Laplacian in non-local Fokker-Planck Equations in
2D. The quadrature rules are trapezoidal rules equipped with correction
weights for points around singularity. We prove the order of convergence is
$2p+4-\alpha$, where $p\in\mathbb{N}_{0}$ is associated with total number of
correction weights. Although we work in 2D setting, we formulate definitions
and theorems in $n\in\mathbb{N}$ dimensions when appropriate for the sake of
generality. We present numerical experiments to validate the order of
convergence of the proposed modified quadrature rules.

##  1 Introduction

In this article we design and analyze arbitrarily high-order quadrature rules
for a class of weakly singular integrals in two dimensions, with an isolated
singularity at the origin. Such weakly singular integrals arises naturally in
non-local partial differential equations involving the fractional Laplacian.

There are different methods of numerical integration for weakly singular
integrals, such as singularity subtraction and singularity removing by change
of co-ordinates, see for example [6] [5] [4] [11] [17] [15]. The approach we
use here is to modify the trapezoidal rule around singularity in order to
achieve high-order convergence rate. This method is efficient as it introduces
correction weights for the trapezoidal rule near the singularity. The
correction weights are pre-computed. Once being stored, the modified
trapezoidal rule can be easily set up. It is well-known that the trapezoidal
rule is spectrally accurate for smooth functions with compact support or
periodic functions, while it is second order accurate for other functions.
There are two main sources of errors: one comes from singularities and the
other is from the boundaries. We focus on the error caused by singularity in
this paper.

Rokhlin [16] was first to propose singularity-corrected trapezoidal rule and
Alpert [3], Kapur and Rokhlin [9] further improves the method. Their
quadrature rules are designed for functions of the form
$f(x)=\phi(x)s(x)+\Psi(x)$ or $f(x)=\phi(x)s(x)$ in 1D, where
$\phi(x),\Psi(x)$ are regular functions and $s(x)$ is a singular function with
isolated singularity such as $s(x)=|x|^{\gamma},\gamma>-1$ or $s(x)=\log(x)$.
Aguilar and Chen [1, 2] designed singularity-boundary-corrected trapezoidal
rule for singular kernel $s(x)=\log(x)$ in 2D and $s(x)=\nicefrac{{1}}{{|x|}}$
in 3D.

The method we used for constructing increasingly accurate modified trapezoidal
rule is to require the modified rule exactly evaluate the integral of the
product of some monomials with increasing degree and the singular kernel. This
gives rise to a linear system of equations for the correction weights. The
more weights one is using, the larger the linear system is. In general, the
weights obtained from solving the linear systems depend on the mesh size $h$,
as in Keast and Lyness [10], Duan and Rokhlin [9]. In the latter paper, the
authors are able to derive analytical formula for the weights. Analytic
formula are not feasible for the weights relying on mesh size, especially when
the number of weights becomes large. In studying weakly singular integral
involving singular kernels of the forms $s(x)=|x|^{\gamma},\gamma>-1$ in 1D
and $s(x)=\nicefrac{{1}}{{|x|}}$ in 2D, Marin et al [14] designed a modified
trapezoidal rule and managed to obtain the weights independent of the mesh
size, by letting $h\to 0$. Those weights can then be pre-computed and stored.
They provided convergence analysis for the weights in 1D and the numerical
evidence of the convergence in 2D. They also gave a rigorous proof for the
rate of convergence of the quadrature rule in 1D and conditional proof for the
two dimensional counterpart. In fact, the difficulty for proving the linear
systems of the correction weights has a well-defined solution when $h\to 0$
lies in showing the linear systems are non-singular as the size of the systems
varies with the number of the correction weights. The same authors also
applied the quadrature rule to study the fundamental solution of the Stokes
equations in [13].

In [8], Ha investigates numerically non-local Fokker-Planck equation in two
dimensions

$\displaystyle u_{t}(x,t)$
$\displaystyle=-\partial_{1}(b_{1}(x)u(x,t))-\partial_{2}(b_{2}(x)u(x,t))+\frac{1}{2}\sigma_{11}^{2}\partial_{1}^{2}u(x,t)+\frac{1}{2}\sigma_{22}^{2}\partial_{2}^{2}u(x,t)$
(1)
$\displaystyle+C_{2,\alpha}\int_{\mathbb{R}^{2}\setminus\\{0\\}}\frac{u(x,t)-u(x+y,t)}{|y|^{2+\alpha}}\operatorname{d}\\!{y},\
(t,x)\in[0,1]\times\Omega$ (2) $\displaystyle u(x,t)$ $\displaystyle\equiv 0,\
x\not\in\Omega.$ (3)

Here,

$\displaystyle(-\Delta)^{\frac{\alpha}{2}}u(x)$ $\displaystyle\coloneqq
C_{2,\alpha}\int_{\mathbb{R}^{2}\setminus\\{0\\}}\frac{u(x)-u(x+y)}{|y|^{2+\alpha}}\operatorname{d}\\!{y},\quad\alpha\in(0,2),$
(4)

is the fractional Laplacian and

$\displaystyle C_{d,\alpha}$
$\displaystyle\coloneqq\frac{2^{\alpha}\Gamma(\frac{\alpha+d}{2})}{\pi^{\frac{d}{2}}|\Gamma(-\frac{\alpha}{2})|},\
d\in\mathbb{N}.$ (6)

Numerical discretization of fractional Laplacian boils down to designing
accurate quadrature rules for the weakly singular integrals
$I_{i,j}=\int_{\mathbb{R}^{2}}\phi(x)\frac{x_{i}x_{j}}{|x|^{2+\alpha}}\operatorname{d}\\!{x}$
, where $\phi$ is a regular function and $i,j\in\\{1,2\\}$. Ha [8] adapted the
simplest form of the modified trapezoidal rule of Marin et al [14] for the
issue. The quadrature rule is effective in evaluating the fractional
Laplacian, however there is no rigorous convergence analysis or comprehensive
numerical investigation for the rule itself. For more details please refer to
[8]. For a comprehensive review of fractional Laplacian please refer to [12]
and references therein.

In this work, we propose the complete version of the modified trapezoidal rule
for the two dimensional weakly singular integrals of the forms:

$\displaystyle I_{i,j}$
$\displaystyle=\int_{\mathbb{R}^{2}}\phi(x)\frac{x_{i}x_{j}}{|x|^{2+\alpha}}\operatorname{d}\\!{x},\quad\alpha\in(0,2).$
(7)

Here, $\phi(x)$ is a function with compact support satisfying some smoothness
criteria and $i,j\in\\{1,2\\}$. We call $I_{11},I_{22}$ the _on-diagonal_
weakly singular integrals and $I_{12},I_{21}$ the _off-diagonal_ weakly
singular integrals. We provide theoretical analysis of the proposed quadrature
rules and verify them by numerical results. We prove the quadrature rule has
the order of convergence $2p+4-\alpha,\ p\geq 0$ for on-diagonal case and
$2p+2-\alpha,\ p\geq 2$ for off-diagonal case, where $p$ is associated with
total number of correction weights. Our correction weights are independent
with the mesh size. We are able to show the linear systems of the correction
weights are non-singular for arbitrary $p$. By using our method one could give
unconditional affirmation to the conditional proof of convergence rate in
[14]. In the proof of the order of convergence of our quadrature rule, we
generalize the ideas used in the proofs in [14].

We organize the paper as follows. In Section 2, we introduce the modified
trapezoidal rules for on/off-diagonal weakly singular integrals and their
associated matrix formulations for determining correction weights. We prove
some preliminary results and the main theorems in Section 3. We address the
technical non-singularity problem arising from matrix formulation of section
Section 2 in Section 4. In final Section 5 numerical results for order of
convergence and associated correction weights are presented for both cases. In
our presentation, we elaborate on on-diagonal case and outline those off-
diagonal counterparts.

##  2 The Corrected Trapezoidal Rules

###  2.1 Notations

Through-out this paper, we use $h$ for mesh size, lowercase letters such as
$x,y,z$ for either scalar or vector and bold uppercase letters such as
$\boldsymbol{A},\boldsymbol{B}$ for matrices. The Euclidean norms
$|x|\coloneqq\sqrt{\sum_{i=1}^{n}x_{i}^{2}}$ and
$|x|_{1}:=\sum_{i=1}^{n}|x_{i}|$. The natural number with zero
$\mathbb{N}_{0}\coloneqq\\{0,1,2,\cdots\\}$. Multi-index notations appeared
frequently in this paper. In particular, we use $p\leq q$ for
$p=(p_{1},\cdots,p_{n}),\ q=(q_{1},\cdots,q_{n})$ if and only if $p_{i}\leq
q_{i}$ for all $i$ and $x^{\gamma}\coloneqq\prod_{j}x_{j}^{\gamma_{j}}$ for
any multi-index $\gamma=(\gamma_{1},\cdots,\gamma_{n})$.

###  2.2 On-diagonal Corrected Trapezoidal Rule

For any compactly supported function $f$ on $\mathbb{R}^{n}$, the punctured-
hole trapezoidal rule is defined by

$\displaystyle T_{h}^{0}[f]\coloneqq
h^{n}\sum_{\beta\in\mathbb{Z}^{n}\setminus\\{0\\}}f(\beta h).$ (8)

Let $f=\phi\cdot s$, where $\phi\in C_{c}(\mathbb{R}^{n})$ and $s$ is a
singular function of the form

$\displaystyle s(x)\coloneqq\frac{x_{1}^{2}}{|x|^{n+\alpha}},\quad\forall
x\in\mathbb{R}^{n},\ \text{with}\ 0<\alpha<2.$ (9)

We have chosen $s$ to be associated with $x_{1}$, the results associated with
this $s$ can be easily extended to other $x_{i}$’s. We call such $s$ on-
diagonal singular kernel. Let $p\in\mathbb{N}_{0}$, we introduce a corrected
trapezoidal rule in $\mathbb{R}^{2}$

$\displaystyle Q_{h}^{p}[f]\coloneqq
T_{h}^{0}[f]+a(h)\sum_{\beta\in\mathcal{L}_{p}}\bar{\omega}_{\beta}\phi(\beta
h),$ (10)

where $\bar{\omega}_{\beta}$’s are the correction weights to be defined,
$\mathcal{L}_{p}=\\{\beta\in\mathbb{Z}^{2}:|\beta|_{1}\leq p\\}$. $a(h)$ is
the leading error term between the singular integral and punctured-hole
trapezoidal rule and it turns out to be $h^{2-\alpha}$. To see this, we can
roughly evaluate

$\displaystyle\int_{(-h,h)^{2}}|\phi(x)|s(x)\ \operatorname{d}\\!{x}$
$\displaystyle\sim\int_{B(0,h)}|\phi(x)|s(x)\ \operatorname{d}\\!{x}$
$\displaystyle\lesssim\int_{B(0,h)}s(x)\
\operatorname{d}\\!{x}\lesssim\int_{0}^{h}r^{1-\alpha}\
\operatorname{d}\\!{r}\sim h^{2-\alpha},$

where $\sim$ means on the same order as $h\to 0$. In Section 3 we will show
rigorously that $a(h)=h^{2-\alpha}$.

We now elaborate on how to determine the correction weights
$\bar{\omega}_{\beta}$. For each $h>0$ we require the corrected trapezoidal
rule $Q_{h}^{p}$ in Eq. 10 with weights $\omega_{\beta}(h)$ evaluate the
following integrals exactly

$\displaystyle\int_{\mathbb{R}^{2}}g(x)s(x)x^{2\xi}\operatorname{d}\\!{x}$
$\displaystyle=T_{h}^{0}[g\cdot s\cdot x^{2\xi}]$
$\displaystyle+a(h)\sum_{\beta\in\mathcal{L}_{p}}\omega_{\beta}(h)(\beta
h)^{2\xi}g(\beta h),\quad\forall\xi\in\mathbb{N}_{0}^{2}:|\xi|_{1}\leq p,$
(11)

where $g$ is a radially symmetric, smooth function with compact support such
that $g(x)\not\equiv 0$. By symmetry of the integrand $g\cdot s$, we impose
same symmetry on weights

$\displaystyle\bar{\omega}_{\beta_{1},\beta_{2}}=\bar{\omega}_{|\beta_{1}|,|\beta_{2}|},\quad\forall\beta=(\beta_{1},\beta_{2})\in\mathcal{L}_{p}.$
(12)

Figure 1 illustrates the grid of the correction weights for $p=6$ after
imposing symmetry. We say that the grid of the on-diagonal correction weights
for $p=N,N\in\mathbb{N}_{0}$ is made up of $N+1$ correction layers, each layer
is the set of points $\beta\in\mathbb{N}^{2}$ such that $|\beta|_{1}=k,\
k=0,\cdots,N$.

Figure 1: The location of the on-diagonal grid of correction weights after
imposing symmetry when $p=6$.

To facilitate further analysis, we define

$\displaystyle\mathcal{I}_{p}$ $\displaystyle\coloneqq\\{(q-m,m):0\leq q\leq
p,\ 0\leq m\leq q\\},$ (13) $\displaystyle\mathcal{G}_{q-m}^{m}$
$\displaystyle\coloneqq\\{(\pm(q-m),\pm m),(\mp(q-m),\pm m)\\},\quad\text{for
each }(q-m,m)\in\mathcal{I}_{p},$ (14) $\displaystyle\mathcal{G}$
$\displaystyle\coloneqq\\{\mathcal{G}_{q-m}^{m}:(q-m,m)\in\mathcal{I}_{p}\\}.$
(15)

It is clear that $\mathcal{I}_{p}=\\{\xi\in\mathbb{N}_{0}^{2}:|\xi|_{1}\leq
p\\}$ and $|\mathcal{G}|=|\mathcal{I}_{p}|=N_{p}$, where
$N_{p}\coloneqq\frac{(p+1)(p+2)}{2}$. We occasionally use
$\mathcal{G}_{\gamma}$ to denote $\mathcal{G}_{q-m}^{m}$ where
$\gamma=(q-m,m)\in\mathcal{I}_{p}$ for the sake of notational convenience. We
write the system of equations Section 2.2 as $N_{p}\times N_{p}$ linear system
for the weights $\omega_{\gamma}$

$\displaystyle
a(h)\sum_{\gamma\in\mathcal{I}_{p}}\omega_{\gamma}(h)\sum_{\beta\in\mathcal{G}_{\gamma}}(\beta
h)^{2\xi}g(\beta
h)=\int_{\mathbb{R}^{2}}g(x)s(x)x^{2\xi}\,\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot
s\cdot x^{2\xi}],\quad\xi\in\mathcal{I}_{p}.$ (16)

and re-write Eq. 10 into

$\displaystyle Q_{h}^{p}[f]=T_{h}^{0}[f]+A_{h}^{p}[\phi]$ (17)

where

$\displaystyle
A^{p}_{h}[\phi]=a(h)\sum_{q=0}^{p}\sum_{m=0}^{q}\bar{\omega}^{m}_{q-m}\sum_{\beta\in\mathcal{G}^{m}_{q-m}}\phi(\beta
h)=a(h)\sum_{\gamma\in\mathcal{I}_{p}}\bar{\omega}_{\gamma}\sum_{\beta\in\mathcal{G}_{\gamma}}\phi(\beta
h).$ (18)

To formulate Eq. 16 in matrix form, we re-index the set
$\mathcal{I}_{p}=(\gamma_{i})_{1\leq i\leq N_{p}}$. There are more than one
way to index the set $\mathcal{I}_{p}$ and the indexing plays an important
role in proving the linear system has unique solution. We specify the indexing
later. By one-to-one correspondence between $\mathcal{G}$ and
$\mathcal{I}_{p}$, we index $\mathcal{G}$ by
$(\mathcal{G}_{j}\coloneqq\mathcal{G}_{\gamma_{j}})_{1\leq j\leq N_{p}}$, for
each $\gamma_{j}\in\mathcal{I}_{p},1\leq j\leq N_{p}$. We denote
$\boldsymbol{K},\boldsymbol{G}\in\mathbb{R}^{N_{p}\times N_{p}}$ such that for
all $1\leq i,j\leq N_{p}$

$\displaystyle K_{i,j}=\sum_{\beta\in\mathcal{G}_{j}}\beta^{2\xi_{i}},\
G_{i,j}=g(\beta h)\delta_{i,j},\ \beta\in\mathcal{G}_{j}.$ (19)

We call $\boldsymbol{K}$ the coefficient matrix. We note that by radial
symmetry, $g(\beta h)$ is constant over each $\mathcal{G}_{j}$, thus the
choice of $\beta$ within $\mathcal{G}_{j}$ is arbitrary. Now, let

$\displaystyle\omega(h)=(\omega_{1}(h),\dots,\omega_{N_{p}}(h))^{T},$ (20)

the linear system Eq. 16 becomes

$\displaystyle\boldsymbol{K}\boldsymbol{G}(h)\omega(h)=c(h),$ (21)

where the right-hand side of Eq. 21

$\displaystyle c_{i}(h)$
$\displaystyle=\frac{1}{h^{2|\xi_{i}|_{1}}a(h)}\left(\int
g(x)s(x)x^{2\xi_{i}}\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot s\cdot
x^{2\xi_{i}}]\right),\quad 1\leq i\leq N_{p},$ (22) $\displaystyle c(h)$
$\displaystyle=(c_{1}(h),\dots,c_{N_{p}}(h))^{T}.$ (23)

We solve the linear system Eq. 21 for each $h>0$ and the coefficients
$\bar{\omega}^{m}_{q-m}:(q-m,m)\in\mathcal{I}_{p}$ are the limits of the
solution of the linear system Eq. 21: $\bar{\omega}=\lim_{h\to 0}\omega(h)$,
provided the following claims hold

  * •

the limit of right-hand side of equation Eq. 22 exists,

  * •

the limit of $\boldsymbol{G}(h)$ exists, denoted by $\boldsymbol{G}$,

  * •

matrices $\boldsymbol{K}$ and $\boldsymbol{G}$ are non-singular.

The first two claims are answered in next section while the non-singularity of
$\boldsymbol{K}$ is addressed in section Section 4.

###  2.3 Off-diagonal Corrected Trapezoidal Rule

In this subsection, we construct corrected trapezoidal rule for the off-
diagonal singular kernel. In $n$ dimensions, the off-diagonal singular kernel
is defined by

$\displaystyle s_{12}(x)=\frac{x_{1}x_{2}}{|x|^{n+\alpha}},\ 0<\alpha<2.$ (24)

In the absence of ambiguity, we recycle the notations used for the on-diagonal
case. Due to similarity with on-diagonal case, we will simplify the
description for off-diagonal case.

In two dimensions, we assign a correction weight $\omega_{\beta}$ on each grid
points $\beta\coloneqq(\beta_{1},\beta_{2})\in\mathbb{Z}^{2}$ such that
$\beta_{1}\beta_{2}\not=0$. Imposing symmetry of $s_{12}$ on correction
weights, we have

$\displaystyle\bar{\omega}_{\beta_{1},\beta_{2}}=\text{sgn}\\{\beta_{1}\beta_{2}\\}\bar{\omega}_{|\beta_{1}|,|\beta_{2}|},\quad\bar{\omega}_{\beta_{1},\beta_{2}}=\bar{\omega}_{\beta_{2},\beta_{1}}.$
(25)

Figure 2 illustrates the grid of the correction weights for $p=6$ after
imposing symmetry. We say that the grid of the off-diagonal correction weights
for $p=N,N\geq 2$ is made up of $N-1$ correction layers, each layer is the set
of points $\beta\in\mathbb{N}^{2}$ such that $|\beta|_{1}=k,\ k=2,\cdots,N$
and $\beta_{1}\geq\beta_{2}$.

Figure 2: The location of the off-diagonal grid of correction weights after
imposing symmetry when $p=6$.

For all $p\in\mathbb{N},p\geq 2$, we denote

$\displaystyle\mathcal{G}_{q-m,m}=$ $\displaystyle\big{\\{}(\pm(q-m),\pm
m),(\mp(q-m),\pm m),$ $\displaystyle(\pm m,\pm(q-m)),(\mp
m,\pm(q-m))\big{\\}},\ 2\leq q\leq p,\ 1\leq
m\leq\left\lfloor\frac{q}{2}\right\rfloor,$ (26)
$\displaystyle\mathcal{I}_{p}=$ $\displaystyle\left\\{(q-m,m):2\leq q\leq
p,1\leq m\leq\left\lfloor\frac{q}{2}\right\rfloor\right\\}.$ (27)

We define the modified trapezoidal rule for a regular function $\phi$ in
$\mathbb{R}^{2}$ as

$\displaystyle Q_{h}^{1}[\phi\cdot s_{12}]$ $\displaystyle\coloneqq
T_{h}^{0}[\phi\cdot s_{12}],$ $\displaystyle Q_{h}^{p}[\phi\cdot s_{12}]$
$\displaystyle\coloneqq T_{h}^{0}[\phi\cdot
s_{12}]+h^{2-\alpha}\sum_{(q-m,m)\in\mathcal{I}_{p}}\bar{\omega}_{q-m,m}\sum_{\beta\in\mathcal{G}_{q-m,m}}\text{sgn}(\beta_{1}\beta_{2})\phi(\beta
h),$ $\displaystyle\eqqcolon T_{h}^{0}[\phi\cdot s_{12}]+A_{h}^{p}[\phi],\
\text{if }p\geq 2.$ (28)

where $T_{h}^{0}[f]$ means punctured-hole trapezoidal rule applied on function
$f$ defined at Eq. 8. To determine the correction weights, we introduce the
linear system: for each $\xi\in\mathcal{I}_{p}$

$h^{2-\alpha}\sum_{(q-m,m)\in\mathcal{I}_{p}}\omega_{q-m,m}(h)\sum_{\beta\in\mathcal{G}_{q-m,m}}\text{sgn}(\beta_{1}\beta_{2})(\beta
h)^{2\xi-e_{1}-e_{2}}g(\beta h)\\\
=\int_{\mathbb{R}^{2}\setminus\\{0\\}}g(x)s_{12}(x)x^{2\xi-e_{1}-e_{2}}\
\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot s_{12}\cdot x^{2\xi-e_{1}-e_{2}}],$
(29)

where $g$ is a smooth, radially symmetric, compactly supported function on
$\mathbb{R}^{2}$. We index $\mathcal{I}_{p}=(\gamma_{j})_{1\leq j\leq
N_{p}},N_{p}:=|\mathcal{I}_{p}|$, which will be specified later. Then by using
the same indexing we can write $\\{\mathcal{G}_{q-m,m}:2\leq q\leq p,1\leq
m\leq\lfloor\frac{q}{2}\rfloor\\}=\\{\mathcal{G}_{\gamma_{j}}\\}_{j}$. In this
case the matrix $\boldsymbol{K}\in\mathbb{R}^{N_{p}\times N_{p}}$ is re-
defined as

$\displaystyle
K_{i,j}=\sum_{\beta\in\mathcal{G}_{\gamma_{j}}}\text{sgn}(\beta_{1}\beta_{2})\beta^{2\xi_{i}-e_{1}-e_{2}},$
(30)

For the off-diagonal kernel $s_{12}$, the linear system Eq. 21 has a new
$\boldsymbol{K}$ and $c(h)$

$\displaystyle\boldsymbol{K}\boldsymbol{G}(h)\omega(h)$ $\displaystyle=c(h),$
(31) $\displaystyle c_{j}(h)$
$\displaystyle=\frac{1}{h^{2|\xi_{j}|_{1}-\alpha}}\left(\int g\cdot
s_{12}\cdot x^{2\xi_{j}-e_{1}-e_{2}}\ \operatorname{d}\\!{x}-T_{h}^{0}[g\cdot
s_{12}\cdot x^{2\xi_{j}-e_{1}-e_{2}}]\right).$ (32)

In Sections 3, LABEL: and 4 we will address the existence of the limit
$\lim_{h\to 0}c(h)$, the convergence of $\lim_{h\to 0}\omega(h)$ and non-
singularity of coefficient matrix, respectively.

##  3 The Analysis of Orders of Convergence

###  3.1 On-diagonal Analysis

In this subsection, we show the convergence order of corrected trapezoidal
rule $Q^{p}_{h}$ for on-diagonal singular integral in Eq. 17 is $2p+4-\alpha$
via the following main theorem:

######  Theorem 3.1.

Let $p\in\mathbb{N}_{0}$. Given $\phi\in C_{c}^{N}(\mathbb{R}^{2})$ such that
$N\geq 2p+4$ and the on-diagonal singular kernel
$s(x)=\frac{x_{1}^{2}}{|x|^{2+\alpha}},\ x\in\mathbb{R}^{2}$, $0<\alpha<2$.
Denote $f\coloneqq\phi\cdot s$, let $g\in C_{c}^{\infty}(\mathbb{R}^{2})$ such
that $g$ is radially symmetric, $g(0)=1$ and $\partial^{k}g(0)=0$ for all
$|k|_{1}=1,\dots,2p+1$. Let $\omega(h)$ be the solution of Eq. 21 for this
$g$, then $\omega(h)$ converges to some $\bar{\omega}\in\mathbb{R}^{N_{p}}$
such that

$\displaystyle|\omega(h)-\bar{\omega}|_{2}=O(h^{2p+2}).$ (33)

Moreover, there exists $C>0$ such that

$\displaystyle\left|Q_{h}^{p}[f]-\int_{\mathbb{R}^{2}}f(x)\
\operatorname{d}\\!{x}\right|\leq Ch^{2p+4-\alpha}.$ (34)

We will introduce some definitions, lemmas and a theorem before we prove the
main theorem Theorem 3.1.

######  Definition 3.1.

Given $f\in L^{1}(\mathbb{R}^{n})$, we define

$\displaystyle\widehat{f}(\xi)=\int_{\mathbb{R}^{n}}f(x)\exp{(-2\pi i\xi\cdot
x)}\operatorname{d}\\!{x}.$ (35)

$\widehat{f}$ is the Fourier transform of $f$.

We state the Poisson Summation formula for later use, for a proof, see [7].

######  Lemma 3.2 (Poisson Summation Formula).

Let $f$ be a continuous function on $\mathbb{R}^{n}$ which satisfies

$\displaystyle|f(x)|\leq\frac{C}{(1+|x|)^{n+\delta}}\quad\text{for some}\
C,\delta>0\ \text{and for all}\ x\in\mathbb{R}^{n},$ (36)

and whose Fourier transform $\widehat{f}$ restricted on $\mathbb{Z}^{n}$
satisfies

$\displaystyle\sum_{m\in\mathbb{Z}^{n}}|\widehat{f}(m)|<\infty,$ (37)

then

$\displaystyle\sum_{m\in\mathbb{Z}^{n}}\widehat{f}(m)=\sum_{k\in\mathbb{Z}^{n}}f(k).$
(38)

It is not hard to verify that if $f\in C_{c}^{n+1}(\mathbb{R}^{n})$, then $f$
satisfies the hypotheses of Poisson Summation formula, in fact,

$\displaystyle|f(x)|\leq\left(\frac{1+|R|}{1+|x|}\right)^{n+1}\|f\|_{\infty},\
\text{where}\ supp(f)\subset B_{R}(0)\ \text{for some}\ R>0,$ (39)

and for fixed $m\in\mathbb{Z}^{n}\setminus\\{0\\}$, choose
$|m_{j}|=\max_{i=1,\cdots,n}|m_{i}|$, then by using integration by part with
respect to $x_{j}$ $(n+1)$-times and the property of being compactly
supported, we obtain

$\displaystyle\widehat{f}(m)=(-1)^{n+1}\int_{\mathbb{R}^{n}}(\partial^{n+1}_{j}f)(x)\frac{\exp\\{-2\pi
ix\cdot m\\}}{(-2\pi im_{j})^{n+1}}\ \operatorname{d}\\!{x},$ (40)

note that $|m|\leq\sqrt{n}|m_{j}|$, we immediately have

$\displaystyle|\widehat{f}(m)|\leq\left(\frac{\sqrt{n}}{2\pi}\right)^{n+1}\frac{\max_{|l|_{1}=n+1}\|\partial^{l}f\|_{\mathcal{L}_{1}}}{|m|^{n+1}}.$
(41)

The right hand side of above inequality is summable over
$\mathbb{Z}^{n}\setminus\\{0\\}$. Thus the condition Eq. 37 is satisfied.

Now, we let $f\coloneqq g\cdot s$ where $g\in C_{c}^{\infty}(\mathbb{R}^{n})$
and $s(x)$ is defined in Eq. 9. Introducing a radially symmetric cut-off
function $\Psi\in C^{\infty}(\mathbb{R}^{n})$ such that

$\displaystyle\Psi(x)=\left\\{\begin{array}[]{cc}0,&|x|\leq\frac{1}{2}\\\
1,&|x|\geq 1\end{array}\right..$ (44)

A concrete example of such cut-off function is

$\displaystyle\Psi(x)=\left\\{\begin{array}[]{cc}0,&|x|\leq\frac{1}{2}\\\
\exp(\frac{-36}{9-16(1-|x|^{2})^{2}}+4),&\frac{1}{2}\leq|x|\leq 1\\\
1,&|x|\geq 1\end{array}\right..$ (48)

By the property of $\Psi$, $s(x)\Psi\left(\frac{x}{h}\right)$ can be
continuously extend to the origin by letting $s(0)\Psi(0)=\lim_{x\to
0}s(x)\Psi\left(\frac{x}{h}\right)=0$. Therefore,

$\displaystyle T_{h}^{0}[f]$
$\displaystyle=h^{n}\sum_{\beta\in\mathbb{Z}^{n}\setminus\\{0\\}}f(\beta h)$
$\displaystyle=h^{n}\sum_{\beta\in\mathbb{Z}^{n}\setminus\\{0\\}}f(\beta
h)\Psi(\beta)+h^{n}\underbrace{f(0)\Psi(0)}_{=0}=T_{h}\left[f(\cdot)\Psi\left(\frac{\cdot}{h}\right)\right].$
(49)

Hence we can split

$\displaystyle\int f(x)\operatorname{d}\\!{x}-T_{h}^{0}[f]$
$\displaystyle=\int
f(x)\left(1-\Psi\left(\frac{x}{h}\right)\right)\operatorname{d}\\!{x}$
$\displaystyle+\int
f(x)\Psi\left(\frac{x}{h}\right)\operatorname{d}\\!{x}-T_{h}\left[f(\cdot)\Psi\left(\frac{\cdot}{h}\right)\right].$
(50)

We now state and prove Theorem 3.4, which establishes the existence of the
limit $\lim_{h\to 0}c(h)$ from Eq. 22 and plays an important role in the proof
of Theorem 3.1. Before doing so, we need the following lemma.

######  Lemma 3.3.

Let $\eta\in\mathbb{N}_{0}^{n}$, the partial derivatives of order $k$
$(\in\mathbb{N}_{0}^{n})$ of
$x\rightarrow\frac{x^{\eta}x_{1}^{2}}{|x|^{n+\alpha}}$, $x\in\mathbb{R}^{n}$
are given as

$\displaystyle\partial^{k}\frac{x^{\eta}x_{1}^{2}}{|x|^{n+\alpha}}=\frac{1}{|x|^{n-2+\alpha-|\eta|_{1}+|k|_{1}}}P_{k,\eta}\left(\frac{x}{|x|}\right),$
(51)

where $P_{k,\eta}$ is a polynomial in $n$-variables with
$\text{deg}P_{k,\eta}=2+|\eta|_{1}+|k|_{1}$.

###### Proof.

For $k=0$,

$\displaystyle\partial^{0}\frac{x^{\eta}x_{1}^{2}}{|x|^{n+\alpha}}=\frac{1}{|x|^{n-2+\alpha-|\eta|_{1}}}\left(\frac{x}{|x|}\right)^{\eta+2e_{1}},$
(52)

which shows Eq. 52 is true for $k=0$ with $P_{0,\eta}(x)=x^{\eta+2e_{1}}$.
Suppose for $k\in\mathbb{N}_{0}^{n}$ the equation Eq. 51 is satisfied. We
define $P_{k,\eta}^{(j)}\coloneqq\partial_{x_{j}}P_{k,\eta}$ with
$j=1,\cdots,n$. Note that $P_{k,\eta}^{(j)}$ is a polynomial with one degree
less than $P_{k,\eta}$. Now,

$\displaystyle\partial^{k+e_{i}}\left\\{\frac{x^{\eta+2e_{1}}}{|x|^{n+\alpha}}\right\\}$
$\displaystyle=\partial_{i}\left\\{\frac{1}{|x|^{n-2+\alpha-|\eta|_{1}+|k|_{1}}}P_{k,\eta}\left(\frac{x}{|x|}\right)\right\\}$
$\displaystyle=-(n-2+\alpha-|\eta|_{1}+|k|_{1})\frac{1}{|x|^{n-1+\alpha-|\eta|_{1}+|k|_{1}}}\frac{x_{i}}{|x|}P_{k,\eta}\left(\frac{x}{|x|}\right)$
$\displaystyle+\frac{1}{|x|^{n-1+\alpha-|\eta|_{1}+|k|_{1}}}\left[P_{k,\eta}^{(i)}\left(\frac{x}{|x|}\right)\left(1-\frac{x_{i}^{2}}{|x|^{2}}\right)-\sum_{j\not=i}P_{k,\eta}^{(j)}\left(\frac{x}{|x|}\right)\frac{x_{i}x_{j}}{|x|^{2}}\right]$
$\displaystyle=\frac{1}{|x|^{n-1+\alpha-|\eta|_{1}+|k|_{1}}}\Bigg{[}-(n-2+\alpha-|\eta|_{1}+|k|_{1})\frac{x_{i}}{|x|}P_{k,\eta}\left(\frac{x}{|x|}\right)$
$\displaystyle+\left(1-\frac{x_{i}^{2}}{|x|^{2}}\right)P_{k,\eta}^{(i)}\left(\frac{x}{|x|}\right)-\sum_{j\not=i}\frac{x_{i}x_{j}}{|x|^{2}}P^{(j)}_{k,\eta}\left(\frac{x}{|x|}\right)\Bigg{]}$
$\displaystyle\eqqcolon\frac{1}{|x|^{n-2+\alpha-|\eta|_{1}+|k|_{1}+1}}P_{k+e_{i},\eta}\left(\frac{x}{|x|}\right).$
(53)

One can check that $\text{deg}(P_{k+e_{i},\eta})=3+|\eta|_{1}+|k|_{1}$. The
proof is then completed by induction. ∎

In this work we follow the convention that the generic constant $C$ may vary
from line to line in the proof.

######  Theorem 3.4.

Let $n\in\mathbb{N}$, $p\in\mathbb{N}_{0}$ and $M=\max(4p+4,n+1)$. Assume
$g\in C_{c}^{N}(\mathbb{R}^{n}):N\geq M$ such that $\partial^{k}g(0)=0$ for
all $|k|_{1}=1,\dots,2p+1$. For any $\xi\in\mathbb{N}_{0}^{n}$ with
$|\xi|_{1}\leq p$, we have

$\displaystyle\left|\int_{\mathbb{R}^{n}}g\cdot s\cdot x^{2\xi}\
\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot s\cdot
x^{2\xi}]-g(0)h^{2|\xi|_{1}+2-\alpha}c(2\xi,\alpha)\right|\leq
Ch^{2p+4-\alpha+2|\xi|_{1}},$ (54)

where $c(\xi,\alpha)\in\mathbb{R}$ and the constant $C$ depends only on $g$
and $p$.

###### Proof.

Let $f=g\cdot s_{\xi}$ where $s_{\xi}(x)=s(x)x^{2\xi}$. By the property of
$\Psi$ in Eq. 44, the first term in Section 3.1 can be computed by

$\displaystyle\int
f(x)\left(1-\Psi\left(\frac{x}{h}\right)\right)\operatorname{d}\\!{x}=\int_{B(0,h)}g(x)s_{\xi}(x)\left(1-\Psi\left(\frac{x}{h}\right)\right)\operatorname{d}\\!{x}$
$\displaystyle=h^{n}\int_{|x|\leq
1}g(hx)s_{\xi}(hx)\left(1-\Psi(x)\right)\operatorname{d}\\!{x}=h^{2|\xi|_{1}+2-\alpha}\int_{|x|\leq
1}g(hx)s_{\xi}(x)\left(1-\Psi(x)\right)\operatorname{d}\\!{x}.$ (55)

From the assumptions on $g$ and $2p+2|\xi|_{1}+3-\alpha>0$, Taylor’s theorem
and $n$-dimensional spherical co-ordinate transform, with radial direction
denoted by $r$, we have

$\displaystyle\left|\int_{|x|\leq
1}(g(hx)-g(0))s_{\xi}(x)(1-\Psi(x))\operatorname{d}\\!{x}\right|$
$\displaystyle=\left|\int_{|x|\leq
1}\left(g(hx)-\sum_{\begin{subarray}{c}k\in\mathbb{N}_{0}^{n},\\\ |k|_{1}\leq
2p+1\end{subarray}}\frac{\partial^{k}g(0)}{k!\partial
x^{k}}(hx)^{k}\right)s_{\xi}(x)(1-\Psi(x))\operatorname{d}\\!{x}\right|$
$\displaystyle\leq h^{2p+2}\sum_{\begin{subarray}{c}k\in\mathbb{N}_{0}^{n},\\\
|k|_{1}=2p+2\end{subarray}}\int_{|x|\leq
1}\left|\frac{\partial^{k}g(\rho)}{k!\partial
x^{k}}\right||x|^{2p+2|\xi|_{1}-n+4-\alpha}|1-\Psi(x)|\operatorname{d}\\!{x}$
$\displaystyle\leq C\max_{|k|_{1}=2p+2}\|\partial^{k}g\|_{\infty}\
h^{2p+2}\int_{0}^{1}r^{2p+2|\xi|_{1}+3-\alpha}\ \operatorname{d}\\!{r}\leq C\
h^{2p+2},$ (56)

where $\rho$ is some point between $hx$ and $0$ . Combining Section 3.1 and
Section 3.1, we have

$\displaystyle\left|\int
f(x)\left(1-\Psi\left(\frac{x}{h}\right)\right)\operatorname{d}\\!{x}-g(0)h^{2|\xi|_{1}+2-\alpha}\int_{|x|\leq
1}s_{\xi}(x)(1-\Psi(x))\operatorname{d}\\!{x}\right|\leq
Ch^{2p+2|\xi|_{1}+4-\alpha}.$ (57)

We define a dilation operator $(\tau^{a}\theta)(x)\coloneqq\theta(ax)$ where
$a>0$ and $\theta$ is any measurable function on $\mathbb{R}^{n}$. Noting that
$f_{\Psi,h}\coloneqq f(\cdot)\Psi(\frac{\cdot}{h})\in
C_{c}^{N}(\mathbb{R}^{n})$, then so is $\tau^{h}f_{\Psi,h}$. Hence, we apply
Poisson Summation formula Lemma 3.2 to $\tau^{h}f_{\Psi,h}$ and get

$\displaystyle\sum_{\beta\in\mathbb{Z}^{n}}(\tau^{h}f_{\Psi,h})(\beta)=\sum_{k\in\mathbb{Z}^{n}}\widehat{\tau^{h}f_{\Psi,h}}(k),$
(58)

thus

$\displaystyle T_{h}[f_{\Psi,h}]$
$\displaystyle=T_{h}\left[f(\cdot)\Psi(\frac{\cdot}{h})\right]=h^{n}\sum_{\beta\in\mathbb{Z}^{n}}f(h\beta)\Psi(\beta)=h^{n}\sum_{\beta\in\mathbb{Z}^{n}}(\tau^{h}f_{\Psi,h})(\beta)$
$\displaystyle=h^{n}\sum_{k\in\mathbb{Z}^{n}}\widehat{\tau^{h}f_{\Psi,h}}(k)=h^{n}\sum_{k\in\mathbb{Z}^{n}}h^{-n}\widehat{f_{\Psi,h}}\left(\frac{k}{h}\right)$
$\displaystyle=\int_{\mathbb{R}^{n}}f(x)\Psi\left(\frac{x}{h}\right)\operatorname{d}\\!{x}+\sum_{k\in\mathbb{Z}^{n}\setminus\\{0\\}}\widehat{f_{\Psi,h}}\left(\frac{k}{h}\right).$
(59)

Denote $I_{\xi}(h,k)\coloneqq\widehat{f_{\Psi,h}}\left(\frac{k}{h}\right)$,
from Eq. 58 and Section 3.1 we know

$\displaystyle\sum_{k\in\mathbb{Z}^{n}\setminus\\{0\\}}|I_{\xi}(h,k)|$
$\displaystyle<\infty,$ (60) $\displaystyle
T_{h}[f_{\Psi,h}]-\int_{\mathbb{R}^{n}}f(x)\Psi\left(\frac{x}{h}\right)\
\operatorname{d}\\!{x}$
$\displaystyle=\sum_{k\in\mathbb{Z}^{n}\setminus\\{0\\}}I_{\xi}(h,k).$ (61)

The following lemma provides an error estimate for each $I_{\xi}(h,k)$.

######  Lemma 3.5.

Assume all conditions stated in Theorem 3.4. For each
$k\in\mathbb{Z}^{n}\setminus\\{0\\}$, let $k_{j}=\max_{l=1,\cdots,n}|k_{l}|$.
There exists a function $W(\xi,k)$, independent of $h\ \text{and}\ g$, and
constant $C$, depending on $g$ and $M$, such that

$\displaystyle\left|I_{\xi}(h,k)-g(0)W(\xi,k)h^{2|\xi|_{1}+2-\alpha}(2\pi
k_{j})^{-M}\right|\leq C|k|^{-M}(h^{2p+4-\alpha+2|\xi|_{1}}+h^{M}).$ (62)

###### Proof of Lemma 3.5.

Denote $\hat{p}=2p+1$. By using $\exp(i\lambda
z)=\frac{1}{(i\lambda)^{M}}\frac{\operatorname{d}\\!{}^{M}}{\operatorname{d}\\!{z}^{M}}(\exp(i\lambda
z)),\ \lambda,z\in\mathbb{R}\setminus\\{0\\}$ and integration by parts
repeatedly,

$\displaystyle I_{\xi}(h,k)$
$\displaystyle=\int_{\mathbb{R}^{n}}\exp\left(-\frac{2\pi i}{h}k\cdot
x\right)g(x)s_{\xi}(x)\Psi\left(\frac{x}{h}\right)\operatorname{d}\\!{x}$
$\displaystyle=\left(\frac{-ih}{2\pi
k_{j}}\right)^{M}\int_{\mathbb{R}^{n}}\partial^{M}_{j}\left(\Psi\left(\frac{x}{h}\right)g(x)s_{\xi}(x)\right)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}.$ (63)

Here, $i$ is the imaginary unit and $\partial_{j}$ means taking derivative
with respect to $x_{j}$. Define

$\displaystyle
W(\xi,k)\coloneqq(-i)^{M}\int_{\mathbb{R}^{n}}\partial^{M}_{j}\left(\Psi\cdot
s_{\xi}\right)(x)\exp\left(-2\pi ik\cdot x\right)\operatorname{d}\\!{x}.$ (64)

In the rest of the proof we drop the $\xi$ in first argument of $W$ for
convenience. We note that $W$ does not depend on $h$ or $g$ and
$\partial_{j}^{l}\Psi\in C_{c}^{\infty}$ for any $l>0$. By the condition on
$M$ we know $M-(2|\xi|_{1}+2-\alpha)>M-(2p+2)\geq 2p+2>1$. It follows from
Lemma 3.3 that

$\displaystyle\int_{|x|>1}|\partial^{M}_{j}s_{\xi}(x)|\operatorname{d}\\!{x}$
$\displaystyle=\int_{|x|>1}\frac{1}{|x|^{n-2+\alpha-2|\xi|_{1}+M}}\left|P_{M,\xi}\left(\frac{x}{|x|}\right)\right|\operatorname{d}\\!{x}$
$\displaystyle\leq
C\int_{1}^{\infty}r^{1-\alpha+2|\xi|_{1}-M}\operatorname{d}\\!{r}<\infty.$
(65)

Therefore,

$\displaystyle|W(k)|\leq
C\left(\int_{\mathbb{R}^{n}}\frac{\Psi(x)}{|x|^{n-2+\alpha-2|\xi|_{1}+M}}\operatorname{d}\\!{x}+\sum_{0<l\leq
M}\int_{\mathbb{R}^{n}}\partial^{l}_{j}\Psi(x)\partial^{M-l}_{j}s_{\xi}(x)\operatorname{d}\\!{x}\right)<\infty.$
(66)

The estimate Eq. 66 shows the boundedness of $W$ is independent of $k$. Now we
re-scale the integral

$W(k)=(-i)^{M}h^{M-n}\int_{\mathbb{R}^{n}}\partial^{M}_{j}\left(\Psi\left(\frac{x}{h}\right)s_{\xi}\left(\frac{x}{h}\right)\right)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}\\\
=(-i)^{M}h^{M+\alpha-2-2|\xi|_{1}}\int_{\mathbb{R}^{n}}\partial^{M}_{j}\left(\Psi\left(\frac{x}{h}\right)s_{\xi}(x)\right)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}.$ (67)

Then, denote $r(x)\coloneqq\Psi(\nicefrac{{x}}{{h}})(g(x)-g(0))$, from Eq. 67

$\displaystyle I_{\xi}(h,k)-g(0)W(k)h^{2|\xi|_{1}+2-\alpha}(2\pi k_{j})^{-M}$
$\displaystyle=\left(\frac{-ih}{2\pi
k_{j}}\right)^{M}\int_{\mathbb{R}^{n}}\partial^{M}_{j}\left\\{r(x)s_{\xi}(x)\right\\}\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\ \operatorname{d}\\!{x}.$ (68)

It can be seen that for all $0\leq m\leq M$,

$\displaystyle\partial^{m}_{j}r(x)=\left\\{\begin{array}[]{ccc}0&,&0\leq|x|\leq\frac{h}{2},\\\
\sum_{0\leq l\leq
m}\binom{m}{l}\frac{1}{h^{m-l}}\partial^{m-l}_{j}\Psi\left(\frac{x}{h}\right)(\partial^{l}_{j}g(x)-\delta_{0,l}g(0))&,&\frac{h}{2}\leq|x|\leq
h\\\ \partial^{m}_{j}g(x)-\delta_{0,m}g(0)&,&|x|>h\end{array}\right..$ (72)

Here we used Kronecker delta
$\delta_{k_{1},k_{2}}=\mathbbm{1}_{k_{1}}(k_{2})$. Let
$\text{supp}(g)\subseteq B(0,R)$ for some $R>1$, then we have

$\displaystyle\int_{\mathbb{R}^{n}}\partial^{m}_{j}r(x)\partial_{j}^{M-m}s_{\xi}(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}$ $\displaystyle=\sum_{0\leq l\leq
m}\frac{\binom{m}{l}}{h^{m-l}}\int_{\frac{h}{2}\leq|x|\leq
h}\partial^{m-l}_{j}\Psi\left(\frac{x}{h}\right)\Big{[}\partial^{l}_{j}g(x)-\delta_{0,l}g(0)\Big{]}\partial^{M-m}_{j}s_{\xi}(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}$ $\displaystyle+\int_{h<|x|\leq
R}\Big{[}\partial^{m}_{j}g(x)-\delta_{0,m}g(0)\Big{]}\partial^{M-m}_{j}s_{\xi}(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}$
$\displaystyle-\delta_{0,m}g(0)\int_{|x|>R}\partial^{M-m}_{j}s_{\xi}(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}.$ (73)

Let $T_{m,1},T_{m,2},T_{m,3}$ be the first, second and third terms in RHS of
Section 3.1. Since $(\partial^{l}_{j}g)(0)=0$ for all $l=1,\cdots,\hat{p}$, we
have for all $l\leq\hat{p}$

$\displaystyle|\partial^{l}_{j}g(x)-(\partial^{l}_{j}g)(0)|\leq
C\max_{\beta\in\mathbb{N}_{0}^{n}:|\beta|_{1}=\hat{p}+1}\|\partial^{\beta}g\|_{\infty}|x|^{\hat{p}+1-l}.$
(74)

It follows that for all $0\leq m\leq M$

$\displaystyle|\partial^{m}_{j}g(x)-\delta_{0,m}(\partial^{m}_{j}g)(0)|\leq
C|x|^{\max\\{\hat{p}+1-m,0\\}}.$ (75)

For all $0\leq m\leq M$ and $0\leq l\leq m$, we estimate each term in
$T_{m,1}$, from Eq. 75 and Lemma 3.3,

$\displaystyle\frac{1}{h^{m-l}}\left|\int_{\frac{h}{2}\leq|x|\leq
h}\partial^{m-l}_{j}\Psi\left(\frac{x}{h}\right)\Big{[}\partial^{l}_{j}g(x)-\delta_{0,l}g(0)\Big{]}\partial^{M-m}_{j}s_{\xi}(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}\right|$
$\displaystyle\leq\frac{C}{h^{m-l}}h^{n}\int_{\frac{1}{2}\leq|x|\leq
1}|hx|^{-(n-2+\alpha-2|\xi|_{1}+M-m)+\max\\{\hat{p}+1-l,0\\}}\operatorname{d}\\!{x}$
$\displaystyle\leq
Ch^{\hat{p}+3-\alpha+2|\xi|_{1}-M}\int_{\frac{1}{2}}^{1}r^{\hat{p}+2-\alpha+2|\xi|_{1}+m+l-M}\operatorname{d}\\!{r}$
$\displaystyle\leq Ch^{\hat{p}+3-\alpha+2|\xi|_{1}-M},$ (76)

where we used $h^{\max\\{\gamma,0\\}}\leq h^{\gamma}$ if $0<h\leq 1$. It
follows that $|T_{m,1}|\leq Ch^{\hat{p}+3-\alpha+2|\xi|_{1}-M}$. By a similar
computation as in Section 3.1, we have

$\displaystyle|T_{m,2}|$ $\displaystyle=\left|\int_{h<|x|\leq
R}\Big{[}\partial^{m}_{j}g(x)-\delta_{0,m}g(0)\Big{]}\partial^{M-m}_{j}s_{\xi}(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}\right|$ $\displaystyle\leq
C\int_{h<|x|\leq
R}|x|^{\max\\{\hat{p}+1-m,0\\}-(n-2+\alpha-2|\xi|_{1}+M-m)}\operatorname{d}\\!{x}$
$\displaystyle\leq
C\int_{h}^{R}r^{\max(\hat{p}+1-m,0)+1-\alpha+2|\xi|_{1}+m-M}\operatorname{d}\\!{r}$
$\displaystyle\leq
C\int_{h}^{1}r^{\hat{p}+2-\alpha+2|\xi|_{1}-M}\operatorname{d}\\!{r}+C$
$\displaystyle\leq C(h^{\hat{p}+3-\alpha+2|\xi|_{1}-M}+1).$ (77)

We used $M-(\hat{p}+2-\alpha+2|\xi|_{1})\geq 1+\alpha$ in the last inequality.
Lastly, we estimate

$\displaystyle|T_{m,3}|=\left|\int_{|x|\geq
R}\partial^{M}_{j}s_{\xi}(x)\exp\left(\frac{-2\pi i}{h}k\cdot x\right)\
\operatorname{d}\\!{x}\right|\leq C\int_{|x|\geq
R}\frac{\operatorname{d}\\!{x}}{|x|^{n-2+\alpha-2|\xi|_{1}+M}}<\infty.$ (78)

Putting together the results from Section 3.1 to Eq. 78 and using
$|k|\leq\sqrt{n}|k_{j}|$, we obtain from Section 3.1

$\displaystyle|I_{\xi}(h,k)-g(0)W(\xi,k)h^{2|\xi|_{1}+2-\alpha}(2\pi
k_{j})^{-M}|$ $\displaystyle\leq
C\left(\frac{h}{|k_{j}|}\right)^{M}\sum_{m=0}^{M}(|T_{m,1}|+|T_{m,2}|+|T_{m,3}|)$
$\displaystyle\leq
C\left(\frac{h}{|k_{j}|}\right)^{M}(h^{\hat{p}+3-\alpha+2|\xi|_{1}-M}+1)$
$\displaystyle\leq C|k_{j}|^{-M}(h^{2p+4-\alpha+2|\xi|_{1}}+h^{M})\leq
C|k|^{-M}(h^{2p+4-\alpha+2|\xi|_{1}}+h^{M}).$ (79)

∎

We define

$\displaystyle c(\xi,\alpha)$
$\displaystyle=c_{1}(\xi,\alpha)-c_{2}(\xi,\alpha),$ (80)
$\displaystyle\text{where}\quad c_{1}(\xi,\alpha)$
$\displaystyle\coloneqq\int_{|x|\leq
1}s_{\xi}(x)(1-\Psi(x))\operatorname{d}\\!{x}<\infty,$ (81) $\displaystyle
c_{2}(\xi,\alpha)$
$\displaystyle\coloneqq\sum_{\begin{subarray}{c}k\in\mathbb{Z}^{n}\setminus\\{0\\},\\\
k_{j}=\max_{l}|k_{l}|\end{subarray}}W(\xi,k)(2\pi k_{j})^{-M}<\infty.$ (82)

Equation 82 holds since for each $k\in\mathbb{Z}^{n}\setminus\\{0\\}$,
$|k|\leq\sqrt{n}|k_{j}|$ and so

$\displaystyle\sum_{k\not=0}|k_{j}|^{-M}\lesssim\sum_{k\not=0}|k|^{-M}\leq\sum_{k\not=0}|k|^{-(n+1)}<\infty.$

Using Sections 3.1, 57, 3.1 and 61, Eqs. 80 to 82 and Lemma 3.5, we obtain

$\displaystyle\left|\int_{\mathbb{R}^{n}}f(x)\operatorname{d}\\!{x}-T_{h}^{0}[f]-g(0)h^{2-\alpha+2|\xi|_{1}}c(\xi,\alpha)\right|$
$\displaystyle\leq\left|\int_{\mathbb{R}^{n}}f(x)\Big{(}1-\Psi\Big{(}\frac{x}{h}\Big{)}\Big{)}\
\operatorname{d}\\!{x}-g(0)h^{2-\alpha+2|\xi|_{1}}c_{1}(\xi,\alpha)\right|$
$\displaystyle+\left|\int_{\mathbb{R}^{n}}f(x)\Psi\Big{(}\frac{x}{h}\Big{)}\operatorname{d}\\!{x}-T_{h}\Big{[}f(\cdot)\Psi\Big{(}\frac{\cdot}{h}\Big{)}\Big{]}+g(0)h^{2-\alpha+2|\xi|_{1}}c_{2}(\xi,\alpha)\right|$
$\displaystyle\overset{\lx@cref{creftype~refnum}{TrapzIntegralDiff}}{\leq}\left|\int_{\mathbb{R}^{n}}f(x)\Big{(}1-\Psi\Big{(}\frac{x}{h}\Big{)}\Big{)}dx-g(0)h^{2-\alpha+2|\xi|_{1}}c_{1}(\xi,\alpha)\right|+\left|g(0)h^{2-\alpha+2|\xi|_{1}}c_{2}(\xi,\alpha)-\sum_{k\in\mathbb{Z}^{n}\setminus\\{0\\}}I_{\xi}(h,k)\right|$
$\displaystyle\overset{\lx@cref{creftype~refnum}{FirstPartEstimate},\lx@cref{creftype~refnum}{Weights2}}{\leq}Ch^{2p+4-\alpha+2|\xi|_{1}}+\sum_{\begin{subarray}{c}k\in\mathbb{Z}^{n},k\not=0,\\\
k_{j}=\max_{l}|k_{l}|\end{subarray}}|I_{\xi}(h,k)-g(0)h^{2-\alpha+2|\xi|_{1}}W(\xi,k)(2\pi
k_{j})^{-M}|$ $\displaystyle\leq
Ch^{2p+4-\alpha+2|\xi|_{1}}+C(h^{2p+4-\alpha+2|\xi|_{1}}+h^{M})\sum_{\begin{subarray}{c}k\in\mathbb{Z}^{n},k\not=0,\\\
k_{j}=\max_{l}|k_{l}|\end{subarray}}|k_{j}|^{-M}$ $\displaystyle\leq
Ch^{2p+4-\alpha+2|\xi|_{1}}+Ch^{2p+4-\alpha+2|\xi|_{1}}\sum_{k\in\mathbb{Z}^{n},k\not=0}|k|^{-(n+1)}\leq
Ch^{2p+4-\alpha+2|\xi|_{1}}$ (83)

This concludes the proof. ∎

Before we prove the Theorem 3.1, we need to know the order of convergence of
the punctured-hole trapezoidal rule Eq. 8.

######  Corollary 3.6.

Let $p\in\mathbb{N}_{0}$. Assume $g\in C_{c}^{N}(\mathbb{R}^{n})$ such that
$N\geq\max\\{2p+4,n+1\\}$ and $(\partial^{k}g)(0)=0$ for all
$k\in\mathbb{N}_{0}^{n}:|k|_{1}\leq 2p+1$. Then

$\displaystyle\left|\int_{\mathbb{R}^{n}}g\cdot
s\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot s]\right|\leq Ch^{2p+4-\alpha}.$

###### Proof.

The proof is rather similar to that of Theorem 3.4. We only present the
essential parts. Let $\Psi$ be the cut-off function introduced in Eq. 44 and
$\hat{p}=2p+1$. By Sections 3.1, LABEL: and 61 for $\xi=0$, we have

$\displaystyle\int_{\mathbb{R}^{n}}g(x)s(x)\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot
s]$
$\displaystyle=\int_{\mathbb{R}^{n}}g(x)s(x)\left(1-\Psi\left(\frac{x}{h}\right)\right)\operatorname{d}\\!{x}-\sum_{k\in\mathbb{Z}^{n}\setminus\\{0\\}}I_{0}(h,k),$
(84) $\displaystyle\text{where }I_{0}(h,k)$
$\displaystyle=\int_{\mathbb{R}^{n}}\exp\left(\frac{-2\pi i}{h}k\cdot
x\right)g(x)s(x)\Psi\left(\frac{x}{h}\right)\operatorname{d}\\!{x}.$

By using $g(0)=0$ and property of $\Psi$, we have $|g(x)|\leq C|x|^{2p+2}$ and

$\displaystyle\left|\int_{\mathbb{R}^{n}}g(x)s(x)\left(1-\Psi\left(\frac{x}{h}\right)\right)\operatorname{d}\\!{x}\right|\leq
C\int_{|x|\leq h}|x|^{2p+4-n-\alpha}\operatorname{d}\\!{x}$ $\displaystyle\leq
C\int_{0}^{h}r^{2p+3-\alpha}\operatorname{d}\\!{r}\leq Ch^{2p+4-\alpha}.$ (85)

We now consider $I_{0}(h,k)$. Fix $k\in\mathbb{Z}^{n},k\not=0$, let
$k_{j}=\max_{i=1}^{n}|k_{i}|$. Repeated integration by parts with respect to
$x_{j}$ yields

$\displaystyle
I_{0}(h,k)=\left(\frac{-ih}{2\pi}\right)^{N}\frac{1}{k_{j}^{N}}\int_{\mathbb{R}^{n}}\partial_{j}^{N}\left(\Psi\left(\frac{x}{h}\right)g(x)s(x)\right)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}.$

For a fixed $m$ with $0\leq m\leq N$, let $r(x)=\Psi(\nicefrac{{x}}{{h}})g(x)$
and

$\displaystyle
I_{m}=\int_{\mathbb{R}^{n}}\partial_{j}^{m}r(x)\partial_{j}^{N-m}s(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}.$

Since

$\displaystyle\partial^{m}_{j}r(x)=\left\\{\begin{array}[]{ccc}0&,&0\leq|x|\leq\frac{1}{2}\\\
\sum_{0\leq l\leq
m}\binom{m}{l}\frac{1}{h^{m-l}}\partial^{m-l}_{j}\Psi\left(\frac{x}{h}\right)\partial^{l}_{j}g(x)&,&\frac{h}{2}\leq|x|\leq
h\\\ \partial^{m}_{j}g(x)&,&|x|>h\end{array}\right.,$

we have

$\displaystyle|I_{m}|$ $\displaystyle\leq C\sum_{0\leq l\leq
m}\frac{1}{h^{m-l}}\left|\int_{\frac{h}{2}\leq|x|\leq
h}\partial^{m-l}_{j}\Psi\left(\frac{x}{h}\right)\partial^{l}_{j}g(x)\partial^{N-m}_{j}s(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}\right|$
$\displaystyle+\left|\int_{h<|x|<R}\partial^{m}_{j}g(x)\partial^{N-m}_{j}s(x)\exp\left(\frac{-2\pi
i}{h}k\cdot x\right)\operatorname{d}\\!{x}\right|,$

where $supp(g)\subset B(0,R)$. The argument of $|I_{m}|\leq
C(h^{2p+4-\alpha-N}+1)$ follows exactly Sections 3.1, LABEL: and 3.1 by
replacing $M$ with $N$ and $\xi$ with $0$. Since $N\geq\max(2p+4,n+1)$, we
obtain

$\displaystyle|I(h,k)|$ $\displaystyle\leq
C\frac{1}{|k_{j}|^{N}}(h^{2p+4-\alpha}+h^{N})\leq
C\frac{1}{|k_{j}|^{N}}h^{2p+4-\alpha},$
$\displaystyle\sum_{k\in\mathbb{Z}^{n},k\not=0}|I(h,k)|$ $\displaystyle\leq
Ch^{2p+4-\alpha}\sum_{\begin{subarray}{c}k\in\mathbb{Z}^{n},k\not=0,\\\
k_{j}=\max_{l}|k_{l}|\end{subarray}}\frac{1}{|k_{j}|^{n+1}}\leq
Ch^{2p+4-\alpha}.$ (86)

Combining Eqs. 85, LABEL: and 86, by Eq. 84 we have

$\displaystyle\left|\int_{\mathbb{R}^{n}}g\cdot
s\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot s]\right|\leq Ch^{2p+4-\alpha}.$ (87)

∎

We are in the position to prove our Theorem 3.1.

###### Proof of Theorem 3.1.

We first show the limit of the solution $\omega(h)$ of the linear system Eq.
21 exists, i.e. $\lim_{h\to 0}\omega(h)=\bar{\omega}$. To this end, we apply
Theorem 3.4 to each $c_{i}(h)$ in Eq. 22 and let $a(h)=h^{2-\alpha}$, yielding

$\displaystyle c_{i}(h)$
$\displaystyle=\frac{1}{h^{2|\xi_{i}|_{1}+2-\alpha}}\left(\int
g(x)s(x)x^{2\xi_{i}}\ \operatorname{d}\\!{x}-T_{h}^{0}[g\cdot s\cdot
x^{2\xi_{i}}]\right)$ (88)
$\displaystyle=\frac{1}{h^{2|\xi_{i}|_{1}+2-\alpha}}\left(g(0)h^{2|\xi_{i}|_{1}+2-\alpha}c(\xi_{i},\alpha)+O(h^{2p+4-\alpha+2|\xi_{i}|_{1}})\right)$
$\displaystyle=c(\xi_{i},\alpha)+O(h^{2p+2})\underset{h\to
0}{\longrightarrow}c(\xi_{i},\alpha).$

Define

$\displaystyle\mathfrak{C}(\alpha)=(c(\xi_{1},\alpha),\dots,c(\xi_{N_{p}},\alpha)).$
(89)

It is also clear that

$\displaystyle\boldsymbol{G}(h)=\text{diag}(g(\xi_{1}h),\dots,g(\xi_{N_{p}}h)))\underset{h\to
0}{\longrightarrow}\boldsymbol{I}_{N_{p}\times N_{p}},$ (90)

and so

$\displaystyle\boldsymbol{G}(h)^{-1}=\text{diag}(g(\xi_{1}h)^{-1},\dots,g(\xi_{N_{p}}h)^{-1}))\underset{h\to
0}{\longrightarrow}\boldsymbol{I}_{N_{p}\times N_{p}}.$ (91)

In addition, $\|\boldsymbol{G}(h)^{-1}\|_{2}=\max\\{|g(\xi_{k}h)|^{-1}:1\leq
k\leq N_{p}\\}$ is bounded in a neighbourhood of origin. Define $\bar{\omega}$
to be the solution of the linear system

$\displaystyle\boldsymbol{K}\bar{\omega}=\mathfrak{C}(\alpha).$ (92)

Assuming $\boldsymbol{K}$ is non-singular, which will be proved in next
section, from Eq. 21 and Eq. 92, we obtain

$\displaystyle|\omega(h)-\bar{\omega}|_{2}$
$\displaystyle=|\boldsymbol{G}(h)^{-1}\boldsymbol{K}^{-1}c(h)-\boldsymbol{K}^{-1}\mathfrak{C}(\alpha)|_{2}$
$\displaystyle=|\boldsymbol{G}(h)^{-1}(\boldsymbol{I}-\boldsymbol{G}(h))\boldsymbol{K}^{-1}c(h)+\boldsymbol{K}^{-1}(c(h)-\mathfrak{C}(\alpha))|_{2}$
$\displaystyle\leq\|\boldsymbol{G}(h)^{-1}\|_{2}\
\|\boldsymbol{I}-\boldsymbol{G}(h)\|_{2}\ \|\boldsymbol{K}^{-1}\|_{2}\
|c(h)|_{2}+\|\boldsymbol{K}^{-1}\|_{2}\ |c(h)-\mathfrak{C}(\alpha)|_{2}.$ (93)

A Taylor expansion of $g$ around the origin gives

$\displaystyle|1-g(\xi_{k}h)|\leq Ch^{2p+2},\quad k=1,\dots,N_{p},$ (94)

and hence

$\displaystyle\|\boldsymbol{I}-\boldsymbol{G}(h)\|_{2}=\max\\{|1-g(\xi_{k}h)|:k=1,\dots,N_{p}\\}\leq
Ch^{2p+2}.$ (95)

From Eq. 88, Eq. 93 and Eq. 95, we therefore have

$\displaystyle|\omega(h)-\bar{\omega}|_{2}\leq Ch^{2p+2},$ (96)

where the constant $C$ depends only on $g$ and $p$.

We now prove the accuracy of the corrected trapezoidal rule Eq. 34. We denote
the Taylor polynomial of $\phi$ with degree $2p+1$ at the origin by

$\displaystyle
P_{\phi}(x)=\sum_{\begin{subarray}{c}\xi\in\mathbb{N}_{0}^{2},\\\
0\leq|\xi|_{1}\leq
2p+1\end{subarray}}\frac{x^{\xi}}{\xi!}\frac{\partial^{\xi}\phi(0)}{\partial
x^{\xi}},$ (97)

and define $\tilde{\phi}=P_{\phi}\cdot g$. It is clear that

$\displaystyle\partial^{k}[\phi(x)-\tilde{\phi}(x)]|_{x=0}=0,\quad\forall\
k\in\mathbb{N}_{0}^{2}:|k|_{1}=0,1,\dots,2p+1.$ (98)

Let $f\coloneqq\phi\cdot s$ and
$I[f]\coloneqq\int_{\mathbb{R}^{2}}f(x)\operatorname{d}\\!{x}$, we split
$Q_{h}^{p}[f]-I[f]$ as

$\displaystyle Q_{h}^{p}[f]-I[f]$
$\displaystyle=A_{h}^{p}[\phi-\tilde{\phi}]+(T_{h}^{0}-I)[(\phi-\tilde{\phi})\cdot
s]+(Q_{h}^{p}-I)[\tilde{\phi}\cdot s]$ $\displaystyle\eqqcolon
E_{1}+E_{2}+E_{3}.$ (99)

To estimate $E_{1}$, from Eq. 98 and Taylor’s theorem we have

$\displaystyle|\phi(x)-\tilde{\phi}(x)|\leq\max_{\begin{subarray}{c}x\in\mathbb{R}^{2},\\\
k\in\mathbb{N}_{0}^{2},\\\
|k|_{1}=2p+2\end{subarray}}\left(\frac{|\partial^{k}(\phi-\tilde{\phi})(x)|}{k!}\right)\
|x|^{2p+2}.$

Then from Eq. 18 we get

$\displaystyle|E_{1}|=|A_{h}^{p}[\phi-\tilde{\phi}]|\leq
Ch^{2-\alpha}\max_{1\leq j\leq
N_{p}}|\bar{\omega}_{j}|\max_{\xi\in\mathcal{G}_{j}:1\leq j\leq
N_{p}}|\phi(\xi h)-\tilde{\phi}(\xi h)|\leq Ch^{2p+4-\alpha}.$ (100)

In order to bound $E_{2}$, we apply Corollary 3.6 to $\phi-\tilde{\phi}$,
giving rise to

$\displaystyle|(T_{h}^{0}-I)[(\phi-\tilde{\phi})\cdot s]|\leq
Ch^{2p+4-\alpha}.$ (101)

Finally, it remains to estimate $E_{3}$. Since

$\displaystyle|E_{3}|$ $\displaystyle=|Q^{p}_{h}[\tilde{\phi}\cdot
s]-I[\tilde{\phi}\cdot s]|$
$\displaystyle\leq\sum_{\begin{subarray}{c}\xi\in\mathbb{N}_{0}^{2},\\\
0\leq|\xi|_{1}\leq
2p+1\end{subarray}}\frac{1}{\xi!}\left|\frac{\partial^{\xi}\phi(0)}{\partial
x^{\xi}}\right||Q_{h}^{p}[g\cdot s\cdot x^{\xi}]-I[g\cdot s\cdot x^{\xi}]|$
$\displaystyle\leq C\sum_{\begin{subarray}{c}\xi\in\mathbb{N}_{0}^{2},\\\
0\leq|\xi|_{1}\leq 2p+1\end{subarray}}|T_{h}^{0}[g\cdot s\cdot
x^{\xi}]-I[g\cdot s\cdot x^{\xi}]+A^{p}_{h}[g\cdot x^{\xi}]|.$ (102)

We claim that $T_{h}^{0}[g\cdot s\cdot x^{\xi}]$, $I[g\cdot s\cdot x^{\xi}]$
and $A^{p}_{h}[g\cdot x^{\xi}]$ vanish when at least one of $\xi_{1}$ or
$\xi_{2}$ is odd. To see this, by radial symmetry of $g$ and
$\frac{1}{|x|^{2+\alpha}}$,

$\displaystyle I[g\cdot s\cdot
x^{\xi}]=\int_{\mathbb{R}^{2}}\frac{g(x)}{|x|^{2+\alpha}}x^{\xi+2e_{1}}\operatorname{d}\\!{x}$
(103)
$\displaystyle=\int_{[0,\infty)^{2}}\frac{g(x)}{|x|^{2+\alpha}}[x^{\xi+2e_{1}}+(-x_{1})^{\xi_{1}+2}(x_{2})^{\xi_{2}}+(x_{1})^{\xi_{1}+2}(-x_{2})^{\xi_{2}}+(-x_{1})^{\xi_{1}+2}(-x_{2})^{\xi_{2}}]\operatorname{d}\\!{x}$
$\displaystyle=\int_{[0,\infty)^{2}}g(x)s(x)x^{\xi}[1+(-1)^{\xi_{1}}+(-1)^{\xi_{2}}+(-1)^{|\xi|_{1}}]\operatorname{d}\\!{x}.$

It is easy to see that Eq. 103 vanish if at least one of $\xi_{1}$ or
$\xi_{2}$ is odd. The argument for $T_{h}^{0}[g\cdot s\cdot x^{\xi}]=0$ under
same condition for $\xi$ closely resembles above. Recall that $g(\xi h)$ is
constant over each $\mathcal{G}_{q-m}^{m}$ and if
$(a_{1},a_{2})\in\mathcal{G}_{q-m}^{m}$ then $(\pm a_{1},\pm a_{2}),(\mp
a_{1},\pm a_{2})\in\mathcal{G}_{q-m}^{m}$ for all $(q-m,m)\in\mathcal{I}_{p}$.
Denote $g_{q-m,m}=g(\xi h),\ \forall\xi\in\mathcal{G}_{q-m}^{m}$. By Eq. 18 we
have

$\displaystyle A^{p}_{h}[g\cdot
x^{\xi}]=a(h)\sum_{(q-m,m)\in\mathcal{I}_{p}}\omega_{q-m}^{m}g_{q-m,m}\sum_{\beta\in\mathcal{G}_{q-m}^{m}}(\beta
h)^{\xi}$ (104)
$\displaystyle=a(h)\sum_{(q-m,m)\in\mathcal{I}_{p}}\omega_{q-m}^{m}g_{q-m,m}h^{|\xi|_{1}}\sum_{\beta\in\mathcal{G}_{q-m}^{m}}\beta^{\xi}\frac{1+(-1)^{\xi_{1}}+(-1)^{\xi_{2}}+(-1)^{|\xi|_{1}}}{4}$
$\displaystyle=0,$

if at least one of $\xi_{1}$ or $\xi_{2}$ is odd. Therefore, we can re-write
Section 3.1 as

$\displaystyle|E_{3}|$ $\displaystyle\leq
C\sum_{\begin{subarray}{c}0\leq|\xi|_{1}\leq p,\\\
\xi\in\mathbb{N}_{0}^{2}\end{subarray}}|T_{h}^{0}[g\cdot s\cdot
x^{2\xi}]-I[g\cdot s\cdot x^{2\xi}]+A^{p}_{h}[g\cdot x^{2\xi}]|.$ (105)

By Eqs. 16, LABEL: and 96 and $a(h)=h^{2-\alpha}$, $E_{3}$ can be estimated by

$\displaystyle|E_{3}|$ $\displaystyle\leq
Ca(h)\sum_{(q-m,m)\in\mathcal{I}_{p}}|\omega^{m}_{q-m}(h)-\bar{\omega}_{q-m}^{m}|\left(\sum_{\beta\in\mathcal{G}_{q-m}^{m}}|g(\beta
h)||\beta h|^{2|\xi|_{1}}\right)$ (106) $\displaystyle\leq
Ch^{2-\alpha}|\omega(h)-\bar{\omega}|_{2}\leq Ch^{2p+4-\alpha}.$

Putting together Eqs. 100, LABEL:, 101, LABEL: and 106, we have proved

$\displaystyle|Q_{h}^{p}[f]-I[f]|\leq Ch^{2p+4-\alpha}.$ (107)

∎

The following corollary provides a Taylor’s expansion of $c(h)$ as in Eq. 22,
it is used in Richardson extrapolation for numerical computation of correction
weights. Recall that $\mathfrak{C}(\alpha)=\lim_{h\to 0}c(h)$ in Eq. 89.

######  Corollary 3.7.

Let $p,k\in\mathbb{N}_{0}$. For any $g\in C_{c}^{\infty}(\mathbb{R}^{2})$ such
that $g(0)=1$, $\partial^{l}g(0)=0,\ l\in\mathbb{N}_{0}^{2}$ with
$0\not=|l|_{1}\leq 2p+1$, we have

$\displaystyle c(h)=\mathfrak{C}(\alpha)+h^{2p+2}\sum_{j=0}^{k}C_{j}h^{2j},$
(108)

where $C_{j}$’s are constants independent of $h$.

###### Proof.

If $k=0$, this is proved in the Theorem 3.1. Assume corollary holds for some
$k\in\mathbb{N}$. We prove the corollary holds for $k+1$. Fix
$\xi_{i}\in\mathcal{I}_{p},\ i\in\\{1,\dots,N_{p}\\}$. It is suffices to show
$c_{i}(h)=C(\xi_{i},\alpha)+h^{2p+2}O(\sum_{j=0}^{k+1}h^{2j})$. To this end,
let $g^{*}\in C_{c}^{\infty}(\mathbb{R}^{2})$ such that $g^{*}(0)=2$ and
$\partial^{l}g^{*}(0)=0,\ l\in\mathbb{N}_{0}^{2}:0\not=|l|_{1}\leq 2(p+k)+3$,
then by Theorem 3.4 we have

$\displaystyle\left|\frac{1}{h^{2|\xi_{i}|_{1}+2-\alpha}}\left[\int_{\mathbb{R}^{2}}g^{*}(x)s(x)x^{2\xi_{i}}\operatorname{d}\\!{x}-T_{h}^{0}[g^{*}\cdot
s\cdot x^{2\xi_{i}}]\right]-2C(\xi_{i},\alpha)\right|\lesssim h^{2p+2+2k+2}.$
(109)

Let $\bar{g}=g^{*}-g$, then $\bar{g}(0)=1,\ \partial^{l}\bar{g}(0)=0$ for all
$l\in\mathbb{N}^{2}_{0}:0\not=|l|_{1}\leq 2p+1$. Applying inductive hypothesis
to $\bar{g}$ we have

$\displaystyle\left|\frac{1}{h^{2|\xi_{i}|_{1}+2-\alpha}}\left[\int_{\mathbb{R}^{2}}\bar{g}(x)s(x)x^{2\xi_{i}}\operatorname{d}\\!{x}-T_{h}^{0}[\bar{g}\cdot
s\cdot x^{2\xi_{i}}]\right]-C(\xi_{i},\alpha)\right|\leq
h^{2p+2}O\left(\sum_{j=0}^{k}h^{2j}\right).$ (110)

Summing Eq. 109 with Eq. 110,

$\displaystyle|c_{i}(h)-C(\xi_{i},\alpha)|$
$\displaystyle=\left|\frac{1}{h^{2|\xi_{i}|_{1}+2-\alpha}}\left[\int_{\mathbb{R}^{2}}g(x)s(x)x^{2\xi_{i}}\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot
s\cdot x^{2\xi_{i}}]\right]-C(\xi_{i},\alpha)\right|$ (111) $\displaystyle\leq
h^{2p+2}O\left(\sum_{j=0}^{k+1}h^{2j}\right).$

This completes the induction. ∎

###  3.2 Off-diagonal Analysis

In this subsection we prove order of convergence of corrected trapezoidal rule
for off-diagonal singular integral Eq. 28 is $2p+2-\alpha,\ p\geq 1$. Note
that in off-diagonal case, it turns out that the correction weights for
$\xi=(0,0),(0,1)$ are zero and hence the correction for punctured-hole
trapezoidal rule is automatic at these grid points. As a consequence, the
number of correction layers is $p-1$, compared with $p+1$ in the on-diagonal
case. We omit the proofs of Lemmas 3.8, LABEL:, 3.9, LABEL: and 3.10, as they
are almost identical with that of Lemmas 3.3, LABEL:, 3.4, LABEL: and 3.6,
respectively.

######  Lemma 3.8.

For all $k,\eta\in\mathbb{N}_{0}^{n}$, we have

$\displaystyle\partial^{k}\left\\{\frac{x^{\eta+e_{1}+e_{2}}}{|x|^{n+\alpha}}\right\\}=\frac{1}{|x|^{n+\alpha-2-|\eta|_{1}+|k|_{1}}}P_{k,\eta}\left(\frac{x}{|x|}\right),$
(112)

where $P_{k,\eta}$ is a polynomial with
$\deg{P_{k,\eta}}=2+|k|_{1}+|\eta|_{1}$.

######  Theorem 3.9.

Let $n\in\mathbb{N}$, $p\in\mathbb{N}$ and $M=\max(4p,n+1)$. Assume $g\in
C_{c}^{N}(\mathbb{R}^{n}):N\geq M$ with $\partial^{k}g(0)=0$ for all
$|k|_{1}=1,2,\dots,2p-1$. For any $\xi\in\mathbb{N}_{0}^{n}:|\xi|_{1}\leq
2p-2$, we have

$\displaystyle\Big{|}\int_{\mathbb{R}^{n}}g\cdot s_{12}\cdot
x^{\xi}\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot s_{12}\cdot
x^{\xi}]-g(0)h^{|\xi|_{1}+2-\alpha}c(\xi,\alpha)\Big{|}\leq
Ch^{2p+2-\alpha+|\xi|_{1}},$ (113)

where $c(\xi,\alpha)\in\mathbb{R}$ is a constant.

######  Remark.

By replacing $M,\xi$ in Theorem 3.4 with the ones given in Theorem 3.9 and
replacing Lemma 3.3 with Lemma 3.8, the proof of Theorem 3.9 is verbatim proof
of Theorem 3.4.

######  Corollary 3.10.

Let $n\in\mathbb{N}$, $p\in\mathbb{N}$. Assume $g\in
C_{c}^{N}(\mathbb{R}^{n})$ with $N\geq\max(2p+2,n+1)$ such that
$\partial^{k}g(0)=0$ for all $|k|_{1}=0,1,\dots,2p-1$.Then

$\displaystyle\Big{|}\int_{\mathbb{R}^{n}}g\cdot
s_{12}\operatorname{d}\\!{x}-T_{h}^{0}[g\cdot s_{12}]\Big{|}\leq
Ch^{2p+2-\alpha}.$ (114)

We now state and prove the main theorem for the off-diagonal case.

######  Theorem 3.11.

Let $p\in\mathbb{N}$. Given $\phi\in C_{c}^{N}(\mathbb{R}^{2})$ with $N\geq
2p+2$ and the singular function $s_{12}(x)=\frac{x_{1}x_{2}}{|x|^{2+\alpha}},\
x\in\mathbb{R}^{2}$, $0<\alpha<2$. Denote $f\coloneqq\phi\cdot s_{12}$ and let
$g\in C_{c}^{\infty}(\mathbb{R}^{2})$ such that $g$ is radially symmetric,
$g(0)=1$ and $\partial^{k}g(0)=0$ for all $|k|_{1}=1,\dots,2p-1$. Let
$\omega(h)$ be the solution of Eq. 31 for this $g$, then $\omega(h)$ converges
to some $\bar{\omega}\in\mathbb{R}^{N_{p}}$ and the convergence rate is given
by

$\displaystyle|\omega(h)-\bar{\omega}|_{2}=O(h^{2p}).$ (115)

Moreover, there exists $C>0$ such that

$\displaystyle\left|Q_{h}^{p}[f]-\int_{\mathbb{R}^{2}}f(x)\ dx\right|\leq
Ch^{2p+2-\alpha}.$ (116)

###### Proof.

We first show the limit of the solution $\omega(h)$ of the linear system Eq.
31 exists. To this end, we apply Theorem 3.9 to each $c_{i}(h)$ in Eq. 32,
yielding

$\displaystyle c_{i}(h)$
$\displaystyle=\frac{1}{h^{2|\xi_{i}|_{1}-\alpha}}\left(\int
g(x)s(x)x^{2\xi_{i}-e_{1}-e_{2}}\ dx-T_{h}^{0}[g\cdot s\cdot
x^{2\xi_{i}-e_{1}-e_{2}}]\right)$
$\displaystyle=\frac{1}{h^{2|\xi_{i}|_{1}-\alpha}}\left(g(0)h^{2|\xi_{i}|_{1}-\alpha}c(\xi_{i},\alpha)+O(h^{2p-\alpha+2|\xi_{i}|_{1}})\right)$
$\displaystyle=c(\xi_{i},\alpha)+O(h^{2p}).$ (117)

Note that we write $c(2\xi_{i}-e_{1}-e_{2},\alpha)$ as $c(\xi_{i},\alpha)$ for
convenience. Define

$\displaystyle\mathfrak{C}(\alpha)$
$\displaystyle\coloneqq(c(\xi_{1},\alpha),\dots,c(\xi_{N_{p}},\alpha)),$ (118)

and let $\bar{\omega}$ be the solution of the linear system

$\displaystyle\boldsymbol{K}\bar{\omega}=\mathfrak{C}(\alpha).$ (119)

Assuming $\boldsymbol{K}$ is non-singular and using the arguments from Eqs. 89
to 93, we obtain

$\displaystyle|\omega(h)-\bar{\omega}|_{2}\leq\|\boldsymbol{G}(h)^{-1}\|_{2}\
\|\boldsymbol{I}-\boldsymbol{G}(h)\|_{2}\ \|\boldsymbol{K}^{-1}\|_{2}\
|c(h)|_{2}+\|\boldsymbol{K}^{-1}\|_{2}\ |c(h)-\mathfrak{C}(\alpha)|_{2}.$
(120)

A Taylor expansion of $g$ around the origin gives

$\displaystyle|1-g(\xi_{k}h)|\leq Ch^{2p},\quad k=1,\dots,N_{p},$ (121)

and hence

$\displaystyle\|\boldsymbol{I}-\boldsymbol{G}(h)\|_{2}=\max\\{|1-g(\xi_{k}h)|:k=1,\dots,N_{p}\\}\lesssim
h^{2p}.$ (122)

From Eq. 120, Section 3.2 and Eq. 121 we therefore have

$\displaystyle|\omega(h)-\bar{\omega}|_{2}\leq Ch^{2p}.$ (123)

We now prove the the accuracy of the corrected trapezoidal rule Eq. 28. We
denote the Taylor’s polynomial of $\phi$ at the origin by

$\displaystyle
P_{\phi}(x)=\sum_{\begin{subarray}{c}\xi\in\mathbb{N}_{0}^{2},\\\
0\leq|\xi|_{1}\leq
2p-1\end{subarray}}\frac{x^{\xi}}{\xi!}\frac{\partial^{\xi}\phi(0)}{\partial
x^{\xi}},$ (124)

and define $\tilde{\phi}=P_{\phi}\cdot g$. It is clear that

$\displaystyle\partial^{k}[\phi(x)-\tilde{\phi}(x)]|_{x=0}=0,\quad
k\in\mathbb{N}_{0}^{2}:|k|_{1}=0,1,\dots,2p-1.$ (125)

Writing $I[f]\coloneqq\int_{\mathbb{R}^{2}}f(x)dx$, we can split
$Q_{h}^{p}[f]-I[f]$ as

$\displaystyle Q_{h}^{p}[f]-I[f]$
$\displaystyle=A_{h}^{p}[\phi-\tilde{\phi}]+(T_{h}^{0}-I)[(\phi-\tilde{\phi})\cdot
s_{12}]+(Q_{h}^{p}-I)[\tilde{\phi}\cdot s_{12}]$ $\displaystyle\eqqcolon
E_{1}+E_{2}+E_{3}.$

To estimate $E_{1}$, from Eq. 125 and Taylor’s theorem we have

$\displaystyle|\phi(x)-\tilde{\phi}(x)|\leq\max_{\begin{subarray}{c}x\in\mathbb{R}^{2},\\\
k\in\mathbb{N}_{0}^{2},\\\
|k|_{1}=2p\end{subarray}}\left(\frac{|\partial^{k}(\phi-\tilde{\phi})(x)|}{k!}\right)\
|x|^{2p},$

From definition of $A_{h}^{p}$ in Eq. 28 we obtain

$\displaystyle|E_{1}|=|A_{h}^{p}[\phi-\tilde{\phi}]|\leq
h^{2-\alpha}\max_{1\leq j\leq
N_{p}}|\bar{\omega}_{j}|\max_{\xi\in\mathcal{G}_{j}:1\leq j\leq
N_{p}}|\phi(\xi h)-\tilde{\phi}(\xi h)|\leq h^{2p+2-\alpha}.$ (126)

To bound $E_{2}$, we apply Corollary 3.10 to $\phi-\tilde{\phi}$, giving rise
to

$\displaystyle|E_{2}|=|(T_{h}^{0}-I)[(\phi-\tilde{\phi})\cdot s_{12}]|\leq
Ch^{2p+2-\alpha}.$ (127)

Finally, it remains to estimate $E_{3}$. Since

$\displaystyle|E_{3}|$ $\displaystyle=|Q^{p}_{h}[\tilde{\phi}\cdot
s_{12}]-I[\tilde{\phi}\cdot s_{12}]|$
$\displaystyle\leq\sum_{\begin{subarray}{c}\xi\in\mathbb{N}_{0}^{2},\\\
0\leq|\xi|_{1}\leq
2p-1\end{subarray}}\frac{1}{\xi!}\left|\frac{\partial^{\xi}\phi(0)}{\partial
x^{\xi}}\right||Q_{h}^{p}[g\cdot s_{12}\cdot x^{\xi}]-I[g\cdot s_{12}\cdot
x^{\xi}]|$
$\displaystyle\leq\sum_{\begin{subarray}{c}\xi\in\mathbb{N}_{0}^{2},\\\
0\leq|\xi|_{1}\leq 2p-1\end{subarray}}|T_{h}^{0}[g\cdot s_{12}\cdot
x^{\xi}]-I[g\cdot s_{12}\cdot x^{\xi}]+A^{p}_{h}[g\cdot x^{\xi}]|.$ (128)

We claim $T_{h}^{0}[g\cdot s_{12}\cdot x^{\xi}]$, $I[g\cdot s_{12}\cdot
x^{\xi}]$ and $A^{p}_{h}[g\cdot x^{\xi}]$ vanish when at least one of
$\xi_{1}$ or $\xi_{2}$ is even. To see this, by radial symmetry of $g$ and
$\frac{1}{|x|^{2+\alpha}}$, we have

$\displaystyle I[g\cdot s_{12}\cdot
x^{\xi}]=\int_{\mathbb{R}^{2}}\frac{g(x)}{|x|^{2+\alpha}}x^{\xi+e_{1}+e_{2}}\
\operatorname{d}\\!{x}$
$\displaystyle=\int_{[0,\infty)^{2}}\frac{g(x)}{|x|^{2+\alpha}}[x^{\xi+e_{1}+e_{2}}+(-x_{1})^{\xi_{1}+1}(x_{2})^{\xi_{2}+1}+(x_{1})^{\xi_{1}+1}(-x_{2})^{\xi_{2}+1}+(-x_{1})^{\xi_{1}+1}(-x_{2})^{\xi_{2}+1}]\
\operatorname{d}\\!{x}$
$\displaystyle=\int_{[0,\infty)^{2}}g(x)s_{12}(x)x^{\xi}[1-(-1)^{\xi_{1}}-(-1)^{\xi_{2}}+(-1)^{|\xi|_{1}}]\
\operatorname{d}\\!{x}.$ (129)

It is easy to see that Section 3.2 vanish when at least one of $\xi_{1}$ or
$\xi_{2}$ is even. Argument for $T_{h}^{0}[g\cdot s_{12}\cdot x^{\xi}]=0$
under same condition for $\xi$ closely resembles above. Recall that $g(\xi h)$
is constant over each $\mathcal{G}_{q-m,m}$ and if
$(a_{1},a_{2})\in\mathcal{G}_{q-m,m}$ then $(\pm a_{1},\pm a_{2}),(\mp
a_{1},\pm a_{2}),(\pm a_{2},\pm a_{1}),(\mp a_{2},\pm
a_{1})\in\mathcal{G}_{q-m,m}$ for all $(q-m,m)\in\mathcal{I}_{p}$. Denote
$g_{q-m,m}=g(\xi h),\ \text{for any }\xi\in\mathcal{G}_{q-m,m}$, we have

$\displaystyle A^{p}_{h}[g\cdot x^{\xi}]$
$\displaystyle=h^{2-\alpha}\sum_{(q-m,m)\in\mathcal{I}_{p}}\bar{\omega}_{q-m,m}g_{q-m,m}\sum_{\beta\in\mathcal{G}_{q-m,m}}\text{sgn}(\beta_{1}\beta_{2})(\beta
h)^{\xi},$
$\displaystyle\sum_{\beta\in\mathcal{G}_{q-m,m}}\text{sgn}(\beta_{1}\beta_{2})(\beta
h)^{\xi}$
$\displaystyle=(q-m)^{\xi_{1}}m^{\xi_{2}}\Big{(}1-(-1)^{\xi_{1}}-(-1)^{\xi_{2}}+(-1)^{|\xi|_{1}}\Big{)}$
$\displaystyle+(1-\delta_{q-m,m})m^{\xi_{1}}(q-m)^{\xi_{2}}\Big{(}1-(-1)^{\xi_{1}}-(-1)^{\xi_{2}}+(-1)^{|\xi|_{1}}\Big{)}$
$\displaystyle=0,$ (130)

if at least one of $\xi_{1}$ or $\xi_{2}$ is even.

By another round of symmetry argument, it is also clear that for all
$\xi\in\mathbb{N}_{0}^{2}$,

$\displaystyle I(g\cdot s_{12}\cdot x^{(\xi_{1},\xi_{2})})$
$\displaystyle=I(g\cdot s_{12}\cdot x^{(\xi_{2},\xi_{1})}),$ (131)
$\displaystyle T_{h}^{0}(g\cdot s_{12}\cdot x^{(\xi_{1},\xi_{2})})$
$\displaystyle=T_{h}^{0}(g\cdot s_{12}\cdot x^{(\xi_{2},\xi_{1})}),$ (132)
$\displaystyle A^{p}_{h}(g\cdot x^{(\xi_{1},\xi_{2})})$
$\displaystyle=A^{p}_{h}(g\cdot x^{(\xi_{2},\xi_{1})}).$ (133)

When $\xi=0$, using the odd symmetry of $s_{12}$ and property of $A_{h}^{p}$,
we have

$\displaystyle I(g\cdot s_{12})$ $\displaystyle=T_{h}^{0}(g\cdot s_{12})=0,$
(134) $\displaystyle A_{h}^{p}(g)$
$\displaystyle=h^{2-\alpha}\sum_{(q-m,m)\in\mathcal{I}_{p}}\omega_{q-m,m}g_{q-m,m}\underbrace{\sum_{\beta\in\mathcal{G}_{q-m,m}}\text{sgn}(\beta_{1}\beta_{2})}_{=0}=0.$
(135)

Therefore, we can now re-write Section 3.2 as

$\displaystyle|E_{3}|$
$\displaystyle\leq\sum_{|\xi|_{1}=0,2,\dots,2p-2}|T_{h}^{0}[g\cdot s_{12}\cdot
x^{\xi}]-I[g\cdot s_{12}\cdot x^{\xi}]+A^{p}_{h}[g\cdot x^{\xi}]|$
$\displaystyle=\sum_{\xi\in\mathcal{I}_{p}}|T_{h}^{0}[g\cdot s_{12}\cdot
x^{2\xi-e_{1}-e_{2}}]-I[g\cdot s_{12}\cdot x^{2\xi-
e_{1}-e_{2}}]+A^{p}_{h}[g\cdot x^{2\xi-e_{1}-e_{2}}]|.$ (136)

By Eq. 31 and Eq. 123, $E_{3}$ can be estimated by

$\displaystyle|E_{3}|$ $\displaystyle\leq
Ch^{2-\alpha}\sum_{\xi\in\mathcal{I}_{p}}\sum_{(q-m,m)\in\mathcal{I}_{p}}|\omega_{q-m,m}(h)-\bar{\omega}_{q-m,m}|\left(\sum_{\beta\in\mathcal{G}_{q-m}^{m}}|g(\beta
h)||\beta h|^{2|\xi|_{1}-e_{1}-e_{2}}\right)$ $\displaystyle\leq
Ch^{2-\alpha}|\omega(h)-\bar{\omega}|_{2}\leq Ch^{2p+2-\alpha}.$ (137)

Putting together Eq. 126, Eq. 127 and Section 3.2, we proved

$\displaystyle|Q_{h}^{p}[f]-I[f]|\leq Ch^{2p+2-\alpha}.$ (138)

∎

The following corollary provides a Taylor’s expansion of $c(h)$ as in Eq. 31.
It is used in Richardson extrapolation for numerical computation of correction
weights. Recall that $\mathfrak{C}(\alpha)=\lim_{h\to 0}c(h)$ in Eq. 118. We
skip the proof that runs almost identically with that of Corollary 3.7.

######  Corollary 3.12.

Let $p,k\in\mathbb{N}_{0}:p\geq 1$. For any $g\in
C_{c}^{\infty}(\mathbb{R}^{2})$ such that $g(0)=1$, $\partial^{l}g(0)=0,\
l\in\mathbb{N}_{0}^{2}$ with $0\not=|l|_{1}\leq 2p-1$, we have

$\displaystyle c(h)=\mathfrak{C}(\alpha)+h^{2p}\sum_{j=0}^{k}C_{j}h^{2j},$
(139)

where $C_{j}$’s are constants independent of $h$.

##  4 Non-singularity of Coefficient Matrices $\boldsymbol{K}$

###  4.1 The On-diagonal Case

In this subsection we prove the on-diagonal coefficient matrix
$\boldsymbol{K}$ defined in Eq. 19 is non-singular. We order all
$\xi\in\mathcal{I}_{p}$ as follows, $\xi\in\mathcal{I}_{p}$ with $\xi_{1}=0$
first, followed by all $\xi\in\mathcal{I}_{p}$ with $\xi_{2}=0$ and then
followed by the rest $\xi\in\mathcal{I}_{p}$ with $\xi_{1}\xi_{2}\not=0$,
where we denote the three subsets by
$\mathcal{C}_{1},\mathcal{C}_{2},\mathcal{C}_{3}$ respectively. Namely,
$\mathcal{C}_{1}$ is arranged by

$\displaystyle\mathcal{C}_{1}=\\{(0,r):0\leq r\leq
p\\}=\\{(0,0),(0,1),\dots,(0,p)\\},$ (140)

and similarly,

$\displaystyle\mathcal{C}_{2}$ $\displaystyle=\\{(s,0):1\leq s\leq
p\\}=\\{(1,0),(2,0),\dots,(p,0)\\}.$ (141)

The arrangement for $\mathcal{C}_{3}$ can in fact be arbitrary, thus we simply
choose a natural one as follows, pictorially, we can list all elements in
$\mathcal{C}_{3}$ as a triangle

$\displaystyle\begin{array}[]{cccccc}(1,1)&(1,2)&\dots&(1,k)&\cdots&(1,p-1)\\\
\downarrow&(k-1,2)\nearrow&\vdots&\nearrow\\\ (k,1)\nearrow&\dots&(k,p-k)&\\\
\vdots&\nearrow&&\\\ (p-1,1)&&&\end{array}\ $ (147)

We arrange all elements in $\mathcal{C}_{3}$ from $(1,1)$, followed by
diagonal $(2,1),(1,2)$, $\dots$, followed by diagonal $(k,1),\dots,(1,k)$,
$\dots$, i.e.

$\displaystyle\mathcal{C}_{3}=\\{(1,1),(2,1),(1,2),\dots,(p-1,1),(p-2,2),\dots,(1,p-1)\\}.$
(148)

We also arrange $\mathcal{G}$ in the same way. We first show matrix
$\boldsymbol{K}$ organized under this arrangement has a special structure.

######  Lemma 4.1.

Let $p\geq 1$, there exist an ordering on $\mathcal{I}_{p}$ such that

  1. 1.

$K_{1,1}=1$ while $K_{j,1}=0$ for all $2\leq j\leq N_{p}$.

  2. 2.

$\boldsymbol{K}$ can be expressed as a block matrix of the form

$\displaystyle\boldsymbol{K}=\begin{pmatrix}\boldsymbol{A}&\boldsymbol{B}\\\
\boldsymbol{C}&\boldsymbol{D}\end{pmatrix},$ (149)

such that $\boldsymbol{A}\in\mathbb{R}^{(2p+1)\times(2p+1)}$,
$\boldsymbol{C}=\boldsymbol{0}\in\mathbb{R}^{\frac{p(p-1)}{2}\times(2p+1)}$,
$\boldsymbol{D}\in\mathbb{R}^{\frac{p(p-1)}{2}\times\frac{p(p-1)}{2}}$. As a
result,

$\displaystyle\det\boldsymbol{K}=\det\boldsymbol{A}\cdot\det\boldsymbol{D}.$
(150)

Moreover, the block $\boldsymbol{A}$ is non-singular.

###### Proof.

We use the ordering on $\mathcal{I}_{p}$ introduced in the beginning of this
subsection. It is easy to verify

$\displaystyle K_{1,1}=0^{2\cdot 0}\cdot 0^{2\cdot 0}=1,$

and for all $\xi_{j}\in\mathcal{I}_{p}\setminus\\{0\\}$, either
$\xi_{j,1}\not=0$ or $\xi_{j,2}\not=0$, then

$\displaystyle K_{j,1}=0^{2\xi_{j,1}}\cdot 0^{2\xi_{j,2}}=0,\quad\forall\
2\leq j\leq N_{p}.$ (151)

The first part of the lemma is proved.

It is obvious that $|\mathcal{C}_{1}|=p+1$, $|\mathcal{C}_{2}|=p$ and
$|\mathcal{C}_{3}|=N_{p}-2p-1=\frac{p(p-1)}{2}$. $\boldsymbol{A}$ is generated
by all $\xi\in\mathcal{C}_{1}\cup\mathcal{C}_{2}$ and
$\mathcal{G}_{\gamma}:\gamma\in\mathcal{C}_{1}\cup\mathcal{C}_{2}$, thus it
has size $(2p+1)\times(2p+1)$. Similarly, $\boldsymbol{D}$ is generated by all
$\xi\in\mathcal{C}_{3}$ and $\mathcal{G}_{\gamma}:\gamma\in\mathcal{C}_{3}$,
hence its dimension is $\frac{p(p-1)}{2}\times\frac{p(p-1)}{2}$. To show
$\boldsymbol{C}$ is a _zero_ matrix, we note that all elements in
$\boldsymbol{C}$ are of the form
$\sum_{\beta\in\mathcal{G}_{\gamma}}\beta^{2\xi},\quad\forall\xi\in\mathcal{C}_{3}$,
$\forall\gamma\in\mathcal{C}_{1}\cup\mathcal{C}_{2}$. If
$\gamma\in\mathcal{C}_{1}$, then $\gamma=(0,r)$ for some $0\leq r\leq p$.
Since $\xi_{1}$ and $\xi_{2}$ are both non-zeros, we have
$\sum_{\beta\in\mathcal{G}_{\gamma}}\beta^{2\xi}=2r^{2\xi_{1}}\cdot
0^{2\xi_{2}}=0$. The case for $\gamma\in\mathcal{C}_{2}$ also holds by
identical argument.

We are left to show $\boldsymbol{A}$ is non-singular. We claim that
$\boldsymbol{A}$ has the form

$\displaystyle\boldsymbol{A}=\begin{pmatrix}1&*&*\\\ 0_{p\times
1}&2\boldsymbol{A}_{1}&\boldsymbol{A}_{2}=\boldsymbol{0}_{p\times p}\\\
0_{p\times 1}&\boldsymbol{A}_{3}=\boldsymbol{0}_{p\times
p}&2\boldsymbol{A}_{4}\end{pmatrix},$ (152)

where both $\boldsymbol{A}_{1},\boldsymbol{A}_{4}\in\mathbb{R}^{p\times p}$
are Vandermonde matrices. To see this, all elements in $2\boldsymbol{A}_{1}$
can be expressed by

$\displaystyle\sum_{\beta\in\mathcal{G}_{\gamma}}\beta^{2\xi},\quad\text{where}\
\xi,\gamma\in\mathcal{C}_{1}\setminus\\{(0,0)\\},$ (153)

let $\xi=(0,r)$ and $\gamma=(0,s)$ for some $1\leq r,s\leq p$, so

$\displaystyle
2\boldsymbol{A}_{1}(r,s)=\boldsymbol{A}(r+1,s+1)=\sum_{\beta\in\mathcal{G}_{\gamma}}\beta^{2\xi}=2(s^{2})^{r},$
(154)

which shows $\boldsymbol{A}_{1}$ is a Vandermonde matrix. Repeating the above
argument, elements in $2\boldsymbol{A}_{4}$ can be expressed by

$\displaystyle\sum_{\beta\in\mathcal{G}_{\gamma}}\beta^{2\xi},\quad\text{where}\
\xi,\gamma\in\mathcal{C}_{2},$ (155)

let $\xi=(r,0)$ and $\gamma=(s,0)$ for some $1\leq r,s\leq p$, we have

$\displaystyle
2\boldsymbol{A}_{4}(r,s)=\boldsymbol{A}(p+1+r,p+1+s)=\sum_{\beta\in\mathcal{G}_{\gamma}}\beta^{2\xi}=2(s^{2})^{r},$
(156)

which show $\boldsymbol{A}_{4}$ is a Vandermonde matrix. Finally, let
$\xi=(0,r)\in\mathcal{C}_{1}$ and $\gamma=(s,0)\in\mathcal{C}_{2}$, we compute

$\displaystyle\sum_{\beta\in\mathcal{G}_{\gamma}}\beta^{2\xi}=0.$ (157)

Thus $\boldsymbol{A}_{3}$ is the $p\times p$ zeros matrix.
$\boldsymbol{A}_{2}=\boldsymbol{0}_{p\times p}$ can then be proved
analogously. The non-singularity of matrix $\boldsymbol{A}$ is now obvious. ∎

We now focus on showing the sub-matrix $\boldsymbol{D}$ in Eq. 149 is non-
singular. We use algebraic combinatorial method to show the determinant of
$\boldsymbol{D}$ is non-zero. To do this, let $x_{i}\in\mathbb{R}$, $1\leq
i\leq p-1$, we define matrix
$\boldsymbol{E}\in\mathbb{R}^{\frac{p(p-1)}{2}\times\frac{p(p-1)}{2}}$ as
follows: for all $1\leq i,j\leq\frac{p(p-1)}{2}$

$\displaystyle E_{i,j}=x_{\xi_{j,1}}^{\xi_{i,1}}\cdot
x_{\xi_{j,2}}^{\xi_{i,2}}=(x_{\xi_{j,1}},x_{\xi_{j,2}})^{(\xi_{i,1},\xi_{i,2})},$
(158)

where $\xi_{i}=(\xi_{i,1},\xi_{i,2})$ and $\xi_{j}=(\xi_{j,1},\xi_{j,2})$ are
the $i$-th and $j$-th elements in $\mathcal{C}_{3}$, respectively. We now
prove

######  Theorem 4.2.

For all $p\geq 2$,

$\displaystyle\det\boldsymbol{E}=C\prod_{1\leq j\leq
p-1}x_{j}^{2(p-j)}\prod_{1\leq i<j\leq p-1}(x_{j}-x_{i})^{2(p-j)},$ (159)

where $C$ is a non-zero constant.

###### Proof.

Denote $\det\boldsymbol{E}\eqqcolon F(x_{1},\dots,x_{p-1})$. We observe that
the multivariate polynomial $F$ is homogeneous, i.e., all terms in $F$ have
the same degree. To see this, apply Leibniz’s formula, we have

$\displaystyle F(x_{1},\dots,x_{p-1})=\sum_{\sigma\in
S_{\frac{p(p-1)}{2}}}\text{sgn}(\sigma)\prod_{1\leq
i\leq\frac{p(p-1)}{2}}E_{i,\sigma(i)},$ (160)

where $S_{m}$ is the permutation group of the set $\\{1,\dots,m\\}$,
$m\in\mathbb{N}$. By construction of $\boldsymbol{E}$, every elements in
$i$-th row of $\boldsymbol{E}$ share the same degree, namely
$\xi_{i,1}+\xi_{i,2}$, for any $1\leq i\leq\frac{p(p-1)}{2}$, so $\prod_{1\leq
i\leq\frac{p(p-1)}{2}}K_{i,\sigma(i)}$ has same degree for all $\sigma\in
S_{\frac{p(p-1)}{2}}$ and this implies homogeneity of $F$. We can then compute
the degree of $F$ by

$\displaystyle\deg F=\sum_{q=2}^{p}q(q-1)=\frac{p^{3}-p}{3}.$ (161)

Consider $j\in\mathbb{N}$ such that $1\leq j\leq p-1$, we claim that
$x_{j}^{2(p-j)}$ divides $F$. By construction of $\boldsymbol{E}$, there are
$(p-j)$’s $x_{j}$ appearing in $(x_{j},\cdot)$ and $(p-j)$’s $x_{j}$ appearing
in $(\cdot,x_{j})$, so there are $2(p-j)$’s $x_{j}$ in total. Hence the
algebraic multiplicity of $x_{j}$ is at least $2(p-j)$.

For arbitrary pair of integers $(i,j)$ such that $1\leq i<j\leq p-1$, if we
substitute $x_{j}$ with $x_{i}$, then there would be at least two identical
columns in $\boldsymbol{E}$ and hence $F$ will vanish. Furthermore, there are
$(p-j)$’s $x_{j}$ appearing in $(x_{j},\cdot)$ and $(p-j)$’s $x_{j}$ appearing
in $(\cdot,x_{j})$. As a result, there are $2(p-j)$ instances. Therefore,
$(x_{j}-x_{i})^{2(p-j)}$ divides $F$.

Now, by unique factorization property of multivariate polynomial ring, there
exists a multivariate polynomial $C(x_{1},\dots,x_{p-1})$ such that

$\displaystyle F(x_{1},\dots,x_{p-1})$
$\displaystyle=C(x_{1},\dots,x_{p-1})\prod_{1\leq j\leq
p-1}x_{j}^{2(p-j)}\prod_{1\leq i<j\leq p-1}(x_{j}-x_{i})^{2(p-j)},$ (162)
$\displaystyle\eqqcolon C(x_{1},\cdots,x_{p-1})H(x_{1},\dots,x_{p-1}).$ (163)

$C(x_{1},\dots,x_{p-1})$ will be a constant polynomial if

$\displaystyle\deg H=\deg F.$ (164)

We compute

$\displaystyle\deg
H=\sum_{j=1}^{p-1}2(p-j)+\sum_{j=2}^{p-1}\sum_{i=1}^{j-1}2(p-j)=p(p-1)+\frac{1}{3}p(p-1)(p-2)=\frac{p^{3}-p}{3}=\deg
F.$ (165)

This concludes the proof of Theorem 4.2. ∎

######  Corollary 4.3.

The coefficient matrix $\boldsymbol{K}$ as in Eq. 19 is non-singular for all
$p\geq 1$.

###### Proof.

By Lemma 4.1, it suffices to show $\boldsymbol{D}$ is non-singular. Note that
$\boldsymbol{D}=4\boldsymbol{E}$ where $\boldsymbol{E}$ defined in Eq. 158
with $x_{i}=i^{2}$ for all $1\leq i\leq p-1$ and all $x_{i}$’s are non-zeros
and mutually distinct. Hence by Theorem 4.2

$\displaystyle\det\boldsymbol{D}=2^{p(p-1)}\det\boldsymbol{E}\not=0.$ (166)

This completes the proof. ∎

###  4.2 The Off-diagonal Case

In this subsection, we will prove off-diagonal coefficient matrix
$\boldsymbol{K}$ defined at Eq. 30 is non-singular, therefore the linear
system Eq. 31 always admits unique solution.

Let $p\in\mathbb{N},\ p\geq 2$. We recall the notations for
$\mathcal{G}_{q-m,m}$ Eq. 26 and $\mathcal{I}_{p}$ Eq. 27 and define:

$\displaystyle\mathcal{G}\coloneqq\left\\{\mathcal{G}_{q-m,m}:2\leq q\leq p,\
1\leq m\leq\left\lfloor\frac{q}{2}\right\rfloor\right\\}.$ (167)

The value of $N_{p}\coloneqq|\mathcal{I}_{p}|=|\mathcal{G}|$ is given

$\displaystyle N_{p}=\left\\{\begin{array}[]{cc}\frac{p^{2}}{4}&,\text{if
}p\text{ is even}\\\ \frac{p^{2}-1}{4}&,\text{if }p\text{ is
odd}\end{array}\right..$ (170)

We arrange $\mathcal{I}_{p}$ first by the increasing order of magnitude of
1-norm and followed by, within each group of same 1-norm, the reverse
dictionary order:

$\displaystyle\mathcal{I}_{p}$
$\displaystyle=\left\\{(1,1),(2,1),(3,1),(2,2),\dots,(p-1,1),(p-2,2),\dots,\left(\left\lceil\frac{p}{2}\right\rceil,\left\lfloor\frac{p}{2}\right\rfloor\right)\right\\}\eqqcolon\\{\beta_{j}\\}_{1\leq
j\leq N_{p}}.$ (171)

We arrange $\mathcal{G}$ in the same fashion with $\mathcal{I}_{p}$ so that

$\displaystyle\mathcal{G}$
$\displaystyle=\left\\{\mathcal{G}_{1,1},\mathcal{G}_{2,1},\dots,\mathcal{G}_{p-1,1},\mathcal{G}_{p-2,2},\dots,\mathcal{G}_{\lceil\frac{p}{2}\rceil,\lfloor\frac{p}{2}\rfloor}\right\\}\eqqcolon\\{\mathcal{G}_{j}\\}_{1\leq
j\leq N_{p}}.$ (172)

The following lemma reveals the structure of the coefficient matrix
$\boldsymbol{K}$ given in Eq. 30.

######  Lemma 4.4.

For all $p\in\mathbb{N},\ p\geq 2$

$\displaystyle\boldsymbol{K}=\boldsymbol{E}\boldsymbol{H},$ (173)

where $\boldsymbol{H}\in\mathbb{R}^{N_{p}\times N_{p}}$ is a diagonal matrix
with

$\displaystyle\boldsymbol{H}_{jj}=\left\\{\begin{array}[]{ccc}\nicefrac{{2}}{{|\xi_{j,1}\xi_{j,2}|}},&\text{if}&|\xi_{j,1}|=|\xi_{j,2}|,\
\forall\xi_{j}=(\xi_{j,1},\xi_{j,2})\in\mathcal{G}_{j}\\\
\nicefrac{{4}}{{|\xi_{j,1}\xi_{j,2}|}},&\text{if}&\text{otherwise},\end{array}\right.$
(176)

and $\boldsymbol{E}$ is given by

$\displaystyle\boldsymbol{E}_{i,j}=|\xi_{j,1}|^{2\xi_{i,1}}|\xi_{j,2}|^{2\xi_{i,2}}+|\xi_{j,2}|^{2\xi_{i,1}}|\xi_{j,1}|^{2\xi_{i,2}},$
(177)

for all $1\leq i,j\leq N_{p}$. Here, $\xi_{i},\ \xi_{j}$ is $i$-th and $j$-th
element in $\mathcal{I}_{p}$, respectively.

###### Proof.

From the definition of coefficient matrix $\boldsymbol{K}$ in Eq. 30, for any
$1\leq i,j\leq N_{p}$ we have

$\displaystyle K_{i,j}$
$\displaystyle=\sum_{\xi_{j}\in\mathcal{G}_{j}}\text{sgn}(\xi_{j,1}\xi_{j,2})\xi_{j}^{2\xi_{i}-(1,1)}$
$\displaystyle=\left\\{\begin{array}[]{ccc}4\nicefrac{{(\xi_{j,1},\xi_{j,2})^{2\xi_{i}}}}{{|\xi_{j,1}\xi_{j,2}|}},&\text{if}&|\xi_{j,1}|=|\xi_{j,2}|,\\\
4\nicefrac{{[(\xi_{j,1},\xi_{j,2})^{2\xi_{i}}+(\xi_{j,2},\xi_{j,1})^{2\xi_{i}}]}}{{|\xi_{j,1}\xi_{j,2}|}},&\text{if}&\text{otherwise}.\end{array}\right.$
(180)

The lemma follows immediately. ∎

From Lemma 4.4 we know that
$\det{\boldsymbol{K}}=\det{\boldsymbol{E}}\det{\boldsymbol{H}}$ and it is
clear that $\boldsymbol{H}$ is non-singular, therefore non-singularity of
$\boldsymbol{K}$ follows once we show $\boldsymbol{E}$ is non-singular. We
prove this fact by a algebraic combinatorial argument. Let
$x_{i}\in\mathbb{R},1\leq i\leq p-1$, we define matrix
$\boldsymbol{T}(x_{1},\dots,x_{p-1})\in\mathbb{R}^{N_{p}\times N_{p}}$ by

$\displaystyle
T_{i,j}=x_{\xi_{j,1}}^{\xi_{i,1}}x_{\xi_{j,2}}^{\xi_{i,2}}+x_{\xi_{j,1}}^{\xi_{i,2}}x_{\xi_{j,2}}^{\xi_{i,1}}=(x_{\xi_{j,1}},x_{\xi_{j,2}})^{(\xi_{i,1},\xi_{i,2})}+(x_{\xi_{j,2}},x_{\xi_{j,1}})^{(\xi_{i,1},\xi_{i,2})},$
(181)

where $\xi_{i}=(\xi_{i,1},\xi_{i,2})$ and $\xi_{j}=(\xi_{j,1},\xi_{j,2})$ are
the $i$-th and $j$-th element in $\mathcal{I}_{p}$, respectively.

######  Theorem 4.5.

Let $p\in\mathbb{N},\ p\geq 2$

$\displaystyle\det\boldsymbol{T}=C\prod_{1\leq
j\leq\lfloor\frac{p}{2}\rfloor}x_{j}^{p+1-j}\prod_{\lfloor\frac{p}{2}\rfloor<j\leq
p-1}x_{j}^{p-j}\prod_{1\leq
i<j\leq\lfloor\frac{p}{2}\rfloor}(x_{j}-x_{i})^{p+1-j}\prod_{\begin{subarray}{c}1\leq
i<j,\\\ \lfloor\frac{p}{2}\rfloor<j\leq
p-1\end{subarray}}(x_{j}-x_{i})^{p-j},$ (182)

where $C$ is a non-zero constant.

###### Proof.

Let us prove when $p=2k$, $k\geq 1$, the proof for the $p$ odd is very
similar. From Eq. 170 we note that $N_{2k}=k^{2}$. Denote
$\det\boldsymbol{T}\eqqcolon F(x_{1},\dots,x_{2k-1})$. We observe that the
multivariate polynomial $F$ is homogeneous. To see this, apply Leibniz’s
formula, we have

$\displaystyle F(x_{1},\dots,x_{2k-1})=\sum_{\sigma\in
S_{k^{2}}}\text{sgn}(\sigma)\prod_{1\leq j\leq k^{2}}T_{j,\sigma(j)},$ (183)

where $S_{n}$ is the permutation group of the set $\\{1,\dots,n\\}$,
$n\in\mathbb{N}$. By the construction of $\boldsymbol{T}$, all elements in
$i$-th row of $\boldsymbol{T}$ share the degree $|\xi_{i}|_{1}$, where
$\xi_{i}\in\mathcal{I}_{p}$ for any $1\leq i\leq k^{2}$. Thus, $\prod_{1\leq
j\leq k^{2}}\boldsymbol{T}_{j,\sigma(j)}$ has same degree for all $\sigma\in
S_{k^{2}}$ and this implies homogeneity of $T$. We can then compute the degree
of $F$ by

$\displaystyle\deg F$
$\displaystyle=\sum_{q=2}^{2k}q\lfloor\nicefrac{{q}}{{2}}\rfloor=\frac{4}{3}k^{3}+\frac{1}{2}k^{2}+\frac{1}{6}k.$
(184)

Consider $j\in\mathbb{N}$ such that $1\leq j\leq 2k-1$, suppose we substitute
$0$ into $x_{j}$, then there would be at least one _zero_ -column in
$\boldsymbol{T}$ and hence $F$ would be zero. Therefore, each $x_{j}$ divides
$F$. We further prove that, if $1\leq j\leq k$, then $x_{j}^{2k+1-j}$ divides
$F$ and if $k<j\leq 2k-1$, then $x_{j}^{2k-j}$ divides $F$. By construction of
$\boldsymbol{T}$, for all $1\leq j\leq k$, there are $j$’s $x_{j}$ appearing
in $(x_{j},\cdot)$ and $2k-2j+1$’s $x_{j}$ appearing in $(\cdot,x_{j})$, thus
there are $2k-j+1$’s $x_{j}$ in total. On the other hand, for all $k<j\leq
2k-1$, there are $2k-j$’s $x_{j}$ appearing in $(x_{j},\cdot)$ and there is no
$(\cdot,x_{j})$, so there are $2k-j$’s $x_{j}$ in total. Hence the algebraic
multiplicity of $x_{j}$ is at least $2k+1-j$ for all $1\leq j\leq k$ and
$2k-j$ for all $k<j\leq 2k-1$.

For arbitrary pair of integers $(i,j)$ such that $1\leq i<j\leq 2k-1$, if we
substituted $x_{j}$ with $x_{i}$, then there would be at least two identical
columns in $\boldsymbol{T}$ and hence $F$ would vanish. Furthermore, for all
integer pair $(i,j)$ such that $1\leq i<j\leq k$, there are $j$’s $x_{j}$
appearing in $(x_{j},\cdot)$ and $2k-2j+1$’s $x_{j}$ appearing in
$(\cdot,x_{j})$, and for all integer pair $(i,j)$ such that $1\leq i<j,k<j\leq
2k-1$, there are $2k-j$’s $x_{j}$ appearing in $(x_{j},\cdot)$. As a result,
there are $2k-j+1$’s $x_{j}$ for all $2\leq j\leq k$ and $2k-j$’s $x_{j}$ for
all $k<j\leq 2k-1$. Therefore, $(x_{j}-x_{i})^{2k-j+1}$ divides $F$ if $2\leq
j\leq k$ and $(x_{j}-x_{i})^{2k-j}$ divides $F$ if $k<j\leq 2k-1$.

Now, by the unique factorization property of multivariate polynomial ring,
there exists a multivariate polynomial $C(x_{1},\dots,x_{2k-1})$ such that

$\displaystyle F(x_{1}\dots,x_{2k-1})$
$\displaystyle=C(x_{1},\dots,x_{2k-1})\prod_{1\leq j\leq
k}x_{j}^{2k-j+1}\prod_{k<j\leq 2k-1}x_{j}^{2k-j}$ $\displaystyle\prod_{1\leq
i<j\leq k}(x_{j}-x_{i})^{2k-j+1}\prod_{\begin{subarray}{c}1\leq i<j,\\\
k<j\leq 2k-1\end{subarray}}(x_{j}-x_{i})^{2k-j},$ (185)
$\displaystyle\eqqcolon C(x_{1},\cdots,x_{2k-1})H(x_{1},\cdots,x_{2k-1}).$

$C(x_{1},\cdots,x_{2k-1})$ will be a constant if we can prove

$\displaystyle\deg H=\deg F.$ (186)

We compute

$\displaystyle\deg H$
$\displaystyle=\sum_{j=1}^{k}(2k-j+1)+\sum_{j=k+1}^{2k-1}(2k-j)+\sum_{j=2}^{k}\sum_{i=1}^{j-1}(2k-j+1)+\sum_{j=k+1}^{2k-1}\sum_{i=1}^{j-1}(2k-j)$
$\displaystyle=\left(\frac{3k^{2}+k}{2}\right)+\left(\frac{k^{2}-k}{2}\right)+\left(\frac{2}{3}k^{3}-\frac{1}{2}k^{2}-\frac{1}{6}k\right)+\left(\frac{2}{3}k^{3}-k^{2}+\frac{1}{3}k\right)$
$\displaystyle=\frac{4}{3}k^{3}+\frac{1}{2}k^{2}+\frac{1}{6}k=\deg F.$ (187)

This concludes the proof. ∎

We are now ready to prove the non-singularity of the coefficient matrix
$\boldsymbol{K}$ as in Eq. 30.

######  Corollary 4.6.

Let $p\in\mathbb{N},\ p\geq 2$. The matrix $\boldsymbol{K}$ define in Eq. 30
is non-singular. As a result, the linear system Eq. 31 always admits a unique
solution.

###### Proof.

We only need to prove $\boldsymbol{E}$ is non-singular, according to Lemma
4.4. Letting $x_{i}=i^{2},\ 1\leq i\leq p-1$, we have

$\displaystyle\boldsymbol{E}=\boldsymbol{F}(1^{2},\dots,(p-1)^{2}).$ (188)

Since all $x_{i}$’s are non-zeros and mutually distinct, therefore Theorem 4.5
implies $\boldsymbol{E}$ is non-singular and the proof is complete. ∎

######  Remark.

By combining Theorem 4.2 and Corollary 4.6, with a few modifications one can
show the linear systems (2.13) in [14] corresponding to 2D weakly singular
integral with singular kernel $s(x)=\nicefrac{{1}}{{|x|}}$ is non-singular for
arbitrary $p$. Consequently, this leads to a unconditional proof of Theorem
3.7 in [14].

##  5 Numerical Results

###  5.1 Computation of Correction Weights

In order to use the modified trapezoidal rules Eqs. 10, LABEL: and 28, one
needs accurate values of weights $\bar{\omega}_{\beta}$’s. In this section, we
describe numerical method for computing $\bar{\omega}_{\beta}$’s.

####  5.1.1 On-diagonal case

We present the computation of the correction weights for the grids
$\mathcal{I}_{0}$, $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ in Eq. 13. We solve
the linear system Eq. 92 in the proof of Theorem 3.1 to compute correction
weight. Compared with the linear system Eq. 21, it has the advantage of
evaluating less limits as $h\to 0$. To find out $\mathfrak{C}(\alpha)$, which
is the RHS of the system Eq. 92, we need $c(h)$ in Eq. 22. We choose a
radially symmetric, practically compactly supported function

$\displaystyle g(x)=\exp(-|x|^{k}),$ (189)

with $k=6$. The reason that we choose $g$ as in Eq. 189, rather than truly
compactly supported one is that we can analytically evaluate the integral

$\displaystyle\int_{\mathbb{R}^{2}\setminus\\{0\\}}g(x)s(x)x^{2\xi}\
\operatorname{d}\\!{x}=\frac{2\Gamma(1.5+\xi_{1})\Gamma(0.5+\xi_{2})\Gamma(\nicefrac{{(2-\alpha+2\xi_{1}+2\xi_{2})}}{{k}})}{k\Gamma(2+\xi_{1}+\xi_{2})},$

for all $0<\alpha<2$ and $\xi\in\mathcal{I}_{2}$. Combining Eq. 92 with
Corollary 3.7 we have

$\displaystyle\boldsymbol{K}\bar{\omega}=\mathfrak{C}(\alpha)=c(h)+h^{6}\sum_{j=0}^{J}C_{j}h^{2j}.$
(190)

Equation 190 enables us to apply twice the Richardson extrapolation to find
out $\mathfrak{C}(\alpha)$:

$\displaystyle
c_{1}(h)=c(\nicefrac{{h}}{{2}})+\frac{c(\nicefrac{{h}}{{2}})-c(h)}{4^{3}-1},\quad
c_{2}(h)=c_{1}(\nicefrac{{h}}{{2}})+\frac{c_{1}(\nicefrac{{h}}{{2}})-c_{1}(h)}{4^{4}-1}.$

We obtain the weights with 20 correct digits by solving
$\bar{\omega}\approx\boldsymbol{K}^{-1}c_{2}(h)$ and ensuring
$|c_{2}(h)-c_{2}(\nicefrac{{h}}{{2}})|<10^{-21}$ for some small enough $h$. In
practice, we found $h=\frac{1}{32}$ sufficiently small to give more than 20
correct digits. Table 1 and Table 2 present the correction weights
corresponding to the singular kernel $s$ in Eq. 9 when $\alpha=0.5,1.5$ and
$p=0,1,2$, respectively.

Table 1: The on-diagonal correction weights for $\alpha=0.5$. $p$ | Grid points | Correction weights
---|---|---
0 | $\mathcal{G}_{0,0}=\\{(0,0)\\}$ | $\omega_{0,0}=9.608446105899650591\times 10^{-1}$
1 | $\mathcal{G}_{0,0}=\\{(0,0)\\}$ | $\omega_{0,0}=9.2275199269460481567\times 10^{-1}$
$\mathcal{G}_{1,0}=\\{(\pm 1,0)\\}$ | $\omega_{1,0}=-3.8305792599451481531\times 10^{-2}$
$\mathcal{G}_{0,1}=\\{(0,\pm 1)\\}$ | $\omega_{0,1}=5.7352101547131603247\times 10^{-2}$
2 | $\mathcal{G}_{0,0}=\\{(0,0)\\}$ | $\omega_{0,0}=0.91354757991861649779\times 10^{-1}$
$\mathcal{G}_{1,0}=\\{(\pm 1,0)\\}$ | $\omega_{1,0}=-4.9714459296827069288\times 10^{-2}$
$\mathcal{G}_{0,1}=\\{(0,\pm 1)\\}$ | $\omega_{0,1}=7.3324618127490001511\times 10^{-2}$
$\mathcal{G}_{2,0}=\\{(\pm 2,0)\\}$ | $\omega_{2,0}=2.2625071864653714109\times 10^{-3}$
$\mathcal{G}_{1,1}=\\{(\pm 1,\pm 1),(\mp 1,\pm 1)\\}$ | $\omega_{1,1}=1.1793189757570510571\times 10^{-3}$
$\mathcal{G}_{0,2}=\\{(0,\pm 2)\\}$ | $\omega_{0,2}=-4.5827886329681250944\times 10^{-3}$
Table 2: The on-diagonal correction weights for $\alpha=1.5$. $p$ | Grid points | Correction weights
---|---|---
0 | $\mathcal{G}_{0,0}=\\{(0,0)\\}$ | $\omega_{0,0}=5.0387797393965760507\times 10^{0}$
1 | $\mathcal{G}_{0,0}=\\{(0,0)\\}$ | $\omega_{0,0}=4.7857569346819649328\times 10^{0}$
$\mathcal{G}_{1,0}=\\{(\pm 1,0)\\}$ | $\omega_{1,0}=1.0971059048869895449\times 10^{-2}$
$\mathcal{G}_{0,1}=\\{(0,\pm 1)\\}$ | $\omega_{0,1}=1.1554034330843566347\times 10^{-1}$
2 | $\mathcal{G}_{0,0}=\\{(0,0)\\}$ | $\omega_{0,0}=4.7305900462046469972\times 10^{0}$
$\mathcal{G}_{1,0}=\\{(\pm 1,0)\\}$ | $\omega_{1,0}=1.7018648395611181367\times 10^{-2}$
$\mathcal{G}_{0,1}=\\{(0,\pm 1)\\}$ | $\omega_{0,1}=1.3848756814856511801\times 10^{-1}$
$\mathcal{G}_{2,0}=\\{(\pm 2,0)\\}$ | $\omega_{2,0}=-4.4305641359382777203\times 10^{-3}$
$\mathcal{G}_{1,1}=\\{(\pm 1,\pm 1),(\mp 1,\pm 1)\\}$ | $\omega_{1,1}=5.8373335985059124819\times 10^{-3}$
$\mathcal{G}_{0,2}=\\{(0,\pm 2)\\}$ | $\omega_{0,2}=-8.6554730092853198753\times 10^{-3}$

####  5.1.2 Off-diagonal Case

We present correction weights for grids $\mathcal{I}_{2}$, $\mathcal{I}_{3}$
and $\mathcal{I}_{4}$ in Eq. 27. In off-diagonal case, the procedure for
computing correction weights is similar with that of the on-diagonal case
described above. We choose $g$ as in Eq. 189 with $k=8$. The analytical values
of the integrals are

$\displaystyle\int_{\mathbb{R}^{2}\setminus\\{0\\}}g(x)s(x)x^{2\xi-(1,1)}\
\operatorname{d}\\!{x}=\frac{2\Gamma(0.5+\xi_{1})\Gamma(0.5+\xi_{2})\Gamma(\nicefrac{{(2\xi_{1}+2\xi_{2}-\alpha)}}{{k}})}{k\Gamma(1+\xi_{1}+\xi_{2})},$

for all $0<\alpha<2$ and $\xi\in\mathcal{I}_{4}$. Table 3 and Table 4 present
the correction weights corresponding to singular kernel in Eq. 24 for
$\alpha=0.5,1.5$ and $p=2,3,4$, respectively.

Table 3: The off-diagonal corrections weights for $\alpha=0.5$. $p$ | Discretization points | Correction weights
---|---|---
2 | $\mathcal{G}_{1,1}=\\{(\pm 1,\pm 1),(\mp 1,\pm 1)\\}$ | $\omega_{1,1}=0.0286760507735658016236634025724$
3 | $\mathcal{G}_{1,1}=\\{(\pm 1,\pm 1),(\mp 1,\pm 1)\\}$ | $\omega_{1,1}=0.0470072053054383020013851917611$
$\mathcal{G}_{2,1}=\\{(\pm\mp 2,\pm\pm 1),(\pm\mp 1,\pm\pm 2)\\}$ | $\omega_{2,1}=-0.00458278863296812509443044729718$
4 | $\mathcal{G}_{1,1}=\\{(\pm 1,\pm 1),(\mp 1,\pm 1)\\}$ | $\omega_{1,1}=0.058498692309201978109$
$\mathcal{G}_{2,1}=\\{(\pm\mp 2,\pm\pm 1),(\pm\mp 1,\pm\pm 2)\\}$ | $\omega_{2,1}=-0.0092844902620645196084$
$\mathcal{G}_{3,1}=\\{(\pm\mp 3,\pm\pm 1),(\pm\mp 1,\pm\pm 3)\\}$ | $\omega_{3,1}=0.0010440418727854435399$
$\mathcal{G}_{2,2}=\\{(\pm 2,\pm 2),(\mp 2,\pm 2)\\}$ | $\omega_{2,2}=0.00026276706897731017725$
Table 4: The off-diagonal corrections weights for $\alpha=1.5$. $p$ | Discretization points | Correction weights
---|---|---
2 | $\mathcal{G}_{1,1}=\\{(\pm 1,\pm 1)\\}$ | $\omega_{1,1}=0.0577701716542178317339761161235$
3 | $\mathcal{G}_{1,1}=\\{(\pm\mp 1,\pm 1)\\}$ | $\omega_{1,1}=0.0923920636913591112353501723599$
$\mathcal{G}_{2,1}=\\{(\pm 2,\pm 1),(\pm 1,\pm 2)\\}$ | $\omega_{2,1}=-0.0086554730092853198753435140591$
4 | $\mathcal{G}_{1,1}=\\{(\pm 1,\pm 1)\\}$ | $\omega_{1,1}=0.11372612810258708544$
$\mathcal{G}_{2,1}=\\{(\pm 2,\pm 1),(\pm 1,\pm 2)\\}$ | $\omega_{2,1}=-0.017474957624915655234$
$\mathcal{G}_{3,1}=\\{(\pm 3,\pm 1),(\pm 1,\pm 3)\\}$ | $\omega_{3,1}=0.0018475475899836517452$
$\mathcal{G}_{2,2}=\\{(\pm 2,\pm 2),(\mp 2,\pm 2)\\}$ | $\omega_{2,2}=0.00071464712784786418872$

###  5.2 Order of Convergence of $Q_{h}^{p}$

####  5.2.1 On-diagonal case

We verify the order of convergence of the corrected trapezoidal rule Eq. 17
for the on-diagonal weakly singular integral $I_{11}$ in Eq. 7. We choose a
compactly supported function

$\displaystyle\phi(x)$
$\displaystyle=(1+x_{1}+x_{1}^{2})(1+x_{2}+x_{2}^{2})\big{(}(1-x_{1}^{2})_{+}(1-x_{2}^{2})_{+}\big{)}^{7},$
(191)

where $x_{+}\coloneqq\max\\{x,0\\}$. Note that $\phi\in
C_{c}^{6}(\mathbb{R}^{2})$ such that $\partial^{k}\phi(0)\not=0,\ \forall
k\in\mathbb{N}_{0}^{2}:|k|_{1}\leq 6$ and $supp(\phi)=[-1,1]^{2}$. We use
$\phi$ to test Theorem 3.1 for $p=0,1,2$. Note that due to limited smoothness
of this $\phi$, Theorem 3.1 guarantees the order of convergence $2p+4-\alpha$
for $p=0,1$, but is inconclusive for $p=2$. However we still observe the high
order of accuracy. This implies there is room for further relaxation on
smoothness criteria on $\phi$. The analytical integral
$I=\int_{\mathbb{R}^{2}\setminus\\{0\\}}\phi(x)s(x)\ \operatorname{d}\\!{x}$
for arbitrary $0<\alpha<2$ can be found by Mathematica. However the resulting
formula is extremely long, therefore we only evaluate the integral for
$\alpha=0.5,\ 1.5$. We evaluate $|I-Q_{h}^{p}|$ for different values of $h$,
where $I$ denotes the true values, and perform linear regressions in log-log
plots to find out the order of convergence, as shown in Fig. 3. We find that
the numerical results in general match very well with theoretically predicted
order of accuracy $2p+4-\alpha$ for all $p=0,1,2$ and $\alpha=0.5,1.5$. When
the value of $p$ becomes large and $h$ becomes small, the round-off errors
dominate. In this case, multi-precision will help when one chooses large $p$.

(a) $\alpha=0.5$

(b) $\alpha=1.5$

Figure 3: Numerical orders of accuracy: Log-log plots of the quadrature error
against $h$ for the modified trapezoidal rule Eq. 17 corresponding to the on-
diagonal singular kernel $s(x)=\frac{x_{1}^{2}}{|x|^{2+\alpha}}$ and
$p=0,1,2$, (a) $\alpha=0.5$, (b) $\alpha=1.5$.

####  5.2.2 Off-diagonal case

We verify the order of convergence of the corrected trapezoidal rule Eq. 28
for the off-diagonal weakly singular integral $I_{12}$ in Eq. 7, we choose a
compact supported function

$\displaystyle\phi(x)$
$\displaystyle=(1+x_{1})(1+x_{2})\big{(}(1-x_{1}^{2})_{+}(1-x_{2}^{2})_{+}\big{)}^{7},$
(192)

Note that $\phi\in C_{c}^{6}(\mathbb{R}^{2})$ and $supp(\phi)=[-1,1]^{2}$. We
use $\phi$ to test Theorem 3.11 for $p=1,2,3$. Note that due to limited
smoothness of this $\phi$, Theorem 3.11 guarantees the order of accuracy of
$2p+2-\alpha$ for $p=1,2$, but is inconclusive for $p=3$. However we still
observe the high order of accuracy. This implies there are room for further
relaxation on smoothness criteria on $\phi$. The true value of the integral
$I=\int_{\mathbb{R}^{2}\setminus\\{0\\}}\phi(x)s(x)\ \operatorname{d}\\!{x}$
for $\alpha=0.5,1.5$ can be found by Mathematica. We evaluate $|I-Q_{h}^{p}|$
for different values of $h$, where $I$ denotes the true values and perform
linear regressions in log-log plots to find out the order of convergence, as
shown in Fig. 4. We find that the numerical results in general match very well
with the theoretically predicted order of accuracy $2p+2-\alpha$ for $p=1,2,3$
and $\alpha=0.5,1.5$.

(a) $\alpha=0.5$

(b) $\alpha=1.5$

Figure 4: Numerical orders of accuracy: Log-log plots of the quadrature error
against $h$ for the modified trapezoidal rule Eq. 28 corresponding to the off-
diagonal singular kernel $s(x)=\frac{x_{1}x_{2}}{|x|^{2+\alpha}}$ and
$p=1,2,3$, (a) $\alpha=0.5$, (b) $\alpha=1.5$.

##  6 Conclusion

We proposed an arbitrarily high-order modified trapezoidal rules for both on-
diagonal and off-diagonal weakly singular integrals
$I_{i,j}=\int_{\mathbb{R}^{2}}\phi(x)\frac{x_{i}x_{j}}{|x|^{2+\alpha}}\operatorname{d}\\!{x}$,
$i,j\in\\{1,2\\}$ with some sufficient smooth function $\phi$ with compact
support. The quadrature rule is a punctured-hole trapezoidal rule with
correction terms. We proved the order of accuracy of the modified quadrature
given the number of correction layers. We focused on error due to singularity
in this paper. The quadrature rule can be combined with any boundary error
correction for regular functions $\phi$ without compact support to attain
high-order convergence. The correction weights can be pre-computed and stored
for future use. We tabulated correction weights in more-than-17 digits in
scientific form for $\alpha=0.5,1.5$ and for both on-diagonal and off-diagonal
case. We provided theoretical guarantee that the quadrature rule works with
arbitrary $p$, the number associated with correction layers, though we do not
recommend to use large $p$ since numerical evidence suggested that linear
systems Eq. 21 and Eq. 31 become increasingly ill-conditioned for large $p$.
By using the method suggested in Remark Remark, one can prove the non-
singularity left in the work of [14]. For future works, one can further relax
the smoothness criteria for $\phi$ in Theorems 3.4, LABEL: and 3.11.
Alternatively, the modified trapezoidal rule is of course not confined to
singular kernels described here. With proper modification, one can easily
adapt this rule to other singular kernels in 2D such as $x_{1}\log{(|x|)}$ or
$\frac{P(x_{1},x_{2})}{|x|^{\nicefrac{{(2k+1)}}{{2}}}}$ where $P(x_{1},x_{2})$
is a irreducible polynomial in $x_{1},x_{2}$ with degree $k$. Another possible
direction is to generalize the quadrature rule to arbitrary $n$ dimensions, as
some results in this paper are already set up in $n$ dimensions.

## References

  * [1] J.C Aguilar and Y Chen.  High-order corrected trapezoidal quadrature rules for functions with a logarithmic singularity in 2-d.  Computers & Mathematics with Applications, 44(8):1031–1039, 2002\.
  * [2] J.C. Aguilar and Yu Chen.  High-order corrected trapezoidal quadrature rules for the coulomb potential in three dimensions.  Computers & Mathematics with Applications, 49(4):625–631, 2005\.
  * [3] BK Alpert.  Rapidly-convergent quadratures for integral operators with singular kernels.  Lawrence Berkeley National Laboratory. LBNL Report: LBL-30091, 1990\.
  * [4] Kendall Atkinson.  Quadrature of singular integrands over surfaces.  ETNA. Electronic Transactions on Numerical Analysis [electronic only], 17, 01 2004.
  * [5] Oscar P. Bruno and Leonid A. Kunyansky.  A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications.  Journal of Computational Physics, 169:80–110, 2001.
  * [6] Michael G. Duffy.  Quadrature over a pyramid or cube of integrands with a singularity at a vertex.  SIAM Journal on Numerical Analysis, 19(6):1260–1262, 1982.
  * [7] L. Grafakos.  Classical Fourier Analysis.  Graduate Texts in Mathematics. Springer New York, 2014.
  * [8] Hansen Ha.  Non-Gaussian Stochastic Dynamics with Deterministic Numerical Tools.  PhD thesis, 2018.  Copyright - Database copyright ProQuest LLC; ProQuest does not claim copyright in the individual underlying works; Last updated - 2021-08-03.
  * [9] Sharad Kapur and Vladimir Rokhlin.  High-order corrected trapezoidal quadrature rules for singular functions.  SIAM Journal on Numerical Analysis, 34(4):1331–1356, 1997.
  * [10] Patrick Keast and James N. Lyness.  On the structure of fully symmetric multidimensional quadrature rules.  SIAM Journal on Numerical Analysis, 16:11–29, 1979.
  * [11] M.A. Khayat and D.R. Wilton.  Numerical evaluation of singular and near-singular potential integrals.  IEEE Transactions on Antennas and Propagation, 53(10):3180–3190, 2005.
  * [12] Anna Lischke, Guofei Pang, Mamikon Gulian, Fangying Song, Christian Glusa, Xiaoning Zheng, Zhiping Mao, Wei Cai, Mark M. Meerschaert, Mark Ainsworth, and George Em Karniadakis.  What is the fractional laplacian? a comparative review with new results.  Journal of Computational Physics, 404:109009, 2020.
  * [13] Oana Marin, Katarina Gustavsson, and Anna-Karin Tornberg.  A highly accurate boundary treatment for confined stokes flow.  Computers & Fluids, 66:215 – 230, 2012.
  * [14] Oana Marin, Olof Runborg, and Anna-Karin Tornberg.  Corrected trapezoidal rules for a class of singular functions.  IMA Journal of Numerical Analysis, 34(4):1509–1540, 2014.
  * [15] S. E. Mousavi and N. Sukumar.  Generalized duffy transformation for integrating vertex singularities.  Computational Mechanics, 45:127–140, 2009.
  * [16] V. Rokhlin.  End-point corrected trapezoidal quadrature rules for singular functions.  Computers & Mathematics with Applications, 20(7):51–62, 1990.
  * [17] Avram Sidi.  Application of class i m variable transformations to numerical integration over surfaces of spheres.  Journal of Computational and Applied Mathematics, 184:475–492, 2005\.