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

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.

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) Yu et al. [3] extended the analysis to FIDEs of neutral-type,𝑑[𝑥(𝑡)((𝑡)(),𝑑𝑡−𝑁𝑥𝑡−𝜏)]=𝐹𝑡,𝑥,𝑥𝑡−𝜏^𝑡_𝑡−𝜏𝑔(𝑡,𝜉,𝑥(𝜉))𝑑𝜉,𝑡≥𝑡₀,𝑥(𝑡)=𝜑(𝑡),𝑡₀−𝜏≤𝑡≤𝑡₀, (1.2) 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) 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) 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 𝑥 on the interval [𝑡−𝜏,𝑡]. 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 𝑑-dimensional complex space ℂ^𝑑, respectively. The functions 𝜑∶[𝑡₀−𝜏,𝑡₀]→ℂ^𝑑,𝑓∶[𝑡₀,+∞)×ℂ^𝑑×ℂ^𝑑→ℂ^𝑑and𝑔∶𝔻×ℂ^𝑑→ℂ^𝑑(with𝔻={(𝑡,𝑠)∶𝑡∈[𝑡₀,+∞),𝑠∈[𝑡−𝜏,𝑡]}) are assumed to be continuous and possess the propertiesℜ𝑓()()(‖‖𝑡,𝑢,𝑣−𝑓𝑡,𝑢,𝑣,𝑢−𝑢−𝑤−𝑤)≤𝛼𝑢−𝑢²+𝛽‖𝑣−𝑣‖²+𝛾‖𝑤−𝑤‖²,𝑢,𝑢,𝑣,𝑣,𝑤,𝑤∈ℂ^𝑑, (2.1)𝑔(𝑡,𝜉,𝑢)−𝑔(𝑡,𝜉,𝑢)≤𝜂‖𝑢−𝑢‖,(𝑡,𝜉)∈𝔻,𝑢,𝑢∈ℂ^𝑑, (2.2) 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 𝜓(𝑡),𝑑(𝑡)−𝑑𝑡𝑥^𝑡_𝑡−𝜏𝑔(𝑡,𝜉,𝑥(𝜉))𝑑𝜉=𝑓(𝑡,𝑥(𝑡),𝑥(𝑡−𝜏)),𝑡≥𝑡₀,x(𝑡)=𝜓(𝑡),𝑡₀−𝜏≤𝑡≤𝑡₀, (2.3) which also belongs to the class 𝔽𝕀𝔻(𝛼,𝛽,𝛾,𝜂).

Definition 2.1. System (1.4) is called globally stable if there exists a constant 𝐶>0 such that ‖𝑥(𝑡)−𝑥(𝑡)‖≤𝐶max_{𝜉∈[𝑡0−𝜏,𝑡0]}𝜑(𝜉)−𝜓(𝜉),𝑡≥𝑡₀. (2.4) Moreover, system (1.4) is called asymptotically stable if lim_𝑡→+∞‖𝑥(𝑡)−𝑥(𝑡)‖=0. (2.5)
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 functions 𝑢,𝑣∶[𝑡₀,+∞)→ℝ satisfy the inequalities 𝑢^′(𝑡)≤−𝑎(𝑡)𝑢(𝑡)+𝑏(𝑡)max_{[]𝜉∈𝑡−𝜏,𝑡}𝑢(𝜉)+𝑐(𝑡)max_{[]𝜉∈𝑡−𝜏,𝑡}𝑣(𝜉),𝑡≥𝑡₀,𝑣(𝑡)≤𝑑(𝑡)max_{[]𝜉∈𝑡−𝜏,𝑡}𝑢(𝜉)+𝑒(𝑡)max_{[]𝜉∈𝑡−𝜏,𝑡}𝑣(𝜉),𝑡≥𝑡₀. (2.6)Here, 𝑎,𝑏,𝑐,𝑑,𝑒 are given nonnegative continuous functions on [𝑡₀,+∞) for which there exist constants 𝑎,𝑏̂,𝑐̂ such that 𝑎(𝑡)(𝑡)≤̂𝑏(𝑡)≥𝑎>0,𝑒𝑏<1,𝑎(𝑡)+𝑐(𝑡)𝑑(𝑡)𝑎(𝑡)[1−𝑒(𝑡)]≤𝑐̂<1,∀𝑡≥𝑡₀. (2.7) Then the following inequalities hold: 𝑢(𝑡)≤max_{𝑡𝜉∈0−𝜏,𝑡0}𝑢(𝜉)exp𝜎𝑡−𝑡₀,∀𝑡≥𝑡₀,𝑣(𝑡)≤max_{𝑡𝜉∈0−𝜏,𝑡0}𝑣(𝜉)exp𝜎𝑡−𝑡₀,∀𝑡≥𝑡₀, (2.8) where 𝜎∶=sup_𝑡≥𝑡0𝜎(𝑡)∶𝜎(𝑡)+𝑎(𝑡)−𝑏(𝑡)𝑒^{−𝜎(𝑡)𝜏}−𝑐(𝑡)𝑑(𝑡)𝑒^{−2𝜎(𝑡)𝜏}1−𝑒(𝑡)𝑒^{−𝜎(𝑡)𝜏}=0<0. (2.9)

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)

