A Fitted Operator Finite Difference Approximation for Singularly Perturbed Volterra–Fredholm Integro-Differential Equations

: This paper presents a ε -uniform and reliable numerical scheme to solve second-order singularly perturbed Volterra–Fredholm integro-differential equations. Some properties of the analytical solution are given, and the ﬁnite difference scheme is established on a non-uniform mesh by using interpolating quadrature rules and the linear basis functions. An error analysis is successfully carried out on the Boglaev–Bakhvalov-type mesh. Some numerical experiments are included to authenticate the theoretical ﬁndings. In this regard, the main advantage of the suggested method is to yield stable results on layer-adapted meshes.


Introduction
Volterra-Fredholm integro-differential equations (VFIDEs) have led to many scientific computings. They play important roles in different branches of science involving aerodynamics, the economy, electricity and electronics, industrial networks, hydrodynamics, oceanography and chemistry [1][2][3] (see the references detailed within). Particularly, VFIDEs have been used widely for population growth, medicine processes and pandemic research. For example, the dissipation of tumor cells and the response of immune system were modeled in [4]. The effect of the COVID-19 pandemic was investigated in Italy, Germany and France with the help of modeling by some integro-differential equations [5,6].
This article concerns with boundary-value problem of second-order Volterra-Fredholm integro-differential equation in the form subject to boundary conditions u(0) = A, u(l) = B, (2) in which the differential operator Lu, the Volterra integral operator T and the Fredholm integral operator S are given as follow, respectively: and Additionally, ε ∈ (0, 1] is a small parameter,Ī = [0, l], λ is a real parameter, the functions a(x) ≥ α > 0, f (x) (x ∈Ī), u(x) (x ∈Ī), K 1 (x, t) and K 2 (x, t) ((x, t) ∈Ī ×Ī) are sufficiently smooth. Under these conditions, the problem (1) and (2) has unique solution u. As ε tends toward zero, the boundary layers appear in neighborhood of x = 0 and x = l.
In recent times, a lot of papers have been published about singularly perturbed integro-differential equations and various numerical schemes have been suggested. Iragi and Munyakazi have considered fitted operator finite difference method by using right-side rectangle rule and trapezoidal integration on Shishkin mesh for Volterra integro-differential equations with singularity [27,28]. In [29][30][31], Volterra delay integro-differential equations with initial layer have been investigated on uniform mesh. Mbroh et. al. have proposed non-standard finite difference scheme by using composite Simpson's rule. Additionally, they have improved the order of convergence by applying Richardson extrapolation [32]. In [33], second-order discretization have been presented on piecewise uniform mesh. Tao and Zhang have introduced the coupled method involving local discontinuous Galerkin technique and continuous finite element method in [34]. In [35], using composite trapezoidal rule, fitted mesh finite difference schemes have been established on Shishkin type mesh. Moreover, almost second-order accuracies for the presented method have been obtained. Exponentially fitted difference schemes have been suggested for singularly perturbed Fredholm integro-differential equations in [36][37][38]. In [39], Durmaz and Amiraliyev have constructed fitted second-order homogeneous difference scheme on Shishkin mesh for Fredholm integro differential equations with layer behavior. Authors in [40,41] have presented a new discrete scheme for singularly perturbed Volterra-Fredholm integrodifferential equations.
To the best of our knowledge, the problems in (1) and (2) have not been investigated using the finite difference schemes. Therefore, this study aims to fill this gap. This paper introduces the new difference scheme for the boundary value problems of second-order singularly perturbed Volterra-Fredholm integro-differential equations as the major novelty of this work. The second contribution is the convergence analysis of the presented scheme on Boglaev-Bakhvalov-type mesh. Last but not least, the proposed algorithm is easy to construct, and it provides stable results in a short time in terms of computation. From these objectives, the theory and applications of the presented method have been extensively studied.
The remainder of this article is organized is as follows. In Section 2, first, some preliminary results are given. Then, using composite numerical quadrature rules and implicit difference rules, the finite difference scheme is constructed on Boglaev-Bakhvalovtype mesh in Section 3. Section 4 is devoted to error approximations and stability analysis. In Section 5, three numerical examples are solved by the proposed method. Furthermore, the corresponding algorithm and the computational results are presented. Finally, the paper ends with "Concluding Remarks".

Asymptotic Properties
This section is devoted to some a priori bounds. For this aim, the following lemma is expressed.
Then, the solution u(x) of the problems in (1) and (2) holds Proof. By considering the maximum principle for the problems in (1)-(2), we obtain Then, it follows that Finally, we have By taking into consideration (3), we reach the proof of the (4). Now, we prove the relation (5). Since |K 1 | ≤K 1 , |K 2 | ≤K 2 and |u(x)| ≤ C, we can write for the second derivative of u(x) For any function g ∈ C 2 [0, l], the following formula is used to estimate |u (0)| and |u (l)| : where In (6), by taking g(x) = u(x), x = 0, α 0 = 0 and α 1 = ε, we find the estimation of |u (0)| : In the same way, rewriting g(x) = u(x), x = l, α 0 = l − ε and α 1 = l in (6), we obtain Differentiating (1), we have From (4), it is clear that We investigate the solution of the problems in (7) and (8) in the following form: Here, the functions v 1 (x) and v 2 (x) are the solutions of the following problems, respectively: By using the maximum principle, we obtain and Here, the function w(x) is the solution of the following problem: For the solution of the problems in (16) and (17), it is obvious that From here, we can write Thus, we obtain which hints at the proof of the relation (5) [42]. Therefore, the proof of the lemma is completed.

