Composite trapezoidal quadrature for computing hypersingular integrals on interval

1.   Introduction

Consider the hypersingular integrals

where $N = 1, 2, \cdots$, $\int_{a}^{b} \hskip-.22in =$ denotes a hypersingular integrals with $p+1$ is singular order and $y$ is the singular point.

Hypersingular integrals (1.1) originated in the 1970s and 1980s, and have gradually gained attention as the boundary element method (BEM) developed. In the boundary element method, hypersingular integrals often appear when dealing with problems in elasticity, fracture mechanics, fluid mechanics, electromagnetic fields, heat conduction, and computational biology.

The definition of hypersingular integral can be found by the derivative of the Cauchy principal integral,

subtraction of the singularity

where $r > p$ and Hadamard finite-part integrals

where

and so on. These definitions can be proved mathematically equal under certain conditions. For the case $p = 1, 2$, the Hadamard finite-part integrals definition is expressed as:

Because of the hypersingular kernel, classical numerical methods cannot be used directly to compute the hypersingular integral. Different numerical methods can be obtained by different definitions. In 1975, original work of Kutt ^[22] stressed Gaussian quadrature, then more papers ^[2,3,7,8,9] on the numerical computation hypersingular based on the definition of the derivative of Cauchy principal integral and subtraction of the singularity. The Newton-Cote methods were presented by Linz ^[19], while the density function was approximated by the trapezoidal quadrature and Simpson quadrature, and the hypersingular kernel is calculated analysis. Then, the superconvergence phenomenon ^[14,15,16] of hypersingular integrals ^[23,24,28] are investigated based on the definition of Hadamard finite-part integrals ^{[4,11,12,13,27]}. There are also the extrapolation methods ^[6,10,17,20] and some other methods ^[21] to numerical evaluation hypersingular integrals, which are also based on the Hadamard finite-part integral definition. The Gaussian integral formula can attain complete accuracy for low-order polynomials. For the integration of nonstandard intervals or nonstandard weight functions, it is requisite to recalculate the nodes and weights, which might be more complex. The accuracy of extrapolation algorithms hinges on the regularity of the function and the design of the algorithm. If the regularity of the function is not high, the extrapolation results might be unreliable. During multiple extrapolation processes, errors may accumulate, resulting in inaccurate final results.

Now, we recall the Newton-Cote rule first. Let $a = x_0 < x_1 < \cdots < x_{n-1} < x_{n} = b$ be a uniform partition of $[a, b]$ with $h = (b-a)/n$. We also define the piecewise Lagrangian polynomial interpolation of degree $k$ by

where

and a linear transformation

changes the subinterval $[x_{i}, x_{i+1}]$ into $[-1, 1]$.

Replacing $u(x)$ in (1.1) by ${u}_{kn}(x)$ gives the composite Newton-Cotes quadrature

where $\mathcal{Q}_{kn}^{p}(u; y)$ denotes quadrature rule for $p+1$ singular order and ${E}_{kn}^p(u; y)$ denotes error estimate, and

is the Cotes coefficient. Different from the Cotes coefficient of the Riemann integral, for $k\ge2$ it is not easy to get the Cotes coefficient for hypersingular integral as we have to compute equation (1.11). There are lots of works that investigate the hypersingular integral with $p = 1$. The error estimate is obtained by Linz ^[19]

and

is the distance of singular point to the mesh point.

A more precise estimate was given in ^[25], where

with $u(x) \in C^{k+\alpha}[a, b], \alpha\in[0, 1)$. For the case $y = x_{i}, i = 0, 1, \cdots, n$, Yu^[27] has presented the error estimate with

while in the recent paper, Wu^[23] has given the modify rule and the $O(h^{k}), k = 1, 2$ convergence rate is obtained. Superconvergence phenomenon of (1.1) has been investigated by ^[24,25] and the error estimations are presented as

for odd $k$, and

for even $k$, and

is the distance of singular point $s$ to the boundary of $a$ and $b$.

In the following, let $y = x_{m+1}+\frac{(\xi+1)h}{2}$. We present the error expansion of the second order finite-part integral as below:

then we get

where $Q_{1}(\tau) = -1+\frac{\tau}{2}\log\left|\frac{1+\tau}{1-\tau}\right|$ is the Legendre function of second kind and $Q'_{1}(\tau), Q''_{1}(\tau)$ are the first and second derivatives of $Q_{1}(\tau)$, then we have

