Method of mechanical quadratures for solving singular integral equations of various types

The method of mechanical quadratures is proposed as a common approach intended for solving the integral equations defined on finite intervals and containing Cauchy-type singular integrals. This method can be used to solve singular integral equations of the first and second kind, equations with generalized kernel, weakly singular equations, and integro-differential equations. The quadrature rules for several different integrals represented through the same coefficients are presented. This allows one to reduce the integral equations containing integrals of different types to a system of linear algebraic equations.


Introduction
The rapid development of computer technologies in the second half of the last century revived the interest of researchers in computational methods and their effective use, in particular, in methods for solving singular integral equations with Cauchy kernel. The importance of "bringing to a number" the solutions of singular integral equations had been repeatedly emphasized by well-known mathematicians, such as N. I. Mukhelishvili, M. A. Lavrentiev, S. G. Mikhlin, H. Multhopp, I. N. Vekua.
The first important results in this direction were obtained by M. A. Lavrentiev [1]. Noting the importance of this work, Academician N. I. Muskhelishvili noted in all editions of his monograph [2]: "The further development of this and similar methods for solving singular integral equations approximately is, I believe, one of the most important problems in the theory of these equations." From the middle of the twentieth century, approximate methods for solving singular integral equations began to develop intensively. There are a huge number of papers pertaining to this problem, but the vast majority of them are related to the most widespread particular case, when the behavior of the singular integral density at the ends of the segment of integration is described by a root function. The number of papers devoted to approximate calculations of singular integrals whose densities contain the weight function of Jacobi polynomials with arbitrary admissible real exponents is much smaller. The works, where these exponents can be complex, are not known to the authors.
Of particular interest is the work of F. Erdogan, G. D. Gupta and T. S. Cook [3]. This work is important because of a wider range of the considered material and good presentation of the methods for investigating the behavior of solutions at singular points and for solving singular integral equations of various types. However, unfortunately, in this work, there is no 2 1234567890 ''"" single approach to the problem of reducing singular integral equations to a system of algebraic equations. Despite such a large number of works devoted to approximate methods for solving singular integral equations, interest in this problem does not decrease and continues to be on the agenda. This is testified by a considerable number of publications in the recent years, for example, [4][5][6][7][8][9][10][11].
In this paper, we give a series of quadrature formulas for various integrals arising in singular integral or integro-differential equations of the first and second kind, and describe the method of mechanical quadratures for solving these equations.

Quadrature formulas
In this section, several quadrature formulas of Gauss type for calculating integrals containing the weight function of Jacobi polynomials are given. Some of them were obtained in [12][13][14][15][16].
Suppose we have a function Φ(x) defined on the interval (−1, 1) and representable in the form where ϕ(x) is a function that satisfies the Hölder condition on the interval [−1, 1]. The application of these quadrature formulas for complex values of exponents α and β is explained in detail in [17].
Replacing the function ϕ(x) by an interpolation polynomial we obtain the quadrature formulas Here the points ξ i (i = 1, n) are the roots of the Jacobi polynomial P The weighting factors w i (i = 1, n) are defined as Note that the functions R (α,β) n (z) and L (α,β) n (z) are uniquely determined in the complex plane which is cut respectively along the segment [−1, 1] and the ray [−1, ∞], when for the cut points, the half-sum of the values from above (y + i0) and below (y − i0) is used.
It is known that quadrature formula (2) is the formula of the highest algebraic accuracy, i.e., it is satisfied exactly if the function ϕ(x) is a polynomial up to order 2n − 1. In the general case, i.e., at an arbitrary value z, formula (3) is not a formula of the highest algebraic accuracy, but it becomes so if z coincides with the roots of the function R (α,β) n (z). The remaining formulas are exact only in the case when the function ϕ(x) is a polynomial of order less than n.
The main advantage of quadrature formulas (2)-(7) is that they are all expressed in terms of the same quantities ϕ(ξ j ), which are the unknowns in the integral equations. This allows one to reduce an integral equation containing arbitrary combinations of the integrals shown above to a system of linear algebraic equations with respect to the unknowns ϕ(ξ j ).