Discrete Scheme
In this section, the finite difference discretization is presented for the problem (1) and (2). First, we give the definition of the mesh. Let ω N be a non-uniform mesh on [0, l]: Here, we use the non-uniform mesh called Boglaev-Bakhvalov-type mesh in [43]. The transition point is taken as For an even number N, we divide each of the subintervals [0, σ 1 ], [σ 1 , σ 2 ] and [σ 2 , l]. Here, σ 2 = l − σ 1 . x i node points are specified as Before constructing the difference scheme, we define some notation for the mesh functions. For any mesh function v(x) defined onω N , we use the following implicit difference rules: Here,h i is defined ash and the discrete maximum norm is denoted by To establish the difference scheme for the problems in (1) and (2), we use the following integral identity: where the basis function is defined as follows: Moreover, it can be easily seen that For the differential operator Lu in (19), after using interpolating quadrature rules in [44] and some manipulations, we find where and Here, for s = 0, T 0 is computed as Moreover, for the right-side of the relation (19), we obtain with the remainder term given by For the Volterra operator in the relation (19), using interpolating quadrature rules in [44], we obtain After applying the composite right side rectangle rule, we have Then, we can write that Here, Similarly, for the Fredholm operator in the relation (19), applying the interpolating quadrature rules in [44] and the composite right side rectangle rule, it is found that where Combining (20), (22), (26) and (27), the following difference scheme is written for the problems in (1) and (2): where By omitting the error term R i in (28), we present the following difference problem for the approximate solution: where

The Stability and Convergence
Let the error function z i = y i − u i , i = 0, 1, 2, ..., N be the solution of the following problem: Here, and the remainder term R i is denoted by (29).
the solution of the problem (32) and (33) satisfies that Proof. Applying the discrete maximum principle to the discrete problem (32) and (33), we get From here, we can write Therefore, we arrive at the proof of the lemma.

Lemma 3.
For the remainder term R i , the following estimate is satisfied: Proof. By applying the mean value theorem to function a(x) in R Thus, taking into account a ∈ C 1 [0, l] and (35), it is found that Similarly, we can write Additionally, because of the boundedness of T 0 and |ϕ i (x)| ≤ 1, it is obvious that For the remainder term R (4) i , using the Leibnitz rule for the integral term in (24), we have Since ∂K 1 ∂x ≤ C and |u(x)| ≤ C, we can find Thus, it is seen that R In a similar way, we can show For the remainder term R i , applying the Leibnitz rule to the integral term in (25), we obtain Analogously, we have By substituting (36), (37), (38), (39), (40), (41) and (42) into (29), we estimate that Now, we consider the node points of adaptive mesh. The mesh stepsizes hold Then, we estimate the remainder terms for each sub-intervals separately. For the interval [0, σ 1 ], if σ 1 < l 4 , it is found that Thus, we obtain If σ 1 = l 4 , it can be written that Now, we consider the interval [σ 1 , σ 2 ]. We have For σ 1 < l 4 , we obtain h i ≤ lN −1 , and for σ 1 = l 4 , we obtain h i = lN −1 . Performing similar operations on the interval [σ 2 , l], we find which concludes the proof of the lemma. (1) and (2), and let y be the solution of the discrete problems in (30) and (31). Then, the following estimate is satisfied that

Theorem 1. Let u be the solution of the problems in
Proof. The proof of the theorem can be derived from the previous two lemmas.

Results and Discussion
In this section, we test the numerical method on several examples. For this, the elimination method is used to obtain maximum pointwise errors and convergence rates. Then, the discretization (30) and (31) can be written as the following form: where Here, the coefficients of the elimination method are as follow [45]: The corresponding Algorithm 1 is given by Algorithm 1: To compute the numerical solution y i .

Example 1. Consider a particular problem
in which the exact solution is u(x) = e −x ε . The nodal maximum errors are specified as where u i is the exact solution and y i is approximate solution. Furthermore, the convergence rates are computed as follows p N = ln(e N /e 2N ) ln 2 .  From Figures 1 and 2, it can be seen that the maximal errors are concentrated within the boundary layers.

Example 2.
We take into account another problem −ε 2 u + (x + 1)u + 1 2 The exact solution of this equation is unknown. Since the exact solution is unknown, we use the double-mesh technique. The error approximations are indicated by i and the order of convergence is defined as follows Error approximations are shown in Table 2.    In Figures 3 and 4, from ε = 10 −2 to ε = 10 −12 , it is seen that numerical solution behaves stable.

Example 3. Examine the last problem
The exact solution of this problem is unknown. Thus, we apply the double-mesh principle again. The obtained results are summarized in Table 3. The improvement in the numerical solution within the boundary layers is illustrated in Figures 5 and 6.  In summary, from Tables 1-3, as N increases, maximum pointwise errors decrease and the order of uniform convergence of the presented difference scheme is almost 1. It can be concluded that the presented difference scheme yields stable and uniform numerical results. Moreover, from Figures 1-6, it is seen that the numerical solution curves converge to the coordinate axes for smaller values of ε. Thus, the proposed method seems to be effective for solving such problems.

Concluding Remarks
In this paper, we suggested a new and stable difference scheme for solving boundary value problems of singularly perturbed second-order Volterra-Fredholm integrodifferential equations. The stability and convergence of the method were examined in the discrete maximum norm. The effectiveness and reliability of the presented scheme are demonstrated by three numerical examples. The computed results reveal that the order of convergence was found as O N −1 . Namely, the scheme is first-order accurate. Although the proposed method is reliable, the order of convergence of the numerical scheme may need to be improved. In order to expand numerical research, experiments in which improvements in the order and different boundary conditions (i.e., nonlocal, integral, Robin, and mixed) can be considered. Data Availability Statement: The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

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

Abbreviations
The abbreviations in this paper are given below: