# New Phase Fitted and Amplification Fitted Numerov-Type Methods for Periodic IVPs with Two Frequencies

Authors: Yonglei Fang; Xiong You; Zhaoxia Chen
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

Phase fitted and amplification fitted Numerov-type methods for periodic initial value problems with two frequencies are investigated. A one-frequency method and a two-frequency method are constructed. The two new methods both have algebraic order five and are dispersive of order six and dissipative of order five. The two-dimensional absolute stability region for the one-frequency method and the three-dimensional absolute stability region for the two-frequency method are plotted. Numerical experiments are reported to show the efficiency and competence of the two new methods.


---

## Body

## 1. Introduction

In this paper we are interested in the numerical integration of the initial value problem (IVP) of ordinary differential equations in the form(1.1)y′′(x)=f(x,y),y(x0)=y0,y′(x0)=y0′,x∈[x0,xend],
whose solutions exhibit a pronounced oscillatory character. Such problems are frequently encountered in a variety of scientific fields and engineering applications. Since the analytical solutions of these equations are usually not available, the numerical solutions become very important. To our knowledge, most existing methods of fitted type (trigonometrically/exponentially fitted, phase fitted and amplification fitted, etc.) are adapted to one frequency (see [1–3]). We call them one-frequency methods. For example, by phase fitting, Van de Vyver [4] constructed an adapted explicit Numerov-type method for oscillatory problems with one frequency. Other related work can be found in [5–20]. However, for oscillatory problems with many frequencies in the solution, when we apply a one-frequency method, it is difficult to choose the most suitable fitting frequency. Contrast to the one-frequency methods, papers on integrators that are adapted to more than one frequency are few. In [21], Wang obtained a trigonometrically fitted Numerov-type method for periodic IVPs with two frequencies and provided some numerical examples to show that his method associated with two frequencies is more accurate and more efficient than the trigonometrically fitted Numerov methods with one frequency. More recently, Fang et al. [22] gave a trigonometrically fitted explicit Numerov-type method for the numerical integration of oscillatory problems with two frequencies. Inspired by Van de Vyver’s idea in [4], we investigate, in this paper, explicit Numerov-type methods for the numerical integration of periodic IVPs with two frequencies. The paper is organized as follows. In Section 2, the preliminaries of stability and phase properties of explicit Numerov-type methods are recalled. In Section 3, we construct two new phase fitted and amplification fitted explicit Numerov-type methods for periodic initial value problems with two frequencies. In Section 4, the stability and phase properties of the newly derived methods are analyzed. In Section 5, we present some numerical results. Section 6 is devoted to the conclusions.


## 2. Numerov-Type Methods

In this section, we recall some elementary notions about Numerov-type methods. The explicit Numerov-type methods considered here have the form(2.1)Yi=(1+ci)yn-ciyn-1+h2∑j=1s-1‍aijf(xn+cjh,Yj),i=1,…,s.yn+1=2yn-yn-1+h2∑i=1s‍bif(xn+cih,Yi),
with c1=-1, c2=0, and can be expressed in the Butcher tableau as
(2.2)cAbT=-1000⋯00000⋯0c3a31a320⋯0⋮⋮⋮⋱⋱⋮csas1as2⋯ass-10b1b2⋯bs-1bs
or equivalently by the triplet (c,A,b). The order conditions for this kind of Numerov-type methods can be found in [23]. A favorable property of these methods is that they need only s-1 function evaluations at each step. According to the idea of Lambert and Watson [24], the stability of the Numerov-type methods can be analyzed by the test equation
(2.3)y′′(x)=-λ2y(x).
Applying the explicit Numerov-type method (2.1) to (2.3), we obtain the following relation
(2.4)yn+1-S(H2)yn+P(H2)yn-1=0,H=λh
with
(2.5)S(H2)=2-H2bT(I+H2A)-1(e+c),P(H2)=1-H2bT(I+H2A)-1c.
Since the coefficient matrix A satisfies As-1=0, we have
(2.6)S(H2)=2-H2bT(e+c)+H4bTA(e+c)-⋯+(-1)s-1H2s-2bTAs-2(e+c),P(H2)=1-H2bTc+H4bTAc-⋯+(-1)s-1H2s-2bTAs-2c,
where e is the s-dimensional vector of units, and the vectors c,b, and matrix A are defined by the Butcher tableau of the scheme (2.1). The polynomial
(2.7)ζ2-S(H2)ζ+P(H2)
is called the stability polynomial of the difference equation (2.4) associated with the method (2.1).Definition 2.1.
The method (2.1) is said to have a periodicity intervalIp=(0,H02) if the roots of its stability polynomial, ξ1,2(H) are conjugate complex and satisy |ξ1,2(H)|=1, for all H∈(0,H02). The method is called P-stable if the periodicity interval Ip=(0,∞).The existence of a nonempty periodicity intervalIp=(0,H02) is equivalent to the fact that there exists a nonzero real number H0 such that the coefficients of the stability polynomial satisfy the conditions
(2.8)P(H2)≡1,|S(H2)|<2,∀H∈(0,H02).
Therefore, the P-stable conditions are given by
(2.9)P(H2)≡1,|S(H2)|<2,∀H>0.
An interval Is=(0,H12) is called the interval of absolute stability of the method (2.1) if (see [4])
(2.10)P(H2)<1,|S(H2)|<P(H2)+1,H∈(0,H12).Definition 2.2.
The quantities (see [25])
(2.11)ϕ(H)=H-arccos(S(H2)2P(H2)),d(H)=1-P(H2),
are called the phase-lag and the dissipation, respectively. Accordingly, the method is said to be dispersive of orderq and dissipative of orderp if ϕ(H)=𝒪(Hq+1), and d(H)=𝒪(Hp+1), respectively. If ϕ(H)=0, the method is said to be phase fitted or zero-dispersive, and if d(H)=0, the method is said to be amplification fitted or zero-dissipative.