and

By similar calculation, we also have

Because of singularity of hypersingular integral, trapezoidal quadrature has not been used to compute (1.1). From the error expansion above, we know that the trapezoidal quadrature is divergence $O(h^{-p}), p = 1, 2$ in general, which is the reason that the trapezoidal rule cannot be used to compute the singular or hypersingular integral. In order to get the simple calculation scheme, trapezoidal quadrature for numerical evaluation of the hypersingular integral is presented. Compared with the numerical scheme just compute the density function presented by Linz, (see Eq (1.10)) and we compute the density function and singular kernel at the same time.

In this paper, we pay our attentions to trapezoidal quadrature for the numerical evaluation hypersingular integral. Furthermore, the relationship of hypersingular integral, supersingular integral, Cauchy principal value integrals, and Riemann integral is also illustrated. In the following, the error functional of trapezoidal quadrature is presented with the first part $h^{-1}u(y)$ which is divergence; see Eq (1.19). From the error expansion, the special function is studied which is related with the function $\frac{1}{\sin^2 x}$ and $\cot x$. Then, we prove that the convergence rate can reach $O(h^2)$ when singular point $s$ is located at the center of interval. Trapezoidal quadrature is presented to approximate the hypersingular integral and the modified trapezoidal quadrature are also given to illustrate the asymptotic expansion which can be used to compute the hypersingular integral. With the help of the modified trapezoidal quadrature, we can compute the hypersingular integral more easy and this rule can be used to solve the hypersingular integral equation. At last, numerical results show the relationship between the Riemann integral and the hypersingular integral.

2.   Main results

In this part, we present our quadrature, different from $u_{1n}(x)$ in the (1.7) substituted density function $u(x)$, then we have

when $x \in (x_{i}, x_{i+1})$. Replacing $u(x)$ in (1.1) with $u_{1n}(x)$ gives

where $E_{1n}^p(u; y)$ is the error functional for $k = 1$ as the trapezoidal rule, and

is the Cote coefficients.

Before, we present the error functional of trapezoidal quadrature for hypersingular integrals. We define the special function

which comes from the error functional by linear transformation, and we also define

and

Theorem 2.1. Assume $u(x)\in C^{p+2}[a, b], p = 1, 2$. Let $Q_{1n}^p(u; y)$ be computed by (2.2). Assume that $y = x_{m}+(1+\xi)h/2, s\neq {x}_{m}$, and $S''_{1}(\xi), S'_{1}(\xi), S_{1}(\xi)$ are defined as (2.10), (2.9), and (2.8), and there hold that

where

$\gamma(\xi)$ is defined as (1.13) and $\eta(y)$ is defined as (1.18).

3.   Proof of the Theorem 2.1

The asymptotic expansion of the composite trapezoidal quadrature to compute hypersingular integrals is presented.

3.1.   The case of $p=1$

In order to prove Theorem 2.1, we give the lemmas as below.

Lemma 3.1. With the same condition of Theorem 2.1 as $p = 1$, it holds that

where

and

where $\alpha_{1i}, \alpha_{2i}, \beta_{1i}, \beta_{2i}, \beta_{3i}\in (x_{i}, x_{i+1})$.

Proof: By Taylor expansion of $u_{1n}(x)$ at $x$, we get

and

Similarly, we have

and

and

Combining (3.7)–(3.10) together, we get the results.

Lemma 3.2. Let $c_{i} = 2(y- x_i)/h -1, 0 \leq i \leq n$, for $y \in (x_{m}, x_{m+1})$. Then, we have

and

and

Proof: By the Eqs (1.6) and (1.9), we have

and

For $i\ne m$, we can obtain it similarly.

Let

and

Lemma 3.3. Let $u(x)\in C^{p+2}[a, b]$ and $E_m^{p}(x), p = 1, 2$ be defined by (3.18) and (3.19) with $y\neq {x}_m, m = 0, 1, 2, \cdots, n$, and there holds

and

Proof: For the case $p = 1$, as $u(x) \in C^{3}[a, b]$, we get $E_{m}^{1}(x) \in C^{3}[a, b]$, and

where we have used

Then, we have

The proof of this lemma is finished.