Theorem 2.3. Assume that the system (1.4) belongs to the class 𝔽𝔻𝕀(𝛼,𝛽,𝛾,𝜂) with 𝛼<0,2𝜂²𝜏²<1,4𝛽+𝛾−𝛼𝜂²𝜏²<−𝛼1−2𝜂²𝜏². (2.11) Then this system is globally and asymptotically stable.

Proof. It follows from (2.1), (1.4), and (2.3) that (𝑡)𝑓((𝑡)(((𝑡)((𝑡)(𝑡(𝑡)()𝑧(𝑡𝑦′=2ℜ𝑡,𝑥,𝑥𝑡−𝜏))−𝑓𝑡,𝑥,𝑥𝑡−𝜏)),𝑦−𝑦)≤2𝛼𝑥+𝛽𝑥𝑡−𝜏+𝛾),∀𝑡≥𝑡₀. (2.12) By (2.2), it holds that 𝑧(𝑡)≤𝜂^𝑡_𝑡−𝜏𝑥(𝜉)𝑑𝜉²≤𝜂²𝜏²max_{[]𝜉∈𝑡−𝜏,𝑡}𝑥(𝜉),∀𝑡≥𝑡₀. (2.13) Substituting (2.13) into (2.12) yields 𝑦′(𝑡)≤2𝛼𝑥(𝑡)+2𝛽+𝛾𝜂²𝜏²max_{[]𝜉∈𝑡−𝜏,𝑡}𝑥(𝜉),∀𝑡≥𝑡₀. (2.14) Note that 𝑦(𝑡)=‖𝑥(𝑡)−𝑥(𝑡)−[𝑧(𝑡)−𝑧(𝑡)]‖²≤2[𝑥(𝑡)+𝑧(𝑡)],∀𝑡≥𝑡₀. (2.15) This, together with 𝛼<0 and (2.13), implies that 2𝛼𝑥(𝑡)≤𝛼𝑦(𝑡)−2𝑧(𝑡)≤𝛼𝑦(𝑡)−2𝛼𝜂²𝜏²max_{[]𝜉∈𝑡−𝜏,𝑡}𝑥(𝜉),∀𝑡≥𝑡₀. (2.16) Hence, by combining (2.14) and (2.16) we are led to 𝑦′(𝑡)≤𝛼𝑦(𝑡)+2𝛽+𝛾𝜂²𝜏²−𝛼𝜂²𝜏²max_{[]𝜉∈𝑡−𝜏,𝑡}𝑥(𝜉),∀𝑡≥𝑡₀. (2.17) On the other hand, we have 𝑥(𝑡)=𝑦(𝑡)−𝑦(𝑡)+[𝑧(𝑡)−𝑧(𝑡)]²(𝑡)+𝑧(𝑡(𝑡)≤2𝑦)≤2𝑦+2𝜂²𝜏²max_{[]𝜉∈𝑡−𝜏,𝑡}𝑥(𝜉),∀𝑡≥𝑡₀, (2.18) 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 Theorem 2.3, we present several examples as follows.

Example 2.4. Consider the 𝑑-dimensional system of linear functional-integrodifferential equation 𝑑𝑥(𝑡)−𝑑𝑡^𝑡_𝑡−𝜏𝑁(𝑡,𝜉)𝑥(𝜉)𝑑𝜉=𝐿(𝑡)𝑥(𝑡)+𝑀(𝑡)𝑥(𝑡−𝜏)+𝐺(𝑡),𝑡≥𝑡₀,𝑥(𝑡)=𝜑(𝑡),𝑡₀−𝜏≤𝑡≤𝑡₀, (2.19) where 𝐺∶[𝑡₀,+∞)→ℂ^𝑑 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 𝑡≥𝑡₀ and (𝑡,𝜉)∈𝔻, 𝜇(𝐿(𝑡()))≤𝜇≤−𝜄+𝜅2,‖𝐿(𝑡)‖≤𝜄,‖𝑀(𝑡)‖≤𝜅,𝑁(𝑡,𝜉)≤𝜚, (2.20) 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𝜚²𝜏²<1,4𝜅−𝜇𝜚²𝜏²()<−𝜇+𝜄+𝜅21−2𝜚²𝜏². (2.21)