## 3. Construction of the New Methods

In this section, we construct explicit Numerov-type methods which are both phase fitted and amplification fitted. The construction is based on the following Butcher tableau of coefficients(3.1)-10000000002528132543904357754390400-235167443312538311115625-138666088281250b1b2b3b4.
If we choose {b1,b2,b3,b4}={173/1908,2791/3450,307328/3056775,125/636732}, the classical Numerov-type method is recovered in [26]. In the following, we derive the frequency-dependent weights bi,i=1,…,4 of a one-frequency method and a two-frequency method by phase fitting and amplification fitting.


### 3.1. A One-Frequency Numerov-Type Method

Applying the method (3.1) to the test equation
(3.2)y′′(x)=-ω12y(x),
we obtain the following expressions for S and P(3.3)S(u)=2-b2u2+53b3u2(-1568+675u2)4390+3b4u2(42000-83720u2+159183u4)3500,P(u)=1+b1u2-25b3u2(1568+53u2)43904-23b4u2(-7000+23520u2+769u4),u=ω1h.
Nullifying the dispersion ϕ(u) and the dissipation d(u) in Definition 2.2 for the test equation (3.2) yields
(3.4)S(u)=2cos(u),P(u)=1.
Choosing b1=173/1908,b2=2791/3450, we solve the linear systems (3.4) for the b-values as follows
(3.5)b3=-196(5119800000-19455100000u2+6757940320u4-460960901u6+731400(23520u2+769u4-7000)cos(u))/3056775Mu2,b4=625(-64915200+30302600u2-1798287u4+41400(1568+53u2)cos(u))636732Mu2,
where M=(-19000+952935u2+11707u4).It is easy to verify the above coefficients satisfy the algebraic fifth-order conditions for the Numerov-type methods with one frequency (see [23]):
(3.6)bTe=1-173940800u6+𝒪(u8),bTc=0+110080u4+𝒪(u6),bTc2=16-173470400u2+𝒪(u4),bT(Ae)=112-3792822400u2+𝒪(u4),bTc3=0+349861197568000u2+𝒪(u4),bT(c·(Ae))=112-277201395136000u2+𝒪(u4),bT(Ac)=0+110087395136000u2+𝒪(u4),bTc4=115-744524539219840000u2+𝒪(u4),bT(c2·Ae)=130-17437681955319040000u2+𝒪(u4),bT(c·(Ac))=-160-7106165355319040000u2+𝒪(u4),bT(Ae·Ae)=7120-135568679110638080000u2+𝒪(u4),bT(Ac2)=1180+153119658560000u2+𝒪(u4),bT(A2e)=1360+2390000u2+𝒪(u4).
This one-frequency phase fitted and amplification fitted fifth-order method is denoted as PFAFI.


### 3.2. A Two-Frequency Numerov-Type Method

