PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 96 (110) (2014), 143–157 DOI: 10.2298/PIM1410143F

A FAST ALGORITHM FOR THE NUMERICAL SOLUTION OF AN INTEGRAL EQUATION WITH LOGARITHMIC KERNEL

Katharina Flemming and Peter Junghanns

Abstract. We describe an algorithm for the numerical solution of an integral equation of the form

− \frac{1}{π} ∫_{−1}^{1} [(y − x) ln |y − x| − h(x, y)] \frac{u(y) dy}{\sqrt{1 − y^{2}}} = f(x), −1 < x < 1,

which is based on a collocation-quadrature method and which has the same convergence rate as this method, but only O(n log n) complexity. This integral equation turns out to be an ill-posed problem in (the best possible choice of) a pair of non-periodic Sobolev-like spaces. The present paper presents the technique, how to overcome this peculiarity in the investigation of the fast algorithm.

1. Introduction

Integral equations in which the kernel function of the integral operator contains a logarithmic part like ln|y − x| occur in various situations, in particular in applying boundary integral methods to boundary value problems for partial differential equations. For example, the single layer potential with density f(x) for the Laplace equation in two dimensions is defined by −\frac{1}{2π} ∫_{Γ} ln|y −x|f(y) dΓ_{y}, where Γ denotes the boundary of the region under consideration. In [7] kernel functions of the form (y − x)^{κ} ln|y − x|, (x, y) ∈ [−1, 1]^{2} with κ = 1 and κ = 2 are considered (cf. (2.1)), motivated by the problem of determining the shear force when a two-dimensional plate is supported along a line segment.

In [10, eq. (11.1)] kernels of the form

(1.1)  b_{0}(x, y) ln|y − x| + b_{1}(x, y), (x, y) ∈ [−1, 1]^{2}

are investigated, where the essential assumption is that the C^{∞} -function b_{0}(x, y) satisfies b_{0}(x, x) ≠ 0 for all x ∈ [−1, 1]. This is not the case under consideration here, which leads to the fact that the respective operator is not longer a Fredholm operator in pairs of Banach or Hilbert spaces which seem to be suitable chosen for the problem under consideration (cf. [10, p. 326] with Proposition 2.1 of the present paper). Here, the image of the operator is not closed (and not dense), which means that the associated operator equation is an ill-posed problem. In [6] it is shown how it is possible to handle this peculiarity in proving the applicability of a collocation-quadrature method to such an equation. It is well known that discretization can imply also a kind of regularization (see, for example, [9, Chapter 17] and the references given there). The aim of the present paper is to use the results of [6] in order to realize a known idea for the construction of a fast algorithm also in the present situation of an operator with not closed and not dense image, where we concentrate on the case κ = 1.

2010 Mathematics Subject Classification: Primary 65R20; Secondary 45B05, 45E99.

Key words and phrases: first kind integral equation, ill-posed problem, collocation method, quadrature method.

Dedicated to Giuseppe Mastroianni on the occasion of his retirement.

The outline of the paper is as follows. In Section 2 we collect the results from [7] on the mapping properties needed here for the involved operator. In Section 3 the collocation and the quadrature method and respective results from [6] are presented, which are necessary for designing and investigating the fast algorithm studied in Section 4. Finally, in Section 5 we present some numerical results confirming the theoretical ones.

2. Preliminary results

The paper [7] is devoted to solvability properties of integral equations of the form

(2.1)  − \frac{1}{π} ∫_{−1}^{1} [(y − x)^{κ} ln|y − x| − h(x, y)] \frac{u(y) dy}{\sqrt{1 − y^{2}}} = f(x), −1 < x < 1,

where κ = 0, 1, 2, while [6] deals with their numerical solution by collocation and collocation-quadrature methods. Of course, the case κ = 0 is well understood in both analytic and numerical sense (see, for example, [1, 4, 8, 11, 12]). To examine the existence and uniqueness of solutions of equations like (2.1), in [7] integral operators of the form

(A_{κ}u) (x) = −\frac{1}{π} ∫_{−1}^{1} (y − x)^{κ} ln|y − x| \frac{u(y) dy}{\sqrt{1 − y^{2}}} , −1 < x < 1,

are considered in weighted Sobolev spaces L_{σ}^{2,s}. Here, for s ⩾ 0, L_{σ}^{2,s} denotes the Hilbert space (a Sobolev-like space) of all f : (−1, 1) → ℂ satisfying

∥f∥_{σ,s} := ( ∑_{n=0}^{∞} (n + 1)^{2s} |⟨f, p_{n}^{σ}⟩_{σ}|^{2} ) ^{1/2} < ∞

equipped with the inner product ⟨f, g⟩_{σ,s} = ∑ _{n=0}^{∞} (n + 1)^{2s} ⟨f, p_{n}^{σ}⟩_{σ} \overline{⟨g, p_{n}^{σ}⟩_{σ}}, where σ(x) = (1 − x^{2})^{− \frac{1}{2}} , ⟨f, g⟩_{σ} = ∫_{−1}^{1} f(x)\overline{g(x)}σ(x) dx, and p_{n}^{σ}(x) denotes the nth normalized Chebyshev polynomial of first kind, p_{0}^{σ}(x)= \frac{1}{\sqrt{π}} , p_{n}^{σ}(coss)=\sqrt{2/π} cos(ns). Spaces of this kind were introduced in [1,3,5,11] for investigating numerical methods for different classes of singular integral equations. So, in [11] there is shown that the operator V : L_{φ}^{2,s} → L_{σ}^{2,s+1} with φ(x) = \sqrt{1 − x^{2}} and

(V u)(x) = −\frac{1}{π} ∫_{−1}^{1} ln|y − x| u(y) dy, −1 < x < 1

is an isomorphism (cf. also [12]). In [10, Chapter 11] there are considered operators with kernels of the form (1.1), provided that b_{0}(x, x) ≠ 0 for all x ∈ [−1, 1], in the pair of spaces (L_{φ}^{2} , H_{φ}^{1} ) with L_{φ}^{2} := L_{φ}^{2,0}, where the respective Sobolevlike space is constructed in another way. Namely, the space H_{φ}^{1} is defined as the set of all continuous functions f : [−1, 1] → ℂ having a distributional derivative f′ ∈ L_{φ}^{2} equipped with the norm ∥f∥_{H_{φ}^{1}} = \sqrt{∥f∥_{L_{σ}^{2}}^{2} + ∥f′∥_{L_{φ}^{2}}^{2}} . That H_{φ}^{1} equals L_{σ}^{2,1}

is proved in [1, Theorems 2.8 and 2.17]. We prefer to use the above description of the spaces L_{σ}^{2,s}, since in this way one can easily see that these spaces are isometrically isomorphic to weighted sequence spaces (see below) with the help of which we can describe the mapping properties of the involved operators very precisely (cf. Proposition 2.1 together with Remark 2.1, Lemma 2.1 and Remark 2.2).