Example 2.5. Consider the system of partial functional-integrodifferential equations 𝜕𝑢()−1𝜕𝑡𝑣,𝑡4^𝑡_{√𝑡−3/4}()𝑢()1exp𝜉−𝑡𝑣,𝜉𝑑𝜉=−2𝑢()+𝑢√𝑣,𝑡𝑣,𝑡−3/4241+𝑢²𝑣,𝑡−√()()(),𝑢()()−√3/4+𝑔𝑣,𝑡,𝑡≥0,𝑣∈0,2𝜋,𝑢𝑣,𝑡𝑣,𝑡=sin𝑣exp−𝑣𝑡,𝑡∈34[],𝑢()(),0,𝑣∈0,2𝜋0,𝑡=𝑢2𝜋,𝑡=0,𝑡≥0, (2.22) where 𝑔(𝑣,𝑡) is a continuous function chosen such that this system has the exact solution 𝑢(𝑣,𝑡)=sin𝑣exp(−𝑣t). By discretizing the spatial variable 𝑣 by a uniform mesh 𝑣_𝑖=𝑖Δ𝑣(𝑖=0,1,…,𝑙,Δ𝑣=2𝜋/𝑙), the system (2.22) can be transformed into a system of ordinary functional-integrodifferential equations, 𝑑𝑥(𝑡)−1𝑑𝑡4^𝑡_{√𝑡−3/4}()𝑥(𝜉)1exp𝜉−𝑡𝑑𝜉=−2𝑥(𝑡)+124𝑥√𝑡−34(𝑡)𝑥(𝑡)=+𝐺,𝑡≥0,sin𝑣₁exp−𝑣₁𝑡,sin𝑣₂exp−𝑣₂𝑡,…,sin𝑣_𝑙−1exp−𝑣_𝑙−1𝑡^𝑇−√,𝑡∈34,0, (2.23) where 𝑥(𝑡)=𝑥₁(𝑡),𝑥₂(𝑡),…,𝑥_𝑙−1(𝑡)^𝑇,𝐺(𝑡)=𝑔𝑣₁,𝑡,𝑔𝑣₂,𝑡,…,𝑔𝑣_𝑙−1,𝑡^𝑇,𝑥(𝑡)=𝑥₁(𝑡)1+𝑥²₁(𝑡),𝑥₂(𝑡)1+𝑥²₂(𝑡)𝑥,…,^𝑙−1(𝑡)1+𝑥²_𝑙−1(𝑡)^𝑇,𝑥_𝑖(𝑡)≈𝑢𝑣_𝑖,𝑡. (2.24) 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 𝑠-stage Runge-Kutta method is described by the Butcher tableau𝑐𝐴𝑏^𝑇 (3.1) where𝐴∶=𝑎_𝑖𝑗∈ℝ^𝑠×𝑠,𝑏∶=𝑏₁,𝑏₂,…,𝑏_𝑠^𝑇∈ℝ^𝑠,^𝑠_𝑖=1𝑏_𝑖𝑐=1.𝑐=₁,𝑐₂,…,𝑐_𝑠^𝑇,0≤𝑐_𝑖≤1(𝑖=1,2,…,𝑠). (3.2)

A Runge-Kutta method (3.1) is called algebraically stable if𝐷∶=diag𝑏₁,𝑏₂,…,𝑏_𝑠≥0,𝑀∶=𝐷𝐴+𝐴^𝑇𝐷−𝑏𝑏^𝑇≥0, (3.3)

where the notation “≥” means that a matrix is nonnegative definite. It is said to be strictly stable at infinity if 𝑅(∞)∶=lim_𝑧→∞𝑅(𝑧) exists and satisfies |𝑅(∞)|<1, where𝑅(𝑧)=1+𝑧𝑏^𝑇𝐼_𝑠−𝑧𝐴⁻¹𝑒(𝑧∈ℂ),𝑒=(1,1,…,1)^𝑇∈ℝ^𝑠, (3.4)

and 𝐼_𝑠 denotes the 𝑠×𝑠 identity matrix; it is said DJ-irreducible if there is no nonempty index set 𝕃⊂{1,2,…,𝑠} such that𝑏_𝑗=0for𝑗∈𝕃,𝑎_𝑖𝑗=0for𝑖∉𝕃,𝑗∈𝕃, (3.5)

S-irreducible if there is no partition (𝕊₁,𝕊₂,…,𝕊_𝑟) of {1,2,…,𝑠} with 𝑟<𝑠 such that for all 𝑙 and 𝑚_{𝑘∈𝕊𝑚}𝑎_𝑖𝑘=_{𝑘∈𝕊𝑚}𝑎_𝑗𝑘,for𝑖,𝑗∈𝕊_𝑙, (3.6)

and irreducible if it is both DJ-irreducible and S-irreducible.

Proposition 3.1 (cf. [15]). A DJ-irreducible, algebraically stable Runge-Kutta methods satisfies 𝑏_𝑖>0  for all 𝑖.

Proposition 3.2 (cf. [14]). Assume that a Runge-Kutta method (3.1) with distinct 𝑐_𝑖 and positive 𝑏_𝑖 satisfies the simplifying condition 𝐵(2𝑠−2),𝐶(𝑠−1),𝐷(𝑠−1). Then this method is algebraically stable if and only if  |𝑅(∞)|≤1.

The class of extended Pouzet-Runge-Kutta methods for FIDEs (1.4) with underlying Runge-Kutta method (3.1) is given by𝑥_𝑖^(𝑛)−𝑧_𝑖^(𝑛)=𝑥_𝑛−𝑧_𝑛+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡_𝑗^(𝑛),𝑥_𝑗^(𝑛),𝑥_𝑗^{(𝑛−𝑚)}𝑥,𝑖=1,2,…,𝑠,_𝑛+1−𝑧_𝑛+1=𝑥_𝑛−𝑧_𝑛+ℎ^𝑠_𝑗=1𝑏_𝑗𝑓𝑡_𝑗^(𝑛),𝑥_𝑗^(𝑛),𝑥_𝑗^{(𝑛−𝑚)},𝑛≥0. (3.7) Here, the stepsize ℎ is chosen as ℎ=𝜏/𝑚(𝑚∈ℕ),  𝑡_𝑛∶=𝑡₀+𝑛ℎ,𝑡_𝑗^(𝑛)∶=𝑡_𝑛+𝑐_𝑗ℎ, 𝑥_𝑛, 𝑧_𝑛, 𝑥_𝑖^(𝑛) and 𝑧_𝑖^(𝑛) are approximations to 𝑥(𝑡_𝑛), 𝑧(𝑡_𝑛), 𝑥(𝑡_𝑖^(𝑛)) and 𝑧(𝑡_𝑖^(𝑛)), respectively, with 𝑧(𝑡) denoting the memory term𝑧(𝑡)∶=^𝑡_𝑡−𝜏𝑔(𝑡,𝜉,x(𝜉))𝑑𝜉. (3.8)