For the method (3.1), we only need the weights bi(i=1,…,4). In order to solve the oscillatory problems with two frequencies, we consider another test equation
(3.7)y′′(x)=-ω22y(x),ω2≠ω1.
With the same idea stated above, we find the following relation
(3.8)S(v)=2cos(v),P(v)=1,v=ω2h.
Solving (3.4) and (3.8) gives
(3.9)b1=(56(865+28v2)+u2(1568+53v2))(u2-v2+v2cos(u)-u2cos(v))M,b2=(-1107120(u2-v2)+446936(u4-v4)+15525u2v2(u2-v2)+v2(23u2(-1568+675v2)+56(-19770+7981v2))cos(u)-u2(23u2(19432+675v2)-112(9885+322v2))cos(v))/N,b3=-1568(23520+769(u2+v2))(v2-u2-v2cos(u)+u2cos(v))T,b4=250(v2-u2-v2cos(u)+u2cos(v))53061u2v2(u2-v2).
In which M=22260u2v2(u2-v2),N=241500u2v2(u2-v2),T=15283875u2v2(u2-v2).It is also easy to verify the algebraic fifth order conditions for the Numerov-type methods with two frequencies (see [23] for reference)
(3.10)bTe=1-13(ω12ω22(ω12+ω22))h6302400+𝒪(h8),bTc=0-(ω12ω22)h410080+𝒪(h6),bTc2=16-13(ω12+ω22)2h4151200+𝒪(h6),bT(Ae)=112-(13ω14+41ω12ω22+13ω24)h4302400+𝒪(h6),bTc3=0-(ω12+ω22)h21680+𝒪(h4),bT(c·(Ae))=112-(ω12+ω22)h23360+𝒪(h4),bT(Ac)=0-(ω12+ω22)h210080+𝒪(h4),bTc4=115+659(ω12+ω22)h2235200+𝒪(h4),bT(c2·Ae)=130+173(ω12+ω22)h2156800+𝒪(h4),bT(c·(Ac))=-160+659(ω12+ω22)h21411200+𝒪(h4),bT(Ae·Ae)=7120+379(ω12+ω22)h2940800+𝒪(h4),bT(Ac2)=1180-13(ω12+ω22)h2151200+𝒪(h4),bT(A2e)=1360-(ω12+ω22)h210800+𝒪(h4).
The two-frequency fifth-order Numerov-type method with the weights in (3.9) is denoted as PFAFII.


## 3.1. A One-Frequency Numerov-Type Method

Applying the method (3.1) to the test equation
(3.2)y′′(x)=-ω12y(x),
we obtain the following expressions for S and P(3.3)S(u)=2-b2u2+53b3u2(-1568+675u2)4390+3b4u2(42000-83720u2+159183u4)3500,P(u)=1+b1u2-25b3u2(1568+53u2)43904-23b4u2(-7000+23520u2+769u4),u=ω1h.
Nullifying the dispersion ϕ(u) and the dissipation d(u) in Definition 2.2 for the test equation (3.2) yields
(3.4)S(u)=2cos(u),P(u)=1.
Choosing b1=173/1908,b2=2791/3450, we solve the linear systems (3.4) for the b-values as follows
(3.5)b3=-196(5119800000-19455100000u2+6757940320u4-460960901u6+731400(23520u2+769u4-7000)cos(u))/3056775Mu2,b4=625(-64915200+30302600u2-1798287u4+41400(1568+53u2)cos(u))636732Mu2,
where M=(-19000+952935u2+11707u4).It is easy to verify the above coefficients satisfy the algebraic fifth-order conditions for the Numerov-type methods with one frequency (see [23]):
(3.6)bTe=1-173940800u6+𝒪(u8),bTc=0+110080u4+𝒪(u6),bTc2=16-173470400u2+𝒪(u4),bT(Ae)=112-3792822400u2+𝒪(u4),bTc3=0+349861197568000u2+𝒪(u4),bT(c·(Ae))=112-277201395136000u2+𝒪(u4),bT(Ac)=0+110087395136000u2+𝒪(u4),bTc4=115-744524539219840000u2+𝒪(u4),bT(c2·Ae)=130-17437681955319040000u2+𝒪(u4),bT(c·(Ac))=-160-7106165355319040000u2+𝒪(u4),bT(Ae·Ae)=7120-135568679110638080000u2+𝒪(u4),bT(Ac2)=1180+153119658560000u2+𝒪(u4),bT(A2e)=1360+2390000u2+𝒪(u4).
This one-frequency phase fitted and amplification fitted fifth-order method is denoted as PFAFI.


## 3.2. A Two-Frequency Numerov-Type Method

