A New Family of Phase-Fitted and Amplification-Fitted Runge-Kutta Type Methods for Oscillators

Abstract

In order to solve initial value problems of differential equations with oscillatory solutions, this paper improves traditional Runge-Kutta (RK) methods by introducing frequency-depending weights in the update. New practical RK integrators are obtained with the phase-fitting and amplification-fitting conditions and algebraic order conditions. Two of the new methods have updates that are also phase-fitted and amplification-fitted. The linear stability and phase properties of the new methods are examined. The results of numerical experiments on physical and biological problems show the robustness and competence of the new methods compared to some highly efficient integrators in the literature.

1. Introduction

A vast of problems in applied fields, such as elastics, mechanics, astrophysics, electronics, molecular dynamics, ecology, and biochemistry, can be cast into the form of an initial value problem of the system of first-order ordinary differential equations as follows: 𝑦^′=𝑓𝑥,𝑦,𝑦𝑥₀=𝑦₀, (1.1) where 𝑦∈ℝ^𝑑,𝑓∶ℝ×ℝ^𝑑→ℝ^𝑑 is a sufficiently smooth function. Traditionally, the problem (1.1) has been solved numerically by Runge-Kutta methods (RK) or linear multistep methods (LMMs) (see [1–3]). In many applications, the solution to the problem (1.1) is oscillatory or periodic. But the general-purpose methods do not take into account the oscillatory feature of the problem and the numerical results they produce are generally not satisfactory.

Recently, some authors have proposed to adapt traditional integrators to the oscillatory character of the solution to the problem (1.1) (see [4–7]). Bettis [8] constructs a three-stage method and a four-stage method which solve the equation 𝑦^′=i𝜔𝑦(i²=−1) without truncation error. Paternoster [9] develops a class of implicit methods of Runge-Kutta (RK) and Runge-Kutta-Nyström (RKN) types by the trigonometric fitting technique. For oscillatory problems which can be put in the form of a second-order equation 𝑦^′′=𝑓(𝑥,𝑦), Franco [10] improves the update of the classical RKN methods and proposes a family of explicit RKN methods adapted to perturbed oscillators (ARKN) and a class of explicit adapted RK methods in [11]. Anastassi and Simos [12] construct a phase-fitted and amplification-fitted RK method of “almost” order five. Van de Vyver [13] investigates phase-fitted and amplification-fitted two step hybrid methods (FTSH). For other important work on frequency dependent integrators for general second-order oscillatory equations, the reader is referred to [14–18]. Many authors focus on effective numerical integration of specific categories of oscillatory problems. For example, Vigo-Aguiar and Simos [18] construct an exponentially fitted and trigonometrically fitted method for orbital problems. The papers [19–21] have designed highly efficient integrators for the Schrödinger equation.

Based on the previous work, we consider, in this paper, phase-fitted and amplification-fitted RK type integrators whose coefficients in the update depend on the product of the fitting frequency and the step size. In Section 2, we present a result on order conditions for frequency-depending modified Runge-Kutta type methods. Section 3 introduces the notion of phase-fitted and amplification-fitted RK type methods (FRK) and derives the phase-fitting and amplification-fitting conditions. In Section 4 two FRK methods of order four and two FRK methods of order five are constructed. Their error coefficients and error constants are also calculated. The linear stability and phase properties of the new FRK methods are analyzed in Section 5. In Section 6, numerical experiments are carried out to illustrate the effectiveness and superiority of our new methods compared to two well-known highly efficient integrators we have chosen from the recent literature. Section 7 is devoted to conclusive comments.

2. Order Conditions for RK Type Methods with Frequency-Dependent Weights

Assume that the principal frequency of the problem (1.1) is known or can be accurately estimated in advance. This estimated frequency is denoted by 𝜔, which is also called the fitting frequency. We consider the following 𝑠-stage modified Runge-Kutta method: 𝑘_𝑖=𝑓𝑥₀+𝑐_𝑖ℎ,𝑦₀+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑘_𝑗𝑦,𝑖=1,…,𝑠,₁=𝑦₀+ℎ^𝑠_𝑖=1𝑏_𝑖(𝜈)𝑘_𝑖, (2.1) where ℎ is the step size, 𝑐_𝑖,𝑎_𝑖𝑗,𝑖,𝑗=1,…,𝑠 are real numbers, and 𝑏_𝑖(𝜈),𝑖=1,…,𝑠 are even functions of 𝜈=ℎ𝜔. The scheme (2.1) can be represented by the Butcher tableau as follows: 𝑐𝐴𝑏^𝑇=𝑐₁𝑎₁₁…𝑎_1𝑠⋮𝑐⋮⋱⋮_𝑠𝑎_𝑠1…𝑎_𝑠𝑠𝑏₁(𝜈)…𝑏_𝑠(𝜈), (2.2) or simply by (𝐴,𝑏,𝑐). Conventionally, we assume that ^𝑠_𝑗=1𝑎_𝑖𝑗=𝑐_𝑖,𝑖=1,…,𝑠. (2.3)

As explained in [3], the order conditions for the modified RK method (2.1) can be derived by just considering the autonomous equation 𝑦^′=𝑓(𝑦). Using Butcher's rooted tree theory in [1], the exact solution to the problem (1.1) and the numerical solution produced by the modified RK type method (2.1) have the 𝐵-series expressions as follows: 𝑦𝑥₀+ℎ=^∞_𝑗=0ℎ^𝑗𝑗!_{𝜌(𝑡)=𝑗}𝛼(𝑡)ℱ(𝑡)𝑦₀,𝑦₁=^∞_𝑗=0ℎ^𝑗𝑗!_{𝜌(𝑡)=𝑗}𝛼(𝑡)𝛾(𝑡)𝑏^𝑇Φ(𝑡)ℱ(𝑡)𝑦₀, (2.4) where the trees 𝑡, the functions 𝜌(𝑡) (order), 𝛼(𝑡) (number of equivalent trees), 𝛾(𝑡) (density), Φ(𝑡)=(Φ₁(𝑡),…,Φ_𝑠(𝑡))^𝑇 (the 𝑠×1 vector of elementary weights), and the elementary differentials ℱ(𝑡)(𝑦₀) are defined in [3]. Then the local truncation error has the following series expansion: LTE=^∞_𝑗=0ℎ^𝑗𝑗!_{𝜌(𝑡)=𝑗}𝛼(𝑡)1−𝛾(𝑡)𝑏^𝑇Φ(𝑡)ℱ(𝑡)𝑦₀. (2.5) If for any (𝑝+1)th differentiable function 𝑓(𝑦), when the scheme (2.1) is applied to the problem (1.1), the local truncation error LTE=𝑦(𝑥₀+ℎ)−𝑦₁=𝒪(ℎ^𝑝+1), then the method (2.1) is said to have (algebraic) order 𝑝.

Theorem 2.1 (see Franco [11]). The modified RK type method (2.1) has order 𝑝 if and only if the following conditions are satisfied: 𝑏^𝑇Φ(𝑡)=1𝛾(𝑡)+𝒪𝜈^{𝑝−𝜌(𝑡)+1},𝜌(𝑡)=1,2,…,𝑝. (2.6)

