# A Class of New Pouzet-Runge-Kutta-Type Methods for Nonlinear Functional Integro-Differential Equations

Authors: Chengjian Zhang
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

This paper presents a class of new numerical methods for nonlinear functional-integrodifferential equations, which are derived by an adaptation of Pouzet-Runge-Kutta methods originally introduced for standard Volterra integrodifferential equations. Based on the nonclassical Lipschitz condition, analytical and numerical stability is studied and some novel stability criteria are obtained. Numerical experiments further illustrate the theoretical results and the effectiveness of the methods. In the end, a comparison between the presented methods and the existed related methods is given.


---

## Body

## 1. Introduction

In the last ten years the numerical analysis and computational solution of various types of functional-integrodifferential equations (FIDEs) have received considerable attention. Many of these numerical schemes were derived by suitably adapting classical numerical methods for ordinary differential equations (ODEs) or integrodifferential equations (IDEs) to FIDEs, and there is a growing literature on their convergence and stability properties. Of these papers, Zhang and Vandewalle [1, 2] dealt with nonlinear numerical stability of FIDEs of the form(1.1)x′(t)=F(t,x(t),x(t-τ),∫t-τtg(t,ξ,x(ξ))dξ),t≥t0,x(t)=φ(t),t0-τ≤t≤t0.
Yu et al. [3] extended the analysis to FIDEs of neutral-type,(1.2)ddt[x(t)-Nx(t-τ)]=F(t,x(t),x(t-τ),∫t-τtg(t,ξ,x(ξ))dξ),t≥t0,x(t)=φ(t),t0-τ≤t≤t0,
while Zhang et al. [4, 5] derived several improved numerical stability results for neutral FIDEs of the form (1.2). For the class of neutral FIDEs:(1.3)ddt[x(t)-∫0ta(t-ξ)G(ξ,x(ξ-τ))dξ]=F(t,x(t)),t≥t0,x(t)=φ(t),t0-τ≤t≤t0,
Brunner and Vermiglio [6] made an insight into the analytical and numerical stability of continuous Runge-Kutta methods. Moreover, in the papers [7, 8], Brunner presented superconvergence results of collocation methods for several classes of FIDEs with constant or variable (vanishing) delays.The reader may wish to consult Baker’s survey paper [9] and Brunner’s monograph [10] for details on related earlier work and for additional references.However, up to now, no numerical investigation appears to have been carried out for general nonlinear FIDEs of the form:(1.4)ddt[x(t)-∫t-τtg(t,ξ,x(ξ))dξ]=f(t,x(t),x(t-τ)),t≥t0,x(t)=φ(t),t0-τ≤t≤t0,
in which the integral term on the left-hand side is no longer pure delay type: in contrast to the FIDEs (1.1)–(1.3), it contains information on the solution x on the interval [t-τ,t]. Hence, the numerical analysis and computational solution of (1.4) is rather more complex than is the case for (1.1)–(1.3).In the present paper, with an adaptation of the underlying Pouzet-Runge-Kutta methods (cf. Brunner and van der Houwen [11]), we obtain a class of new numerical methods for nonlinear FIDEs (1.4) and study analytical and numerical stability of the equations. The paper is organized as follows. In Section 2 we derive results on the asymptotic stability of analytical solutions, under the assumption of nonclassical Lipchitz conditions. Section 3 describes the adaptation of the Pouzet-Runge-Kutta method to the FIDE (1.4). In Section 4 some lemmas are given which will play a key role in the analysis of the global and asymptotical stability properties of the Pouzet-Runge-Kutta solutions (Section 5). Here, we also state stability results for a number of concrete methods. Some numerical experiments are given in Section 6 to illustrate the theoretical results and the effectiveness of the methods. Finally, in Section 7, a comparison between the presented methods and the existed related methods is given.


## 2. Stability Results for Exact Solutions