For s < 0, L_{σ}^{2,s} is defined as the dual space (L_{σ}^{2,−s} )^{∗}. In addition to these spaces we use the sequence spaces ℓ_{s}^{2} = {ξ = (ξ_{n})_{n=0}^{∞} : ∑ _{n=0}^{∞} |ξ_{n}|^{2}(n + 1)^{2s} < ∞}, s ∈ ℝ, equipped with the inner product ⟨ξ, η⟩_{s} = ∑ _{n=0}^{∞} (n+1)^{2s} ξ_{n} \overline{η_{n}} and the norm ∥ξ∥_{s} = \sqrt{⟨ξ, ξ⟩_{s}}. Clearly, for s ⩾ 0, the operator J : L_{σ}^{2,s} → ℓ_{s}^{2}, f ↦ (⟨f, p_{n}^{σ}⟩_{σ})_{n=0}^{∞} is an isometrical isomorphism. Moreover, the dual space (L_{σ}^{2,s} )^{∗} is isometric isomorphic to ℓ_{−s}^{2} via the map J : (L_{σ}^{2,s} )^{∗} → ℓ_{−s}^{2} , v ↦ (v(p_{n}^{σ}))_{n=0}^{∞} (cf. [7, Section 5]). This allows us to consider the operators A_{κ} in all spaces L_{σ}^{2,s}, s ∈ ℝ by identifying them with the operators J^{−1} A_{κ}J, where A_{κ}ξ is defined for any complex number sequence ξ = (ξ_{n})_{n=0}^{∞} by the condition A_{κ}ξ = J A_{κ}J^{−1} ξ for all ξ ∈ ℓ^{2} (see also Proposition 2.1).

In what follows, for simplicity we will restrict to the case κ = 1 and set B := A_{1} as well as B := A_{1}. Hence, we are interested in designing a fast algorithm for the numerical solution of an integral equation of the form

(2.2)  − \frac{1}{π} ∫_{−1}^{1} [(y − x) ln|y − x| − h(x, y)] \frac{u(y) dy}{\sqrt{1 − y^{2}}} = f(x), −1 < x < 1.

Let us summarize some results from [6,7] which are important for our investigations here.

Proposition 2.1. [7, Corollaries 5.2, 5.3 and Remark 5.4] The operator B : L_{σ}^{2,s} → L_{σ}^{2,s+2} is linear and bounded. Its image is not closed in L_{σ}^{2,s+2} . Moreover,

(2.3)  Bξ = (β_{n}ξ_{n+1} − β_{n−1} ξ_{n−1} ) _{n=0}^{∞}

with ξ_{−1} := 0, β_{−1} := 0, β_{0} = \frac{ln 2−1}{\sqrt{2}} , β_{n} = \frac{1}{2n(n+1)} , n ⩾ 1. If f ∈ L_{σ}^{2,s} for some s > τ + 3 and τ ⩾ −1, then the equation

(2.4)  Bu = f + ζ_{0}p_{0}^{σ}

has a unique solution (u, ζ_{0}) ∈ L_{σ}^{2,τ} × ℂ. Furthermore, for f ≡ 0, equation (2.4) has only the trivial solution in L_{σ}^{2,− \frac{3}{2}} × ℂ.

Let us note that the assumption f ∈ L_{σ}^{2,s} with s > τ + 3 is sharp in the sense that in case s < τ + 3 the existence of a solution in L_{σ}^{2,τ} × ℂ is not guaranteed (cf. [7, Cor. 5.7]).

Remark 2.1. [7, p. 2819] For η ∈ ℓ_{s}^{2}, s > τ + 3, τ ⩾ −1, the solution of Bξ = η + (ζ_{0}, 0, …) is given by

ζ_{0} = −4β_{0} ∑_{m=1}^{∞} m η_{2m} − η_{0}, ξ_{0} = −\frac{1}{β_{0}} ∑_{m=0}^{∞} (2m + 1)η_{2m+1} ,

and, for j = 0, 1, …, by

ξ_{2j+1} = −4(2j + 1) ∑_{m=j+1}^{∞} m η_{2m} , ξ_{2(j+1)} = −4(j + 1) ∑_{m=j+1}^{∞} (2m + 1)η_{2m+1} .

For given η ∈ ℓ_{s}^{2} with s > 2, define B_{0}η = (ξ, ζ_{0}), where ξ = [ξ_{j}]_{j=0}^{∞} and ζ_{0} ∈ ℂ are given by Remark 2.1.

Lemma 2.1. [6, Lemma 3.1] Let 2 < t + 3 < s. Then the operator B_{0} : L_{σ}^{2,s} → L_{σ}^{2,t} × ℂ defined by Ĵ^{−1} B_{0}J, where Ĵ(u, ζ_{0}) = (J u, ζ_{0}), is linear and bounded and satisfies

B_{0}(Bu − ζ_{0}p_{0}^{σ}) = (u, ζ_{0}) ∀ (u, ζ_{0}) ∈ L_{σ}^{2,−1} × ℂ : Bu ∈ L_{σ}^{2,s} .

As norm of (u, ζ) in L̃_{σ}^{2,t} := L_{σ}^{2,t} × ℂ we take ∥(u, ζ)∥_{σ,t,∼} := \sqrt{∥u∥_{σ,t}^{2} + |ζ|^{2}} .

Remark 2.2. From the proof of [6, Lemma 3.1] one can see that there is a

constant c such that ∥B_{0}∥_{L_{σ}^{2,s} →_{L_{σ}^{2,t}}} ⩽ c\sqrt{∑ _{n=0}^{∞} (n + 1)^{2(t−s)+5}} =: ρ_{s−t} .

3. Collocation for the dominant equation and a quadrature method

In this section we describe the numerical method which is the basis for the construction of the fast algorithm in Section 4. Having in mind Prop. 2.1 we are interested in the numerical approximation of the solution (u, ζ_{0}) ∈ L_{σ}^{2} × ℂ of the modified (in comparison to (2.2)) equation

−\frac{1}{π} ∫_{−1}^{1} [(y − x) ln|y − x| − h(x, y)] \frac{u(y) dy}{\sqrt{1 − y^{2}}} = f(x) + \frac{ζ_{0}}{\sqrt{π}} , −1 < x < 1,

which we write shortly as

(3.1)  Bu = f + ζ_{0}p_{0}^{σ}

and which is equivalent to the equation

(3.2)  Bξ = η + (ζ_{0}, 0, 0, …), η = J f, ξ = J u.

We look for u_{n} ∈ imP_{n−1} , where P_{n} : L_{σ}^{2} → L_{σ}^{2} denotes the orthoprojection

P_{n}f = ∑_{k=0}^{n−1} ⟨f, p_{k}^{σ}⟩_{σ} p_{k}^{σ},

