# Fast and Accurate Solvers for Weakly Singular Integral Equations

Laurence Grammont, Rekha Panditra Kulkarni, Paulo Vasconcelos

► To cite this version:

Laurence Grammont, Rekha Panditra Kulkarni, Paulo Vasconcelos. Fast and Accurate Solvers for Weakly Singular Integral Equations. Numerical Algorithms, 2022, 10.1007/s11075-022-01376-x. hal-04053854

HAL Id: hal-04053854

https://hal.science/hal-04053854

Submitted on 31 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fast and Accurate Solvers for Weakly Singular Integral Equations

Laurence Grammont<fnref n="1" />, Rekha P. Kulkarni<fnref n="2" /> and Paulo B. Vasconcelos<fnref n="3" />*

*Corresponding author(s). E-mail(s): pjv@fep.up.pt;
Contributing authors: laurence.grammont@univ-st-etienne.fr;
rpkiitb@gmail.com;

## Abstract

Consider an integral equation $\lambda u - Tu = f$, where $T$ is an integral operator, defined on $C[0, 1]$, with a kernel having an algebraic or a logarithmic singularity. Let $\pi_m$ denote an interpolatory projection onto a space of piecewise polynomials of degree $\le r-1$ with respect to a graded partition of $[0, 1]$ consisting of $m$ subintervals. In the product integration method, an approximate solution is obtained by solving $\lambda u_m - T\pi_m u_m = f$. As in order to achieve a desired accuracy, one may have to choose $m$ large, we find approximations of $u_m$ using a discrete modified projection method and its iterative version. We define a two-grid iteration scheme based on this method and show that it needs less number of iterates than the two-grid iteration scheme associated with the discrete collocation method. Numerical results are given which validate the theoretical results.

**Keywords:** Algebraic singularity, Logarithmic singularity, Collocation method, Interpolatory projection, Graded mesh, Piecewise polynomial approximation, Product integration, Weakly singular integral operator.

<footnote n="1">¹*Université de Lyon, Institut Camille Jordan, UMR5208, 23, rue du Dr Paul Michelon, Saint-Étienne Cedex 2, 42023, France.</footnote>

<footnote n="2">²Department of Mathematics, I.I.T. Bombay, Powai, Mumbai, 400076, India.</footnote>

<footnote n="3">³Faculdade de Economia da Universidade do Porto and Centro Mathemática da Universidade do Porto, Street, City, 610101, State, Portugal.</footnote>
2 Fast and Accurate Solvers for Weakly Singular Integral Equations

MSC Classification: 45L05, 65J15, 65R20

# 1 Introduction

Consider the integral operator

$$ (Tu)(s) = \int_{0}^{1} g_{\alpha}(|s-t|) u(t) dt, \quad 0 \le s \le 1, $$

with

$$ g_{\alpha}(s) = s^{\alpha-1}, \quad 0 < \alpha < 1 \quad \text{and} \quad g_{1}(s) = \log(s). $$

The former corresponds to the so called algebraic singularity and the latter, $\alpha = 1$ following the convention in [1], to the logarithm singularity. Then $T$ is a compact operator from $C[0, 1]$ to $C[0, 1]$. We are interested in an approximate solution of the following integral equation:

$$ \lambda u - Tu = f, \quad (1) $$

where $\lambda$ is not an eigenvalue of $T$ and $f \in C[0, 1]$.

One of the important methods for finding an approximate solution of weakly singular equations is the Product Integration method. Let $r \ge 1$ and let $X_m$ denote the space of piecewise polynomial functions of degree $\le r - 1$ with respect to a partition of $[0, 1]$ with $m$ subintervals. Choose $r$ distinct points in each of the subinterval of the partition and let $\pi_m: C([0, 1]) \to X_m$ be the map which interpolates at these $mr$ points. In the product integration method, (1) is approximated by

$$ \lambda u_m - T\pi_m u_m = f. \quad (2) $$

The solution $u_m$ of the above equation is obtained by solving a linear system of equations of size $mr$. In order to achieve a desired accuracy, it may be necessary to choose $m$ large. We would like to find an approximation of $u_m$ by solving a system of a relatively small size. For this purpose, we approximate $T\pi_m$ by some discrete operators defined on a coarser grid.

In the collocation method, (1) is approximated by $\lambda u_m^C - \pi_m T \pi_m u_m^C = \pi_m f$.

Vainikko and Uba [2] obtain the error bounds in the collocation method. Hakk and Pedas [3] prove that if the interpolation points are chosen appropriately, then the collocation solution exhibits superconvergence at the interpolation points. A two-grid method based on the collocation method is analysed in Hakk and Pedas [4]. Kaneko et al [5] consider approximate solution of weakly singular Hammerstein integral equations by the collocation methods.
In [6], Graham chooses approximate space consisting of spline functions and investigates the Galerkin and the iterated Galerkin methods. Cao and Xu [7] choose approximating space consisting of piecewise polynomials and of mildly singular functions which reflect the singularities of the exact solution. The Galerkin method is then considered with an orthogonal projection onto this approximating space. A disadvantage of this approach is that it is problem specific.

In order to reduce the size of the linear system to solve, Hakk and Pedas analyse a two-grid method based on the collocation method in [4].

The aim of this paper is to investigate the performance of the discrete version of the modified projection method of Kulkarni [8] for approximate solution of (2). In the case of the algebraic singularity, $0 < \alpha < 1$, for an appropriate choice of the coarse grid (value of the grading parameter $q$), we prove that the convergence rate of the discrete modified projection solution is better than those of the discrete collocation and the discrete iterated collocation.

There are not many computational results available for an approximate solution of (1) in research literature and we try to fill this gap. We report some computational results for the weakly singular integral equations with both the algebraic as well as logarithmic singularities. As we considered a graded mesh, a careful computation of the weights in the product integration method was necessary. We followed the approach developed by Atkinson and Shampine in [9]. For the algebraic singularity, numerical results confirm the superiority of the discrete modified projection method over the collocation ones. In the case of the logarithmic singularity, $\alpha = 1$, numerical results show that the discrete iterated modified projection solution converges faster than the other three approaches.

Then, we define 2 two-grid iteration schemes associated with the discrete modified projection method and compare their performance with the two-grid iteration schemes associated with the collocation method. We show that the iterates in the two-grid iteration scheme associated with the discrete modified projection method converge faster than the corresponding iterates in the discrete collocation method.

The paper is arranged as follows. In Section 2 we describe the product integration method, quote a result about the smoothness of the exact solution and state the results about the order of convergence of the approximate solution $u_m$ in the product integration method. In Section 3 the discrete collocation, the discrete iterated collocation, the discrete modified projection and the discrete iterated modified projection methods are defined. Then, we obtain orders of convergence of approximate solutions in various discrete methods. Numerical results related to these discrete methods for the integral operator with both an algebraic singularity and a logarithmic singularity are given in Section 4. Section 5 is devoted to two-grid methods followed by its numerical illustration in Section 6.
4 Fast and Accurate Solvers for Weakly Singular Integral Equations

## 2 Product Integration

In this section we recall the product integration method and state some convergence results which are available in the research literature.

The exact solution of (1) has a singular behavior at the two end points of the interval $[0, 1]$. Hence a graded mesh is needed to define a piecewise polynomial approximation.

The concept of a graded mesh was developed by Rice [10]. In [1] Schneider considers a graded mesh and obtains the orders of convergence for $||u - u_m||_\infty$. He shows that an appropriate choice of interpolation points, such as the Gauss points, improves the order of convergence. Atkinson and Shampine [9] discuss the details of the computation of weights in the product integration method and give computational results.

Let $q \ge 1$, $m \in \mathbb{N}$ and $m$ be even. Consider the following graded partition of $[0, 1]$:

$$ 
\begin{aligned}
y_i &= \frac{1}{2} \left( \frac{2(i-1)}{m} \right)^q, \quad 1 \le i \le \frac{m}{2} + 1, \\
y_i &= 1 - y_{m+2-i}, \quad \frac{m}{2} + 1 \le i \le m + 1.
\end{aligned} 
\quad (3) 
$$

Then

$$ \max_{1 \le i \le m} (y_{i+1} - y_i) \le \frac{q}{m}. \quad (4) $$

Let $r \ge 1$ and let $X_m$ denote the space of piecewise polynomial functions of degree $\le r - 1$ with respect to the above partition. Choose

$$ 0 \le \mu_1 < \mu_2 < \dots < \mu_r \le 1. \quad (5) $$

The mr collocation points are chosen as $s_{i,j} = y_i + \mu_j(y_{i+1} - y_i)$, $i = 1, \dots, m$, $j = 1, \dots, r$.

An interpolatory map $\pi_m : C[0, 1] \to X_m$ is defined by

$$ (\pi_m u)(s_{i,j}) = u(s_{i,j}), \quad i = 1, \dots, m, \quad j = 1, \dots, r. $$

If $\mu_1 = 0$ and $\mu_r = 1$, then $X_m \subset C[0, 1]$ and the $m(r - 1) + 1$ collocation points are given by

$$ s_{i,j} = y_i + \mu_j(y_{i+1} - y_i), \quad i = 1, \dots, m, \quad j = 1, \dots, r-1, s_{m,r} = y_{m+1}. $$

Then $\pi_m$ converges to the identity operator pointwise on $C[0, 1]$. As in Atkinson et al [11], $\pi_m$ can be extended to $L^\infty([0, 1])$ and

$$ \sup_m ||\pi_m|| \le C_1. \quad (6) $$
Let {$\eta_{i,j} : i = 1, \dots, m, j = 1, \dots, r$} be a basis for $X_m$ such that

$$ \eta_{i,j}(s_{k,l}) = \delta_{i,k}\delta_{j,l}. $$

Then

$$ (\pi_m u)(t) = \sum_{i=1}^{m} \sum_{j=1}^{r} u(s_{i,j}) \eta_{i,j}(t), \quad t \in [0, 1]. $$

It follows that

$$ \begin{aligned} T\pi_m u(s) &= \sum_{i=1}^{m} \sum_{j=1}^{r} u(s_{i,j}) \int_{0}^{1} g_{\alpha}(|s-t|) \eta_{i,j}(t) dt \\ &= \sum_{i=1}^{m} \sum_{j=1}^{r} \omega_{i,j}(s) u(s_{i,j}). \end{aligned} \quad (7) $$

The exact equation (1) is approximated by $\lambda u_m - T\pi_m u_m = f$ (equation (2)).

Note that the operators {$T\pi_m$} are collectively compact and converge pointwise to $T$ on $C[0, 1]$. Hence for $m$ large enough, $(\lambda I - T\pi_m)$ is invertible,

$$ \|(\lambda I - T\pi_m)^{-1}\| \le C_2 \quad (8) $$

and (2) has a unique solution. Note that

