Estimation of Longest Stability Interval for a Kind of Explicit Linear Multistep Methods

Abstract

The new explicit linear three-order four-step methods with longest interval of absolute stability are proposed. Some numerical experiments are made for comparing different kinds of linear multistep methods. It is shown that the stability intervals of proposed methods can be longer than that of known explicit linear multistep methods.

1. Introduction

For the initial value problem of the ordinary differential equation (ODE)𝑦^′(𝑡)=𝑓(𝑡,𝑦(𝑡)), (1.1) where 𝑓∶[𝑡₀,𝑡_end]×ℛ^𝑚→ℛ^𝑚 and 𝑦(𝑡₀)=𝑦₀ with 𝑦₀∈ℛ^𝑚, there are a lot of numerical methods to be proposed for the numerical integration. Among them, linear multistep methods (LMMs) are a class of the most prominent and most widely used methods, see [1, 2] and the references therein.

Adams methods are among the oldest of LMMs, dating back to the nineteenth century. Nevertheless, they continue to play a key role in efficient modern algorithms. The first to use such a method was Adams in solving a problem of Bashforth in connection with capillary action, see [3]. In contrast to one-step methods, where the numerical solution is obtained solely from the differential equation and the initial value 𝑦₀, a linear multistep (𝑘-step) method requires (𝑘) starting values 𝑦₀,𝑦₁,…,𝑦_𝑘−1 and a multistep (𝑘-step) formula to obtain an approximation to the exact solution, see [4].

So far as we know, explicit linear multistep methods (ELMMs) have some advantages such as simple calculation formulae, and small error constants. However, due to the famous Dahlquist barrier in [5], an explicit linear multistep method cannot be A-stable. Therefore, we try to find the new explicit linear multistep methods with the longest interval of stability region in this paper. And some numerical experiments are given to compare the proposed methods with existing methods such as Adams-Bashforth method, Adams-Moulton methods, and BDF methods. Practical calculations have shown that these proposed methods are adaptive.

2. Linear Multistep Methods

Applying the linear multistep (𝑘-step) methods to the initial value problem (1.1), we obtain the recurrence relation^𝑘_𝑗=0𝛼_𝑗𝑦_𝑛+𝑗=ℎ^𝑘_𝑗=0𝛽_𝑗𝑓_𝑛+𝑗, (2.1) where 𝑦_𝑛 denotes an approximation to the solution 𝑦(𝑡_𝑛), 𝑓_𝑛+𝑗∶=𝑓(𝑡_𝑛+𝑗,𝑦_𝑛+𝑗), 𝑡_𝑛+𝑗=𝑡_𝑛+𝑗ℎ for 𝑗=0,1,…,𝑘, the constant steplength ℎ>0, and 𝑘 starting conditions 𝑦₀,𝑦₁,…,𝑦_𝑘−1 are required. Here, 𝛼_𝑗 and 𝛽_𝑗(𝑗=0,1,…,𝑘) are constants subject to the condition 𝛼_𝑘=1,|𝛼₀|+|𝛽₀|≠0. If 𝛽_𝑘=0, then the corresponding methods (2.1) are explicit, and implicit otherwise.

Then, we define the first and second generating polynomials by𝜌(𝜁)∶=^𝑘_𝑗=0𝛼_𝑗𝜁^𝑗,𝜎(𝜁)∶=^𝑘_𝑗=0𝛽_𝑗𝜁^𝑗, (2.2) where 𝜁∈𝒞 is a dummy variable.

Consider the scalar test equation𝑦^′(𝑡)=𝜆𝑦(𝑡), (2.3) where 𝜆∈𝒞 and ℜ𝜆<0. Its characteristic polynomial can be written asΠ(𝜉;ℎ)∶=𝜌(𝜉)−ℎ𝜎(𝜉), (2.4) where ℎ=𝜆ℎ.

Here, we quote some important definitions (see Sections 3.2, 3.8, and 6.3 in the reference [2]).

Definition 2.1. The set Ω={ℎ∈𝒞∣allrootsofΠ(𝜉;ℎ)=0havemoduluslessthan1} (2.5) is called the region of absolute stability of the linear multistep method. The corresponding numerical method is said to be absolutely stable.

Definition 2.2. The intersection of Ω with the real axis in the complex ℎ-plane is called the interval of absolute stability.

Definition 2.3. The multistep method is said to be A-stable if Ω⊃𝒞⁻, where 𝒞⁻∶={𝑧∈𝒞∣ℜ𝑧<0}.

Definition 2.4. The linear multistep methods (2.1) are said to be of order 𝑝 if 𝐶₀=𝐶₁=⋯=𝐶_𝑝=0 and 𝐶_𝑝+1≠0, where 𝐶₀=^𝑘_𝑗=0𝛼_𝑗≡𝜌(1),𝐶₁=^𝑘_𝑗=0(𝑗𝛼_𝑗−𝛽_𝑗)≡𝜌^′𝐶(1)−𝜎(1),_𝑞=^𝑘_𝑗=0[1𝑗𝑞!^𝑞𝛼_𝑗−1𝑗(𝑞−1)!^𝑞−1𝛽_𝑗],𝑞=2,3,…. (2.6)

3. Boundary Locus Technique

Since an explicit linear multistep method cannot be A-stable, we focus our attention on its absolute stability. The most convenient method for finding regions of absolute stability is the boundary locus technique (BLT).

Let the contour 𝜕Ω in the complex ℎ-plane be defined by the requirement that for all ℎ∈𝜕Ω, one of the roots of Π(𝜉;ℎ) has modulus 1, that is, is of the form 𝜉=exp(𝑖𝜃). Thus, we have Π(exp(𝑖𝜃);ℎ)=𝜌(exp(𝑖𝜃))−ℎ𝜎(exp(𝑖𝜃))=0.

This identity is readily solved for ℎ, and we introduce the root locus curveℎ∶𝜃→ℎ(𝜃)=𝜌(exp(𝑖𝜃)),𝜎(exp(𝑖𝜃)) (3.1) which maps [0,2𝜋] onto Γ, where Γ∶={ℎ∈𝒞∣ there exsist 𝜉,|𝜉|=1,Π(𝜉;ℎ)=0}. If this map is one to one, then Γ is a Jordan curve, which can constitute the boundary of Ω, see [6].

Theorem 3.1. If 𝜎(−1)≠0, then two intersected points of curve Γ with the real axis in complex ℎ-plane are ℎ₁=0, ℎ₂=𝜌(−1)/𝜎(−1). Furthermore, if the real interval Δ=(𝜌(−1)𝜎(−1),0) (3.2) does not intersect with the curve Γ and there exists 𝜀₀∈Δ such that all roots of Π(𝜉;𝜀₀)=0 have modulus less than 1 in complex ℎ-plane, then Δ is just the interval of absolute stability of linear multistep methods.

Proof. This proof is referred to [7].

Although this criterion is not easy to prove, we always use condition (3.2) to estimate the maximal length of stability interval for certain linear multistep method.

Fortunately, the true stability interval is just the same as that estimated by (3.2) for most cases. Therefore, the longest interval of absolute stability can be evaluated by condition (3.2) for given stepnumber, order and other conditions.

For example, consider the two-order explicit linear two-step methods, whose first and second generating polynomials can be written as𝜌(𝜉)=²_𝑗=0𝛼_𝑗𝜉^𝑗,𝜎(𝜉)=¹_𝑗=0𝛽_𝑗𝜉^𝑗. (3.3)

Without loss of generality, we assume that −1<𝛼₀<1.

Corollary 3.2. The estimated interval of absolute stability of two-order explicit linear two-step methods (3.3) is Δ=(−1−𝛼₀,0), where −1<𝛼₀<1.

Corollary 3.3. The maximal length of stability interval of two-order explicit linear two-step methods (3.3) can not exceed 2.

Corollary 3.4. Let Δ_{𝐸;𝑝,𝑠} or Δ_{𝐼;𝑝,𝑠} be the estimated interval of absolute stability of 𝑝-order explicit or implicit linear 𝑠-step methods. Then(I)Δ_𝐸;3,3=(6(𝛼₀+𝛼₂)/(𝛼₀−𝛼₂+10),0), where −1<𝛼₀<1, 𝛼₀+𝛼₂<0 and 𝛼₀−𝛼₂<2;(II)Δ_𝐼;3,2=(−6(𝛼₀+1)/(1−𝛼₀),0), where −1<𝛼₀<1;(III)Δ_𝐼;4,3=(−3(𝛼₀+𝛼₂)/(1+𝛼₀),0), where −1<𝛼₀<1, 𝛼₀+𝛼₂<0 and 𝛼₀−𝛼₂<2.

In fact, we can also study the interval of absolute stability of linear multistep methods from definition directly.

For example, consider the two-order explicit linear two-step methods𝑦_𝑛+2−𝛼𝑦_𝑛+1−(1−𝛼)𝑦_𝑛=ℎ2[(4−𝛼)𝑓_𝑛+1−𝛼𝑓_𝑛], (3.4) where 0<𝛼<2.

Applying to the test equation (2.3), we have the stability polynomial Π(𝜉;ℎ)=𝜉²−𝛼𝜉−(1−𝛼)−ℎ2[(4−𝛼)𝜉−𝛼]. (3.5)

Define 𝜉=exp(𝑖𝜃), then the roots of above characteristic polynomial can be written as ℎ=ℜℎ+𝑖ℑℎ, where ℜℎ=(2𝛼/𝐷(𝜃))(4cos𝜃−cos2𝜃−3), ℑℎ=(2/𝐷(𝜃))[2(𝛼²−3𝛼+4)sin𝜃−𝛼sin2𝜃] and 𝐷(𝜃)=(4−𝛼)²+𝛼²−2𝛼(4−𝛼)cos𝜃.

Let ℑℎ=0, we have sin𝜃(𝛼²−3𝛼+4−𝛼cos𝜃)=0. It is easy to see that sin𝜃=0. Otherwise, |cos𝜃|=|(𝛼²−3𝛼+4)/𝛼|>1 for 0<𝛼<2, which makes a contradiction. Thus, 𝜃=0 or 𝜋. In particular, ℜℎ=−𝛼 for 𝜃=𝜋.

Therefore, the interval of absolute stability of methods (3.4) is (−𝛼,0), where 0<𝛼<2. It means that the maximal length of stability interval can not exceed 2, which is in agreement with Corollary 3.3.

4. LMMs with Longest Stability Interval

In this section, we mainly consider the three-order explicit linear four-step methods with the first and second generating polynomials𝜌(𝜉)=(𝜉−1)(𝜉−𝑎)(𝜉−𝑏)(𝜉−𝑐),𝜎(𝜉)=𝛽₃𝜉³+𝛽₂𝜉²+𝛽₁𝜉+𝛽₀, (4.1) where 𝑎,𝑏,𝑐 are real numbers or one is a real number and the other two are a pair of conjugated complex numbers. For simplicity, we define 𝑚=𝑎+𝑏+𝑐 and 𝑙=𝑎𝑏+𝑏𝑐+𝑎𝑐. Due to (2.6), we have𝛽₃=112(−5𝑚−𝑙−5𝑎𝑏𝑐+23)−𝛽₀,𝛽₂=112(−8𝑚+8𝑙+16𝑎𝑏𝑐−16)+3𝛽₀,𝛽₁=112(𝑚+5𝑙−23𝑎𝑏𝑐+5)−3𝛽₀, (4.2) then the above three-order four-step methods can be written as𝑦_𝑛+4−(𝑚+1)𝑦_𝑛+3+(𝑚+𝑙)𝑦_𝑛+2−(𝑙+𝑎𝑏𝑐)𝑦_𝑛+1+𝑎𝑏𝑐𝑦_𝑛=ℎ(𝛽₃𝑓_𝑛+3+𝛽₂𝑓_𝑛+2+𝛽₁𝑓_𝑛+1+𝛽₀𝑓_𝑛), (4.3) where the additional starting values 𝑦₁, 𝑦₂, and 𝑦₃ should be calculated to an accuracy at least as high as the local accuracy of (4.3). The method used to evaluate 𝑦_𝑘(𝑘=1,2,3) must be a one-step method such as Runge-Kutta method.

Here, the error constant is𝐶₄=172(3𝑚+3𝑙+27𝑎𝑏𝑐+27)+𝛽₀, (4.4) with 𝑎,𝑏,𝑐, and 𝛽₀ free parameters.

To guarantee the convergence, according to [2, Theorem 2.2], we require methods (4.3) to be zero stable; that is, the first generating polynomial 𝜌(𝜉) satisfies the root condition. Thus,|𝑎|≤1,|𝑏|≤1,|𝑐|≤1,𝑎≠𝑏≠𝑐,𝑎≠1,𝑏≠1,𝑐≠1,when|𝑎|=|𝑏|=|𝑐|=1. (4.5)

It is obvious that all methods satisfying conditions 𝐶₄≠0 and (4.5) are convergent three-order linear four-step methods. So, these two conditions are always assumed to be held in the following.

Making the transformation 𝜉=(1+𝑧)/(1−𝑧) to characteristic equation Π(𝜉;ℎ)=0, then we have𝑎₀𝑧⁴+𝑎₁𝑧³+𝑎₂𝑧²+𝑎₃𝑧+𝑎₄=0, (4.6) where 𝑎₀=2(1+𝑎)(1+𝑏)(1+𝑐)−1ℎ[3(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)+8𝛽₀𝑎],₁=−2[(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)]−ℎ[(1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)−8𝛽₀𝑎],₂=−2[(1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)]+4ℎ[3𝑎(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)],₃=2(1−𝑎)(1−𝑏)(1−𝑐)+𝑎ℎ[(1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)],₄=−ℎ(1−𝑎)(1−𝑏)(1−𝑐). (4.7)

