An Approximation Method to Compute Highly Oscillatory Singular Fredholm Integro-Differential Equations

: This paper appertains the presentation of a Clenshaw–Curtis rule to evaluate highly oscillatory Fredholm integro-differential equations (FIDEs) with Cauchy and weak singularities. To calculate the singular integral, the unknown function approximated by an interpolation polynomial is rewritten as a Taylor series expansion. A system of linear equations of FIDEs obtained by using equally spaced points as collocation points is solved to obtain the unknown function. The proposed method attains higher accuracy rates, which are proven by error analysis and some numerical examples as well.


Introduction
Integro-differential equations (IDEs) are studied in physics, biology, and engineering applications, as well as in advanced integral equations literature.Applications of the IDEs in electromagnetic theory, dispersive waves, and ocean circulations are enormous [1][2][3][4][5][6].The addressed paper considers the highly oscillatory singular FIDEs of the following kind: along with the initial condition: The given highly oscillatory kernel function k(x, t) possesses weak and Cauchy singularities, i.e., k(x, t) = w(t)e ikt t−x , w(t) = (1 + t) α (1 − t) β .P(x), f (x), Q(x) are smooth functions on [−1, 1], whereas y(x) is the unknown function that needs to be determined.The FIDEs having weak or strong singularities are considered in [7].To gain analytical approximation, the author used the power series expansion technique.However, for the convergence of the method, the author applied a ratio test and proved that the proposed method gives exact solutions if the solutions of the equations are finite-degree polynomials.Otherwise, by increasing the number of polynomials, better accuracies are obtained.
For a FIDE of the kind a reproducing kernel method is introduced in [8], where the kernel function in the equation is a weakly singular kernel.The method converts the weakly singular kernel into a logarithmic kernel to a Kalman kernel.Furthermore, a smooth transformation helps to remove the weak singularity of the Kalman kernel.The authors claimed that the reproducing kernel method is not restricted by the order of the equation.
A quadrature formula is applied to discretize the FIDEs.Based on the behavior of the exact solution, the special graded grid points are used for piecewise polynomial collocation method [9].Smooth parts of the integrands are approximated by the piecewise polynomial interpolation, whereas the singular parts are integrated exactly.By the same authors in [10], two approaches, an integral equation reformation and discrete Galerkin method, are used to find the approximations for the solutions and derivatives of the nth order weakly singular FIDEs.The approximations to the solutions are piecewise polynomial functions.In another research work, the Taylor series expansion along with the Galerkin method is considered for FIDEs with Cauchy singularity.The Legendre polynomials are used as basis to approximate the solution of the FIDEs [11].The traditional piecewise homotopy perturbation method is extended for FIDEs with weak singular kernels.The accuracy and calculation speed with the Gauss quadrature rule and piecewise low-order interpolation is significantly improved in [12].
For an FIDE with weak singular or other non-smooth kernels, the regularity properties of the solutions are briefly studied by the author in [13].These obtained results are further used in the analysis to solve such problems by piecewise polynomial collocation method.As compared to highly oscillatory Volterra integral equations, Fredholm integral equations with high oscillations have received less attention by researchers.Nevertheless, there is still a gap to be filled to propose the approximate methods to obtain the numerical solutions of highly oscillatory FIDEs along with the singular kernel functions.The purpose of this paper is to represent an efficient algorithm for such highly oscillatory singular FIDEs and try to fill this gap.The main novelty is the simplicity and accuracy of the proposed method.This research work aims at introducing an approximation method for a highly oscillatory Fredholm integro-differential equation.The general form of the integral term in Equation ( 1) with highly oscillatory function, weak singularities and Cauchy singularity in . Along with weak and Cauchy singularities, the highly oscillatory FIDE (1) is impossible to solve by the classical methods discussed above.To overcome such adversity, an approximated method following the steps of the Clenshaw-Curtis rule is presented, which interpolates the unknown function with an interpolation polynomial H N (x) of degree N. In the integral term, the interpolation polynomial extended to a sum of Chebyshev series of the first kind is rewritten in the form of truncated Taylor series expansion for x ∈ (−1, 1) .The pivotal point of this method is that efficient results with higher accuracy rates can be achieved for smaller values of N and m, where m defines the degree of the Taylor series expansion of the Chebyshev polynomial.With the evaluation of the highly oscillatory singular integral, a system of linear equations is constructed for Equation (1) by using the equally spaced points as the collocation points.This system is eventually solved to obtain the unknown coefficients for the unknown function y(x).
The rest of this paper is organized as follows: Section 2 illustrates the methodology of a newly proposed method briefly, and obtains a system of linear equations.Section 3 gives an idea of error estimation.The efficiency and accuracy of the method is demonstrated in Section 4 by some numerical examples.