$$ \begin{aligned} u - u_m &= (\lambda I - T)^{-1}f - (\lambda I - T\pi_m)^{-1}f \\ &= (\lambda I - T\pi_m)^{-1}(T - T\pi_m)(\lambda I - T)^{-1}f \\ &= (\lambda I - T\pi_m)^{-1}T(I - \pi_m)u. \end{aligned} $$

Hence

$$ \|u - u_m\|_{\infty} \le C_2 \|T(I - \pi_m)u\|_{\infty}. \quad (9) $$

We quote a result about the smoothness of the solution $u$ of (1) from Atkinson [12].

The Hölder spaces are defined as follows. For $0 < \alpha < 1$,

$$ C^{(0,\alpha)}[0,1] = \left\{ g \in C[0,1] : d_{\alpha}(g) = \sup_{0 \le t, \tau \le 1} \frac{|g(t) - g(\tau)|}{|t - \tau|^{\alpha}} < \infty \right\} \quad (10) $$

and

$$ C^{(0,1)}[0,1] = \left\{ g \in C[0,1] : \sup_{0 \le t, \tau \le 1} \frac{|g(t) - g(\tau)|}{|t - \tau \log B/(t - \tau)|} < \infty \right\} $$

for some $B > 1$.
6 Fast and Accurate Solvers for Weakly Singular Integral Equations

**Theorem 1** (Atkinson [12]) Let $k \ge 0$ be an integer, and let $0 < \alpha \le 1$. Assume that $f \in C^{(0,\alpha)}[0,1]$, $f \in C^k(0,1)$ and $t^i(1-t)^i f^{(i)}(t) \in C^{(0,\alpha)}[0,1]$, $i=1, \dots, k$. Then

1. The solution $u$ of (1) satisfies $u \in C^{(0,\alpha)}[0,1]$, $u \in C^k(0,1)$ and

$$u_i(t) = t^i(1-t)^i u^{(i)}(t) \in C^{(0,\alpha)}[0,1], \quad i=1, \dots, k.$$

$$\text{Further } u_i(0) = u_i(1) = 0, \quad i=1, \dots, k.$$

2. For $0 < \alpha < 1$, $|u^{(i)}(t)| \le c_i t^{\alpha-i}$, $0 < t \le \frac{1}{2}$, $i=1, \dots, k$,

$$|u^{(i)}(t)| \le c_i (1-t)^{\alpha-i}, \quad \frac{1}{2} \le t < 1, \quad i=1, \dots, k.$$

With $\alpha = 1$, for any $\epsilon \in (0,1)$,

$$|u^{(i)}(t)| \le c_i t^{1-\epsilon-i}, \quad 0 < t \le \frac{1}{2}, \quad i=1, \dots, k,$$

$$|u^{(i)}(t)| \le c_i (1-t)^{1-\epsilon-i}, \quad \frac{1}{2} \le t < 1, \quad i=1, \dots, k.$$

Throughout this paper we assume that $f \in C^{(0,\alpha)}[0,1] \cap C^r(0,1)$ and for $i=1, \dots, r$, $t^i(1-t)^i f^{(i)}(t) \in C^{(0,\alpha)}[0,1]$.

In [10], Rice has proved the following estimate: If $0 < \alpha < 1$ and $q \ge \frac{r}{\alpha}$, then

$$\|u - \pi_m u\|_{\infty} = O(m^{-r}). \qquad (11)$$

It then follows from (9) that

$$\|u - u_m\|_{\infty} = O(m^{-r}). \qquad (12)$$

If $\alpha = 1$ and $q > r$, then

$$\|u - u_m\|_{\infty} = O(m^{-r}). \qquad (13)$$

The above results can be improved by a suitable choice of collocation points.

We describe them below. Let $M = \int_0^1 \prod_{j=1}^r (\mu_j - t)dt$. Note that if $\mu_j$, $j = 1, \dots, r$, are the Gauss points, then $M = 0$. In [1] Schneider has proved the following error estimates.

Let $M = 0$ and $0 < \alpha < 1$. If $q \ge \frac{\alpha + r + 1}{2\alpha}$, then

$$\|T(I - \pi_m)u\|_{\infty} = O(m^{-r-\alpha}). \qquad (14)$$

It then follows from (9) that

$$\|u - u_m\|_{\infty} = O(m^{-r-\alpha}). \qquad (15)$$
Let $M = 0$ and $\alpha = 1$. If $q > \frac{r+2}{2}$, then

$$ \|T(I - \pi_m)u\|_{\infty} = O(m^{-r-1} \log(m)). \qquad (16) $$

It follows that

$$ \|u - u_m\|_{\infty} = O(m^{-r-1} \log(m)). \qquad (17) $$

In the next Section, some discrete operators which approximate $T\pi_m$ are defined and then two-grid methods associated with these discrete operators are described.

## 3 Projection based Discrete Methods

In the case of an integral operator

$$ (Tu)(s) = \int_{0}^{1} k(s, t)u(t)dt, \quad 0 \le s \le 1, $$

where the kernel $k(\cdot, \cdot)$ is either smooth or is of the type of Green's function, the discrete versions of various projection methods are defined by replacing the integral in the definition of $T$ by a convergent quadrature formula.

The discrete version of the iterated collocation method for an Urysohn integral equation with a smooth kernel is considered in Atkinson-Flores [13], whereas the discrete versions of the modified projection and the iterated modified projection methods are investigated in Kulkarni-Rakshit [14].

In the present case, we replace $T$ by

$$ T\pi_m(u)(s) = \sum_{i=1}^{m} \sum_{j=1}^{r} \omega_{i,j}(s)u(s_{i,j}). $$

Choose $n \in \mathbb{N}$ such that $n$ is even and $m = np$, $p \in \mathbb{N}$. Let $q \ge 1$ and consider the following coarse partition of $[0, 1]$:

$$ \begin{aligned} x_i &= \frac{1}{2} \left( \frac{2(i-1)}{n} \right)^q, & 1 \le i \le \frac{n}{2} + 1, \\ x_i &= 1 - x_{n+2-i}, & \frac{n}{2} + 1 < i \le n + 1. \end{aligned} \qquad (18) $$

Let $r \ge 1$ and $X_n$ denote the space of piecewise polynomial functions of degree $\le r - 1$ with respect to the above partition. The $nr$ collocation points are chosen as $t_{i,j} = x_i + \mu_j(x_{i+1} - x_i)$, $i = 1, \dots, n$, $j = 1, \dots, r$, where $\mu_j$'s are from (5).
8 Fast and Accurate Solvers for Weakly Singular Integral Equations

Let $\pi_n : C[0, 1] \to X_n$ be defined as

$$ (\pi_n u)(t_{i,j}) = u(t_{i,j}), \quad i = 1, \dots, n, \quad j = 1, \dots, r. \qquad (19) $$

Note that

$$ y_{ip+1} = \frac{1}{2} \left( \frac{2(ip)}{m} \right)^q = \frac{1}{2} \left( \frac{2(i+1-1)}{n} \right)^q = x_{i+1}. $$

### Discrete Collocation Method

We define $T_n^C = \pi_n T \pi_m \pi_n$ and approximate (2) by

$$ \lambda u_n^C - \pi_n T \pi_m \pi_n u_n^C = \pi_n f. $$

Note that $\pi_n u_n^C = u_n^C$. Hence

$$ (\lambda - \pi_n T \pi_m) u_n^C = \pi_n f. \qquad (20) $$

Furthermore $\pi_n$ converges to the Identity operator on $C[0, 1]$, $\{T\pi_m : m \ge 1\}$ is collectively compact and is pointwise convergent on $C[0, 1]$. Hence $\|(\pi_n - I)T\pi_m\| \to 0$ as $n \to \infty$. It follows that $\{\pi_n T \pi_m : n \ge 1, m \ge 1\}$ is collectively compact. For $n$ large enough, $(\lambda I - \pi_n T \pi_m)$ is invertible,

$$ \|(\lambda I - \pi_n T \pi_m)^{-1}\| \le C_3 \qquad (21) $$

and (20) has a unique solution. Apply $\pi_n$ to (2) to obtain

$$ \lambda \pi_n u_m - \pi_n T \pi_m u_m = \pi_n f. $$

Thus,

$$ (\lambda - \pi_n T \pi_m) u_m = \pi_n f + \lambda (I - \pi_n) u_m. \qquad (22) $$

From (20) and (22) we obtain

$$ u_m - u_n^C = \lambda(\lambda - \pi_n T \pi_m)^{-1}(I - \pi_n)u_m. $$

Hence using the estimates (6) and (21), we obtain

$$ \begin{aligned} \|u_m - u_n^C\|_{\infty} &\le |\lambda| C_3 \|(I - \pi_n)u_m\|_{\infty} \\ &\le |\lambda| C_3 \|(I - \pi_n)u\|_{\infty} + |\lambda| C_3 (1 + C_1) \|u - u_m\|_{\infty}. \end{aligned} \qquad (23) $$

Bearing in mind that $\pi_m \pi_n = \pi_n$ and $T \pi_m \pi_n = T \pi_n$, $u_n^C$ is also the solution of $(\lambda - \pi_n T \pi_n)u_n^C = \pi_n f$, so it is the collocation approximation at the coarse grid.
## Discrete Iterated Collocation Method

We define $T_n^S = T\pi_m\pi_n = T\pi_n$ and approximate (2) by

$$ \lambda u_n^S - T\pi_n u_n^S = f. \qquad (24) $$

Indeed, we recover the product integration on the coarse grid.

Then $u - u_n^S = (\lambda I - T\pi_n)^{-1}T(I - \pi_n)u$. Hence

$$ \|u - u_n^S\|_{\infty} \le C_2 \|T(I - \pi_n)u\|_{\infty} $$

and from (9)

$$ \begin{aligned} \|u_m - u_n^S\|_{\infty} &\le \|u_m - u\|_{\infty} + \|u - u_n^S\|_{\infty} \\ &\le C_2(\|T(I - \pi_m)u\|_{\infty} + \|T(I - \pi_n)u\|_{\infty}). \end{aligned} \qquad (25) $$

## Discrete Modified Projection Method

To build this method, and contrary to the preceding ones, two grids are needed (a coarse and a fine).

We define $T_n^M = \pi_n T \pi_m + T \pi_m \pi_n - \pi_n T \pi_m \pi_n$. Then

$$ T\pi_m - T_n^M = (I - \pi_n)T\pi_m(I - \pi_n) = (I - \pi_n)T(\pi_m - \pi_n). $$

We approximate (2) by

$$ \lambda u_n^M - T_n^M u_n^M = f. \qquad (26) $$

Since $\|T\pi_m - T_n^M\| \le 2C_1 \|(I - \pi_n)T\| \to 0$ as $n \to \infty$, for n large enough, $\lambda I - T_n^M$ is invertible and

$$ \|(\lambda I - T_n^M)^{-1}\| \le C_4. \qquad (27) $$