Since the transformation 𝜉=(1+𝑧)/(1−𝑧) maps the circle |𝜉|=1 into the imaginary axis ℜ𝑧=0, the interior of the circle into the half-plane ℜ𝑧<0, and the point 𝜉=1 into 𝑧=0 (see Section 3.7 in [8]). Appeal to the well-known Routh-Hurwitz criterion (see Section 1.9 in [2]), the necessary and sufficient conditions that Π(𝜉;ℎ) is Schur polynomial are equivalent to the conditions as follows:𝑎_𝑗>0(𝑗=0,1,2,3,4),𝑎₁𝑎₂𝑎₃−𝑎₀𝑎²₃−𝑎₄𝑎²₁>0. (4.8)

Define 𝜙_𝑖∶={ℎ∈ℛ∣𝑎_𝑖>0} for 𝑖=0,1,2,3,4 and 𝜑(𝑎,𝑏,𝑐,𝛽₀)∶={ℎ∈ℛ∣allrootsofΠ(𝜉;ℎ)=0havemoduluslessthan1forfixed𝑎,𝑏,𝑐,𝛽₀𝜑},₁(𝑎,𝑏,𝑐,𝛽₀)∶={ℎ∈ℛ∣𝑎_𝑗>0,𝑗=0,1,2,3,4forfixed𝑎,𝑏,𝑐,𝛽₀}, (4.9) then it is obvious that 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)=𝜙₀∩𝜙₁∩𝜙₂∩𝜙₃∩𝜙₄.

Lemma 4.1. If (1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)≥0, then 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)=∅ for any convergent three-order four-step methods (4.3).

Proof. For convergent methods (4.3), we have (4.5). So, it is easy to see that ℎ<0 from 𝑎₄>0, that is, 𝜙₄=(−∞,0) and 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)={ℎ∈ℛ∣𝑎₀>0,𝑎₁>0,𝑎₂>0,𝑎₃>0,ℎ<0forfixed𝑎,𝑏,𝑐,𝛽₀}. (4.10)(I)If (1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)>0, then 𝜙₃=(−2(1−𝑎)(1−𝑏)(1−𝑐)(1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐),0). (4.11) However, since |𝑎|<1, |𝑏|<1, |𝑐|<1, we also have (1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)>0. Thus, 𝜙₂=(2[(1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)](4/3)(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐),+∞)⊂(0,+∞). (4.12) It is easy to see that 𝜙₂∩𝜙₃=∅, that is, 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)=∅.(II)If (1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)=0, then (1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)=0. Hence, 𝜙₂=(0,+∞), which means that 𝜙₂∩𝜙₄=∅, that is, 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)=∅.

Lemma 4.2. For convergent three-order four-step methods (4.3), one has 𝜙₂∩𝜙₃∩𝜙₄=(−∞,0), if (1−𝑎)(1−𝑏)(1−𝑐)−3(1−𝑎𝑏𝑐)<0. (4.13)

Proof. From condition (4.13), we have (1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)<0. Similarly, we also have (1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)<0.
After simple calculations, we can obtain 𝜙₀𝜙=(−∞,0),₁=(−∞,−2(1−𝑎)(1−𝑏)(1−𝑐)𝜙(1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)),₂=(−∞,2(1+𝑎)(1+𝑏)(1+𝑐)−8(1+𝑎𝑏𝑐)(4/3)(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)). (4.14) It is easy to see that 𝜙₀∩𝜙₁∩𝜙₂=(−∞,0).

Theorem 4.3. The longest interval of absolute stability for convergent LMMs (4.3) with (4.13) is (6(𝑎+𝑏+𝑐+2−𝑎𝑏𝑐)−10+𝑎+𝑏+𝑐+2𝑎𝑏+2𝑏𝑐+2𝑎𝑐+𝑎𝑏𝑐,0), (4.15) where certain fixed parameters 𝑎,𝑏,𝑐 satisfying 𝛽₀=3(1+𝑎)²(1+𝑏)²(1+𝑐)²+(1−𝑎)²(1−𝑏)²(1−𝑐)²−48(𝑎+𝑏+𝑐−𝑎𝑏𝑐+2)(1+𝑎𝑏𝑐)(1+𝑎)(1+𝑏)(1+𝑐)−𝑎+𝑏+𝑐−𝑎𝑏𝑐+2(1−𝑎𝑏𝑐)(1−𝑎)(1−𝑏)(1−𝑐)−3(1−𝑎𝑏𝑐)².3(𝑎+𝑏+𝑐−𝑎𝑏𝑐+2) (4.16)