For the method (3.1), we only need the weights bi(i=1,…,4). In order to solve the oscillatory problems with two frequencies, we consider another test equation
(3.7)y′′(x)=-ω22y(x),ω2≠ω1.
With the same idea stated above, we find the following relation
(3.8)S(v)=2cos(v),P(v)=1,v=ω2h.
Solving (3.4) and (3.8) gives
(3.9)b1=(56(865+28v2)+u2(1568+53v2))(u2-v2+v2cos(u)-u2cos(v))M,b2=(-1107120(u2-v2)+446936(u4-v4)+15525u2v2(u2-v2)+v2(23u2(-1568+675v2)+56(-19770+7981v2))cos(u)-u2(23u2(19432+675v2)-112(9885+322v2))cos(v))/N,b3=-1568(23520+769(u2+v2))(v2-u2-v2cos(u)+u2cos(v))T,b4=250(v2-u2-v2cos(u)+u2cos(v))53061u2v2(u2-v2).
In which M=22260u2v2(u2-v2),N=241500u2v2(u2-v2),T=15283875u2v2(u2-v2).It is also easy to verify the algebraic fifth order conditions for the Numerov-type methods with two frequencies (see [23] for reference)
(3.10)bTe=1-13(ω12ω22(ω12+ω22))h6302400+𝒪(h8),bTc=0-(ω12ω22)h410080+𝒪(h6),bTc2=16-13(ω12+ω22)2h4151200+𝒪(h6),bT(Ae)=112-(13ω14+41ω12ω22+13ω24)h4302400+𝒪(h6),bTc3=0-(ω12+ω22)h21680+𝒪(h4),bT(c·(Ae))=112-(ω12+ω22)h23360+𝒪(h4),bT(Ac)=0-(ω12+ω22)h210080+𝒪(h4),bTc4=115+659(ω12+ω22)h2235200+𝒪(h4),bT(c2·Ae)=130+173(ω12+ω22)h2156800+𝒪(h4),bT(c·(Ac))=-160+659(ω12+ω22)h21411200+𝒪(h4),bT(Ae·Ae)=7120+379(ω12+ω22)h2940800+𝒪(h4),bT(Ac2)=1180-13(ω12+ω22)h2151200+𝒪(h4),bT(A2e)=1360-(ω12+ω22)h210800+𝒪(h4).
The two-frequency fifth-order Numerov-type method with the weights in (3.9) is denoted as PFAFII.


## 4. Stability and Phase Analysis

In this section, we analyze the stability and phase properties of the two new methods derived in Section3. The theory of Lambert and Watson was reconsidered by Coleman and Ixaru [27].


### 4.1. Stability and Phase Analysis of PFAFI

We apply a four-stage one-frequency Numerov-type method to the test equation (2.3) and obtain the functions P and S in (2.4) as follows
(4.1)P(H2,u)=1-H2b(u)Tc+H4b(u)TAc-H6b(u)TA2c,S(H2,u)=2-H2b(u)T(e+c)+H4b(u)TA(e+c)-H6b(u)TA2(e+c).
The stability and the phase properties depend on the bivariate functions S(H2,u) and P(H2,u). In this situation, the interval of absolute stability is replaced by a two-dimensional region.Definition 4.1.
For the Numerov-type method withS(H2,u) and P(H2,u) given by (4.1), the quantities
(4.2)ϕ(H,u)=H-arccos(S(H2,u)2P(H2,u)),d(H,u)=1-P(H2,u)
are called the phase-lag and the dissipation, respectively. If
(4.3)ϕ(H,u)=𝒪(Hq+1),d(H,u)=𝒪(Hp+1),
then the method is said to be dissipative of orderq and dissipative of orderp, respectively.For convenience of analyzing the phase-lag and the dissipation, we denote the ratior=u/H. Then the phase-lag and the dissipation of the new method PFAFI are given by
(4.4)ϕ=-13(r2(r2-1))604800H7+𝒪(H9),d=r2-120160H6+𝒪(H8).
Thus the method PFAFI is dispersive of order six and dissipative of order five.Definition 4.2.
In theH-u plane, the region
(4.5)Ω={(H,u)∣H>0,u>0,P(H2,u)<1,|S(H2,u)|<P(H2,u)+1}
is said to be the region of absolute stability of the Numerov-type method PFAFI. Any closed curve defined by
(4.6)P(H2,u)=1or|S(H2,u)|=P(H2,u)+1
is a stability boundary of the method.In Figure1, we plot the region of the absolute stability of the new method PFAFI derived in Section 3.1 (a similar discussion can be found in [27, 28]).Figure 1
Stability region for the method PFAFI.It should be noted that when the exact value of the frequency of the problem is known (i.e., the fitted frequencyω equals the test frequency λ), then u=H and we have that S(H2,H)=2cos(H) and P(H2,H)=1. Consequently, when the main frequency is known the interval of periodicity is (0,∞) except for H=2π,4π,…. So the new method PFAFI is almost P-stable.