Let〈·,·〉 and ∥·∥ denote a given inner product and its induced norm on the d-dimensional complex space ℂd, respectively. The functions φ:[t0-τ,t0]→ℂd,f:[t0,+∞)×ℂd×ℂd→ℂdandg:𝔻×ℂd→ℂd(with𝔻={(t,s):t∈[t0,+∞),s∈[t-τ,t]}) are assumed to be continuous and possess the properties(2.1)R〈f(t,u,v)-f(t,ũ,ṽ),u-ũ-(w-w̃)〉≤α‖u-ũ‖2+β‖v-ṽ‖2+γ‖w-w̃‖2,u,ũ,v,ṽ,w,w̃∈Cd,(2.2)‖g(t,ξ,u)-g(t,ξ,ũ)‖≤η‖u-ũ‖,(t,ξ)∈D,u,ũ∈Cd,
where -α,β,γ,η are nonnegative constants. We will refer to the class of FIDEs of the form (1.4) which satisfies (2.1)-(2.2) as FIDEs of class𝔽𝕀𝔻(α,β,γ,η).In order to study the stability of solutions to (1.4), we need to consider the system with a different initial function ψ(t),(2.3)ddt[x̃(t)-∫t-τtg(t,ξ,x̃(ξ))dξ]=f(t,x̃(t),x̃(t-τ)),t≥t0,x̃(t)=ψ(t),t0-τ≤t≤t0,
which also belongs to the class 𝔽𝕀𝔻(α,β,γ,η).Definition 2.1.
System (1.4) is called globally stable if there exists a constant C>0 such that
(2.4)‖x(t)-x̃(t)‖≤Cmax⁡ξ∈[t0-τ,t0]‖φ(ξ)-ψ(ξ)‖,t≥t0.
Moreover, system (1.4) is called asymptotically stable if
(2.5)lim⁡t→+∞‖x(t)-x̃(t)‖=0.
In order to gain insight into the global and asymptotical stability of system (1.4), we will use the generalized Halanay inequality as presented in Wang [12], compare also to [13].Lemma 2.2 (see [12]).
Assume that the functionsu,v:[t0,+∞)→ℝ satisfy the inequalities
(2.6)u′(t)≤-a(t)u(t)+b(t)max⁡ξ∈[t-τ,t]u(ξ)+c(t)max⁡ξ∈[t-τ,t]v(ξ),t≥t0,v(t)≤d(t)max⁡ξ∈[t-τ,t]u(ξ)+e(t)max⁡ξ∈[t-τ,t]v(ξ),t≥t0.Here, a,b,c,d,e are given nonnegative continuous functions on [t0,+∞) for which there exist constants â,b̂,ĉ such that
(2.7)a(t)≥â>0,e(t)≤b̂<1,b(t)a(t)+c(t)d(t)a(t)[1-e(t)]≤ĉ<1,∀t≥t0.
Then the following inequalities hold:
(2.8)u(t)≤max⁡ξ∈[t0-τ,t0]u(ξ)exp⁡[σ̂(t-t0)],∀t≥t0,v(t)≤max⁡ξ∈[t0-τ,t0]v(ξ)exp⁡[σ̂(t-t0)],∀t≥t0,
where
(2.9)σ̂:=sup⁡t≥t0{σ(t):σ(t)+a(t)-b(t)e-σ(t)τ-c(t)d(t)e-2σ(t)τ1-e(t)e-σ(t)τ=0}<0.With this lemma, we will be able to obtain an analytical stability result for the system (1.4). In order to do so, we introduce some notations:(2.10)z(t)=∫t-τtg(t,ξ,x(ξ))dξ,z̃(t)=∫t-τtg(t,ξ,x̃(ξ))dξ,y(t)=x(t)-z(t),ỹ(t)=x̃(t)-z̃(t),x̂(t)=‖x(t)-x̃(t)‖2,ŷ(t)=‖y(t)-ỹ(t)‖2,ẑ(t)=‖z(t)-z̃(t)‖2.Theorem 2.3.
Assume that the system (1.4) belongs to the class 𝔽𝔻𝕀(α,β,γ,η) with
(2.11)α<0,2η2τ2<1,4[β+(γ-α)η2τ2]<-α(1-2η2τ2).
Then this system is globally and asymptotically stable.Proof.
It follows from (2.1), (1.4), and (2.3) that
(2.12)ŷ′(t)=2R〈f(t,x(t),x(t-τ))-f(t,x̃(t),x̃(t-τ)),y(t)-ỹ(t)〉≤2[αx̂(t)+βx̂(t-τ)+γẑ(t)],∀t≥t0.
By (2.2), it holds that
(2.13)ẑ(t)≤[η∫t-τtx̂(ξ)dξ]2≤η2τ2max⁡ξ∈[t-τ,t]x̂(ξ),∀t≥t0.
Substituting (2.13) into (2.12) yields
(2.14)ŷ′(t)≤2αx̂(t)+2(β+γη2τ2)max⁡ξ∈[t-τ,t]x̂(ξ),∀t≥t0.
Note that
(2.15)ŷ(t)=∥x(t)-x̃(t)-[z(t)-z̃(t)]∥2≤2[x̂(t)+ẑ(t)],∀t≥t0.
This, together with α<0 and (2.13), implies that
(2.16)2αx̂(t)≤α[ŷ(t)-2ẑ(t)]≤αŷ(t)-2αη2τ2max⁡ξ∈[t-τ,t]x̂(ξ),∀t≥t0.
Hence, by combining (2.14) and (2.16) we are led to
(2.17)ŷ′(t)≤αŷ(t)+2(β+γη2τ2-αη2τ2)max⁡ξ∈[t-τ,t]x̂(ξ),∀t≥t0.
On the other hand, we have
(2.18)x̂(t)=‖y(t)-ỹ(t)+[z(t)-z̃(t)]‖2≤2[ŷ(t)+ẑ(t)]≤2ŷ(t)+2η2τ2max⁡ξ∈[t-τ,t]x̂(ξ),∀t≥t0,
where we have used (2.13). Therefore, under the condition (2.11), an application of Lemma 2.2 to (2.17)-(2.18) yields the conclusion.As an application of Theorem2.3, we present several examples as follows.Example 2.4.
Consider thed-dimensional system of linear functional-integrodifferential equation
(2.19)ddt[x(t)-∫t-τtN(t,ξ)x(ξ)dξ]=L(t)x(t)+M(t)x(t-τ)+G(t),t≥t0,x(t)=φ(t),t0-τ≤t≤t0,
where G:[t0,+∞)→ℂd is a known continuous function such that (2.19) has a unique solution. It is easy to check that system (2.19) belongs to the class 𝔽𝕀𝔻(μ+(ι+κ)/2,κ,(ι+κ)/2,ϱ) if, provided there exist constants μ,ι,κ,ϱ such that, for all t≥t0 and (t,ξ)∈𝔻,
(2.20)μ(L(t))≤μ≤-(ι+κ)2,‖L(t)‖≤ι,‖M(t)‖≤κ,‖N(t,ξ)‖≤ϱ,
in which the matrix norm ∥·∥ and the logarithmic norm μ(·) are induced by the vector inner-product norm. A direct application of Theorem 2.3 shows that the system (2.19) is globally and asymptotically stable whenever condition (2.20) and the following condition hold:
(2.21)2ϱ2τ2<1,4(κ-μϱ2τ2)<-[μ+(ι+κ)2](1-2ϱ2τ2).Example 2.5.
Consider the system of partial functional-integrodifferential equations(2.22)∂∂t[u(v,t)-14∫t-3/4texp⁡(ξ-t)u(v,ξ)dξ]=-12u(v,t)+u(v,t-3/4)24[1+u2(v,t-3/4)]+g(v,t),t≥0,v∈(0,2π),u(v,t),u(v,t)=sin⁡⁡vexp⁡(-vt),t∈[-34,0],v∈[0,2π],u(0,t)=u(2π,t)=0,t≥0,
where g(v,t) is a continuous function chosen such that this system has the exact solution u(v,t)=sinvexp(-vt). By discretizing the spatial variable v by a uniform mesh vi=iΔv(i=0,1,…,l,Δv=2π/l), the system (2.22) can be transformed into a system of ordinary functional-integrodifferential equations,
(2.23)ddt[x(t)-14∫t-3/4texp⁡⁡(ξ-t)x(ξ)dξ]=-12x(t)+124x¯(t-34)+G(t),t≥0,x(t)=(sin⁡v1exp⁡(-v1t),sin⁡v2exp⁡(-v2t),…,sin⁡vl-1exp⁡(-vl-1t))T,t∈[-34,0],
where
(2.24)x(t)=(x1(t),x2(t),…,xl-1(t))T,G(t)=(g(v1,t),g(v2,t),…,g(vl-1,t))T,x¯(t)=(x1(t)1+x12(t),x2(t)1+x22(t),…,xl-1(t)1+xl-12(t))T,xi(t)≈u(vi,t).
When the standard inner product and its induced norm are used, one easily verifies that the system (2.23) belongs to the class 𝔽𝕀𝔻(-11/48,1/24,13/48,1/4). Moreover, in light of Theorem 2.3, we know that the system (2.23) is globally and asymptotically stable.


