# A Mechanical Quadrature Method for Solving Delay Volterra Integral Equation with Weakly Singular Kernels

Authors: Li Zhang; Jin Huang; Yubin Pan; Xiaoxia Wen
Journal: Complexity
 (2019)

---

## Abstract

In this work, a mechanical quadrature method based on modified trapezoid formula is used for solving weakly singular Volterra integral equation with proportional delays. An improved Gronwall inequality is testified and adopted to prove the existence and uniqueness of the solution of the original equation. Then, we study the convergence and the error estimation of the mechanical quadrature method. Moreover, Richardson extrapolation based on the asymptotic expansion of error not only possesses a high accuracy but also has the posterior error estimate which can be used to design self-adaptive algorithm. Numerical experiments demonstrate the efficiency and applicability of the proposed method.


---

## Body

## 1. Introduction

In recent years, Volterra integral equation with delay has received a considerable amount of attention. This paper considers the following weakly singular Volterra integral equations with proportional delays:(1)ut=gt+Iut+Iθut,t∈I=0,T,with(2)Iut≔∫0tK1t,susds,Iθut≔∫0θtK2t,susds,and(3)θt≔qt,0<q<1,K1t,s≔sαk1t,s,K2t,s≔sγk2t,s,-1<α<0,-1<γ≤0,where k1(t,s) and k2(t,s) are known continuous functions defined on the domains D≔{(t,s):0≤s≤t≤T} and Dθ≔{(t,s):0≤s≤θ(t),t∈I=[0,T]}, respectively. g(t) is a known function and u(t) is an unknown function. In practice, the delay arguments are consistent with the real phenomena which make the models more realistic for simulation. Delay integral equation and partial differential equation have been widely used in many population growth and relevant phenomena in mathematical biology [1–3].Numerical algorithms for implementing delay models should be designed specially according to the nature of the equations. There are many numerical techniques for the delay differential equations [4–8] and integral-differential equations [9–11]. Simultaneously, vast researchers also focused their interests on the numerical techniques of delay integral equations with continuous kernels, such as least squares approximation method [12], spectral method [13], Bernoulli wavelet method [14], and collocation method [15]. Xie et al. in [16] handled the Volterra integral equation with the delay function vanishing at the initial point in the given interval; they found that the iterated collocation solution possessed local superconvergence at the mesh points. In [17], the authors adopted multistep method based on Hermite collocation method and it turned out that the numerical method had uniform order 2m+2r with m collocation points and r previous time steps. In [18], an hp-spectral collocation method was used for nonlinear Volterra integral equations with vanishing variable delays.There are a few researches on the delay integral equation with weakly singular kernels, such as [19]. In this paper, we concentrate on (1) whose upper limit of integral is a delay function and the integral kernels are weakly singular functions, which increase the computational complexity and the theoretical difficulty. To the best of our knowledge, there are no studies in (1) by the mechanical quadrature method in recent years. A further advantage of this method for (1) is that the error has an asymptotic expansion. Thus, we can improve the accuracy order of approximation. Simultaneously, the theoretical analysis is complete and the calculation is simplified.In this paper, we firstly deduce the improved Gronwall inequality; the existence and uniqueness of the solution for (1) are testified via the improved Gronwall inequality. Then the equation is approximated by employing the floor technique to the delay argument θ(t) and by adopting the quadrature formula [20] to the weakly integrals. Next, the approximate equation for (1) is constructed by combining the mechanical quadrature method and the interpolation technique; then, we use the iterative method for solving the approximate equation. The existence and uniqueness of the solution for the approximate equation are testified via the discrete Gronwall inequality. Finally, we prove that the convergent order is O(h2+min(α,γ)). In order to achieve a higher accuracy order O(h2), the Richardson-h2+min(α,γ) extrapolation based on the asymptotic expansion of error is adopted; moreover, a posterior error estimate is realized conveniently.The layouts of this paper are as follows. In Section2, we prove the existence and uniqueness of the solution for (1). In Section 3, we introduce the quadrature method, iterative method, and interpolation technique. In Section 4, the existence and uniqueness of the solution for the approximate equation are discussed. In Section 5, the convergence and the error estimation are obtained to ensure the reliability of the method. In Section 6, the asymptotic expansion of the error is achieved, a higher accuracy order is realized by extrapolation, and a posterior error estimate is derived. In Section 7, some numerical examples are demonstrated to illustrate the theoretical results. Some concluding remarks are provided in Section 8.


