*Article* **Approximate Methods for Solving Linear and Nonlinear Hypersingular Integral Equations**

#### **Ilya Boykov 1,\*, Vladimir Roudnev <sup>2</sup> and Alla Boykova <sup>1</sup>**


Received: 2 June 2020; Accepted: 25 June 2020 ; Published: 1 July 2020

**Abstract:** We propose an iterative projection method for solving linear and nonlinear hypersingular integral equations with non-Riemann integrable functions on the right-hand sides. We investigate hypersingular integral equations with second order singularities. Today,hypersingular integral equations of this type are widely used in physics and technology. The convergence of the proposed method is based on the Lyapunov stability theory of solutions of ordinary differential equation systems. The advantage of the method for linear equations is in simplicity of unique solvability verification for the approximate equations system in terms of the operator logarithmic norm. This makes it possible to estimate the norm of the inverse matrix for an approximating system. The advantage of the method for nonlinear equations is that neither the existence or reversibility of the nonlinear operator derivative is required. Examples are given illustrating the effectiveness of the proposed method.

**Keywords:** hypersingular integral equations; iterative projection method; Lyapunov stability theory

#### **1. Introduction**

The importance of developing analytical and numerical methods for solving hypersingular integral equations is determined by a variety of fields of mathematics and by applications that use hypersingular integral equations.

Hadamard introduced the concept of a finite part of an integral, or the hypersingular integral in modern terminology, when studying hyperbolic equations. The Riemann boundary problem leads to hypersingular integral equations in exceptional cases. The boundary integral equations method reduces the dimensions of partial differential equations; that leads to hypersingular integral equations.

Hypersingular integral equations, singular integral equations and Riemann boundary problem are widely used in aerodynamics, electrodynamics, quantum physics, antennae theory and many other fields of physics and engineering [1–5].

Analytical methods for solving singular and hypersingular integral equations are known only for certain particular types of equations [6–8]. Thus, the importance of constructing numerical solutions is clear.

Developing numerical methods for solving singular integral equations began in the middle of the last century. By now, exhaustive results have been obtained for many types of equations. A detailed account of numerical methods for solving singular integral equations as well as numerous bibliography references can be found in [9–14].

Numerical methods for solving hypersingular integral equations have been developed to a much lesser extent. Mostly numerical methods to solve hypersingular integral equations of the first kind have been developed. Numerical methods for solving hypersingular integral equations of the second kind have been much less developed. Apparently, hypersingular integral equations of the first kind are more common. Naturally, the equations of the first kind are widely used in aerodynamics (one-dimensional [15] and multi-dimensional [5,16] Prandtl equation), electrodynamics, antennae theory, etc.

The following methods are used in solving hypersingular integral equations of the first kind.

Collocations, mechanical quadratures and Galerkin methods were employed to solve equations with *p* = 2 singularity [6,17–19].

Approximate methods for solving hypersingular integral equations having singularities of order *p* = 2, 3, . . . , and defined on closed smooth integration contours are constructed in [20].

In [21,22] spline-collocation methods for solving hypersingular and polyhypersingular integral equations of the second kind with odd and even singularities have been developed and justified. The spline-collocation methods for solving nonlinear hypersingular and polyhypersingular integral equations have been developed and justified in [23].

An iterative projection method for solving linear and nonlinear hypersingular integral equations, and polyhyperpersingular and multidimensional hypersingular equations, was proposed in [24].

In [22] the unique solvability of hypersingular integral equations with even singularities (*p* = 2, 4, . . . ) was proven. Meanwhile the convergence of approximate solution to the exact one was not justified. In [24] a unique solvability of the spline-collocation method was proven. In addition, for hypersingular integral equations with bounded right-hand sides the convergence of an approximate solution sequence to the exact solution was proven under certain additional conditions.

The iterative projection method proposed here overcomes these limitations. It was shown that if the exact equation has a solution for large enough N, where N is the dimension of an approximate system of equations, an approximate solution converges to the exact one.

Hypersingular integral equations with bounded right-hand sides are a small subset of the hypersingular integral equations family. Therefore, the problem arises of constructing and justifying approximate methods for solutions for hypersingular integral equations with non-Riemann integrable functions on the right-hand sides. This paper is devoted to those issues.

A large number of works are devoted to approximate methods for solving hypersingular integral equations of the first kind

$$\int\_{-1}^{1} \frac{\mathbf{x}(\tau)d\tau}{(\tau - t)^2} + \int\_{-1}^{1} h(t, \tau)d\tau = f(t). \tag{1}$$

To solve the Equation (1), collocation and mechanical quadrature methods [17,18], the method of orthogonal polynomials [25], the method of discrete vortices [19], the method of homotopy [26] and others are used.

In the works [27–29] computational schemes for the approximate solution of the Equation (1) are constructed and their justification is carried out under the assumption that the solution has the forms *x*(*t*)=(1 − *t* <sup>2</sup>)±1/2*ω*(*t*) or *<sup>x</sup>*(*t*) = ((<sup>1</sup> <sup>−</sup> *<sup>t</sup>*)/(<sup>1</sup> <sup>+</sup> *<sup>t</sup>*))±1/2*ω*(*t*), where *<sup>ω</sup>*(*t*) is a smooth function.

The hypersingular integral equations

$$\frac{1}{\pi} \int\_{-1}^{1} \frac{\mathbf{x}(\tau) d\tau}{(\tau - t)^2} = f(t), \quad -1 < t < 1,\tag{2}$$

are widely used in aerodinamical problems and in the theory of antennae [30,31]. In the works [30,31] the Equation (2) is investigated under the assumption that the right-hand side has the form *f*(*t*) = 1/(*t* − *c*) or *f*(*t*) = *δ*(*t* − *c*), where *δ*(*t*) is the delta-function. An analytical solution of the Equation (2) with the indicated right-hand sides is obtained under the assumption that it has the form *<sup>x</sup>*(*t*) = <sup>√</sup> 1 − *t*2*ϕ*(*t*).

A fairly detailed review of analytical and numerical methods for solving hypersingular integral equations is given in [32].

In this paper, we propose an approach to solving linear and nonlinear hypersingular integral equations, the right parts of which contain functions with power features.

In particular, the right-hand sides of the form

$$f(t) = g(t)\frac{1}{t - c\_1}\frac{1}{t - c\_2}\cdots\frac{1}{t - c\_l}, \\ l = 1, 2, \ldots, \ -1 < c\_1 < \cdots < c\_l < 1,\tag{3}$$

are considered. Here *g*(*t*) is a smooth function.

Below, for simplicity of notation, we put *l* = 1 in (3).

**Remark 1.** *It can be shown that if in the hypersingular integral Equation (1) of the first kind the right side f*(*t*) ∈ *H, H is a Holder class, then the solution to this equation has the form x*(*t*)=(1 − *t* <sup>2</sup>)±1/2 *or <sup>x</sup>*(*t*) = ((<sup>1</sup> <sup>+</sup> *<sup>t</sup>*)/(<sup>1</sup> <sup>−</sup> *<sup>t</sup>*))±1/2*. For singular right-hand sides, the classes of solutions of (1) are unknown.*

*Below, when constructing and justifying the computational method, we assume that the Equation (1) with a given right-hand side has a unique solution.*

The proposed method has the following advantages:


The paper is organized as follows. The continuous method for linear and nonlinear operator equations is explained in Section 2. The numerical method for solving hypersingular integral equations is presented in Section 2.

#### **2. Continuous Method and Its Convergence Properties**

The method we employ in the next section for solving hypersingular integral equations is based on the continuous method introduced in [33].

#### *Continuous Method for Solving Operator Equations*

The continuous method for solving operator equations is based on the Lyapunov theory of stability.

Let *x*(*t*) be a solution of the differential equation in a Banach space *B*

$$\frac{d\mathbf{x}}{dt} = F(\mathbf{t}, \mathbf{x})\tag{4}$$

which is defined for all *t* ≥ *t*0. The solution *x*(*t*) is said to be stable if (i) for each *ε* > 0 there is a corresponding *δ* = *δ*(*ε*) > 0 such that any solution *x*˜(*t*) of (4) which satisfies the inequality |*x*˜(*t*0) − *x*(*t*0)| < *δ* exists and satisfies the inequality |*x*˜(*t*) − *x*(*t*)| < *ε* for all *t* ≥ *t*0.

It is said to be asymptotically stable if in addition (ii) |*x*˜(*t*) − *x*(*t*)| → 0 if *t* → ∞ whenever |*x*˜(*t*0) − *x*(*t*0)| is sufficiently small.

We will use the following notation:

*B*(*a*,*r*) = {*z* ∈ *B* : *z* − *a* ≤ *r*}, *S*(*a*,*r*) = {*z* ∈ *B* : *z* − *a* = *r*}, *Re*(*K*) = (*K*)=(*K* + *K*∗)/2, Λ(*K*) = lim *h*↓0 ( *<sup>I</sup>* <sup>+</sup> *hK*<sup>−</sup> <sup>1</sup>)*h*−1.

Here *B* is a Banach space, *a* ∈ *B*, *K* is a linear operator on *B*, Λ(*K*) is the logarithmic norm [34] of the operator *K*, *K*∗ is the conjugate operator to *K* and *I* is the identity operator.

The analytical expressions for logarithmic norms are known for operators in many spaces. We restrict ourselves to a description of the three norms.

Let *A* = {*aij*}, *i*, *j* = 1, 2, . . . , *n*, be a matrix.

In the *n*-dimensional space *Rn* of vectors *x* = (*x*1,..., *xn*) the following norms are often used:

$$\begin{aligned} \text{totalhedral- } \|\|\mathbf{x}\|\|\_{1} &= \sum\_{i=1}^{n} |\mathbf{x}\_{i}|; \text{ cubic- } \|\|\mathbf{x}\|\|\_{2} = \max\_{1 \le i \le n} |\mathbf{x}\_{i}|;\\ \text{spherical (Euclidean)- } \|\|\mathbf{x}\|\|\_{3} &= (\sum\_{i=1}^{n} \mathbf{x}\_{i}^{2})^{1/2}. \end{aligned}$$