The integral approximations 𝑧_𝑛 and 𝑧_𝑗^(𝑛) are given by the Pouzet quadrature rules𝑧_𝑛=ℎ^𝑛−1𝑠_{𝑞=𝑛−𝑚𝑗=1}𝑏_𝑗𝑔𝑡_𝑛,𝑡_𝑗^(𝑞),𝑥_𝑗^(𝑞)𝑧,𝑛≥0,_𝑖(𝑛)=ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑔𝑡_𝑖^(𝑛),𝑡_𝑗^(𝑛),𝑥_𝑗^(𝑛)+ℎ^𝑛−1𝑠_{𝑞=𝑛−𝑚𝑗=1}𝑏_𝑗𝑔𝑡_𝑖^(𝑛),𝑡_𝑗^(𝑞),𝑥_𝑗^(𝑞)−ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑔𝑡_𝑖^(𝑛),𝑡_𝑗^{(𝑛−𝑚)},𝑥_𝑗^{(𝑛−𝑚)},𝑖=1,2,…,𝑠;𝑛≥0. (3.9)(3.10) Moreover, it is assumed that𝑥_𝑛=𝜑𝑡_𝑛,𝑥_𝑖^(𝑛)=𝜑𝑡_𝑖^(𝑛)for𝑡₀−𝜏≤𝑡_𝑛,𝑡_𝑖^(𝑛)≤𝑡₀. (3.11)

4. Some Basic Lemmas

In the subsequent numerical analysis we will imply the following notations:𝑥_𝑛∶=𝑥_𝑛−𝑥_𝑛,𝑧_𝑛∶=𝑧_𝑛−𝑧_𝑛,𝑦_𝑛∶=𝑥_𝑛−𝑧_𝑛,𝑦_𝑛∶=𝑥_𝑛−𝑧_𝑛,𝑦_𝑛∶=𝑦_𝑛−𝑦_𝑛,𝑥_𝑖^(𝑛)∶=𝑥_𝑖^(𝑛)−𝑥_𝑖^(𝑛),𝑧_𝑖^(𝑛)∶=𝑧_𝑖^(𝑛)−𝑧_𝑖^(𝑛),𝑦_𝑖^(𝑛)∶=𝑥_𝑖^(𝑛)−𝑧_𝑖^(𝑛),𝑦_𝑖^(𝑛)∶=𝑥_𝑖^(𝑛)−𝑧_𝑖^(𝑛),𝑦_𝑖^(𝑛)∶=𝑦_𝑖^(𝑛)−𝑦_𝑖^(𝑛),𝑓_𝑖^(𝑛)∶=𝑓𝑡_𝑖^(𝑛),𝑥_𝑖^(𝑛),𝑥_𝑖^{(𝑛−𝑚)}−𝑓𝑡_𝑖^(𝑛),𝑥_𝑖^(𝑛),𝑥_𝑖^{(𝑛−𝑚)}. (4.1)

Moreover, when the method (3.7)–(3.10) is applied to the problem (2.3), we set𝑧(𝑡)∶=^𝑡_𝑡−𝜏𝑔(𝑡,𝜉,𝑥(𝜉))𝑑𝜉 (4.2)

and denote the corresponding approximations of 𝑥(𝑡_𝑛),𝑧(𝑡_𝑛),𝑥(𝑡_𝑖^(𝑛)), and 𝑧(𝑡_𝑖^(𝑛)), by 𝑥_𝑛,𝑧_𝑛,𝑥_𝑖^(𝑛), and 𝑧_𝑖^(𝑛), respectively. Similarly, it is assumed that𝑥_𝑛=𝜓𝑡_𝑛,𝑥_𝑖^(𝑛)=𝜓𝑡_𝑖^(𝑛),for𝑡₀−𝜏≤𝑡_𝑛,𝑡_𝑖^(𝑛)≤𝑡₀. (4.3)

As mentioned in Section 1, 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 𝑦_𝑛+1²≤2𝛽𝜏+𝜂²𝜏²max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)²+2ℎ^𝑛𝑠_𝑖=0𝑗=1𝑏_𝑗𝛼+𝛽𝑥_𝑗^(𝑖)2+𝛾𝑧_𝑗^(𝑖)2,𝑛≥0. (4.4)