## 3. The Pouzet-Runge-Kutta Discretization

In order to obtain a class of effective numerical schemes for FIDEs of the form (1.4), we first recall some related concepts and results on the underlying Runge-Kutta methods for ordinary differential equations (see, e.g., [14]). An s-stage Runge-Kutta method is described by the Butcher tableau(3.1)cAbT
where(3.2)A:=(aij)∈Rs×s,b:=(b1,b2,…,bs)T∈Rs,∑i=1sbi=1.c=(c1,c2,…,cs)T,0≤ci≤1(i=1,2,…,s).A Runge-Kutta method (3.1) is called algebraically stable if(3.3)D:=diag⁡(b1,b2,…,bs)≥0,M:=DA+ATD-bbT≥0,where the notation “≥” means that a matrix is nonnegative definite. It is said to bestrictly stable at infinity if R(∞):=limz→∞R(z) exists and satisfies |R(∞)|<1, where(3.4)R(z)=1+zbT(Is-zA)-1e(z∈C),e=(1,1,…,1)T∈Rs,andIs denotes the s×s identity matrix; it is said DJ-irreducible if there is no nonempty index set 𝕃⊂{1,2,…,s} such that(3.5)bj=0forj∈L,aij=0fori∉L,j∈L,S-irreducible if there is no partition (𝕊1,𝕊2,…,𝕊r) of {1,2,…,s} with r<s such that for all l and m(3.6)∑k∈Smaik=∑k∈Smajk,fori,j∈Sl,andirreducible if it is both DJ-irreducible and S-irreducible.Proposition 3.1 (cf. [15]).
A DJ-irreducible, algebraically stable Runge-Kutta methods satisfiesbi>0  for all i.Proposition 3.2 (cf. [14]).
Assume that a Runge-Kutta method (3.1) with distinct ci and positive bi satisfies the simplifying condition B(2s-2),C(s-1),D(s-1). Then this method is algebraically stable if and only if  |R(∞)|≤1.The class ofextended Pouzet-Runge-Kutta methods for FIDEs (1.4) with underlying Runge-Kutta method (3.1) is given by(3.7)xi(n)-zi(n)=xn-zn+h∑j=1saijf(tj(n),xj(n),xj(n-m)),i=1,2,…,s,xn+1-zn+1=xn-zn+h∑j=1sbjf(tj(n),xj(n),xj(n-m)),n≥0.
Here, the stepsize h is chosen as h=τ/m(m∈ℕ),  tn:=t0+nh,tj(n):=tn+cjh, xn, zn, xi(n) and zi(n) are approximations to x(tn), z(tn), x(ti(n)) and z(ti(n)), respectively, with z(t) denoting the memory term(3.8)z(t):=∫t-τtg(t,ξ,x(ξ))dξ.The integral approximationszn and zj(n) are given by the Pouzet quadrature rules(3.9)zn=h∑q=n-mn-1∑j=1sbjg(tn,tj(q),xj(q)),n≥0,(3.10)zi(n)=h∑j=1saijg(ti(n),tj(n),xj(n))+h∑q=n-mn-1∑j=1sbjg(ti(n),tj(q),xj(q))-h∑j=1saijg(ti(n),tj(n-m),xj(n-m)),i=1,2,…,s;n≥0.
Moreover, it is assumed that(3.11)xn=φ(tn),xi(n)=φ(ti(n))fort0-τ≤tn,ti(n)≤t0.