We have

$$ \begin{aligned} u_m - u_n^M &= (\lambda I - T\pi_m)^{-1}f - (\lambda I - T_n^M)^{-1}f \\ &= (\lambda I - T_n^M)^{-1}(T\pi_m - T_n^M)(\lambda I - T\pi_m)^{-1}f \\ &= (\lambda I - T_n^M)^{-1}(I - \pi_n)T(\pi_m - \pi_n)u_m. \end{aligned} $$

Note that

$$ \|(I - \pi_n)T(\pi_m - \pi_n)u_m\|_{\infty} \le 2C_1 \|(I - \pi_n)T\| \|u_m - u\|_{\infty} + $$
10 Fast and Accurate Solvers for Weakly Singular Integral Equations

$$ \|(I - \pi_n)T(I - \pi_m)u\|_{\infty} + \|(I - \pi_n)T(I - \pi_n)u\|_{\infty}. $$

Hence

$$ \|u_m - u_n^M\|_{\infty} \le 2C_1C_4\|(I - \pi_n)T\|\|u_m - u\|_{\infty} + C_4[\|(I - \pi_n)T(I - \pi_m)u\|_{\infty} + \|(I - \pi_n)T(I - \pi_n)u\|_{\infty}] \quad (28) $$

## Discrete Iterated Modified Projection Method

We define

$$ \tilde{u}_n^M = \frac{1}{\lambda} (T\pi_m u_n^M + f). \quad (29) $$

It can be proved that

$$ \|u_m - \tilde{u}_n^M\|_{\infty} \le C_5(\|T(I - \pi_n)T\|\|u_n^M - u_m\|_{\infty} + \|T(I - \pi_n)T(I - \pi_n)u\|_{\infty}). \quad (30) $$

We skip the details.

## Error Estimates

We consider an integral operator with an algebraic singularity and obtain orders of convergence for the discrete collocation solution, the discrete iterated collocation solution and the discrete modified projection solution.

**Proposition 1** Let $(Tu)(s) = \int_0^1 \frac{1}{|s-t|^{1-\alpha}} u(t) dt, \quad 0 < \alpha < 1. \text{ If } u \in L^\infty[0, 1], \\ \text{then } Tu \in C^{(0,\alpha)}[0, 1] \text{ and}$

$$ d_{\alpha}(Tu) = \sup_{0 \le s, s' \le 1} \frac{(\mathbf{T} u)(s) - (\mathbf{T} u)(s')}{|s - s'|^{\alpha}} \le \frac{2^{2-\alpha} + 1}{\alpha} \|u\|_{\infty}. \quad (31) $$

*Proof* Note that

$$ \begin{aligned} |(Tu)(s) - (Tu)(s')| & \le \left( \int_0^1 \left| \frac{1}{|s-t|^{1-\alpha}} - \frac{1}{|s'-t|^{1-\alpha}} \right| dt \right) \|u\|_{\infty} \\ & \le \left[ \frac{2^{2-\alpha} + 1}{\alpha} \|u\|_{\infty} \right] |s' - s|^{\alpha}. \end{aligned} $$

It follows that $d_{\alpha}(Tu) \le \frac{2^{2-\alpha} + 1}{\alpha} \|u\|_{\infty}$, which completes the proof. $\square$
**Theorem 2** Let $(Tu)(s) = \int_0^1 \frac{1}{|s-t|^{1-\alpha}} u(t) dt$, $0 < \alpha < 1$, and $\pi_n$ be the interpolatory projection defined by (19). Then

$$ \|(I - \pi_n)T\| = O(n^{-\alpha}). \qquad (32) $$

*Proof* Let $g \in C^{(0,\alpha)}[0,1]$. Then $|g(t) - g(s)| \le d_\alpha(g)|t - s|^\alpha$, $s, t \in [0,1]$. We define a piecewise constant function $g_0$ as follows. Let

$$ t_i = \frac{x_i + x_{i+1}}{2}, \quad g_0(x_1) = g(t_1), \quad g_0(t) = g(t_i), \quad t \in (x_i, x_{i+1}], \quad i = 1, \dots, n. $$

Then $\pi_n(g_0) = g_0$. From (4) we recall that $\max_{1 \le i \le n} (x_{i+1} - x_i) \le \frac{q}{n}$. Hence for $t \in [x_i, x_{i+1}]$,

$$ \begin{aligned} |g(t) - (\pi_n g)(t)| &= |(I - \pi_n)(g - g_0)(t)| \\ &\le (1 + \|\pi_n\|)|g(t) - g(t_i)| \\ &\le (1 + C_1)q^\alpha d_\alpha(g)n^{-\alpha}. \end{aligned} $$

Hence if $g \in C^{(0,\alpha)}[0,1]$, then

$$ \|g - \pi_n g\|_{\infty} \le (1 + C_1)q^{\alpha} d_{\alpha}(g)n^{-\alpha}. $$

Let $u \in L^\infty[0,1]$. Then from (31), $d_\alpha(Tu) \le \frac{2^{2-\alpha}+1}{\alpha} \|u\|_\infty$. It follows that

$$ \|(I - \pi_n)Tu\|_{\infty} \le (1 + C_1)q^{\alpha}d_{\alpha}(Tu)n^{-\alpha} \le (1 + C_1)q^{\alpha}\frac{2^{2-\alpha}+1}{\alpha}\|u\|_{\infty}n^{-\alpha}. $$

The required estimate follows by taking the supremum over the set $\{u \in L^\infty[0,1] : \|u\|_\infty \le 1\}$. $\square$

**Theorem 3** If $q \ge \frac{r}{\alpha}$ for $0 < \alpha < 1$ and $q > r$ for $\alpha = 1$, then

$$ \|u_m - u_n^C\|_{\infty} = O(n^{-r}), \quad \|u_m - u_n^S\|_{\infty} = O(n^{-r}). \qquad (33) $$

If $0 < \alpha < 1$ and $q \ge \min\left\{\frac{r}{\alpha}, \frac{\alpha + r + 1}{2\alpha}\right\}$, then

$$ \|u_m - u_n^M\|_{\infty} = O(n^{-r-\alpha}). \qquad (34) $$

If $M = \int_0^1 \prod_{j=1}^r (\mu_j - t)dt = 0$, $0 < \alpha < 1$ and $q \ge \frac{\alpha + r + 1}{2\alpha}$, then

$$ \|u_m - u_n^S\|_{\infty} = O(n^{-r-\alpha}). \qquad (35) $$

If $M = 0$, $\alpha = 1$ and $q > \frac{r+2}{2}$, then

$$ \|u_m - u_n^S\|_{\infty} = O(n^{-r-1} \log(n)). \qquad (36) $$
12 Fast and Accurate Solvers for Weakly Singular Integral Equations

Proof Recall from (23) that

$$
\|u_m - u_n^C\|_{\infty} \le |\lambda|C_3\|(I - \pi_n)u\|_{\infty} + |\lambda|C_3(1 + C_1)\|u - u_m\|_{\infty}.
$$

Since $m \ge n$, from (11) and (12) we obtain $\|u_m - u_n^C\|_\infty = O(n^{-r})$.

Recall from (25) that

$$
\|u_m - u_n^S\|_{\infty} \le C_2(\|T(I - \pi_m)u\|_{\infty} + \|T(I - \pi_n)u\|_{\infty}) \quad (37)
$$

$$
\le C_2 \|T\|(\|(I - \pi_m)u\|_{\infty} + \|(I - \pi_n)u\|_{\infty}). \quad (38)
$$

Hence from (11), $\|u_m - u_n^S\|_\infty = O(n^{-r})$.

We recall from (28) the following error estimate in the discrete modified
projection solution:

$$
\|u_m - u_n^M\|_{\infty} \le 2C_1 C_4 \|(I - \pi_n)T\| \|u_m - u_n\| + \\
C_4 [\|(I - \pi_n)T(I - \pi_m)u\|_{\infty} + \|(I - \pi_n)T(I - \pi_n)u\|_{\infty}]
$$

Using (11), (12), (14) and (32), we conclude that if $0 < \alpha < 1$ and $q \ge$
$\min \left\{ \frac{r}{\alpha}, \frac{\alpha + r + 1}{2\alpha} \right\}$, then

$$
\|u_m - u_n^M\|_{\infty} = O\left(n^{-r-\alpha}\right).
$$

If $M = 0, 0 < \alpha < 1$ and $q \ge \frac{\alpha + r + 1}{2\alpha}$, then from (14) and (37),

$$
\|u_m - u_n^S\|_{\infty} = O(n^{-r-\alpha}).
$$

If $M = 0, \alpha = 1$ and $q > \frac{r+2}{2}$, then from (16) and (37),

$$
\|u_m - u_n^S\|_{\infty} = O(n^{-r-1} \log(n)).
$$

This completes the proof.

<figure><img src="image_2.png" /></figure>

We now deduce an error estimates from Theorem 3 in the case of an
algebraic singularity with $\alpha = \frac{1}{2}$ and $u(t) = \sqrt{t} + \sqrt{1-t}$ as exact solution.

**Corollary 1** Let $\alpha = \frac{1}{2}$ and $u(t) = \sqrt{t} + \sqrt{1-t}$.

1. Let $X_n$ be the space of piecewise linear polynomials, that is, $r=2$, and let $\mu_1 = 0$, $\mu_2 = 1$. Then $M \neq 0$. If $q=4$, then $\|u_m - u_n^C\| = O(n^{-2})$, $\|u_m - u_n^S\| = O(n^{-2.5})$.

2. Let $X_n$ be the space of piecewise constant polynomials, that is, $r = 1$, and
let $\mu_1 = \frac{1}{2}$. Then $M = 0$. If $q \ge 2$, then $\|u_m - u_n^C\| = O(n^{-1})$, $\|u_m - u_n^M\| = O(n^{-1.5})$. If $q \ge 2.5$, then $\|u_m - u_n^S\| = O(n^{-1.5})$.

3. Let $X_n$ be the space of piecewise linear polynomials, that is, $r = 2$, and let $\mu_1, \mu_2$ be the Gauss 2 points. Then $M = 0$. If $q \ge 4$, then $\|u_m - u_n^C\| = O(n^{-2})$. If $q \ge 3.5$, then $\|u_m - u_n^S\| = O(n^{-2.5})$, $\|u_m - u_n^M\| = O(n^{-2.5})$.

We also deduce error estimates for logarithmic singularity ($\alpha = 1$) with the
exact solution chosen to be $u(t) = t \log(t) + (1-t) \log(1-t)$, from Theorem 3.
**Corollary 2** Let $\alpha = 1$ and $u(t) = t \log(t) + (1-t) \log(1-t)$.

1. Let $X_n$ be the space of piecewise linear polynomials, that is, $r = 2$, and let $\mu_1 = 0$, $\mu_2 = 1$. Then $M \neq 0$. If $q > 2$, then $\|u_m - u_n^C\| = O(n^{-2})$, $\|u_m - u_n^S\| = O(n^{-2})$.