Lemma 3.4. For $\xi \in(-1, 1)$, and $m \geq1$, we have

and

Proof: By (2.5) of $\phi_{1}(\xi)$, we have

then, we have

and

Noting that $y = a+(m+\frac{\xi+1}{2})h$, we have $\frac{2(y-a)}{h} = \xi+2m+1$, and

Since $b = a+nh$, we have $\frac{2(b-y)}{h} = 2(n-m)-1-\xi$, and

Thus, we finished the proof of lemma.

Lemma 3.5. With the same condition of Theorem 2.1, we have

and

where

and

Proof: By Taylor expansion of $u(x), u(x_{i}), u(x_{i+1})$ at the point $y$, we have

and

which completed the the proof of (3.33). The proof (3.34) can be similarly obtained.

Proof of Theorem 2.1: By Lemma 3.1, we get

where we have used the (3.18) of $E_m^1(x)$

and

and

Now, we estimate $\mathcal{R}^1_n(u; y)$. By using $(x-y)^2 = (x-\hat{x}_{i}+\hat{x}_{i}-y)^2 = (x-\hat{x}_{i})^2+2(x-\hat{x}_{i})(\hat{x}_{i}-y)^2+(\hat{x}_{i}-y)^2$ and $x-y = x-\hat{x}_{i}+\hat{x}_{i}-y$, we have

For the first part of $R^1_{1}(u; y)$, we have

For the second part of $R^1_{1}(u; y)$, we similarly have

For the third part of $R^1_{1}(u; y)$,

For the fourth part of $R^1_{1}(u; y)$,

For the last part of $R^1_{1}(u; y)$, we have

where $\xi_{2i}, \xi_{3i}, \beta_{2i}, \beta_{3i} \in (x_{i}, x_{i+1})$.

For $R^1_{2}(u; y)$ and $R^1_{3}(u; y)$, by Lemmas 3.3 and 3.4, then, we get

and we finished the proof of Theorem 2.1.

3.2.   The case of $p=2$

By Lemma 3.1, we have

where we have used the identity of (3.19),

and

For $R^2_{2}(u; y)$ and $R^3_{2}(u; y)$, by Lemmas 3.3 and 3.4, we get

and the proof Theorem 2.1 is completed.

Corollary 3.1. Under same assumption of Theorem 2.1, we have

and $\eta(y)$ is defined as (1.18).

From Theorem 2.1, modify trapezoidal quadrature Ⅰ is presented

as well as modify trapezoidal quadrature Ⅱ

and error functional

Then, we have the following.

Corollary 3.2. Under the same assumption of Theorem 2.1, for the modify trapezoidal quadrature I, we have

where $\gamma(\xi)$ is defined as (1.13).

In the year of 2013, we have presented the results ^[18]

where $E_{1n}^0(f; y) = \int_{a}^{b} \hskip-.20in- \ \frac{u(x)}{x-y}dx-\int_{a}^{b} \hskip-.20in- \ \frac{u_{1n}(x)}{x_{i}-y}dx$, and

are the Cote coefficients.

Theorem 3.1. Assume $u(x)\in C^{2}[a, b]$. For the trapezoidal quadrature $Q_{0n}^0(u; y)$ defined as (3.66), assume that $y = x_{[m]}+(1+\xi)h/2$, and there holds

Now, we present the following theorem:

Theorem 3.2. Assume $u(x)\in C^{p+2}[a, b]$. Let $Q_{1n}^p(u; y)$ be computed by (2.2), assume that $y = x_{m}+(1+\xi)h/2, y\neq \hat{x}_{m}$, and there holds

where

and $\gamma(\xi)$ is defined as (1.13).

This theorem can be proved similarly as Theorem 1, and here we omit it.

4.   The calculation of $S_{1}(\xi), S'_{1}(\xi), S''_{1}(\xi)$

Let

be the Legendre function of the second kind and $Q_{n}(x)$ associated with the $P_{n}(x)$, defined by (cf.^[1]).

We also define

Then, by (4.2) of $W$,

and it follows that

Then, we have

where we have used the identify

5.   Numerical examples

In this part, some examples are reported to illustrate our theorem.

Example 5.1. Consider the hypersingular integral with $u(x) = x^4+1$ and with $a = 0, b = 1$, and