### 4.2. Stability and Phase Analysis of Phase Fitted Numerov-Type Method with Two Frequencies

We apply a four-stage two-frequency Numerov-type method to the test equation (2.3) and obtain the functions P and S in (2.4) as follows
(4.7)P(H2,u,v)=1-H2b(u,v)Tc+H4b(u,v)TAc-H6b(u,v)TA2c,S(H2,u,v)=2-H2b(u,v)T(e+c)+H4b(u,v)TA(e+c)-H6b(u,v)TA2(e+c).Definition 4.3.
For the Numerov-type method withS(H2,u,v) and P(H2,u,v) given by (4.7), the quantities
(4.8)ϕ(H,u,v)=H-arccos(S(H2,u,v)2P(H2,u,v)),d(H,u,v)=1-P(H2,u,v)
are called the phase-lag and the dissipation, respectively. If
(4.9)ϕ(H,u,v)=𝒪(Hq+1),d(H,u,v)=𝒪(Hp+1),
then this method is said to be dispersive of orderq and dissipative of orderp, respectively.Denoter=u/H and k=v/u. Then the phase-lag and the dissipation of the new method PFAFII are given by
(4.10)ϕ=13r2(r2-1)(1+k2)(k2r2-1)604800H7+𝒪(H9),d=-(r2-1)(r2k2-1)20160H6+𝒪(H8).
Thus the new method PFAFII is dispersive of order six and dissipative of order five.Following the main idea of Coleman and Ixaru in [27], the absolute stability of the two-frequency method PFAFII can be described by a three-dimensional region.Definition 4.4.
A three-dimensional regionΩ in the H-u-v space
(4.11)V={(H,u,v)∣P(H2,u,v)<1,|S(H2,u,v)|<P(H2,u,v)+1}
is called the region of absolute stability of a Numerov-type method with the S(H2,u,v) and P(H2,u,v) given by (4.7), and any closed curve defined by
(4.12)P(H2,u,v)=1or|S(H2,u,v)|=P(H2,u,v)+1
is a stability boundary of the method.In Figure2, we plot the region of absolute stability for the method PFAFII, and Figure 3 shows the projection of the three-dimensional region onto the u-H plane. (A similar discussion can be found in [22, 29].)Figure 2
Stability region for the method PFAFII.Figure 3
The projection of the three-dimensional region of PFAFII on theu-H plane.It should be noted that, for the new method PFAFII and for a problem with one frequency, when the exact value of the frequency of the problem is known, that is,ω1=ω2=λ, then u=v=H and we have that S(H2,H,H)=2cos(H) and P(H2,H,H)=1. Consequently, the interval of periodicity is (0,∞) with the points H=2π,4π,… deleted. That is, the new method PFAFII is also almost P-stable.


## 4.1. Stability and Phase Analysis of PFAFI

We apply a four-stage one-frequency Numerov-type method to the test equation (2.3) and obtain the functions P and S in (2.4) as follows
(4.1)P(H2,u)=1-H2b(u)Tc+H4b(u)TAc-H6b(u)TA2c,S(H2,u)=2-H2b(u)T(e+c)+H4b(u)TA(e+c)-H6b(u)TA2(e+c).
The stability and the phase properties depend on the bivariate functions S(H2,u) and P(H2,u). In this situation, the interval of absolute stability is replaced by a two-dimensional region.Definition 4.1.
For the Numerov-type method withS(H2,u) and P(H2,u) given by (4.1), the quantities
(4.2)ϕ(H,u)=H-arccos(S(H2,u)2P(H2,u)),d(H,u)=1-P(H2,u)
are called the phase-lag and the dissipation, respectively. If
(4.3)ϕ(H,u)=𝒪(Hq+1),d(H,u)=𝒪(Hp+1),
then the method is said to be dissipative of orderq and dissipative of orderp, respectively.For convenience of analyzing the phase-lag and the dissipation, we denote the ratior=u/H. Then the phase-lag and the dissipation of the new method PFAFI are given by
(4.4)ϕ=-13(r2(r2-1))604800H7+𝒪(H9),d=r2-120160H6+𝒪(H8).
Thus the method PFAFI is dispersive of order six and dissipative of order five.Definition 4.2.
In theH-u plane, the region
(4.5)Ω={(H,u)∣H>0,u>0,P(H2,u)<1,|S(H2,u)|<P(H2,u)+1}
is said to be the region of absolute stability of the Numerov-type method PFAFI. Any closed curve defined by
(4.6)P(H2,u)=1or|S(H2,u)|=P(H2,u)+1
is a stability boundary of the method.In Figure1, we plot the region of the absolute stability of the new method PFAFI derived in Section 3.1 (a similar discussion can be found in [27, 28]).Figure 1
Stability region for the method PFAFI.It should be noted that when the exact value of the frequency of the problem is known (i.e., the fitted frequencyω equals the test frequency λ), then u=H and we have that S(H2,H)=2cos(H) and P(H2,H)=1. Consequently, when the main frequency is known the interval of periodicity is (0,∞) except for H=2π,4π,…. So the new method PFAFI is almost P-stable.