Proof. According to Lemmas 4.1 and 4.2, we have 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)={ℎ∈ℛ∣𝑎₀>0,𝑎₁>0,ℎ<0forfixed𝑎,𝑏,𝑐,𝛽₀}. (4.17)
Dividing the real axis into three parts as 𝑅₁1=(−∞,81(1+𝑎)(1+𝑏)(1+𝑐)−2𝑅(1+𝑎𝑏𝑐)],₂1=(81(1+𝑎)(1+𝑏)(1+𝑐)−21(1+𝑎𝑏𝑐),−124(1−𝑎)(1−𝑏)(1−𝑐)+2𝑅(1−𝑎𝑏𝑐)),₃1=[−124(1−𝑎)(1−𝑏)(1−𝑐)+2(1−𝑎𝑏𝑐),+∞). (4.18)
Let 𝑑₁=2(1+𝑎)(1+𝑏)(1+𝑐)(1/3)(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)+8𝛽₀,𝑑₂=−2[(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)](1+𝑎)(1+𝑏)(1+𝑐)−4(1−𝑎𝑏𝑐)−8𝛽₀. (4.19)(I)If 𝛽₀∈𝑅₁, then 𝜙₁=(−∞,2[(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)]−[(1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)−8𝛽₀]𝜙),₀=(2(1+𝑎)(1+𝑏)(1+𝑐)(1/3)(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)+8𝛽₀,+∞), (4.20) that is, 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)=(𝑑₁,0). Hence, 𝜑₁(𝑎,𝑏,𝑐,𝛽_0,1)=(3(1+𝑎)(1+𝑏)(1+𝑐)−10+𝑚+2𝑙+𝑎𝑏𝑐,0)⊂𝜑₁(𝑎,𝑏,𝑐,𝛽₀) (4.21) for any 𝛽₀∈𝑅₁ and 𝛽_0,1=(1/8)(1+𝑎)(1+𝑏)(1+𝑐)−(1/2)(1+𝑎𝑏𝑐).(II)If 𝛽₀∈𝑅₂, then 𝜙₁=(−2[(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)](1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)−8𝛽₀𝜙,+∞),₀=(2(1+𝑎)(1+𝑏)(1+𝑐)(1/3)(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)+8𝛽₀,+∞), (4.22) that is, 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)=(𝑑₁,0)∩(𝑑₂,0) for fixed 𝑎,𝑏,𝑐.
In fact, 𝑑₁ decreases and 𝑑₂ increases monotonously with 𝛽₀ increasing on 𝑅₂. So, 𝜑₁(𝑎,𝑏,𝑐,𝛽_0,2)=(6(𝑎+𝑏+𝑐+2−𝑎𝑏𝑐)−10+𝑚+2𝑙+𝑎𝑏𝑐,0)⊂𝜑₁(𝑎,𝑏,𝑐,𝛽₀), (4.23) for any 𝛽₀∈𝑅₂ and 𝑑₁=𝑑₂, that is, 𝛽_0,2=3(1+𝑎)²(1+𝑏)²(1+𝑐)²+(1−𝑎)²(1−𝑏)²(1−𝑐)²−48(𝑎+𝑏+𝑐−𝑎𝑏𝑐+2)(1+𝑎𝑏𝑐)(1+𝑎)(1+𝑏)(1+𝑐)−𝑎+𝑏+𝑐−𝑎𝑏𝑐+216(1−𝑎𝑏𝑐)(1−𝑎)(1−𝑏)(1−𝑐)−48(1−𝑎𝑏𝑐)².48(𝑎+𝑏+𝑐−𝑎𝑏𝑐+2) (4.24)(III) If 𝛽₀∈𝑅₃, then 𝜙₁=(−2[(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)](1+𝑎)(1+𝑏)(1+𝑐)−4(1+𝑎𝑏𝑐)−8𝛽₀𝜙,+∞),₀=(−∞,2(1+𝑎)(1+𝑏)(1+𝑐)(1/3)(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)+8𝛽₀), (4.25) that is, 𝜑₁(𝑎,𝑏,𝑐,𝛽₀)=(𝑑₂,0). Therefore, 𝜑₁(𝑎,𝑏,𝑐,𝛽_0,3)=(−3[(1−𝑎)(1−𝑏)(1−𝑐)−4(1−𝑎𝑏𝑐)]−10+𝑚+2𝑙+𝑎𝑏𝑐,0)⊂𝜑₁(𝑎,𝑏,𝑐,𝛽₀), (4.26) for any 𝛽₀∈𝑅₃ and 𝛽_0,3=−(1/24)(1−𝑎)(1−𝑏)(1−𝑐)+(1/2)(1−𝑎𝑏𝑐).
Comparing above three cases by Matlab software, we obtain 𝜑₁(𝑎,𝑏,𝑐,𝛽_0,2)⊃𝜑₁(𝑎,𝑏,𝑐,𝛽_0,1) and 𝜑₁(𝑎,𝑏,𝑐,𝛽_0,2)⊃𝜑₁(𝑎,𝑏,𝑐,𝛽_0,3). Hence, 𝜑₁(𝑎,𝑏,𝑐,𝛽_0,2)⊃𝜑₁(𝑎,𝑏,𝑐,𝛽₀) for any 𝛽₀∈𝑅.
If 𝛽₀=𝛽_0,2 and ℎ∈(6(𝑎+𝑏+𝑐+2−𝑎𝑏𝑐)/(−10+𝑚+2𝑙+𝑎𝑏𝑐),0), after direct calculations by Matlab software, we arrive at 𝑎₁𝑎₂𝑎₃−𝑎₀𝑎²₃−𝑎₄𝑎²₁>0, which means that the second condition of (4.8) holds.

