Product Quasi-Interpolation in Logarithmically Singular Integral Equations

Abstract A discrete high order method is constructed and justified for a class of Fredholm integral equations of the second kind with kernels that may have boundary and logarithmic diagonal singularities. The method is based on the improving the boundary behaviour of the kernel with the help of a change of variables, and on the product integration using quasi-interpolation by smooth splines of order m. Properties of different proposed calculation schemes are compared through numerical experiments using, in particular, variable precision interval arithmetics.


Introduction
In the present paper we treat a fully discrete method of accuracy O(h m ) (see Section 4) for the integral equation with the logarithmic diagonal singularity in the kernel. The coefficient functions a, b ∈ C m ([0, 1] × (0, 1)) and the free term f ∈ C[0, 1] ∩ C m (0, 1) may have certain boundary singularities described below in detail. Due to diagonal and boundary singularities of the kernel, the derivatives of the solution to equation (1.1), as a rule, have certain boundary singularities. In Section 2 we reduce (1.1) with the help of a smoothing change of variables to a similar problem v(t) = 1 0 A(t, s) log |t − s| + B(t, s) v(s) ds + g(t), 0 ≤ t ≤ 1, in which the coefficients A(t, s) and B(t, s) have no singularities and vanish for s = 0 and s = 1. Simultaneously the singularities of the solution will be milder or disappear at all for suitable parameters of the change of variables. On the contrast, the logarithmic diagonal singularity of the kernel of equation (1.1) still remains to be present. In Section 3 we recall some results about quasiinterpolation of functions by smooth splines of order m (of degree m − 1) on the uniform grid of the step size h = 1/n. In Section 4 we introduce and justify our method based on quasi-interpolation of the products A(t, s)v(s) and B(t, s)v(s) by smooth splines. The use of smooth splines of order m enables an m-fold reduction of degrees of freedom, compared with the traditional use of discontinuous splines in the product integration methods [1,3,10]. The use of quasi-interpolation instead of real interpolation leads to some simplification of the numerical algorithm. Section 5 is devoted to the computation of the coefficients in the matrix form of the method. An important advantage of our method is that we can present simple exact formulae for the coefficients occurring in the discretized problem. Unfortunately, for big n, far from the main diagonal, these formulae become numerically unstable, so we complete them by traditional numerically stable quadrature approximation; the error of quadrature approximation rapidly decreases as we move away from the main diagonal. The use of only exact formulae in double precision arithmetics for all coefficients is usually quite acceptable in engineer computations. A numerical check and illustration of the method is performed in Section 6; also the interval arithmetics approach is used for the analysis of the numerical accuracy of exact formulae and the quadrature computation of the coefficients. Linear integral equations with a logarithmic diagonal singularity often occur modelling physical processes. For instance, one of the main equations in radiative transfer theory, the Milne integral equation [2], has the form is the integral exponent function, c = lim n→∞ n k=1 1 k − log n ≈ 0.5772 is the Euler constant.
The present paper is a continuation of our work [13], where the convergence, error and numerical analysis of the same method has been presented for the integral equation where ν ∈ (0, 1). Since the convergence analysis for equation (1.1) differs from that for (1.2) presented in [13] only in small details, we omit detailed proofs and concentrate mainly on the specific numerical aspects of the method for equation (1.1). Denote by T the integral operator of equation (1.1), The following lemma can be proved by standard argument, cf. [3,5].

The father B-spline
The father B-spline B m of order m (or of degree m − 1) can be defined by the formula Here, as usual, 0! = 1, 0 0 := lim x↓0 x x = 1, Let us recall some properties of B m :