## 4. Some Basic Lemmas

In the subsequent numerical analysis we will imply the following notations:(4.1)x̂n:=xn-x̃n,ẑn:=zn-z̃n,yn:=xn-zn,ỹn:=x̃n-z̃n,ŷn:=yn-ỹn,x̂i(n):=xi(n)-x̃i(n),ẑi(n):=zi(n)-z̃i(n),yi(n):=xi(n)-zi(n),ỹi(n):=x̃i(n)-z̃i(n),ŷi(n):=yi(n)-ỹi(n),f̂i(n):=f(ti(n),xi(n),xi(n-m))-f(ti(n),x̃i(n),x̃i(n-m)).Moreover, when the method (3.7)–(3.10) is applied to the problem (2.3), we set(4.2)z̃(t):=∫t-τtg(t,ξ,x̃(ξ))dξand denote the corresponding approximations ofx̃(tn),z̃(tn),x̃(ti(n)), and z̃(ti(n)), by x̃n,z̃n,x̃i(n), and z̃i(n), respectively. Similarly, it is assumed that(4.3)x̃n=ψ(tn),x̃i(n)=ψ(ti(n)),fort0-τ≤tn,ti(n)≤t0.As mentioned in Section1, the following lemmas will play key roles in derivation of our main results.Lemma 4.1.
Assume that the underlying Runge-Kutta method (3.1) is algebraically stable, and the conditions (2.1)-(2.2) hold. Then the extended Pouzet-Runge-Kutta scheme (3.7) satisfies
(4.4)‖ŷn+1‖2≤(2βτ+η2τ2)max⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2+2h∑i=0n∑j=1sbj[(α+β)‖x̂j(i)‖2+γ‖ẑj(i)‖2],n≥0.Proof.
A straightforward computation and the assumption of algebraic stability lead to(4.5)‖ŷn+1‖2-‖ŷn‖2-2h∑j=1sbjR〈ŷj(n),f̂j(n)〉=-h2∑i=1s∑j=1smij〈f̂i,f̂j〉≤0,∀n≥0,
where M=(mij). Hence, it holds that
(4.6)‖ŷn+1‖2≤‖ŷn‖2+2h∑j=1sbjR〈ŷj(n),f̂j(n)〉,∀n≥0.
By (2.1), we further have for n≥0,
(4.7)‖ŷn+1‖2≤‖ŷn‖2+2h∑j=1sbj[α‖x̂j(n)‖2+β‖x̂j(n-m)‖2+γ‖ẑj(n)‖2].
An induction argument shows that the inequality (4.7) implies
(4.8)‖ŷn+1‖2≤‖ŷ0‖2+2h∑i=0n∑j=1sbj[α‖x̂j(i)‖2+β‖x̂j(i-m)‖2+γ‖ẑj(i)‖2],n≥0,
in which
(4.9)‖ŷ0‖2=‖x0-z0-(x̃0-z̃0)‖2=‖∫t0-τt0[g(t0,ξ,φ(ξ))-g(t0,ξ,ψ(ξ))]dξ‖2≤η2τ2max⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2.
Also, the assumptions mh=τ and ∑j=1sbj=1 allow us to write
(4.10)h∑i=0n∑j=1sbj‖x̂j(i-m)‖2=h∑j=1sbj(∑i=-m-1‖x̂j(i)‖2+∑i=0n-m‖x̂j(i)‖2)≤τmax⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2+h∑i=0n∑j=1sbj‖x̂j(i)‖2,n≥0.
A combination of (4.8), (4.9), and (4.10) yields (4.4). Hence the lemma is proven.Lemma 4.2.
Under the condition (2.2), the Pouzet quadrature rule (3.9) satisfies
(4.11)‖ẑn‖2≤hη2τ(∑j=1s|bj|2)(∑i=0n-1∑j=1s‖x̂j(i)‖2)+η2τ2s(∑j=1s|bj|)max⁡t0-τ≤ξ≤t0⁡‖φ(ξ)-ψ(ξ)‖2,n≥0.Proof.
By (2.2), mh=τ and the Cauchy-Schwarz inequality, we have
(4.12)‖ẑn‖2≤(hη∑i=n-mn-1∑j=1s|bj|‖x̂j(i)‖)2≤h2η2(∑i=n-mn-1∑j=1s|bj|2)(∑i=n-mn-1∑j=1s‖x̂j(i)‖2)=hη2τ(∑j=1s|bj|2)(∑i=n-mn-1∑j=1s‖x̂j(i)‖2).
Also, it follows from (3.11), (4.3), and mh=τ that
(4.13)h∑i=n-mn-1∑j=1s‖x̂j(i)‖2≤τsmax⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2,forn=0,(4.14)h∑i=n-mn-1∑j=1s‖x̂j(i)‖2≤h(∑i=-m-1∑j=1s‖x̂j(i)‖2+∑i=0n-1∑j=1s‖x̂j(i)‖2)≤τsmax⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2+h∑i=0n-1∑j=1s‖x̂j(i)‖2,for1≤n<m,(4.15)h∑i=n-mn-1∑j=1s‖x̂j(i)‖2≤h∑i=0n-1∑j=1s‖x̂j(i)‖2,forn≥m.
A combination of (4.13)–(4.15) yields that for all n≥0,
(4.16)h∑i=n-mn-1∑j=1s‖x̂j(i)‖2≤τsmax⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2+h∑i=0n-1∑j=1s‖x̂j(i)‖2.
Inserting (4.16) into (4.12) generates (4.11). This completes the proof.Lemma 4.3 (cf. [5]).
Under the condition (2.2), the Pouzet quadrature rule (3.10) satisfies
(4.17)‖ẑj(i)‖2≤3h2η2[(∑r=1s|ajr|2)(∑r=1s‖x̂r(i)‖2)+m∑q=1m(∑r=1s|br|2)(∑r=1s‖x̂r(i-q)‖2)+(∑r=1s|ajr|2)(∑r=1s‖x̂r(i-m)‖2)],i≥0,j=1,2,…,s.


