Abstract

Based on the linear multistep methods for ordinary differential equations (ODEs) and the canonical interpolation theory that was presented by Shoufu Li who is exactly the second author of this paper, we propose the linear multistep methods for general Volterra functional differential equations (VFDEs) and build the classical stability, consistency, and convergence theories of the methods. The methods and theories presented in this paper are applicable to nonneutral, nonstiff, and nonlinear initial value problems in ODEs, Volterra delay differential equations (VDDEs), Volterra integro-differential equations (VIDEs), Volterra delay integro-differential equations (VDIDEs), etc. At last, some numerical experiments verify the correctness of our theories.

1. Introduction

VFDEs contain many subtypes, such as VDDEs, VIDEs, VDIDEs, etc. Certainly ODEs are also a subtype of them. VFDEs are widely applied in many fields of science and technology (see [15] and their references), for which there have been remarkable theoretical and numerical analysis research results. As for VDDEs, please refer to [619], as for VIDEs, please refer to [2023], and as for VDIDEs, please refer to [2429]. In recent decades, Li [3035] has carried on systematic research for stiff general VFDEs and the numerical methods for them. In 2014, Li [36] established the classical stability and convergence theories of Runge–Kutta methods for nonstiff, nonlinear general VFDEs. It is well known that in solving ODEs, linear multistep methods have significant advantages in the aspects of format simplification and computation cost, and experts have presented many famous linear multistep methods such as Backward Differentiation Formula (BDF) methods and Adams methods. Based on the linear multistep methods for ODEs and the canonical interpolation theory presented by Li [34], we propose the linear multistep methods for general VFDEs and build the classical stability, consistency, and convergence theories of the methods. The methods and theories presented in this paper are applicable to nonneutral, nonstiff, and nonlinear initial value problems in ODEs, VDDEs, VIDEs, VDIDEs, etc., and are interesting companions to the methods and theories in [36]. Furthermore, the paper may have some value for the study of numerical methods for fractional order VFDEs because many numerical methods for fractional-order differential equations are presented by extending the methods for integer-order differential equations.

This paper is organized as follows. In Section 2, the linear multistep methods for general VFDEs are introduced. The classical stability of the linear multistep methods is discussed in Section 3. The classical consistency and convergence analyses are carried out in Section 4. Some numerical experiments are carried out in Section 5, which verifies the correctness of the theories presented in this paper.

2. Derivation of the Numerical Methods

Let be the -dimensional Euclidean space with standard inner product and the corresponding norm . For any given closed interval , let denote the Banach space consisting of all continuous mappings with the maximum norm .

Consider VFDEs of the form [3032, 34, 36]where and are constants, is an initial function, and is a given continuous mapping which satisfies the classical Lipschitz conditions [36]:where and are classical Lipschitz constants. Equations (1) has a unique solution [37], and we further make the same assumption as [36] that the system in (1) is stable for the initial function , which is to say there exists a function which satisfies has only moderate size such thatwhere is the solution of the perturbed problem:

We use to denote the problem class consisting of all the equations of form (1) which satisfy condition (2) and all the other assumption conditions made above. The problem class contains nonneutral initial value problems in ODEs, DDEs, IDEs, DIDEs, etc. [36].

Combine the linear -step method:for ODEs with piecewise Lagrangian interpolation operators which are constructed based on canonical interpolation theory [34, 36], and we get the linear -step method:for the VFDEs (1), where , all and are real constants, , , is called the characteristic polynomial of the method (6), is the stepsize, is a given appropriate positive integer, are grid points, is an approximation to the initial function , , is an approximation to the true value , and the piecewise interpolation operators : based on -degree piecewise Lagrangian interpolation polynomials are defined as follows: if ,and if ,withwhere the integer can be freely determined under the conditions

According to [34, 36], we know a piecewise interpolation operator defined by (7) or by (8) and (9) which satisfies the canonical conditionwhere the constant only depends on the polynomial order and does not depend on and . In this paper, we always assume that is of moderate size. In fact, Li [34, 36] has given the computational formula of and found , respectively, equals , and 1.631 when , respectively, equals , and 3. We call method (5) as the mother method of method (6).

If in method (6), in order to reduce the amount of calculation, we defineand now the method is an explicit linear -step method which is a special case of (6).

In particular, we point out that when the linear -step method (6) for VFDEs starts, on the one hand, like method (5) for ODEs, it needs starting values, and on the other hand, it may need several grid points as interpolation nodes to construct interpolation operator . Based on the above statements, in this paper, we always assume the values , which are called additional starting values, should be provided in advance by other ways, and method (6) starts with . By the Banach contraction mapping principle, we conclude that, for any given starting data , method (6) can uniquely determine the sequence when , where can be any number under the condition

3. Stability Analysis

Firstly, we list some concepts and elementary facts about scalar linear difference equations as follows [35, 38, 39].