In Tables 1 and 2, error and a posterior estimation of the trapezoidal rule are are presented, and the extrapolation method can be found in reference ^[12].

The uniform meshes are chosen to test convergence rate of the trapezoidal quadrature $Q_{1n}^{1}(u; y)$, modified Ⅰ trapezoidal quadrature ${\tilde{Q}}_{1n}^{1}(u; y)$, and modified Ⅱ trapezoidal quadrature $\tilde{\tilde{Q}}_{1n}^{1}(u; y)$ with the dynamic point with $y = x_{[n/4]}+(1+\xi)h/2$.

Table 3 shows that when the local coordinate of singular point is $\xi$, the modified Ⅱ trapezoidal quadrature $\tilde{\tilde{Q}}_{1n}^{1}(u; y)$ reaches $O(h^2)$. From the Table 4, the modified Ⅰ trapezoidal quadrature ${\tilde{Q}}_{1n}^{1}(u; y)$ can reach $O(h^2)$ at $\xi = 0$; as for $\xi\ne 0$, there are no convergence rates. From the Table 5, for the trapezoidal quadrature ${Q}_{1n}^{1}(u; y)$, there are no convergence rates.

Example 5.2. Consider the hypersingular integral $u(x) = x^3$, and

The uniform meshes are chosen to test convergence rate of the trapezoidal quadrature $Q_{1n}^{2}(u; y)$, modified Ⅰ trapezoidal quadrature ${\tilde{Q}}_{1n}^{2}(u; y)$, and modified Ⅱ trapezoidal quadrature $\tilde{\tilde{Q}}_{1n}^{2}(u; y)$ with $y = x_{[n/4]}+(1+\xi)h/2$ and $y = a+(1+\xi)h/2$.

For the case $p = 2$ with $y = x_{[n/4]}+(1+\xi)h/2$, the Table 6 shows that when the local coordinate is $\xi$ of the singular point, modified Ⅱ trapezoidal quadrature $\tilde{\tilde{Q}}_{1n}^{2}(u; y)$ reaches $O(h^2)$. Table 7 shows that the modified Ⅰ trapezoidal quadrature ${\tilde{Q}}_{1n}^{2}(u; y)$ can reach $O(h^2)$ at the $\xi = 0$; as for $\xi\ne 0$, there are no convergence rates. Table 8 shows the trapezoidal quadrature ${Q}_{1n}^{2}(u; y)$, and there are no convergence rates.

Table 9 shows $y = a+(\xi+1) h/2$, and modified Ⅱ trapezoidal quadrature $\tilde{\tilde{Q}}_{1n}^{2}(u; y)$ reaches $O(h)$. Table 10 shows that the modified Ⅰ trapezoidal quadrature ${\tilde{Q}}_{1n}^{2}(u; y)$ can reach $O(h)$ at the $\xi = 0$; as for $\xi\ne 0$, there are no convergence rates which agree with our Theorem 2.1. Table 11 shows the trapezoidal quadrature ${Q}_{1n}^{2}(u; y)$, and there are also convergence rates $O(h)$ because of $f(x) = 0$ at the end of intervals.

For the case $p = 1$ with $y = x_{[n/4]}+(1+\xi)h/2$, Table 12 shows that when the local coordinate is $\xi$ of singular point, modified Ⅱ trapezoidal quadrature $\tilde{\tilde{Q}}_{1n}^{1}(u; y)$ reaches $O(h^2)$. Table 13 shows that modified Ⅰ trapezoidal quadrature ${\tilde{Q}}_{1n}^{1}(u; y)$ can reach $O(h^2)$ at the $\xi = 0$; as for $\xi\ne 0$, there are no convergence rates which agree with our Theorem 2.1. Table 14 shows the trapezoidal quadrature ${Q}_{1n}^{1}(u; y)$, and there are no convergence rates.

Table 15 shows $y = a+(\xi+1) h/2$ and that modified Ⅱ trapezoidal quadrature $\tilde{\tilde{Q}}_{1n}^{1}(u; y)$ reaches $O(h^2)$. Table 16 shows that modified Ⅰ trapezoidal quadrature ${\tilde{Q}}_{1n}^{1}(u; y)$ can reach $O(h^2)$ at the $\xi = 0$; as for $\xi\ne 0$, there is also convergence rate $O(h^2)$. Table 17 shows the trapezoidal quadrature ${Q}_{1n}^{1}(u; y)$, and there are also convergence rate $O(h^2)$ because of $f(x) = 0$ at the end of intervals.