If the method (2.1) is of order 𝑝, then we have ℎLTE=^𝑝+1𝑝+1!_{𝜌(𝑡)=𝑝+1}𝛼(𝑡)1−𝛾(𝑡)𝑏^𝑇Φ(𝑡)ℱ(𝑡)𝑦₀+𝒪𝜈^𝑝+2ℎ∶=^𝑝+1𝑝+1!_{𝜌(𝑡)=𝑝+1}𝛼(𝑡)𝜖(𝑡)ℱ(𝑡)𝑦₀+𝒪𝜈^𝑝+2, (2.7) where 𝜖(𝑡)=1−𝛾(𝑡)𝑏^𝑇Φ(𝑡) is called the error coefficient associated with the tree 𝑡 of order 𝑝+1 which is involved in the leading term of the local truncation error. Denote 𝐸𝐶_𝑝+1(𝜈)=_{𝜌(𝑡)=𝑝+1}𝜖(𝑡)^21/2. (2.8) The positive number 𝐸𝐶_𝑝+1=lim_𝜈→0𝐸𝐶_𝑝+1(𝜈) (2.9) is called the error constant of the method.

3. Phase-Fitted and Amplification-Fitted Conditions

For oscillatory problems, as suggested by Paternoster [9] and Van der Houwen and Sommeijer [22], apart from the algebraic order, the analysis of phase lag and dissipation is important. To this end, we consider the following linear scalar equation: 𝑦^′=i𝜔𝑦. (3.1) The exact solution of this equation with the initial value 𝑦(𝑥₀)=𝑦₀ satisfies 𝑦𝑥₀+ℎ=𝑅₀(𝑧)𝑦₀, (3.2) where 𝑅₀(𝑧)=exp(𝑧),𝑧=i𝜈. This means that after a period of time ℎ, the exact solution experiences a phase advance 𝜈=ℎ𝜔 and the amplification remains constant.

An application of the modified RK method (2.1) to (3.1) yields 𝑦₁=𝑅(𝑧)𝑦₀, (3.3) where 𝑅(𝑧)=1+𝑧𝑏^𝑇(𝐼−𝑧𝐴)⁻¹𝑒,𝑧=i𝜈,𝑒=(1,…,1)^𝑇. (3.4) Thus, after a time step ℎ, the numerical solution obtains a phase advance arg𝑅(𝑧) and the amplification factor |𝑅(i𝜈)|. 𝑅(𝑧) is called the stability function of the method (2.1). Denote the real and imaginary part of 𝑅(𝑧) by 𝑈(𝜈) and 𝑉(𝜈), respectively. Then, for small ℎ, we have 𝑈(𝜈)=1−𝜈²𝑏^𝑇𝐴𝑒+𝜈⁴𝑏^𝑇𝐴³𝑒−𝜈⁶𝑏^𝑇𝐴⁵𝑒+𝜈⁸𝑏^𝑇𝐴⁷𝑒𝑉(𝜈)𝑏−⋯,=𝜈^𝑇𝑒−𝜈³𝑏^𝑇𝐴²𝑒+𝜈⁵𝑏^𝑇𝐴⁴𝑒−𝜈⁷𝑏^𝑇𝐴⁶𝑒+⋯. (3.5) For small ℎ, arg𝑅(𝑧)=tan⁻¹(𝑉(𝜈)/𝑈(𝜈)) and |𝑅(𝑧)|=𝑈²(𝜈)+𝑉²(𝜈). The above analysis leads to the following definition.

Definition 3.1 (see [22]). The quantities 𝑃(𝜈)=𝜈−arg𝑅(i𝜈),𝐷(𝜈)=1−|𝑅(i𝜈)| (3.6) are called the phase lag (or dispersion) and the error of amplification factor (or dissipation) of the method, respectively. If 𝑃(𝜈)=𝒪𝜈^𝑞+1,𝐷(𝜈)=𝒪𝜈^𝑝+1, (3.7) then the method is called dispersive of order 𝑞 and dissipative of order 𝑝, respectively. If 𝑃(𝜈)=0,𝐷(𝜈)=0, (3.8) the method is called phase-fitted (or zero-dispersive) and amplification-fitted (or zero-dissipative), respectively. A modified RK method which is phase-fitted and amplification-fitted is called in short an FRK method.

It is interesting to consider the phase properties of the update of the scheme (2.1). Suppose that the internal stages have been exact for the linear equation (3.1), that is, 𝑌_𝑖=exp(i𝑐_𝑖𝜈)𝑦₀, then the update gives 𝑦₁=𝑅_𝑢(𝑧)𝑦₀, (3.9) where 𝑅_𝑢(𝑧)=1+𝑧^𝑠_𝑖=1𝑏_𝑖(𝜈)exp𝑐_𝑖𝑧,𝑧=i𝜈. (3.10) Denote the real and imaginary part of 𝑅_𝑢(𝑧) by 𝑈_𝑢(𝜈) and 𝑉_𝑢(𝜈), respectively. Then, for small ℎ, 𝑈_𝑢(𝜈)=1−𝜈^𝑠_𝑖=1𝑏_𝑖(𝜈)sin𝑐_𝑖𝜈,𝑉_𝑢(𝜈)=𝜈^𝑠_𝑖=1𝑏_𝑖(𝜈)cos𝑐_𝑖𝜈. (3.11)

Definition 3.2. The quantities 𝑃_𝑢(𝜈)=𝜈−arg𝑅_𝑢(i𝜈),𝐷_𝑢(𝜈)=1−𝑅_𝑢(i𝜈) (3.12) are called the phase lag (or dispersion) and the error of amplification factor (or dissipation) of the update of the method, respectively. If 𝑃_𝑢(𝜈)=𝒪𝜈^𝑞+1,𝐷_𝑢(𝜈)=𝒪𝜈^𝑝+1, (3.13) then the update of the method is called dispersive of order 𝑞 and dissipative of order 𝑝, respectively. If 𝑃_𝑢(𝜈)=0,𝐷_𝑢(𝜈)=0, (3.14) the update of the method is called phase-fitted (or zero-dispersive) and amplification-fitted (or zero-dissipative) in the update, respectively.

Generally speaking, a traditional RK method with constant coefficients inevitably carries a nonzero phase lag and a nonzero error of amplification factor when applied to the linear oscillatory equation (3.1). Therefore, they are neither phase-fitted nor amplification-fitted. So does the update. For example, the classical RK method of order four with constant coefficients (denoted as RK4, see [3]) given by 01/21/21/211/211/62/62/61/6 (3.15) has a phase lag and an error of amplification factor as follows: 𝑃(𝜈)=1𝜈120⁵+𝒪𝜈⁷(𝜈)=1,𝐷𝜈144⁶+𝒪𝜈⁸. (3.16) For the update, 𝑃_𝑢(𝜈)1=−𝜈2880⁵+𝒪𝜈⁷,𝐷_𝑢(𝜈)1=−𝜈5760⁶+𝒪𝜈⁸. (3.17) Therefore, both the method and its update are dispersive of order four and dissipative of order five.

The following theorem gives the necessary and sufficient conditions for a modified RK method and its update to be phase-fitted and amplification-fitted, respectively.

Theorem 3.3. (i) The method (2.1) is phase-fitted and amplification-fitted if and only if 𝑈(𝜈)=cos(𝜈),𝑉(𝜈)=sin(𝜈). (3.18)
(ii) The update of the method (2.1) is phase-fitted and amplification-fitted if and only if 𝑈_𝑢(𝜈)=cos(𝜈),𝑉_𝑢(𝜈)=sin(𝜈). (3.19)

The proof of this theorem is immediate.

4. Construction of New Methods

Now we proceed to construct modified RK type methods that are both phase-fitted and amplification-fitted based on the internal coefficients of two classical RK methods. For convenience we restrict ourselves to explicit methods.

4.1. Fourth-Order Methods