Proof. A straightforward computation and the assumption of algebraic stability lead to 𝑦_𝑛+1²−𝑦_𝑛²−2ℎ^𝑠_𝑗=1𝑏_𝑗ℜ𝑦_𝑗^(𝑛),𝑓_𝑗^(𝑛)=−ℎ^2𝑠𝑠_𝑖=1𝑗=1𝑚_𝑖𝑗𝑓_𝑖,𝑓_𝑗≤0,∀𝑛≥0, (4.5) where 𝑀=(𝑚_𝑖𝑗). Hence, it holds that 𝑦_𝑛+1²≤𝑦_𝑛²+2ℎ^𝑠_𝑗=1𝑏_𝑗ℜ𝑦_𝑗^(𝑛),𝑓_𝑗^(𝑛),∀𝑛≥0. (4.6) By (2.1), we further have for 𝑛≥0, 𝑦_𝑛+1²≤𝑦_𝑛²+2ℎ^𝑠_𝑗=1𝑏_𝑗𝛼𝑥_𝑗^(𝑛)2+𝛽𝑥_𝑗^{(𝑛−𝑚)2}+𝛾𝑧_𝑗^(𝑛)2. (4.7) An induction argument shows that the inequality (4.7) implies 𝑦_𝑛+1²≤𝑦₀²+2ℎ^𝑛𝑠_𝑖=0𝑗=1𝑏_𝑗𝛼𝑥_𝑗^(𝑖)2+𝛽𝑥_𝑗^{(𝑖−𝑚)2}+𝛾𝑧_𝑗^(𝑖)2,𝑛≥0, (4.8) in which 𝑦₀²=𝑥₀−𝑧₀−(𝑥₀−𝑧₀)²=^𝑡0_𝑡0−𝜏𝑔𝑡₀,𝜉,𝜑(𝜉)−𝑔𝑡₀,𝜉,𝜓(𝜉)𝑑𝜉²≤𝜂²𝜏²max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)². (4.9) Also, the assumptions 𝑚ℎ=𝜏 and ∑^𝑠_𝑗=1𝑏_𝑗=1 allow us to write ℎ^𝑛𝑠_𝑖=0𝑗=1𝑏_𝑗𝑥_𝑗^{(𝑖−𝑚)2}=ℎ^𝑠_𝑗=1𝑏_𝑗⁻¹_𝑖=−𝑚𝑥_𝑗^(𝑖)2+^𝑛−𝑚_𝑖=0𝑥_𝑗^(𝑖)2≤𝜏max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)²+ℎ^𝑛𝑠_𝑖=0𝑗=1𝑏_𝑗𝑥_𝑗^(𝑖)2,𝑛≥0. (4.10) 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 𝑧_𝑛²≤ℎ𝜂²𝜏^𝑠_𝑗=1𝑏_𝑗^{2𝑛−1𝑠}_𝑖=0𝑗=1𝑥_𝑗^(𝑖)2+𝜂²𝜏²𝑠^𝑠_𝑗=1𝑏_𝑗max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)²,𝑛≥0. (4.11)

Proof. By (2.2), 𝑚ℎ=𝜏 and the Cauchy-Schwarz inequality, we have 𝑧_𝑛²≤ℎ𝜂^𝑛−1𝑠_{𝑖=𝑛−𝑚𝑗=1}𝑏_𝑗𝑥_𝑗^(𝑖)2≤ℎ²𝜂^{2𝑛−1𝑠}_{𝑖=𝑛−𝑚𝑗=1}𝑏_𝑗^{2𝑛−1𝑠}_{𝑖=𝑛−𝑚𝑗=1}𝑥_𝑗^(𝑖)2=ℎ𝜂²𝜏^𝑠_𝑗=1𝑏_𝑗^{2𝑛−1𝑠}_{𝑖=𝑛−𝑚𝑗=1}𝑥_𝑗^(𝑖)2. (4.12) Also, it follows from (3.11), (4.3), and 𝑚ℎ=𝜏 that ℎ^𝑛−1𝑠_{𝑖=𝑛−𝑚𝑗=1}𝑥_𝑗^(𝑖)2≤𝜏𝑠max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)²ℎ,for𝑛=0,_{𝑛−1𝑠𝑖=𝑛−𝑚𝑗=1}𝑥_𝑗^(𝑖)2≤ℎ^−1𝑠_{𝑖=−𝑚𝑗=1}𝑥_𝑗^(𝑖)2+^𝑛−1𝑠_𝑖=0𝑗=1𝑥_𝑗^(𝑖)2≤𝜏𝑠max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)²+ℎ^𝑛−1𝑠_𝑖=0𝑗=1𝑥_𝑗^(𝑖)2ℎ,for1≤𝑛<𝑚,_{𝑛−1𝑠𝑖=𝑛−𝑚𝑗=1}𝑥_𝑗^(𝑖)2≤ℎ^𝑛−1𝑠_𝑖=0𝑗=1𝑥_𝑗^(𝑖)2,for𝑛≥𝑚. (4.13)(4.14)(4.15) A combination of (4.13)–(4.15) yields that for all 𝑛≥0, ℎ^𝑛−1𝑠_{𝑖=𝑛−𝑚𝑗=1}𝑥_𝑗^(𝑖)2≤𝜏𝑠max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)²+ℎ^𝑛−1𝑠_𝑖=0𝑗=1𝑥_𝑗^(𝑖)2. (4.16) 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 𝑧_𝑗^(𝑖)2≤3ℎ²𝜂^2𝑠_𝑟=1𝑎_𝑗𝑟^2𝑠_𝑟=1𝑥_𝑟^(𝑖)2+𝑚^𝑚_𝑞=1^𝑠_𝑟=1𝑏_𝑟^2𝑠_𝑟=1𝑥_𝑟^{(𝑖−𝑞)2}+^𝑠_𝑟=1𝑎_𝑗𝑟^2𝑠_𝑟=1𝑥_𝑟^{(𝑖−𝑚)2},𝑖≥0,𝑗=1,2,…,𝑠. (4.17)

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 𝑥_𝑛−𝑥_𝑛≤ℋmax_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉),∀𝑛≥1. (5.1) Moreover, the method (3.7)–(3.10) is called asymptotically stable for the problems of class 𝔽𝕀𝔻(𝛼,𝛽,𝛾,𝜂) if lim_𝑛→∞𝑥_𝑛−𝑥_𝑛=0. (5.2)
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 Theorem 5.7.