## 4.2. Stability and Phase Analysis of Phase Fitted Numerov-Type Method with Two Frequencies

We apply a four-stage two-frequency Numerov-type method to the test equation (2.3) and obtain the functions P and S in (2.4) as follows
(4.7)P(H2,u,v)=1-H2b(u,v)Tc+H4b(u,v)TAc-H6b(u,v)TA2c,S(H2,u,v)=2-H2b(u,v)T(e+c)+H4b(u,v)TA(e+c)-H6b(u,v)TA2(e+c).Definition 4.3.
For the Numerov-type method withS(H2,u,v) and P(H2,u,v) given by (4.7), the quantities
(4.8)ϕ(H,u,v)=H-arccos(S(H2,u,v)2P(H2,u,v)),d(H,u,v)=1-P(H2,u,v)
are called the phase-lag and the dissipation, respectively. If
(4.9)ϕ(H,u,v)=𝒪(Hq+1),d(H,u,v)=𝒪(Hp+1),
then this method is said to be dispersive of orderq and dissipative of orderp, respectively.Denoter=u/H and k=v/u. Then the phase-lag and the dissipation of the new method PFAFII are given by
(4.10)ϕ=13r2(r2-1)(1+k2)(k2r2-1)604800H7+𝒪(H9),d=-(r2-1)(r2k2-1)20160H6+𝒪(H8).
Thus the new method PFAFII is dispersive of order six and dissipative of order five.Following the main idea of Coleman and Ixaru in [27], the absolute stability of the two-frequency method PFAFII can be described by a three-dimensional region.Definition 4.4.
A three-dimensional regionΩ in the H-u-v space
(4.11)V={(H,u,v)∣P(H2,u,v)<1,|S(H2,u,v)|<P(H2,u,v)+1}
is called the region of absolute stability of a Numerov-type method with the S(H2,u,v) and P(H2,u,v) given by (4.7), and any closed curve defined by
(4.12)P(H2,u,v)=1or|S(H2,u,v)|=P(H2,u,v)+1
is a stability boundary of the method.In Figure2, we plot the region of absolute stability for the method PFAFII, and Figure 3 shows the projection of the three-dimensional region onto the u-H plane. (A similar discussion can be found in [22, 29].)Figure 2
Stability region for the method PFAFII.Figure 3
The projection of the three-dimensional region of PFAFII on theu-H plane.It should be noted that, for the new method PFAFII and for a problem with one frequency, when the exact value of the frequency of the problem is known, that is,ω1=ω2=λ, then u=v=H and we have that S(H2,H,H)=2cos(H) and P(H2,H,H)=1. Consequently, the interval of periodicity is (0,∞) with the points H=2π,4π,… deleted. That is, the new method PFAFII is also almost P-stable.


## 5.  Numerical Experiments