## 5. Numerical Stability

This section will involve the analysis of the global and asymptotical stability properties of the Pouzet-Runge-Kutta method (3.7)–(3.10). The relevant numerical stability concepts are defined as follows.Definition 5.1.
The Pouzet-Runge-Kutta method (3.7)–(3.10) is called globally stable for the problems of class 𝔽𝕀𝔻(α,β,γ,η) if there exists a stability constant ℋ>0, which depends only on α,β,γ,η,τ and the method, such that
(5.1)‖xn-x̃n‖≤Hmax⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖,∀n≥1.
Moreover, the method (3.7)–(3.10) is called asymptotically stable for the problems of class 𝔽𝕀𝔻(α,β,γ,η) if
(5.2)lim⁡n→∞‖xn-x̃n‖=0.
We first establish a result on the global stability of Pouzet-Runge-Kutta methods. An analogous result on their asymptotical stability will be given in Theorem5.7.Theorem 5.2.
Assume that the underlying Runge-Kutta method (3.1) with positive bi is algebraically stable. Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is globally stable for the class 𝔽𝕀𝔻(α,β,γ,η), with stability constant
(5.3)H=2τ[2β+η2τ(1+s)+6sγη2τ2∑j=1sbj(bj+∑r=1s|ajr|2)],
whenever
(5.4)2(α+β)+η2τσ∑j=1sbj[(1+6γτ)bj+12γτ∑r=1s|ajr|2]≤0,
where σ=min1≤i≤s{bi}>0.Proof.
It follows from Lemma4.3 and ∑j=1sbj=1(bj>0) that
(5.5)h∑i=0n∑j=1sbj‖ẑj(i)‖2≤3h3η2[(∑j=1sbj∑r=1s|ajr|2)(∑i=0n∑r=1s‖x̂r(i)‖2)+m(∑r=1sbr2)(∑r=1s∑i=0n∑q=1m‖x̂r(i-q)‖2)+(∑j=1sbj∑r=1s|ajr|2)(∑i=0n∑r=1s‖x̂r(i-m)‖2)].
Also, we have that
(5.6)∑i=0n∑r=1s‖x̂r(i-m)‖2=∑i=-mn-m∑r=1s‖x̂r(i)‖2≤∑i=0n∑r=1s‖x̂r(i)‖2+m∑r=1smax⁡-m≤i≤-1‖x̂r(i)‖2
and, by the inequality (3.13) in [1],
(5.7)∑i=0n∑q=1m‖x̂r(i-q)‖2≤m∑i=0n‖x̂r(i)‖2+m(m+1)2max⁡-m≤i≤-1‖x̂r(i)‖2.
Inserting both (5.6) and (5.7) into (5.5) yields
(5.8)h∑i=0n∑j=1sbj‖ẑj(i)‖2≤3h3η2[∑j=1sbj(2∑r=1s|ajr|2+m2bj)∑i=0n∑r=1s‖x̂r(i)‖2+m∑j=1sbj(m(m+1)2bj+∑r=1s|ajr|2)∑r=1smax⁡-m≤i≤-1‖xr(i)‖2].
With h=τ/m≤τ, (3.11) and (4.3), the left side of (5.8) can be bounded by
(5.9)3η2τ2[hσ∑j=1sbj(2∑r=1s|ajr|2+bj)∑i=0n∑r=1sbr‖x̂r(i)‖2+τs∑j=1sbj(bj+∑r=1s|ajr|2)max⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2].
Substituting this bound into (4.4) yields
(5.10)‖ŷn+1‖2≤τ[2β+η2τ+6sγη2τ2∑j=1sbj(bj+∑r=1s|ajr|2)]max⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2+2h[α+β+3γη2τ2σ∑j=1sbj(bj+2∑r=1s|ajr|2)]∑i=0n∑r=1sbr‖x̂r(i)‖2,n≥0,
which, together with (4.11), (5.10), and ∑j=1sbj=1(bj>0), implies that for n≥1,
(5.11)‖x̂n‖2=‖ŷn+ẑn‖2≤2(‖ŷn‖2+‖ẑn‖2)≤2τ[2β+η2τ(1+s)+6sγη2τ2∑j=1sbj(bj+∑r=1s|ajr|2)]max⁡t0-τ≤ξ≤t0‖φ(ξ)-ψ(ξ)‖2+2h[2(α+β)+η2τσ∑j=1sbj((1+6γτ)bj+12γτ∑r=1s|ajr|2)]∑i=0n-1∑r=1sbr‖x̂r(i)‖2.
Observing the condition (5.4) to (5.11), we readily derive the global stability conclusion. Hence, the theorem is proven.Applying Proposition3.1 and Proposition 3.2 to Theorem 5.2, respectively, we can obtain the following corollaries.Corollary 5.3.
Assume that an underlying Runge-Kutta method is DJ-irreducible and algebraically stable. Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is globally stable for the class 𝔽𝕀𝔻(α,β,γ,η) whenever (5.4) holds.Corollary 5.4.
Assume that an underlying Runge-Kutta method, with distinctci and positive bi, satisfies the simplifying condition B(2s-2),C(s-1),D(s-1), and |R(∞)|≤1. Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is globally stable for the class 𝔽𝕀𝔻(α,β,γ,η) whenever (5.4) holds.Remark 5.5.
In [14, Chapters IV.5 and IV.12], Hairer and Wanner have shown that the Runge-Kutta methods of type Gauss, Radau IA, Radau IIA, and Lobatto IIIC are all algebraically stable and have invertible matrices A, distinct ci, and bi>0 for all i. Hence, by Theorem 5.2, we have a more concrete stability result for the Pouzet-Runge-Kutta schemes based on these important Runge-Kutta methods.Corollary 5.6.
Assume that an underlying Runge-Kutta method (3.1) is of type Gauss, Radau IA, Radau IIA, or Lobatto IIIC. Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is globally stable for the class 𝔽𝕀𝔻(α,β,γ,η) whenever (5.4) holds.Next, we will address the asymptotic stability of the extended Pouzet-Runge-Kutta methods. The notation,(5.12)X̂n:=(x̂1(n)T,…,x̂s(n)T)T,Ẑn:=(ẑ1(n)T,…,ẑs(n)T)T,F̂n:=(f̂1(n)T,…,f̂s(n)T)T,will subsequently be used, and we define a vector norm on the spaceℂds by(5.13)‖U‖=∑r=1s‖ur‖2,∀U=(u1T,u2T,…,usT)T∈Cds.Moreover, we will employ the Kronecker product⊗ and its well-known properties (cf. [16]).Theorem 5.7.
Assume that an underlying Runge-Kutta method (3.1) is irreducible, algebraically stable, and strictly stable at infinity. Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is asymptotically stable for the class 𝔽𝔻𝕀(α,β,γ,η) whenever (5.4) holds.Proof.
By (3.7) we have
(5.14)X̂n-Ẑn=(e⊗Id)(x̂n-ẑn)+h(A⊗Id)F̂n,x̂n+1-ẑn+1=x̂n-ẑn+h(bT⊗Id)F̂n,n≥0.
When the matrix A is invertible, it follows from the first equation of (5.14) that
(5.15)hF̂n=(A-1⊗Id)(X̂n-Ẑn)-(A-1e⊗Id)(x̂n-ẑn).
Insertion of (5.15) into the second equation of (5.14) yields
(5.16)x̂n+1-ẑn+1=(1-bTA-1e)(x̂n-ẑn)+(bTA-1⊗Id)(X̂n-Ẑn).
When A is singular, we adopt a technique, proposed by Hairer and Wanner [14]. It consists in replacing A by the regular matrix A+ɛIs everywhere, followed by considering the limit limɛ→0(A+ɛIs)-1, whose existence is assured by irreducibility and algebraic stability of the method (see Lemma  3.2 in [14]). Thus, we set
(5.17)A-1:={A-1,ifAisinvertible,lim⁡ϵ→0(A+ϵIs)-1,ifAissingular.
Hence, it holds that
(5.18)x̂n+1-ẑn+1=R(∞)(x̂n-ẑn)+(bTA-1⊗Id)(X̂n-Ẑn),
where R(∞)=1-bT𝒜-1e. Applying Lemma 4.6 in [2] to (5.18) we obtain
(5.19)lim⁡n→∞‖x̂n+1-ẑn+1‖=0,
if and only if
(5.20)|R(∞)|<1,lim⁡n→∞‖X̂n-Ẑn‖=0.
Next, we prove (5.19). We need only to show that(5.21)lim⁡n→∞‖X̂n-Ẑn‖=0,
since |R(∞)|<1 is a known condition. In fact, by Proposition 3.1 we know that bi>0 for all i. According to this and (5.4), we derive from (5.11) that
(5.22)lim⁡n→∞‖x̂r(n)‖=0,r=1,2,…,s.
This implies that
(5.23)lim⁡n→∞‖X̂n‖=lim⁡n→∞∑r=1s‖x̂r(n)‖2=0.
Combining (4.17) and (5.22) yields limn→∞∥ẑj(n)∥=0 for all j, which leads to
(5.24)lim⁡n→∞‖Ẑn‖=lim⁡n→∞∑j=1s‖ẑj(n)‖2=0.
A combination of (5.23) and (5.24) gives (5.21). Hence, (5.19) is true.
Finally, we prove thatlimn→∞∥x̂n∥=0. By (3.9), (2.2), and (5.22), we have(5.25)‖ẑn‖≤hη∑q=n-mn-1∑j=1sbj‖x̂j(q)‖⟶0,asn⟶∞,
which means that limn→∞∥ẑn∥=0. This, together with (5.19), implies that
(5.26)‖x̂n‖≤‖x̂n-ẑn‖+‖ẑn‖⟶0,asn⟶∞.
Accordingly, the theorem is proven.The proof of the above theorem reveals that the irreducibility of the method is used only for the case where the Runge-Kutta matrixA is singular. Hence, for Pouzet-Runge-Kutta methods whose underlying Runge-Kutta has an invertible the matrix A, the irreducibility condition in Theorem 5.7 can be dropped. This is made precise in the following theorem.Theorem 5.8.
Assume that an underlying Runge-Kutta method (3.1) with invertible matrix A and positive bi is algebraically stable and strictly stable at infinity. Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is asymptotically stable for the class 𝔽𝔻𝕀(α,β,γ,η) whenever the condition (5.4) holds.In light of Theorem5.7, Propositions 3.1, and Propositions 3.2, we obtain the following analogues of Corollaries 5.3 and 5.4.Corollary 5.9.
Assume that an underlying Runge-Kutta method (3.1) with distinct ci is irreducible, algebraically stable, and satisfies |R(∞)|≠1 and the simplifying condition B(2s-2),C(s-1),D(s-1). Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is asymptotically stable for the class 𝔽𝔻𝕀(α,β,γ,η) whenever (5.4) holds.Corollary 5.10.
Assume that an underlying Runge-Kutta method (3.1) with distinct ci and positive bi is irreducible, strictly stable at infinity, and satisfies the simplifying condition B(2s-2),C(s-1),D(s-1). Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is asymptotically stable for the class 𝔽𝔻𝕀(α,β,γ,η) whenever (5.4) holds.Similar to Theorem5.8, the irreducibility condition can be dropped in Corollary 5.10 when the matrix A is invertible. Moreover, by Remark 5.5, Theorem 5.8, and the fact that the stability functions of Radau IA, Radau IIA, and Lobatto IIIC methods all satisfy R(∞)=0 (cf. [14]), one can establish an asymptotic stability result analogous to the one in Corollary 5.6.Corollary 5.11.
Assume that the underlying Runge-Kutta method (3.1) is of type Radau IA, Radau IIA, or Lobatto IIIC. Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is asymptotically stable for the class 𝔽𝔻𝕀(α,β,γ,η) whenever (5.4) holds.


## 6. Numerical Illustration

In order to illustrate the effectiveness of the extended Pouzet-Runge-Kutta methods (3.7)–(3.10), we will apply the two-stage methods of type Gauss, Radau IA, Radau IIA, or Lobatto IIIC to the system (2.23), respectively, where the solution domain is chosen as [0,43]. These methods produce a series of high-precision numerical solutions for the (2.22) on [0,43;0,2π].Let(6.1)Γ:=2(α+β)+η2τmin⁡1≤i≤s⁡{bi}∑j=1sbj[(1+6γτ)bj+12γτ∑r=1s|ajr|2].The correspondingΓ-values of the above methods applied to the system (2.23) are listed in the Table 1(a), indicating that the methods satisfy the condition (5.4). Hence, the methods are globally stable by Corollary 5.6, and asymptotically stable by Corollary 5.11, except for the Gauss-type method.Table 1
(a) TheΓ-values of the methods for (2.23), (b) The errors of the numerical solutions of (2.22) by methods (3.7)–(3.10).
(a)GaussRadau IARadau IIALobatto IIIC-3.1302e-001-2.2800e-001-2.1530e-001-2.9081e-001(b)(Δv,h)GaussRadau IARadau IIALobatto IIIC(π/30,3/16)9.1064e-0063.7388e-0043.8044e-0041.3926e-002(π/60,3/32)5.7295e-0074.6647e-0054.8352e-0053.5011e-003(π/120,3/64)3.5879e-0085.8133e-0066.0818e-0068.7739e-004(π/240,3/128)2.2433e-0097.2537e-0077.6236e-0072.1957e-004The excellent stability properties of the methods lead us to expect good numerical results. In order to confirm this, we use the Newton-Raphson iteration technique to implement the above numerical schemes. Taking the following four groups of space-time stepsizes:(6.2)(Δv,h)=(π30,316),(π60,332),(π120,364),(π240,3128),
and then applying the above Pouzet-Runge-Kutta schemes to the system (2.23) on [0,43], respectively, we obtain sixteen sets of numerical solutions. The numerical solution generated by the extended two-stage Gauss-type method with space-time stepsizes (π/240,3/128) is plotted in Figure 1. The solution figures for the other methods are quite similar to Figure 1, and hence we omit them here. In order to show the computational precision of the obtained numerical solutions, we use(6.3)err⁡:=max⁡1≤n≤N‖xn-û(tn)‖∞Figure 1
Numerical solution of (2.22) by 2-stage Gauss-type method.to characterize the errors of the methods, where(6.4)û(tn):=(sin⁡v1exp⁡(-v1tn),sin⁡v2exp⁡(-v2tn),…,sin⁡vl-1exp⁡(-vl-1tn))Tis the vector whose entries consist of the system (2.22) true solutions at the meshpoints (vi,tn), i=1,2,…,l-1. The errors of the above four methods with different stepsizes are displayed in Table 1(b); they confirm again that the methods are rather effective.


## 7. Concluding Remarks

In paper [1], with a combination of Runge-Kutta methods and Pouzet quadrature rules, the authors also obtained an alternative type of numerical methods for (1.1), namely,(7.1)xi(n)=xn+h∑j=1saijF(tj(n),xj(n),xj(n-m),zj(n)),i=1,2,…,s,xn+1=xn+h∑j=1sbjF(tj(n),xj(n),xj(n-m),zj(n)),n≥0,
where the meanings of the notations are similar to those indicated in method (3.7) and {zj(n)}j=1s are determined by (3.10). But such methods cannot be applied directly to systems (1.4) unless the integrated function g is continuously differentiable on its domain. When the latter holds, the system (1.4) can be transformed into(7.2)x′(t)=F(t,x(t),x(t-τ),∫t-τtgt(t,ξ,x(ξ))dξ),t≥t0,x(t)=φ(t),t0-τ≤t≤t0,
where(7.3)F(t,x,y,z)=f(t,x,y)+g(t,t,x)-g(t,t-τ,y)+z.This is just of the form (1.1) and implies, under the condition that g is continuously differentiable, that systems (1.4) can be solved by methods {(7.1), (3.10)}. Thus, the numerical stability theory in [1] is applicable and hence a series of global and asymptotical stability results can be followed immediately. However, it is difficult to compare the theoretical results in this way with those in previous sections since both discretization schemes are different and there is no direct relationship between the condition (2.1)-(2.2) and the condition imposed on (1.1) (see [1]). Moreover, we have noted that a system (1.4) can be solved by scheme {(7.1), (3.10)} only when g continuously differentiable, which shows that schemes (3.7)–(3.10) have a wider applicable range than schemes {(7.1), (3.10)} do.In the following, we give a comparison between methods (3.7)–(3.10) and methods {(7.1), (3.10)} with some numerical experiments. It is evident that the system (2.23) can be changed into the form (1.1). Hence this system also can be solved by methods {(7.1), (3.10)}. Similarly, we take the space-time stepsizes in (6.2) and apply the two-stage methods {(7.1), (3.10)} of type Gauss, Radau IA, Radau IIA, and Lobatto IIIC to the system (2.23) on [0,43], respectively, then a series of high-precision numerical solutions can be worked out, whose errors are displayed in Table 2(a). It follows from Tables 1(b) and 2(a)–2(c) that the numerical precisions and the computational times of the both methods based on the same type of underlying Runge-Kutta method are almost similar under the same stepsize. This implies that methods (3.7)–(3.10) are comparable.Table 2
(a) The errors of the numerical solutions of (2.22) by methods {(7.1), (3.10)}, (b) The CPU times (in second) of methods (3.7)–(3.10) for (2.23) and (c) The CPU times (in second) of methods {(7.1), (3.10)} for (2.23).
(a)(Δv,h)GaussRadau IARadau IIALobatto IIIC(π/30,3/16)9.0193e-0063.7273e-0043.8966e-0041.4083e-002(π/60,3/32)5.6739e-0074.6491e-0054.9558e-0053.5453e-003(π/120,3/64)3.5534e-0085.7933e-0066.2360e-0068.8912e-004(π/240,3/128)2.2217e-0097.2284e-0077.8186e-0072.2261e-004(b)(Δv,h)GaussRadau IARadau IIALobatto IIIC(π/30,3/16)7.0780e+0006.8750e+0006.9060e+0006.9380e+000(π/60,3/32)4.4468e+0014.4500e+0014.4766e+0014.5172e+001(π/120,3/64)3.3039e+0023.3156e+0023.3281e+0023.3264e+002(π/240,3/128)2.7117e+0032.7042e+0032.7109e+0032.7291e+003(c)(Δv,h)GaussRadau IARadau IIALobatto IIIC(π/30,3/16)4.2180e+0004.2350e+0004.2650e+0004.2820e+000(π/60,3/32)3.0953e+0013.0234e+0013.1031e+0013.1125e+001(π/120,3/64)2.7055e+0022.8400e+0022.7213e+0022.7105e+002(π/240,3/128)2.7821e+0032.8164e+0032.7950e+0032.8006e+003



---

*Source: 642318-2012-04-23.xml*