Consider the th order scalar linear difference equation with constant coefficients:here and later, is an integer constant, and where is an integer constant, every is a complex constant, , any is a complex number depending on . If a complex number sequence , satisfies equation (14), it is denoted by , in this paper, and is called a solution of (14). It is easily known that equation (14) has a unique solution when given initial values . If every is zero, equation (14) is a so-called homogeneous equation:

A system of linearly independent solutions of equation (15) is called a fundamental system of (15), and the general solution of equation (15) can be written bywhere each is an arbitrary complex constant.

Call the algebraic equationwhich corresponds to (14) or (15) the characteristic equation.

Proposition 1. (see [35, 38]). Suppose are all different roots of the characteristic equation (17), and let denote the multiple number of . Then, when , the sequences with the elements which can be equal toform a fundamental system of equation (15), and when , which means equation (17) has zero root and we further suppose is a -tuple root of the characteristic equation (17), the sequences with the elements which can be equal toform a fundamental system of equation (15), where is the Kronecker function:

Remark 1. Some explanations about the proof of Proposition 1 are as follows:(i)When and , P. Henrici has given the proof in the famous book [38].(ii)When and , from the conclusion in the case that and , the proof is easily completed.(iii)When and , the corresponding part of Proposition 1 is based on the proof [38] in the case that and mentioned in the above (i) and some conclusions in the case that and in [35], and we offer the explanations about the proof as follows. Firstly, is a -tuple root of the characteristic equation (17) which means that , and it is easily further known as the sequences with the elements are the solutions of equation (15). Secondly, based on the proof presented by P. Henrici mentioned in the above (i), we can know the sequences with the elements shown in (20) are solutions of equation (15) and that the solutions with the elements shown in (19) and (20) are linearly independent.(iv)When and , from the conclusion in the case that and , the proof is easily completed.

Proposition 2. (see [35]). When the initial values , the solution of equation (14) is given bywhere each satisfies the homogeneous equation:

Remark 2. In this paper, when , we define and , where each is any a given expression.

Remark 3. By Theorem 5.2 in [38], Proposition 2 is easily shown.

Proposition 3. (see [35]). Suppose is a fundamental system of equation (15) and is the solution of the homogeneous equation (23). Then, the general solution of equation (14) can be expressed bywhere each is an arbitrary constant.

Proposition 4. (discrete Bellman inequality). Let , be a sequence of nonnegative real numbers which satisfyThen,

Definition 1. We say method (6) is zero-stable if the sequence determined by method (6) applied to equation (1) with starting data (the method (6) starts with ) and the sequence determined by the perturbation equationswith starting data , satisfywhere and are both perturbations, , the constant depends on method (6), and , and the constant is independent of .

Definition 2. We say method (6) satisfies the root condition if for all roots of the characteristic polynomial , the moduli of them are all no larger than 1, and among them any root whose modulus equals 1 is simple root.

Theorem 1. Method (6) is zero-stable if and only if it satisfies the root condition.

Proof. Firstly, we prove method (6) satisfies the root condition if it is zero-stable. Consider the ODEs:which are a special case of (1), where , , and the continuous mapping satisfies the classical Lipschitz condition:where is a classical Lipschitz constant. For the ODEs (29), method (6) degrades into its mother method (5). Method (6) is zero-stable implies that its mother method (5) is zero-stable, which means method (6) should satisfy the root condition if it is zero-stable [35, 40].
We assume that method (6) satisfies the root condition, and now we prove it is zero-stable.
Denote . From (6) and (27), we obtainwhereUsing Proposition 1 and considering that all αi are real constants, we can find a real fundamental system [38] of the following scalar homogeneous equation:and by applying Proposition 3 to (31), we obtainwhere each is a constant vector and each satisfies the scalar homogeneous equation:From the assumption that method (6) satisfies the root condition, we can conclude for all and in (34), there exists a constant which is only dependent on the mother method (5) such thatIt follows from (34) and (36) thatSetting in (34), we thus findwhere the matrixBy (38), we knowwhere . Combining (37) with (40), we haveIt is known from (2), (6), (11), (27), and (32) thatThen, (42) yieldsSetting and noticing (43), we haveDenoteCombining (44) with (45), we concludeFrom inequalities (41) and (46), noticing Remark 2, we obtainBy (45), we haveIt can be concluded from (47) and (48) thatLetFrom (49) and (50), we know that, for ,For any given , setwhich depends on method (6), and , where is defined in Section 2. Noticing , when , from (52), we obtainBecause does not decrease as increases, when and , from (54), we knowIt is concluded from (54) and (55) thatwhich impliesFrom (57) and Proposition 4, we havesowherewhich is independent of . This completes the proof.

4. Consistency and Convergence Analyses

In this section, for equation (1), we assume the mapping possesses sufficiently high-order continuous partial derivatives, and its true solution possesses sufficiently high-order continuous derivatives on the interval and the constants used later which are defined asare all of moderate size except for some special cases such as in the transient phase of a stiff problem.

Definition 3. It is said the piecewise Lagrangian interpolation operator in (6) is consistent of order if, for any given function which is sufficiently differentiable on the subinterval , there existswhere the function is a restriction of the function on the subinterval and the constant only depends on some defined as .

