Singularity Subtraction for Nonlinear Weakly Singular Integral Equations of the Second Kind

The singularity subtraction technique for computing an approximate solution of a linear weakly singular Fredholm integral equation of the second kind is generalized to the case of a nonlinear integral equation of the same kind. Convergence of the sequence of approximate solutions to the exact one is proved under mild standard hypotheses on the nonlinear factor, and on the sequence of quadrature rules used to build an approximate equation. A numerical example is provided with a Hammerstein operator to illustrate some practical aspects of effective computations. 1.


Introduction
To be consistent with [An81] and [AhEtAl01], we assume that g is a decreasing function on ]0, b − a]. The factor N, containing the values x(t) ∈ R of the functional variable x ∈ X for t ∈ [a, b], is a continous function with continuous partial derivative with respect to the third variable. Then the operator K maps X into itself, and it is Fréchet-differentiable over X. When N(s,t, x(t)) := κ(s,t) x(t) for some continuous function κ : [a, b]×[a, b] → R, then K is a linear bounded operator from X into itself. In this paper, we are interested in the general, possibly nonlinear, case.
The main idea of the singularity subtraction method is to compensate the singularity of g(|s−t|) along the diagonal s = t, by multiplying g(|s−t|) by a factor which tends to 0 as t → s. If K is linear, this factor is κ(s,t)(x(t) − x(s)). In the general case, the factor is N(s,t, x(t)) − N(s, s, x(s)). This leads to rewrite K as (1.1) The singularity subtraction method builds an approximation of K as it is written in (1.1), and, as described in [An81] for the linear case, it is a double approximation scheme consisting of truncation and numerical integration.
The ideas worked out in [An81] and [AhEtAl01] for the linear case, are extended here to the nonlinear case. Truncation: Given δ ∈ ]0, b − a[, we replace g with a truncated approximation g δ in a δ -right-neighborhood of 0. This function coincides with g outside a small interval [0, δ ], and is constantly equal to g(δ ) in [0, δ ]. Hence g δ is a continuous function. In the sequence of singularity subtraction approximations, the role of δ is played by a sequence (a n ) n≥2 in ]0, b − a[ leading to the function g n defined by g n (s) := g(a n ) for s ∈ [0, a n ], g(s) for s ∈ ]a n , b − a].
Numerical integration: To proceed with the singularity subtraction idea -like in the linear case -we define a general grid with n ≥ 2 points on [a, b]: a ≤ τ n,1 < τ n,2 < . . . < τ n,n ≤ b. (

1.2)
This grid is called the basic grid, and it determines n − 1 subintervals of [a, b]. The integrals in the first line of (1.1), after replacing g with g n , are approximated by some quadrature rule Q n with p(n) nodes depending on the nodes of the basic grid. For instance, if Q n is the composite trapezoidal rule, then the quadrature grid is the basic grid, so p(n) = n; if Q n is the composite Simpson rule, then its nodes are the points of the basic grid and the mid-points of the corresponding subintervals, and hence p(n) = 2n − 1; for some other rules Q n , the nodes are the so-called Gaussian points, which are obtained by shifting to each subinterval of the basic grid the zeros of a polynomial of a given degree m belonging to a complete sequence of orthogonal polynomials in some particular Hilbert space, and hence p(n) = m(n − 1).
In this paper, we consider a sequence (Q n ) n≥2 of quadrature rules built upon the basic grid. The nodes of Q n are denoted by t p(n), j , j = 1, . . . , p(n), and are numbered so that a ≤ t p(n),1 < · · · < t p(n),p(n) ≤ b. The weights of Q n are denoted by w p(n), j , j = 1, . . . , p(n). We suppose that they are all positive, and that there exists a constant γ > 0 such that