Corollary 4.4. The interval of absolute stability of convergent LMMs (4.3) with (4.13) tends to the whole negative real axis, that is, 𝐶₄→1 and 𝛽₀→0, when 𝑎→1, 𝑏→1 and 𝑐→1.

As we known, explicit linear multistep methods cannot be A-stable. In fact, LMMs (4.3) are no longer convergent when 𝑎=𝑏=𝑐=1. Corollary 4.4 implies that the stability interval can increase infinitely, so we can find convergent methods with sufficiently large stability interval.

Corollary 4.5. For the convergent three-order linear four-step methods 𝑦_𝑛+4−(𝑎+1)𝑦_𝑛+3+𝑎𝑦_𝑛+2=ℎ(𝛽₃𝑓_𝑛+3+𝛽₂𝑓_𝑛+2+𝛽₁𝑓_𝑛+1+𝛽₀𝑓_𝑛), (4.27) the stability interval is (−6(𝑎+2)/(10−𝑎),0) if 𝛽₀=(𝑎²−7𝑎−3)/12(𝑎+2); that is, its maximal length cannot exceed 2.

Corollary 4.6. If 𝑎=𝑥+𝑖𝑧, 𝑏=𝑥−𝑖𝑧 with 𝑖²=−1 and 𝑐=0 in (4.1), then one has the convergent three-order linear four-step methods 𝑦_𝑛+4−(2𝑥+1)𝑦_𝑛+3+(𝑥²+2𝑥+𝑧²)𝑦_𝑛+2−(𝑥²+𝑧²)𝑦_𝑛+1=ℎ(𝛽₃𝑓_𝑛+3+𝛽₂𝑓_𝑛+2+𝛽₁𝑓_𝑛+1+𝛽₀𝑓_𝑛). (4.28) Furthermore, if we choose 𝛽₀=6+4𝑥−𝑥²+2𝑥³+𝑥⁴+2𝑥𝑧²+2𝑥²𝑧²−5𝑧²+𝑧⁴,24(𝑥+1) (4.29) then the stability interval of (4.28) is (12(𝑥+1)/(𝑧²+2𝑥²+2𝑥−10),0), that is, its maximum length cannot exceed 4.

At the end of this section, we use BLT to study the stability region of (4.3).

In Figure 1, we plot the region 𝑆, the interval length 𝛼 and the maximal height along the positive imaginary axis Imℎ of absolute stability.

Figure 1 

Region of absolute stability.

It is obvious that 𝛼, 𝐼ℎ and the area of 𝑆 are quite different for different 𝛽₀ of LMMs (4.3) with 𝑎=𝑏=𝑐=0 in Figure 2. Furthermore, when 𝛽₀>0.25, for example, 𝛽₀=0.3, there is a loop on the left of its curve. It demonstrates that the map ℎ(𝜃) is no longer one to one on this interval, which means that the interior of this loop does not belong to the stability region. Similarly, when 𝛽₀<−0.05, for example, 𝛽₀=−0.25, there are two loops on the right part of its curve. Therefore, the curve constructs one single connected region of corresponding numerical method only when −0.05<𝛽₀<0.25.

Figure 2 

Regions of absolute stability for LMMs (4.3) with 𝑎=𝑏=𝑐=0.

In Table 1, the upper bounds 𝐵_𝑈 of 𝛽₀ for curve with no loop on the left and the lower bounds 𝐵_𝐿 of 𝛽₀ for curve with no loop on the right are given. Moreover, 𝛽₀(𝛼) and 𝛽₀(𝐼ℎ) are evaluated on [𝐵_𝐿,𝐵_𝑈] such that 𝛼 and 𝐼ℎ are maximal, respectively. The corresponding maximums 𝛼 and 𝐼ℎ are also given out. Especially, several stability regions of LMMs (4.3) with different 𝑎,𝑏,𝑐 are shown in Figure 3.

Table 1 

Stability results of LMMs (4.3).

(a)

(b)

(c)

(d)

Figure 3 

Regions of absolute stability for LMMs (4.3).

Finally, some stability intervals are calculated for special three-order four-step methods.

It is shown that the maximal length of stability interval for methods (4.27) cannot exceed 2 and the maximal length for methods (4.28) cannot exceed 4 in Table 2. These results are in agreement with the conclusions of Corollaries 4.5 and 4.6.