Consider a four-stage modified RK method with the following Butcher tableau: 01/21/21/211/21𝑏₁𝑏₂𝑏₃𝑏₄. (4.1) For 𝑠=4, the phase-fitting and amplification-fitting conditions (3.18) have the following form: 11−2𝑏₂(𝜈)+12𝑏₃(𝜈)+𝑏₄(𝜈)𝜈²+14𝑏₄(𝜈)𝜈⁴(𝜈),𝑏=cos₁(𝜈)+𝑏₂(𝜈)+𝑏₃(𝜈)+𝑏₄(𝜈)1𝜈−4𝑏₃(𝜈)+12𝑏₄(𝜈)𝜈³=sin(𝜈). (4.2) On the other hand, the first-order and second-order conditions become 𝑏^𝑇𝑒=𝑏₁(𝜈)+𝑏₂(𝜈)+𝑏₃(𝜈)+𝑏₄(𝜈)𝑏=1,_𝑇1𝑐=2𝑏₂(𝜈)+12𝑏₃(𝜈)+𝑏₄(𝜈)=12. (4.3) Here, for the purpose of deriving the weights of the method, the higher order terms in the order conditions are omitted. We will follow this convention in the sequel. Solving (4.2) and (4.3), we obtain 𝑏₁(𝜈)=2−2+𝜈²+2cos(𝜈)𝜈⁴,𝑏₂(𝜈)=𝜈³−4𝜈+4sin(𝜈)𝜈³,𝑏₃(𝜈)=−4(−2+2cos(𝜈)+𝜈sin(𝜈))𝜈⁴,𝑏₄(𝜈)=𝑏₁(𝜈). (4.4) As 𝜈→0, 𝑏_𝑖(𝜈) have the following Taylor expansions: 𝑏₁(𝜈)=16−1𝜈180²+1𝜈10080⁴−1𝜈907200⁶+1𝜈119750400⁸𝑏+⋯,₂(𝜈)=13+1𝜈30²−1𝜈1260⁴+1𝜈907200⁶−1𝜈9979200⁸𝑏+⋯,₃(𝜈)=13−1𝜈45²+1𝜈1680⁴−1𝜈134400⁶+1𝜈11975040⁸𝑏+⋯,₄(𝜈)=16−1𝜈180²+1𝜈10080⁴−1𝜈907200⁶+1𝜈119750400⁸+⋯. (4.5) By direct calculation, we obtain the following expansions: 𝑏^𝑇𝑐²=−8+4𝜈²+𝜈⁴+8cos(𝜈)4𝜈⁴=13+𝒪𝜈²,𝑏^𝑇(𝜈)𝐴𝑐=𝜈−sin𝜈³=16+𝒪𝜈²,𝑏^𝑇𝑐³=−24+12𝜈²+𝜈⁴+24cos(𝜈)8𝜈⁴=14(𝜈),𝑏+𝒪_𝑇(𝑐⋅𝐴𝑐)=−2+2𝜈²+2cos(𝜈)−𝜈sin(𝜈)2𝜈⁴=18(𝜈),𝑏+𝒪_𝑇𝐴𝑐²=𝜈−sin(𝜈)2𝜈³=1(𝜈),𝑏12+𝒪_𝑇𝐴²𝑐=−2+𝜈²+2cos(𝜈)2𝜈⁴=1(𝜈),24+𝒪 (4.6) where the dot “⋅” between two vectors indicates componentwise product, and 𝑐² and 𝑐³ indicate the componentwise square of the vector. We will follow the convention in the sequel. Thus, the coefficients given by (4.1) and (4.4) satisfy all the conditions of order four in Theorem 2.1. But they cannot satisfy all the conditions of order five. For instance, 𝑏^𝑇𝑐⁴=5−724𝜈1440²+1𝜈11520⁴−1𝜈1036800⁶≠15+𝒪(𝜈). (4.7) Therefore, this method is of order four. The method was originally obtained by Simos [23] and we denote it by Simos4.

Corresponding to the nine fifth-order rooted trees 𝑡_5𝑗,𝑗=1,…,9, the error coefficients of Simos4 are given by 𝑒𝑡₅₁=280−140𝜈²+11𝜈⁴−280cos(𝜈)16𝜈⁴,𝑒𝑡₅₂=30−20𝜈²+2𝜈⁴−30cos(𝜈)+5𝜈sin(𝜈)2𝜈⁴,𝑒𝑡₅₃=30−30𝜈²+4𝜈⁴−30cos(𝜈)+15𝜈sin(𝜈)4𝜈⁴,𝑒𝑡₅₄=1−15−2+𝜈²+2cos(𝜈)𝜈⁴,𝑒𝑡₅₅=10−10𝜈²+𝜈⁴−10cos(𝜈)+5𝜈sin(𝜈)𝜈⁴,𝑒𝑡₅₆=−5𝜈+𝜈³+5sin(𝜈)𝜈³,𝑒𝑡₅₇=1−10−2+𝜈²+2cos(𝜈)𝜈⁴,𝑒𝑡₅₈=1−15−2+𝜈²+2cos(𝜈)𝜈⁴,𝑒𝑡₅₉=1. (4.8) Then for Simos4, 𝐸𝐶₅(𝜈)=√𝐾16𝜈⁴, (4.9) where 𝐾=739200−798400𝜈²+296160𝜈⁴−40200𝜈⁶+2169𝜈⁸+739200cos²(𝜈)+320𝜈185−205𝜈²+26𝜈⁴sin(𝜈)+18000𝜈²sin²(𝜈)−80cos(𝜈)18480−9980𝜈²+797𝜈⁴+740𝜈sin(𝜈). (4.10)

Now we require that the update of the method (4.1) is phase-fitted and amplification-fitted. By Theorem 3.3 (ii), we have 𝑏₁(𝜈)+𝑏₂(𝜈)1cos2𝜈+𝑏₃(𝜈)1cos2𝜈+𝑏₄(𝜈)(𝜈)=(𝜈)cossin𝜈,𝑏₂(𝜈)1sin2𝜈+𝑏₃(𝜈)1sin2𝜈+𝑏₄(𝜈)(𝜈)=(𝜈)sin1−cos𝜈. (4.11) Solving (4.2) and (4.11), we obtain 𝑏₁(𝜈)=4(𝜈−2sin(𝜈/2))sin(𝜈/2)𝜈²−4+𝜈²(),𝑏+4cos𝜈/2₂(𝜈)=𝐿𝜈⁴−4+𝜈²(),𝑏+4cos𝜈/2₃(𝜈)8(()(()=−𝜈cos𝜈/2−2sin𝜈/2))sin𝜈/2𝜈⁴,𝑏₄(𝜈)=𝑏₁(𝜈), (4.12) where 𝜈𝐿=2sin28𝜈−4𝜈³+𝜈⁵+4𝜈−4+𝜈²𝜈cos2(𝜈)𝜈+8𝜈cos+32sin2−8𝜈²𝜈sin2−16sin(𝜈)+4𝜈²sin(𝜈). (4.13) It can be verified that the method given by (4.1) and (4.12) is of order four. We denote the method by FRK4.