Singularity subtraction
Consider the basic grid (1.2) and define h n, j := τ n, j+1 − τ n, j for j = 1, . . . , n − 1, and h n := max j=1,...,n−1 h n, j . The singularity subtraction technique, as presented in [An81], relates truncation and numerical integration through the following condition on the sequences (a n ) n≥2 and (Q n ) n≥2 : There exist constants α 1 > 0 and β 1 > 0 such that α 1 h n ≤ a n ≤ β 1 h n for all n ≥ 2, i.e. the width of truncation must tend to zero at the same rate as the mesh sizes.
These considerations lead to approximate K, as written in (1.1), by the following operator K n : For x ∈ X, and s ∈ [a, b], The exact equation, to be solved numerically, is: For y ∈ X, find ψ ∈ X such that We assume that 1 is not in the spectrum of the Fréchet-derivative of K at ψ, so ψ is an isolated solution of (1.4).
The approximate equation, to be solved exactly, is: Find ψ n ∈ X such that If we take the values of (1.5) at t p(n),i , i = 1, . . . , p(n), we get the following, possibly nonlinear, system of order p(n) with unknowns x n (i) := ψ n (t p(n),i ), i = 1, . . . , p(n): This p(n)-dimensional system can be written as where, for all x ∈ R p(n)×1 , and i = 1, . . . , p(n), The system (1.6) must be solved accurately by some numerical method like, for instance, Gauss' method in the linear case, and Newton's method -as described in the sequel -in the nonlinear case.
The Jacobian matrix of F n : R p(n)×1 → R p(n)×1 at x ∈ R p(n)×1 is given by where δ i, j is the Kronecker delta, and i, j = 1, . . . , p(n).
The Newton's sequence x n k≥0 in R p(n)×1 is defined, for a given starting column x n is the unknown. Equivalently, where I is the identity matrix of order p(n), and, for i, j = 1, . . . , p(n), Theorem 1 Let (Q n ) n≥2 be a sequence of composite quadrature rules with nodes t p(n), j and weights w p(n), j , j = 1, . . . , p(n), satisfying (1.3). Then K n n≥2 is pointwise convergent to K, F n n≥2 is pointwise convergent to F , and ψ n n≥2 is convergent with limit ψ.
Proof. The Fréchet-derivatives T := K and T n := (K n ) at ψ are given by: Let us consider the decomposition T n (ψ) = T A n (ψ) + T B n (ψ), where The proof is done in 5 steps: 1. We show that T n (ψ) p → T (ψ) and that T n (ψ) ν → T (ψ) too: As T n (ψ) and T (ψ) are bounded linear operators, we use the results of [An81].
Since (s,t, u) → ∂ N(s,t, u) ∂ u is a continuous function, and since (1.3) holds, then T n (ψ) and T A n (ψ) satisfy the hypotheses of Proposition 4.18 in [AhEtAl01], page 227, and T A n (ψ) .
3. We prove that F n is locally invertible with continuous inverse in a neighborhood of 0: I − K n is a continuously differentiable operator from the Banach space X into itself. By the Inverse Function Theorem, I − T n (ψ), being invertible, there is a neighborhood of ψ where I − K n is invertible with continuous inverse in some neighborhood of y. Hence F −1 n exists and is continuous in some neighborhood of 0. 4. We prove that (K n ) n≥2 is pointwise convergent to K, and (F n ) n≥2 is pointwise convergent to F : An auxiliary operator K n is used in the proof. For x ∈ X, and s ∈ [a, b], define w p(n), j g n (|s − t p(n), j |)N(s,t p(n), j , x(t p(n), j )).
K n can be rewritten as K n (x)(s) = K n (x)(s) + N(s, s, x(s))(U −U n )e(s). (1.9) which is finite because of the continuity of N in its three variables, and of x in its single one. In the linear case, σ (x) = ρ x for some constant ρ > 0. Now, since U n p → U as stated in the beginning of this proof. Following the ideas of the proof of Proposition 4.18 in [AhEtAl01], we decompose where λ δ , µ n and η n are defined as follows. Let γ > 0 be the constant introduced in Then the following upper bounds hold for all n greater than some integer n 0 (x): we conclude that K n p → K. K n p → K, and F n p → F .