Remark 4. By Definition 3 and the Lagrangian interpolation theorem, we know the piecewise interpolation operators in (6) which based on -degree Lagrangian interpolation polynomial are all consistent of order .

Definition 4. It is said method (6) is consistent of order if its mother method (5) and the piecewise interpolation operators are all consistent of order .

Remark 5. As for the consistency of the mother method (5), please see references [35, 38, 40].

Definition 5. Method (6) is said to be convergent of order if, for the sequence determined by method (6) applied to any given equations (1) in with the starting data (the method (6) starts with ), there exists a sufficiently small positive number , such thatwhere the continuous functions and depend on , method (6), , , and some .

Theorem 2. If method (6) is consistent of order and satisfies the root condition, it is convergent of order .

Proof. Assume method (6) is consistent of order and satisfies the root condition. In this proving course, we always assume , where is defined by (53). Using method (6) to solve any given equations (1) , we easily know the true solution obviously satisfies the perturbation equations:whereFrom (1) and (67), we concludeSince the mother method (5) is consistent of order , equation (68) implies [35, 38, 40]where the constant only depends on the mother method (5) and defined above. Because all the interpolation operators are consistent of order , by Definition 3, we obtainwhere the constant only depends on some defined above. By Theorem 1 and assumption that method (6) satisfies the root condition, we know method (6) is zero-stable, which tells us that inequality (59) holds. From (6), (27), (50), (59), (65), (69), and (70), we further conclude thatwhereand defined in the course of proving Theorem 1 in Section 3. It can be seen that the continuous functions and depend on , method (6), , , and some . This completes the proof.

Remark 6. In this paper, the strict condition (2) can be generally weakened as the inequalities in (2) hold only in some neighborhood of the true solution of the equations [35, 36, 38, 41].

5. Numerical Examples

In this section, the linear multistep methods formed as (6) are used to solve VFDEs formed as (1) , and for convenience, we give the starting data according to true solution, i.e., let

For a method formed as (6) whose stepsize is , we denote the maximum global error:where and are the th components of the true solution and the approximation solution , respectively, and furthermore, according to [36, 41], the observing order of the numerical method is defined as

In order to adequately verify the correctness of the theory presented in this paper, we use various different zero-stable mother methods including the Adams–Bashforth methods [38, 42], the Adams–Moulton method [38, 43], the Backward Differentiation Formula (BDF) method, and the Milne-Simpson method [44, 45] with corresponding piecewise Lagrangian interpolation operator to build the methods formed as (6), where each mother method and its corresponding interpolation operator are both consistent of the same order . According to Definition 4, these methods are consistent of order , and by Theorem 2, we know they are convergent of order . For convenience, these methods are named the same as their corresponding mother methods. For simplicity, let AB denote the Adams–Bashforth method which is consistent of order , let AM denote the Adams–Moulton method which is consistent of order , let B denote the BDF method which is consistent of order , and let MS4 denote the Milne-Simpson method which is consistent of order 4.

Remark 7. We list the mother methods used in this section as follows:(i)The Adams–Bashforth method which is consistent of order 2 takes the form(ii)The Adams–Bashforth method which is consistent of order 3 takes the form(iii)The Adams–Moulton method which is consistent of order 2 takes the form(iv)The BDF method which is consistent of order 2 takes the form(v)The Milne-Simpson method takes the form

Example 1. Consider the pantograph equations [36]:whose unique true solution is . Equation (81) is a special case , where [36]. We use the methods AB2, AM2, and B2 with to solve (81), respectively. The maximum global errors and the observing orders are listed in Table 1.
From Table 1, we can see for the pantograph equation (81) the observing orders of the three methods AB2, AM2, and B2 are all very close to the theoretical convergent order 2.

Example 2. Consider the nonlinear VIDEs [46]:The unique true solution of equation (82) is . Letwhere , . Then, equation (82) can be written in form (1), where and are both of moderate size. We use the method AB2 with to solve (82), respectively. The maximum global errors and the observing orders are listed in Table 2.
From Table 2, we can see, for the nonlinear VIDEs (82), the observing orders of the method AB2 are all very close to the theoretical convergent order 2.

Example 3. Consider the nonlinear VDIDEs [36]:The unique true solution of equation (84) is . Equation (84) in the sense of Remark 6, where and are both of moderate size [36]. We use MS4 with to solve (84), respectively, and list the maximum global errors and the observing orders in Table 3.
From Table 3, we can see, for the nonlinear VDIDEs (84), the observing orders of the method MS4 are all very close to the theoretical convergent order 4.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare there are no conflicts of interest.

Acknowledgments

This work was supported by Hunan Provincial Natural Science Foundation of China (Grant no. 2017JJ2194), the Scientific Research Fund of Hunan Provincial Education Department (Grant no. 20C1261), and the PhD Start-Up Fund from Hunan University of Arts and Science (Grant no. 16BSQD24).