Theorem 5.2. Assume that the underlying Runge-Kutta method (3.1) with positive 𝑏_𝑖 is algebraically stable. Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is globally stable for the class 𝔽𝕀𝔻(𝛼,𝛽,𝛾,𝜂), with stability constant ℋ=2𝜏2𝛽+𝜂²𝜏(1+𝑠)+6𝑠𝛾𝜂²𝜏^2𝑠_j=1𝑏_𝑗𝑏_𝑗+^𝑠_𝑟=1𝑎_𝑗𝑟², (5.3) whenever 2+𝜂𝛼+𝛽²𝜏𝜎^𝑠_𝑗=1𝑏_𝑗1+6𝛾𝜏𝑏_𝑗+12𝛾𝜏^𝑠_𝑟=1𝑎_𝑗𝑟²≤0, (5.4) where 𝜎=min_{1≤𝑖≤𝑠}{𝑏_𝑖}>0.

Proof. It follows from Lemma 4.3 and ∑^𝑠_𝑗=1𝑏_𝑗=1(𝑏_𝑗>0) that ℎ^𝑛𝑠_𝑖=0𝑗=1𝑏_𝑗𝑧_𝑗^(𝑖)2≤3ℎ³𝜂^2𝑠_𝑗=1𝑏_𝑗^𝑠_𝑟=1𝑎_𝑗𝑟^2𝑛𝑠_𝑖=0𝑟=1𝑥_𝑟^(𝑖)2+𝑚^𝑠_𝑟=1𝑏²_𝑟^𝑠𝑛_𝑟=1^𝑚_𝑖=0𝑞=1𝑥_𝑟^{(𝑖−𝑞)2}+^𝑠_𝑗=1𝑏_𝑗^𝑠_𝑟=1𝑎_𝑗𝑟^2𝑛𝑠_𝑖=0𝑟=1𝑥_𝑟^{(𝑖−𝑚)2}. (5.5) Also, we have that ^𝑛𝑠_𝑖=0𝑟=1𝑥_𝑟^{(𝑖−𝑚)2}=^{𝑛−𝑚𝑠}_{𝑖=−𝑚𝑟=1}𝑥_𝑟^(𝑖)2≤^𝑛𝑠_𝑖=0𝑟=1𝑥_𝑟^(𝑖)2+𝑚^𝑠_𝑟=1max_{−𝑚≤𝑖≤−1}𝑥_𝑟^(𝑖)2 (5.6) and, by the inequality (3.13) in [1], ^𝑛𝑚_𝑖=0𝑞=1𝑥_𝑟^{(𝑖−𝑞)2}≤𝑚^𝑛_𝑖=0𝑥_𝑟^(𝑖)2+𝑚(𝑚+1)2max_{−𝑚≤𝑖≤−1}𝑥_𝑟^(𝑖)2. (5.7) Inserting both (5.6) and (5.7) into (5.5) yields ℎ^𝑛𝑠_𝑖=0𝑗=1𝑏_𝑗𝑧_𝑗^(𝑖)2≤3ℎ³𝜂^2𝑠_𝑗=1𝑏_𝑗2^𝑠_𝑟=1𝑎_𝑗𝑟²+𝑚²𝑏_𝑗^𝑛𝑠_𝑖=0𝑟=1𝑥_𝑟^(𝑖)2+𝑚^𝑠_𝑗=1𝑏_𝑗𝑚(𝑚+1)2𝑏_𝑗+^𝑠_𝑟=1𝑎_𝑗𝑟^2𝑠_𝑟=1max_{−𝑚≤𝑖≤−1}𝑥_𝑟^(𝑖)2. (5.8) With ℎ=𝜏/𝑚≤𝜏, (3.11) and (4.3), the left side of (5.8) can be bounded by 3𝜂²𝜏²ℎ𝜎^𝑠_𝑗=1𝑏_𝑗2^𝑠_𝑟=1𝑎_𝑗𝑟²+𝑏_𝑗^𝑛𝑠_𝑖=0𝑟=1𝑏_𝑟𝑥_𝑟^(𝑖)2+𝜏𝑠^𝑠_𝑗=1𝑏_𝑗𝑏_𝑗+^𝑠_𝑟=1𝑎_𝑗𝑟²max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)². (5.9) Substituting this bound into (4.4) yields 𝑦_𝑛+1²≤𝜏2𝛽+𝜂²𝜏+6𝑠𝛾𝜂²𝜏^2𝑠_𝑗=1𝑏_𝑗𝑏_𝑗+^𝑠_𝑟=1𝑎_𝑗𝑟²max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)²+2ℎ𝛼+𝛽+3𝛾𝜂²𝜏²𝜎^𝑠_𝑗=1𝑏_𝑗𝑏_𝑗+2^𝑠_𝑟=1𝑎_𝑗𝑟^2𝑛𝑠_𝑖=0𝑟=1𝑏_𝑟𝑥_𝑟^(𝑖)2,𝑛≥0, (5.10) which, together with (4.11), (5.10), and ∑^𝑠_𝑗=1𝑏_𝑗=1(𝑏_𝑗>0), implies that for 𝑛≥1, 𝑥_𝑛²=𝑦_𝑛+𝑧_𝑛²≤2𝑦_𝑛²+𝑧_𝑛²≤2𝜏2𝛽+𝜂²𝜏(1+𝑠)+6𝑠𝛾𝜂²𝜏^2𝑠_𝑗=1𝑏_𝑗𝑏_𝑗+^𝑠_𝑟=1𝑎_𝑗𝑟²max_{𝑡0−𝜏≤𝜉≤𝑡0}𝜑(𝜉)−𝜓(𝜉)²2+𝜂+2ℎ𝛼+𝛽²𝜏𝜎^𝑠_𝑗=1𝑏_𝑗1+6𝛾𝜏𝑏_𝑗+12𝛾𝜏^𝑠_𝑟=1𝑎_𝑗𝑟^{2𝑛−1𝑠}_𝑖=0𝑟=1𝑏_𝑟𝑥_𝑟^(𝑖)2. (5.11) Observing the condition (5.4) to (5.11), we readily derive the global stability conclusion. Hence, the theorem is proven.