2. Let $X_n$ be the space of piecewise constant polynomials, that is, $r = 1$, and let $\mu_1 = \frac{1}{2}$. Then $M = 0$. If $q > 1.5$, then $\|u_m - u_n^C\| = O(n^{-1})$, $\|u_m - u_n^S\| = O(n^{-2} \log(n))$.

3. Let $X_n$ be the space of piecewise linear polynomials, that is, $r = 2$, and let $\mu_1, \mu_2$ be the Gauss 2 points. Then $M = 0$. If $q > 2$, then $\|u_m - u_n^C\| = O(n^{-2})$, $\|u_m - u_n^S\| = O(n^{-3} \log(n))$.

# 4 Numerical Experiments for Projection based Discrete Methods

Recall that $X_m$ is the space of piecewise polynomials with respect to the partition (3) and $\{\eta_{i,j} : i = 1, \dots, m, j = 1, \dots, r\}$ is a basis for $X_m$ such that $\eta_{i,j}(s_{k,l}) = \delta_{i,k}\delta_{j,l}$. Then

$$ (\pi_m u)(t) = \sum_{i=1}^{m} \sum_{j=1}^{r} u(s_{i,j}) \eta_{i,j}(t), \quad t \in [0, 1] $$

and

$$ T\pi_m u(s) = \sum_{i=1}^{m} \sum_{j=1}^{r} u(s_{i,j}) \int_{0}^{1} g_{\alpha}(|s-t|) \eta_{i,j}(t) dt = \sum_{i=1}^{m} \sum_{j=1}^{r} \omega_{i,j}(s) u(s_{i,j}). $$

Consider the case of the piecewise constant polynomials which interpolate at the mid-points. In this case, $r = 1$ and $s_{i,1} = \frac{y_i + y_{i+1}}{2}$.

Note that $\eta_{i,1}(t) = \begin{cases} 1, & \text{if } t \in [y_i, y_{i+1}), \\ 0, & \text{otherwise.} \end{cases}$

Hence

$$ T\pi_m u(s) = \sum_{i=1}^{m} u(s_{i,1}) \int_{y_i}^{y_{i+1}} g_{\alpha}(|s-t|) dt = \sum_{i=1}^{m} \omega_i(s) u(s_{i,1}). $$

In the case of the continuous piecewise linear polynomials which interpolate at the partition points, the dimension of $X_m$ is $m + 1$. The interpolation points are $y_i, i = 1, \dots, m + 1$. We have

$$ (\pi_m u)(t) = \sum_{i=1}^{m+1} u(y_i) \eta_i(t), \quad t \in [0, 1], $$
14 Fast and Accurate Solvers for Weakly Singular Integral Equations

where the basis functions are $\eta_i(t) = \begin{cases} \frac{t-y_{i-1}}{y_i-y_{i-1}}, & \text{if } t \in [y_{i-1}, y_i], \\ \frac{y_{i+1}-t}{y_{i+1}-y_i}, & \text{if } t \in [y_i, y_{i+1}], \\ 0, & \text{otherwise.} \end{cases}$

It follows that

$$ 
\begin{aligned}
T\pi_m u(s) &= \left( \frac{1}{y_2 - y_1} \int_{y_1}^{y_2} g_\alpha(|s-t|)(y_2-t)dt \right) u(y_1) + \\
&\quad \sum_{i=2}^{m} \left( \frac{1}{y_i - y_{i-1}} \int_{y_{i-1}}^{y_i} g_\alpha(|s-t|)(t - y_{i-1})dt + \right. \\
&\qquad \left. \frac{1}{y_{i+1} - y_i} \int_{y_i}^{y_{i+1}} g_\alpha(|s-t|)(y_{i+1} - t)dt \right) u(y_i) + \\
&\quad \left( \frac{1}{y_{m+1} - y_m} \int_{y_m}^{y_{m+1}} g_\alpha(|s-t|)(t - y_m)dt \right) u(y_{m+1}) \\
&= \sum_{i=1}^{m+1} \omega_i(s)u(y_i)
\end{aligned}
$$

In the case of the piecewise linear polynomials which interpolate at Gauss 2 points, the dimension of $X_m$ is $2m$. The interpolation points are given by

$$ s_{i,1} = \frac{y_i + y_{i+1}}{2} - \frac{1}{\sqrt{3}} \frac{y_{i+1} - y_i}{2}, \quad s_{i,2} = \frac{y_i + y_{i+1}}{2} + \frac{1}{\sqrt{3}} \frac{y_{i+1} - y_i}{2}, \quad i = 1, \dots, m. $$

The basis functions are

$$ \eta_{i,1}(t) = \begin{cases} \frac{s_{i,2}-t}{s_{i,2}-s_{i,1}}, & \text{if } t \in [y_i, y_{i+1}], \\ 0, & \text{otherwise,} \end{cases} \quad \eta_{i,2}(t) = \begin{cases} \frac{t-s_{i,1}}{s_{i,2}-s_{i,1}}, & \text{if } t \in [y_i, y_{i+1}], \\ 0, & \text{otherwise.} \end{cases} $$

It follows that

$$ 
\begin{aligned}
T\pi_m u(s) &= \sum_{i=1}^{m} \left( \frac{1}{s_{i,2} - s_{i,1}} \int_{y_i}^{y_{i+1}} g_\alpha(|s-t|)(s_{i,2}-t)dt \right) u(s_{i,1}) + \\
&\quad \sum_{i=1}^{m} \left( \frac{1}{s_{i,2} - s_{i,1}} \int_{y_i}^{y_{i+1}} g_\alpha(|s-t|)(t-s_{i,1})dt \right) u(s_{i,2}) \\
&= \sum_{i=1}^{m} [\omega_{i,1}(s)u(s_{i,1}) + \omega_{i,2}(s)u(s_{i,2})].
\end{aligned}
$$

Thus, in the above three cases, the computation of the weights $\omega_{i,j}(s)$ reduces to the integrals of the form $\int_a^b g_\alpha(|s-t|)dt$ or $\int_a^b g_\alpha(|s-t|)t dt$. These integrals can be evaluated analytically. However, since we are using a graded
mesh, there is a loss of significance. Hence we follow the approach of Atkinson and Shampine [9] and evaluate them numerically with a sufficient accuracy.

### Example 1: Algebraic Singularity

Let

$$T(u)(s) = \int_{0}^{1} \frac{1}{\sqrt{|s-t|}} u(t) dt, \quad \text{that is,} \quad \alpha = \frac{1}{2}.$$

We solve the integral equation

$$10 u - T(u) = f,$$

where the right hand side $f$ is so chosen that $u(t) = \sqrt{t} + \sqrt{1-t}$, $t \in [0, 1]$. This solution satisfies the properties listed in Theorem 1.

We choose $n = 4$ and $m = 256$ (higher values lead to similar results). In the following, we illustrate and confirm the bounds described in Corollary 1. The error and the computed orders of convergence achieved with (i) continuous piecewise linear interpolation at the partition points, (ii) piecewise constant interpolation at the mid-points and (iii) piecewise linear interpolation at Gauss 2 points, are reported. The notation is $u_n^C$ for collocation solution, $u_n^S$ for iterated collocation, $u_n^M$ for modified projection, and $\tilde{u}_n^M$ for iterated modified projection. The integer $q$ is the parameter for the graded mesh and its value has impact on the convergence rate.

<table><caption>**Table 1:** Errors and computed orders of convergence: continuous piecewise linear interpolation at the partition points (Algebraic Singularity).</caption><thead><tr><th>q</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>C</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>S</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u&#771;<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th></tr></thead><tbody><tr><td>1</td><td>3.18 × 10<sup>-02</sup></td><td>0.52</td><td>8.84 × 10<sup>-04</sup></td><td>1.19</td><td>3.98 × 10<sup>-04</sup></td><td>1.11</td><td>1.19 × 10<sup>-05</sup></td><td>1.64</td></tr><tr><td>2</td><td>2.11 × 10<sup>-04</sup></td><td>1.95</td><td>6.77 × 10<sup>-05</sup></td><td>1.98</td><td>1.27 × 10<sup>-05</sup></td><td>2.04</td><td>4.61 × 10<sup>-07</sup></td><td>2.79</td></tr><tr><td>3</td><td>5.33 × 10<sup>-04</sup></td><td>1.58</td><td>8.07 × 10<sup>-05</sup></td><td>2.03</td><td>1.04 × 10<sup>-06</sup></td><td>2.06</td><td>1.38 × 10<sup>-06</sup></td><td>2.31</td></tr><tr><td>4</td><td>7.17 × 10<sup>-04</sup></td><td>1.84</td><td>1.29 × 10<sup>-04</sup></td><td>2.01</td><td>1.70 × 10<sup>-05</sup></td><td>2.34</td><td>2.30 × 10<sup>-06</sup></td><td>2.50</td></tr></tbody></table>

Table 1 shows that for continuous piecewise linear interpolation at the partition points the error in the discrete iterated modified projection method is less than that in the discrete iterated collocation method, as well as less than the two non-iterated versions.

From Corollary 1, if $q = 4$, then the expected orders of convergence are 2 for the collocation and the iterated collocation solutions, whereas 2.5 for the two modified projection solutions. These estimates are confirmed by the numerical computations.

For the piecewise constant interpolation at the mid-points, Table 2, the error in the discrete iterated collocation method is less than the error in the
<table><caption>**Table 2:** Errors and computed orders of convergence: piecewise constant interpolation at the mid-points (Algebraic Singularity).</caption><thead><tr><th>q</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>C</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>S</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u&#773;<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th></tr></thead><tbody><tr><td>1</td><td>4.14 × 10<sup>-02</sup></td><td>0.84</td><td>7.29 × 10<sup>-04</sup></td><td>1.15</td><td>5.55 × 10<sup>-04</sup></td><td>1.17</td><td>9.82 × 10<sup>-06</sup></td><td>1.52</td></tr><tr><td>2</td><td>1.17 × 10<sup>-02</sup></td><td>1.21</td><td>9.23 × 10<sup>-05</sup></td><td>1.68</td><td>6.97 × 10<sup>-05</sup></td><td>1.37</td><td>2.38 × 10<sup>-06</sup></td><td>1.38</td></tr><tr><td>3</td><td>6.84 × 10<sup>-03</sup></td><td>1.27</td><td>1.20 × 10<sup>-04</sup></td><td>1.60</td><td>9.52 × 10<sup>-05</sup></td><td>1.33</td><td>4.15 × 10<sup>-06</sup></td><td>1.26</td></tr><tr><td>4</td><td>7.79 × 10<sup>-03</sup></td><td>1.28</td><td>1.66 × 10<sup>-04</sup></td><td>1.63</td><td>1.30 × 10<sup>-04</sup></td><td>1.36</td><td>6.07 × 10<sup>-06</sup></td><td>1.34</td></tr></tbody></table>