## 2. The Existence and Uniqueness of the Solution for the Original Equation

In this section, we will verify the existence and uniqueness of the solution for (1). We first prove the improved Gronwall inequality.Lemma 1.
Suppose thatg(t), h(t), and u(t) are nonnegative integrable functions, t∈[0,T], 0<q<1, and A≥0; based on the inequality (4)ut≤A+∫0tgsusds+∫0qthsusds,we have(5)ut≤Ae∫0tgs+hsds.Proof.
Due to the fact that0<q<1, qt<t and(6)ut≤A+∫0tgsusds+∫0thsusds≤A+∫0tgs+hsusds.Let H(s)=g(s)+h(s); we can deduce(7)HtutA+∫0tHsusds≤Ht.Integrate on both sides (8)ln⁡A+∫0tHsusds-ln⁡A≤∫0tHsds,and then (9)A+∫0tHsusds≤Ae∫0tHsds;namely, (10)ut≤Ae∫0tHsds.The proof of the Lemma 1 is completed.Theorem 2.
Assume thatki(t,s)(i=1,2) are known continuous functions defined on the domains D≔{(t,s):0≤s≤t≤T} and Dθ≔(t,s):0≤s≤θ(t),t∈I=[0,T], respectively; then the solution of (1) is existent uniquely.Proof.
We construct the sequenceuk(t)(k∈N), satisfying (11)u0t=u0=gt,ukt=gt+∫0tK1t,suk-1sds+∫0qtK2t,suk-1sds,where K1(t,s) and K2(t,s) are defined in (3), with k∈N. ki(t,s)(i=1,2) are continuous functions; then there exists a constant C such that kit,s≤C. We have(12)u2t-u1t≤C∫0tsαu1s-u0sds+C∫0qtsγu1s-u0sds≤C∫0tsαu1s-u0s+sγu1s-u0sds≤Ca∫0tsα+sγds≤2Ca∫0tsμds=2Catμ+1μ+1,with a=max0≤t≤Tu1t-u0t; μ=min⁡α,γ. Now, we can deduce(13)ukt-uk-1t≤a2Ck-1k-1!μ+1k-1tk-1μ+1.By means of the mathematical induction, when n=k+1, we obtain(14)uk+1t-ukt≤C∫0tsαuks-uk-1sds+C∫0qtsγuks-uk-1sds≤C∫0tsαuks-uk-1s+sγuks-uk-1sds≤2C∫0tsμuks-uk-1sds≤a2Ckk!μ+1ktkμ+1.
Next, we prove thatuk(t) is the basic sequence in C[0,T]; in fact, (15)ukt-uk+mt≤uk+1t-ukt+uk+2t-uk+1t+⋯+uk+mt-uk+m-1t≤a2Ckk!μ+1ktkμ+1+⋯+a2Ck+m-1k+m-1!μ+1k+m-1tk+m-1μ+1.For sufficiently small ϵ>0, there exists a positive integer N such that when n>N and any m>0, we have(16)unt-un+mt≤ϵ.According to Cauchy’s test for convergence, the sequence un(t)(n∈N) is convergent uniformly to u(t) which is the solution of (1).
Suppose that bothu and v are the solutions of (1); let w=u-v; then we get(17)wt=∫0tsαk1t,swsds+∫0qtsγk2t,swsds≤∫0tsαk1t,swsds+∫0qtsγk2t,swsds≤∫0tsαk1t,sws+sγk2t,swsds≤∫0tsαk1t,s+sγk2t,swsds.Next, we verify that sαk1t,s+sγk2t,s is integrable;(18)∫0tsαk1t,s+sγk2t,sds≤C∫0tsα+sγds≤Ct1+α1+α+t1+γ1+γ. According to Lemma 1, we can derive that wt=0; the solution is unique.