Applying Proposition 3.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 distinct 𝑐_𝑖 and positive 𝑏_𝑖, satisfies the simplifying condition 𝐵(2𝑠−2),𝐶(𝑠−1),𝐷(𝑠−1), and |𝑅(∞)|≤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 𝐴, distinct 𝑐_𝑖, and 𝑏_𝑖>0 for all 𝑖. 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)

will subsequently be used, and we define a vector norm on the space ℂ^𝑑𝑠 by‖𝑈‖=^𝑠_𝑟=1𝑢_𝑟²,∀𝑈=𝑢₁^𝑇,𝑢₂^𝑇,…,𝑢_𝑠^𝑇𝑇∈ℂ^𝑑𝑠. (5.13)

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 𝑋_𝑛−𝑍_𝑛=𝑒⊗𝐼_𝑑𝑥_𝑛−𝑧_𝑛+ℎ𝐴⊗𝐼_𝑑𝐹_𝑛,𝑥_𝑛+1−𝑧_𝑛+1=𝑥_𝑛−𝑧_𝑛+ℎ𝑏^𝑇⊗𝐼_𝑑𝐹_𝑛,𝑛≥0. (5.14) When the matrix 𝐴 is invertible, it follows from the first equation of (5.14) that ℎ𝐹_𝑛=𝐴⁻¹⊗𝐼_𝑑𝑋_𝑛−𝑍_𝑛−𝐴⁻¹𝑒⊗𝐼_𝑑𝑥_𝑛−𝑧_𝑛. (5.15) Insertion of (5.15) into the second equation of (5.14) yields 𝑥_𝑛+1−𝑧_𝑛+1=1−𝑏^𝑇𝐴⁻¹𝑒𝑥_𝑛−𝑧_𝑛+𝑏^𝑇𝐴⁻¹⊗𝐼_𝑑𝑋_𝑛−𝑍_𝑛. (5.16) When 𝐴 is singular, we adopt a technique, proposed by Hairer and Wanner [14]. It consists in replacing 𝐴 by the regular matrix 𝐴+𝜀𝐼_𝑠 everywhere, followed by considering the limit lim_𝜀→0(𝐴+𝜀𝐼_𝑠)⁻¹, whose existence is assured by irreducibility and algebraic stability of the method (see Lemma  3.2 in [14]). Thus, we set 𝒜⁻¹𝐴∶=⁻¹,if𝐴isinvertible,lim_𝜖→0𝐴+𝜖𝐼_𝑠⁻¹,if𝐴issingular. (5.17) Hence, it holds that 𝑥_𝑛+1−𝑧_𝑛+1=𝑅(∞)𝑥_𝑛−𝑧_𝑛+𝑏^𝑇𝒜⁻¹⊗𝐼_𝑑𝑋_𝑛−𝑍_𝑛, (5.18) where 𝑅(∞)=1−𝑏^𝑇𝒜⁻¹𝑒. Applying Lemma 4.6 in [2] to (5.18) we obtain lim_𝑛→∞𝑥_𝑛+1−𝑧_𝑛+1=0, (5.19) if and only if |𝑅(∞)|<1,lim_𝑛→∞𝑋_𝑛−𝑍_𝑛=0. (5.20)
Next, we prove (5.19). We need only to show thatlim_𝑛→∞𝑋_𝑛−𝑍_𝑛=0, (5.21) since |𝑅(∞)|<1 is a known condition. In fact, by Proposition 3.1 we know that 𝑏_𝑖>0 for all 𝑖. According to this and (5.4), we derive from (5.11) that lim_𝑛→∞𝑥_𝑟^(𝑛)=0,𝑟=1,2,…,𝑠. (5.22) This implies that lim_𝑛→∞𝑋_𝑛=lim_𝑛→∞^𝑠_𝑟=1𝑥_𝑟^(𝑛)2=0. (5.23) Combining (4.17) and (5.22) yields lim_𝑛→∞‖𝑧_𝑗^(𝑛)‖=0 for all 𝑗, which leads to lim_𝑛→∞𝑍_𝑛=lim_𝑛→∞^𝑠_𝑗=1𝑧_𝑗^(𝑛)2=0. (5.24) A combination of (5.23) and (5.24) gives (5.21). Hence, (5.19) is true.
Finally, we prove that lim_𝑛→∞‖𝑥_𝑛‖=0. By (3.9), (2.2), and (5.22), we have𝑧_𝑛≤ℎ𝜂^𝑛−1𝑠_{𝑞=𝑛−𝑚𝑗=1}𝑏_𝑗𝑥_𝑗^(𝑞)⟶0,as𝑛⟶∞, (5.25) which means that lim_𝑛→∞‖𝑧_𝑛‖=0. This, together with (5.19), implies that 𝑥_𝑛≤𝑥_𝑛−𝑧_𝑛+𝑧_𝑛⟶0,as𝑛⟶∞. (5.26) 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 matrix 𝐴 is singular. Hence, for Pouzet-Runge-Kutta methods whose underlying Runge-Kutta has an invertible the matrix 𝐴, 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 𝐴 and positive 𝑏_𝑖 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 Theorem 5.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 𝑐_𝑖 is irreducible, algebraically stable, and satisfies |𝑅(∞)|≠1 and the simplifying condition 𝐵(2𝑠−2),𝐶(𝑠−1),𝐷(𝑠−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 𝑐_𝑖 and positive 𝑏_𝑖 is irreducible, strictly stable at infinity, and satisfies the simplifying condition 𝐵(2𝑠−2),𝐶(𝑠−1),𝐷(𝑠−1). Then the corresponding Pouzet-Runge-Kutta method (3.7)–(3.10) is asymptotically stable for the class 𝔽𝔻𝕀(𝛼,𝛽,𝛾,𝜂) whenever (5.4) holds.

Similar to Theorem 5.8, the irreducibility condition can be dropped in Corollary 5.10 when the matrix 𝐴 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 𝑅(∞)=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,4√3]. These methods produce a series of high-precision numerical solutions for the (2.22) on [0,4√3;0,2𝜋].