by solving the collocation equations

(3.3)  (Bu_{n})(x_{nj}^{σ} ) = f(x_{nj}^{σ} ) + \frac{ζ_{0}^{n}}{\sqrt{π}} , j = 1, …, n,

with the zeros x_{nj}^{σ} = cos \frac{2j−1}{2n} π of p_{n}^{σ}. Denoting by L_{n}^{σ} the Lagrange interpolating operator with respect to the nodes x_{nj}^{σ} , j = 1, …, n, the system (3.3) can be written equivalently as

(3.4)  Bu_{n} = L_{n}^{σ}f + ζ_{0}^{n}p_{0}^{σ},

taking into account Bu_{n} ∈ imP_{n} for all u_{n} ∈ imP_{n−1} (cf., (2.3)), or as

(3.5)  Bξ^{n} = (η_{0}^{n} + ζ_{0}^{n}, η_{1}^{n}, …, η_{n−1}^{n} , 0, 0, …), ξ^{n} = J u_{n} ∈ imP_{n−1} ,

with η_{k}^{n} given by

(3.6) η_{k}^{n} = \frac{π}{n} ∑_{j=1}^{n} f(x_{nj}^{σ} )p_{k}^{σ}(x_{nj}^{σ} ) = \{ \begin{matrix}\frac{\sqrt{π}}{n} ∑_{j=1}^{n} f(x_{nj}^{σ} ) & k = 0, \\ \frac{\sqrt{2π}}{n} ∑_{j=1}^{n} cos \frac{k(2j−1)π}{2n} f(x_{nj}^{σ} ) & k = 1, …, n−1. \end{matrix}

Moreover, P_{n} : ℓ^{2} → ℓ^{2}, ξ = (ξ_{n})_{n=0}^{∞} ↦ (ξ_{0}, …, ξ_{n−1} , 0, …).

By c ≠ c(n, f, …) we will mean that the constant c does not depend on n, f, …, where c can have different values at different places. We remark that (see, for example, [1, Theorem 3.4]), for s > \frac{1}{2},

(3.7)  ∥f − L_{n}^{σ}f∥ _{σ,t} ⩽ c n^{t−s} ∥f∥_{σ,s} , 0 ⩽ t ⩽ s, f ∈ L_{σ}^{2,s} ,

with c ≠ c(n, f, t). Note that (3.7) implies the uniform boundedness of L_{n}^{σ} in L_{σ}^{2,s} for s > \frac{1}{2}.

Proposition 3.1. [6, Prop. 2.1] Assume that f ∈ L_{σ}^{2,s} for some s > τ + 3, τ ⩾ 0. Then the unique solution (u_{n}, ζ_{0}^{n}) of equation (3.4) converges in L_{σ}^{2} × ℂ to the unique solution (u, ζ_{0}) ∈ L_{σ}^{2} × ℂ of equation (3.1), where

|ζ_{0}^{n} − ζ_{0}| ⩽ c n^{\frac{3}{2} −s} ∥f∥_{σ,s} and ∥u_{n} − u∥ _{σ,t} ⩽ c n^{t−τ} ∥f∥_{σ,s} , 0 ⩽ t ⩽ τ,

with c ≠ c(n, f, t).

Now, assume that s_{0} > 3 and that the continuous function h : [−1, 1]^{2} → ℂ satisfies

(A)  h(·, y) ∈ L_{σ}^{2,s_{0}} uniformly w.r.t. y ∈ [−1, 1].

Define H : L_{σ}^{2} → L_{σ}^{2} by

(Hu)(x) := \frac{1}{π} ∫_{−1}^{1} h(x, y) \frac{u(y) dy}{\sqrt{1 − y^{2}}} , −1 < x < 1.

Then condition (A) ensures that H ∈ L(L_{σ}^{2}, L_{σ}^{2,s_{0}} ) (cf. [1, Lemma 4.2]). We look for an approximate solution (u_{n}, ζ_{0}^{n}) ∈ imP_{n−1} × ℂ of the equation

(3.8)  (B + H)u = f + ζ_{0}p_{0}^{σ},

where we assume that

(B)  Equation (3.8) with f ≡ 0 has only the trivial solution in L_{σ}^{2} × ℂ.

For u ∈ L_{σ}^{2} we define H_{n}^{0}u by (H_{n}^{0}u) (x) = (L_{n}^{σ}g_{n}) (x) with

(3.9)  g_{n}(x) = \frac{1}{π} ∫_{−1}^{1} [L_{n−1}^{σ} h(x, .)] (y) \frac{u(y) dy}{\sqrt{1 − y^{2}}} , −1 < x < 1.

We look for an approximate solution (u_{n}, ζ_{0}^{n}) ∈ imP_{n−1} × ℂ of equation (3.8) by solving

(3.10)  (B + H_{n}^{0})u_{n} = L_{n}^{σ}f + ζ_{0}^{n}p_{0}^{σ}.

Taking into account the algebraic accuracy of the Gauss–Chebyshev quadrature, we can write equation (3.10) as the system

−\frac{1}{π} ∫_{−1}^{1} (y −x_{nj}^{σ} ) ln|y −x_{nj}^{σ} | \frac{u_{n}(y) dy}{\sqrt{1 − y^{2}}} + \frac{1}{n − 1} ∑_{k=1}^{n−1} h(x_{nj}^{σ} , x_{n−1,k}^{σ} )χ_{k}^{n} = f(x_{nj}^{σ} )+ \frac{ζ_{0}^{n}}{\sqrt{π}} ,

j = 1, …, n, with χ_{k}^{n} = u_{n}(x_{n−1,k}^{σ} ). Assume that

(C)  h(x, .) ∈ L_{σ}^{2,s_{0}} uniformly w.r.t. x ∈ [−1, 1].

Lemma 3.1. [6, Lemma 3.2] Assume that 3 ⩽ t + 3 < s < s_{0} and that conditions (A) and (B) are fulfilled. There exists a constant γ_{ts} > 0 such that, for all sufficiently large n,

∥(B + L_{n}^{σ}H)u − ζ_{0}p_{0}^{σ}∥ _{σ,s} ⩾ γ_{ts} ∥(u, ζ_{0})∥_{σ,t,∼} ∀ (u, ζ_{0}) ∈ L_{σ}^{2} × ℂ : Bu ∈ L_{σ}^{2,s} .

Lemma 3.2. Let (C) be fulfilled with s_{0} > \frac{1}{2}. Then, for t ⩾ 0,

∥(L_{n}^{σ}H − H_{n}^{0})u∥ _{σ,t} ⩽ c n^{t−s_{0}} ∥u∥_{σ,0} ∀ u ∈ L_{σ}^{2},

where c ≠ c(n, t, u).

Proof. Since (L_{n}^{σ}H − H_{n}^{0})u ∈ imP_{n} we have
∥(L_{n}^{σ}H − H_{n}^{0}) u∥_{σ,t}^{2}

⩽ n^{2t} ∥(L_{n}^{σ}H − H_{n}^{0})u∥ _{σ,0}^{2}

= n^{2t} \frac{π}{n} ∑_{j=1}^{n} | \frac{1}{π} ∫_{−1}^{1} u(y) [h(x_{nj}^{σ} , .) − L_{n−1}^{σ} h(x_{nj}^{σ} , .)] (y)σ(y) dy|^{2}

⩽ n^{2t}∥u∥_{σ,0}^{2} \frac{1}{πn} ∑_{j=1}^{n} ∥h(x_{nj}^{σ} , .) − L_{n−1}^{σ} h(x_{nj}^{σ} , .)∥ _{σ,0}^{2}

⩽^{(3.7)} c n^{2(t−s_{0} )} ∥u∥_{σ,0}^{2} \frac{1}{n} ∑_{j=1}^{n} ∥h(x_{nj}^{σ} , .)∥_{σ,s_{0}}^{2} ⩽ c n^{2(t−s_{0} )} ∥u∥_{σ,0}^{2} ,

where we used the algebraic accuracy of the Gaussian rule and (3.7) together with assumption (C). □

Corollary 3.1. Additionally to the assumptions of Lemma (3.1), let condition (C) be satisfied. Then, for all sufficiently large n,

∥(B + H_{n}^{0})u − ζ_{0}p_{0}^{σ} ∥ _{σ,s} ⩾ \frac{γ_{ts}}{2} ∥(u, ζ_{0})∥_{σ,t,∼} ∀ (u, ζ_{0}) ∈ L_{σ}^{2} × ℂ : Bu ∈ L_{σ}^{2,s} .

Proof. From Lemma 3.1 and Lemma 3.2 we get

γ_{ts} ∥(u, ζ_{0})∥_{σ,t,∼} ⩽ ∥(B + H_{n}^{0})u − ζ_{0}p_{0}^{σ} ∥ _{σ,s} + ∥(L_{n}^{σ}H − H_{n}^{0}) u∥_{σ,s}
⩽ ∥(B + H_{n}^{0})u − ζ_{0}p_{0}^{σ} ∥ _{σ,s} + c n^{s−s_{0}} ∥(u, ζ_{0})∥_{σ,s,∼} ,

which proves the corollary. □

Proposition 3.2. [6, Prop. 4.1] Let 3 ⩽ 3 + τ < s < s_{0}, 0 ⩽ t ⩽ τ,

f ∈ L_{σ}^{2,s}, and assume that conditions (A), (B), and (C) are fulfilled. Then, for all sufficiently large n, equation (3.10) has a unique solution (u_{n}, ζ_{0}^{n}). Moreover, there is a constant c ≠ c(n, f, t) such that

|ζ_{0}^{n} − ζ_{0}^{∗}| ⩽ c n^{−τ} ∥f∥_{σ,s} and ∥u_{n} − u^{∗}∥ _{σ,t} ⩽ c n^{t−τ} ∥f∥_{σ,s}

where (u^{∗}, ζ_{0}^{∗}) ∈ L_{σ}^{2,τ} × ℂ is the unique solution of (3.8).

The proof of the unique solvability of (3.8) in L_{σ}^{2,τ} × ℂ is part of the proof of [6, Prop. 4.1].

Remark 3.1. Since χ_{i}^{n} = ∑ _{k=0}^{n−2} p_{k}^{σ}(x_{n−1,i}^{σ} )ξ_{k}^{n} holds for ξ^{n} = [ξ_{k}^{n}]_{k=0}^{n−2} (cf. (3.5)), a matrix version of equation (3.10) is given by

(B_{n} + D_{n}C_{n}^{2}H_{n}C_{n−1}^{3} D_{n−1} )ξ^{n} = η^{n} + [ζ_{0}^{n} 0 ⋯ 0]^{T} ,

where

B_{n} = [ \begin{matrix}0 & β_{0} & 0 & ⋯ & 0 \\ −β_{0} & 0 & β_{1} & ⋯ & 0 \\ & ... & ... & ... \\ 0 & ⋯ & −β_{n−4} & 0 & β_{n−3} \\ 0 & ⋯ & 0 & −β_{n−3} & 0 \\ 0 & ⋯ & 0 & 0 & −β_{n−2} \end{matrix} ] , (3.11)

H_{n} = [h(x_{nj}^{σ} , x_{n−1,k}^{σ} )]_{j=1, k=1}^{n, n−1} , C_{n}^{2} = [cos \frac{j(2k−1)π}{2n} ]_{j=0,k=1}^{n−1, n} , (3.12)

C_{n−1}^{3} = (C_{n−1}^{2} )^{T} , D_{n} = \frac{1}{n} diag[1 \sqrt{2} ⋯ \sqrt{2}], (3.13)

and the vector η^{n} = [η_{k}^{n}]_{k=0}^{n−1} is defined by (3.6).

4. A fast algorithm

Now, our aim is to design a fast algorithm according to the pattern given in [1], which was also used, for example, in [2]. So, on one hand, we want to reduce the complexity of the method to O(n logn) and on the other hand, we want to keep the convergence rate stated in Propositon 3.2, at least in a smaller range for the parameter t (see Prop. 4.3). For this, we use the structure properties of that part of the operator, which is given by the (y − x) ln|y − x|-kernel (cf. (3.11)), and the smoothing properties of the remaining part, which are due to the assumptions (A) and (C).

In what follows we assume that the conditions (A), (B), and (C) are fulfilled. Let us remark that the main difference of the present situation to the situations in [1, 2] is that the main part B : L_{σ}^{2,s} → L^{2,s+2} of the operator B + H of the original equation does not have a closed image and, consequently, is not bounded invertible.

We search for an approximate solution (u_{n}, ζ_{0}^{n}) ∈ imP_{n−1} ×ℂ of equation (3.8) and write u_{n} in the form

(4.1)  u_{n} = ∑_{k=0}^{m−2} ξ_{k}^{n}p_{k}^{σ} + ∑_{k=m−1}^{n−2} ξ_{k}^{n}p_{k}^{σ},

where 2 < m < n.

For k = m − 1, …, n − 2, we set ξ_{k}^{n} := ⟨v_{n}, p_{k}^{σ}⟩_{σ}, where (v_{n}, ζ̃_{0}^{n}) ∈ imP_{n−1} × ℂ is the unique solution of (cf. (3.4) and (3.5))

(4.2)  Bv_{n} = L_{n}^{σ}f + ζ̃_{0}^{n}p_{0}^{σ}.

Equation (4.2) is equivalent to

B_{n}ξ^{n} = η^{n} + [ζ̃_{0}^{n} 0 … 0]^{T} with ξ^{n} = [ξ_{k}^{n}]_{k=0}^{n−2} ,

where the matrix B_{n} is defined in (3.11) and where the entries of

η^{n} = [η_{k}^{n}]_{k=0}^{n−1} = \frac{\sqrt{π}}{n} D_{n}C_{n}^{2}[f(x_{nj}^{σ} )]_{j=1}^{n}

(cf. (3.12), (3.13)) can be computed with O(n logn) complexity. Due to the structure of the matrix B_{n} we have

ξ_{n−2}^{n} = −\frac{η_{n−1}^{n}}{β_{n−2}} , ξ_{n−3}^{n} = −\frac{η_{n−2}^{n}}{β_{n−3}} , ξ_{k}^{n} = \frac{β_{k+1} ξ_{k+2}^{n} − η_{k+1}^{n}}{β_{k}} , k = n − 4, …, m − 1.

which can be realized with O(n) complexity. Moreover, due to Lemma 2.1,

(4.3)  (v_{n}, ζ̃_{0}^{n}) = B_{0}L_{n}^{σ}f.

Define Q_{n} := I − P_{n} : L_{σ}^{2} → L_{σ}^{2}, i.e., Q_{n}u = ∑ _{k=n}^{∞} ⟨u, p_{k}^{σ}⟩_{σ} p_{k}^{σ}. Note that

(4.4)  ∥Q_{n}u∥_{σ,t} ⩽ (n + 1)^{t−s} ∥u∥_{σ,s} , 0 ⩽ t ⩽ s, u ∈ L_{σ}^{2,s} .
Proposition 4.1. Let 3 ⩽ τ + 3 < s < s_{0} and f ∈ L_{σ}^{2,s}. Furthermore, let

(u^{∗}, ζ_{0}^{∗}) ∈ L_{σ}^{2,τ} × ℂ be the unique solution of (3.8) and (v_{n}, ζ̃_{0}^{n}) ∈ imP_{n−1} × ℂ be the unique solution of (4.2). Then, for every ε ∈ (0, s − τ − 3) there is a constant c ≠ c(m, n, f, t) such that, for 0 ⩽ t ⩽ τ + s_{0} − s + ε,

∥Q_{m−1} u^{∗} − Q_{m−1} v_{n}∥ _{σ,t} ⩽ c(m^{t−τ−s_{0} +s−ε} ∥u^{∗}∥_{σ,0} + n^{t−τ−ε} ∥f∥_{σ,s} ).
Proof. Set θ := s_{0} − s and define R : L̃_{σ}^{2,t} → L_{σ}^{2,t}, (u, ζ) ↦ u. Since Bu^{∗} =

f + ζ_{0}^{∗}p_{0}^{σ} − Hu^{∗} ∈ L_{σ}^{2,s} we have, in view of Lemma 2.1,

B_{0}f = B_{0}(Bu^{∗} + Hu^{∗} − ζ_{0}^{∗}p_{0}^{σ}) = (u^{∗}, ζ_{0}^{∗}) + B_{0}Hu^{∗}.

Consequently, for 0 ⩽ t ⩽ τ,

∥Q_{m−1} u^{∗} − Q_{m−1} v_{n}∥ _{σ,t}
=  ∥Q_{m−1} R(u^{∗}, ζ_{0}^{∗}) − Q_{m−1} R(v_{n}, ζ̃_{0}^{n})∥_{σ,t}
=^{(4.3)}  ∥Q_{m−1} R(B_{0}f − B_{0}Hu^{∗}) − Q_{m−1} RB_{0}L_{n}^{σ}f∥_{σ,t}
≤  ∥Q_{m−1} RB_{0}Hu^{∗}∥_{σ,t} + ∥Q_{m−1} RB_{0}(f − L_{n}^{σ}f)∥ _{σ,t}
≤^{(4.4)}  m^{t−(τ+θ+ε)} ∥B_{0}Hu^{∗}∥_{σ,τ+θ+ε,∼} + ∥B_{0}(f − L_{n}^{σ}f)∥ _{σ,t,∼}
≤^{Lemma 2.1} m^{t−τ−θ−ε} ∥B_{0}∥_{L_{σ}^{2,s_{0}} →_{L_{σ}^{2,τ+θ+ε}}} ∥Hu^{∗}∥_{σ,s_{0}}

+ ∥B_{0}∥_{L_{σ}^{2,t+3+ε_{1}} →_{L_{σ}^{2,t}}} ∥f − L_{n}^{σ}f∥ _{σ,t+3+ε_{1}}

⩽^{(3.7)}  c(m^{t−τ−θ−ε} ∥u^{∗}∥_{σ,0} + n^{t+3+ε_{1} −s} ∥f∥_{σ,s} ),

where ε_{1} := s − τ − 3 − ε > 0, so that t + 3 + ε_{1} − s = t − τ − ε. Moreover, we took into account that, due to Remark 2.2, ∥B_{0}∥_{L_{σ}^{2,t+3+ε_{1}} →_{L_{σ}^{2,t}}} ⩽ ρ_{3+ε_{1}} independently of t. □

The numbers ξ_{k}^{n}, k = 0, …, m − 2 in (4.1) we take from the solution w_{m} ∈ imP_{m−1} of (cf. (3.10))

(4.5)  (B + H_{n}^{0}) w_{m} = L_{n}^{σ}f + ζ̂_{0}^{m}p_{0}^{σ} − L_{m}^{σ}BQ_{m−1} v_{n},

where v_{n} ∈ imP_{n−1} is the solution of (4.2), i.e., ξ_{k}^{n} := ⟨w_{m}, p_{k}^{σ}⟩_{σ}, k = 0, …, m − 2. Moreover, we set ζ_{0}^{n} := ζ̂_{0}^{m}.

In view of Remark 3.1 the solving (4.5) is equivalent to solving the system

(B_{m} + D_{m}C_{m}^{2}H_{m}C_{m−1}^{3} D_{m−1} )ξ̂^{m} = η^{m} + [ζ̂_{0}^{m} 0 ⋯0]^{T},

which can be done with O(m^{3}) complexity if the numbers

η_{k}^{m} = ⟨L_{m}^{σ}f − L_{m}^{σ}BQ_{m−1} v_{n}, p_{k}^{σ}⟩ _{σ} , k = 0, …, m − 1,

are computed. In order to do this analogously to (3.6), we need the function values (BQ_{m−1} v_{n})(x_{mj}^{σ} ), j = 1, …, m. To get these values effectively we can choose m and n in such a way that \frac{m}{n} is an odd integer, since in this case we have {x_{mj}^{σ} : j = 1, …, m} ⊂ {x_{nj}^{σ} : j = 1, …, n}. Let us write

BQ_{m−1} v_{n} = ∑_{k=0}^{n−1} χ_{k}^{n}p_{k}^{σ}.

In view of Proposition 2.1 we get, setting ξ_{n−1}^{n} = ξ_{n}^{n} = 0,

BQ_{m−1} v_{n} = ∑_{k=m−1}^{n−2} ξ_{k}^{n} (β_{k−1} p_{k−1}^{σ} − β_{k}p_{k+1}^{σ} )

= β_{m−2} ξ_{m−1}^{n} p_{m−1}^{σ} + β_{m−1} ξ_{m}^{n}p_{m}^{σ} + ∑_{k=m}^{n−1} (β_{k}ξ_{k+1}^{n} − β_{k−1} ξ_{k−1}^{n} ) p_{k}^{σ},

such that, due to (4.2),

χ_{k}^{n} := \{ \begin{matrix}0 & 0 ⩽ k ⩽ m − 3, \\ β_{k}ξ_{k+1}^{n} & k = m − 2, m − 1, \\ η_{k}^{n} & m ⩽ k ⩽ n − 1. \end{matrix}

Finally,

(BQ_{m−1} v_{n})(x_{nj}^{σ} ) = ∑_{k=0}^{n−1} χ_{k}p_{k}^{σ}(x_{nj}^{σ} ), j = 1, …, n,

which can be realized with O(n logn) complexity using the discrete cosine transform C_{n−1}^{2} . Consequently, (4.5) can be solved with O(m^{3}) + O(n logn) = O(n logn) complexity if \frac{m^{3}}{n} ⩽ γ for some constant γ.

Proposition 4.2. Let 3 ⩽ τ + 3 < s < s_{0} and f ∈ L_{σ}^{2,s}. Assume that

there is an ε ∈ (0, s − τ − 3) such that τ + s_{0} − s − 1 + ε > 0. Then, for every δ ∈ (0, τ + s_{0} − s − 1 + ε], there exists a constant c ≠ c(m, n, f) such that, for all sufficiently large m and n with 0 < m < n,

\sqrt{∥w_{m} − P_{m−1} u^{∗}∥ _{σ,0}^{2} + |ζ̂_{0}^{m} − ζ_{0}^{∗} |^{2}} ⩽ c(m^{1+δ−τ−s_{0} +s−ε} ∥u^{∗}∥_{σ,0} + n^{1+δ−τ−ε} ∥f∥_{σ,s} )

where (u^{∗}, ζ_{0}^{∗}) and (w_{m}, ζ̂_{0}^{m}) are the solutions of (3.8) and (4.5), respectively.

Proof. We compute

(B + H_{m}^{0})(w_{m} − P_{m−1} u^{∗})
= L_{m}^{σ}f + ζ̂_{0}^{m}p_{0}^{σ} − L_{m}^{σ}BQ_{m−1} v_{n} − (B + H_{m}^{0})P_{m−1} u^{∗}
= L_{m}^{σ}(B + H)u^{∗} − ζ_{0}^{∗}p_{0}^{σ} + ζ̂_{0}^{m}p_{0}^{σ} − L_{m}^{σ}BQ_{m−1} v_{n} − (B + H_{m}^{0})P_{m−1} u^{∗}
= (L_{m}^{σ} − I)BP_{m−1} u^{∗} + (ζ̂_{0}^{m} − ζ_{0}^{∗})p_{0}^{σ}
+ L_{m}^{σ}BQ_{m−1} (u^{∗} − v_{n}) + (L_{m}^{σ}H − H_{m}^{0})u^{∗} + H_{m}^{0}Q_{m−1} u^{∗}
= L_{m}^{σ}BQ_{m−1} (u^{∗} − v_{n}) + (ζ̂_{0}^{m} − ζ_{0}^{∗})p_{0}^{σ} + (L_{m}^{σ}H − H_{m}^{0})u^{∗},

where we took into account that (L_{m}^{σ} − I) BP_{m−1} u^{∗} = 0 because of BP_{m−1} u^{∗} ∈ imP_{m} and that H_{m}^{0}Q_{m−1} u^{∗} = L_{m}^{σ}g_{m} with (cf. (3.9))

g_{m}(x) = \frac{1}{π} ∫_{−1}^{1} [L_{m−1} h(x, .)](y) (Q_{m−1} u^{∗})(y) σ(y) dy

= \frac{1}{π} ∫_{−1}^{1} ∑_{k=1}^{m−1} h(x, x_{m−1,k}^{σ} ) ℓ_{m−1,k}^{σ} (y)(Q_{m−1} u^{∗})(y) σ(y) dy

= \frac{1}{π} ∑_{k=1}^{m−1} h(x, x_{m−1,k}^{σ} ) ⟨Q_{m−1} u^{∗}, ℓ_{m−1,k}^{σ} ⟩_{σ} = 0,

where ℓ_{m−1,k}^{σ} denote the respective fundamental Lagrange interpolation polynomials of degree m−2. With the help of Cor. 3.1, Prop. 2.1, Prop. 4.1, and Lemma 3.2 we get

∥(w_{m} − P_{m−1} u^{∗}, ζ̂_{0}^{m} − ζ_{0}^{∗})∥_{σ,0,∼}

⩽ \frac{2}{γ_{0,3+δ}} ∥(B + H_{m}^{0})(w_{m} − P_{m−1} u^{∗}) − (ζ̂_{0}^{m} − ζ_{0}^{∗})p_{0}^{σ} ∥_{σ,3+δ}

⩽ c( ∥L_{m}^{σ}BQ_{m−1} (u^{∗} − v_{n})∥ _{σ,3+δ} + ∥(L_{m}^{σ}H − H_{m}^{0})u^{∗} ∥ _{σ,3+δ} )

⩽ c( ∥Q_{m−1} (u^{∗} − v_{n})∥ _{σ,1+δ} + ∥(L_{m}^{σ}H − H_{m}^{0})u^{∗} ∥ _{σ,3+δ} )

⩽ c( m^{1+δ−τ−s_{0} +s−ε} ∥u^{∗}∥_{σ,0} + n^{1+δ−τ−ε} ∥f∥_{σ,s} + m^{3+δ−s_{0}} ∥u^{∗}∥_{σ,0} )

⩽ c( m^{1+δ−τ−s_{0} +s−ε} ∥u^{∗}∥_{σ,0} + n^{1+δ−τ−ε} ∥f∥_{σ,s} ) ,

and the assertion is proved. □

The following proposition shows that, for the O(n logn)-algorithm, we can attain the same convergence rate (in a restricted interval for t) as in Prop. 3.2 under the assumptions that s_{0} > s > 4.

Proposition 4.3. Let u_{n} in (4.1) together with ζ_{0}^{n} = ζ̂_{0}^{m} be determined by the above described algorithm and let (u^{∗}, ζ_{0}) be the solution of (3.8), where we assume that there is a constant γ > 0 such that \frac{n}{m^{3}} ⩽ γ and that 4 ⩽ τ + 4 < s < s_{0}, f ∈ L_{σ}^{2,s}. Then, there is a constant c ≠ c(m, n, t) such that

∥u_{n} − u^{∗}∥ _{σ,t} ⩽ c n^{t−τ} , max{0, τ − (s_{0} − s)/2} ⩽ t ⩽ τ.

Moreover, if additionally τ ⩽ (s_{0} − s)/2, then |ζ_{0}^{n} − ζ_{0}^{∗}| ⩽ c n^{−τ} .

Proof. Choose ε ∈ (1, s−τ −3) and δ = ε−1. Then 0 < δ < τ +s_{0}−s−1+ε,

and we can apply Prop. 4.1 and Prop. 4.2 to estimate

∥u^{∗} − u_{n}∥ _{σ,t} ⩽ ∥P_{m−1} u^{∗} − P_{m−1} u_{n}∥ _{σ,t} + ∥Q_{m−1} u^{∗} − Q_{m−1} u_{n}∥ _{σ,t}
= ∥P_{m−1} u^{∗} − w_{m}∥ _{σ,t} + ∥Q_{m−1} u^{∗} − Q_{m−1} v_{n}∥ _{σ,t}

⩽ c( m^{t} ∥P_{m−1} u^{∗} − w_{m}∥ _{σ,0} + m^{t−τ−s_{0} +s−ε} ∥u^{∗}∥_{σ,0} + n^{t−τ−ε} ∥f∥_{σ,s} )

⩽ c( m^{t−τ−s_{0} +s} ∥u^{∗}∥_{σ,0} + m^{t}n^{−τ} ∥f∥_{σ,s} + n^{t−τ−ε} ∥f∥_{σ,s} ) ⩽ cn^{t−τ}

if τ − (s_{0} − s)/2 ⩽ t. Furthermore, due to Prop. 4.2,

|ζ_{0}^{n} − ζ_{0}^{∗}| = |ζ̂_{0}^{m} − ζ_{0}^{∗} | ⩽ c(m^{−τ−s_{0} +s} ∥u^{∗}∥_{σ,0} + n^{−τ} ∥f∥_{σ,s} ) ⩽ cn^{−τ}

if t ⩽ (s_{0} − s)/2. □

5. Numerical examples

As an example we take the equation

(5.1)  − \frac{1}{π} ∫_{−1}^{1} [(y − x) ln|y − x| − cos(xy + 2y)] \frac{u(y) dy}{\sqrt{1 − y^{2}}} = f(x) + \frac{ζ_{0}}{\sqrt{π}}

which was already considered in [6, Section 5]. For the right-hand side function f(x) we choose f(x) = x^{3}|x|, f(x) = x^{4}|x|, and

(5.2)  f(x) = \frac{1}{π}[ \frac{2 sin(x + 2)}{x + 2} − \frac{1}{2}(1 − x)^{2}( ln(1 − x) − \frac{1}{2})

+ \frac{1}{2}(1 + x)^{2}( ln(1 + x) − \frac{1}{2}) ].

In case of (5.2), u^{∗}(x) = \sqrt{1 − x^{2}} with ζ_{0}^{∗} = 0 is the exact solution. In the following tables one can see the numerical results (we restrict to 9 decimal digits) of the described fast algorithm (FA) for various choices of n and m in comparison with the results of the quadrature method (QM) presented in [6]. The computations are performed with double precision.

Equation (5.1), Fourier coefficients, f(x) = x^{3}|x|

n m  ξ_{0}^{n}  ξ_{1}^{n}  ξ_{2}^{n}  ξ_{3}^{n}  ξ_{4}^{n}
FA 256 8 -8.47904757 -1.12502515 -7.02640549 0.14626959 -1.38702022
FA 256 16 -8.47874416 -1.12489088 -7.02624013 0.14632797 -1.38700379
QM 256  -8.47874416 -1.12489088 -7.02624013 0.14632797 -1.38700379
FA 512 8 -8.47904761 -1.12502516 -7.02640553 0.14626959 -1.38702028
FA 512 16 -8.47874420 -1.12489089 -7.02624017 0.14632797 -1.38700386
QM 512  -8.47874420 -1.12489089 -7.02624017 0.14632797 -1.38700386
FA 1024 8 -8.47904761 -1.12502516 -7.02640554 0.14626959 -1.38702029
FA 1024 16 -8.47874420 -1.12489089 -7.02624017 0.14632797 -1.38700386

Equation (5.1), ζ_{0}^{n} and function values, f(x) = x^{3}|x|

n m  ζ_{0}^{n}  u_{n}(−0.9)  u_{n}(0)  u_{n}(0.5)  u_{n}(0.95)
FA 256 8 1.36088509 -7.34653331 -0.64725315 -1.79908541 -10.4375062
FA 256 16 1.36088532 -7.34636880 -0.64720210 -1.79895947 -10.4371124
QM 256  1.36088532 -7.34636881 -0.64720206 -1.79895946 -10.4371124
FA 512 8 1.36088510 -7.34653709 -0.64731563 -1.79908137 -10.4375044
FA 512 16 1.36088533 -7.34637260 -0.64726455 -1.79895542 -10.4371106
QM 512  1.36088533 -7.34637260 -0.64726455 -1.79895542 -10.4371106
FA 1024 8 1.36088510 -7.34653652 -0.64733128 -1.79908188 -10.4375040
FA 1024 16 1.36088533 -7.34637202 -0.64728020 -1.79895593 -10.4371102

In cases n = 256, 512, for the fast algorithm we get the same results (with respect to 9 decimal digits) as in the quadrature method of [6] already for m = 16, which remains true if we increase n. In the following two tables we consider a more smooth right-hand side which results in the fact, that the fast algorithm gives the same results as the quadrature method already for m = 8.

Equation (5.1), Fourier coefficients, f(x) = x^{4}|x|

n m  ξ_{0}^{n}  ξ_{1}^{n}  ξ_{2}^{n}  ξ_{3}^{n}  ξ_{4}^{n}
FA 256 8 0.00000000 -4.25538432 0.00000000 -5.47120842 0.00000000
QM 256  0.00000000 -4.25538432 0.00000000 -5.47120842 0.00000000
FA 512 8 0.00000000 -4.25538432 0.00000000 -5.47120842 0.00000000
QM 512  0.00000000 -4.25538432 0.00000000 -5.47120842 0.00000000
FA 1024 8 0.00000000 -4.25538432 0.00000000 -5.47120842 0.00000000

Equation (5.1), ζ_{0}^{n} and function values, f(x) = x^{4}|x|

n m  ζ_{0}^{n}  u_{n}(−0.9)  u_{n}(0)  u_{n}(0.5)  u_{n}(0.95)
FA 256 8 0.32152132 3.57217748 0.00000000 2.34570884 -5.78339205
QM 256  0.32152132 3.57217748 0.00000000 2.34570884 -5.78339205
FA 512 8 0.32152132 3.57217748 0.00000000 2.34570883 -5.78339205
QM 512  0.32152132 3.57217748 0.00000000 2.34570883 -5.78339205
FA 1024 8 0.32152132 3.57217748 0.00000000 2.34570883 -5.78339205

In the last example (see the following two tables), the right-hand side is less smooth. But, beginning with n = 1024 the results are exact (w.r.t. 9 or 8 decimal digits) already for m = 16.

Equation (5.1), Fourier coefficients, f(x) from (5.2)

n m  ξ_{0}^{n}  ξ_{1}^{n}  ξ_{2}^{n}  ξ_{3}^{n}  ξ_{4}^{n}
FA 256 8 1.12848920 0.00004609 -0.53186634 0.00001504 -0.10638087
FA 256 16 1.12837951 0.00000008 -0.53192266 -0.00000001 -0.10638403
QM 256  1.12837951 0.00000008 -0.53192266 -0.00000001 -0.10638403
FA 512 8 1.12848890 0.00004602 -0.53186667 0.00001505 -0.10638138
FA 512 16 1.12837921 0.00000001 -0.53192299 0.00000000 -0.10638454
QM 512  1.12837921 0.00000001 -0.53192299 0.00000000 -0.10638454
FA 1024 8 1.12848886 0.00004601 -0.53186671 0.00001505 -0.10638144
FA 1024 16 1.12837917 0.00000000 -0.53192303 0.00000000 -0.10638460
u^{∗}  1.12837917 0.00000000 -0.53192304 0.00000000 -0.10638461

Equation (5.1), ζ_{0}^{n} and function values, f(x) from (5.2)

n m  ζ_{0}^{n}  u_{n}(−0.9)  u_{n}(0)  u_{n}(0.5)  u_{n}(0.95)
FA  256 8 0.00000237 0.43594383 1.00001952 0.86607051 0.31238988
FA  256 16 -0.00000007 0.43588980 1.00000006 0.86602553 0.31224956
QM 256 -0.00000007 0.43588980 1.00000006 0.86602553 0.31224956
FA  512 8 0.00000243 0.43594390 1.00001947 0.86607040 0.31239018
FA  512 16 -0.00000001 0.43588988 1.00000001 0.86602542 0.31224986
QM 512 -0.00000001 0.43588988 1.00000001 0.86602542 0.31224986
FA 1024 8 0.00000244 0.43594392 1.00001947 0.86607039 0.31239022
FA 1024 16 0.00000000 0.43588989 1.00000001 0.86602541 0.31224990
u^{∗}(x)  0.43588989 1.00000000 0.86602540 0.31224990

In the following two tables, let us show some results with 15 decimal digits for the third example. One can see that, due to rounding errors, the precision of the results decreases when going from n = 16384 to n = 32768.

Equation (5.1), Fourier coefficients, f(x) from (5.2)

n m  ξ_{0}^{n}  ξ_{4}^{n}
FA 8192 8 1.12848885640610 -0.10638144804001
FA 8192 16 1.12837916711157 -0.10638460808187
FA 8192 32 1.12837916711008 -0.10638460808242
FA 16384 8 1.12848885640011 -0.10638144805023
FA 16384 16 1.12837916710552 -0.10638460809210
FA 16384 32 1.12837916710401 -0.10638460809269
FA 32768 8 1.12848885641034 -0.10638144803278
FA 32768 16 1.12837916711585 -0.10638460807464
FA 32768 32 1.12837916711433 -0.10638460807525
u^{∗}  1.12837916709551 -0.10638460810705

Equation (5.1), ζ_{0}^{n} and function values, f(x) from (5.2)

n m  u_{n}(−0.9)  u_{n}(0.5)
FA  8192 8 0.43594391928916 0.86607039879272
FA  8192 16 0.43588989107899 0.86602541682333
FA  8192 32 0.43588989107818 0.86602541682336
FA 16384 8 0.43594401611548 0.86607034614225
FA 16384 16 0.43588998790529 0.86602536417281
FA 16384 32 0.43588998790464 0.86602536417291
FA 32768 8 0.43594421947381 0.86607061508862
FA 32768 16 0.43589019126367 0.86602563311924
FA 32768 32 0.43589019126272 0.86602563311909
u^{∗}(x)  0.43588989435407 0.86602540378444

The last table presents the comparison of the CPU-time used for establishing and solving the system by the quadrature method and by the fast algorithm (in case of m = 16).

CPU-time in seconds for QM and FA (m = 16)

n  64 128 256 512 1024 2048 4096 8192 16384 32768
QM 0.02 0.16 1.25 9.52
FA 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.02 0.03 0.06

References

1. D. Berthold, W. Hoppe, B. Silbermann, A fast algorithm for solving the generalized airfoil equation, J. Comput. Appl. Math. 43 (1992), 185–219.

2. M. R. Capobianco, G. Criscuolo, P. Junghanns, A fast algorithm for Prandtl’s integrodifferential equation, J. Comput. Appl. Math. 77 (1997), 103–128.

3. V. J. Ervin, E. P. Stephan, Collocation with Chebyshev polynomials for a hypersingular integral equation on an interval, J. Comput. Appl. Math. 43 (1992), 221–229.

4. R. Estrada, R. P. Kanwal, Integral equations with logarithmic kernels, IMA J. Appl. Math. 43 (1989), 133–155.

5. T. Hartmann, E. P. Stephan, Rates of convergence for collocation with Jacobi polynomials for the airfoil equation, J. Comput. Appl. Math. 51 (1994), 179–191.

6. P. Junghanns, G. Monegato, On the numerical solution of integral equations with logarithmic kernels, Rend. Circ. Mat. Palermo (2) 82 (2010), 1–20.

7. P. Junghanns, G. Monegato, A. Strozzi, On the integral equation formulations of some 2D contact problems, J. Comput. Appl. Math. 234 (2010), 2808–2825.

8. R. P. Kanwal, Linear Integral Equations, Theory and Technique, Birkhäuser-Verlag, Boston, Basel, Berlin, 1997.

9. R. Kress, Linear Integral Equations, Springer-Verlag, New York, 1999.

10. J. Saranen, G. Vainikko, Periodic Integral and Pseudodifferential Equations with Numerical Approximation, Springer-Verlag, Berlin, Heidelberg, New York, 2002.

11. I. H. Sloan, E. P. Stephan, Collocation with Chebyshev polynomials for Symm’s integral euqtion on an interval, J. Austral. Math. Soc., Ser. B 34 (1992), 199–211.

12. Y. Yan, I. H. Sloan, On integral equations of the first kind with logarithmic kernels, J. Integral Equations Appl. 1 (1988), 517–548.

Fakultät für Mathematik, Technische Universität Chemnitz, Germany

 katharina.flemming@mathematik.tu-chemnitz.de
 peter.junghanns@mathematik.tu-chemnitz.de