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

Authors: Y. Xu; J. J. Zhao
Journal: Discrete Dynamics in Nature and Society
 (2010)

---

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


---

## Body

## 1. Introduction

For the initial value problem of the ordinary differential equation (ODE)(1.1)y'(t)=f(t,y(t)),
where f:[t0,tend]×ℛm→ℛm and y(t0)=y0 with y0∈ℛm, 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 y0, a linear multistep (k-step) method requires (k) starting values y0,y1,…,yk-1 and a multistep (k-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 (k-step) methods to the initial value problem (1.1), we obtain the recurrence relation(2.1)∑j=0kαjyn+j=h∑j=0kβjfn+j,
where yn denotes an approximation to the solution y(tn), fn+j:=f(tn+j,yn+j), tn+j=tn+jh for j=0,1,…,k, the constant steplength h>0, and k starting conditions y0,y1,…,yk-1 are required. Here, αj and βj(j=0,1,…,k) are constants subject to the condition αk=1,|α0|+|β0|≠0. If βk=0, then the corresponding methods (2.1) are explicit, and implicit otherwise.Then, we define the first and second generating polynomials by(2.2)ρ(ζ):=∑j=0kαjζj,σ(ζ):=∑j=0kβjζj,
where ζ∈𝒞 is a dummy variable.Consider the scalar test equation(2.3)y′(t)=λy(t),
where λ∈𝒞 and ℜλ<0. Its characteristic polynomial can be written as(2.4)Π(ξ;h¯):=ρ(ξ)-h¯σ(ξ),
where h¯=λh.Here, we quote some important definitions (see Sections3.2, 3.8, and 6.3 in the reference [2]).Definition 2.1.
The set(2.5)Ω={h¯∈𝒞∣allrootsofΠ(ξ;h¯)=0havemoduluslessthan1}
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 h¯-plane is called the interval of absolute stability.Definition 2.3.
The multistep method is said to be A-stable ifΩ⊃𝒞-, where 𝒞-:={z∈𝒞∣ℜz<0}.Definition 2.4.
The linear multistep methods (2.1) are said to be of order p if C0=C1=⋯=Cp=0 and Cp+1≠0, where
(2.6)C0=∑j=0kαj≡ρ(1),C1=∑j=0k(jαj-βj)≡ρ'(1)-σ(1),Cq=∑j=0k[1q!jqαj-1(q-1)!jq-1βj],q=2,3,….


## 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 h¯-plane be defined by the requirement that for all h¯∈∂Ω, one of the roots of Π(ξ;h¯) has modulus 1, that is, is of the form ξ=exp(iθ). Thus, we have Π(exp(iθ);h¯)=ρ(exp(iθ))-h¯σ(exp(iθ))=0.This identity is readily solved forh¯, and we introduce the root locus curve(3.1)h¯:θ→h¯(θ)=ρ(exp(iθ))σ(exp(iθ)),
which maps [0,2π] onto Γ, where Γ:={h¯∈𝒞∣ there exsist ξ,|ξ|=1,Π(ξ;h¯)=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 h¯-plane are h¯1=0, h¯2=ρ(-1)/σ(-1). Furthermore, if the real interval
(3.2)Δ=(ρ(-1)σ(-1),0)
does not intersect with the curve Γ and there exists ɛ0∈Δ such that all roots of Π(ξ;ɛ0)=0 have modulus less than 1 in complex h¯-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(3.3)ρ(ξ)=∑j=02αjξj,σ(ξ)=∑j=01βjξj.Without loss of generality, we assume that-1<α0<1.Corollary 3.2.
The estimated interval of absolute stability of two-order explicit linear two-step methods (3.3) is Δ=(-1-α0,0), where -1<α0<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ΔE;p,s or ΔI;p,s be the estimated interval of absolute stability of p-order explicit or implicit linear s-step methods. Then(I)
ΔE;3,3=(6(α0+α2)/(α0-α2+10),0), where -1<α0<1, α0+α2<0 and α0-α2<2;(II)
ΔI;3,2=(-6(α0+1)/(1-α0),0), where -1<α0<1;(III)
ΔI;4,3=(-3(α0+α2)/(1+α0),0), where -1<α0<1, α0+α2<0 and α0-α2<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(3.4)yn+2-αyn+1-(1-α)yn=h2[(4-α)fn+1-αfn],
where 0<α<2.Applying to the test equation (2.3), we have the stability polynomial (3.5)Π(ξ;h¯)=ξ2-αξ-(1-α)-h¯2[(4-α)ξ-α].Defineξ=exp(iθ), then the roots of above characteristic polynomial can be written as h¯=ℜh¯+iℑh¯, where ℜh¯=(2α/D(θ))(4cosθ-cos2θ-3), ℑh¯=(2/D(θ))[2(α2-3α+4)sinθ-αsin2θ] and D(θ)=(4-α)2+α2-2α(4-α)cosθ.Letℑh¯=0, we have sinθ(α2-3α+4-αcosθ)=0. It is easy to see that sinθ=0. Otherwise, |cosθ|=|(α2-3α+4)/α|>1 for 0<α<2, which makes a contradiction. Thus, θ=0 or π. In particular, ℜh¯=-α 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(4.1)ρ(ξ)=(ξ-1)(ξ-a)(ξ-b)(ξ-c),σ(ξ)=β3ξ3+β2ξ2+β1ξ+β0,
where a,b,c are real numbers or one is a real number and the other two are a pair of conjugated complex numbers. For simplicity, we define m=a+b+c and l=ab+bc+ac. Due to (2.6), we have(4.2)β3=112(-5m-l-5abc+23)-β0,β2=112(-8m+8l+16abc-16)+3β0,β1=112(m+5l-23abc+5)-3β0,
then the above three-order four-step methods can be written as(4.3)yn+4-(m+1)yn+3+(m+l)yn+2-(l+abc)yn+1+abcyn=h(β3fn+3+β2fn+2+β1fn+1+β0fn),
where the additional starting values y1, y2, and y3 should be calculated to an accuracy at least as high as the local accuracy of (4.3). The method used to evaluate yk(k=1,2,3) must be a one-step method such as Runge-Kutta method.Here, the error constant is(4.4)C4=172(3m+3l+27abc+27)+β0,
with a,b,c, and β0 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,(4.5)|a|≤1,|b|≤1,|c|≤1,a≠b≠c,a≠1,b≠1,c≠1,when|a|=|b|=|c|=1.It is obvious that all methods satisfying conditionsC4≠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+z)/(1-z) to characteristic equation Π(ξ;h¯)=0, then we have(4.6)a0z4+a1z3+a2z2+a3z+a4=0,
where (4.7)a0=2(1+a)(1+b)(1+c)-h¯[13(1-a)(1-b)(1-c)-4(1-abc)+8β0],a1=-2[(1-a)(1-b)(1-c)-4(1-abc)]-h¯[(1+a)(1+b)(1+c)-4(1+abc)-8β0],a2=-2[(1+a)(1+b)(1+c)-4(1+abc)]+h¯[43(1-a)(1-b)(1-c)-4(1-abc)],a3=2(1-a)(1-b)(1-c)+h¯[(1+a)(1+b)(1+c)-4(1+abc)],a4=-h¯(1-a)(1-b)(1-c).Since the transformationξ=(1+z)/(1-z) maps the circle |ξ|=1 into the imaginary axis ℜz=0, the interior of the circle into the half-plane ℜz<0, and the point ξ=1 into z=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 Π(ξ;h¯) is Schur polynomial are equivalent to the conditions as follows:(4.8)aj>0(j=0,1,2,3,4),a1a2a3-a0a32-a4a12>0.Defineϕi:={h¯∈ℛ∣ai>0} for i=0,1,2,3,4 and (4.9)φ(a,b,c,β0):={h¯∈ℛ∣allrootsofΠ(ξ;h¯)=0havemoduluslessthan1forfixeda,b,c,β0},φ1(a,b,c,β0):={h¯∈ℛ∣aj>0,j=0,1,2,3,4forfixeda,b,c,β0},
then it is obvious that φ1(a,b,c,β0)=ϕ0∩ϕ1∩ϕ2∩ϕ3∩ϕ4.Lemma 4.1.
If(1+a)(1+b)(1+c)-4(1+abc)≥0, then φ1(a,b,c,β0)=∅ 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 h¯<0 from a4>0, that is, ϕ4=(-∞,0) and
(4.10)φ1(a,b,c,β0)={h¯∈ℛ∣a0>0,a1>0,a2>0,a3>0,h¯<0forfixeda,b,c,β0}.(I)
If(1+a)(1+b)(1+c)-4(1+abc)>0, then
(4.11)ϕ3=(-2(1-a)(1-b)(1-c)(1+a)(1+b)(1+c)-4(1+abc),0).
However, since |a|<1, |b|<1, |c|<1, we also have (1-a)(1-b)(1-c)-4(1-abc)>0. Thus,
(4.12)ϕ2=(2[(1+a)(1+b)(1+c)-4(1+abc)](4/3)(1-a)(1-b)(1-c)-4(1-abc),+∞)⊂(0,+∞).
It is easy to see that ϕ2∩ϕ3=∅, that is, φ1(a,b,c,β0)=∅.(II)
If(1+a)(1+b)(1+c)-4(1+abc)=0, then (1-a)(1-b)(1-c)-4(1-abc)=0. Hence, ϕ2=(0,+∞), which means that ϕ2∩ϕ4=∅, that is, φ1(a,b,c,β0)=∅.Lemma 4.2.
For convergent three-order four-step methods (4.3), one has ϕ2∩ϕ3∩ϕ4=(-∞,0), if
(4.13)(1-a)(1-b)(1-c)-3(1-abc)<0.Proof.
From condition (4.13), we have (1-a)(1-b)(1-c)-4(1-abc)<0. Similarly, we also have (1+a)(1+b)(1+c)-4(1+abc)<0.
After simple calculations, we can obtain(4.14)ϕ0=(-∞,0),ϕ1=(-∞,-2(1-a)(1-b)(1-c)(1+a)(1+b)(1+c)-4(1+abc)),ϕ2=(-∞,2(1+a)(1+b)(1+c)-8(1+abc)(4/3)(1-a)(1-b)(1-c)-4(1-abc)).
It is easy to see that ϕ0∩ϕ1∩ϕ2=(-∞,0).Theorem 4.3.
The longest interval of absolute stability for convergent LMMs (4.3) with (4.13) is
(4.15)(6(a+b+c+2-abc)-10+a+b+c+2ab+2bc+2ac+abc,0),
where certain fixed parameters a,b,c satisfying
(4.16)β0=3(1+a)2(1+b)2(1+c)2+(1-a)2(1-b)2(1-c)248(a+b+c-abc+2)-(1+abc)(1+a)(1+b)(1+c)a+b+c-abc+2-(1-abc)(1-a)(1-b)(1-c)-3(1-abc)23(a+b+c-abc+2).Proof.
According to Lemmas4.1 and 4.2, we have
(4.17)φ1(a,b,c,β0)={h¯∈ℛ∣a0>0,a1>0,h¯<0forfixeda,b,c,β0}.
Dividing the real axis into three parts as(4.18)R1=(-∞,18(1+a)(1+b)(1+c)-12(1+abc)],R2=(18(1+a)(1+b)(1+c)-12(1+abc),-124(1-a)(1-b)(1-c)+12(1-abc)),R3=[-124(1-a)(1-b)(1-c)+12(1-abc),+∞).
Let(4.19)d1=2(1+a)(1+b)(1+c)(1/3)(1-a)(1-b)(1-c)-4(1-abc)+8β0,d2=-2[(1-a)(1-b)(1-c)-4(1-abc)](1+a)(1+b)(1+c)-4(1-abc)-8β0.(I)
Ifβ0∈R1, then
(4.20)ϕ1=(-∞,2[(1-a)(1-b)(1-c)-4(1-abc)]-[(1+a)(1+b)(1+c)-4(1+abc)-8β0]),ϕ0=(2(1+a)(1+b)(1+c)(1/3)(1-a)(1-b)(1-c)-4(1-abc)+8β0,+∞),
that is, φ1(a,b,c,β0)=(d1,0). Hence,
(4.21)φ1(a,b,c,β0,1)=(3(1+a)(1+b)(1+c)-10+m+2l+abc,0)⊂φ1(a,b,c,β0)
for any β0∈R1 and β0,1=(1/8)(1+a)(1+b)(1+c)-(1/2)(1+abc).(II)
Ifβ0∈R2, then
(4.22)ϕ1=(-2[(1-a)(1-b)(1-c)-4(1-abc)](1+a)(1+b)(1+c)-4(1+abc)-8β0,+∞),ϕ0=(2(1+a)(1+b)(1+c)(1/3)(1-a)(1-b)(1-c)-4(1-abc)+8β0,+∞),
that is, φ1(a,b,c,β0)=(d1,0)∩(d2,0) for fixed a,b,c.
In fact,d1 decreases and d2 increases monotonously with β0 increasing on R2. So,
(4.23)φ1(a,b,c,β0,2)=(6(a+b+c+2-abc)-10+m+2l+abc,0)⊂φ1(a,b,c,β0),
for any β0∈R2 and d1=d2, that is,
(4.24)β0,2=3(1+a)2(1+b)2(1+c)2+(1-a)2(1-b)2(1-c)248(a+b+c-abc+2)-(1+abc)(1+a)(1+b)(1+c)a+b+c-abc+2-16(1-abc)(1-a)(1-b)(1-c)-48(1-abc)248(a+b+c-abc+2).(III)
Ifβ0∈R3, then
(4.25)ϕ1=(-2[(1-a)(1-b)(1-c)-4(1-abc)](1+a)(1+b)(1+c)-4(1+abc)-8β0,+∞),ϕ0=(-∞,2(1+a)(1+b)(1+c)(1/3)(1-a)(1-b)(1-c)-4(1-abc)+8β0),
that is, φ1(a,b,c,β0)=(d2,0). Therefore,
(4.26)φ1(a,b,c,β0,3)=(-3[(1-a)(1-b)(1-c)-4(1-abc)]-10+m+2l+abc,0)⊂φ1(a,b,c,β0),
for any β0∈R3 and β0,3=-(1/24)(1-a)(1-b)(1-c)+(1/2)(1-abc).
Comparing above three cases by Matlab software, we obtainφ1(a,b,c,β0,2)⊃φ1(a,b,c,β0,1) and φ1(a,b,c,β0,2)⊃φ1(a,b,c,β0,3). Hence, φ1(a,b,c,β0,2)⊃φ1(a,b,c,β0) for any β0∈R.
Ifβ0=β0,2 and h¯∈(6(a+b+c+2-abc)/(-10+m+2l+abc),0), after direct calculations by Matlab software, we arrive at a1a2a3-a0a32-a4a12>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, C4→1 and β0→0, when a→1, b→1 and c→1.As we known, explicit linear multistep methods cannot be A-stable. In fact, LMMs (4.3) are no longer convergent when a=b=c=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.27)yn+4-(a+1)yn+3+ayn+2=h(β3fn+3+β2fn+2+β1fn+1+β0fn),
the stability interval is (-6(a+2)/(10-a),0) if β0=(a2-7a-3)/12(a+2); that is, its maximal length cannot exceed 2.Corollary 4.6.
Ifa=x+iz, b=x-iz with i2=-1 and c=0 in (4.1), then one has the convergent three-order linear four-step methods
(4.28)yn+4-(2x+1)yn+3+(x2+2x+z2)yn+2-(x2+z2)yn+1=h(β3fn+3+β2fn+2+β1fn+1+β0fn).
Furthermore, if we choose
(4.29)β0=6+4x-x2+2x3+x4+2xz2+2x2z2-5z2+z424(x+1),
then the stability interval of (4.28) is (12(x+1)/(z2+2x2+2x-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 Figure1, we plot the region S, the interval length α and the maximal height along the positive imaginary axis Imh¯ of absolute stability.Figure 1
Region of absolute stability.It is obvious thatα, Ih and the area of S are quite different for different β0 of LMMs (4.3) with a=b=c=0 in Figure 2. Furthermore, when β0>0.25, for example, β0=0.3, there is a loop on the left of its curve. It demonstrates that the map h¯(θ) 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<-0.05, for example, β0=-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<0.25.Figure 2
Regions of absolute stability for LMMs (4.3) with a=b=c=0.In Table1, the upper bounds BU of β0 for curve with no loop on the left and the lower bounds BL of β0 for curve with no loop on the right are given. Moreover, β0(α) and β0(Ih) are evaluated on [BL,BU] such that α and Ih are maximal, respectively. The corresponding maximums α and Ih are also given out. Especially, several stability regions of LMMs (4.3) with different a,b,c are shown in Figure 3.Table 1
Stability results of LMMs (4.3).
abcBUBLβ0(α)αβ0(Ih)Ih0.00000.00000.00000.2500−0.05000.25001.20000.00000.72310.25000.25000.25000.2109−0.35000.21091.8519−0.20000.70900.30000.30000.30000.2003−0.40000.20032.0202−0.20000.70470.75000.75000.75000.0578−1.35000.05785.8103−0.20000.96930.90000.90000.90000.0123−1.95000.012313.9252−0.20001.74600.25000.50000.75000.1670−0.85000.16702.9067−0.20000.7620−0.25000.50000.50000.2854−0.40000.28541.8121−0.10000.7795−0.1000−0.1000−0.10000.26220.05000.26220.99660.10000.7092Regions of absolute stability for LMMs (4.3).
(a)(b)(c)(d)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).
a(4.27)b(4.27)c(4.27)β(4.27)α(4.27)a(4.28)b(4.28)c(4.28)β(4.28)α(4.28)0.99990.00000.00000.25001.99990.99990.99990.00000.25003.99960.90000.00000.00000.24741.91200.90000.90000.00000.23913.46470.75000.00000.00000.24431.78370.75000.75000.00000.22852.84730.50000.00000.00000.24171.57920.50000.50000.00000.22402.11800.25000.00000.00000.24311.38490.25000.25000.00000.23241.5999−0.50000.00000.00000.26491.0245−0.5000−0.50000.00000.29690.5716−0.99990.00000.00000.41660.3529−0.9999−0.99990.00000.33330.500E-4


## 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, Section8.5 in [8])
(5.1)y′(t)=Ay(t),
where y(0)=(1,0,-1)T, A=(-2119-2019-212040-40-40) and its exact solutions are
(5.2)y1(t)=12exp(-2t)+12exp(-40t)(cos40t+sin40t),y2(t)=12exp(-2t)-12exp(-40t)(cos40t+sin40t),y3(t)=-exp(-40t)(cos40t-sin40t).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)(5.3)yn+1=yn+h12(23fn-16fn-1+5fn-2).(II)
Linear four-step methods (4.3) of order three with a=b=c=0 and β0=0.25 (LMM1)
(5.4)yn+4-yn+3=h(53fn+3-712fn+2-13fn+1+14fn).(III)
Linear four-step methods (4.3) of order three with a=b=c=0.9, β0=0.01 (LMM2)
(5.5)yn+4-3.7yn+3+5.13yn+2-3.159yn+1+0.729yn=h(0.2754fn+3-0.5113fn+2+0.2269fn+1+0.01fn).Here, we let the global errorEn:=Y(tn)-Yn, which Y(tn) and Yn denote the exact solutions and numerical solutions, respectively. Then, we attempt to solve (5.1) by these three methods with |En|≤10-2. It can be found that the steplengths h taken to reach the point tn=0.1 are as follows.In Table3, 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 a→1, b→1, and c→1, which is predicted by Corollary 4.4.Table 3
Numerical solutions and errors of AB3, LMM1 and LMM2 for (5.1).
tnY(tn)Yn(AB3)En(AB3)Yn(LMM1)En(LMM1)Yn(LMM2)En(LMM2)h=0.0001h=0.004h=0.010.79810.7803-0.01780.7981-0.105E-40.7957−0.241E-20.02000.16270.18010.01740.16270.174E-40.16390.122E-20.00930.20860.19930.01040.109E-20.00790.145E-20.55950.5217-0.03780.5595-0.250E-50.5593−0.146E-30.04000.36360.40100.03740.3636-0.420E-50.36430.709E-30.20770.25880.05110.20760.100E-30.20850.816E-30.44070.4205-0.02020.44070.100E-50.4407−0.551E-40.06000.44630.46610.01980.44630.242E-50.44630.550E-40.12820.1199-0.00830.1280-0.220E-30.12820.854E-40.40450.3994-0.00510.4045-0.122E-60.40450.214E-50.08000.44760.45240.00480.44760.274E-60.44760.324E-50.03830.0228-0.01550.0383-0.517E-40.0383−0.859E-50.39640.39720.800E-30.39640.815E-40.39640.662E-50.10000.42230.4212-0.100E-20.4223-0.512E-40.4223−0.160E-5−0.0019−0.00900.800E-2−0.00190.541E-4−0.00190.221E-5Example 5.2.
Consider the problem (see page 213, Section6.1 in [2])
(5.6)y′(t)=By(t)+f(t),
where y(0)=(0,0)T, B=(-21998-999), f(t)=(2sint999(cost-sint)) and its exact solution is
(5.7)(y1(t)y2(t))=K1exp(-t)(11)+K2exp(-1000t)(1-998)+(sintcost).Here,K1=-1/999 and K2=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) (5.8)(1-α)(yn+2-yn-h3(fn+2+4fn+1+fn))+α(yn+2-yn+1-h12(5fn+2+8fn+1-fn))=0,
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, h=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.Comparison of exact solutions with numerical solutions for (5.6).
(a)(b)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)
(5.9)yn+3-1811yn+2+911yn+1-211yn=611hfn+3.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 Table4, the numerical solutions for (5.1) are all calculated with the same steplength h=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).
tnY(tn)Yn(LMM1)En(LMM1)Yn(BDF3)En(BDF3)Yn(AM3)En(AM3)h=0.0025h=0.0025h=0.00257.9808E-17.9826E-11.7435E-47.9789E-1-1.9274E-47.9778E-1-3.0096E-40.02001.6271E-11.6253E-1-1.7435E-41.6290E-11.9274E-41.6301E-13.0097E-49.2784E-39.2644E-3-1.3993E-59.3240E-34.5652E-59.4174E-31.3908E-45.5952E-15.5966E-11.4495E-45.5935E-1-1.6673E-45.5872E-1-8.0064E-40.04003.6360E-13.6346E-1-1.4496E-43.6377E-11.6673E-43.6440E-18.0064E-42.0771E-12.0798E-12.7493E-42.0743E-1-2.7512E-42.0737E-1-3.3240E-44.4065E-14.4065E-11.7550E-64.4064E-1-1.2265E-54.3975E-1-9.0526E-40.06004.4627E-14.4627E-1-1.7579E-64.4628E-11.2268E-54.4717E-19.0527E-41.2817E-11.2846E-12.8861E-41.2786E-1-3.0907E-41.2702E-1-1.1478E-34.0454E-14.0447E-1-6.3347E-54.0460E-16.3188E-54.0382E-1-7.1785E-40.08004.4761E-14.4767E-16.3343E-54.4754E-1-6.3184E-54.4833E-17.1787E-43.8313E-23.8439E-21.2611E-43.8170E-2-1.4361E-43.6721E-2-1.5922E-33.9645E-13.9640E-1-5.1188E-53.9650E-15.4584E-53.9595E-1-5.0282E-40.10004.2228E-14.2233E-15.1183E-54.2223E-1-5.4580E-54.2278E-15.0285E-41.8894E-3-1.8911E-3-1.7216E-6-1.8947E-3-5.3187E-6-3.4785E-3-1.5891E-3Example 5.4.
It is well known that, fors-stage explicit Runge-Kutta methods of order s(s=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.Letx=0.99 and z=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 1000h=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<1000h=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.Errors of Runge-Kutta method and LMM (4.28) with h=0.002.
(a)(b)Errors of Runge-Kutta method and LMM (4.28) with h=0.003.
(a)(b)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.



---

*Source: 912691-2010-11-07.xml*