Let+𝜂Γ∶=2𝛼+𝛽²𝜏min_{1≤𝑖≤𝑠}𝑏_𝑖^𝑠_𝑗=1𝑏_𝑗1+6𝛾𝜏𝑏_𝑗+12𝛾𝜏^𝑠_𝑟=1𝑎_𝑗𝑟². (6.1)

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).

The 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:()=𝜋Δ𝑣,ℎ,√303,𝜋16,√603,𝜋32,√1203,𝜋64,√2403,128 (6.2) and then applying the above Pouzet-Runge-Kutta schemes to the system (2.23) on [0,4√3], 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 useerr∶=max_{1≤𝑛≤𝑁}𝑥_𝑛−𝑢𝑡_𝑛∞ (6.3)

Figure 1 

Numerical solution of (2.22) by 2-stage Gauss-type method.

to characterize the errors of the methods, where𝑢𝑡_𝑛∶=sin𝑣₁exp−𝑣₁𝑡_𝑛,sin𝑣₂exp−𝑣₂𝑡_𝑛,…,sin𝑣_𝑙−1exp−𝑣_𝑙−1𝑡_𝑛^𝑇 (6.4)

is the vector whose entries consist of the system (2.22) true solutions at the meshpoints (𝑣_𝑖,𝑡_𝑛), 𝑖=1,2,…,𝑙−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,𝑥_𝑖^(𝑛)=𝑥_𝑛+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝐹𝑡_𝑗^(𝑛),𝑥_𝑗^(𝑛),𝑥_𝑗^{(𝑛−𝑚)},𝑧_𝑗^(𝑛)𝑥,𝑖=1,2,…,𝑠,_𝑛+1=𝑥_𝑛+ℎ^𝑠_𝑗=1𝑏_𝑗𝐹𝑡_𝑗^(𝑛),𝑥_𝑗^(𝑛),𝑥_𝑗^{(𝑛−𝑚)},𝑧_𝑗^(𝑛),𝑛≥0, (7.1) where the meanings of the notations are similar to those indicated in method (3.7) and {𝑧_𝑗^(𝑛)}^𝑠_𝑗=1 are determined by (3.10). But such methods cannot be applied directly to systems (1.4) unless the integrated function 𝑔 is continuously differentiable on its domain. When the latter holds, the system (1.4) can be transformed into𝑥^′(𝑡)=𝐹𝑡,𝑥(𝑡),𝑥(𝑡−𝜏),^𝑡_𝑡−𝜏𝑔_𝑡(𝑡,𝜉,𝑥(𝜉))𝑑𝜉,𝑡≥𝑡₀,𝑥(𝑡)=𝜑(𝑡),𝑡₀−𝜏≤𝑡≤𝑡₀, (7.2) where𝐹𝑡,𝑥,𝑦,𝑧=𝑓𝑡,𝑥,𝑦+𝑔(𝑡,𝑡,𝑥)−𝑔𝑡,𝑡−𝜏,𝑦+𝑧. (7.3)

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 𝑔 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,4√3], 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).

Acknowledgments

This work is supported by NSFC (Grant no. 11171125, 91130003) and NSFH (Grant no. 2011CDB289), and its initial version was completed at Hong Kong Baptist University when Co Zhang visited Professor Hermann Brunner. The author is grateful to Hermann Brunner for his valuable discussion and aborative correction to the paper.