In this section, we compare the numerical performance of the new methods with some well-known methods from the literature. We select three periodic IVPs with two frequencies. The methods we choose for comparison are listed as follows.(i)
EFRKN4F: the trigonometrically fitted RKN method of order 4 with the FSAL property given in [2].(ii)
VPHA: the phase fitted and amplification fitted Numerov-type method with one frequency given in [4].(iii)
TFNTONE: the trigonometrically fitted Numerov-type method obtained by Fang and Wu in [3].(iv)
PFAFI: the phase fitted and amplification fitted Numerov-type method with one frequency obtained in this paper.(v)
PFAFII: the phase fitted and amplification fitted Numerov-type method with two frequencies obtained in this paper.Problem 1.
We consider the inhomogeneous equation(5.1)y′′(x)+100y(x)=99sin(x),y(0)=1,y′(0)=11,
whose exact solution is given by
(5.2)y(x)=cos(10x)+sin(10x)+sin(x).
In our test we choose ω1=10 for VPHA, TFNTONE, EFRKN4F and PFAFI, and  ω1=10,ω2=1 for PFAFII. The problem is integrated on the interval [0,1000] and the numerical results are presented in Table 1.Table 1
Comparison of the end-point global errors and the consumed CPU time in seconds (in the brackets) produced by the five methods for Problem1.
h
EFRKN4F
VPHA
TFNTONE
PFAFI
PFAFII
1
NaN (0.09)
674.5 (0.08)
63.2 (0.06)
0.27 (0.06)
7.9e − 3 (0.06)
1
/
2
0.17 (1.19)
12.7 (0.13)
1.99 (0.11)
0.04 (0.13)
1.3e − 3 (0.11)
1
/
4
4.6e − 3 (0.34)
2.3e − 2 (0.28)
2.7e − 3 (0.23)
1.5e − 4 (0.23)
1.1e − 5 (0.23)
1
/
8
1.5e − 4 (0.71)
2.4e − 4 (0.53)
2.8e − 5 (0.47)
2.3e − 6 (0.48)
1.9e − 7 (0.48)
1
/
16
0.7e − 6 (1.39)
3.4e − 6 (1.13)
1.43e − 7 (0.95)
1.3
×
10
-
7 (0.95)
4.9e − 9 (0.95)Problem 2.
Consider the famous Kramarz system (see [30])
(5.3)y′′=2498y+4998z,y(0)=2,y′(0)=0,z′′=-2499y-4999z,z(0)=-1,z′(0)=0,
whose exact solution is
(5.4)y(x)=2cos(x),z(x)=-cos(x).
The eigenvalues of the matrix of coefficients of the equations for y′′ and z′′ are −1 and −2500. For this stiff equation, we choose ω1=1 for VPHA, TFNTONE, EFRKN4F and PFAFI, and ω1=50,ω2=1 for PFAFII. The problem is integrated on the interval [0,1000] and the numerical results are presented in Table 2.Table 2
Comparison of the end-point global errors and the consumed CPU time in seconds (in the brackets) produced by the five methods for Problem2.
h
EFRKN4F
VPHA
TFNTONE
PFAFI
PFAFII
8
NaN (0.02)
NaN (0.03)
NaN (0.02)
NaN (0.02)
6.7e − 12 (0.02)
4
NaN (0.05)
NaN (0.05)
NaN (0.03)
NaN (0.03)
2.5e − 12 (0.03)
2
NaN (0.09)
NaN (0.09)
NaN (0.06)
NaN (0.05)
2.0e − 11 (0.06)
1
NaN (0.17)
NaN (0.19)
NaN (0.13)
NaN (0.11)
1.8e − 11 (0.09)
1
/
2
NaN (0.34)
NaN (0.38)
NaN (0.25)
NaN (0.23)
1.0e − 11 (0.20)
1
/
4
NaN (0.70)
NaN (0.77)
NaN (0.5)
NaN (0.47)
8.0e − 12 (0.40)
1
/
8
NaN (1.42)
NaN (1.53)
NaN (1.0)
NaN (0.92)
1.7e − 11 (0.82)
1
/
16
1.2e − 12 (2.56)
1.6e − 11 (2.70)
2.3e − 12(1.75)
NaN (1.86)
1.8e − 11 (1.65)Problem 3.
Consider the linear stiff system in [31]
(5.5)y′′(t)=-12(ω2+1ω2-1ω2-1ω2+1)y(t),y(0)=(1+ϵ,ϵ-1)T,y′(0)=(ϵω+1,ϵω-1)T,
whose analytic solution is given by
(5.6)y(t)=(cos(t)+sin(t)-cos(t)-sin(t))+ϵ(cos(ωt)+sin(ωt)cos(ωt)+sin(ωt)).
In our test we choose ω1=ω for VPHA, TFNTONE, EFRKN4F and PFAFI, and ω1=ω,ω2=1 for PFAFII. The problem is integrated on the interval [0,10] and the numerical results are presented in Table 3 with ϵ=10-7,ω=103, and in Table 4 with ϵ=10-7,ω=105.Table 3
Comparison of the end-point global errors and the consumed CPU time in seconds (in the brackets) produced by the five methods for Problem3 with ϵ=10-7,ω=103.
h
EFRKN4F
VPHA
TFNTONE
PFAFI
PFAFII
2
NaN (0)
NaN (0)
NaN (0)
NaN (0)
2.6e − 11 (0)
1
NaN (0)
NaN (0)
0.01 (0)
NaN (0)
3.0e − 11 (0)
1
/
2
NaN (0.015)
NaN (0.015)
0.35 (0.015)
1.38 (0.015)
2.9e − 12 (0.015)
1
/
4
NaN (0.016)
NaN (0.016)
0.025 (0.016)
1.43 (0.016)
4.7e − 11 (0.016)
1
/
8
NaN (0.016)
NaN (0.016)
3.3e − 3 (0.016)
1.56 (0.016)
1.4e − 11 (0.015)
1
/
16
NaN (0.031)
NaN (0.016)
2.6e − 3 (0.016)
1.68 (0.016)
1.5e − 11 (0.031)
1
/
32
NaN (0.047)
NaN (0.031)
8.6e − 4 (0.031)
1.42 (0.031)
1.3e − 12 (0.031)Table 4
Comparison of the end-point global errors and the consumed CPU time in seconds (in the brackets) produced by the five methods for Problem3 with ϵ=10-7,ω=105.
h
EFRKN4F
VPHA
TFNTONE
PFAFI
PFAFII
1/8
NaN (0.02)
NaN (0.02)
0.29 (0.01)
NaN (0.01)
5.7e − 7 (0.02)
1/16
NaN (0.03)
NaN (0.03)
0.05 (0.02)
1.76 (0.02)
8.3e − 8 (0.03)
1/32
NaN (0.06)
NaN (0.06)
8.8e − 3 (0.03)
1.91 (0.03)
7.5e − 8 (0.05)
1/64
NaN (0.11)
NaN (0.11)
8.6e − 3 (0.06)
2.01 (0.05)
4.1e − 8 (0.08)
1/128
NaN (0.20)
NaN (0.23)
1.8e − 3 (0.14)
2.07 (0.08)
3.3e − 8 (0.14)
1/256
NaN (0.43)
NaN (0.45)
9.0e − 5 (0.27)
2.10 (0.17)
1.1e − 8 (0.27)
1/512
NaN (0.89)
NaN (0.93)
1.8e − 6 (0.53)
2.11 (0.36)
4.6e − 9 (0.53)Problem 4.
Consider the nonlinear system in [32]
(5.7)y′′(t)=(-100-λ2)y(t)+k(11)(∥y(t)∥22+cos2(t)-1),y(0)=(0,0)T,y′(0)=(1,0)T,
where λ=103,k=0.01. The analytic solution of the problem is
(5.8)y(t)=(sin(t),0)T.In our test we chooseω=λ for VPHA, TFNTONE, EFRKN4F and PFAFI, and ω1=λ,ω2=1 for PFAFII. The problem is integrated on the integral interval [0,1000] and the numerical results are stated in Table 5.Table 5
Comparison of the end-point global errors and the consumed CPU time in seconds (in the brackets) produced by the five methods for Problem4.
h
EFRKN4F
VPHA
TFNTONE
PFAFI
PFAFII
1/4
NaN (0.76)
NaN (0.61)
NaN (0.53)
NaN (0.42)
3.45e − 7 (0.43)
1/8
NaN (1.52)
NaN (1.24)
NaN (0.89)
NaN (0.81)
2.14e − 8 (0.89)
1/16
NaN (3.06)
NaN (2.45)
NaN (1.79)
NaN (1.60)
1.33e − 9 (1.78)
1/32
NaN (6.04)
NaN (4.89)
0.008 (3.56)
0.82 (2.5)
3.98e − 11 (3.59)
1/64
NaN (12.07)
NaN (9.73)
0.08 (7.1)
0.82 (5)
4.39e − 12 (7.1)


## 6. Conclusions

In this paper, we investigate Numerov-type methods adapted to oscillatory problems with two frequencies obtained by phase fitting and amplification fitting. The numerical stability and the phase properties of the new method are analyzed. Unlike the case of one-frequency methods, the absolute stability region of a two-frequency method is a region in the three-dimensional space. The results of numerical experiments on four oscillatory test problems with two frequencies, including non-stiff and stiff problems, show that the new methods are superior to some well-known high quality codes proposed in the recent scientific literature.



---

*Source: 742585-2012-10-12.xml*