Method of mechanical quadratures
Below we step by step describe the basic procedures required to apply the method of mechanical quadratures for solving the singular integral equations.  (3), also contains integrals of some type from (2)- (7). We note that, in the integral equation, formula (2) appears in the form where R(x, y) is a regular function in both variables. It is obvious that, in the case of an arbitrary interval of integration, the equation can be reduced to the interval (−1, 1) by a simple linear transformation. Further, the solution of the equation is sought in the form (1) with unknown coefficients α and β.
Step 1. We assume that the solution of the given equation has representation (1), where the unknown function ϕ(x) is sufficiently smooth and bounded on the closed interval [−1, 1]. We will call it the regular part of the solution of the given equation.
Using the results of N. I. Mukhelishvili on the behavior of the Cauchy-type integral at the ends of the line of integration [2], we write the principal terms of the asymptotics of the behavior of each of the integrals in the equation, and of the right-hand side, the known part of the equation in a neighborhood of each of the ends of the interval of integration. Equating the coefficients of the corresponding terms of the asymptotics to each other, we obtain equations for determining the exponents α and β. Usually, these roots are in the interval (−1, 1). In the case when the equations obtained have two roots in (−1, 1), usually of opposite signs, the choice of the necessary root is dictated by the statement of the problem.
Step 2. Having found the exponents α and β and choosing a certain order of approximation n, we can replace each of the integrals of the equation, according to formulas (2)-(7), by a quadrature sum. As a result, we obtain a functional equation that is linear with respect to unknown coefficients, which are the values of the regular part of the solution of the integral equation in the roots of the Jacobi polynomial corresponding to the found exponents.
Step 3. To determine the unknown coefficients, it is necessary to equate the right-and lefthand sides of the obtained functional equation at a certain number of points, called collocation points. The effective choice of these points is dictated by the following considerations.
Suppose we have a singular integral equation of the first kind, which contains only the singular integral (3) and the regular integral (9). As was mentioned above, because, at the roots of the function R (α,β) n (z), the quadrature formulas for both integrals are quadrature formulas of the highest algebraic accuracy, it is obvious that these roots should be chosen as collocation points. This also fixes the number of collocation points necessary to obtain a closed system of linear algebraic equations. In the case under consideration, when the regular kernel and the righthand side of the equation do not have singularities at the ends of the interval of integration, the exponents α and β can take the values α = ±β = ±0.5. As was indicated above, the choice of adequate values is dictated by the formulation of the problem. For clarity, we illustrate this with an example of the classical problem of indenting a rigid stamp into an elastic half-plane.
Let us have an elastic strip or an inhomogeneous half-plane into which a rigid stamp of a certain shape is pressed under the action of a force P (figure 1). It is assumed that there is a smooth contact and only normal stresses appear in the contact zone.
After reducing the contact zone intervals (−a, a), (−c, a), and (−c, b) to the interval (−1, 1) in all three cases, the problem reduces to solving the equation under the condition of the stamp equilibrium Here Φ(ξ) = p(x(ξ))/µ is the dimensionless contact pressure, f (η) is a known dimensionless function determined by the shape of the base of the stamp, µ is the modulus of shear of the material of the base, and P * = P/(aµ) is the dimensionless value of the pressing force. The solution of equation (10) is sought in the form (1). According to step 1, we find that the exponents α and β can take the values α = ±β = ±0.5.
In case a) in figure 1, at the ends of the contact zone, the contact pressure has a singularity and, therefore, the values α = β = −0.5 should be selected. The function R (−0.5,−0.5) n (y) becomes the Chebyshev polynomial of the second kind U n−1 (y) which has n − 1 roots, and therefore the number of collocation points is n − 1. This is in complete agreement with the well-known fact that equation (10) does not have a unique solution and is determined up to the summand A(1 − ξ 2 ) −0.5 whose coefficient A must be found from the condition (11).
In case b) in figure 1, at the left end of the contact zone, due to continuity, the contact pressure must go to zero and, therefore, the values α = −0.5, β = 0.5 should be selected. The function R (−0.5,0.5) n (y) is the Jacobi polynomial P (0.5,−0.5) n (y) and the number of collocation points is n. Therefore, the resulting system of linear algebraic equations will contain n + 1 equations. However, it should be noted that in addition to unknown coefficients ϕ(ξ j ), in this case the location of the left end of the contact zone c/a is also unknown. Since P * appears only in condition (11), it is more expedient to solve the inverse problem, i.e., to set the value c/a, find the coefficients ϕ(ξ j ) from the system of collocation equations, and, further, from equilibrium condition (11), to determine for what P * this is the case. The system of linear algebraic equations in this case becomes n j=1 w j ϕ(ξ j ) ξ j − t k + n j=1 w j K(ξ j , t k )ϕ(ξ j ) = f (t k ) (k = 1, 2, . . . , n), In case c) in figure 1, the contact pressure is zero at both ends of the contact zone. The function R (0.5,0.5) n (y) is a Chebyshev polynomial of the first kind T n+1 (y), which has n + 1 roots, that is, the number of collocation points is n + 1. It is known that, in this class of functions, equation (10) has a solution only if the following condition is satisfied: Since, unlike the preceding case, the location of the right end of the contact zone b/a is also unknown here, condition (14) is necessary for its determination.
In [18] it is shown that if an additional unknown γ 0n , with a unit coefficient, is introduced into the resulting system of linear algebraic equations, then lim n→∞ γ 0n = 0 if and only if condition (14) is satisfied. The unknown γ 0n is called the regularizing factor, and its behavior serves as a convenient indicator of the existence of the solution. As in the preceding case, it is convenient here to solve the inverse problem, i.e., to set, for example, the value c/a, then select a value b/a so that when the collocation system is solved, then the regularizing factor γ 0n is close to zero. Then the corresponding value P * is determined from condition (11).
It is worth to note that, in this case, the unknown γ 0n can be interpreted as the angle of rotation of the stamp. Really, if both values c/a and b/a are set, and the value P * is found, then γ 0n means the angle of the stamp rotation needed for the contact zone to coincide with the interval (−c/a, b/a). When γ 0n = 0, the stamp is pressed vertically into the half-plane.
Unlike the equations of the first kind, other singular integral equations contain at least one term for which the approximation formula is exact only for polynomials of order at most n − 1. Although the above justification for choosing the roots of R (α,β) n (z) as collocation points is not valid here, nevertheless, the experience shows that it is also acceptable in this case. Indeed, the number of roots of the function R (α,β) n (z), if any, coincides with the necessary number of collocation points. The choice of collocation points themselves can be arbitrary, since there are no restrictions on the external variable in formulas (2)- (7).
When solving the singular integral equations of the second kind, it is more convenient to use a slightly different representation of quadrature formula (3).
The function R (α,β) n (z) for the points of the interval of integration can be written as Since, in the case of a singular integral equation of the second kind, the index of the equation κ is equal to κ = α + β = {−1, 0, 1}, the last expression becomes  (3), after simple transformations, we obtain the quadrature formula Suppose we have a singular integral equation of the second kind Under the assumption that the regular kernel and the right-hand side of the equation have no singularity at the ends of the integration interval, the exponents α and β are determined from the equations cot(απ) − λ = 0, cot(βπ) + λ = 0.
Here, depending on the value of the index κ, all three cases mentioned above can occur. We note that, using the quadrature formula (17) and choosing the zeros of the polynomial P (−α,−β) n+κ (y) as collocation points, equation (18) reduces to a system of linear algebraic equations of the form (12), (13) or (15).
A similar approach can be used for more complex singular integral equations as well.