Table 2 

Stability intervals of LMMs (4.27) and LMMs (4.28).

5. Numerical Experiments

In this section, some examples are given to demonstrate the validity of our proposed methods.

Example 5.1. Consider the ordinary differential system (see page 229, Section 8.5 in [8]) 𝑦^′(𝑡)=𝐴𝑦(𝑡), (5.1) where 𝑦(0)=(1,0,−1)^𝑇, 𝐴=(^{−2119−2019−2120}_40−40−40) and its exact solutions are 𝑦₁1(𝑡)=21exp(−2𝑡)+2𝑦exp(−40𝑡)(cos40𝑡+sin40𝑡),₂1(𝑡)=21exp(−2𝑡)−2𝑦exp(−40𝑡)(cos40𝑡+sin40𝑡),₃(𝑡)=−exp(−40𝑡)(cos40𝑡−sin40𝑡). (5.2)

As we all know, this equation is stiff because its solutions contain a component which decays much more rapidly than the other.

Then, we apply the following three numerical methods to problem (5.1).(I)Explicit three-step Adams-Bashforth methods of order three (AB3) 𝑦_𝑛+1=𝑦_𝑛+ℎ12(23𝑓_𝑛−16𝑓_𝑛−1+5𝑓_𝑛−2). (5.3)(II)Linear four-step methods (4.3) of order three with 𝑎=𝑏=𝑐=0 and 𝛽₀=0.25 (LMM1) 𝑦_𝑛+4−𝑦_𝑛+35=ℎ(3𝑓_𝑛+3−7𝑓12^𝑛+2−13𝑓_𝑛+1+14𝑓_𝑛). (5.4)(III)Linear four-step methods (4.3) of order three with 𝑎=𝑏=𝑐=0.9, 𝛽₀=0.01 (LMM2) 𝑦_𝑛+4−3.7𝑦_𝑛+3+5.13𝑦_𝑛+2−3.159𝑦_𝑛+1+0.729𝑦_𝑛=ℎ(0.2754𝑓_𝑛+3−0.5113𝑓_𝑛+2+0.2269𝑓_𝑛+1+0.01𝑓_𝑛). (5.5)

Here, we let the global error 𝐸_𝑛∶=𝑌(𝑡_𝑛)−𝑌_𝑛, which 𝑌(𝑡_𝑛) and 𝑌_𝑛 denote the exact solutions and numerical solutions, respectively. Then, we attempt to solve (5.1) by these three methods with |𝐸_𝑛|≤10⁻². It can be found that the steplengths ℎ taken to reach the point 𝑡_𝑛=0.1 are as follows.

In Table 3, the steplength of the linear four-step methods (LMM1 and LMM2) can be chosen much larger than that of the classical Adams-Bashforth methods (AB3) with the same order. This implies that the proposed methods (4.3) have much longer absolute stability interval. On the other hand, the selection of steplength as indicted in Table 3 for LMM1 and LMM2 gives another suggested fact that the steplength can increase gradually with 𝑎→1, 𝑏→1, and 𝑐→1, which is predicted by Corollary 4.4.

Table 3 

Numerical solutions and errors of AB3, LMM1 and LMM2 for (5.1).

Example 5.2. Consider the problem (see page 213, Section 6.1 in [2]) 𝑦^′(𝑡)=𝐵𝑦(𝑡)+𝑓(𝑡), (5.6) where 𝑦(0)=(0,0)^𝑇, 𝐵=(⁻²¹_998−999), 𝑓(𝑡)=(^2sin𝑡_{999(cos𝑡−sin𝑡)}) and its exact solution is 𝑦₁𝑦(𝑡)₂(𝑡)=𝐾₁11exp(−𝑡)+𝐾₂1exp(−1000𝑡)+−998sin𝑡cos𝑡. (5.7)

Here, 𝐾₁=−1/999 and 𝐾₂=1/999. Then, system (5.6) is also stiff because it has a stiffness ratio of 1000.

To compare the proposed methods (4.3) with other well-known stiff methods, we consider the Adams-Moulton methods of order three (AM3) (1−𝛼)(𝑦_𝑛+2−𝑦_𝑛−ℎ3(𝑓_𝑛+2+4𝑓_𝑛+1+𝑓_𝑛))+𝛼(𝑦_𝑛+2−𝑦_𝑛+1−ℎ12(5𝑓_𝑛+2+8𝑓_𝑛+1−𝑓_𝑛))=0, (5.8) where 𝛼≠1 and 0≤𝛼<2. Such numerical methods have good stability because their absolute stability interval is (6𝛼/(𝛼−2),0) and they are appropriate for moderately stiff systems (see [8]).