The error coefficients of FRK4 are given by 𝑒𝑡₅₁=𝑊8𝜈²𝑄,𝑒𝑡₅₂=1+5sin(𝜈/2)𝜈−4+𝜈²cos(𝜈/2)+2𝜈−2𝜈³+𝜈cos(𝜈)+4+3𝜈²sin(𝜈/2)−2sin(𝜈)/𝜈⁴𝑄,𝑒𝑡₅₃()()𝜈=1−15sin𝜈/2−𝜈cos𝜈/2𝑄+2³−4+𝜈²sin(𝜈/2)+2sin(𝜈)2𝜈⁴𝑄,𝑒𝑡₅₄((()=1−30𝜈−2sin𝜈/2))sin𝜈/2𝜈²𝑄,𝑒𝑡₅₅=1+10sin(𝜈/2)𝜈−4+𝜈²cos(𝜈/2)+2𝜈−𝜈³+𝜈cos(𝜈)+4+𝜈²sin(𝜈/2)−2sin(𝜈)/𝜈⁴𝑄,𝑒𝑡₅₆()=1−40cos𝜈/4sin³(𝜈/4)𝜈³+4𝜈cos(𝜈/2)−8sin(𝜈/2)𝜈⁴𝑄,𝑒𝑡₅₇((()=1−20𝜈−2sin𝜈/2))sin𝜈/2𝜈²𝑄,𝑒𝑡₅₈((()=1−30𝜈−2sin𝜈/2))sin𝜈/2𝜈²𝑄,𝑒𝑡₅₉=1, (4.14) where 𝑊=160−32𝜈²+8𝜈⁴+2−5+16𝜈²𝜈cos2−160cos(𝜈)+10cos3𝜈2𝜈−140𝜈sin2−5𝜈³𝜈sin2,𝑄=−4+𝜈²𝜈+4cos2. (4.15) Then for FRK4, we have 𝐸𝐶₅(𝜈)=√2881−48697288√𝜈2881²+12901376871115154432√𝜈2881⁴−13289682175603385531190231040√𝜈2881⁶+𝒪𝜈⁸. (4.16)

It can be seen that as 𝜈→0, both the methods Simos4 and FRK4 reduce to a classical RK method of order four (on page 138 of [3]), which we denote by RK4. The method RK4 is called the prototype method or the limit method of the methods Simos4 and FRK4. Moreover, as 𝜈→0, both Simos4 and FRK4 have the same error constant as that of RK4: 𝐸𝐶₅=√2881/48.

4.2. Fifth-Order Methods

Consider the following modified RK type method with FSAL property, the prototype of which can be found on page 167 of [2] or on page 178 of [3], 01/51/53/103/409/404/544/45−56/1532/98/9119372/6561−25360/218764448/6561−212/72919017/3168−355/3346732/524749/176−5103/18656𝑏₁𝑏₂𝑏₃𝑏₄𝑏₅𝑏₆0𝑏₁𝑏₂𝑏₃𝑏₄𝑏₅𝑏₆0. (4.17)

For the method (4.17), the phase-fitting and amplification-fitting conditions (3.18) become 11−5𝑏₂(𝜈)+3𝑏10³(𝜈)+45𝑏₄(𝜈)+89𝑏₅(𝜈)+𝑏₆(𝜈)𝜈²+4𝑏25⁴(𝜈)+424𝑏1215⁵(𝜈)+21𝑏55⁶(𝜈)𝜈⁴−7𝑏550⁶(𝜈)𝜈⁶(𝜈),𝑏=cos₁(𝜈)+𝑏₂(𝜈)+𝑏₃(𝜈)+𝑏₄(𝜈)+𝑏₅(𝜈)+𝑏₆(𝜈)𝜈−9𝑏200³(𝜈)+8𝑏25⁴(𝜈)+32𝑏81⁵(𝜈)+12𝑏₆(𝜈)𝜈³−848𝑏18225⁵(𝜈)+14𝑏275⁶(𝜈)𝜈⁵=sin(𝜈). (4.18) Combining these two equations and the following sufficient conditions in Theorem 2.1 corresponding the trees of orders one, two, and three, we have 𝑏^𝑇𝑒=𝑏₁(𝜈)+𝑏₂(𝜈)+𝑏₃(𝜈)+𝑏₄(𝜈)+𝑏₅(𝜈)+𝑏₆(𝜈)𝑏=1,_𝑇1𝑐=5𝑏₂(𝜈)+3𝑏10³(𝜈)+45𝑏₄(𝜈)+89𝑏₅(𝜈)+𝑏₆(𝜈)=12,𝑏^𝑇𝑐²=1𝑏25²(𝜈)+9𝑏100³(𝜈)+16𝑏25⁴(𝜈)+64𝑏81⁵(𝜈)+𝑏₆(𝜈)=13,𝑏^𝑇9𝐴𝑐=𝑏200³(𝜈)+8𝑏25⁴(𝜈)+32𝑏81⁵(𝜈)+12𝑏₆(𝜈)=16. (4.19) Solving (4.18) and (4.19), we obtain 𝑏₁(𝜈)=28𝜈⁷−235𝜈⁵+28800sin(𝜈)−36600𝜈+7350𝜈³+7800𝜈cos(𝜈)+1350𝜈²sin(𝜈)2884+𝜈²𝜈⁵,𝑏₂(𝜈)𝑏=0,₃(𝜈)=43550𝜈⁵+371𝜈⁷−186750sin(𝜈)+236400𝜈−46500𝜈³−49650𝜈cos(𝜈)−9450𝜈²sin(𝜈)/33394+𝜈²𝜈⁵,𝑏₄(𝜈)=225𝜈⁵+22𝜈⁷+9000sin(𝜈)−10200𝜈+750𝜈³+1200𝜈cos(𝜈)+1350𝜈²sin(𝜈)484+𝜈²𝜈⁵,𝑏₅(𝜈)(𝜈)=−2431800𝜈−1200sin−650𝜈³−600𝜈cos(𝜈)+69𝜈⁵+150𝜈²sin(𝜈)16964+𝜈²𝜈⁵,𝑏₆(𝜈)=11600𝜈−450sin(𝜈)−150𝜈³−150𝜈cos(𝜈)+11𝜈⁵214+𝜈²𝜈⁵. (4.20) As 𝜈→0, we obtain the following expansions: 𝑏₁(𝜈)=35−384149𝜈48384²+899𝜈1161216⁴−2111𝜈10948608⁶+52379𝜈1086898176⁸𝑏+⋯,₂(𝜈)𝑏=0,₃(𝜈)=500+1113977𝜈140238²−1433𝜈841428⁴+94067𝜈222136992⁶−333413𝜈3150306432⁸𝑏+⋯,₄(𝜈)=125−19261𝜈8064²+151𝜈193536⁴−2281𝜈12773376⁶+8059𝜈181149696⁸𝑏+⋯,₅(𝜈)=−2187−6784243𝜈94976²+1791𝜈759808⁴−10215𝜈16715776⁶+12087𝜈79020032⁸𝑏+⋯,₆(𝜈)=11+8411𝜈1764²−187𝜈84672⁴+71𝜈127008⁶−1385𝜈9906624⁸+⋯. (4.21)