discrete collocation method. As expected the errors in discrete iterated collocation method and the discrete modified projection method are similar. The error in the discrete iterated modified projection method is the least.

From Corollary 1, if $q \ge 2.5$, then the expected orders of convergence for the collocation solution is 1, whereas for the iterated collocation and the modified projection solutions it is 1.5.

<table><caption>**Table 3:** Errors and computed orders of convergence: Piecewise linear interpolation at Gauss 2 points (Algebraic Singularity).</caption><thead><tr><th>q</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>C</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>S</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u&#773;<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th></tr></thead><tbody><tr><td>1</td><td>1.99 × 10<sup>-02</sup></td><td>1.01</td><td>3.18 × 10<sup>-04</sup></td><td>1.27</td><td>2.51 × 10<sup>-04</sup></td><td>1.24</td><td>3.60 × 10<sup>-06</sup></td><td>1.47</td></tr><tr><td>2</td><td>8.75 × 10<sup>-03</sup></td><td>1.33</td><td>3.05 × 10<sup>-05</sup></td><td>1.78</td><td>2.06 × 10<sup>-05</sup></td><td>1.97</td><td>3.49 × 10<sup>-07</sup></td><td>2.36</td></tr><tr><td>3</td><td>2.74 × 10<sup>-03</sup></td><td>1.69</td><td>2.57 × 10<sup>-05</sup></td><td>1.86</td><td>1.75 × 10<sup>-05</sup></td><td>2.49</td><td>7.66 × 10<sup>-07</sup></td><td>2.76</td></tr><tr><td>4</td><td>1.14 × 10<sup>-03</sup></td><td>2.17</td><td>3.99 × 10<sup>-05</sup></td><td>2.10</td><td>3.51 × 10<sup>-05</sup></td><td>2.34</td><td>1.57 × 10<sup>-06</sup></td><td>2.56</td></tr></tbody></table>

Observe that for piecewise linear interpolation at Gauss 2 points (Table 3), again, the error in the discrete iterated collocation method is less than the error in the discrete collocation method. As expected the errors in discrete iterated collocation method and the discrete modified projection method are similar. The error in the discrete iterated modified projection method is the least and there is an improvement in the order of convergence.

The orders of convergence estimated by Corollary 1 are also verified: for $q \ge 3.5$ the expected orders of convergence for the collocation solution is 2, whereas for the modified projection solution it is 2.5.

**Example 2: Logarithmic Singularity**

Let $T(u)(s) = \int_{0}^{1} \log|s - t| u(t) dt$. We solve the integral equation

$$u - T(u) = f,$$

where the right hand side $f$ is so chosen that $u(t) = t \log(t) + (1 - t) \log(1 - t)$, $t \in [0, 1]$. This solution satisfies the properties listed in Theorem 1.
In Table 4 we observe that the discrete modified projection method improves over the discrete iterated collocation method. The discrete iterated modified projection method is still the best.

From Corollary 2, if $q > 2$, then the expected orders of convergence in the discrete collocation method and the discrete iterated collocation method are 2, which are inline with the numerical computations for $q = 4$. Numerical results indicate that the order of convergence for the discrete modified projection methods is 3, as well as for iterated modified projection.

<table><caption>**Table 4:** Errors and computed orders of convergence: continuous piecewise linear interpolation at the partition points (Logarithmic Singularity).</caption><thead><tr><th>q</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>C</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>S</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u&#770;<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th></tr></thead><tbody><tr><td>1</td><td>3.01 × 10<sup>-02</sup></td><td>0.45</td><td>1.49 × 10<sup>-03</sup></td><td>1.40</td><td>3.12 × 10<sup>-04</sup></td><td>1.54</td><td>1.47 × 10<sup>-05</sup></td><td>2.37</td></tr><tr><td>2</td><td>5.43 × 10<sup>-03</sup></td><td>0.97</td><td>9.67 × 10<sup>-05</sup></td><td>2.08</td><td>1.66 × 10<sup>-06</sup></td><td>2.96</td><td>6.81 × 10<sup>-07</sup></td><td>3.08</td></tr><tr><td>3</td><td>8.11 × 10<sup>-04</sup></td><td>1.58</td><td>1.33 × 10<sup>-04</sup></td><td>2.09</td><td>5.85 × 10<sup>-06</sup></td><td>2.92</td><td>2.10 × 10<sup>-06</sup></td><td>3.12</td></tr><tr><td>4</td><td>4.26 × 10<sup>-04</sup></td><td>1.84</td><td>2.24 × 10<sup>-04</sup></td><td>2.09</td><td>1.39 × 10<sup>-05</sup></td><td>2.88</td><td>4.88 × 10<sup>-06</sup></td><td>3.13</td></tr></tbody></table>

For piecewise constant interpolation at the mid-points, Table 5, it can be seen that the error in the discrete iterated collocation method is less than the error in the discrete collocation method. Whatever the value $q$ of the graded mesh can be, the error in the discrete iterated modified projection method is always the smallest.

From Corollary 2, if $q > 1.5$, then the expected order of convergence for the collocation solution is 1, whereas $||u_m - u_n^S||_\infty = O(n^{-2} \log(n))$. Numerical results indicate that the order of convergence in the discrete modified projection method is 2 and in the discrete iterated modified projection method it is 2.5 (even higher).

<table><caption>**Table 5:** Errors and computed orders of convergence: piecewise constant interpolation at the mid-points (Logarithmic Singularity).</caption><thead><tr><th>q</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>C</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>S</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u&#770;<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th></tr></thead><tbody><tr><td>1</td><td>3.13 × 10<sup>-02</sup></td><td>1.07</td><td>6.35 × 10<sup>-04</sup></td><td>1.73</td><td>3.15 × 10<sup>-04</sup></td><td>1.85</td><td>5.71 × 10<sup>-06</sup></td><td>2.69</td></tr><tr><td>2</td><td>1.17 × 10<sup>-02</sup></td><td>1.24</td><td>1.74 × 10<sup>-04</sup></td><td>2.06</td><td>9.07 × 10<sup>-05</sup></td><td>2.09</td><td>2.68 × 10<sup>-06</sup></td><td>2.53</td></tr><tr><td>3</td><td>1.39 × 10<sup>-02</sup></td><td>1.24</td><td>2.91 × 10<sup>-04</sup></td><td>2.06</td><td>1.65 × 10<sup>-04</sup></td><td>2.11</td><td>6.58 × 10<sup>-06</sup></td><td>2.70</td></tr><tr><td>4</td><td>1.69 × 10<sup>-02</sup></td><td>1.25</td><td>4.66 × 10<sup>-04</sup></td><td>2.06</td><td>2.73 × 10<sup>-05</sup></td><td>2.13</td><td>1.35 × 10<sup>-05</sup></td><td>2.98</td></tr></tbody></table>

Considering now piecewise linear interpolation at Gauss 2 points, Table 6, from Corollary 2, if $q > 2$, then the expected order of convergence for the collocation solution is 2, whereas $||u_m - u_n^S||_\infty = O(n^{-3} \log(n))$.

While the performance of the discrete iterated collocation method and the discrete modified projection method is similar, the discrete iterated modified projection method remains the best, both in terms of the error and the order of convergence.
<table><caption>**Table 6:** Errors and computed orders of convergence: piecewise linear interpolation at Gauss 2 points (Logarithmic Singularity).</caption><thead><tr><th>q</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>C</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>S</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th><th>||u<sub>m</sub> - u&#773;<sub>n</sub><sup>M</sup>||<sub>∞</sub></th><th>order</th></tr></thead><tbody><tr><td>1</td><td>6.59 × 10<sup>-03</sup></td><td>1.39</td><td>9.20 × 10<sup>-05</sup></td><td>2.00</td><td>5.90 × 10<sup>-05</sup></td><td>2.06</td><td>1.84 × 10<sup>-06</sup></td><td>2.87</td></tr><tr><td>2</td><td>6.24 × 10<sup>-04</sup></td><td>2.17</td><td>1.16 × 10<sup>-05</sup></td><td>3.03</td><td>1.20 × 10<sup>-05</sup></td><td>3.06</td><td>1.02 × 10<sup>-06</sup></td><td>3.48</td></tr><tr><td>3</td><td>1.04 × 10<sup>-03</sup></td><td>2.27</td><td>3.76 × 10<sup>-05</sup></td><td>2.90</td><td>3.75 × 10<sup>-05</sup></td><td>2.92</td><td>3.67 × 10<sup>-06</sup></td><td>3.45</td></tr><tr><td>4</td><td>1.83 × 10<sup>-03</sup></td><td>2.24</td><td>8.49 × 10<sup>-05</sup></td><td>2.78</td><td>8.30 × 10<sup>-05</sup></td><td>2.80</td><td>9.33 × 10<sup>-06</sup></td><td>3.36</td></tr></tbody></table>

Overall, the discrete iterated modified method is clearly superior to the others, and the discrete modified method can also provide good performance. All theoretical bounds deduced are illustrated by the numerical results.

# 5 Two Grid Methods

In order to reduce the size of the linear systems to be solved, two-grid methods are proposed for approximating the solution of

$$ (\lambda I - T\pi_m)u_m = f. \qquad (39) $$

These methods were considered in Kulkarni [15] for the multivariable Fredholm integral equations with smooth kernels. In [15], T was replaced by a Nyström approximation, whereas in this paper it is replaced by $T\pi_m$. Previously, Hakk and Pedas [4], have proposed a two-grid method for collocation and analyzed the convergence order, but without numerical experiments. The aim of this section is also to fill this gap.

## Method 1

Let $u_m^{(0)}$ be an approximate solution of (39) and $T_n$ be an approximation of $T\pi_m$. The iteration is defined as

$$ \begin{aligned} r^{(k)} &= f - (\lambda I - T\pi_m)u_m^{(k)} = (\lambda I - T\pi_m)(u_m - u_m^{(k)}). \\ u_m^{(k+1)} &= u_m^{(k)} + (\lambda I - T_n)^{-1}r^{(k)}, \quad k = 0, 1, 2, \dots \end{aligned} \qquad (40) $$

Then

$$ \begin{aligned} u_m - u_m^{(k+1)} &= (I - (\lambda I - T_n)^{-1}(\lambda I - T\pi_m))(u_m - u_m^{(k)}) \\ &= (\lambda I - T_n)^{-1}(T\pi_m - T_n)(u_m - u_m^{(k)}). \end{aligned} $$

Let

$$ \tau_n = \sup_{m \ge n} \|(\lambda I - T_n)^{-1}(T\pi_m - T_n)\|. $$

Then