Applying LMM1 and AM3 to system (5.6), we plot the exact solutions and numerical solutions with 𝛼=1.99, ℎ=0.0002 in Figure 4, respectively. It can be noted that the accuracy of AM3 is indeed better than that of LMM1 in the beginning interval. However, LMM1 can repair the initial error and obtain the same accuracy as AM3 eventually, which can be seen from the second component (the bottom figure). This phenomenon reveals that the proposed methods (LMM1) have good stability.

(a)

(b)

Figure 4 

Comparison of exact solutions with numerical solutions for (5.6).

Example 5.3. To illustrate the efficiency of the proposed methods (4.3) further, we will do some comparisons with other famous stiff methods, such as, BDF methods of order three (BDF3) 𝑦_𝑛+3−18𝑦11^𝑛+2+9𝑦11^𝑛+1−2𝑦11^𝑛=611ℎ𝑓_𝑛+3. (5.9)

It is well known that the BDF methods are central to the construction of efficient algorithms for handling stiff systems. In fact, they play the same role in stiff problems as the Adams methods do in nonstiff ones.

In Table 4, the numerical solutions for (5.1) are all calculated with the same steplength ℎ=0.0025. It is noted that the accuracy of LMM1 is no worse than that of the other two stiff methods (BDF3 and AM3) and LMM1 almost behave as well as BDF3.

Table 4 

Numerical solutions and errors of LMM1, BDF3, and AM3 for (5.1).

Example 5.4. It is well known that, for 𝑠-stage explicit Runge-Kutta methods of order 𝑠(𝑠=1,2,3,4), the intervals of absolute stability become larger as the order increases. Meanwhile, the opposite happened for the explicit linear multistep methods. Therefore, people are willing to use the explicit Runge-Kutta methods in the physical problems. However, the maximal length of stability interval for methods (4.28) approach 4 in Table 2, while that stability interval length for three-stage explicit Runge-Kutta methods with the same order is about 2.51. Furthermore, Corollary 4.4 implies that the stability interval can be infinite for LMMs (4.3), while the explicit Runge-Kutta methods only have infinite intervals of absolute stability. In order to illustrate the efficiency of the proposed methods further, we will give a comparison between the explicit three-order linear four-step methods (4.28) and the three-stage explicit Runge-Kutta methods of order three.

Let 𝑥=0.99 and 𝑧=0.1 in methods (4.28), then 𝛼=3.9472. In Figure 5, we plot the error curves of Runge-Kutta method and the proposed method (4.28) with 1000ℎ=2<2.51 (also less than 𝛼) for problem (5.6). It can be noted that both numerical methods behave well in stability, which is in agreement with the theory in [1]. Furthermore, we plot the error curves of both methods with 2.51<1000ℎ=3<𝛼 in Figure 6. It can be found that the proposed method is stable, while the Runge-Kutta method is unstable. This implies that the stability interval of the proposed methods (4.28) is greatly improved.

(a)

(b)

Figure 5 

Errors of Runge-Kutta method and LMM (4.28) with ℎ=0.002.

(a)

(b)

Figure 6 

Errors of Runge-Kutta method and LMM (4.28) with ℎ=0.003.

Remark 5.5. In the above three examples, there are two kinds of comparison between the proposed methods (4.3) and the other existing methods. In Example 5.1, it is shown that our methods in this paper have longer absolute stability interval than the classical Adams-Bashforth methods. In Examples 5.2 and 5.3, other comparisons between the proposed methods (4.3) and some well-known stiff methods are given. Here, we do not claim that the proposed methods in this paper are better than AM3 and BDF3, after all, the proposed methods (4.3) are explicit methods. However, it can give a comparative results at least. These comparisons show in depth that the proposed methods (4.3) are favorable in applications to the stiff systems and improved in stability for the classic Adams-Bashforth methods.

Remark 5.6. Although implicit methods are so favoured that explicit linear multistep methods are seldom used on their own, they do, however, play an important ancillary role in predictor-corrector pairs (see [8, Section 3.8]). The common software for ODEs is not based on Adams-Bashforth methods alone, but on predictor-corrector methods with Adams-Bashforth predictor and Adams-Moulton corrector. However, the research in our paper can supply a kind of tools for practical computation when Adams-Bashforth methods or predictor-corrector methods are applied.

6. Conclusions

In this paper, several three-order explicit linear four-step methods are proposed, which possess far longer intervals of absolute stability than the classical Adams-Bashforth methods with the same order. Because the steplength of the proposed methods can be chosen much larger, these kinds of explicit linear multistep methods are found to be more adaptive.

Acknowledgments

The authors wish to thank the anonymous referees for their valuable comments which helped us to improve the present paper. This paper was supported by projects from Science Research Foundation (HITC200710) and Natural Scientific Research Innovation Foundation (HIT.NSRIF.2009053) in Harbin Institute of Technology.