By simple computation, we have 𝑏^𝑇𝑐³=𝜈4800−1000𝜈²+200𝜈⁴+37𝜈⁶−1050𝜈cos(𝜈)−15025+𝜈²sin(𝜈)150𝜈⁵4+𝜈²=14+𝒪𝜈²,𝑏^𝑇()=1𝑐⋅𝐴𝑐2𝑏^𝑇𝑐³=18+𝒪𝜈²,𝑏^𝑇𝐴𝑐²=13𝑏^𝑇𝑐³=1𝜈12+𝒪²,𝑏^𝑇𝐴²𝑐=−600+750𝜈²−75𝜈⁴+11𝜈⁶+600cos(𝜈)−450𝜈sin(𝜈)150𝜈⁴4+𝜈²=1𝜈24+𝒪²,𝑏^𝑇𝑐⁴=𝜈925200−195500𝜈²+20850𝜈⁴+2583𝜈⁶−204450𝜈cos(𝜈)−1504805+179𝜈²sin(𝜈)13500𝜈⁵4+𝜈²=15(𝜈),𝑏+𝒪_𝑇𝑐²=1⋅𝐴𝑐2𝑏^𝑇𝑐⁴=1(𝜈)10+𝒪,𝑏^𝑇𝑐⋅𝐴𝑐²=13𝑏^𝑇𝑐⁴=1(𝜈),𝑏15+𝒪_𝑇𝑐⋅𝐴²𝑐2=−−𝜈7200+500𝜈²−150𝜈⁴+33𝜈⁶+75𝜈cos(𝜈)+7595+23𝜈²sin(𝜈)1125𝜈⁵4+𝜈²=1(𝜈),𝑏30+𝒪_𝑇()=1𝐴𝑐⋅𝐴𝑐4𝑏^𝑇𝑐⁴=1(𝜈),𝑏20+𝒪_𝑇𝐴𝑐³=𝜈3000−1150𝜈²+155𝜈⁴+7𝜈⁶−1050𝜈cos(𝜈)+150−13+2𝜈²sin(𝜈)225𝜈⁵4+𝜈²=1(𝜈),𝑏20+𝒪_𝑇𝐴()=𝜈𝑐⋅𝐴𝑐7200−2000𝜈²+300𝜈⁴+33𝜈⁶−1950𝜈cos(𝜈)+150−35+𝜈²sin(𝜈)1500𝜈⁵4+𝜈²=1(𝜈),𝑏40+𝒪_𝑇𝐴²𝑐²=23𝑏^𝑇𝐴()=1𝑐⋅𝐴𝑐(𝜈),𝑏60+𝒪_𝑇𝐴³𝑐=−6𝜈+𝜈³+6sin(𝜈)6𝜈⁵=1(𝜈).120+𝒪 (4.22)

By Theorem 2.1, the coefficients given by (4.17) and (4.20) satisfy all the conditions of order five. But it cannot satisfy all the conditions of order six. For example, 𝑏^𝑇𝑐⁵=899+540022883𝜈9720000²+𝒪𝜈⁴≠16+𝒪(𝜈). (4.23) Therefore, this method is of order five and we denote it as FRK5a.

Corresponding to the twenty sixth-order rooted trees 𝑡_6𝑗,𝑗=1,…,20, the error coefficients of FRK5a are given by 𝑒𝑡₆₁=𝜈−124203600+26615500𝜈²−1403850𝜈⁴+18711𝜈⁶[phone]𝜈cos(𝜈)[phone]+22177𝜈²sin(𝜈)/202500𝜈⁵4+𝜈²,𝑒𝑡₆₂=𝑒𝑡₆₃=𝑒𝑡₆₅=𝑒𝑡_6,10=𝑒𝑡₆₁,𝑒𝑡₆₄=−5493600𝜈+551000𝜈³+31800𝜈⁵−3879𝜈⁷+756600𝜈cos(𝜈)+6007895+1028𝜈²sin(𝜈)/5625𝜈⁵4+𝜈²,𝑒𝑡₆₆=𝜈−6667200+2224000𝜈²−194700𝜈⁴+6957𝜈⁶+2084700𝜈cos(𝜈)−300−15275+1393𝜈²sin(𝜈)/16875𝜈⁵4+𝜈²,𝑒𝑡₆₇=𝜈−2296800+563000𝜈²−36600𝜈⁴+873𝜈⁶+565800𝜈cos(𝜈)+6002885+14𝜈²sin(𝜈)/5625𝜈⁵4+𝜈²,𝑒𝑡₆₈=𝑒𝑡₆₇,𝑒𝑡₆₉=𝜈57600−1600𝜈²−948𝜈⁴+75𝜈⁶−4800𝜈cos(𝜈)−480011+2𝜈²sin(𝜈)/75𝜈⁵4+𝜈²𝑒𝑡_6,11=𝑒𝑡₆₄,𝑒𝑡_6,12=𝜈−652800+333000𝜈²−36400𝜈⁴+1533𝜈⁶+290550𝜈cos(𝜈)−450−805+283𝜈²sin(𝜈)/2250𝜈⁵4+𝜈²,𝑒𝑡_6,13=𝜈−602400+300000𝜈²−32500𝜈⁴+1359𝜈⁶+262650𝜈cos(𝜈)−450−755+249𝜈²sin(𝜈)/2250𝜈⁵4+𝜈²,𝑒𝑡_6,14=𝑒𝑡_6,16=𝑒𝑡_6,13,𝑒𝑡_6,15=𝜈5280−360𝜈²−56𝜈⁴+5𝜈⁶−600𝜈cos(𝜈)−36013+2𝜈²sin(𝜈)5𝜈⁵4+𝜈²,𝑒𝑡_6,17=𝜈−352800+101000𝜈²−7800𝜈⁴+243𝜈⁶+97800𝜈cos(𝜈)−600−425+16𝜈²sin(𝜈)/375𝜈⁵4+𝜈²,𝑒𝑡_6,18=𝜈−3360+1320𝜈²−128𝜈⁴+5𝜈⁶+1200𝜈cos(𝜈)−360−6+𝜈²sin(𝜈)5𝜈⁵4+𝜈²,𝑒𝑡_6,19=𝑒𝑡_6,18,𝑒𝑡_6,20=1+24−600𝜈+150𝜈³−11𝜈⁵+150𝜈cos(𝜈)+450sin(𝜈)5𝜈⁵4+𝜈² (4.24) Then for FRK5a, we have 𝐸𝐶₆(𝜈)=√𝑁202500𝜈⁵4+𝜈², (4.25) where 𝑁=200598105893906250+704003855089083750𝜈²−297809752687323750𝜈⁴+56298277495370000𝜈⁶−5785233488749200𝜈⁸[phone]𝜈¹⁰−15770219564700𝜈¹²[phone]𝜈¹⁴+300𝜈²−1028483389544400+258042555395500𝜈²−17441989606650𝜈⁴[phone]𝜈⁶cos(𝜈)−1125017830942746125−509988761071𝜈²[phone]𝜈⁴()(𝜈)cos2𝜈−1033593801921000000𝜈sin+195145999800990000𝜈³sin(𝜈)−14464008183675000𝜈⁵sin(𝜈)+1060750694812500𝜈⁷sin(𝜈)−58620919851900𝜈⁹sin(𝜈)+115600689172687500𝜈sin(2𝜈)−103692151057500𝜈³sin(2𝜈). (4.26)

Now we require that the update of the method (4.17) is phase-fitted and amplification-fitted. By Theorem 3.3 (ii), we have 𝑏₁(𝜈)+𝑏₂(𝜈)1cos5𝜈+𝑏₃(𝜈)3cos𝜈10+𝑏₄(𝜈)4cos5𝜈8+cos9𝜈(𝜈)=(𝜈)+cossin𝜈,𝑏₂(𝜈)1sin5𝜈+𝑏₃(𝜈)3sin𝜈10+𝑏₄(𝜈)4sin5𝜈8+sin9𝜈(𝜈)=(𝜈)+sin1−cos𝜈. (4.27) We can solve the algebraic systems (4.18), (4.27), and the last two equations in (4.19) for 𝑏_𝑖(𝜈). The closed expressions for 𝑏_𝑖(𝜈) is too complicated. As 𝜈→0, they have the series expressions as follows: 𝑏₁(𝜈)=35−38420437583𝜈11904122880²−119367195208354081𝜈968707998823219200000⁴+95011884964518844932145979𝜈2252269346255956397260800000000⁶−2870157856674050512150034054708268901𝜈403216681907224475027036857958400000000000⁸𝑏+⋯,₂(𝑣)𝑏=0,₃(𝜈)=500+111333964523𝜈8625839040²−7661109532360639𝜈140386979516958720000⁴−22790810005346940380832667𝜈326403096664437431009280000000⁶+172880608491993732432855997019686657𝜈11686983514654709393361771429888000000000⁸𝑏+⋯,₄(𝜈)=125−19210830931𝜈1984020480²[phone]𝜈32290266627440640000⁴+1736990439801521504752259𝜈75075644875198546575360000000⁶−11178170615372984995156901626580537𝜈2688111212714829833513579053056000000000⁸𝑏+⋯,₅(𝜈)=−2187+67841224279𝜈2596372480²+18753211139151121𝜈23475776834764800000⁴−1958127900513001510763993𝜈18193914853157398118400000000⁶+126793156608794251730624145537232661𝜈9771597686823455499054848409600000000000⁸𝑏+⋯,₆(𝜈)=11+842401663𝜈868008960²−44892540738708751𝜈70634958247526400000⁴+12565766033461779855286229𝜈164227973164496820633600000000⁶−23848033853016506130445461060833921𝜈2672837853551677391277706444800000000000⁸+⋯. (4.28)