6.   Concluding remarks

In this paper, hypersingular integrals are studied by the composite trapezoidal quadrature on the interval. With the linear transformation of subintervals to the identity intervals, the error functional is obtained related with the special function. The modify trapezoidal rule is presented, with the help of second term equal zeros which gives the superconvergence phenomenon. In this paper of ^[18], we have proved the results that the error functional has related with the special function

for the Cauchy principal integrals. In this paper, we have presented the results that the error functional has related with the special function

for the hypersingular integrals. Based on the above results, we conjecture that there are certain relationship between the hypersingular integral $\int_{a}^{b} \hskip-.25in = \frac{u(x)}{(x-y)^{p+1}}dx, p = 0, 1, \cdots$, and the special function

We will give further investigation to illustrate the relationship.

The research prospects regarding the calculation methods of hypersingular integrals that might center on the following aspects in the future. We construct efficient numerical approaches to handle two-dimensional and three-dimensional problems in practical applications. For fractional order hyper singular integrals emerging in fractional order partial differential equations, corresponding numerical methods are developed for numerical computations and theoretical investigations. We establish cost-effective program packages suitable for the engineering community to facilitate the application of super singular integration calculation methods. With the advancement of computing technology and the introduction of novel algorithms, the calculation method of hypersingular integrals will become more efficient and precise, thereby exerting a greater impact in various application fields.

Symbol list

$\int_{a}^{b} \hskip-.22in =$ — —- —- hypersingular integrals

$\int_{a}^{b} \hskip-.22in -$ — —- —- Cauchy Principal integral

${I}_{p}(u, y)$ — —- —- denotes a hypersingular integrals with $p+1$ singular order

${u}_{kn}(x)$ — —- —- Lagrangian polynomial interpolation of degree $k$

$\mathcal{Q}_{kn}^{p}(u; y)$ — —- —- $k$ degree composite Newton-Cotes quadrature for ${I}_{p}(u, y)$

${E}_{kn}^{p}(u; y)$ — —- —- error functional of $k$ degree composite Newton-Cotes quadrature for ${I}_{p}(u, y)$

$\gamma(\xi)$ — —- —- distance of singular point to the mesh point

$\eta(y)$ — —- —- distance of singular point $s$ to the boundary point.

$\omega_i^p(y)$ — —- —- Cote coefficients of ${I}_{p}(u, y)$ with ${u}_{kn}(x)$

$\phi_{1}(x)$ — —- —- special function of ${I}_{p}(u, y)$

$S_{1}(\xi)$ — —- —- linear operator

$S'_{1}(\xi)$ — —- —- first order derivative of linear operator

$S''_{1}(\xi)$ — —- —- second order derivative of linear operator

$\mathcal{R}_n^p(u; y)$ — —- —- remain part of ${E}_{kn}^{p}(u; y)$

$E_{m}^1(x)$ — —- —- errors of

subinterval $[x_{m}, x_{m+1}]$

$E_{m}^2(x)$ — —- —- errors of

subinterval $[x_{m}, x_{m+1}]$

$R^1_{1}(u; y),$ — —- —- first part of $\mathcal{R}_n^1(u; y)$

$R^1_{2}(u; y)$ — —- —- second part of $\mathcal{R}_n^1(u; y)$

$R^1_{3}(u; y)$ — —- —- third part of $\mathcal{R}_n^1(u; y)$

$R^2_{1}(u; y),$ — —- —- first part of $\mathcal{R}_n^2(u; y)$

$R^2_{2}(u; y)$ — —- —- second part of $\mathcal{R}_n^2(u; y)$

$R^2_{3}(u; y)$ — —- —- third part of $\mathcal{R}_n^2(u; y)$

Author contributions

Xiaoping Zhang: conceptualization, methodology, investigation, funding acquisition, writing-original draft; Jin Li: conceptualization, investigation, writing-review and editing methodology, software, writing-original draft. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

The work of Jin Li was supported by Natural Science Foundation of Shandong Province (Grant No. ZR2022MA003). The authors also gratefully acknowledge the helpful comments and suggestions of the reviewers, which have improved the presentation.

Conflict of interest

The authors declare there is no conflict of interest.