## 3. The Quadrature Method and the Iterative Algorithm

Let the delay functionθ(t)=t-τ(t)≥0, t∈[0,T], satisfy the following conditions [21]:(1)τ(0)=0, and τ(t)>0 with t∈(0,T] (vanishing delay);(2)θ(t)≤q1t on I=[0,T] for some q1∈(0,1) and θ′(t)≥q0>0; qi(i=0,1) are constants;(3)θ(t)∈C1(I).For0<q<1, the special case is τ(t)=(1-q)t; we get θ(t)=t-τ(t)=qt; (1) is a weakly integral equation with vanishing delay. In order to solve (1), the modified trapezoid quadrature formula is considered to deal with the integrals (Iu)(t) and (Iθu)(t).It is challenging because the integral upper limit is a delay functionθ(t) and the integral kernels are weakly singular at the endpoint.Lemma 3 (see [22]).
LetI(G)=∫abG(x)dx=∫abx-aαg(x)dx, -1<α<0, and g(x)∈C2m[a,b]; then the modified trapezoidal rule is(19)QNG=-h1+αgaζ-α+h∑i=1N-1Gxi+h2Gb,where ζ is the zeta function. Further, I(G) has the following asymptotic expansion:(20)ENG=QNG-IG=∑j=1m-1B2j2j!G2j-1bh2j+∑j=12m-1-1jζ-α-jgjahj+α+1j!+Oh2m,where B(2j),j=1,…,m-1, are the Bernoulli numbers.Lemma 4.
Letu∈C3[0,1], and z=βx+(1-β)y with β∈[0,1],x,y∈[0,T]; then there is(21)uz=βux+1-βuy-β1-β2x-y2u′′z+Ox-y3.Proof.
The Taylor expansion of functionu(x) at the point z is(22)ux=uβx+1-βx=uβx+1-βy+1-βx-y=uz+1-βx-y=uz+1-βx-yu′z+1-β22x-y2u′′z+Ox-y3.Similarly, the Taylor expansion of function u(y) at the point z is(23)uy=uz-βx-y=uz-βx-yu′z+β22x-y2u′′z+Ox-y3.Combining (22) with (23), we derive (21).Now, the concrete approximate scheme of (1) will be constructed on the basis of the above lemmas. We divide I∈[0,T] into N subintervals with the uniform step size h=T/N; ti=ih,i=0,…,N. Let t be replaced by ti in (1); that is to say,(24)uti=gti+∫0tisαk1ti,susds+∫0qtisγk2ti,susds=gti+∫0tisαk1ti,susds+∫0tqisγk2ti,susds+∫tqiqtisγk2ti,susds=gti+I1+I2+I3,where [qi] denotes the maximum integer less than qi. It is obvious that u(qti), 0<q<1, are not the node values, and we can turn u(qti) to the adjacent points u(t[qi]) and u(t[qi]+1) by the linear interpolation approximation with t[qi]≤qti≤t[qi]+1. Then, we have that qti=βit[qi]+(1-βi)t[qi]+1 are established for βi∈[0,1] with βi=1+[qi]-qi, and(25)uqti=βiutqi+1-βiutqi+1.Now, we structure the quadrature algorithm; by Lemma3 and the trapezoidal quadrature rule, we obtain(26)I1=∫0tisαk1ti,susds=h∑k=1i-1tkαk1ti,tkutk+h2tiαk1ti,tiuti-ζ-αh1+αk1ti,t0ut0+E1,isαk1t,s,(27)I2=∫0tqisγk2ti,susds=h∑k=1qi-1tkγk2ti,tkutk+h2tqiγk2ti,tqiutqi-ζ-γh1+γk2ti,t0ut0+E2,isγk2t,s,(28)I3=∫tqiqtisγk2ti,susds=qti-tqi2tqiγk2ti,tqiutqi+qtiγk2ti,qtiuqti+E3,i(sγk2t,s=qti-tqi2tqiγk2ti,tqiutqi+qtiγk2ti,qtiβiutqi+1-βiutqi+1+E3,isγk2t,s,where E1,i, E2,i, and E3,i are error functions, which are, respectively, equal to(29)E1,i=k1ti,sus′s=0ζ-α-1h2+α+k1ti,sus′s=02!ζ-α-2h3+α+oh4+α=T1tih2+α+Oh3+α,E2,i=k2ti,sus′s=0ζ-γ-1h2+γ+k2ti,sus′s=02!ζ-γ-2h3+γ+oh4+γ=T2tih2+γ+Oh3+γ,E3,i=-β1-β2h2u′′qtiqtiγk2ti,qtiβiutqi+1-βiutqi+1+qti-tqi212∫tqiqti∂2∂s2k2ti,sussγds+Oh3=T3tih2+qti-tqi2-h212∫tqiqti∂2∂s2k2ti,sussγds+Oh3=T3tih2+Oh3,where (30)T1ti=k1ti,sus′s=0ζ-α-1,T2ti=k2ti,sus′s=0ζ-γ-1,T3ti=-β1-β2u′′qtiqtiγk2ti,qtiβiutqi+1-βiutqi+1+112∫tqiqti∂2∂s2k2ti,sussγds.The discrete forms of (1) are obtained:(31)ut0=gt0;(32)uti=gti+h∑k=1i-1tkαk1ti,tkutk+h2tiαk1ti,tiuti-ζ-αh1+αk1ti,t0ut0+h∑k=1qi-1tkγk2ti,tkutk+h2tqiγk2ti,tqiutqi-ζ-γh1+γk2ti,t0ut0+qti-tqi2tqiγk2ti,tqiutqi+qtiγk2ti,qtiβiutqi+1-βiutqi+1+E1,i+E2,i+E3,i,i=1,…,N.Let ui be the approximate solution of u(ti) and ignore the error function; (32) becomes(33)u0=gt0;ui=gti+h∑k=1i-1tkαk1ti,tkuk+h2tiαk1ti,tiui-ζ-αh1+αk1ti,t0u0+h∑k=1qi-1tkγk2ti,tkuk+h2tqiγk2ti,tqiuqi-ζ-γh1+γk2ti,t0u0+qti-tqi2tqiγk2ti,tqiuqi+qtiγk2ti.qti(βiuqi+1-βiuqi+1.The iterative algorithm is built to solve (33).Iterative AlgorithmStep 1.
Take sufficiently smallϵ>0 and set u0=g(t0); i≔1;Step 2.
Letui0=ui-1,m≔0; then compute uim+1(i≤N) as follows:(34)uim+1=gti+h∑k=1i-1tkαk1ti,tkukm+h2tiαk1ti,tiuim-ζ-αh1+αk1ti,t0u0+h∑k=1qi-1tkγk2ti,tkukm+h2tqiγk2ti,tqiuqim-ζ-γh1+γk2ti,t0u0+qti-tqi2tqiγk2ti,tqiuqim+qtiγk2ti,qti(βiuqim+1-βiuqi+1m;Step 3.
Ifuim+1-uim≤ϵ, set ui≔uim+1 and i=i+1 and return to Step 2; else, let m=m+1 and return to Step 2.


## 4. The Existence and Uniqueness of the Solution for the Approximate Equation

Now, we prove the existence and uniqueness of the solution for the approximate equation. We first introduce the following lemma.Lemma 5 (see [23]).
Suppose that the sequence{wn}, n=0,…,N, satisfies(35)wn≤h∑k=1nBkwk+A,0≤n≤N,where A and Bk(k=1,…,N) are nonnegative constants. Let hmax1≤k≤N⁡Bk≤1/2 with h=1/N; then we can derive (36)max0≤n≤Nwn≤Aexp⁡2h∑k=1NBk.Theorem 6.
Assume thath is sufficiently small; then the solution of (34) is existent uniquely, and the algorithm converges at a geometrical rate.Proof.
From the nature of the delay functionθ(t), we discuss the existence and uniqueness of the solution for the approximate equation under two situations.
First, we prove that the solution of (34) is existent under two situations.
(1) One situation is[qi]+1=i; that is, when i<1/(1-q), we can easily obtain (37)uim+1-uim≤h2tiαk1ti,tiuim-uim-1+qti-tqi2qtiγk2ti,qti1-βiuqi+1m-uqi+1m-1≤Ch2tiα+qtiγ1-βiuim-uim-1.Let Ch/2(tiα+(qti)γ(1-βi))≤1/2 for a sufficiently small h; then uim+1-uim≤1/2uim-uim-1 holds.
(2) The other situation is[qi]+1<i; that is, when i>1/(1-q), then(38)uim+1-uim≤h2tiαk1ti,tiuim-uim-1≤Ch2tiαuim-uim-1.Let Ch/2tiα≤1/2 for a sufficiently small h; then uim+1-uim≤1/2uim-uim-1 holds.
With the discussion of the above two situations, one can conclude that the iterative algorithm is convergent geometrically, and the limit is the solution of (34); therefore, the solution of (34) is existent.
Next, we prove that the solution of (34) is unique. If ui and vi are solutions of (34), the difference can be represented as zi=ui-vi, 1≤i≤N, and(39)z0=0;zi=∑j=1i-1htjαk1ti,tjzj+h2tiαk1ti,tizi-ζ-αh1+αk1ti,t0z0+∑j=1qi-1htjγk2ti,tjzj+h2tqiγk2ti,tqizqi-ζ-γh1+γk2ti,t0z0+qti-tqi2tqiγk2ti,tqizqi+qtiγk2ti,qti(βizqi+1-βizqi+1,and kit,s≤C(i=1,2), because ki(t,s) are continuous on bounded domains, with qti-t[qi]≤h; then(40)z0=0;zi≤C∑j=1i-1htjαzj+Ch2tiαzi+∑j=1qi-1htjγCzj+h2tqiγCzqi+Cqti-tqi2tqiγzqi+qtiγ(βizqi+1-βizqi+1,≤C∑j=1i-1htjαzj+Ch2tiαzi+∑j=1qi-1htjγCzj+h2tqiγCzqi+Ch2tqiγzqi+tqiγ(βizqi+1-βizqi+1,with ti=ih,i=1,…,N. Let μ=min{α,γ}.
(1) The first situation is[qi]+1=i, that is, when i<1/(1-q). Let hBj≤1/2 for a sufficiently small h; then (40) can be written as follows:(41)zi≤C∑j=1qi-1htjμzj+Chtqiμzqi+Ch2tiμzi+C∑j=1qi-1htjμzj+Ch2tqiμzqi+Ch2tqiμzqi+tqiμβizqi+1-βizqi+1=2C∑j=1qi-1htjμzj+2Chtqiμ+Ch2tqiμβizqi+Ch2tiμ+Ch2tqiμ1-βizqi+1;namely,(42)zi≤∑j=1iBjzj=2C∑j=1qi-1htjμzj+2Chtqiμ+Ch2tqiμβizqi+Ch2tiμ+Ch2tqiμ1-βizqi+1,where(43)Bj=2C∑j=1qi-1tjμ,j=1,…,qi-1,2Ctqiμ+C12tqiμβi,j=qi,C12tiμ+C12tqiμ1-βi,j=i.
(2) The second situation is[qi]+1<i, that is, when i>1/(1-q). Let hBj≤1/2 for a sufficiently small h; then (40) can be written as(44)zi≤C∑j=1qi-1htjμzj+Chtqiμzqi+Chtqi+1μzqi+1+∑j=qi+1i-1htjμzj+Ch2tiμzi+C∑j=1qi-1htjμzj+Ch2tqiμzqi+Ch2tqiμzqi+tqiμ(βizqi+1-βizqi+1=2Ch∑j=1qi-1tjμzj+2Chtqiμ+Ch2tqiμβizqi+Chtqi+1μ+Ch2tqiμ1-βizqi+1)+∑j=qi+1i-1htjμzj+Ch2tiμzi;namely,(45)zi≤∑j=1iBjzj=2Ch∑j=1qi-1tjμzj+2Chtqiμ+Ch2tqiμβizqi+Chtqi+1μ+Ch2tqiμ1-βizqi+1+∑j=qi+1i-1htjμzj+Ch2tiμzi.And(46)Bj=2C∑j=1qi-1tjμ,j=1,…,qi-1,2Ctqiμ+C12tqiμβi,j=qi,Ctqi+1μ+C12tqiμ1-βi,j=qi+1,∑j=qi+1i-1tjμ,j=qi+1,…,i-1,C12tiμ,j=i.Based on Lemma 5, zi=0 with A=0 and the solution of the discrete equation (34) is unique. The proof of Theorem 6 is completed.


## 5. The Error Estimation

In this section, we give the error estimate between the approximation solution and the exact solution of (1).Theorem 7.
Letu(t) be the exact solution of (1); the kernel functions K1(t,s)=sαk1(t,s), K2(t,s)=sγk2(t,s), -1<α,γ≤0, and the functions k1(t,s) and k2(t,s) are continuous in the domains D≔{(t,s):0≤s≤t≤T} and Dθ≔{(t,s):0≤s≤θ(t),t∈I}, respectively. Then there is a positive constant b independent of h such that kit,s≤C and ei=uti-ui(i=0,…,N) have the following estimation: (47)max1≤i≤N⁡ei≤bh2+μ.Proof.
From (26) and (27), we have(48)E1,i=k1ti,sus′s=01!ζ-α-1h2+α+k1ti,sus′s=02!ζ-α-2h3+α+oh4+α=T1tih2+α+Oh3+α,E2,i=k2ti,sus′s=01!ζ-γ-1h2+γ+k2ti,sus′s=02!ζ-γ-2h3+γ+oh4+γ=T2tih2+γ+Oh3+γ.By the trapezoidal formula and Lemma 4, we have(49)E3,i=-β1-β2h2u′′qtiqtiγk2ti,qtiβiutqi+1-βiutqi+1+qti-tqi212∫tqiqti∂2∂s2k2ti,sussγds+Oh3=T3tih2+qti-tqi2-h212∫tqiqti∂2∂s2k2ti,sussγds+Oh3=T3tih2+Oh3,where(50)T3ti=-β1-β2u′′qtiqtiγk2ti,qtiβiutqi+1-βiutqi+1+112∫tqiqti∂2∂s2k2ti,sussγds.Then, we have(51)uti=gti+QN1+QN2+QN3+T1tih2+α+T2tih2+γ+T3tih2+Oh2,QN1=∑j=1i-1htjαk1ti,tjuj+h2tiαk1ti,tiui-ζ-αh1+αu0kti,t0,QN2=∑j=1qi-1htjγk2ti,tjuj+h2tqiγk2ti,tqiuqi-ζ-γh1+γu0k2ti,t0,QN3=qti-tqi2tqiγk2ti,tqiuqi+qtiγk2ti,qti(βiuqi+1-βiuqi+1, where βi=1+[qi]-qi. Subtracting (34) from (51), we get(52)e0=0,ei=∑j=1i-1htjαk1ti,tjej+h2tiαk1ti,tiei+∑j=1qi-1htjγk2ti,tjej+h2tqiγk2ti,tqieqi+qti-tqi2tqiγk2ti,tqieqi+qtiγk2ti,qti(βieqi+1-βieqi+1+T1tih2+α+T2tih2+β+T3tih2+Oh3+μ≤C∑j=1i-1htjαej+Ch2tiαei+∑j=1qi-1htjγCej+h2tqiγCeqi+Cqti-tqi2tqiγeqi+qtiγ(βieqi+1-βieqi+1+T1tih2+α+T2tih2+γ+T3tih2+Oh3+μ.Let(53)A=T1tih2+α+T2tih2+γ+T3tih2+Oh3+μ=Oh2+μ,and the analysis is the same as (41) and (44); we have(54)ei≤A+∑j=1iBjej,1≤i≤N.From Lemma 5, there is a positive constant b independent of h such that(55)max1≤i≤Nei≤bh2+μ.The proof of Theorem 7 is completed.


## 6. Error Asymptotic Expansion and Extrapolation Algorithm

In this section, we present the main theoretical result of the error asymptotic expansions and the relevant extrapolation algorithm.Theorem 8.
Based on the conditions of Theorem7, there exist continue functions T^i(t),(i=1,2,3) satisfying the asymptotic expansion(56)ui=uti+T^1tih2+α+T^2tih2+γ+T^3tih2+Oh3+μ,-1<α<0,-1<γ≤0.Proof.
Suppose that{T^k(t),k=1,2,3} satisfies the auxiliary delay equations: (57)T^kt=Tkt+∫0tsαk1t,susT^ksds+∫0qtsγk2t,susT^ksds,and T^k(ti), i=1,…,N, satisfy the approximation equations:(58)T^kti=Tkti+∑j=1i-1htjαk1ti,tjujT^ktj+h2tiαk1ti,tiuiT^kti-ζ-αh1+αu0kti,t0T^kt0+∑j=1qi-1htjγk2ti,tjujT^ktj+T^ktjh2tqiγk2ti,tqiuqiT^ktqi-ζ-γh1+γu0k2ti,t0T^kt0+qti-tqi2tqiγk2ti,tqiuqiT^qi+qtiγk2ti,qtiβuqtiT^qi+1-βuqi+1T^qi+1.Similar to the proof of Theorem 7, we obtain(59)max1≤i≤N⁡T^kti-Tti≤Ch2+μ.Let(60)Ei=ei-T1tih2+α+T2tih2+γ+T3tih2.We have(61)Ei=∑j=1i-1htjαk1ti,tjujEj+h2tiαk1ti,tiuiEi-ζ-αh1+αu0kti,t0E0+∑j=1qi-1htjγk2ti,tjujEj+h2tqiγk2ti,tqiuqiEqi-ζ-γh1+γu0k2ti,t0E0+qti-tqi2tqiγk2ti,tqiuqiT^qi+qtiγk2ti,qtiβuqti+1-βuqi+1βEqi+1-βEqi+1.From Lemma 1, there exists a constant d such that(62)max1≤i≤N⁡Ei≤dh3+μ.The asymptotic expansion is (56).Based on Theorem8, we adopt the Richardson extrapolation to improve the accuracy.Extrapolation AlgorithmStep 1.
Assume thatμ=min(α,γ)=α, and halve the step length to obtain(63)uih/2=uti+T^1tih22+α+T^2ti)h22+γ+T^3tih22+Oh23+α.By combining (56) with (63), we get(64)ui1,h=22+αuih/2-uih22+α-1=uti+T^2ti)h2+γ+T^3tih2+Oh3+α.Step 2.
We implement Richardsonh2+γ extrapolation:(65)ui1,h/2=uti+T^2ti)h22+γ+T^3tih22+Oh23+α.Combining (64) with (65), we get(66)ui2,h=22+γui1,h/2-ui1,h22+γ-1=uti+T^3tih2+Oh3+α.Moreover, a posterior asymptotic error estimate(67)uih/2-ui=22+αuih/2-uih22+α-1-ui+uih-uih/222+α-1≤22+αuih/2-uih22+α-1-ui+uih-uih/222+α-1=22+αuih/2-uih22+α-1-ui+uih-uih/222+α-1+Oh2.The error uih/2-u(i) is bounded by (uih-uih/2)/22+α-1, which is essential to construct adaptable algorithms.


## 7. Numerical Experiments

In this section, three examples will be presented to show the efficiency of the quadrature method. We design a set of grids on the intervalI; the absolute error is denoted by(68)ehti=uti-uiwith i=1,…,N; u(ti) and ui are the exact solution and the approximate solution at t=ti, respectively. Set h=1/N, and the convergence order is defined by(69)Rate=log⁡eh/eh/2log⁡2.Example 1.
Consider the following equation:(70)ut=gt-∫0tt-sαusds,t∈0,T,where T=1,-1<α<0, and the analytical solution is u(t)=cos⁡(t).
The numerical results at the pointt=1 with the partitions N=24,25,26,27,28 are addressed in Table 1. By the mechanical quadrature method and iterative method, the obtained absolute errors eh with more refined partitions show a more accurate approximate solution, and the convergent rate is adjacent to 3/2-order. Based on the Richardson h2+α extrapolation, the errors are closer to the exact solution, and the convergence order is improved to 2-order. The posteriori error is also achieved. The two kinds of convergence orders are consistent with the theoretical analysis.Table 1
Numerical results at pointt=1 of Example 1.
N e h Rate R-h2+α extrapolation Rate Posteriori errors 24 9.27e-04 — — — — 25 3.22e-04 2 1.53 9.65e-06 — 3.31e-04 26 1.12e-04 2 1.52 2.24e-06 2 2.10 1.14e-04 27 3.93e-05 2 1.51 5.32e-07 2 2.08 3.99e-05 28 1.38e-05 2 1.51 1.28e-07 2 2.06 1.39e-05Example 2.
Consider the following delay Volterra integral equation with weakly singular kernel as(71)ut=gt-∫0ts-1/2usds+∫0qts-1/2usds,t∈0,T,with q=0.9, T=1, and the initial value u(0)=g(0)=0, and g(t) is determined by the analytical solution u(t)=t+1.
We denote the approximate solution byuh with h=1/32. The error results at some interior points in the interval t∈[0,1] with different partition are listed in Table 2. The comparison of the exact solution and the approximate solution with partition N=23 is shown in Figure 1. It is obvious that this paper provides a high accuracy algorithm for the weakly singular Volterra integral equation with proportional delays.Table 2
Numerical results of Example2.
t e h e h / 2 e h / 4 Posteriori error 0.2 3.89e-06 7.84e-07 1.77e-07 1.69e-06 0.4 5.05e-06 1.01e-06 2.28e-07 2.20e-06 0.6 5.67e-06 1.14e-06 2.55e-07 2.48e-06 0.8 5.98e-06 1.19e-06 2.68e-07 2.61e-06 1.0 6.07e-06 1.21e-06 2.71e-07 2.66e-06Figure 1
Comparison of exact solution and approximate solution withN=23.Example 3.
Consider the following delay Volterra integral equation with weakly singular kernel as(72)ut=gt-∫0ts-1/4usds+∫0qtt2+susds,t∈0,T,with q=0.9,T=1, and the initial value u(0)=g(0)=0, and g(t) is determined by the analytical solution u(t)=t+1.
The absolute erroreht=ut-uht,t∈[0,1], and the posteriori errors at some interior points are listed in Table 3 when N=25,26,27, which can observe that the error results decay quickly with the increasing of N. The absolute errors between the exact solution and the approximate solution with partition N=22,23 are shown in Figure 2, which indicate that our algorithm is effective.Table 3
Numerical results of Example3.
t e h e h / 2 e h / 4 Posteriori errors 0.2 1.43e-06 9.29e-07 7.18e-07 2.75e-07 0.4 1.23e-05 4.52e-07 9.09e-07 6.48e-06 0.6 4.79e-05 8.27e-06 5.23e-07 2.17e-05 0.8 1.34e-04 2.90e-05 5.10e-06 5.74e-05 1.0 3.25e-04 7.63e-05 1.16e-05 1.36e-04Figure 2
Comparison of exact solution and approximate solution withN=23.


## 8. Conclusion

In this paper, we use the mechanical quadrature method and Romberg extrapolation for weakly Volterra integral equation with proportional delays. Most papers analyze the delay Volterra integral equation with continuous kernels; the study for the delay Volterra integral equation with weakly singular kernels still faces a real challenge for the relevant researchers in both numerical computation and theoretical analysis. The improved Gronwall inequality is adopted to prove the existence and uniqueness of the solution of the original equation. At the same time, the discrete Gronwall equality and iterative method are adopted to prove the existence and uniqueness of the solution of the approximate equation. Moreover, according to the error asymptotic expansion of the mechanical quadrature method and the extrapolation, a high order of accuracy can be achieved and a posterior error estimation can be obtained. Both the theoretical analysis and the numerical examples show that the presented method is efficient.



---

*Source: 4813802-2019-06-16.xml*