It can be verified that the method given by (4.17) and (4.28) is of order five. We denote the method by FRK5b.

The error coefficients of FRK5b are given by 𝑒𝑡₆₁=1−90015684728327𝜈2092521600000²+7668220591385576111𝜈6811228116725760000000⁴−3832599437422925637715480451𝜈26393781401436989030400000000000⁶+37040225493178191951169407554740296497𝜈2025083781900122921452752076800000000000000⁸𝑒𝑡+⋯,₆₂=𝑒𝑡₆₃=𝑒𝑡₆₅=𝑒𝑡_6,10=𝑒𝑡₆₁,𝑒𝑡₆₄=−1312943077𝜈58125600000²+101770174398983783𝜈135143415014400000000⁴−50361593135536115530775449𝜈733160594484360806400000000000⁶+722243608228427223857425462329823121𝜈78753258185004780278718136320000000000000⁸𝑒𝑡+⋯,₆₆=−55576679𝜈174376800000²+3349151183983384651𝜈2838011715302400000000⁴−10489794742946357634833857𝜈62842336670088069120000000000⁶+24830088427729996581426427673789858999𝜈1181298872775071704180772044800000000000000⁸𝑒𝑡+⋯,₆₇=−220019101𝜈58125600000²+854131951646627753𝜈946003905100800000000⁴−92480070367754551230905137𝜈733160594484360806400000000000⁶+1287242651211484482736127020303275113𝜈78753258185004780278718136320000000000000⁸𝑒𝑡+⋯,₆₈=𝑒𝑡₆₇,𝑒𝑡₆₉=1+252237059𝜈96876000²[phone]𝜈7883365875840000000⁴−489228700990972151904469𝜈6109671620703006720000000000⁶+38954008540439948287239141526960033𝜈3281385757708532511613255680000000000000⁸𝑒𝑡+⋯,_6,11=𝑒𝑡₆₄,𝑒𝑡_6,121=−[phone]𝜈46500480000²+1020783650901433441𝜈756803124080640000000⁴−349391329566451655705258747𝜈1759585426762465935360000000000⁶+223430128326558463436293654635080087𝜈9000372364000546317567787008000000000000⁸𝑒𝑡+⋯,_6,131=−[phone]𝜈46500480000²+930238779819453217𝜈756803124080640000000⁴−322126891973817461296228667𝜈1759585426762465935360000000000⁶+7281345860082294443565287484305877349𝜈315013032740019121114872545280000000000000⁸𝑒𝑡+⋯,_6,14=𝑒𝑡_6,16=𝑒𝑡_6,13,𝑒𝑡_6,15=7147843𝜈258336000²−5833795112598163𝜈21022309002240000000⁴−2118373575562062120724543𝜈48877372965624053760000000000⁶+69165957915286832600826502459212257𝜈8750362020556086697635348480000000000000⁸𝑒𝑡+⋯,_6,17=−1495267𝜈298080000²+1575285628973489𝜈693043153920000000⁴−221743464671062174816067𝜈751959584086523904000000000⁶+71862921159172247356437619228883827𝜈2019314312436020007146618880000000000000⁸𝑒𝑡+⋯,_6,18=844759𝜈258336000²+6871050469020383𝜈3003187000320000000⁴−15338466371260326350460259𝜈48877372965624053760000000000⁶+332683402120349742930981514490030141𝜈8750362020556086697635348480000000000000⁸𝑒𝑡+⋯,_6,19=𝑒𝑡_6,18,𝑒𝑡_6,201=−5−218333𝜈8611200²+4081140067155341𝜈700743633408000000⁴−1142342366678343623207839𝜈1629245765520801792000000000⁶+23848033853016506130445461060833921𝜈291678734018536223254511616000000000000⁸+⋯. (4.29) Then for FRK5b, we have 𝐸𝐶₆(𝜈)=√33801+9002207344629853418504320000√𝜈33801²+1035095647490481214191033791151131607867237068800000000√𝜈33801⁴−696604907215699160267620098314707264612190465287610822576869356860211200000000000√𝜈33801⁶+24352927165428290711932956191142918847286589328298244559547431372665305140913187085379096747809177600000000000000√𝜈33801⁸+⋯. (4.30)

It can be seen that as 𝜈→0, both the methods FRK5a and FRK5b reduce to a classical RK method of order five (on page 178 of [3]), which we denote by RK5. The method RK5 is called the prototype method or the limit method of the methods FRK5a and FRK5b. Moreover, as 𝜈→0, both FRK5a and FRK5b have the same error constant as that of RK5: 𝐸𝐶₆=√33801/900.

5. Analysis of Stability and Phase Properties

In order to investigate the linear stability of the methods (2.1) or any other frequency-depending method, we applied it to the linear test equation as 𝑦^′=i𝜆𝑦,𝜆>0. (5.1) One step of computation yields 𝑦₁=𝑀i𝜇,i𝜈𝑦₀,𝜇=ℎ𝜆, (5.2) where 𝑀i𝜇,i𝜈=1+i𝜇𝑏^𝑇𝐼−i𝜇𝐴⁻¹𝑒=1−𝜇²𝑏^𝑇(𝜈)𝐴𝑒+𝜇⁴𝑏^𝑇(𝜈)𝐴³𝑒−𝜇⁶𝑏^𝑇(𝜈)𝐴⁵𝑒+𝜇⁸𝑏^𝑇(𝜈)𝐴⁷𝑒𝜇𝑏−⋯+i^𝑇(𝜈)𝑒−𝜇³𝑏^𝑇(𝜈)𝐴²𝑒+𝜇⁵𝑏^𝑇(𝜈)𝐴⁴𝑒−𝜇⁷𝑏^𝑇(𝜈)𝐴⁶𝑒+⋯ (5.3) is called the stability matrix of the method (2.1).

Definition 5.1. The region in the 𝜇-𝜈 plane Ω={(𝜇,𝜈)∣𝜇>0,𝜈>0,|𝑅(i𝜇,i𝜈)∶≤1} is called the imaginary stability region of the method (2.1) and the closed curve defined by |𝑅(i𝜇,i𝜈)|=1 is called the stability boundary of the method.

The imaginary stability regions of the methods Simos4, FRK4, FRK5a and FRK5b are depicted in Figures 1 and 2.

(a)

(b)

Figure 1 

The imaginary stability region of Simos4 (a) and FRK4a (b).

(a)

(b)

Figure 2 

The imaginary stability region of FRK5a (a) and FRK5b (b).