$$ \|u_m - u_m^{(k+1)}\|_{\infty} \le \tau_n \|u_m - u_m^{(k)}\|_{\infty}. $$
Remark 1 Let $T_n = T_n^C = \pi_n T \pi_n$ and consider the iteration scheme

$$r^{(k)} = (\lambda I - T\pi_m)(u_m - u_m^{(k)}).$$

$$u_m^{(k+1)} = u_m^{(k)} + (\lambda I - T_n)^{-1}\pi_n r^{(k)}, \quad k = 0, 1, 2, \dots$$

We show below that $u_m^{(k)}$ is a constant sequence.

Let $u_m^{(0)} = 0$. Then $r^{(0)} = f \Rightarrow \pi_n r^{(0)} = \pi_n f$.

The term $z_n^{(0)} = (\lambda I - T_n)^{-1}\pi_n r^{(0)}$ is obtained by solving

$$\lambda z_n^{(0)} - \pi_n T \pi_n z_n^{(0)} = \pi_n r^{(0)}.$$

Thus, $z_n^{(0)} \in X_n$.

As a consequence, $u_m^{(1)} = z_n^{(0)} \in X_n$. Now since $\pi_m \pi_n = \pi_n$,

$$\begin{aligned}
\pi_n r^{(1)} &= \pi_n f - \lambda \pi_n u_m^{(1)} + \pi_n T \pi_m \pi_n u_m^{(1)} \\
&= \pi_n r^{(0)} - \lambda z_n^{(0)} + \pi_n T \pi_n z_n^{(0)} \\
&= 0.
\end{aligned}$$

Hence $z_n^{(1)} = (\lambda I - T_n)^{-1}\pi_n r^{(1)} = 0$ and

$$u_m^{(2)} = u_m^{(1)} = z_n^{(0)}.$$

Thus there is no improvement by iteration.

Remark 2 If we choose $T_n = T_n^M = \pi_n T \pi_m + T \pi_m \pi_n - \pi_n T \pi_m \pi_n$, then

$$T\pi_m - T_n = (I - \pi_n)T\pi_m (I - \pi_n) = (I - \pi_n)T(\pi_m - \pi_n).$$

Using (6) and (27) we deduce that

$$\begin{aligned}
\tau_n &\le \|(\lambda I - T_n^M)^{-1}\| \sup_{m \ge n} (\|\pi_m\| + \|\pi_n\|) \|(I - \pi_n)T\| \\
&\le 2C_1 C_4 \|(I - \pi_n)T\|.
\end{aligned}$$

It follows that

$$\tau_n = O(\|(I - \pi_n)T\|). \qquad (41)$$

Since $T$ is a compact operator and $\pi_n$ converges to the Identity operator $I$ pointwise, it follows that $\|(I - \pi_n)T\| \to 0$ as $n \to \infty$. Hence for $n$ large enough $\tau_n < 1$ and

$$\|u_m - u_m^{(k+1)}\|_{\infty} \le \tau_n \|u_m - u_m^{(k)}\|_{\infty} \qquad (42)$$

## Method 2

Let $\tilde{u}_m^{(0)}$ be an approximate solution

$$(\lambda I - T\pi_m)u_m = f.$$

Let $T_n$ be an approximation of $T\pi_m$. The iteration scheme is defined as

$$r^{(k)} = f - (\lambda I - T\pi_m)\tilde{u}_m^{(k)} = (\lambda I - T\pi_m)(u_m - \tilde{u}_m^{(k)})$$
20 Fast and Accurate Solvers for Weakly Singular Integral Equations

$$
\tilde{u}_m^{(k+1)} = \tilde{u}_m^{(k)} + \frac{1}{\lambda} \left[ r^{(k)} + (\lambda I - T_n)^{-1} T \pi_m r^{(k)} \right], \quad k = 0, 1, \dots \quad (43)
$$

In this case

$$
\begin{align*}
u_m - \tilde{u}_m^{(k+1)} &= u_m - \tilde{u}_m^{(k)} - \frac{1}{\lambda} [I + (\lambda I - T_n)^{-1} T \pi_m] r^{(k)} \\
&= -\frac{1}{\lambda} (\lambda I - T_n)^{-1} (T_n - T \pi_m) T \pi_m (u_m - \tilde{u}_m^{(k)}).
\end{align*}
$$

Let

$$
\tilde{\tau}_n = \frac{1}{|\lambda|} \sup_{m \ge n} \|(\lambda I - T_n)^{-1} (T_n - T\pi_m) T\pi_m\|.
$$

Then

$$
\|u_m - \tilde{u}_m^{(k+1)}\|_{\infty} \le \tilde{\tau}_n \|u_m - \tilde{u}_m^{(k)}\|_{\infty} \quad (44)
$$

Remark 3 If we choose $T_n = T_n^C = \pi_n T \pi_n$, then

$$
T\pi_m - T_n = (I - \pi_n)T\pi_m + \pi_n T(\pi_m - I) + \pi_n T(I - \pi_n).
$$

Thus,

$$
\begin{align*}
\|(T\pi_m - T_n)T\pi_m\| &\le \|(I - \pi_n)T\|\|T\|\|\pi_m\|^2 + \|(I - \pi_m)T\|\|T\|\|\pi_n\|\|\pi_m\| \\
&\quad + \|(I - \pi_n)T\|\|T\|\|\pi_n\|\|\pi_m\| \\
&\le 3(C_1)^2 \|T\|\|(I - \pi_n)T\|.
\end{align*}
$$

It follows that

$$
\tilde{\tau}_n = O(\|(I - \pi_n)T\|). \tag{45}
$$

Remark 4 If we choose $T_n = T_n^M = \pi_n T \pi_m + T \pi_m \pi_n - \pi_n T \pi_m \pi_n$, then

$$
\begin{align*}
\|(T\pi_m - T_n)T\pi_m\| &\le \|(I - \pi_n)T(I - \pi_m)T\|\|\pi_m\| + \\
&\quad \|(I - \pi_n)T(I - \pi_n)T\|\|\pi_m\| \\
&\le 2C_1 \|(I - \pi_n)T\|^2.
\end{align*}
$$

Hence

$$
\tilde{\tau}_n = O \left( \|(I - \pi_n)T\|^2 \right). \qquad (46)
$$

A comparison of (41), (45) and (46) shows that while the iterates obtained
in the two-grid method 2 with the discrete collocation approximation and
the iterates obtained in the two-grid method 1 with the discrete modified
projection approximation converge at the same rate, the iterates obtained in
the two-grid method 2 with the discrete modified projection approximation
converge faster than these two iterates. These results are validated in Section
5.3 for specific examples of weakly singular integral equations with an algebraic
or a logarithmic singularity.

The two-grid method 2 with $T_n = T_n^C$ is considered in Hakk and Pedas [4].
However, there are no computational results.
# 6 Numerical Experiments for Two Grid Methods

We now give numerical results for the two grid methods. The results show that the two-grid method 2 based on the discrete modified projection method needs about half the number of iterations to achieve the desired accuracy as compared to the two-grid method 2 based on the discrete collocation method and the two-grid method 1 based on the discrete modified projection method. Thus, the results given in (41), (42), (44), (45) and (46) are validated.

## Algebraic Singularity

We consider the integral equation $10 u - \int_0^1 \frac{1}{\sqrt{|s-t|}} u(t) dt = f$, where the right hand side $f$ is so chosen that the exact solution is $u(t) = \sqrt{t} + \sqrt{1-t}$, $t \in [0, 1]$. We choose $m = 256$, $n = 4$.

The stopping criterium used for the iterative process was

$$ \frac{\|u_m - u_m^{(k)}\|}{\|u_m\|} < 10 \times \text{eps} $$

for method 1, and

$$ \frac{\|u_m - \tilde{u}_m^{(k)}\|}{\|\tilde{u}_m\|} < 10 \times \text{eps} $$

for method 2, where $\text{eps} \approx 2.2 \times 10^{-16}$ is the machine precision.

Tables 7 to 9 present the errors for the three interpolation schemes. As predicted by (46) modified projection method 2 delivers the fastest convergence, requiring half of the number of iterations as the others. Nevertheless, all the three schemes allowed to reach the desired accuracy. The parameter $q$ controlling the graded mesh does not seem to influence the speed of convergence (just for illustration, results with $q \neq 1$ are shown together with the ones obtained with the uniform grid $q = 1$).

## Logarithmic Singularity

We solve the integral equation $u(s) - \int_0^1 \log|s-t| u(t) dt = f(s)$, where the right hand side $f$ is so chosen that $u(t) = t \log(t) + (1-t) \log(1-t)$, $t \in [0, 1]$. Again we choose $m = 256$, $n = 4$, and the stopping criterium is established by the relative error