Here are analytical expressions of the logarithmic norm of *n* × *n* matrix *A* = (*aij*), due to the above norms of the vectors:

octahedral logarithmic norm Λ<sup>1</sup>

$$\Lambda\_1(A) = \max\_{1 \le j \le n} \left( a\_{j\bar{j}} + \sum\_{i \ne j} |a\_{i\bar{j}}| \right);$$

cubic logarithmic norm Λ<sup>2</sup>

$$\Lambda\_2(\mathcal{A}) = \max\_{1 \le i \le n} \left( a\_{ii} + \sum\_{j \ne i} |a\_{ij}| \right);$$

spherical (Euclidean) logarithmic norm Λ<sup>3</sup>

$$
\Lambda\_3(A) = \lambda\_{\text{max}} \left( \frac{A + A^\*}{2} \right) \cdot \Lambda
$$

where *A*∗ is the conjugate matrix for *A*.

Note that the logarithmic norm of the same matrix can be positive in one space and negative in another.

The logarithmic norm has the some properties which are very useful for numerical mathematics.

Let *A*, *B* be *n* × *n* matrices with complex elements; and *x* = (*x*1, ... , *xn*), *y* = (*y*1, ... , *yn*), *ξ* = (*ξ*1,..., *ξn*) and *η* = (*η*1,..., *ηn*) are *n*-dimensional vectors with complex components. Let the systems of algebraic equations *Ax* = *ξ* and *By* = *η* be given. The norm of a vector and its subordinate operator norm of the matrix are agreed upon; the logarithmic norm Λ(*A*) corresponds to the operator norm.

**Theorem 1** ([35])**.** *If* <sup>Λ</sup>(*A*) <sup>&</sup>lt; <sup>0</sup>*, the matrix A is non-singular and <sup>A</sup>*−1<sup>≤</sup> 1/|Λ(*A*)|*.*

**Theorem 2** ([35])**.** *Let Ax* = *ξ, By* = *η and* Λ(*A*) < 0*,* Λ(*B*) < 0*. Then*

$$||x - y|| \le \frac{||\xi - \eta||}{|\Lambda(B)|} + \frac{||A - B||}{|\Lambda(A)\Lambda(B)|}.$$

Some properties of the logarithmic norm in a Banach space, which are useful in numerical mathematics, are given in [34].

Let us consider in a Banach space *B*, the Cauchy problem

$$\frac{d\mathbf{x}(t)}{dt} = A(\mathbf{x}(t)),\tag{5}$$

$$\mathbf{x}(0) = \mathbf{x}\_0.\tag{6}$$

Let us assume that the nonlinear operator *A* has a Frechet derivative and *A*(0) = 0.

The sufficiently satisfying conditions of asymptotically stability for the solution of the Cauchy problem (5), (6) were obtained in [36,37].

**Theorem 3.** *Let the integral <sup>t</sup>* 0 Λ(*A* (*ϕ*(*τ*)))*dτ be non-positive (respectively, be negative and satisfy* lim *t*→+∞ 1 *t t* 0 Λ(*A* (*ϕ*(*τ*)))*dτ* ≤ −*αϕ*, *αϕ* > 0*) for any differentiable curve ϕ*(*t*) *lying in a ball B*(0,*r*) *of some radius r. Then the trivial solution of Equation (5) is stable (respectively, asymptotically stable).*

**Remark 2.** *Additionally, the Theorem is valid under r* = ∞*.*

Let us consider in a Banach space *B* a nonlinear operator equation

$$A(\mathbf{x}) - f = 0,\tag{7}$$

where operator *A* acts from *B* into *B*.

We associate Equation (7) with the Cauchy problem

$$\frac{d\mathbf{x}(t)}{dt} = A(\mathbf{x}(t)) - f\_{\prime} \tag{8}$$

$$
\mathfrak{x}(0) = \mathfrak{x}\_0. \tag{9}
$$

Let *x*∗ be a the solution of Equation (7). Let us make the change of variable *x* = *x*∗ + *v* in Equation (8). This change reduces the Cauchy problem (8), (9) to the form

$$\frac{dv(t)}{dt} = A(\mathbf{x}^\* + v(t)) - A(\mathbf{x}^\*),\tag{10}$$

$$
\upsilon(0) = \mathbf{x}\_0 - \mathbf{x}^\*.\tag{11}
$$

It is easy to see that if the trivial solution of Equation (10) is globally asymptotically stable, then lim*t*→+<sup>∞</sup> *v*(*t*) → 0. So, for any initial value the solution of the Cauchy problem (8), (9) tends to *x*∗. It follows from the next assertions which were proven in [33].

**Theorem 4.** *Let Equation*(7) *have a solution x*∗*, and let inequality*

$$\lim\_{t\_1 \to +\infty} \frac{1}{t} \int\_0^t \Lambda(A'(\mathcal{g}(\tau)))d\tau \le -a\_{\mathcal{S}'}a\_{\mathcal{S}} > 0,\tag{12}$$

*be true on each differentiable curve g*(*t*) *lying in the Banach space B. Then the solution of the Cauchy problem (8), (9) converges to the solution x*∗ *of Equation*(7) *for any initial value.*

**Theorem 5.** *Let Equation (7) have a solution x*∗*, and let the following conditions be satisfied on any differentiable curve g*(*t*) *lying in the ball B*(*x*∗,*r*)*.*

*1. The inequality*