Definition 5.2. The quantities 𝜙𝜇,𝜈=𝜇−arg𝑀i𝜇,i𝜈,𝑑𝜇,𝜈=1−𝑀i𝜇,i𝜈 (5.4) are called the phase lag and the error of amplification factor of the method (2.1), respectively. If 𝜙(𝜇,𝜈)=𝒪(𝜇^𝑞+1), and 𝑑(𝜇,𝜈)=𝒪(𝜇^𝑝+1), then the method (2.1) is said to be dispersive of order 𝑞 and dissipative of order 𝑝, respectively (see Ozawa [24] and Van de Vyver [25]).

By definition, an FRK method has zero phase lag and zero dissipation when applied to the standard linear oscillator (5.1) whose frequency 𝜆 coincides with the fitting frequency 𝜔. That is to say, 𝜙(𝜇,𝜇)=0,𝑑(𝜇,𝜇)=0.

Suppose that the internal stages have been exact for the linear equation (5.1), that is, 𝑌_𝑖=exp(i𝑐_𝑖𝜇)𝑦₀, then the update gives 𝑦₁=𝑀_𝑢i𝜇,i𝜈𝑦₀, (5.5) where 𝑀_𝑢i𝜇,i𝜈=1+i𝜇^𝑠_𝑖=1𝑏_𝑖(𝜈)exp𝑐_𝑖𝜇=1+𝜇^𝑠_𝑖=1𝑏_𝑖(𝜈)sin𝑐_𝑖𝜇+i𝜇^𝑠_𝑖=1𝑏_𝑖(𝜈)cos𝑐_𝑖𝜇. (5.6) Accordingly, the quantities 𝜙_𝑢𝜇,𝜈=𝜇−arg𝑀_𝑢i𝜇,i𝜈,𝑑_𝑢𝜇,𝜈=1−𝑀_𝑢i𝜇,i𝜈 (5.7) are called the phase lag and the dissipation of the update of the method (2.1), respectively. If 𝜙_𝑢(𝜇,𝜈)=𝒪(𝜇^𝑞+1), and 𝑑_𝑢(𝜇,𝜈)=𝒪(𝜇^𝑝+1), then the update of the method (2.1) is said to be dispersive of order 𝑞 and dissipative of order 𝑝, respectively.

Letting 𝜈=𝑟𝜇, we obtain the phase lags and dissipations of the four FRK methods derived in Section 4.(i) Simos4: 𝜙𝜇,𝜈=1−𝑟²𝜇120⁵+𝒪𝜇⁷,𝑑𝜇,𝜈=1−𝑟²𝜇144⁶+𝒪𝜇⁸,𝜙_𝑢𝜇,𝜈=−1+4𝑟²𝜇2880⁵+𝒪𝜇⁷,𝑑_𝑢𝜇,𝜈=−1+4𝑟²𝜇5760⁶+𝒪𝜇⁸,𝜙_𝑢1𝜇,𝜇=−𝜇576⁵+𝒪𝜇⁷,𝑑_𝑢1𝜇,𝜇=−𝜇1152⁶+𝒪𝜇⁸. (5.8)(ii) FRK4: 𝜙𝜇,𝜈=12−17𝑟²+5𝑟⁴𝜇1440⁵+𝒪𝜇⁷,𝑑𝜇,𝜈=4−5𝑟²+𝑟⁴𝜇576⁶+𝒪𝜇⁸,𝜙_𝑢𝜇,𝜈=−1−9𝑟²+10𝑟⁴𝜇2880⁵+𝒪𝜇⁷,𝑑_𝑢𝜇,𝜈=−1−9𝑟²+10𝑟⁴𝜇5760⁶+𝒪𝜇⁸,𝜙_𝑢𝜇,𝜇=0,𝑑_𝑢𝜇,𝜇=0. (5.9)(iii) FRK5a: 𝜙=𝑟𝜇,𝜈²−1𝜇2100⁷+𝒪𝜇⁹,𝑑𝜇,𝜈=1−𝑟²𝜇3600⁶+𝒪𝜇⁸,𝜙_𝑢𝜇,𝜈=4021+144855𝑟²𝜇2449440000⁷+𝒪𝜇⁹,𝑑_𝑢𝜇,𝜈=−7+570𝑟²𝜇4536000⁶+𝒪𝜇⁸,𝜙_𝑢𝜇,𝜇=5317𝜇87480000⁷+𝒪𝜇⁹,𝑑_𝑢𝜇,𝜇=−577𝜇4536000⁶+𝒪𝜇⁸. (5.10)(iv) FRK5b: 𝜙𝜇,𝜈=−20666880+22148840𝑟²−3027619𝑟⁴+1545659𝑟⁶𝜇43400448000⁷+𝒪𝜇⁹,𝑑𝜇,𝜈=1722240−2154757𝑟²+432517𝑟⁴𝜇6200064000⁶+𝒪𝜇⁸,𝜙_𝑢𝜇,𝜈=9618232+190446368𝑟²−408728565𝑟⁴[phone]𝑟⁶𝜇5859060480000⁷+𝒪𝜇⁹,𝑑_𝑢𝜇,𝜈=−9568−422949𝑟²+432517𝑟⁴𝜇6200064000⁶+𝒪𝜇⁸,𝜙_𝑢𝜇,𝜇=0,𝑑_𝑢𝜇,𝜇=0. (5.11)Therefore, when integrating the test equation (5.1), the methods Simos4 and FRK4 are dispersive of order four and dissipative of order five, and FRK5a and FRK5b are dispersive of order six and dissipative of order five. Similar conclusions for their updates can be drawn from the above analysis.

Note that if the update of (2.1) is phase-fitted and amplification-fitted, it must be true that 𝜙_𝑢(𝜇,𝜇)=0 and 𝑑_𝑢(𝜇,𝜇)=0, as is verified above for the two methods FRK4 and FRK5b.

6. Numerical Experiments

In order to examine the effectiveness of the new FRK methods proposed in this paper, we apply them to four test problems. We also employ several highly efficient integrators from the literature for comparison. The numerical methods we choose for experiments are as follow: (i)FRK5a: the seven-stage RK method of order five given by (4.17) and (4.20) in Section 4 of this paper; (ii)FRK5b: the seven-stage RK method of order five given by (4.17) and (4.28) in Section 4 of this paper; (iii)Simos4: the RK method of order four presented in [23], that is, the four-stage FRK method of order four given by (4.1) and (4.12) in Section 4 of this paper; (iv)FRK4: the four-stage RK method of order four given by (4.1) and (4.12) in Section 4 of this paper; (v)ARK4: the second four-stage adapted RK method of order four given in the Subsection 3.2 of Franco [11]; (vi)EFRK4: the four-stage exponentially fitted RK method of order four given in [26]; (vii)RK4: the classical RK method of order four presented in [3] (the prototype method of Simos4 and FRK4); (viii)RK5: the classical RK method of order five presented in [3] (the prototype method of FRK5a and FRK5b). In the figures showing the efficiency of these integrators, the horizontal axis stands for the number of function evaluations required and the vertical axis stands for the digital logarithm of the maximal global error (MGE).

Problem 1. Consider the following orbit problem studied in [27]: 𝑝^′=−𝑞+0.001cos(𝑡),𝑞^′𝑝(0)(0)=𝑝,=0,𝑞=1. (6.1) This equation is the real part of the “almost” periodic orbit problem as 𝑦^′′+𝑦=0.001𝑒^𝑖𝑡,𝑦(0)=1,𝑦^′(0)=0.9995𝑖,𝑦∈ℂ. (6.2) The exact solution to this problem is 𝑞(𝑡)=cos(𝑡)+0.0005𝑡sin(𝑡). We select the fitting frequency 𝜔=1.0007 and integrate the equation on the interval [0,1000] with step sizes ℎ=1/2^𝑖,𝑖=0,1,2,3. The numerical results are presented in Figure 3(a).