Spline interpolation
Introduce in R the uniform grid hZ = {ih : i ∈ Z} of the step size h > 0. Denote by S h,m , m ∈ N, the space of (maximally smooth) splines of order m (of degree m − 1) and defect 1 with the knot set hZ. It consists of all functions It occurs (see [11]) that also for m ≥ 3 conditions (3.4) uniquely determine d j , j ∈ Z, in the space of bounded bisequences, namely, Lemma 3 [see [4,15]

Spline quasi-interpolation
It occurs [6] that for v with uniformly continuous mth derivative v (m) and 2p > m, the quasi-interpolant Q h,m v is asymptotically of the same accuracy as the interpolant Q h,m v. It is reasonable to take the smallest p ∈ N for which 2p > m; denote it by m 1 , Lemma 4 [see [6]].
where the constant c is independent of f , h and i (more about c see in [6]).

Product Quasi-Interpolation Method
Let h = 1/n, n ∈ N, n ≥ m. We approximate equation where the quasi-interpolation operator Q h,m is applied to products A(t, s)v(s) and B(t, s)v(s) as functions of s. This can be done for v given on [0, 1] since A(t, s) = B(t, s) = 0 for s ≤ 0 and for s ≥ 1; recall that under con- and Q h,m (B(t, s)v n (s)) depend on the function v n through its knot values v n ((k + m 2 )h), we obtain a closed system to determine v n ((k + m 2 )h) by collocating equation (4.1) at the points (i + m 2 )h. In this way, the following system of equations can be derived (cf. [13]): recall also formulae (3.5)-(3.7) for α k,m , m 0 and m 1 . The dimension of system (4.2) is n − 1 for even and n for odd n. Having the quantities α k,m , [13] for a numerically stable evaluation of the function Φ(t, s) occurring in the definition of B(t, s) and for some other computational details, in particular, for the computation of u n (x) := v n (ϕ −1 (x)), x ∈ [0, 1], the approximate solution to (1.1). The evaluation of integrals (4.6) defining β 0 j is elementary, in particular, β 0 j = h for j = 0, . . . , n − m, due to (3.3). It remains to explain how to evaluate the integrals (4.7) defining β i,j . This is the most delicate part of the method. In Section 5 we present two exact and one approximate algorithm to compute β i,j ; each of algorithms has its advantages and disadvantages. Exploiting the symmetries, the number of different integrals reduces to O(n), and they can be computed in O(n) flops. System (4.2) can be solved in O(n 2 log n) flops by GMRES (n 2 flops per a matrix to vector multiplication, O(log n) iterations, see [14]). Thus the complexity of the method is O(n 2 log n) flops.
Having found the solution v i,n , i = −m 0 , . . . , n − m 1 , of system (4.2) we can use (4.1) to define the Nystrom extension v n (t) of the grid solution for any t ∈ [0, 1]; this is expensive. An essentially cheaper way is to construct the where v ∈ C[0, 1] is the unique solution of equation ( where v n is the quasi-interpolant constructed onto the solution of system While having in mind the implementation of the proposed methods using finite precision arithmetics, we observe that from the perspective on numerical stability the main difficulty is the computation of quadrature coefficients β i,j defined in (4.7). Therefore, we propose and analyse three different calculation schemes for β i,j in the next section with the numerical verification of the proposed algorithms in Section 6.

Computation of β i,j
The condition number of system (4.2) is uniformly bounded in n, and small perturbations ∆β ij of β ij cause the error |∆v h | ≤ c max i,j |∆β ij | in the spline collocation solution v n of (4.1), where the constant c is independent of n.
We elaborate on three different computation schemes for β i,j ; the first two of them are based on exact arithmetics while the third one uses an approximate method for calculation of underlying integrals. We will see that the exact schemes compared with the inexact one behave differently for distinct values of i and j, namely depending on the value of modulus |i − j|.

Exact algorithm 1
The integrals of type (4.7) can be evaluated exploiting the following result.
Proof. This claim can be easily checked by induction in m.
Applying Lemmas 5 and 6 we obtain the (exact) formula where γ j = γ j+1 − γ j is the forward difference and (the difference m γ i,j is taken with respect to the argument j). Exact formulae for β i,j are delicate. Formula (5.1)-(5.2) is suitable for moderate values of |i − j + m 2 |. On the other hand, for great |i − j|, the quantities γ i,j are large in modulus (comparable with |i − j| m ) in contrast to their differences ∆ m γ i (which are comparable with m!), causing a loss of accuracy in standard floating-point arithmetics such as IEEE single-or doubleprecision floating-point format. Nevertheless, for |i − j| m /m! ≤ 10 8 , using double-precision arithmetics, we obtain seven correct decimal digits of γ i,j by (5.1)-(5.2), which is usually sufficient in engineer calculations. For instance, in case m = 4 (cubic splines), seven correct digits of γ i,j are obtained for |i − j| ≤ 220; according to error estimates (4.8) and (4.9), m = 4, n ≈ 200 usually yields a sufficient accuracy of the approximate solution.
If more accurate approximation is needed, formulae (5.1)-(5.2) can be applied for all |i − j| only if a sufficiently high precision arithmetics is used; in Section 6.1 we illustrate the situation numerically. Alternatively, for great |i − j|, γ i,j can be computed in a numerically stable way by quadratures using standard arithmetics. The accuracy of quadratures rapidly grows as |i − j| increases; see Section 5.3. The quadrature rules are based on another way for an exact computation of β i,j described in Section 5.2.

Exact algorithm 2
In case of equation ( it is a polynomial of degree m. Integration by parts yields .

Product Quasi-Interpolation in Logarithmically Singular Integral Equations707
The integrals in (5.4) can be computed exactly since the integrands are polynomials of degree m − 1:  ,i−j ( + 1) are great, and again a loss of accuracy takes place in standard arithmetics.

A numerically stable approximate algorithm
We obtain a numerically stable algorithm approximating the integrals in (see (5.3)) by the µ-point Gauss quadrature, 2µ ≥ m. Fortunately, as we see below, the accuracy of the quadrature increases as |i − j| increases, so we obtain a suitable complement to the exact formulae for β i,j presented in Sections 5.1 and 5.2. Recall that B m (x) ≥ 0. For even m, also log |(i − j + m 2 − x)h| is signconstant for x ∈ ( , + 1), and the Gauss quadratures are stable numerically for the integrals in (5.6). For odd m, the integrand changes its sign when x = i − j + m 2 ± n = + 1 2 but there is still no danger for the numerical instability.
The accuracy can be even raised if we treat B m (x) in the integrals (5.6) as a weight function and apply corresponding Gauss type quadrature. Unfortunately, this way is labour consuming.
It is reasonable to compute β i,j by exact formulae (5.1)-(5.2) or (5.4)-(5.5) for |i − j| ≤ m or perhaps for |i − j| ≤ 2m and to apply the eight point Gauss quadrature in (5.6) for |i − j| > m, respectively, |i − j| > 2m, all in the framework of a standard arithmetics.

Numerical Experiments
Here we are testing numerically the properties of the methods described above. We start with investigating the calculation schemes for different ways of finding β i,j for the choice of the best strategy from the three different methods. For verifying the accuracy of the computed results of β i,j we use interval arithmetics with a specified number of bits of precision. We start with testing numerical properties of the exact calculation methods of β i,j (algorithms 1 and 2) in interval arithmetics with variable number of bits in precision. The third algorithm, namely the approximate integration method, is performed in standard double precision arithmetics.
We demonstrate that for best numerical accuracy in the case of double precision arithmetics different given algorithms need to be combined. Thereafter we set up a model problem with a known solution for which we observe, that the best way of combining the exact and the approximate scheme, does not depend only on i and j in calculation of β i,j , but also on the size of discretization parameter n and spline rank m. Nevertheless, the rules of thumb formulated in the end of Section 5 give quite acceptable accuracy.

Calculation of β i,j accuracy assessments using variable precision interval arithmetics
Here we compare the three different methods given above, namely, the exact method 1 (E1) for β i,j calculation described in Section 5.1, the exact method 2 (E2) given in Section 5.2 and the approximate integration method (AIM) outlined in Section 5.3. For testing the accuracy hypotheses given above, we perform cross-comparison tests with one of the methods in standard 53-bit interval arithmetics (or double precision) and the other method calculated using higher precision than normal, if needed, to achieve the maximal diameter of the result in interval arithmetics less than 10 −15 . To reach this goal we gradually increase the number of bits in the variable precision, starting from 53, until the goal is achieved. Standard double precision numbers are represented with 53 bits (one bit is reserved for the sign and 52 bits for the mantissa). While using interval field for calculation, the worst possible influence of round-off errors is given by the diameter of the calculated result. This makes it possible to numerically find out the algorithm behavior in case of real calculations with given input data. This technique is applied separately to all combinations of indexes i and j in β i,j calculation resulting in method accuracy across different values of i and j indicated with the minimum number of precision bits needed for the result to achieve the needed accuracy. (Due to the large number of computations needed for such calculations, we perform the tests only with moderate n.)  Fig. 1 we plot the error of the result, which is calculated as distance from β E2(1e−15) i,j . Here, we denote with E2(1e − 15) the algorithm E2, which has been calculated with the  warranty of 1e − 15 maximal diameter in the result. This is achieved through gradually increasing the number of bits in precision of real interval field, that results in diameter of the result being at most 1e − 15, starting over for each combination of indexes i and j.
In addition, in Tab. 1 we give the minimum numbers of precision bits for reducing the interval diameter of the resulting values to 10 −15 . As we can see, the method is most accurate in the case i = j and worst accurate with highly distinct values of i and j.

Exact method E2 for calculation of β i,j accuracy
Similarly, in Fig. 2 we present the results of the accuracy calculation for the method E2. The corresponding numbers of precision bits needed for achieving E2(1e − 15) are given in Tab. 2.

Approximate integration method (AIM) accuracy
Here, we compare the result of AIM with E1(1e − 15). The result is shown in Fig. 3. We observe that AIM behaves numerically very much in an opposite manner to the exact methods given above. While the methods E1 and E2 give the most accurate value for β i,j in cases when indexes i and j are not too distinct from each other, AIM works best in case of |i − j| being not too small. Due to better performance, we will perform further experiments with E1 from the two exact methods.
Yet we do not know, what would be a reasonable size in |i − j| for switching from one of the exact methods to AIM. This we are checking out together with the overall method verification with the following model problem.     in Tab. 3. As we can see, the strategy for choosing the switching value k depends also on n. Two example plots of the error depending on k are presented in Fig. 4. But if we use a higher order spline, say m = 8 (m 0 = 3, m 1 = 5), we have found out that with n = 512 the minimal error 3.49446166e-13 is achieved for k = 8. This means that if one really wants to achieve the absolutely best calculation scheme by combining the given β i,j algorithms, it depends on many parameters. Although, usually a general rule of choosing k somewhere between m and 2m should give quite good results, it follows from the results presented in Fig. 4, that the error is practically constant for k > m. Thus the selection of k seems to be not a very sensitive step for controlling accuracy of the proposed numerical algorithm, at least for moderate n.

Model problem solution
From the numerical examples we discovered that the exact method E1 worked actually better than we were originally expecting.
For all the numerical experiments we used the sage worksheet environment (see e.g. [12]) using numpy and scipy packages. The source code of the programs together with the test results with interactive graphics can be found at http://www.ut.ee/∼eero/LogSingIntEq/.

Conclusion
We examined fully discrete methods for solving linear second kind integral equations containing boundary and logarithmic diagonal singularity in the kernel. The methods are based on the smoothing change of variable and the product quasi-interpolation by smooth splines of order m on the uniform grid of step size h = 1/n. The methods are of optimal accuracy order O(h m ). The matrix form of the methods is given by (4.1)-(4.7). The main difficulty in the application of the methods is the computation of quadrature coefficients β i,j defined in (4.7). We presented a simple exact algorithm (5.1)-(5.2) for the computation of β i,j , which due to some numerical instability can be recommended for small and moderate |i − j|, say for |i − j| ≤ 2m. For greater |i − j| we we recommended to use quadratures described in Section 5.3, although our numerical example encourages us to recommend exact formula (5.1)-(5.2) also for much greater |i − j| than 2m; the preciseness of quadratures rapidly grows as |i − j| increases. In engineer computations with n m /m! ≤ 10 8 , formulae (5.1)- (5.2)