$$\int\_{0}^{t} \Lambda(A'(\mathcal{g}(\tau)))d\tau \le 0$$

*holds for all t*(*t* > 0)*.*

*2. Inequality (12) is satisfied.*

*Then the solution of the Cauchy problem (8), (9) converges to the solution x*∗ *of Equation (7).*

**Remark 3.** *The sufficient condition for convergence of the Cauchy problem (8), (9) solution to the solution of the operator Equation (7) is given above. It was obtained by analysing Lyapunov stability. One of the first basic* *results in accretive operator theory was a relation between the solution of operator equation Au* = 0*, where A is a locally Lipschitzian and accretive operator, and the differential equation du dt* = *Au was obtained in [38].*

*Later, accretive operator theory and its applications for finding fixed points and constructing iterative procedures were studied by many authors. Basic results and a detailed bibliography devoted to the subject may be found in [39–42].*

#### **3. An Solution of Hypersingular Integral Equations with the Continuous Method**

Let us consider the method of mechanical quadrature for solving hypersingular integral equation of the types

$$a(t)\mathbf{x}(t) + \int\_{-1}^{1} \frac{h(t,\tau)\mathbf{x}(\tau)d\tau}{(\tau - t)^2} = f(t). \tag{13}$$

and

$$a(t)x(t) + \int\_{-1}^{1} \frac{h(t, \tau, x(\tau))d\tau}{(\tau - t)^2} = f(t). \tag{14}$$

It is assumed that in the Equations (13) and (14) the right-hand sides have features of the following types

$$f(t) = \sum\_{i=1}^{l} g\_i(t) \frac{1}{t - c\_i}, \quad f(t) = \mathcal{g}(t) \prod\_{i=1}^{l} \frac{1}{t - c\_i}.$$

where −1 < *ci* < 1, *i* = 1, 2, . . . , *l*, *l* = 1, 2, . . . ; *g*(*t*), *gi*(*t*), *i* = 1, 2, ··· , *l*,—are continuous functions.

In what follows, without loss of generality, we set *l* = 1.

Let us recall the Hadamard definition of hypersingular integrals [43].

**Definition 1** ([43])**.** *The integral of the type*

$$\int\_{a}^{b} \frac{A(\mathbf{x}) \, d\mathbf{x}}{(b-\mathbf{x})^{p+\alpha}}$$

*for an integer p and* 0 < *α* < 1*, is defined as*

$$\lim\_{x \to b^{-}} \left[ \int\_{a}^{x} \frac{A(t) \, dt}{(b-t)^{p+\alpha}} + \frac{B(x)}{(b-x)^{p+\alpha-1}} \right],$$

*if A*(*x*) *has p derivatives in the neighborhood of the point b. Here B*(*x*) *is any function that satisfies the following two conditions:*


It is easy to see [43], that the conditions (i) and (ii) are sufficient for the existence of the limit. Chikin in [44] introduced the following definition.

**Definition 2** ([44])**.** *The Cauchy–Hadamard principal value of the integral*

$$\int\_{a}^{b} \frac{\varrho(\tau) \, d\tau}{(\tau - c)^{p}}, \quad a < c < b,\tag{15}$$

*is defined as*

$$\int\_{a}^{b} \frac{\varrho(\tau) \, d\tau}{(\tau - c)^{p}} = \lim\_{v \to 0} \left[ \int\_{a}^{c - v} \frac{\varrho(\tau) \, d\tau}{(\tau - c)^{p}} + \int\_{c + v}^{b} \frac{\varrho(\tau) \, d\tau}{(\tau - c)^{p}} + \frac{\xi(v)}{v^{p - 1}} \right],$$

*where ξ*(*v*) *is a function constructed so that the limit exists.*

#### *3.1. An Approximate Solution of Linear Hypersingular Integral Equations with Second Order Singularity*

Consider a one-dimensional hypersingular integral equation of the type

$$\mathbf{Kx} \equiv a(t)\mathbf{x}(t) + \int\_{-1}^{1} \frac{h(t,\tau)\mathbf{x}(\tau)d\tau}{(\tau - t)^2} = f(t),\tag{16}$$

where *f*(*t*) = *g*(*t*)/(*t* − *c*) or *f*(*t*) = *g*(*t*)/((1 − *t* <sup>2</sup>)(*<sup>t</sup>* <sup>−</sup> *<sup>c</sup>*)), <sup>−</sup><sup>1</sup> <sup>&</sup>lt; *<sup>c</sup>* <sup>&</sup>lt; 1, *<sup>g</sup>*(*t*) <sup>∈</sup> *<sup>C</sup>*[−1, 1]. Divide the interval [−1, 1] into two subintervals [−1, *c*], [*c*, 1].

Let us fix a positive integer *N*0. Put *h* = 2/*N*0, *N*<sup>1</sup> = &(1 + *c*)/*h*', *N*<sup>2</sup> = &(1 − *c*)/*h*', *N* = *N*<sup>1</sup> + *N*2. Divide the interval [−1, *c*] into *N*<sup>1</sup> subintervals at the points *tk* = −1 + (*c* + 1)*k*/*N*1, *k* = 0, 1, . . . , *N*1.

Divide the interval [*c*, 1] into *N*<sup>2</sup> subintervals at the points *τ<sup>k</sup>* = *c* + (1 − *c*)*k*/*N*2, *k* = 0, 1, . . . , *N*2. Let us introduce the nodes ¯*t*<sup>0</sup> <sup>=</sup> *<sup>t</sup>*<sup>0</sup> <sup>+</sup> 1/2(*N*1)2, ¯*tk* <sup>=</sup> *tk*, *<sup>k</sup>* <sup>=</sup> 1, 2, ... , *<sup>N</sup>*<sup>1</sup> <sup>−</sup> 1, ¯*tN*<sup>1</sup> <sup>=</sup> *tN*<sup>1</sup> <sup>−</sup> 1/2(*N*1)2; *<sup>τ</sup>*¯0 <sup>=</sup> *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> 1/2(*N*2)2, *<sup>τ</sup>*¯*<sup>k</sup>* <sup>=</sup> *<sup>τ</sup>k*, *<sup>k</sup>* <sup>=</sup> 1, 2, . . . , *<sup>N</sup>*<sup>2</sup> <sup>−</sup> 1, *<sup>τ</sup>*¯*N*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> 1/2(*N*2)2.

As an approximate solution of (16), we shall seek in the form of a continuous function

$$\alpha\_N(t) = \sum\_{k=0}^{N\_1} \alpha\_k \varphi\_k(t) + \sum\_{k=0}^{N\_2} \beta\_k \psi\_k(t),\tag{17}$$

where *ϕk*(*t*), *k* = 0, 1, . . . , *N*1, *ψk*(*t*), *k* = 0, 1, . . . , *N*<sup>2</sup> is a family of basis functions.

For nodes *tk*, *k* = 1, . . . , *N*<sup>1</sup> − 1, the corresponding basis elements are determined by

$$\rho\_{k}(t) = \begin{cases} 0, & t\_{k-1} \le t \le t\_{k-1} + \frac{1}{N\_{1}^{2}}, \\ \frac{N\_{1}^{2}}{(1+\varepsilon)N\_{1} - 2}(t - t\_{k-1}) - \frac{1}{(1+\varepsilon)N\_{1} - 2}, & t\_{k-1} + \frac{1}{N\_{1}^{2}} \le t \le t\_{k} - \frac{1}{N\_{1}^{2}}, \\ 1, & t\_{k} - \frac{1}{N\_{1}^{2}} \le t \le t\_{k} + \frac{1}{N\_{1}^{2}}, \\ -\frac{N\_{1}^{2}}{(1+\varepsilon)N\_{1} - 2}(t - t\_{k} - \frac{1}{N\_{1}^{2}}) + 1, & t\_{k} + \frac{1}{N\_{1}^{2}} \le t \le t\_{k+1} - \frac{1}{N\_{1}^{2}}, \\ 0, & t\_{k+1} - \frac{1}{N\_{1}^{2}} \le t \le t\_{k+1}, \\ 0, & t \in [-1, 1] \cup [t\_{k-1}, t\_{k+1}]. \end{cases} \tag{18}$$

For boundary nodes *tk*, *k* = 0 and *k* = *N*<sup>1</sup> the corresponding basis elements are defined as

$$\varphi\_0(t) = \begin{cases} 1, & -1 \le t \le -1 + \frac{1}{N\_1^{2^\prime}} \\ -\frac{N\_1^2}{(1+c)N\_1 - 2}(t+1 - \frac{1}{N\_1^2}) + 1, & -1 + \frac{1}{N\_1^2} \le t \le t\_1 - \frac{1}{N\_1^{2^\prime}} \\ 0, & t\_1 - \frac{1}{N\_1^2} \le t \le t\_1 \\ 0, & [-1, 1] \backslash [t\_0, t\_1]; \end{cases} \tag{19}$$

and

$$\varphi\_{N\_{\mathbb{I}}}(t) = \begin{cases} 0, & -1 \le t \le t\_{N\_{\mathbb{I}}-1} + \frac{1}{N\_{\mathbb{I}}^{2}},\\ \frac{N\_{\mathbb{I}}^{2}}{(1+c)N\_{\mathbb{I}}-2}(t - t\_{N\_{\mathbb{I}}-1}) - \frac{1}{(1+c)N\_{\mathbb{I}}-2'}, & t\_{N\_{\mathbb{I}}-1} + \frac{1}{N\_{\mathbb{I}}^{2}} \le t \le c - \frac{1}{N\_{\mathbb{I}}^{2}},\\ 1, & c - \frac{1}{N\_{\mathbb{I}}^{2}} \le t \le c. \end{cases} \tag{20}$$

For nodes *τk*, *k* = 0, 1, ... , *N*2, the corresponding basis elements *ψk*, *k* = 0, 1, ... , *N*2, are determined in a the similar way: For nodes *τk*, *k* = 1, ... , *N*<sup>2</sup> − 1, the corresponding basis elements are determined by

$$\psi\_{k}(t) = \begin{cases} 0, & \mathfrak{r}\_{k-1} \le t \le \mathfrak{r}\_{k-1} + \frac{1}{N\_{2}^{2}}, \\ \frac{N\_{2}^{2}}{(1-\varepsilon)N\_{2}-2}(t-\mathfrak{r}\_{k-1}) - \frac{1}{(1-\varepsilon)N\_{2}-2}, & \mathfrak{r}\_{k-1} + \frac{1}{N\_{2}^{2}} \le t \le \mathfrak{r}\_{k} - \frac{1}{N\_{2}^{2}}, \\ 1, & \mathfrak{r}\_{k} - \frac{1}{N\_{2}^{2}} \le t \le \mathfrak{r}\_{k} + \frac{1}{N\_{2}^{2}}, \\ -\frac{N\_{2}^{2}}{(1-\varepsilon)N\_{2}-2}(t-\mathfrak{r}\_{k} - \frac{1}{N\_{2}^{2}}) + 1, & \mathfrak{r}\_{k} + \frac{1}{N\_{2}^{2}} \le t \le \mathfrak{r}\_{k+1} - \frac{1}{N\_{2}^{2}}, \\ 0, & \mathfrak{r}\_{k+1} - \frac{1}{N\_{2}^{2}} \le t \le \mathfrak{r}\_{k+1}, \\ 0, & t \in [-1,1] \cup [\mathfrak{r}\_{k-1}, \mathfrak{r}\_{k+1}]. \end{cases} \tag{21}$$

For boundary nodes *τk*, *k* = 0 and *k* = *N*<sup>2</sup> the corresponding basis elements are defined as

$$\psi\_0(t) = \begin{cases} 1, & \mathcal{c} \le t \le \mathcal{c} + \frac{1}{N\_1^2}, \\ -\frac{N\_2^2}{(1-\varepsilon)N\_2 - 2}(t - \mathcal{c} - \frac{1}{N\_2^2}) + 1, & \mathcal{c} + \frac{1}{N\_2^2} \le t \le \tau\_1 - \frac{1}{N\_2^2}, \\ 0, & \tau\_1 - \frac{1}{N\_2^2} \le t \le \tau\_1, \\ 0, & [-1, 1] \backslash [\mathcal{c}, \tau\_1]; \end{cases} \tag{22}$$

and

$$\psi\_{\rm N\_2}(t) = \begin{cases} 0, & -1 \le t \le \tau\_{\rm N\_2 - 1} + \frac{1}{N\_2^2} \\ \frac{N\_2^2}{(1 - \varepsilon)N\_2 - 2}(t - \tau\_{\rm N\_2 - 1}) - \frac{1}{(1 - \varepsilon)N\_2 - 2}, & \tau\_{\rm N\_2 - 1} + \frac{1}{N\_2^2} \le t \le 1 - \frac{1}{N\_2^2} \\ 1, & 1 - \frac{1}{N\_2^2} \le t \le 1. \end{cases} \tag{23}$$

To simplify the description of computational scheme, we introduce the following notation:


Here *vi* = *ti*, *i* = 0, 1, . . . , *N*<sup>1</sup> , *vN*1+*<sup>i</sup>* = *τi*, *i* = 1, 2, . . . , *N*2,

$$\begin{array}{l} \gamma\_{i} = \alpha\_{i}, i = 0, 1, \ldots, N\_{1} \text{ }, \gamma\_{N\_{1}+1+i} = \beta\_{i}, i = 0, 1, \ldots, N\_{2}, \\\zeta\_{i} = \wp\_{i}, i = 0, 1, \ldots, N\_{1} \text{ }, \zeta\_{N\_{1}+1+i} = \psi\_{i}, i = 0, 1, \ldots, N\_{2}. \end{array}$$

Let us recall that the points *tN*<sup>1</sup> and *τ*<sup>0</sup> coincide.

Applying the collocation method on the knots *v*¯*k*, *k* = 0, 1, ... , *N*<sup>∗</sup> + 1 to the Equation (16), we obtain the following system of algebraic equations for finding unknown coefficients {*γk*} of the polygon (17)

$$a(\psi\_k)\gamma\_k + \sum\_{l=0}^{N^\*+1} h(\psi\_k, \upsilon\_l)\gamma\_l \int\_{-1}^1 \frac{\zeta\_l(\tau)}{(\tau - \vartheta\_k)^2} d\tau = f(\vartheta\_k), \tag{24}$$

*k* = 0, 1, . . . , *N*∗ + 1.

Using the definition of hypersingular integrals, we receive:

$$\int\_{w\_{k-1}}^{v\_{k+1}} \frac{\zeta\_k(\tau)d\tau}{(\tau - \vartheta\_k)^2} = -2 \frac{N\_1^2}{(1+c)N\_1 - 2} \ln((1+c)N\_1 - 1), k = 1, 2, \dots, N\_1 - 1; \tag{25}$$

$$\int\_{\psi\_{k-1}}^{\psi\_{k+1}} \frac{\zeta\_{k+1}(\tau)d\tau}{(\tau - \vartheta\_k)^2} = -2 \frac{N\_2^2}{(1-c)N\_2 - 2} \ln((1-c)N\_2 - 1), k = N\_1 + 2, \dots, N^\* - 1; \tag{26}$$

$$\int\_{-1}^{v\_1} \frac{\zeta\_0(\tau)d\tau}{(\tau+1-\frac{1}{2N\_1^2})^2} = -2N\_1^2 - N\_1^2 \frac{\ln(2(c+1)N\_1 - 3)}{(c+1)N\_1 - 2},\tag{27}$$

$$\int\_{\tau\_{N\_1-1}}^{\tau\_{N\_1}} \frac{\zeta\_{N\_1}(\tau)d\tau}{(\tau - \upsilon\_{N\_1} + \frac{1}{2N\_1^2})^2} = -2N\_1^2 - N\_1^2 \frac{\ln(2(c+1)N\_1 - 3)}{(c+1)N\_1 - 2},\tag{28}$$

$$\int\_{v\_{N\_1}}^{v\_{N\_1+1}} \frac{\zeta\_{N\_1+1}(\tau)d\tau}{(\tau - v\_{N\_1} - \frac{1}{2N\_2^2})^2} = -2N\_2^2 - N\_2^2 \frac{\ln(2(1-c)N\_2 - 3)}{(1-c)N\_2 - 2},\tag{29}$$

$$\int\_{\tau\_{N^\*-1}}^1 \frac{\zeta\_{N^\*+1}(\tau)d\tau}{(\tau - 1 + \frac{1}{2N^2})^2} = -2N\_2^2 - N\_2^2 \frac{\ln(2(1-c)N\_2 - 3)}{(1-c)N\_2 - 2},\tag{30}$$

$$\int\_{-1}^{1} \left[ \sum\_{l=1}^{N^\*+1} \zeta\_l(\tau) \right] \frac{d\tau}{(\tau + 1 - \frac{1}{2N\_1^2})^2} = -\frac{2N\_1^2}{4N\_1^2 - 1} + N\_1^2 \frac{\ln(2(1+c)N\_1 - 3)}{(1+c)N\_1 - 2},\tag{31}$$

$$\int\_{-1}^{1} \left[ \sum\_{l=0}^{N^\*} \zeta\_l(\tau) \right] \frac{d\tau}{(\tau - 1 + \frac{1}{2N\_2^2})^2} = -\frac{2N\_2^2}{4N\_2^2 - 1} + N\_2^2 \frac{\ln(2(1-c)N\_2 - 3)}{(1-c)N\_2 - 2},\tag{32}$$

$$\int\_{-1}^{1} \left[ \sum\_{l=0}^{N^{\prime}+1} {\ulcorner \zeta\_l(\tau)} \right] \frac{d\tau}{(\tau - (c - \frac{1}{2N\_1^2}))^2} = -\frac{2N\_1^2}{4N\_1^2 - 1} + N\_1^2 \frac{\ln(2(1+c)N\_1 - 3)}{(1+c)N\_1 - 2},\tag{33}$$

$$\int\_{-1}^{1} \left[ \sum\_{l=0}^{N^\*+1} \, \prescript{\prime \chi}{}{\zeta}\_l(\tau) \right] \frac{d\tau}{(\tau - (c + \frac{1}{2N\_2^2}))^2} = -\frac{2N\_2^2}{4N\_2^2 - 1} + N\_2^2 \frac{\ln(2(1+c)N\_2 - 3)}{(1+c)N\_2 - 2},\tag{34}$$

$$\frac{1}{\tau\_{-1}} \left[ \sum\_{l=0}^{N^\*+1} {N^\*} \varphi\_l(\mathbf{r}) \right] \frac{d\tau}{(\tau - v\_k)^2} = -\frac{N\_1}{(c+1)k} - \frac{N\_1}{2N\_1 - (c+1)k} + \frac{2N\_1^2}{(1+c)N\_1 - 2} \ln((c+1)N\_1 - 1);\tag{35}$$

$$\frac{1}{\varepsilon} \int\_{-1}^{\frac{1}{2}} \left[ \sum\_{l=0}^{N\_1^\*+1} \mathcal{W}\_l(\mathbf{r}) \right] \frac{d\mathbf{r}}{(\mathbf{r} - v\_k)^2} = -\frac{N\_2}{(1-\varepsilon)k} - \frac{N\_2}{2N\_2 - (1-\varepsilon)k} + \frac{2N\_2^2}{(1-\varepsilon)N\_2 - 2} \ln((1-\varepsilon)N\_2 - 1). \tag{36}$$

Here ∑ *l* , ∑ *<sup>l</sup>* , ∑ *<sup>l</sup>* , ∑ *<sup>l</sup>* indicates a summation over *l* = *N*1, *l* = *N*<sup>1</sup> + 1, *l* = *k*(1 ≤ *k* ≤ *N*<sup>1</sup> − 1), *l* = *k*(*N*<sup>1</sup> + 2 ≤ *k* ≤ *N*<sup>∗</sup> − 1), respectively. Detailed calculations are given in [23].

We can rewrite the system (24) as

$$\begin{split} a(\vec{v}\_k)\gamma\_k - h(\vec{v}\_k, \vec{v}\_k) 2N\_1^2 \frac{\ln(N\_l - 1)}{(1 + c)N\_l - 2} \gamma\_k + \sum\_{l=0}^{N\_r - 1} {}^{\prime}\gamma\_l h(\vec{v}\_k, \vec{v}\_l) \int\_{-1}^1 \zeta\_l(\tau) \frac{d\tau}{(\tau - \vartheta\_k)^2} \\ = f(\vec{v}\_k), \quad k = 1, \dots, N\_1 - 1; \end{split} \tag{37}$$

$$\begin{split} a(\boldsymbol{\upbeta}\_{k})\gamma\_{\mathbf{k}} - h(\boldsymbol{\upbeta}\_{k}, \boldsymbol{\upbeta}\_{k})2N\_{2}^{2} \frac{\ln(N\_{2}-1)}{(1-c)N\_{2}-2} \gamma\_{\mathbf{k}} + \sum\_{l=0}^{N^{\*}+1} \, ^{\prime} \gamma\_{l} h(\boldsymbol{\upbeta}\_{k}, \boldsymbol{\upbeta}\_{l}) \, \stackrel{\scriptstyle 1}{\int} \boldsymbol{\upzeta}\_{l}(\boldsymbol{\uppi}) \frac{d\boldsymbol{\uppi}}{(\boldsymbol{\uppi} - \boldsymbol{\uppi}\_{k})^{2}} \\ = f(\boldsymbol{\upbeta}\_{k}), \quad k = N\_{1} + 2, \ldots, N^{\*}; \end{split} \tag{38}$$

$$\begin{split} &a(\boldsymbol{\vartheta}\_{0})\gamma\_{0} - h(\boldsymbol{\vartheta}\_{0}, \boldsymbol{\vartheta}\_{0})(2N\_{1}^{2} + N\_{1}^{2}\frac{\ln(2(1+\boldsymbol{c})N\_{1} - 3)}{(1+\boldsymbol{c})N\_{1} - 2})\gamma\_{0} \\ &+\sum\_{l=1}^{N^{\*}+1} \gamma\_{l}h(\boldsymbol{\vartheta}\_{k}, \boldsymbol{\vartheta}\_{l})\int\_{-1}^{1} \boldsymbol{\zeta}\_{l}(\boldsymbol{\tau})\frac{d\boldsymbol{\tau}}{(\boldsymbol{\tau}-\boldsymbol{\tau}\_{0})^{2}} = f(\boldsymbol{\vartheta}\_{0}); \end{split} \tag{39}$$

$$\begin{split} a(\boldsymbol{\vartheta}\_{\mathcal{N}\_{1}})\boldsymbol{\gamma}\_{\mathcal{N}\_{1}} - h(\boldsymbol{\vartheta}\_{\mathcal{N}\_{1}}, \boldsymbol{\vartheta}\_{\mathcal{N}\_{1}}) (2\boldsymbol{N}\_{1}^{2} + \boldsymbol{N}\_{1}^{2} \frac{\ln(2(1+\boldsymbol{c})\boldsymbol{N}\_{1} - 3)}{(1+\boldsymbol{c})\boldsymbol{N}\_{1} - 3}) \boldsymbol{\gamma}\_{\mathcal{N}\_{1}} \\ + \sum\_{l=0}^{N^{\*}\_{1}+1} \prescript{}{\prime\prime}{\boldsymbol{\gamma}}\_{l} h(\boldsymbol{\vartheta}\_{\mathcal{N}\_{1}}, \boldsymbol{\vartheta}\_{l}) \int\_{-1}^{1} \boldsymbol{\zeta}\_{l}(\boldsymbol{\tau}) \frac{d\boldsymbol{\tau}}{(\boldsymbol{\tau} - \boldsymbol{\uprho}\_{\mathcal{N}\_{1}})^{2}} = f(\boldsymbol{\vartheta}\_{\mathcal{N}\_{1}}); \end{split} \tag{40}$$

$$\begin{split} a(\bar{\boldsymbol{v}}\_{\text{N}\_{1}+1})\boldsymbol{\gamma}\_{\text{N}\_{1}+1} - h(\bar{\boldsymbol{v}}\_{\text{N}\_{1}+1}, \bar{\boldsymbol{v}}\_{\text{N}\_{1}+1})(2\mathbf{N}\_{2}^{2} + \mathbf{N}\_{2}^{2}\frac{\ln(2(1-\boldsymbol{c})\boldsymbol{N}\_{2} - 3)}{(1-\boldsymbol{c})\boldsymbol{N}\_{2} - 2})\boldsymbol{\gamma}\_{\text{N}\_{1}+1} \\ + \sum\_{l=0}^{N^{\*}\_{1}+1} \prime \gamma\_{l} h(\bar{\boldsymbol{v}}\_{\text{N}\_{1}+1}, \bar{\boldsymbol{v}}\_{l}) \int\_{-1}^{1} \zeta\_{l}(\boldsymbol{\tau}) \frac{d\boldsymbol{\tau}}{(\boldsymbol{\tau} - \bar{\boldsymbol{v}}\_{\text{N}\_{1}+1})^{2}} = f(\bar{\boldsymbol{v}}\_{\text{N}\_{1}}); \end{split} \tag{41}$$

$$\begin{split} \ln \left( \overline{\upsilon}\_{N^\*+1} \right) \gamma\_{N^\*+1} - h \left( \overline{\upsilon}\_{N^\*+1}, \overline{\upsilon}\_{N^\*+1} \right) \left( 2N\_2^2 + N\_2^2 \frac{\ln(2(1-\varepsilon)N\_2 - 3)}{(1-\varepsilon)N\_2 - 2} \right) \gamma\_{N^\*+1} \\ + \sum\_{l=0}^{N^\*} \gamma\_l h (\overline{\upsilon}\_{N^\*+1}, \overline{\upsilon}\_l) \int\_{-1}^1 \overline{\zeta}\_l(\tau) \frac{d\tau}{(\tau - \overline{\upsilon}\_{N^\*-1})^2} = f(\overline{\upsilon}\_{N^\*+1}). \end{split} \tag{42}$$

Here ∑ , ∑ , ∑ indicates a summation over *l* = *k*, *l* = *N*1, *l* = *N*<sup>1</sup> + 1, respectively. The system (37)–(42) is equivalent to the system

$$\begin{aligned} \left(\operatorname{sgn} h(\overline{v}\_k, \overline{v}\_k)\right) \left(a(\overline{v}\_k)\gamma\_k - h(\overline{v}\_k, \overline{v}\_k)2N\_1^2 \frac{\ln(N\_1 - 1)}{(1 + c)N\_1 - 2} \gamma\_k \\ + \sum\_{l=0}^{N^\* + 1} \,' \gamma\_l h(\overline{v}\_k, \overline{v}\_l) \int\_{-1}^1 \zeta\_l(\tau) \frac{d\tau}{(\tau - \overline{v}\_k)^2} \right) = (\operatorname{sgn} h(t\_k, t\_k)) f(v\_k), k = 1, \dots, N\_1 - 1; \end{aligned} \tag{43}$$

$$\begin{aligned} \left(\operatorname{sgn} h(\boldsymbol{\psi}\_{k}, \boldsymbol{\psi}\_{k})\right) \left(a(\boldsymbol{\psi}\_{k})\gamma\_{k} - h(\boldsymbol{\psi}\_{k}, \boldsymbol{\psi}\_{k})2N\_{2}^{2}\frac{\ln(N\_{2}-1)}{(1-\varepsilon)N\_{2}-2}\gamma\_{k} \\ + \sum\_{l=0}^{N^{\*}+1} \gamma\_{l}h(\boldsymbol{\psi}\_{k}, \boldsymbol{\psi}\_{l}) \int\_{-1}^{1} \zeta\_{l}(\boldsymbol{\tau}) \frac{d\boldsymbol{\tau}}{(\boldsymbol{\tau}-\boldsymbol{\bar{v}}\_{k})^{2}}\right) &= (\operatorname{sgn} h(t\_{k}, t\_{k})) f(\boldsymbol{\psi}\_{k}), \quad k = N\_{1}+2, \ldots, N^{\*}; \end{aligned} \tag{44}$$

$$\begin{aligned} \left(\operatorname{sgn} h(\overline{v}\_{0}, \overline{v}\_{0})\right) \left(a(\overline{v}\_{0})\gamma\_{0} - h(\overline{v}\_{0}, \overline{v}\_{0})\left(2N\_{1}^{2} + N\_{1}^{2} \frac{\ln(2(1+\varepsilon)N\_{1} - 3)}{(1+\varepsilon)N\_{1} - 2}\right)\gamma\_{0} \\ + \sum\_{l=1}^{N^{\*}+1} \gamma\_{l} h(\overline{v}\_{k}, \overline{v}\_{l}) \int\_{-1}^{1} \zeta\_{l}(\tau) \frac{d\overline{\tau}}{(\overline{\tau} - \overline{v}\_{0})^{2}}\right) &= \left(\operatorname{sgn} h(v\_{0}, v\_{0})\right) f(\overline{v}\_{0}); \end{aligned} \tag{45}$$

$$\begin{split} \left( \operatorname{sgn} h(\boldsymbol{\upsilon}\_{\text{N}\_{1}}, \boldsymbol{\upsilon}\_{\text{N}\_{1}}) \right) \left( \operatorname{a} (\boldsymbol{\upsilon}\_{\text{N}\_{1}}) \boldsymbol{\gamma}\_{\text{N}\_{1}} - h(\boldsymbol{\upsilon}\_{\text{N}\_{1}}, \boldsymbol{\upsilon}\_{\text{N}\_{1}}) (2N\_{1}^{2} + N\_{1}^{2} \frac{\ln(2(1+c)N\_{1} - 3)}{(1+c)N\_{1} - 2}) \boldsymbol{\gamma}\_{\text{N}\_{1}} \\ + \sum\_{l=0}^{N^{\*}+1} \boldsymbol{\eta}\_{l} h(\boldsymbol{\upsilon}\_{\text{N}\_{1}}, \boldsymbol{\upsilon}\_{l}) \int\_{-1}^{1} \boldsymbol{\zeta}\_{l} (\boldsymbol{\tau}) \frac{d\boldsymbol{\tau}}{(\boldsymbol{\tau} - \boldsymbol{\uprho}\_{\text{N}\_{1}})^{2}} \right) = (\operatorname{sgn} h(\boldsymbol{\upsilon}\_{\text{N}\_{1}}, \boldsymbol{\upsilon}\_{\text{N}\_{1}})) f(\boldsymbol{\upsilon}\_{\text{N}\_{1}}); \end{split} \tag{46}$$

$$\begin{aligned} \left(\operatorname{sgn} h(\upsilon\_{\text{N}\_1+1}, \upsilon\_{\text{N}\_1+1})\right) \left(a(\overline{\upsilon}\_{\text{N}\_1+1})\gamma\_{\text{N}\_1+1} - h(\overline{\upsilon}\_{\text{N}\_1+1}, \overline{\upsilon}\_{\text{N}\_1+1})\right) \left(2N\_2^2\right) \\ + N\_2^2 \frac{\ln(2(1-\varepsilon)N\_2 - 3)}{(1-\varepsilon)N\_2 - 2} (\gamma\_{\text{N}\_1+1} + \sum\_{l=0}^{N^\*+1} \gamma\_{l} h(\overline{\upsilon}\_{\text{N}\_1}, \overline{\upsilon}\_l) \int \zeta\_l(\tau) \frac{d\tau}{(\tau - \varepsilon\_{\text{N}\_1+1})^2} \\ = (\operatorname{sgn} h(\upsilon\_{\text{N}\_1}, \upsilon\_{\text{N}\_1})) f(\overline{\upsilon}\_{\text{N}\_1}); \end{aligned} \tag{47}$$

$$\begin{split} \left( \operatorname{sgn} h(\boldsymbol{\upsilon}\_{N^\*+1}, \boldsymbol{\upsilon}\_{N^\*+1}) \right) \left( \operatorname{a} (\boldsymbol{\overline{\upsilon}\_{N^\*+1}}) \boldsymbol{\gamma}\_{N^\*+1} - h(\boldsymbol{\overline{\upsilon}\_{N^\*+1}}, \boldsymbol{\overline{\upsilon}\_{N^\*+1}}) \right. \\ \left( 2N\_2^2 + N\_2^2 \frac{\ln(2(1-\boldsymbol{\varepsilon})N\_2 - 3)}{(1-\boldsymbol{\varepsilon})N\_2 - 2} \right) \boldsymbol{\gamma}\_{N^\*+1} + \sum\_{l=0}^{N^\*} \boldsymbol{\gamma}\_l h(\boldsymbol{\overline{\upsilon}\_{N^\*-1}}, \boldsymbol{\upsilon}\_l) \left( \operatorname{\boldsymbol{\overline{\upsilon}}\_l} (\boldsymbol{\tau}) \frac{d\boldsymbol{\tau}}{(\boldsymbol{\tau} - \boldsymbol{\underline{\upsilon}\_{N^\*+1}})^2} \right) \\ = (\operatorname{sgn} h(\boldsymbol{\upsilon}\_{N^\*+1}, \boldsymbol{\upsilon}\_{N^\*+1})) f(\boldsymbol{\overline{\upsilon}\_{N^\*+1}}). \end{split} \tag{48}$$

Here ∑ , ∑ , ∑ indicates a summation over *l* = *k*, *l* = *N*1, *l* = *N*<sup>1</sup> + 1, respectively. Let us write the system (43)–(48) in a matrix form

$$DX = F\_{\prime\prime}$$

where *D* = {*dkl*}, *k*, *l* = 0, 1, ... , *N*<sup>∗</sup> + 1, *X* = (*x*0, *x*1, ... , *xN*∗+1), *F* = (*f*0, *f*1, ... , *fN*∗+1). The values {*dkl*}, {*xk*}, and { *fk*} are obvious.

The diagonal elements in the left-hand side of the system of Equations (43)–(48) have the following forms

$$d\_{kk} = \left(\text{sgn}\,h(\vec{v}\_k, \vec{v}\_k)\right) \left(a(\vec{v}\_k) - h(\vec{v}\_k, \vec{v}\_k)2N\_1^2 \frac{\ln(N\_1 - 1)}{(1 + \varepsilon)N\_1 - 2}\right)^2$$

*k* = 1, 2, . . . , *N*<sup>1</sup> − 1,

$$d\_{kk} = \left(\text{sgn}\,h(\vec{v}\_k, \vec{v}\_k)\right) \left(a(\vec{v}\_k) - h(\vec{v}\_k, \vec{v}\_k) 2N\_1^2 \frac{\ln(N\_1 - 1)}{(1 + c)N\_1 - 2}\right),$$

*k* = *N*<sup>1</sup> + 2, . . . , *N*∗,

$$d\_{00} = \left(\text{sgn}\,h(\vec{v}\_0, \vec{v}\_0)\right)\left(a(\vec{v}\_0) - h(\vec{v}\_0, \vec{v}\_0)(2N\_1^2 + N\_1^2 \frac{\ln(2(1+c)N\_1 - 3)}{(1+c)N\_1 - 2}\right),$$

$$\begin{split} d\_{N\_1, N\_1} = \left( \operatorname{sgn} h(v\_{N\_1, \prime} v\_{N\_1}) \right) \left( a(\vartheta\_{N\_1}) - h(\vartheta\_{N\_1, \prime} \vartheta\_{N\_1}) \right) (2N\_1^2 \\ + N\_1^2 \frac{\ln \left( 2(1 + c) N\_1 - 3 \right)}{(1 + c) N\_1 - 2} \right) \end{split}$$

$$\begin{split}d\_{N\_1+1,N\_1+1} &= \left(\operatorname{sgn} h(v\_{N\_1+1}, v\_{N\_1+1})\right) \left(a(\vec{v}\_{N\_1+1}) - h(\vec{v}\_{N\_1+1}, \vec{v}\_{N\_1+1})\right) \\ &\quad \left(2N\_2^2 + N\_2^2 \frac{\ln(2(1-c)N\_2 - 3)}{(1-c)N\_2 - 2}\right) .\end{split}$$

$$\begin{aligned} d\_{N^\*+1, N^\*+1} &= \left( \operatorname{sgn} h(v\_{N^\*+1, \prime} v\_{N^\*+1}) \right) \left( a(\vartheta\_{N^\*+1}) - h(\vartheta\_{N^\*+1, \prime} \vartheta\_{N^\*+1}) \right) \\ &\quad \left( (2N\_2^2 + N\_2^2 \frac{\ln(2(1-c)N\_2 - 3)}{(1-c)N\_2 - 2}) \right) .\end{aligned}$$

The cubic logarithmic norm of the matrix *D* is equal to

<sup>Λ</sup>2(*D*) = max max <sup>1</sup>≤*k*≤*N*1−<sup>1</sup> *dkk* <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=0 |*h*(*v*¯*k*, *v*¯*l*)| 1 −1 *ζl*(*τ*)*dτ* (*τ*−*v*¯*<sup>k</sup>* )<sup>2</sup> max *<sup>N</sup>*1+2≤*k*≤*N*<sup>∗</sup> *dkk* <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=0 |*h*(*v*¯*k*, *v*¯*l*)| 1 −1 *ζl*(*τ*)*dτ* (*τ*−*v*¯*<sup>k</sup>* )<sup>2</sup> , *<sup>d</sup>*<sup>00</sup> <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=1 |*h*(*v*¯0, *v*¯*l*)| 1 −1 *ζl*(*τ*)*dτ* (*τ*−*v*¯0)<sup>2</sup> , *dN*1*N*<sup>1</sup> <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=0 |*h*(*v*¯*N*<sup>1</sup> , *v*¯*l*)| 1 −1 *ζl*(*τ*)*dτ* (*τ*−*v*¯*N*<sup>1</sup> )<sup>2</sup> , *dN*1+1,*N*1+<sup>1</sup> <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=0 |*h*(*v*¯*N*1+1, *v*¯*l*)| 1 −1 *ζl*(*τ*)*dτ* (*τ*−*v*¯*N*1+1)<sup>2</sup> , *dN*∗+1,*N*∗+<sup>1</sup> <sup>+</sup> *<sup>N</sup>*<sup>∗</sup> ∑ *l*=0 |*h*(*v*¯*N*∗+1, *v*¯*l*)| 1 −1 *ζl*(*τ*)*dτ* (*τ*−*v*¯*N*∗+1)<sup>2</sup> .

,

From (25)–(36) it follows that for sufficiently large *N* Λ2(*D*) < 0 occurs. By Theorem 2 it is clear that the system (43)–(48) (and (37)–(42)) has a unique solution *x*∗ *<sup>N</sup>*(*t*) and *<sup>D</sup>*−1<sup>≤</sup> 1/|Λ2(*D*)|.

Let *x*∗(*t*) and *x*∗ *<sup>N</sup>* be solutions of (16 ) and (37)–(42), respectivety.

We recall the following definitions.

**Definition 3.** *The class <sup>W</sup>r*(*M*, [*a*, *<sup>b</sup>*]), *<sup>r</sup>* <sup>=</sup> 1, 2, ... , *consists of all functions <sup>f</sup>* <sup>∈</sup> *<sup>C</sup>*([*a*, *<sup>b</sup>*]), *which have an absolutely continuous derivative f* (*r*−<sup>1</sup>)(*x*) *and piecewise derivative f* (*r*)(*x*) *with* <sup>|</sup> *<sup>f</sup>* (*r*)(*x*)| ≤ *<sup>M</sup>*.

**Definition 4.** *Denote by <sup>W</sup>r*(*<sup>f</sup>* : *<sup>f</sup>*1, *<sup>f</sup>*2; *<sup>M</sup>*, *<sup>c</sup>*),*<sup>r</sup>* <sup>=</sup> 1, 2, ... , *a set of functions <sup>f</sup>*(*x*), *<sup>x</sup>* <sup>∈</sup> [*a*, *<sup>b</sup>*], *such that <sup>f</sup>*(*x*) = *<sup>f</sup>*1(*x*), *<sup>x</sup>* <sup>∈</sup> [*a*, *<sup>c</sup>*), *<sup>f</sup>*(*x*) = *<sup>f</sup>*2(*x*), *<sup>x</sup>* <sup>∈</sup> (*c*, *<sup>b</sup>*]*, where <sup>f</sup>*1(*x*) <sup>∈</sup> *<sup>W</sup>r*(*M*, [*a*, *<sup>c</sup>*]), *<sup>f</sup>*2(*x*) <sup>∈</sup> *<sup>W</sup>r*(*M*, [*c*, *<sup>b</sup>*]), *<sup>f</sup>*1(*c*) <sup>=</sup> *<sup>f</sup>*2(*c*), *<sup>c</sup>* <sup>∈</sup> (*a*, *<sup>b</sup>*).

Repeating the proof presented in [24] we see that the approximation of *<sup>f</sup>*(*t*) <sup>∈</sup> *<sup>W</sup>*1((*<sup>f</sup>* : *<sup>f</sup>*1, *<sup>f</sup>*2; *<sup>M</sup>*, *<sup>c</sup>*)) by piecewise linear functions constructed on the basis *ζk*(*t*), *k* = 0, 1, ... , *N*<sup>∗</sup> + 1, has the error *C <sup>N</sup>* max(*ω*(*f* (1) <sup>1</sup> , <sup>1</sup> *<sup>N</sup>* ), *ω*(*f* (1) <sup>1</sup> , <sup>1</sup> *<sup>N</sup>* )) for *<sup>f</sup>*(*t*) <sup>∈</sup> *<sup>W</sup>*1((*<sup>f</sup>* : *<sup>f</sup>*1, *<sup>f</sup>*2; *<sup>M</sup>*, *<sup>c</sup>*)), and *<sup>C</sup> <sup>N</sup>*<sup>2</sup> for *<sup>f</sup>*(*t*) <sup>∈</sup> *<sup>W</sup>*2((*<sup>f</sup>* : *f*1, *f*2; *M*, *c*)).

In this paper, we denote the constants that do not depend on *N* by *C*.

Let *<sup>x</sup>*∗(*t*) <sup>∈</sup> *<sup>W</sup>*2((*x*<sup>∗</sup> : *<sup>x</sup>*<sup>∗</sup> <sup>1</sup> , *x*<sup>∗</sup> <sup>2</sup> ; *M*, *c*)), and *x* ∗(1) <sup>1</sup> (*t*) *<sup>C</sup>* ≤ *M*1, *t* ∈ [*a*, *c*], *x* ∗(1) <sup>2</sup> (*t*) *<sup>C</sup>* ≤ *M*2, *t* ∈ [*c*, *b*], *M* = *max*(*M*1, *M*2), 0 < *M* < ∞, where *M* is a bounded constant.

Repeating the arguments given in [24], we arrive at the following statement.

**Theorem 6.** *Let the following conditions be fulfilled:*


*Then the system of Equations (37)–(42) has a unique solution x*∗ *<sup>N</sup>*(*t*) *and the following estimate holds:* ||*x*<sup>∗</sup> − *x*<sup>∗</sup> *<sup>N</sup>*||<sup>1</sup> <sup>≤</sup> *CN*−<sup>1</sup> ln *N.*

#### *3.2. Nonlinear Hypersingular Integral Equations*

Consider the nonlinear hypersingular integral equation:

$$a(t)x(t) + \int\_{-1}^{1} \frac{h(t, \tau, x(\tau))d\tau}{(\tau - t)^2} = f(t) \,. \tag{49}$$

The approximate solution of the Equation (49) we shall seek as a continuous function (17) with the coefficients *γk*. The coefficients *γ<sup>k</sup>* are determined by the following system of nonlinear algebraic equations

$$a(\bar{v}\_k)\gamma\_k + \sum\_{l=0}^{N^\*+1} h(\bar{v}\_k, v\_l, \gamma\_l) \int\_{-1}^1 \frac{\zeta\_l(\tau)}{(\tau - \bar{v}\_k)^2} d\tau = f(\bar{v}\_k), \; k = 0, 1, \ldots, N^\*+1. \tag{50}$$

**Remark 4.** *Note that the set γk*, *k* = 0, 1, ... , *N*<sup>∗</sup> + 1, *is union of sets αk*, *k* = 0, 1, ... , *N*1, *and βk*, *k* = 0, 1, . . . , *N*2.

By computing the hypersingular integrals in (50), we can rewrite the system of Equation (50) as

$$\begin{split} a(\vec{v}\_k)\gamma\_k - h(\vec{v}\_k, \vec{v}\_k, \gamma\_k) 2N\_1^2 \frac{\ln(N\_1 - 1)}{(1 + c)N\_1 - 2} + \sum\_{l = 0}^{N^\* + 1} \prime h(\vec{v}\_k, \vec{v}\_l, \gamma\_l) \int\_{-1}^1 \zeta\_l(\tau) \frac{d\tau}{(\tau - \vec{v}\_k)^2} \\ = f(\vec{v}\_k), \quad k = 1, \dots, N\_1 - 1; \end{split} \tag{51}$$

$$\begin{split} \ln a(\vec{v}\_k)\gamma\_k - h(\vec{v}\_k, \vec{v}\_{k'}, \gamma\_k) 2N\_2^2 \frac{\ln(N\_2 - 1)}{(1 - c)N\_2 - 2} + \sum\_{l = 0}^{N\_\*^\* + 1} {}\_l^{\prime} \gamma\_l h(\vec{v}\_{k'}, \vec{v}\_{l'}, \gamma\_l) \int\_{-1}^1 \mathbb{Z}\_l(\tau) \frac{d\tau}{(\tau - \vartheta\_k)^2} \\ = f(\vec{v}\_k), \quad k = N\_1 + 2, \dots, N^\*; \end{split} \tag{52}$$

$$\begin{split} &a(\vec{v}\_{0})\gamma\_{0} - h(\vec{v}\_{0}, \vec{v}\_{0}, \gamma\_{0})(2N\_{1}^{2} + N\_{1}^{2}\frac{\ln(2(1+c)N\_{1} - 3)}{(1+c)N\_{1} - 2}) \\ &+ \sum\_{l=1}^{N^{\prime}+1} h(\vec{v}\_{k\prime}, \vec{v}\_{l\prime}\gamma\_{l}) \int\_{-1}^{1} \zeta\_{l}(\tau) \frac{d\tau}{(\tau - \vartheta\_{0})^{2}} = f(\vec{v}\_{0}); \end{split} \tag{53}$$

$$\begin{split} &a(\boldsymbol{\upsigma}\_{\mathcal{N}\_{1}})\boldsymbol{\upgamma}\_{\mathcal{N}\_{1}} - h(\boldsymbol{\upsigma}\_{\mathcal{N}\_{1}}, \boldsymbol{\upsigma}\_{\mathcal{N}\_{1}}, \boldsymbol{\upgamma}\_{\mathcal{N}\_{1}}) \big(2\mathcal{N}\_{1}^{2} + \mathcal{N}\_{1}^{2} \frac{\ln(2(1+\varepsilon)N\_{1} - 3)}{(1+\varepsilon)N\_{1} - 2} \big) \\ &+ \sum\_{I=0}^{N^{\*}+1} \prime \prime h(\boldsymbol{\upsigma}\_{\mathcal{N}\_{1}}, \boldsymbol{\upsigma}\_{I}) \int\_{-1}^{1} \zeta\_{I}(\boldsymbol{\upsigma}) \frac{d\boldsymbol{\uppi}}{(\boldsymbol{\uppi} - \boldsymbol{\upsigma}\_{\mathcal{N}\_{1}})^{2}} = f(\boldsymbol{\upsigma}\_{\mathcal{N}\_{1}}); \end{split} \tag{54}$$

$$\begin{split} \mathfrak{a}\left(\bar{\upsilon}\_{N\_{1}+1}\right)\gamma\_{N\_{1}+1} - \mathfrak{h}\left(\bar{\upsilon}\_{N\_{1}+1}, \bar{\upsilon}\_{N\_{1}+1}, \gamma\_{N\_{1}+1}\right) &(2N\_{2}^{2} + N\_{2}^{2} \frac{\ln\left(2(1-\varepsilon)N\_{2}-3\right)}{(1-\varepsilon)N\_{2}-2}) \\ + \sum\_{l=0}^{N^{\*}+1} \mathscr{W}(\bar{\upsilon}\_{N\_{1}+1}, \bar{\upsilon}\_{l}, \gamma\_{l}) \int\_{-1}^{1} \tilde{\zeta}\_{l}(\tau) \frac{d\tau}{(\tau - \bar{\upsilon}\_{N\_{1}+1})^{2}} = f(\bar{\upsilon}\_{N\_{1}+1}); \end{split} \tag{55}$$

$$\begin{split} a(\vec{v}\_{N^\*+1})\gamma\_{N^\*+1} - h(\vec{v}\_{N^\*+1}, \vec{v}\_{N^\*+1}, \gamma\_{N^\*+1})(2N\_2^2 + N\_2^2 \frac{\ln(2(1-c)N\_2 - 3)}{(1-c)N\_2 - 2}) \\ + \sum\_{l=0}^{N^\*} \, ^{\prime\prime\prime}h(\vec{v}\_{N^\*+1}, \vec{v}\_{l}, \gamma\_{l}) \int\_{-1}^{\prime} \zeta\_{l}(\tau) \frac{d\tau}{(\tau - \vartheta\_{N^\*+1})^2} = f(\vec{v}\_{N^\*+1}). \end{split} \tag{56}$$

Here ∑ , ∑ , ∑ indicate summations over *l* = *k*, *l* = *N*1, *l* = *N*<sup>1</sup> + 1, respectively. The Frechet derivative on a vector (*α*0, *α*1, ··· , *αN*∗+1) in the space *RN*∗+<sup>1</sup> is equal to

 *a*(*v*¯*k*)*γ<sup>k</sup>* − *h* <sup>3</sup>(*v*¯*k*, *<sup>v</sup>*¯*k*, *<sup>γ</sup>*¯ *<sup>k</sup>*)2*N*<sup>2</sup> 1 ln(*N*1−1) (1+*c*)*N*1−2*γ<sup>k</sup>* <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=0 *h* <sup>3</sup>(*v*¯*k*, *v*¯*l*, *γ*¯*l*)*γ<sup>l</sup>* 1 −1 *ζl*(*τ*) *<sup>d</sup><sup>τ</sup>* (*τ*−*v*¯*<sup>k</sup>* )<sup>2</sup> , *<sup>k</sup>* <sup>=</sup> 1, . . . , *<sup>N</sup>*<sup>1</sup> <sup>−</sup> 1; *a*(*v*¯*k*)*γ<sup>k</sup>* − *h* <sup>3</sup>(*v*¯*k*, *<sup>v</sup>*¯*k*, *<sup>γ</sup>*¯ *<sup>k</sup>*)2*N*<sup>2</sup> 2 ln(*N*2−1) (1−*c*)*N*2−2*γ<sup>k</sup>* <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=0 *h* <sup>3</sup>(*v*¯*k*, *v*¯*l*, *γ*¯*l*)*γ<sup>l</sup>* 1 −1 *ζl*(*τ*) *<sup>d</sup><sup>τ</sup>* (*τ*−*v*¯*<sup>k</sup>* )<sup>2</sup> , *<sup>k</sup>* <sup>=</sup> *<sup>N</sup>*<sup>1</sup> <sup>+</sup> 2, . . . , *<sup>N</sup>*∗; *a*(*v*¯0)*γ*<sup>0</sup> − *h* <sup>3</sup>(*v*¯0, *<sup>v</sup>*¯0, *<sup>γ</sup>*¯0)*γ*0(2*N*<sup>2</sup> <sup>1</sup> + *<sup>N</sup>*<sup>2</sup> 1 ln(2(1+*c*)*N*1−3) (1+*c*)*N*1−<sup>2</sup> ) <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=1 *h* <sup>3</sup>(*v*¯0, *v*¯*l*, *γ*¯*l*)*γ<sup>l</sup>* 1 −1 *ζl*(*τ*) *<sup>d</sup><sup>τ</sup>* (*τ*−*v*¯0)<sup>2</sup> ; *a*(*v*¯*N*<sup>1</sup> )*γN*<sup>1</sup> − *h* <sup>3</sup>(*v*¯*N*<sup>1</sup> , *<sup>v</sup>*¯*N*<sup>1</sup> , *<sup>γ</sup>*¯ *<sup>N</sup>*<sup>1</sup> )*γN*<sup>1</sup> (2*N*<sup>2</sup> <sup>1</sup> + *<sup>N</sup>*<sup>2</sup> 1 ln(2(1+*c*)*N*1−3) (1+*c*)*N*1−<sup>2</sup> ) <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=0 *h* <sup>3</sup>(*v*¯*N*<sup>1</sup> , *v*¯*l*, *γ*¯*l*)*γ<sup>l</sup>* 1 −1 *ζl*(*τ*) *<sup>d</sup><sup>τ</sup>* (*τ*−*v*¯*N*<sup>1</sup> )<sup>2</sup> ; *a*(*v*¯*N*1+1)*γN*1+<sup>1</sup> − *h* <sup>3</sup>(*v*¯*N*1+1, *<sup>v</sup>*¯*N*1+1, *<sup>γ</sup>*¯ *<sup>N</sup>*1+1)*γN*1+1(2*N*<sup>2</sup> <sup>2</sup> + *<sup>N</sup>*<sup>2</sup> 2 ln(2(1−*c*)*N*2−3) (1−*c*)*N*2−<sup>2</sup> ) <sup>+</sup> *<sup>N</sup>*∗+<sup>1</sup> ∑ *l*=0 *h* <sup>3</sup>(*v*¯*N*1+1, *v*¯*l*, *γ*¯*l*)*γ<sup>l</sup>* 1 −1 *ζl*(*τ*) *<sup>d</sup><sup>τ</sup>* (*τ*−*v*¯*N*1+1)<sup>2</sup> ; *a*(*v*¯*N*∗+1)*γN*∗+<sup>1</sup> − *h* <sup>3</sup>(*v*¯*N*∗+1, *<sup>v</sup>*¯*N*∗+1, *<sup>γ</sup>*¯ *<sup>N</sup>*∗+1)*γN*∗+1(2*N*<sup>2</sup> <sup>2</sup> + *<sup>N</sup>*<sup>2</sup> 2 ln(2(1−*c*)*N*2−3) (1−*c*)*N*2−<sup>2</sup> ) <sup>+</sup> *<sup>N</sup>*<sup>∗</sup> ∑ *l*=0 *h* <sup>3</sup>(*v*¯*N*∗+1, *v*¯*l*, *γ*¯*l*)*γ<sup>l</sup>* 1 −1 *ζl*(*τ*) *<sup>d</sup><sup>τ</sup>* (*τ*−*v*¯*N*∗+1)<sup>2</sup> . (57)

Here ∑ , ∑ ∑ indicate summations over *l* = *k*, *l* = *N*1, *l* = *N*<sup>1</sup> + 1, respectively. The notation *h* <sup>3</sup>(*t*, *<sup>τ</sup>*, *<sup>u</sup>*) = *<sup>δ</sup>h*(*t*,*τ*,*u*) *<sup>δ</sup><sup>u</sup>* is used here.

Let the Equation (49) has the unique solution *x*∗(*t*) inside the ball *B*(*x*∗, *δ*). We shall assume that the Frechet derivative (57) in the ball *RN*∗+1(*x*∗, *δ*) satisfies the conditions of Theorem 5. Thus, according to statements of the Theorem 5, the solution of the system of differential equations

$$\begin{array}{rcl} \frac{da\_l(\sigma)}{d\sigma} & = & a(\overline{\mathfrak{f}}\_l)a\_l(\sigma) - \sum\_{k=0}^{N^\*+1} h(\overline{\mathfrak{f}}\_l, \overline{\mathfrak{f}}\_k, a\_k(\sigma)) N \left( \frac{1}{2l - 2k + 1} - \frac{1}{2l - 2k - 1} \right) \\\ -f(\overline{\mathfrak{f}}\_l), \; l & = & 0, 1, \cdots, N^\* + 1, \end{array} \tag{58}$$

converges to the solution of the Equation (49).

Thus, we have proven the following statement.

**Theorem 7.** *Let the following conditions hold:*

*(1) Equation (49) has a unique solution x*∗(*t*) *inside some ball B*(*x*∗, *<sup>δ</sup>*), *<sup>x</sup>*<sup>∗</sup> <sup>∈</sup> *<sup>W</sup>*2(*x*<sup>∗</sup> : *<sup>x</sup>*<sup>∗</sup> <sup>1</sup>, *x*<sup>∗</sup> <sup>2</sup>; *M*, *c*); *(2) The Frechet derivative (57) in the ball RN*∗+1(*x*∗, *δ*) *satisfies the conditions of Theorem 5.*

*Then the system of Equations (51)–(56) has a unique solution inside the ball B*(*x*∗, *δ*)*, and the solution of Equation (58) converges to this solution.*

The effectiveness of the presented algorithms is illustrated by solving two hypersingular integral equations modeling aerodynamics problems.

**Example 1.** *Let us illustrate the effectiveness of continuous method by solving the following linear hypersingular equation*

$$\int\_{-1}^{1} \frac{\mathbf{x}(\tau)}{\left(\tau - t\right)^{2}} d\tau = f(\gamma\_{1}, \gamma\_{2}, t), \tag{59}$$

*where f*(*γ*1, *γ*2, *t*) *is the given right-hand side of the equation:*

$$\begin{array}{rcl} f(\gamma\_1, \gamma\_2; t) &=& \gamma\_1 - \gamma\_2 + (a\_1 - a\_2)\frac{1}{t} - (a\_1 + \gamma\_1 t)\frac{1}{1 + t} \\ &- (a\_2 + \gamma\_2 t)\frac{1}{1 - t} + \gamma\_1 \ln|\frac{t}{1 + t}| + \gamma\_2 \ln|\frac{1 - t}{t}|. \end{array}$$

*The exact solution of the equation is x*(*t*)=(*x*1(*t*), *x*2(*t*)); *xi*(*t*) = *ai* + *γit*, *i* = 1, 2*.*

To solve the Equation (59) numerically we use the continuous method for solving operator equations and arrive to the following evolution equation

$$\frac{da\_k(\sigma)}{d\sigma} = \sum\_{l=0}^{N^\*+1} Na\_l(\sigma) (\frac{1}{2k+2l-1} - \frac{1}{2k+2l+1}) - f(\gamma\_1, \gamma\_2; \overline{v}\_k), k = 0, 1, \dots, N^\*+1.$$

Nodes *vk*, *v*¯*k*, *k* = 0, 1, . . . , *N*<sup>∗</sup> + 1, have been entered above.

In Figure 1 we show the trajectories of the exact solution of the Equation (59); its approximate solution, received with continuous method; and values of error.

**Figure 1.** Numerical solutions for the linear hypersingular equation with a discontinuous right-hand side example.

Here *a*<sup>1</sup> = 1, *a*<sup>2</sup> = 1.5, *γ*<sup>1</sup> = 0.5, *γ*<sup>2</sup> = 0.3.

**Example 2.** *Let us illustrate the effectiveness of the continuous method for the solutions of nonlinear hypersingular equations*

$$\int\_{-1}^{1} \frac{\mathbf{x}^2(\tau)}{(\tau - t)^2} d\tau = f(\gamma\_1, \gamma\_{2'}, t) \tag{60}$$

*where f*(*γ*1, *γ*2, *t*) *is the given right-hand side of the equation:*

$$\begin{array}{rcl} f(\gamma\_1, \gamma\_2; t) &=& \gamma\_1^2 + 2\gamma\_1 a\_1 + \gamma\_1^2 t + \frac{a\_1^2}{l} + (2\gamma\_1 a\_1 + 2\gamma\_1^2 t) \ln|\frac{t}{1+t}| \\ &- (a\_1^2 + 2\gamma\_1 a\_1 t + \gamma\_1^2 t^2) \frac{1}{1+t} + \gamma\_2^2 - 2a\_2 \gamma\_2 - \gamma\_2^2 t - \frac{a\_2^2}{l} - (a\_2^2 + 2a\_2 \gamma\_2 t) \\ &+ \gamma\_2^2 t^2) \frac{1}{1-t} + (2a\_2 \gamma\_2 + 2\gamma\_2^2 t) \ln|\frac{1-t}{t}|.\\ \end{array}$$

The exact solution of the equation is *x*(*t*)=(*x*1(*t*), *x*2(*t*)); *xi*(*t*) = *ai* + *γit*, *i* = 1, 2.

It easy to see that, if *x*(*t*) is a solution of the Equation (60), then functions −*x*(*t*), |*x*(*t*)| and −|*x*(*t*)| are solutions of this equation too.

To solve the Equation (60) numerically we use the continuous method and receive the following evolution equation

$$\frac{d\alpha\_k(\sigma)}{d\sigma} = \sum\_{l=0}^{N^\*+1} N a\_l^2(\sigma) (\frac{1}{2k+2l-1} - \frac{1}{2k+2l+1}) - f(\gamma\_{1\prime}\gamma\_{2\prime}\vec{v}\_l)\_{\prime\prime}$$

*k* = 0, 1, . . . , *N*∗ + 1.

At first, we take *αk*(0) = 0.0 as an initial condition in order to demonstrate applicability of our method in cases of the Newton–Kantorovich method, the minimal residual method and other numerical methods; using in their construction the derivative of nonlinear operator is not applicable. Indeed, in this case the Frechet derivative (57) is not only degenerate—and, therefore, not invertable—but is an identical zero.

In Figure 2 we put *a*<sup>1</sup> = 1, *a*<sup>2</sup> = 1.4, *γ*<sup>1</sup> = 0.5, *γ*<sup>2</sup> = −0.4.

**Figure 2.** Numerical solution for the nonlinear hypersingular equation with a discontinuous right-hand side example.

In Figure 2 we show the trajectories of the exact solution of the Equation (60), its approximate solution, received with continuous method and values of error.

The exact solution at *t* = 0 has a jump discontinuity of *h* = 0.4. The slopes of the exact solution also change at *t* = 0. In Figure 2 we demonstrate that the numerical solution approximates the exact one at [−1, 0) well. At *<sup>t</sup>* <sup>=</sup> 0 the approximate solution has a jump ˜ *h* = 0.15.

#### **4. Summary and Discussion**

An iterative projection method for solving linear and nonlinear hypersingular integral equations has been proposed. The method is based on the use of sufficient conditions for asymptotic stability of ODE systems. Stability conditions are expressed in terms of the logarithmic norms of the corresponding matrices. In a number of spaces often used in computational mathematics, the calculation of logarithmic norms does not cause difficulties, even for large-dimensional matrices.

What are the advantages of the presented method?


The detailed bibliography of approximation methods of hypersingular integral equations of the first and the second kinds is given in [32]. The bibliography on solving hypersingular integral equations of the first kind is presented in [45].

Mostly, papers devoted to hypersingular integral equations of the first kind focused to seek solutions in the class of functions <sup>√</sup> 1 − *t*2*φ*(*t*), where *φ*(*t*) is a smooth function. The presented method provides solutions in a general form.

The theoretical justification of the method is based on Lyapunov stability theory. It connects convergence of the method to the sign of the approximate system matrix logarithmic norm.

Said justification has advantages that allow us


The major advantage of the method for nonlinear equations is as follows.

The Newton–Kantorovich method requires the Frechet derivative reversibility at each iteration step. Similar conditions are required when using other iteration methods. Our method lacks such a deficiency. It does not put any restrictions on the Frechet derivative of the nonlinear operator.

**Author Contributions:** Conceptualization, I.B.; Data curation, V.R.; Formal analysis, I.B. and V.R.; Funding acquisition, V.R.; Investigation, I.B. and V.R.; Methodology, V.R.; Project administration, I.B.; Resources, V.R.; Software, V.R. and A.B.; Supervision, I.B.; Visualization, A.B.; Writing—original draft, I.B.; Writing—review and editing, V.R. and A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