(a)

(b)

Figure 3 

Efficiency curves of Problem 1 (a) and Problem 2 (b).

Problem 2. Consider the following linear problem studied in [9]: 𝑦^′′+𝜔²𝑦=𝜔²(𝑡),𝑦(0)−1sin=1,𝑦^′(0)=𝜔+1. (6.3) The exact solution to this problem is 𝑦(𝑡)=cos(𝜔𝑡)+sin(𝜔𝑡)+sin(𝑡). In this experiment, we take the parameter 𝜔=20 and integrate the equation on the interval [0,100] with the step sizes ℎ=1/(8𝑖),𝑖=2,3,4,5. The numerical results are presented in Figure 3(b).

Problem 3. Consider the prey-predator system in ecology (see [28]) as ̇,̇,𝑢=𝑢𝛼−𝛽𝑣𝑣=𝑣𝛾𝑢−𝛿 (6.4) where 𝑢(𝑡) is the number of prey at time 𝑡, 𝑣(𝑡) is the number of predators at time 𝑡, 𝑢̇ and 𝑣̇ are first derivatives with respect to time, and 𝛼,𝛽,𝛾,𝛿 are positive constants. For given initial values 𝑢(0)=𝑢₀,𝑣(0)=𝑣₀, the system has unique solution, but analytical solution is not available. However, the system has an invariant 𝐼(𝑢,𝑣)=𝛿ln𝑢+𝛼ln𝑢−𝛾𝑢−𝛽𝑣 in the sense that (d/d𝑡)𝐼(𝑢(𝑡),𝑣(𝑡))≡0. The efficient of the methods will be tested by measuring the error growth of the invariant 𝐼(𝑢(𝑡),𝑣(𝑡)) against the number of function evaluations required.
In this experiment, we take the values of parameters 𝛼=2,𝛽=𝛾=𝛿=1 and take initial data 𝑢₀=1.6,𝑣₀=2.2. We select the fitting frequency 𝜔=1.0075 and integrate the equation on the interval [0,30] with step sizes ℎ=1/2^𝑖,𝑖=2,3,4,5. The numerical results are presented in Figure 4(a).

(a)

(b)

Figure 4 

Efficiency curves of Problem 3 (a) and Problem 4 (b).

Problem 4. We consider the following two-gene regulatory system without self-regulation (see Widder et al. [29] and Polynikis et al. [30]): 𝑟̇₁=𝑚₁𝐻⁺𝑝₂;𝜃₂,𝑛₂−𝛾₁𝑟₁,𝑟̇₂=𝑚₂𝐻⁻𝑝₁;𝜃₁,𝑛₁−𝛾₂𝑟₂,𝑝̇₁=𝑘₁𝑟₁−𝛿₁𝑝₁,𝑝̇₂=𝑘₂𝑟₂−𝛿₂𝑝₂, (6.5) where 𝑖=1,2, 𝑟_𝑖 is the concentration of mRNA 𝑅_𝑖 produced by gene 𝐺_𝑖, 𝑝_𝑖 is the concentration of protein 𝑃_𝑖, 𝑚_𝑖 is the maximal transcription rate of 𝐺_𝑖, 𝑘_𝑖 is the translation rate of 𝑅_𝑖, 𝛾_𝑖 is the degradation rate of 𝑅_𝑖, and 𝛿_𝑖 is the degradation rate of 𝑃_𝑖. The two functions 𝐻⁺𝑝₂;𝜃₂,𝑛₂=𝑝^𝑛2₂𝑝^𝑛2₂+𝜃^𝑛2₂,𝐻⁻𝑝₁;𝜃₁,𝑛₁=𝜃^𝑛1₁𝑝^𝑛1₁+𝜃^𝑛1₁ (6.6) are the Hill functions of activation and repression, respectively. The parameters 𝑛₁,𝑛₂ are the Hill coefficients, 𝜃₁,𝜃₂ are the expression thresholds.
The solution (𝑟^∗₁,𝑟^∗₂,𝑝^∗₁,𝑝^∗₂) of the system 𝑟̇₁=𝑟̇₂=𝑝̇₁=𝑝̇₂=0 is an equilibrium of the system (6.5). Now we take the values of parameters as follows: 𝑛₁=𝑛₂=3,𝑚₁=1.15,𝑚₂=2.35,𝑘₁=𝑘₂=1,𝛾₁=𝛾₂𝛿=1,₁=𝛿₂=1,𝜃₁=𝜃₂=0.21. (6.7) For any initial point (𝑟₁(0),𝑟₂(0),𝑝₁(0),𝑝₂(0)) near the equilibrium, the system (6.5) has a stable limit cycle.
In this experiment, we take initial data (0.6,0.8,0.4,0.6) and select the fitting frequency 𝜔=1.48. The problem is integrated on the interval [0,20] with step sizes ℎ=1/2^𝑖,𝑖=1,2,3,4. The numerical results are presented in Figure 4(b).
It can be seen from Figures 3 and 4 that the FRK methods are more efficient than their prototype RK methods and are more efficient than the other frequency depending methods of the same algebraic order.

7. Conclusions

In this paper, classical Runge-Kutta methods are adapted to the time integration of initial value problems of first order differential equations whose solutions have oscillatory properties. The newly developed phase-fitted and amplification-fitted Runge-Kutta methods (FRK) adopt functions of the product 𝜈=𝜔ℎ of the fitting frequency 𝜔 and the step size ℎ as weight coefficients in the update. FRK methods have zero dispersion error and zero dissipation error when applied to the standard linear oscillator 𝑦^′=i𝜔𝑦. That is to say, they preserve initial phase and amplification with time. Therefore, FRK integrators are a kind of integrators which preserve the oscillation structure of the problem.

As the fitting frequency tends to zero, FRK methods reduce to their classical prototypes methods. Furthermore, an FRK method has the same algebraic order and the same error constant with its prototypes method. Numerical experiments illustrate the high efficiency of FRK methods compared with their prototype methods and some other frequency depending methods like exponentially fitted RK type methods.

Theorem 3.3 gives a pair of sufficient conditions for modified RK type methods to be phase-fitted and amplification-fitted. The coefficients of the methods Simos4 and FRK5a are obtained with these conditions combined with an appropriate number of order conditions associated to low order trees. Since the updates of the methods FRK4 and FRK5b are also phase-fitted and amplification-fitted, they may be more efficient than those methods whose updates do not have this property.

In practical computations of oscillatory problems, the true frequency is, in general, not available. The fitting frequency contained in an FRK method is just an estimate of the true frequency. Sometimes the choice of the value of the fitting frequency 𝜔 for the FRK method affects their effectiveness to some extent. For instance, in the experiment of Problem 1, we find that the value of 𝜔=1.0007 is superior to the usual choice 𝜔=1. For approaches to estimating principal frequencies we refer to the papers [31–36].

This paper provides a convenient approach to constructing FRK methods. Other effective approaches are possible. For example, the fitting frequency can also be incorporated into internal coefficients 𝑎_𝑖𝑗 or into nodes 𝑐_𝑖 just as in [12]. Or frequency-dependent coefficients may be introduced to the terms 𝑦_𝑛 in the modified (2.1) as in [26]. In these cases, the internal stages can also be made phase-fitted and amplification-fitted.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 11171155), the Fundamental Research Fund for the Central Universities (Grants no. Y0201100265 and KYZ201125), and the Research Fund for the Doctoral Program of Higher Education (Grant no. 20100091110033). The authors would like to thank the referees for their valuable comments and suggestions.