$$ \frac{\|u_m - u_m^k\|}{\|u_m\|} < 10 \times \text{eps} $$
<table><caption>**Table 7:** Errors for continuous piecewise linear interpolation at the partition points (Algebraic Singularity)</caption><thead><tr><th rowspan="2">it.</th><th>Collocation<br>method 2</th><th>Modified Proj.<br>method 1</th><th>Modified Proj.<br>method 2</th><th>Collocation<br>method 2</th><th>Modified Proj.<br>method 1</th><th>Modified Proj.<br>method 2</th></tr><tr><th colspan="3">q = 1</th><th colspan="3">q = 4</th></tr></thead><tbody><tr><td>1</td><td>3.47 × 10<sup>-02</sup></td><td>7.16 × 10<sup>-03</sup></td><td>1.81 × 10<sup>-03</sup></td><td>1.50 × 10<sup>-02</sup></td><td>6.94 × 10<sup>-03</sup></td><td>2.26 × 10<sup>-03</sup></td></tr><tr><td>2</td><td>1.90 × 10<sup>-03</sup></td><td>4.03 × 10<sup>-04</sup></td><td>6.24 × 10<sup>-06</sup></td><td>8.16 × 10<sup>-04</sup></td><td>4.90 × 10<sup>-04</sup></td><td>1.25 × 10<sup>-05</sup></td></tr><tr><td>3</td><td>1.11 × 10<sup>-04</sup></td><td>2.37 × 10<sup>-05</sup></td><td>2.34 × 10<sup>-08</sup></td><td>5.55 × 10<sup>-05</sup></td><td>3.52 × 10<sup>-05</sup></td><td>7.10 × 10<sup>-08</sup></td></tr><tr><td>4</td><td>6.74 × 10<sup>-06</sup></td><td>1.43 × 10<sup>-06</sup></td><td>9.09 × 10<sup>-11</sup></td><td>3.87 × 10<sup>-06</sup></td><td>2.55 × 10<sup>-06</sup></td><td>4.04 × 10<sup>-10</sup></td></tr><tr><td>5</td><td>4.15 × 10<sup>-07</sup></td><td>8.73 × 10<sup>-08</sup></td><td>3.63 × 10<sup>-13</sup></td><td>2.72 × 10<sup>-07</sup></td><td>1.85 × 10<sup>-07</sup></td><td>2.30 × 10<sup>-12</sup></td></tr><tr><td>6</td><td>2.58 × 10<sup>-08</sup></td><td>5.39 × 10<sup>-09</sup></td><td>3.11 × 10<sup>-15</sup></td><td>1.92 × 10<sup>-08</sup></td><td>1.35 × 10<sup>-08</sup></td><td>1.47 × 10<sup>-14</sup></td></tr><tr><td>7</td><td>1.62 × 10<sup>-09</sup></td><td>3.36 × 10<sup>-10</sup></td><td>2.66 × 10<sup>-15</sup></td><td>1.36 × 10<sup>-09</sup></td><td>9.76 × 10<sup>-10</sup></td><td>2.66 × 10<sup>-15</sup></td></tr><tr><td>8</td><td>1.02 × 10<sup>-09</sup></td><td>2.10 × 10<sup>-11</sup></td><td>2.44 × 10<sup>-15</sup></td><td>9.66 × 10<sup>-11</sup></td><td>7.08 × 10<sup>-11</sup></td><td></td></tr><tr><td>9</td><td>6.48 × 10<sup>-12</sup></td><td>1.33 × 10<sup>-12</sup></td><td></td><td>6.86 × 10<sup>-12</sup></td><td>5.13 × 10<sup>-12</sup></td><td></td></tr><tr><td>10</td><td>4.13 × 10<sup>-13</sup></td><td>8.44 × 10<sup>-14</sup></td><td></td><td>4.89 × 10<sup>-13</sup></td><td>3.72 × 10<sup>-13</sup></td><td></td></tr><tr><td>11</td><td>2.74 × 10<sup>-14</sup></td><td>6.44 × 10<sup>-15</sup></td><td></td><td>3.62 × 10<sup>-14</sup></td><td>2.82 × 10<sup>-14</sup></td><td></td></tr><tr><td>12</td><td>3.33 × 10<sup>-15</sup></td><td>2.44 × 10<sup>-15</sup></td><td></td><td>4.44 × 10<sup>-15</sup></td><td>4.00 × 10<sup>-15</sup></td><td></td></tr><tr><td>13</td><td>2.44 × 10<sup>-15</sup></td><td></td><td></td><td>2.66 × 10<sup>-15</sup></td><td>2.66 × 10<sup>-15</sup></td><td></td></tr></tbody></table>

<table><caption>**Table 8:** Errors for piecewise constant interpolation at the mid-points (Algebraic Singularity)</caption><thead><tr><th rowspan="2">it.</th><th>Collocation<br>method 2</th><th>Modified Proj.<br>method 1</th><th>Modified Proj.<br>method 2</th><th>Collocation<br>method 2</th><th>Modified Proj.<br>method 1</th><th>Modified Proj.<br>method 2</th></tr><tr><th colspan="3">q = 1</th><th colspan="3">q = 3</th></tr></thead><tbody><tr><td>1</td><td>7.40 × 10<sup>-02</sup></td><td>8.65 × 10<sup>-03</sup></td><td>2.71 × 10<sup>-03</sup></td><td>5.74 × 10<sup>-02</sup></td><td>8.02 × 10<sup>-03</sup></td><td>2.91 × 10<sup>-03</sup></td></tr><tr><td>2</td><td>2.68 × 10<sup>-03</sup></td><td>4.17 × 10<sup>-04</sup></td><td>7.27 × 10<sup>-06</sup></td><td>2.88 × 10<sup>-03</sup></td><td>5.35 × 10<sup>-04</sup></td><td>1.46 × 10<sup>-05</sup></td></tr><tr><td>3</td><td>1.32 × 10<sup>-04</sup></td><td>2.24 × 10<sup>-05</sup></td><td>2.42 × 10<sup>-08</sup></td><td>1.94 × 10<sup>-04</sup></td><td>3.94 × 10<sup>-05</sup></td><td>1.22 × 10<sup>-07</sup></td></tr><tr><td>4</td><td>7.19 × 10<sup>-06</sup></td><td>1.26 × 10<sup>-06</sup></td><td>9.44 × 10<sup>-11</sup></td><td>1.44 × 10<sup>-05</sup></td><td>3.47 × 10<sup>-06</sup></td><td>1.01 × 10<sup>-09</sup></td></tr><tr><td>5</td><td>4.10 × 10<sup>-07</sup></td><td>7.33 × 10<sup>-08</sup></td><td>3.94 × 10<sup>-13</sup></td><td>1.32 × 10<sup>-06</sup></td><td>3.21 × 10<sup>-07</sup></td><td>8.19 × 10<sup>-12</sup></td></tr><tr><td>6</td><td>2.40 × 10<sup>-08</sup></td><td>4.39 × 10<sup>-09</sup></td><td>3.11 × 10<sup>-15</sup></td><td>1.22 × 10<sup>-07</sup></td><td>2.94 × 10<sup>-08</sup></td><td>6.59 × 10<sup>-14</sup></td></tr><tr><td>7</td><td>1.47 × 10<sup>-09</sup></td><td>2.83 × 10<sup>-10</sup></td><td>2.22 × 10<sup>-15</sup></td><td>1.11 × 10<sup>-08</sup></td><td>2.66 × 10<sup>-09</sup></td><td>3.33 × 10<sup>-15</sup></td></tr><tr><td>8</td><td>9.46 × 10<sup>-11</sup></td><td>1.82 × 10<sup>-11</sup></td><td></td><td>1.01 × 10<sup>-09</sup></td><td>2.40 × 10<sup>-10</sup></td><td>2.66 × 10<sup>-15</sup></td></tr><tr><td>9</td><td>6.11 × 10<sup>-12</sup></td><td>1.18 × 10<sup>-12</sup></td><td></td><td>9.09 × 10<sup>-11</sup></td><td>2.15 × 10<sup>-11</sup></td><td></td></tr><tr><td>10</td><td>3.96 × 10<sup>-13</sup></td><td>7.75 × 10<sup>-14</sup></td><td></td><td>8.19 × 10<sup>-08</sup></td><td>1.93 × 10<sup>-12</sup></td><td></td></tr><tr><td>11</td><td>2.64 × 10<sup>-14</sup></td><td>6.44 × 10<sup>-15</sup></td><td></td><td>7.35 × 10<sup>-13</sup></td><td>1.72 × 10<sup>-13</sup></td><td></td></tr><tr><td>12</td><td>3.11 × 10<sup>-15</sup></td><td>2.22 × 10<sup>-15</sup></td><td></td><td>6.59 × 10<sup>-14</sup></td><td>1.55 × 10<sup>-14</sup></td><td></td></tr><tr><td>13</td><td>2.22 × 10<sup>-15</sup></td><td></td><td></td><td>6.88 × 10<sup>-15</sup></td><td>3.33 × 10<sup>-15</sup></td><td></td></tr><tr><td>14</td><td></td><td></td><td></td><td>3.11 × 10<sup>-15</sup></td><td>2.89 × 10<sup>-15</sup></td><td></td></tr><tr><td>15</td><td></td><td></td><td></td><td>2.66 × 10<sup>-15</sup></td><td>2.66 × 10<sup>-15</sup></td><td></td></tr></tbody></table>

for method 1, and

$$ \frac{\|u_m - \tilde{u}_m^k\|}{\|\tilde{u}_m\|} < 10 \times \text{eps} $$

for method 2. To lighten the presentation, only piecewise linear interpolation at Gauss 2 points on a uniform grid are considered. For the other alternatives, the outputs are similar.

Figure 1 shows the relative error against the number of iterations to compute the approximate solution for the three different iterative approaches (Collocation - method 2, Modified Projection - methods 1 and 2) on a semi-logarithmic plot.
<table><caption>**Table 9:** Errors for piecewise linear interpolation at Gauss 2 points (Algebraic Singularity)</caption><thead><tr><th rowspan="2">it.</th><th>Collocation<br>method 2</th><th>Modified Proj.<br>method 1</th><th>Modified Proj.<br>method 2</th><th>Collocation<br>method 2</th><th>Modified Proj.<br>method 1</th><th>Modified Proj.<br>method 2</th></tr><tr><th colspan="3">q = 1</th><th colspan="3">q = 4</th></tr></thead><tbody><tr><td>1</td><td>4.98 × 10<sup>-02</sup></td><td>5.42 × 10<sup>-03</sup></td><td>1.55 × 10<sup>-03</sup></td><td>5.82 × 10<sup>-02</sup></td><td>5.43 × 10<sup>-03</sup></td><td>1.55 × 10<sup>-03</sup></td></tr><tr><td>2</td><td>1.53 × 10<sup>-03</sup></td><td>2.32 × 10<sup>-04</sup></td><td>4.71 × 10<sup>-06</sup></td><td>1.53 × 10<sup>-03</sup></td><td>2.32 × 10<sup>-04</sup></td><td>4.71 × 10<sup>-06</sup></td></tr><tr><td>3</td><td>6.86 × 10<sup>-05</sup></td><td>1.41 × 10<sup>-05</sup></td><td>2.52 × 10<sup>-08</sup></td><td>6.87 × 10<sup>-05</sup></td><td>1.41 × 10<sup>-05</sup></td><td>2.52 × 10<sup>-08</sup></td></tr><tr><td>4</td><td>4.71 × 10<sup>-06</sup></td><td>1.04 × 10<sup>-06</sup></td><td>1.30 × 10<sup>-10</sup></td><td>4.70 × 10<sup>-06</sup></td><td>1.04 × 10<sup>-06</sup></td><td>1.30 × 10<sup>-10</sup></td></tr><tr><td>5</td><td>3.47 × 10<sup>-07</sup></td><td>7.54 × 10<sup>-08</sup></td><td>6.57 × 10<sup>-13</sup></td><td>3.47 × 10<sup>-07</sup></td><td>7.54 × 10<sup>-08</sup></td><td>6.56 × 10<sup>-13</sup></td></tr><tr><td>6</td><td>2.53 × 10<sup>-08</sup></td><td>5.42 × 10<sup>-09</sup></td><td>5.33 × 10<sup>-15</sup></td><td>2.53 × 10<sup>-08</sup></td><td>5.41 × 10<sup>-09</sup></td><td>5.11 × 10<sup>-15</sup></td></tr><tr><td>7</td><td>1.82 × 10<sup>-09</sup></td><td>3.86 × 10<sup>-10</sup></td><td>2.89 × 10<sup>-15</sup></td><td>1.82 × 10<sup>-09</sup></td><td>3.86 × 10<sup>-10</sup></td><td>3.11 × 10<sup>-15</sup></td></tr><tr><td>8</td><td>1.30 × 10<sup>-10</sup></td><td>2.75 × 10<sup>-11</sup></td><td>2.66 × 10<sup>-15</sup></td><td>1.31 × 10<sup>-10</sup></td><td>2.75 × 10<sup>-11</sup></td><td></td></tr><tr><td>9</td><td>9.32 × 10<sup>-12</sup></td><td>1.95 × 10<sup>-12</sup></td><td></td><td>9.33 × 10<sup>-12</sup></td><td>1.95 × 10<sup>-12</sup></td><td></td></tr><tr><td>10</td><td>6.66 × 10<sup>-13</sup></td><td>1.39 × 10<sup>-13</sup></td><td></td><td>6.55 × 10<sup>-13</sup></td><td>1.13 × 10<sup>-13</sup></td><td></td></tr><tr><td>11</td><td>4.91 × 10<sup>-14</sup></td><td>1.13 × 10<sup>-14</sup></td><td></td><td>4.80 × 10<sup>-14</sup></td><td>1.13 × 10<sup>-14</sup></td><td></td></tr><tr><td>12</td><td>5.33 × 10<sup>-15</sup></td><td>2.89 × 10<sup>-15</sup></td><td></td><td>5.33 × 10<sup>-15</sup></td><td>3.11 × 10<sup>-15</sup></td><td></td></tr><tr><td>13</td><td>2.89 × 10<sup>-15</sup></td><td>2.66 × 10<sup>-15</sup></td><td></td><td>2.89 × 10<sup>-15</sup></td><td></td><td></td></tr><tr><td>14</td><td>2.66 × 10<sup>-15</sup></td><td></td><td></td><td>2.53 × 10<sup>-15</sup></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 1:** Two grid convergence history for piecewise linear interpolation at Gauss 2 points (Logarithmic Singularity).</figcaption></figure>

