This work is licensed under the Creative Commons Attribution 4.0 International License.
Introduction
Singularly perturbed integro-differential equations with delay are identified by those problems in which the highest order derivative terms are multiplied by small parameters and involving delay terms. In general, the solutions of such equations exhibit multiscale phenomena. Within certain thin subregions, the scale of some derivatives is significantly large. These thin regions of rapid change are called, boundary or interior layers, as appropriate. Such type of equations occur frequently in fluid dynamics, elasticity, quantum mechanics, chemical reactor theory, and other applied areas [1,2,3,4]. A class of fractional integro-differential equations (FIDEs) under the variable order differentiation has been solved using utilize the Haar collocation method [5]. To investigate the dynamics of TB transmission, a TB transmission model with vaccination and time delay was developed in an article by Zhang et al. [6]. It is well known that standard numerical methods for solving singular perturbation problems fail to give accurate results. Therefore, it is important to develop parameter-uniform numerical methods for these problems.
In this paper, we consider the following singularly perturbed integro-differential equation with delay [7]:
Lu+b(t)u(t-r)+\int_{0}^{t}\left[ K(t,s)u(s)+L(t,s)u(s-r)\right] \mathrm{d} s=f(t),\ \ \ t\in I,u(t)=\varphi (t),\ \ \ t\in I_{0},
where
Lu\equiv \varepsilon u^{\prime }(t)+a(t)u(t),
and
I=(0,T]=\mathop\cup \limits_{p=1}^{m}I_{p}
, Ip = {t|rp−1 < t ≤ rp} with 1 ≤ p ≤ m, I0 = [−r,0] and rp = pr, 0 < ɛ ≪ 1 is a perturbation parameter and r is a constant delay, a(t) ≥ α> 0,b(t), f (t), ϕ (t),K(t,s) and L(t,s) sufficiently smooth functions are given thus satisfying certain regularity conditions to be specified. Under these assumptions, the initial value problem (1)–(2) in general has boundary layers on the right side of each point t = rp (0 ≤ p ≤ m − 1) for small parameter ɛ (see [7]).
Recently, there has been a growing interest in the numerical methods for the singularly perturbed integro-differential equations. The existence and unique solution of functional integro-differential equation with nonlocal condition and finite delay function has been investigated in [8]. In [9], the existence of nonlocal impulsive fractional integro-differential equations of order 1 < r < 2 with infinite delay has been investigated. Approximate controllability of a second-order Volterra-Fredholm stochastic differential integral system including delay and impulses is investigated by Ma et al. [10]. Wu and Gan [11] applied linear multistep methods with convergent quadrature rules to solve singularly perturbed Volterra delay-integro-differential equations and derived global error estimates of A(α)-stable linear multistep methods. Kudu et al. [7] proposed an implicit finite difference method on a piecewise uniform Shishkin-type mesh with appropriate quadrature rules for a singularly perturbed delay integro-differential equation and showed that the scheme is almost first-order convergent. Amiraliyev and Yapman [12] used a fitted mesh method to solve a singularly perturbed Volterra delay-integro-differential equation and indicated that the scheme is also almost first-order convergent. Yapman et al. [13] developed a numerical method of integral identities with the use of exponential basis functions and interpolating quadrature rules for a singularly perturbed nonlinear Volterra integro-differential equation with delay and proved that the scheme is first-order uniform convergent. Yapman and Amiraliyev [14] also design a fitted difference operator on a piecewise-uniform mesh by using the exponential basis functions and interpolating quadrature rules to solve a singularly perturbed Volterra integro-differential equation and showed the convergence order of the scheme is O(N−2 lnN).
In this paper, we develop a finite difference scheme for the singularly perturbed integro-differential equation with delay (1)–(2). The construction of the difference scheme is a combination of hybrid difference scheme on a Shishkin-type mesh and appropriate quadrature rules. Error estimates are obtained by using the truncation error estimate techniques and a discrete analogue of Grönwall’s inequality and it is shown that the scheme is O(N−2 ln2N) order convergent, which improves the numerical results given in [7, 12,13,14]. Numerical experiments are presented to support the theoretical result.
The rest of the paper is organized as follows. The asymptotic behavior of the exact solution is derived in Section 2. The discretization method is described in Section 3. Error analysis for the finite difference scheme is presented in Section 4. Finally, numerical experiments are given in Section 5. Conclusion is given by Section 6.
Note 1
Throughout the paper, C will denote a generic positive constant that is independent of ɛ and of the mesh. Note that C is not necessarily the same at each occurrence. To simplify the notation we set gi = g(ti) and gi−1/2 = g((ti−1 + ti)/2) for any function g(t). We use the (pointwise) maximum norm on the interval Ω by ||·||∞,Ω.
Properties of the exact solution
In this section we show the asymptotic behavior of the exact solution, which is needed in the construction of layer-adapted meshes and the analysis of numerical methods. Since the exact solution u(t) may not be smooth at points t = rp (0 ≤ p ≤ m − 1), the bounds of the exact solution are somewhat different from those in [7].
Lemma 2
Under the assumptions ϕ ∈ C3 (Ī0), a,b, F ∈ C3 (Ī) and K,L ∈ C3 (Ī×Ī), the solution u(t) of problem (1)–(2) satisfies the following estimates\left\vert u^{(k)}(t)\right\vert \leq C\left[ 1+\varepsilon ^{-k}e^{-\frac{\alpha \left( t-r_{p-1}\right) }{\varepsilon }}\right] ,\ \ \ \ \ \ \ \ \ \ t\in \bar{I}_{p},\ 1\leq p\leq m,\ 0\leq k\leq 3.
Proof
Our proof is given step by step. On the first intervalĪ1, equation (1) can be rewritten in the form
\varepsilon {u}^{\prime }(t)+a(t)u(t)=F_{0}(t),\quad \quad \quad \quad t\in I_{1},
where
F_{0}(t)=f(t)-b(t)\varphi (t-r)-\int_{0}^{t}\left[ K(t,s)u(s)+L(t,s)\varphi (s-r)\right] ds.
Then we have
\begin{eqnarray*} \left\vert F_{0}(t)\right\vert &\leq &\left\Vert f\right\Vert _{\infty , \bar{I}_{1}}+\left\Vert b\right\Vert _{\infty ,\bar{I}_{1}}\left\vert \varphi (t-r)\right\vert +\int_{0}^{t}\left[ \bar{K}\left\vert u(s)\right\vert +\bar{L}\left\vert \varphi (s-r)\right\vert \right] ds \\ &\leq &\left\Vert f\right\Vert _{\infty ,\bar{I}_{1}}+\left( \left\Vert b\right\Vert _{\infty ,\bar{I}_{1}}+r\bar{L}\right) \left\Vert \varphi \right\Vert _{\infty ,\bar{I}_{0}}+\bar{K}\int_{0}^{t}\left\vert u(s)\right\vert \mathrm{d}s, \end{eqnarray*}
where
{\bar{K}=\max_{\bar{I}\times \bar{I}}\left\vert K(t,s)\right\vert }
and
{\bar{L}=\max_{\bar{I}\times \bar{I}}\left\vert L(t,s)\right\vert }
. Thus, applying the maximum principle to problem (4) we can get
|u(t)|\leq \delta _{0}+\alpha ^{-1}\bar{K}\int_{0}^{t}\left\vert u(s)\right\vert ds,\quad \quad \quad \quad t\in \bar{I}_{1},
where δ0 = |ϕ (0)| +α −1 [||f||∞,Ī1+(||b||∞,Ī1+rL̄)||ϕ||∞,Ī0]. Hence, applying the Grönwall’s inequalty to (5), we have
|u(t)|\leq \delta _{0}e^{\alpha ^{-1}\bar{K}t},\quad \quad \quad \quad t\in \bar{I}_{1},
which implies the inequality (3) for k = 0 and t ∈Ī1 holds true. Differentiating (4), we have
\varepsilon \left( u^{\prime }(t)\right) ^{\prime }+a(t)u^{\prime }(t)=F_{1}(t),\ \ \ \ \ \ \ \ \ \ \ \ \ \ t\in I_{1},
where
\begin{eqnarray*}F_{1}(t) &=&f^{\prime }(t)-a^{\prime }(t)u(t)-b(t)\varphi ^{\prime}(t-r)-b^{\prime }(t)\varphi (t-r) \\ &&-\left[ K(t,t)u(t)+L(t,t)\varphi (t-r)\right] -\int_{0}^{t}\left[ \frac{ \partial }{\partial t}K(t,s)u(s)+\frac{\partial }{\partial t}L(t,s)\varphi (s-r)\right] \mathrm{d}s.\end{eqnarray*}
By applying the result (3) with k = 0, the assumptions on the functions a,b, f,K,L and ϕ to the above equality and equation (1), respectively, we have
\left\vert F_{1}(t)\right\vert \leq C,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t\in I_{1},
and
\left\vert u^{\prime }(0)\right\vert \leq C\varepsilon ^{-1}.
From (6) we have
u^{\prime }(t)=u^{\prime }(0)\exp \left( -\frac{1}{\varepsilon } \int_{0}^{t}a(s)\mathrm{d}s\right) +\frac{1}{\varepsilon } \int_{0}^{t}F_{1}(\tau )\exp \left( -\frac{1}{\varepsilon }\int_{\tau }^{t}a(s)\mathrm{d}s\right) \mathrm{d}\tau .
Taking into account (7)–(8) in (9) we have
\left\vert u^{\prime }(t)\right\vert \leq C\left( 1+\varepsilon^{-1}e^{-\alpha t/\varepsilon }\right) ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ t\in \bar{I}_{1},
which implies that the result (3) for k = 1 and t ∈Ī1 holds true.
Differentiating (6) both sides once and twice, respectively, we have
\varepsilon \left( u^{\prime \prime }(t)\right) ^{\prime }+a(t)u^{\prime\prime }(t)=F_{2}(t),\ \ \ \ \ \ \ \ \ \ \ \ \ \ t\in I_{1},
and
\varepsilon \left( u^{\prime \prime \prime }(t)\right) ^{\prime}+a(t)u^{\prime \prime \prime }(t)=F_{3}(t),\ \ \ \ \ \ \ \ \ \ \ \ \ \ t\in I_{1},
where
F_{2}(t)=F_{1}^{\prime }(t)-a^{\prime }(t)u^{\prime }(t),\ \ \ \ \ \ F_{3}(t)=F_{2}^{\prime }(t)-a^{\prime }(t)u^{\prime \prime }(t).
Using the same method as that for k = 1, we also can prove that the result (3) for k = 2,3 and t ∈ Ī1 holds true.
Similarly, by using the same method as that for t ∈Ī1 step by step, we also can prove the desired results (3) for Īp with 2 ≤ p ≤ m.
Discretization and mesh
Let ωN0 be any non-uniform mesh on Ī:
\omega_{N_{0}}=\left\{0=t_{0}<t_{1}< \cdots <t_{N_{0}}=T\right\}
with mesh sizes hi = ti − ti−1, which contains N mesh point at each subinterval Ip (1 ≤ p ≤ m):
\omega_{N,p}=\left\{t_{i}:(p-1)N+1\leq i\leq pN \right\},\ \ \ \ \ \ \ \ \ \ \ \ \ 1\leq p\leq m,
and consequently
\omega_{N_{0}}=\bigcup_{p=1}^{m}\omega_{N,p}
. We also denote
\omega_{N,p}^{*} =\left\{t_{i}: 1\leq i\leq pN \right\},\ \ \ \ \ \ \ \ \ \ \ \ \ 1\leq p\leq m.
In order to be ɛ -uniform convergent, a Shishkin-type mesh is used to discretize the interval I. For the even number N, the piecewise uniform mesh ωN,p divides each of the intervals [rp−1, σp] and [σp, rp] into N/2 equidistant subintervals, where the transition point σp, which separates the fine and coarse portions of the mesh, is obtained by
\sigma_{p}=r_{p-1}+\min \left\{r/2,\ \ 2\alpha^{-1}\varepsilon \ln N\right\}.
Hence, if
h _{p}^{(1)}
and
h _{p}^{(2)}
denote the mesh sizes in [rp−1,σp] and [σp,rp] respectively, we have
\begin{eqnarray*}h_{i}=\left\{ \begin{array}{ll}h_{p}^{(1)}=2(\sigma_{p}-r_{p-1})N^{-1}, & i\in I_{p}^{1}, \\ h_{p}^{(2)}=2(r_{p}-\sigma_{p})N^{-1}, & i\in I_{p}^{2}\end{array}\right.\end{eqnarray*}
for 1 ≤ p ≤ m, where
I_{p}^{1}=\left\{\left(p-1\right)N+1,\cdots, \left(p-1/2\right)N\right\}, \ \ \ \ I_{p}^{2}=\left\{\left(p-1/2\right)N+1,\cdots, pN\right\}.
For our analysis we shall assume that
\varepsilon \leq CN^{-1}
as is generally the case for the singularly perturbed problems. On this mesh, we propose a hybrid difference scheme for problem (1)–(2):
T^{N}U_{i}=\tilde{f}_{i},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 1\leq i\leq N_{0}, U_{i-N}=\varphi (t_{i}-r),\ \ \ \ \ \ \ \ \ \ \ \ \ \,0\leq i\leq N,
where
\begin{eqnarray} T^{N}U_{i} &=&\left\{ \begin{array}{l} \varepsilon \frac{U_{i}-U_{i-1}}{h_{i}}+a_{i-1/2}\frac{U_{i-1}+U_{i}}{2} +b_{i-1/2}\frac{U_{i-1-N}+U_{i-N}}{2} \\ +{\sum_{j=1}^{i-1}\left[ K(t_{i-1/2},t_{j-1})U_{j-1}+K(t_{i-1/2},t_{j})U_{j} \right] \frac{h_{j}}{2}} \\ +\left[ K(t_{i-1/2},t_{i-1})U_{i-1}+K(t_{i-1/2},t_{i-1/2})\frac{U_{i-1}+U_{i} }{2}\right] \frac{h_{i}}{4} \\ +{\sum_{j=1}^{i-1}\left[ L(t_{i-1/2},t_{j-1})U_{j-1-N}+L(t_{i-1/2},t_{j})U_{j-N}\right] \frac{h_{j}}{2 }} \\ +\left[ L(t_{i-1/2},t_{i-1})U_{i-1-N}+L(t_{i-1/2},t_{i-1/2})\frac{ U_{i-1-N}+U_{i-N}}{2}\right] \frac{h_{i}}{4}, \ \ \ \ \ i\in I_{p}^{1}, \\ \varepsilon \frac{U_{i}-U_{i-1}}{h_{i}}+a_{i}U_{i}+b_{i}U_{i-N} \\ +{\sum_{j=1}^{i}\left[ K(t_{i},t_{j-1})U_{j-1}+K(t_{i},t_{j})U_{j}\right] \frac{h_{j}}{2}} \\ +{\sum_{j=1}^{i}\left[ L(t_{i},t_{j-1})U_{j-1-N}+L(t_{i},t_{j})U_{j-N}\right] \frac{h_{j}}{2}}, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ i\in I_{p}^{2}, \end{array} \right. \end{eqnarray}
and
\tilde{f}_{i}=\left\{ \begin{array}{ll}f_{i-1/2}, & i\in I_{p}^{1}, \\ f_{i}, & i\in I_{p}^{2}.\end{array}\right.
Analysis of the method
To investigate the convergence of the method, note that the error function zi = Ui − ui for 0 ≤ i ≤ N0 is the solution of the following discrete problem
T^{N}z_{i}=R_{i},\ \ \ \ \ \ \ \ \ \ \ \ \ 1\leq i\leq N_{0},z_{i-N}=0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \,0\leq i\leq N,
where
\begin{eqnarray} R_{i} &=&\left\{ \begin{array}{ll} \varepsilon \left( u_{i-1/2}^{\prime }-\frac{u_{i}-u_{i-1}}{h_{i}}\right) +a_{i-1/2}\left( u_{i-1/2}-\frac{u_{i-1}+u_{i}}{2}\right) & \\ +b_{i-1/2}\left( u_{i-1/2-N}-\frac{u_{i-1-N}+u_{i-N}}{2}\right) +\int_{0}^{t_{i-1/2}}K(t_{i-1/2},s)u(s)\mathrm{d}s & \\ -{\sum_{j=1}^{i-1}}\left[ K(t_{i-1/2},t_{j-1})u_{j-1}+K(t_{i-1/2},t_{j})u_{j} \right] \frac{h_{j}}{2} & \\ -\left[ K(t_{i-1/2},t_{i-1})u_{i-1}+K(t_{i-1/2},t_{i-1/2})\frac{u_{i-1}+u_{i} }{2}\right] \frac{h_{i}}{4} & \\ +\int_{0}^{t_{i-1/2}}L(t_{i-1/2},s)u(s-r)\mathrm{d}s & \\ -{\sum_{j=1}^{i-1}}\left[ L(t_{i-1/2},t_{j-1})u_{j-1-N}+L(t_{i-1/2},t_{j})u_{j-N}\right] \frac{h_{j}}{2 } & \\ -\left[ L(t_{i-1/2},t_{i-1})u_{i-1-N}+L(t_{i-1/2},t_{i-1/2})\frac{ u_{i-1-N}+u_{i-N}}{2}\right] \frac{h_{i}}{4}, & \\ & i\in I_{p}^{1}, \\ \varepsilon \left( u_{i}^{\prime }-\frac{u_{i}-u_{i-1}}{h_{i}}\right) +\int_{0}^{t_{i}}\left[ K(t_{i},s)u(s)+L(t_{i},s)u(s-r)\right] \mathrm{d}s & \\ -{\sum_{j=1}^{i}\left[ K(t_{i},t_{j-1})u_{j-1}+K(t_{i},t_{j})u_{j}\right] \frac{h_{j}}{2}} & \\ -{\sum_{j=1}^{i}\left[ L(t_{i},t_{j-1})u_{j-1-N}+L(t_{i},t_{j})u_{j-N}\right] \frac{h_{j}}{2}}, & i\in I_{p}^{2}. \end{array} \right. \end{eqnarray}
The following lemma gives the estimate of the truncation error.
Lemma 3
Under the assumptions ϕ ∈ C3(Ī0), a,b, f ∈ C3(Ī) and K,L ∈ C3(Ī×Ī), the truncation error of the difference scheme (12)–(13) satisfies:\left\vert R_{i}\right\vert \leq CN^{-2}\ln ^{2}N,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 1\leq i\leq N_{0}.
Proof
By using Taylor expansions for u and u′ around ti we have
\left\vert u_{i-1/2}^{\prime }-\frac{u_{i}-u_{i-1}}{h_{i}}\right\vert \leq \frac{3}{2}\int_{t_{i-1}}^{t_{i}}|u^{\prime \prime \prime }(s)|(s-t_{i-1}) \mathrm{d}s,|u_{i-1/2}-\frac{u_{i-1}+u_{i}}{2}| \leq \frac{3}{2} \int_{t_{i-1}}^{t_{i}}|u^{\prime \prime }(s)|(s-t_{i-1})\mathrm{d}s,
and
\left\vert u_{i}^{\prime }-\frac{u_{i}-u_{i-1}}{h_{i}}\right\vert \leq \int_{t_{i-1}}^{t_{i}}|u^{\prime \prime }(s)|\mathrm{d}s.
By using Taylor expansions for u around tj we also have
\begin{eqnarray} &&\left\vert \int_{t_{j-1}}^{t_{j}}K_{i}^{s}u(s)\mathrm{d}s-\left[ K_{i}^{j-1}(u_{j-1}+K_{i}^{j}u_{j}\right] \frac{h_{j}}{2}\right\vert \leq \left\vert \int_{t_{j-1}}^{t_{j}}K_{i}^{s}u(s)\mathrm{d} s-K_{i}^{j-1/2}u_{j-1/2}h_{j}\right\vert \notag \\ &&+\left\vert K_{i}^{j-1/2}u_{j-1/2}h_{j}-\left[ K_{i}^{j-1}u_{j-1}+K_{i}^{j}u_{j}\right] \frac{h_{j}}{2}\right\vert \leq 4h_{j}\int_{t_{j-1}}^{t_{j}}\left\vert \left( K_{i}^{s}u(s)\right) ^{\prime \prime }\right\vert \left( s-t_{j-1}\right) \mathrm{d}s. \end{eqnarray}
For i ∈
i\in I_{1}^{1}
we have
\begin{eqnarray} \left\vert R_{i}\right\vert &\leq &C\int_{t_{i-1}}^{t_{i}}\left( \varepsilon |u^{\prime \prime \prime }(t)|+|u^{\prime \prime }(s)|\right) (s-t_{i-1}) \mathrm{d}s+\int_{t_{i-1-N}}^{t_{i-N}}\left\vert u^{\prime \prime }(s)\right\vert \left( s-t_{i-1-N}\right) \mathrm{d}s \notag \\ &&+C\sum_{j=1}^{i-1}h_{j}\int_{t_{j-1}}^{t_{j}}\left[ \left\vert \left( K_{i-1/2}(s)u(s)\right) ^{\prime \prime }\right\vert +\left\vert \left( L_{i-1/2}(s)u(s-r)\right) ^{\prime \prime }\right\vert \right] \left( s-t_{j-1}\right) \mathrm{d}s \notag \\ &&+Ch_{i}\int_{t_{i-1}}^{t_{i-1/2}}\left[ \left\vert \left( K_{i-1/2}(s)u(s)\right) ^{\prime \prime }\right\vert +\left\vert \left( L_{i-1/2}(s)u(s-r)\right) ^{\prime \prime }\right\vert \right] \left( s-t_{i-1}\right) \mathrm{d}s \notag \\ &&+Ch_{i}\int_{t_{i-1}}^{t_{i-1/2}}\left[ \left\vert K_{i-1/2}^{i-1/2}\right\vert \left\vert u^{\prime \prime }(s)\right\vert +\left\vert L_{i-1/2}^{i-1/2}\right\vert \left\vert u^{\prime \prime }(s-r)\right\vert \right] \left( s-t_{i-1}\right) \mathrm{d}s \notag \\ &\leq &C\int_{t_{i-1}}^{t_{i}}\left( 1+\varepsilon ^{-2}e^{-\frac{\alpha s}{ \varepsilon }}\right) \left( s-t_{i-1}\right) \mathrm{d}s \notag \\ &&+C\sum_{j=1}^{i}h_{j}\int_{t_{j-1}}^{t_{j}}\left( 1+\varepsilon ^{-2}e^{- \frac{\alpha s}{\varepsilon }}\right) \left( s-t_{j-1}\right) \mathrm{d} s+CN^{-2} \notag \\ &\leq &C\left[ \int_{t_{i-1}}^{t_{i}}\left( 1+\varepsilon ^{-1}e^{-\frac{ \alpha s}{2\varepsilon }}\right) \mathrm{d}s\right] ^{2}+C\sum_{j=1}^{i}h_{j} \left[ \int_{t_{j-1}}^{t_{j}}\left( 1+\varepsilon ^{-1}e^{-\frac{\alpha s}{ 2\varepsilon }}\right) \mathrm{d}s\right] ^{2}+CN^{-2}, \end{eqnarray}
where we have used equations (19), (20) and (22), the assumptions on a,b, f,ϕ,K,L and the inequality
\int_{c_{1}}^{c_{2}}g(x)(x-c_{1})^{(k-1)}\mathrm{d}x\leq \frac{1}{k}\left\{ \int_{c_{1}}^{c_{2}}g(x)^{1/k}\mathrm{d}x\right\} ^{k}
for any positive monotonically decreasing function g on [c1,c2] and for arbitrary k ∈ N+; see [15].
For i ∈
i\in I_{1}^{2}
we have
\begin{eqnarray} \left\vert R_{i}\right\vert &\leq &C\int_{t_{i-1}}^{t_{i}}\varepsilon |u^{\prime \prime }(s)|\mathrm{d}s+C\sum_{j=1}^{i}h_{j}\int_{t_{j-1}}^{t_{j}} \left[ \left\vert \left( K_{i}(s)u(s)\right) ^{\prime \prime }\right\vert +\left\vert \left( L_{i}(s)u(s-r)\right) ^{\prime \prime }\right\vert \right] \left( s-t_{j-1}\right) \mathrm{d}s \notag \\ &\leq &C\int_{t_{i-1}}^{t_{i}}\left( \varepsilon +\varepsilon ^{-1}e^{-\frac{ \alpha s}{\varepsilon }}\right) \mathrm{d}s+C\sum_{j=1}^{i}h_{j} \int_{t_{j-1}}^{t_{j}}\left( 1+\varepsilon ^{-2}e^{-\frac{\alpha s}{ \varepsilon }}\right) \left( s-t_{j-1}\right) \mathrm{d}s+CN^{-2} \notag \\ &\leq &C\int_{t_{i-1}}^{t_{i}}\left( \varepsilon +\varepsilon ^{-1}e^{-\frac{ \alpha s}{\varepsilon }}\right) \mathrm{d}s+C\sum_{j=1}^{i}h_{j}\left[ \int_{t_{j-1}}^{t_{j}}\left( 1+\varepsilon ^{-1}e^{-\frac{\alpha s}{ 2\varepsilon }}\right) \mathrm{d}s\right] ^{2}+CN^{-2}, \end{eqnarray}
where we have used equations (21) and (22), the assumptions on a,b, f,ϕ,K,L and the inequality (24).
For i ∈
i\in I_{p}^{2}
with 2 ≤ p ≤ m we can also obtain
\begin{eqnarray} \left\vert R_{i}\right\vert &\leq &C\int_{t_{i-1}}^{t_{i}}\varepsilon |u^{\prime \prime }(s)|\mathrm{d}s+C\sum_{j=1}^{i}h_{j} \notag \\ &&\times \int_{t_{j-1}}^{t_{j}}\left[ \left\vert \left( K_{i}(s)u(s)\right) ^{\prime \prime }\right\vert +\left\vert \left( L_{i}(s)u(s-r)\right) ^{\prime \prime }\right\vert \right] \left( s-t_{j-1}\right) \mathrm{d}s \notag \\ &\leq &C\int_{t_{i-1}}^{t_{i}}\left( \varepsilon +\varepsilon ^{-1}e^{-\frac{ \alpha s}{\varepsilon }}\right) \mathrm{d}s+C\sum_{j=1}^{i}h_{j} \int_{t_{j-1}}^{t_{j}}\left( 1+\varepsilon ^{-2}e^{-\frac{\alpha s}{ \varepsilon }}+\varepsilon ^{-2}e^{-\frac{\alpha (s-r)}{\varepsilon } }\right) \left( s-t_{j-1}\right) \mathrm{d}s \notag \\ &\leq &C\int_{t_{i-1}}^{t_{i}}\left( \varepsilon +\varepsilon ^{-1}e^{-\frac{ \alpha s}{\varepsilon }}\right) \mathrm{d}s+C\sum_{j=1}^{i}h_{j}\left[ \int_{t_{j-1}}^{t_{j}}\left( 1+\varepsilon ^{-1}e^{-\frac{\alpha s}{ 2\varepsilon }}+\varepsilon ^{-1}e^{-\frac{\alpha (s-r)}{2\varepsilon } }\right) \mathrm{d}s\right] ^{2} \notag \\ &\leq &C\int_{t_{i-1}}^{t_{i}}\left( \varepsilon +\varepsilon ^{-1}e^{-\frac{ \alpha s}{\varepsilon }}\right) \mathrm{d}s \notag \\ &&+C\sum_{j=1}^{i}h_{j}\left[ \int_{t_{j-1}}^{t_{j}}\left( 1+\varepsilon ^{-1}e^{-\frac{\alpha s}{2\varepsilon }}\right) \mathrm{d} s+\int_{t_{j-1-N}}^{t_{j-N}}\left( 1+\varepsilon ^{-1}e^{-\frac{\alpha s}{ 2\varepsilon }}\right) \mathrm{d}s\right] ^{2}, \end{eqnarray}
where we also have used equations (21) and (22), the assumptions on a,b, f,K,L and the inequality (24).
For i ∈
i\in I_{1}^{1}
we have
\begin{eqnarray} \int_{t_{i-1}}^{t_{i}}\left( 1+\varepsilon ^{-1}e^{-\frac{\alpha s}{ 2\varepsilon }}\right) \mathrm{d}s=h_{i}-\frac{2}{\alpha }e^{-\frac{\alpha s }{2\varepsilon }}|_{t_{i-1}}^{t_{i}}=h_{i}-\frac{2}{\alpha }e^{-\frac{\alpha t_{i}}{2\varepsilon }}(1-e^{\frac{\alpha h_{i}}{2\varepsilon }}) && \notag \\ \leq C\varepsilon N^{-1}\ln N+C\frac{h_{i}}{\varepsilon }e^{-\frac{\alpha t_{i}}{2\varepsilon }}\leq CN^{-1}\ln N, && \label{eq4.11} \end{eqnarray}
where we have used mesh widths.
For i ∈
i\in I_{1}^{2}
we have
\begin{equation*} \int_{t_{i-1}}^{t_{i}}\left( \varepsilon +\varepsilon ^{-1}e^{-\frac{\alpha s }{\varepsilon }}\right) \mathrm{d}s\leq \varepsilon h_{i}+\alpha ^{-1}e^{- \frac{\alpha t_{i-1}}{\varepsilon }}\leq \varepsilon h_{i}+Ce^{-\frac{\alpha \sigma _{1}}{\varepsilon }}\leq C\varepsilon N^{-1}+CN^{-2}\leq CN^{-2}, \end{equation*}
where we also have used mesh widths and the assumption (11). Similarly, we can get
\begin{eqnarray} \int_{t_{i-1}}^{t_{i}}\left( 1+\varepsilon ^{-1}e^{-\frac{\alpha s}{ 4\varepsilon }}\right) \mathrm{d}s &=&h_{i}-4\alpha ^{-1}e^{-\frac{\alpha t_{i}}{4\varepsilon }}(1-e^{\frac{\alpha h_{i}}{4\varepsilon }})\leq CN^{-1}\ln N,\ \notag \\\ i &\in &I_{p}^{1},\ \ 2\leq p\leq m,\ \ \ \end{eqnarray}\begin{equation}\int_{t_{i-1}}^{t_{i}}\left( \varepsilon +\varepsilon ^{-1}e^{-\frac{\alpha s}{4\varepsilon }}\right) \mathrm{d}s\leq \varepsilon h_{i}+Ce^{-\frac{\alpha \sigma _{p}}{4\varepsilon }}\leq CN^{-2},\ \ i\in I_{p}^{2},\ 2\leq p\leq m.\end{equation}
Combining equations (23), (25), (26) with (27), we can obtain the desired result. Next we give a useful result for the error estimate.
Lemma 4
Assume that the mesh function\left\{ v_{i}\right\}_{0}^{N_{0}}is the solution of initial value probleml^{N}v_{i}=F_{i},\ \ \ \ \ \ \ \ \ \ \ \ 1\leq i\leq N_{0}, v_{0}=A with |Fi| ≤ ℱi and ℱi is a nondecreasing function, where\begin{equation} l^{N}v_{i}=\left\{ \begin{array}{ll} \varepsilon \frac{v_{i}-v_{i-1}}{h_{i}}+a_{i-1/2}\frac{v_{i-1}+v_{i}}{2} & \\ +\frac{1}{4}h_{i}\left[ K\left( t_{i-1/2},t_{i-1}\right) v_{i-1}+K\left( t_{i-1/2},t_{i-1/2}\right) \frac{v_{i-1}+v_{i}}{2}\right] , & i\in I_{p}^{1}, \\ \varepsilon \frac{v_{i}-v_{i-1}}{h_{i}}+a_{i}v_{i}+\frac{1}{2}h_{i}\left[ K\left( t_{i},t_{i-1}\right) v_{i-1}+K\left( t_{i},t_{i}\right) v_{i}\right] , & i\in I_{p}^{2}. \end{array} \right. \end{equation}Then the solution of problem (28)–(29) satisfies\left\vert v_{i}\right\vert \leq \left\vert A\right\vert +2\alpha ^{-1} \mathcal{F}_{i},\ \ \ \ \ \ \ \ \ \ \ \ 1\leq i\leq mNfor sufficiently large N.
Proof
It is easy to verify that the matrix associated with lN is an M-matrix for sufficiently large N. Hence the operator lN satisfies a discrete maximum principle. By applying the discrete maximum principle with the barrier function ψi = |A| + 2α−1 ℱi the desired result can be obtained.
Now we can obtain the following error estimates for the numerical scheme (12)–(13).
Theorem 5
Let u be the solution of problem (1)–(2) and U be the solution of finite difference scheme (12)–(13) on the Shishkin mesh ωN0. Then, under the assumptions ϕ ∈ C3(Ī0), a,b, f ∈ C3(Ī) and K,L ∈ C3(Ī×Ī) we have the following error estimates\left\Vert U-u\right\Vert _{\infty ,\omega _{N_{0}}}\leq CN^{-2}\ln ^{2}N.
Proof
The error equations (16)–(17) can be rewritten in the form
l^{N}z_{i}=F_{i},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 1\leq i\leq N_{0},z_{0}=0,
where lN is defined in (30) and
\begin{equation*} F_{i}=\left\{ \begin{array}{ll} R_{i}-\frac{1}{2}b_{i-1/2}\left( z_{i-1-N}+z_{i-N}\right) & \\ -{\sum_{j=1}^{i-1}\frac{h_{j}}{2}}\left[ K\left( t_{i-1/2},t_{j-1}\right) z_{j-1}+K\left( t_{i-1/2},t_{j}\right) z_{j}\right] & \\ -{\sum_{j=1}^{i-1}\frac{h_{j}}{2}\left[ L\left( t_{i-1/2},t_{j-1}\right) z_{j-1-N}+L\left( t_{i-1/2},t_{j}\right) z_{j-N}\right] } & \\ -\frac{h_{i}}{4}\left[ L(t_{i-1/2},t_{i-1})z_{i-1-N}+L(t_{i-1/2},t_{i-1/2}) \frac{z_{i-1-N}+z_{i-N}}{2}\right] , & \\ & i\in I_{p}^{1}, \\ {R_{i}-b_{i}z_{i-N}-\sum_{j=1}^{i-1}\frac{h_{j}}{2}\left[ K\left( t_{i},t_{j-1}\right) z_{j-1}+K\left( t_{i},t_{j}\right) z_{j}\right] } & \\ -{\sum_{j=1}^{i}\frac{h_{j}}{2}\left[ L\left( t_{i},t_{j-1}\right) z_{j-1-N}+L\left( t_{i},t_{j}\right) z_{j-N}\right] }, & i\in I_{p}^{2}. \end{array} \right. \end{equation*}
Then we have
\begin{equation*} \left\vert F_{i}\right\vert \leq \left\Vert R\right\Vert _{\infty ,\omega _{N,p}^{\ast }}+\left( \left\Vert b\right\Vert _{\infty ,\bar{I}}+T\bar{L} \right) \left\Vert z\right\Vert _{\infty ,\omega _{N,p-1}^{\ast }}+\bar{K} \sum_{j=1}^{i-1}h_{j}\left\vert z_{j}\right\vert . \end{equation*}
Applying Lemma 4 to (31)–(32) we have
\begin{equation*} \left\vert z_{i}\right\vert \leq \alpha ^{-1}\left[ \left\Vert R\right\Vert _{\infty ,\omega _{N,p}^{\ast }}+\left( \left\Vert b\right\Vert _{\infty , \bar{I}}+T\bar{L}\right) \left\Vert z\right\Vert _{\infty ,\omega _{N,p-1}^{\ast }}+\bar{K}\sum_{j=1}^{i-1}h_{j}\left\vert z_{j}\right\vert \right] . \end{equation*}
Thus, by using the discrete analogue of Grönwall’s inequality to the above inequality we can obtain
\begin{equation*}\left\Vert z\right\Vert _{\infty ,\omega _{N,p}^{\ast }}\leq \alpha ^{-1} \left[ \left\Vert R\right\Vert _{\infty ,\omega _{N,p}^{\ast }}+\left(\left\Vert b\right\Vert _{\infty ,\bar{I}}+T\bar{L}\right) \left\Vert z\right\Vert _{\infty ,\omega _{N,p-1}^{\ast }}\right] e^{\alpha ^{-1}\bar{K}T}.\end{equation*}
Therefore, by recursive calculation we have
\left\Vert z\right\Vert _{\infty ,\omega _{N,p}^{\ast }}\leq C\sum_{k=1}^{p}\left\Vert R\right\Vert _{\infty ,\omega _{N,p}^{\ast }}.
Combining (33) with Lemma 3 we can get the desired result.
Numerical results
In this section, we experimentally verify the theoretical results obtained in the preceding section. The error estimates and convergence rates for the numerical scheme are presented for an example which has been presented in [7].
Error estimates and convergence rates for various N for Example
We define the maximum pointwise error
e_{\varepsilon }^{N}
and the computed parameter-uniform maximum pointwise error eN as follows:
e_{\varepsilon }^{N}=\left\Vert {U-u}\right\Vert _{\infty ,\omega _{N}},\ \ \ \ \ \ \ \ \ \ \ e^{N}=\mathop {\max }\limits_{\varepsilon }e_{\varepsilon}^{N},
where U is the numerical approximation to the exact solution u for various values of N and ɛ. We also define the “Shishkin” convergence rate and the computed parameter-uniform “Shishkin” convergence rate as follows:
\begin{eqnarray*}r_{\varepsilon }^{N} &=&\frac{\ln e_{\varepsilon }^{N}-\ln e_{\varepsilon}^{2N}}{\ln \left( 2\ln N/\ln \left( 2N\right) \right) },\ \ \\ \ \ \ \ \ \ r^{N} &=&\frac{\ln e^{N}-\ln e^{2N}}{\ln \left( 2\ln N/\ln\left( 2N\right) \right) }.\end{eqnarray*}
The numerical results for Example 6 are presented in Table 1. From Table 1, we confirm that our method proposed in this paper is more accurate than those of [7]. It is visible that the results seem to be ɛ -uniform as expected and rN are 2 for sufficiently large N independent of ɛ, which is fully in agreement with our theoretical analysis.
Conclusions
In this study, we have focused on the numerical treatment of singularly perturbed Volterra integro-differential equations with delay. Our main objective was to develop a second-order accurate numerical method for efficiently approximating the solutions of these equations. The novelty of our approach lies in the combination of a hybrid difference scheme, quadrature rules, and a Shishkin-type mesh for solving singularly perturbed Volterra integro-differential equations with delay. Previous studies have primarily focused on either singularly perturbed differential equations or Volterra integro-differential equations, but our work bridges the gap between these two areas. Furthermore, the incorporation of delay in the equations adds an additional layer of complexity, which has not been extensively explored in the existing literature. The numerical experiments conducted in this study provide strong support for the theoretical results obtained. The experiments demonstrate that our proposed method accurately approximates the solutions of singularly perturbed Volterra integro-differential equations with delay. Moreover, the estimates obtained from the experiments highlight the sharpness of our theoretical results. In conclusion, this research contributes to the field by presenting a novel second-order numerical method for singularly perturbed Volterra integro-differential equations with delay. The hybrid difference scheme, combined with appropriate quadrature rules and a Shishkin-type mesh, provides an efficient and accurate approach for solving these challenging equations. Theoretical analysis and numerical experiments validate the effectiveness and sharpness of the proposed method. This work opens up avenues for further research in the numerical treatment of singularly perturbed Volterra integro-differential equations with delay, and it has the potential to find applications in various fields, such as mathematical biology and engineering.
Declarations
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Funding
The authors hereby declare that there is no any funding interest.
Author’s contribution
F.E.-Conceptualization, Data Curation, Writing-Review Editing, Software, Writing-Original Draft, Supervisor, Methodology. All authors read and approved the final submitted version of this manuscript.
Acknowledgement
The authors sincerely appreciate the reviewers’ insightful remarks, which resulted in a significant improvement in the manuscript.
Data availability statement
All data that support the findings of this study are included within the article.
Using of AI tools
The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.