Methodology
A function y(x) can be approximated by its interpolation polynomial H N (x) of degree N at Chebyshev points of the second kind, x j = cos jπ N , j = 0, 1, ⋯N.We rewrite the interpolation polynomial in terms of the Chebyshev series as [14]: where T n (x) is the Chebyshev polynomial of the first kind of degree n, double prime denotes a sum whose first and last term are halved and a n are the coefficients which can be efficiently computed by FFT [15,16], defined as Equation ( 5) can further be written in the matrix form as follows: Furthermore, the derivative of a function y (µ) (x), µ = 0, 1, ⋯, in terms of Chebyshev series and matrix form is defined as [17] y (µ) (x) = T(x)A (µ) , where the matrix M (N+1)×(N+1) for odd and even values of N is defined, respectively, as: .
By considering the Equation ( 5) for integral ⨍ , it transforms into the following: where R n (x, k) are called the modified moments and defined as: The following subsection appertains the Clenshaw-Curtis method with truncated Taylor series to evaluate the moments R n (x, k) accurately.

Computation of the B n (x, k) Moments
For moments B n (x, k), we transform these moments into a singular and non-singular integrals by applying some basic steps as follows: For the Cauchy singular integral ⨍ x dt, an explicit calculation has been performed in [16] as where t dt is the exponential integral, and sgn(x) is the sign function.On the other hand, by applying the truncated Taylor series expansion for T n (t) in non-singular integral, we obtain where To calculate the moments A m , a simple recurrence relation is illustrated by integrating by parts with the initial value A , where m denotes the derivatives of the relative terms.

Computation of the B n (x, k, α, β) Moments
For the moments dt, by following the initial necessary steps as above, we obtain In addition, by applying the truncated Taylor series expansion for the left integral, it becomes where Ãm clearly possesses some weak singularities as well as high oscillation and cannot be calculated easily without any numerical method.The following theorem presents the steepest descent method to evaluate this integral significantly accurately.
Theorem 1. Suppose that a function g(t) is analytic in the half-strip of the complex plane −1 ≤ R(x) ≤ 1 and I(x) ≥ 0, and satisfies that for M and k 0 constants, the integral ikt dt can be evaluated as Proof.The integrand g(t)w(t)(t − x) m e ikt is analytic in the half-strip of the complex plane: −1 ≤ R(x) ≤ 1 and I(x) ≥ 0; then, based on the Cauchy's theorem, we obtain where paths of the integration are taken in counterclockwise direction and shown in the Figure 1.Since where R is a large number, then Similarly, In addition, For Γ 0 , by taking t + 1 = re iθ and 0 For r → 0, there exist By Equations ( 19) and ( 18), we yield Similarly, for all integration paths Γ m → 0, r → 0. For r → 0, R → ∞, Equation ( 17) implies For the calculation of integrals involving weight functions t α e −t and t β e −t , the generalized Gauss-Laguerre quadrature rule can be used to approximate these integrals.For this purpose, let x j and w j be the zeros and weight functions of generalized Gauss-Laguerre quadrature rule, respectively, then these integrals are written as The above theorem is proved for a highly oscillatory weakly singular integral which has an analytic function; however, this can be extended to the integral Ãm of our choice.13) is also explicitly calculated by steepest descent method in [18,19].

The integral ⨍
With the successful evaluation of the moments R n (x, k), the Equation ( 1) is transformed after some necessary substitutions as ν µ=0 Furthermore, to obtain the unknown coefficients A in Equation ( 22), we can apply collocation method for equally spaced points x i as ν µ=0 where For unknown coefficient vector A, the Equation ( 23) is obtained in matrix form ν µ=0 where W is an (n + 1) × (n + 1) matrix and corresponds to a system of linear equations that can be solved to obtain the unknown coefficients A. Similarly for the initial conditions (2), we obtain Consequently, by replacing the rows of W with the rows of Equation ( 25), the new obtained system of equations WA = F is solved for the unknown coefficients A. Once the unknown coefficients are derived, these can be substituted in the Equation ( 5) to obtain the unknown function y(x).It should be noted that the matrix W of the linear system of equations WA = F (24) can be ill-conditioned, which indicates that increasing the parameter N does not guarantee the greater accuracy of the approximated solution.
To avoid such adversity, small values of N are taken as the proposed method, giving more accurate results for smaller N.The value of the N in numerical examples is taken as N ≤ 15.The increment in the N values gives no such improvement in the accuracy.

Error Analysis Theorem 2 ([15]). A Lipschitz continuous function f on [−1,1] has a unique representation as a Chebyshev series
which is absolutely and uniformly convergent.The coefficients are given for j ≥ 1 by the formula for j = 0, the factor 2 π changes to 1 π .
• If f is analytic with | f (z)| ≤ M in the region bounded by the ellipse with foci ±1 and major and minor semiaxis lengths summing to ρ > 1, then for each j ≥ 0, • For an integer κ 0 ≥ 0, let f and its derivatives through f (κ 0 −1) be absolutely continuous on [−1, 1] and suppose the κ 0 th derivative of bounded variation V κ 0 .Then for j ≥ κ 0 + 1, the Chebyshev coefficients of f satisfy: Theorem 3. Suppose that y(t) is an analytic function and satisfies the conditions of the Theorem 1; then, the error bound for the integral Proof.For integral ∫ +∞ 0 f (x)x σ e −x dx, σ > −1, the error formula for the l−point general Gauss-Laguerre quadrature rule is defined as [21] l!Γ(l By applying the above formula for the error E[y], we obtain For k ≫ 1, the error for E[y] decays asymptotically as , k ≫ 1, which proves the theorem. Theorem 4. For an analytic function y(t) satisfying the conditions of the Theorem 2 and Lemma 1, the error bound for integral term ⨍ dt for the proposed method is defined as Proof.The error between a function y(t) and its Taylor series at a point x is defined as For an analytic function y(t) and its interpolation polynomial H N (t), the error bound is considered by the Taylor's series expansion as follows: by applying the Theorem 3 and bounds for ∥y Remark 1.Another practical error estimation e (i) N | can be adduced for Equation (1), substituting its approximated polynomial interpolation r i=0 P i (x)H (i) where D N (x) is the residual function associated with H N (x).Subtracting the Equation ( 35) from (1), we obtain r i=0 P i (x)e (i) by solving this equation, the residual function should be zero D N (x) ≃ 0 or less than the error estimated in the above lemmas.

Numerical Examples
Example 1.In the following example, we have presented an absolute error estimation for the moments R n (x, k).The approximated values are calculated by the proposed method, and the exact values are obtained in Mathematica11 software.The Tables 1 and 2 illustrate the efficiency of the method, as for really small values of N, higher error estimations are derived.The proposed method is simple and gives good results.Table 1 shows the absolute error for the moments T n (t)e ikt t−x dt for x = 0.2.However, Table 2 represents the absolute error for the moments B n (x, k, α, β) = ⨍ where the value of f (x) is taken as f (x) = − sin(x) + x cos(x) − x sin(x) + ⨍   Example 3. A second-order weak singular FIDE dt, the absolute error is presented at equally spaced points in Table 4 and Figure 3.   the error for w(t) = (1 + t) α (1 − t) β , α = −0.5, β = −0.5, is shown in Table 5, and for larger values of k, Figure 4 shows the absolute error for N = 10.dt, the absolute error is presented at equally spaced points in Table 8 and Figure 6 for small values of N, whereas k is taken to be large enough.

Conclusions
In this paper, we have given a really simple but efficient method to solve the highly oscillatory singular FIDEs.For a small number of equally spaced points x i , i = 0, 1, 2, ⋯, N as collocation points, the proposed method provides efficient higher accuracy.We do not need to take the higher mth degree Taylor series expansion of the Chebyshev polynomial, as even for m = 1 we obtain satisfactory approximation to the exact values.However, for weak singularities, the w(t) ≠ 1 values of m and N are taken to be identical.All the exact values of the integrals are calculated in Mathematica11 software, whereas the code for the proposed method is developed in MatlabR2018a.

Figure 1 .
Figure 1.The integration path for integral I[g].

1 − 1
sin(t)e ikt t−x dt.The error for this equation is provided for N = 2, N = 5 in Table3.It is shown that even with small values of m, high accuracy can be achieved for large values of k, i.e., k = 10 4 .Figure2also presents this phenomenon.

Table 4 .
The absolute error for N = 5.

Figure 3 .Example 4 .
Figure 3.The absolute error for N = 5, m = 5.Example 4. For FIDE of the following form

Table 1 .
The absolute error for B n (x, k).

Table 3 .
The absolute error for k = 10