5.
We prove that (ψ n ) n≥2 is convergent with limit ψ: Since F and F n are invertible and Fréchet-differentiable, the derivative of their inverses at 0 is equal to the inverse of the derivative of the direct operators at the inverse image of 0, and the integral form of the Mean Value Theorem for Derivatives gives: Since the sequence (F n ) n≥2 is pointwise convergent to F and F (ψ) = 0, then v n (t) := F n (ψ) + t(F (ψ) − F n (ψ)) tends to 0 uniformly in t ∈ [0, 1] as n → ∞.

Programming:
The numerical implementation was carried on MATLAB 2017b.
Newton: Because of its fast convergence, we have performed 7 iterations of Newton's method for each fixed value of n.
Quadrature: The basic grid was chosen to be uniform. We used a composite quadrature rule built with 3 Gauss points on each one of the n − 1 subintervals of the basic grid: the zeros of a Chebyshev polynomial of degree 3. Hence, we had been led to a nonlinear system of order 3(n − 1) for each value of n, and to a linear system of order 3(n − 1) at each Newton's iteration. The numerical integration of the explicit integrals in (1.7) and (1.8) was performed by adaptive quadrature based on a Gauss-Kronrod method (cf. [Sh08]). The use of strict tolerances for the numerical adaptive quadrature approach is mandatory, otherwise finer discretizations would not deliver the predicted convergence results.

Conclusions
In this paper, we have extended to nonlinear integral operators, the singularity subtraction technique for approching linear weakly singular integral operators in the framework of real valued continuous functions. The singularity subtraction technique cannot be settled in Lebesgue spaces.
The equation to be solved numerically, written as F (ψ) = 0, is discretized as F n (ψ n ) = 0, which leads to a finite-dimensional nonlinear problem which is solved by Newton's method. We have proved that the sequence (F n ) n≥2 is pointwise convergent to F , and that the sequence (ψ n ) n≥2 is convergent with limit ψ. The double bound (1.10) shows that the convergence of the latter is neither slower nor faster than that of the sequence of (F n (ψ)) n≥2 to 0.
Since the values of ψ n are approximated by Newton's method only at the p(n) nodes of the quadrature grid, ψ n could be globally approximated by interpolation, if needed. In fact, oppositely to the linear case, once the unknowns x n (i) = ψ n (t p(n),i ) are approximated by the values x n (i) issued from Newton's method, no natural interpolation formula is available for a closed formula of ψ n since, for each s ∈ [a, b], ψ n (s) is hidden implicitly in the nonlinear expression (1.5). This is a significant difference between the linear case and the nonlinear case.
The expression (1.9) shows that the order of the pointwise convergence to K of the sequence of approximaitons (K n ) n≥2 is dominated by the order of convergence of the sequence of quadrature rules Q n n≥2 . In the case of the numerical computations presented in this paper, the sequence of quadrature rules converges in theory to the exact integral at the same rate as n −2 tends to 0, when n → ∞ (cf. [XiEtAl12]). This is confirmed in practice, as it is shown in the loglog plotting of Figure 1.1, where we observe that the slope of the straight line is around −2.
A major survey on numerical approximation of nonlinear integral equations is [At92]. This paper studies numerical methods for calculating fixed points of nonlinear integral operators, i.e. equations of the form ψ = K(ψ) with the notation of our paper. This corresponds to the case y = 0 and is less general than the work presented here since y cannot be incorporated as a part of the integral operator K. Methods treated in [At92] include a product integration type scheme for weakly singular Hammerstein operators, projection methods and Nyström methods. As in our paper, all those methods require the solution of finite-dimensional systems of non-linear equations. An auxiliary numerical method is needed to solve these nonlinear finite-dimensional systems.