The number of iterates in the two-grid method associated with the modified projection method 2 is about half as compared with those in the collocation method 2 and modified projection method 1. This is inline with the theoretical results given in (41), (45), and (46).

We now consider a larger value for *m* to report on elapsed time. Figure 2 shows the required number of iterations and time in seconds for *m* = 1024 in terms of increasing values of *n*. Overall, computation time is greatly reduced for the 2-grid approach when compared with the time to compute the product integration on the fine grid. Just for note, we mention that for *m* = 1024, 3.50
<figure><img src="image_4.png" /><figcaption>**Fig. 2:** Two grid required number of iterations and time in seconds for $m = 1024$ and several values of $n$ (logarithmic scale), for piecewise linear interpolation at Gauss 2 points (Logarithmic Singularity).</figcaption></figure>

seconds are required while only 0.03 seconds for the two-grid approach with $n = 4$. That is, more than 100 times faster. The elapsed time increases for increasing values of $n$ but the number of iterations decreases. An increase of $2^2$ in the coarse grid induces a reduction in the number of iterations in $2^{-1}$. The modified iteration method 2 consistently requires half of the iterations to reach the desired accuracy when compared with the other two approaches, as already reported.

It is worth to mention that Atkinson, in the nineties in [12], stated that “In our experience these [two-grid methods] are the most efficient numerical methods for solving the linear systems obtained when solving integral equations of the second kind”. All our results clearly highlight Atkinson’s statement.

# 7 Conclusion

In this paper we consider approximate solution of a weakly singular Fredholm integral equations. The kernel of the integral operator is of the form $|s-t|^{\alpha-1}$, $0 < \alpha < 1$ or of the form $\log(|s-t|)$, the case corresponding to $\alpha = 1$. The approximating space is chosen to be the space of piecewise polynomials with respect to a graded partition of $[0, 1]$ with grading parameter $q$. We consider a fine partition of $[0, 1]$ with $m$ subintervals and a coarse partition with $n$ subintervals. We choose $r$ distinct points $0 \le \mu_1 < \mu_2 < \cdots < \mu_r \le 1$ and the interpolation points are obtained by shifting these points to each subinterval of the partition. The value $M = \int_0^1 \prod_{j=1}^r (\mu_j - t)dt$ characterizes the grid type:
if $M = 0$, Gauss points are used. The interpolatory projection associated with the fine partition is denoted by $\pi_m$. In the Product Integration method, $T$ is replaced by $T\pi_m$ and the solution is denoted by $u_m$.

Our aim is to obtain approximations of $u_m$. We prove the following results:

1. If $0 < \alpha < 1$, $M \neq 0$ and $q \ge \frac{r}{\alpha}$, then the discrete collocation and the iterated discrete collocation solution converge to $u_m$ at the same rate, whereas in the discrete modified projection method, the order of convergence is improved by $\alpha$.

2. If $0 < \alpha < 1$, $M = 0$ and $q \ge \frac{\alpha + r + 1}{2\alpha}$, then the iterated discrete collocation solution converge faster than the discrete collocation solution. Both the discrete iterated collocation and the discrete modified projection solutions converge to $u_m$ at the same rate.

3. If $\alpha = 1$, $M \neq 0$, $q > r$, then the collocation solution and the iterated collocation solution converge at the same rate. Numerical results indicate that in the discrete modified projection method, the order of convergence is improved by 1.

4. If $\alpha = 1$, $M = 0$, $q > \frac{r+2}{2}$, then the iterated discrete collocation solution converges faster than the discrete collocation solution. Numerical results indicate that the error in the discrete iterated modified projection method is the least and there is an improvement in the order of convergence as compared to the discrete iterated collocation solution.

5. In the case of Gauss points ($M = 0$), we expect that the discrete iterated modified projection solution would converge faster than the discrete modified projection solution. However, we could not prove it theoretically. If $0 < \alpha \le 1$, then the numerical results show that the error is the least in the discrete iterated modified projection method. If $\alpha = 1$, then numerical results show, in addition, an improvement in the order of convergence of the discrete iterated modified projection solution.

6. In the two-grid method 2 based on the discrete collocation method, we have the following estimate:

$$ \|u_m - \tilde{u}_m^{(k+1)}\|_{\infty} \le M_1 \|(I - \pi_n)T\| \|u_m - \tilde{u}_m^{(k)}\|_{\infty}. $$

On the other hand, in the two-grid method 2 based on the discrete modified projection method, we obtain

$$ \|u_m - \tilde{u}_m^{(k+1)}\|_{\infty} \le M_2 \|(I - \pi_n)T\|^2 \|u_m - \tilde{u}_m^{(k)}\|_{\infty}. $$

Numerical results show that the number of iterations in the two-grid method associated with the modified projection method 2 is half as compared with those in the discrete collocation method 2. Both approaches favor the computational efficiency by reducing dramatically the elapsed time to reach a prescribed accuracy (near machine precision).
26 Fast and Accurate Solvers for Weakly Singular Integral Equations

The results in this paper can be extended to an integral operator of the form

$$ (Tu)(s) = \int_{0}^{1} k(s,t)g_{\alpha}(|s-t|)u(t)dt, \quad 0 \le s \le 1, $$

where $k(s,t)$ is a smooth function and

$$ g_{\alpha}(s) = s^{\alpha-1}, \quad 0 < \alpha < 1 \text{ and } g_{1}(s) = \log(s). $$

In this case the formula (7) will be replaced by

$$ \begin{aligned} T_m u(s) &= \sum_{i=1}^{m} \sum_{j=1}^{r} k(s, s_{i,j}) u(s_{i,j}) \int_{0}^{1} g_{\alpha}(|s-t|) \eta_{i,j}(t) dt \\ &= \sum_{i=1}^{m} \sum_{j=1}^{r} \omega_{i,j}(s) k(s, s_{i,j}) u(s_{i,j}). \end{aligned} $$

**Acknowledgments.** The third author was partially supported by CMUP, which is financed by national funds through FCT - Fundação para a Ciência e Tecnologia, I.P., under the project with reference UIDB/00144/2020.

## Declarations

* Conflict of interest/Competing interests: The authors declare that they have no conflict of interest.

* Availability of data and materials: Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

* Code availability: 'Not applicable'

## References

[1] Schneider, C.: Product integration for weakly singular integral equations. Math Comput **36**(153), 207–213 (1981). https://doi.org/10.2307/2007737

[2] Vainikko, G., Uba, P.: A piecewise polynomial approximation to the solution of an integral equation with weakly singular kernel. Anziam J **22**(4), 431–438 (1981). https://doi.org/10.1017/S0334270000002770

[3] Hakk, K., Pedas, A.: Numerical solutions and their superconvergence for weakly singular integral equations with discontinuous coefficients. Math Model Anal **3**(1), 104–113 (1998). https://doi.org/10.3846/13926292.1998.9637093

[4] Hakk, K., Pedas, A.: Two-grid iteration method for weakly singular integral equations. Math Model Anal **5**(1), 76–85 (2000). https://doi.org/10.1080/13926292.2000.9637130
[5] Kaneko, H., Noren, R.D., Xu, Y.: Numerical solutions for weakly singular hammerstein equations and their superconvergence. J Integral Equ Appl 4(3), 391–407 (1992)

[6] Graham, I.G.: Galerkin methods for second kind integral equations with singularities. Math Comput 39(160), 519–533 (1982). https://doi.org/10.1090/S0025-5718-1982-0669644-3

[7] Cao, Y., Xu, Y.: Singularity preserving galerkin methods for weakly singular fredholm integral equations. J Integral Equ Appl, 303–334 (1994). https://doi.org/10.1216/jiea/1181075816

[8] Kulkarni, R.P.: A superconvergence result for solutions of compact operator equations. B Aust Math Soc 68(3), 517–528 (2003). https://doi.org/10.1017/S0004972700037916

[9] Atkinson, K.E., Shampine, L.F.: Algorithm 876: Solving fredholm integral equations of the second kind in matlab. ACM T Math Software 34(4), 1–20 (2008). https://doi.org/10.1145/1377596.1377601

[10] Rice, J.R.: On the degree of convergence of nonlinear spline approximation. In: Approximations with Special Emphasis on Spline Functions, pp. 349–365 (1969). Academic Press New York

[11] Atkinson, K., Graham, I., Sloan, I.: Piecewise continuous collocation for integral equations. SIAM J Numer Anal 20(1), 172–186 (1983). https://doi.org/10.1137/0720012

[12] Atkinson, K.E.: The Numerical Solution of Integral Equations of the Second Kind. Cambridge, Cambridge Monographs on Applied and Computational Mathematics (1997). https://doi.org/10.1017/CBO9780511626340

[13] Atkinson, K., Flores, J.: The discrete collocation method for nonlinear integral equations. IMA J Numer Anal 13(2), 195–213 (1993). https://doi.org/10.1093/imanum/13.2.195

[14] Kulkarni, R.P., Rakshit, G.: Discrete modified projection method for urysohn integral equations with smooth kernels. Appl Numer Math 126, 180–198 (2018). https://doi.org/10.1016/j.apnum.2017.12.008

[15] Kulkarni, R.P.: Approximate solution of multivariable integral equations of the second kind. J Integral Equ Appl, 343–374 (2004). https://doi.org/10.1216/jiea/1181075296