Approbation of the method
The method of mechanical quadratures in application to singular integral equations of the first kind is well known, and there are many works using it. Below we cite the works where the method of mechanical quadratures is used to solve singular equations of other types. Note that in many of these papers, the method is called the method of discrete singularities. In [19,20], contact problems of theory of elasticity are solved in the case where the defining equations are singular integral equations of the first kind with a right-hand side containing a logarithmic singularity at the ends of the interval of integration.
The contact problems for moving stamps, which reduce to solving singular integral equations of the second kind with real coefficient of the free term and, consequently, real exponents of the Jacobi weight function, were solved in [21][22][23][24][25][26]. The problem of indenting a stamp into an elastic massive body with allowance for the bonding zone in the middle part of the contact zone and slip zones at its ends [27][28][29], as well as the mixed problems for a piecewise homogeneous plane with interphase cracks or rigid inclusions, considered in the formulation of the antiplane problem of the theory of elasticity [30], can also be reduced to a system of such equations.
A large number of papers are devoted to problems which can be reduced to systems of singular integral equations of the second kind with imaginary coefficient of the free term and, consequently, with complex exponents of the Jacobi weight function. In [31][32][33][34][35], problems for homogeneous massive bodies containing thin rigid inclusions are solved, when, on one side, the inclusion is connected to a massive body, and on the other side, a crack is formed. In [36][37][38], problems for piecewise-homogeneous massive bodies with cracks and thin rigid inclusions on the material separation line are solved.
In [39][40][41][42][43][44][45][46][47][48], problems of the theory of elasticity in plane and antiplane formulations for massive bodies containing stress concentrators of the types indicated above are solved, which reduce to solving systems of singular integral equations with a generalized Cauchy kernel. These equations are also known as equations with a fixed singularity. This type of singularities often appears when an end of a stress concentrator reaches the boundary of the homogeneous body or the line of separation of materials of the piecewise homogeneous body or another concentrator.
Contact problems for massive bodies, reinforced by thin elastic plates or containing thin elastic inclusions, reduce to solving singular integro-differential equations. In [49][50][51][52][53][54][55], such problems were solved by the method of mechanical quadratures. In [53] the solution of the Prandtl equation is presented in two ways using the method of mechanical quadratures. In a number of works, the rate of convergence of the approximate solution to the exact one, depending on the order of approximation, is investigated by numerical analysis. The practice shows that, with the correct choice of the exponents of the singularity α and β and the corresponding quadrature formulas, at n = 10 ∼ 12, it is already possible to obtain solutions with accuracy higher than four significant digits.

Conclusions
The paper presents the development of the method of mechanical quadratures, which is well known for solving singular integral equations of the first kind. The quadrature formulas given in the paper, which are written with respect to the same coefficients (values of the unknown functions at the same points), allow us considerably to expand the scope of the method. This is tested by a large number of problems solved by this method whose defining equations qualitatively differ from each other.