A Variable Step-Size Exponentially Fitted Explicit Hybrid Method for Solving Oscillatory Problems

Abstract

An exponentially fitted explicit hybrid method for solving oscillatory problems is obtained. This method has four stages. The first three stages of the method integrate exactly differential systems whose solutions can be expressed as linear combinations of {1,𝑥,exp(𝜇𝑥),exp(−𝜇𝑥)},𝜇∈C, while the last stage of this method integrates exactly systems whose solutions are linear combinations of {1,𝑥,𝑥²,𝑥³,𝑥⁴,exp(𝜇𝑥),exp(−𝜇𝑥)}. This method is implemented in variable step-size code basing on an embedding approach. The stability analysis is given. Numerical experiments that have been carried out show the efficiency of our method.

1. Introduction

There has been great interest in the research of new methods for numerically solving the second-order initial value problems of the form 𝐲^(𝑥)=𝐟(𝑥,𝐲(𝑥)),𝐲𝑥₀=𝐲₀,𝐲^𝑥₀=𝐲^0(1.1) whose solution exhibits specific oscillatory behavior. Such problems arises in celestial mechanics, in quantum mechanical scattering problems, and elsewhere, and they can be solved by using general purpose methods or by using codes specially adapted to the structure or to the solution of the problem. In the case of adapted hybrid methods, particular algorithms have been proposed by several authors (e.g., see [1–3]) to solve these classes of problems.

This paper concerns derivation of a new method for the numerical integration of (1.1) with oscillatory solution through the usage of an exponential-fitting technique proposed by Vanden Berghe and Van Daele [4]. The coefficients of the new method depend on the product of frequency and step size, thus the method can only be applied when a good estimate of the dominant frequency of the solution is known in advance. The exponential-fitting technique is generally described in the context of linear multistep methods as follows.

Consider linear 𝑘-step methods of the form^𝑘_𝑗=1𝛼_𝑗𝑦_𝑛+𝑗=ℎ^2𝑘_𝑗=1𝛽_𝑗𝑦^{𝑛+𝑗},(1.2) where 𝛼_𝑘=1, 𝛼₀, and 𝛽₀ do not both vanish. Associate with the method is the linear operator𝐿[ℎ,𝐚]𝑦(𝑥)=^𝑘_𝑗=0𝛼_𝑗𝑦(𝑥+𝑗ℎ)−ℎ^2𝑘_𝑗=0𝛽_𝑗𝑦^(𝑥+𝑗ℎ),(1.3) where 𝐚=[𝛼₀,𝛼₁,…,𝛼_𝑘,𝛽₀,𝛽₁,…,𝛽_𝑘]. The linear operator 𝐿 is said to integrate exactly the function 𝑦(𝑥) if 𝐿[ℎ,𝐚]𝑦(𝑥)=0. We now impose the requirement that the operator 𝐿 integrates exactly the following functions:1,𝑥,𝑥²,…,𝑥^𝐾,𝑒^±𝜇𝑥,𝑥𝑒^±𝜇𝑥,𝑥²𝑒^±𝜇𝑥,…,𝑥^𝑃𝑒^±𝜇𝑥(1.4) with 𝐾+2𝑃=𝑀−3. Then, we solve the algebraic equations for 𝛼_𝑗 and 𝛽_𝑗, 𝑗=0,1,…,𝑘. In particular, the step-by-step exponential-fitting procedure is described as follows.

Step 1. Find the maximal 𝑀 such that 𝐿 integrates exactly the set of power functions {1,𝑥,𝑥²,…,𝑥^𝑀−1} and that the resulting coefficients have the same value as those of a classical method. The classical method is the method with constant coefficients.

Step 2. Consider this set of functions 1,𝑥,𝑥²,…,𝑥^𝐾,𝑒^±𝜇𝑥,𝑥𝑒^±𝜇𝑥,𝑥²𝑒^±𝜇𝑥,…,𝑥^𝑃𝑒^±𝜇𝑥(1.5) with 𝐾+2𝑃=𝑀−3. Using 𝑀 obtained from Step 1, we find suitable 𝑃 and 𝐾 which correspond to a finite set consisting of power functions and exponentials. Then, solve equations obtained by assuming that 𝐿 integrates exactly elements of this set.

Considering each formula stage of a hybrid method as a linear multistep method with nonequidistant grid points; we derive our exponentially fitted hybrid method. This approach has also been employed by Simos and Vigo-Aguiar [1] and D’Ambrosio et al. [5] in their derivation of exponential-fitting methods.

2. Embedded Pairs of Hybrid Methods

For the numerical solution of (1.1), we consider the explicit hybrid method which has been established by Franco [6]: 𝑌_𝑖𝑌₁=𝑦_𝑛−1,𝑌₂=𝑦_𝑛,=1+𝑐_𝑖𝑦_𝑛−𝑐_𝑖𝑦_𝑛−1+ℎ^2𝑖−1_𝑗=1𝑦𝑎_𝑖𝑗𝑓𝑥_𝑛+𝑐_𝑗ℎ,𝑌_𝑗,𝑖=3,…,𝑠,_𝑛+1=2𝑦_𝑛−𝑦_𝑛−1+ℎ²𝑏¹𝑓_𝑛−1+𝑏₂𝑓_𝑛+^𝑠_𝑖=3𝑏_𝑖𝑓𝑥_𝑛+𝑐_𝑖ℎ,𝑌_𝑖,(2.1) where 𝑓_𝑛−1 and 𝑓_𝑛 represent 𝑓(𝑥_𝑛−1,𝑦_𝑛−1) and 𝑓(𝑥_𝑛,𝑦_𝑛), respectively. The method requires 𝑠−1 function evaluations or stages at each step of integration. The explicit hybrid method can be represented by the Butcher tableau:−01000⋯0𝑐000⋯0₃⋮𝑎₃₁𝑎₃₂0⋯0𝑐⋮⋮⋱⋮⋮_𝑠𝑎_𝑠1𝑎_𝑠2⋯𝑎_{𝑠,𝑠−1}0𝑏₁𝑏₂⋯𝑏_𝑠−1𝑏_𝑠=𝐜𝐀𝐛^𝑇(2.2) Adopting the concept of embedded 𝑝(𝑞) pair of Runge-Kutta-Nystrom (RKN) methods (see [7, 8]), we define an embedded 𝑝(𝑞) pair of hybrid methods to be based on the hybrid method (𝐜,𝐀,𝐛) of order 𝑝 and another hybrid method (𝐜,𝐀,𝐛) of order 𝑞<𝑝, represented by the following tableau: 𝐜𝐀𝐛^𝑇𝐛^𝑇(2.3) Embedded pairs of explicit hybrid methods are used in variable step-size algorithm because they provide cheap error estimation. A local error estimation is determined by the formula LTE=‖‖𝐲_𝑛+1−𝐲_𝑛+1‖‖,(2.4) where 𝐲_𝑛+1 and 𝐲_𝑛+1 are solutions obtained using the higher-order and the lower-order formula, respectively. The LTE is used to control the step-size of which the procedure is given by

(i) if tol/div<LTE<div⋅tol, then ℎ_𝑛+1=ℎ_𝑛,

(ii) if LTE≤tol/div, then ℎ_𝑛+1=2ℎ_𝑛,

(iii) if LTE≥div⋅tol, then ℎ_𝑛+1=(1/2)ℎ_𝑛 and repeat the step,

where, from numerical experiments, div is chosen to be 2¹⁷. We do not allow step-size change after each step because it would contribute to unnecessary rounding errors. If the step is acceptable (i.e., tol/div < LTE < div · tol and LTE ≤ tol/div), then we adopt the widely used of performing local extrapolation: although the LTE is the error estimation for the lower-order formula, the solutions obtained by using the higher-order formula are actually accepted at each point.

3. Construction of Hybrid Methods

3.1. The Sixth-Order Hybrid Method of Franco

In this section, we consider the sixth-order hybrid method which is zero dissipative and dispersive of order six derived in [6]. The interval of periodicity of the method is (0, 2.75). In order for the method to be implemented in variable step-size code as described in Section 2, we derive a three-stage fourth-order explicit hybrid method. Order conditions for the fourth-order hybrid methods with 𝑠=4 as given in [9] are listed as follows.

Order conditions:⁴_𝑖=1𝑏_𝑖=1,₄_𝑖=1𝑏_𝑖𝑐_𝑖=0,₄_𝑖=1𝑏_𝑖𝑐^2𝑖=16,₄_4𝑖=1_𝑗=1𝑏_𝑖𝑎_𝑖𝑗=112₄,_𝑖=1𝑏_𝑖𝑐^3𝑖=0,₄_4𝑖=1_𝑗=1𝑏_𝑖𝑐_𝑖𝑎_𝑖𝑗=112₄,_4𝑖=1_𝑗=1𝑏_𝑖𝑎_𝑖𝑗𝑐_𝑗=0.(3.1) The two-step explicit hybrid methods (2.1) are a subclass of the two-step hybrid methods considered by Coleman [9] by taking 𝑐₁=−1,𝑐₂=0,  𝑎₂₁=0, and 𝑎_𝑖𝑗=0  (𝑗≥𝑖). Substituting 𝑐₁=−1, 𝑐₂=0,  𝑎₂₁=0, and 𝑎_𝑖𝑗=0  (𝑗≥𝑖) together with coefficients 𝑐_𝑖 and 𝑎_𝑖𝑗 of the sixth-order formula into (3.1), we get the following values of coefficients 𝐛 for the three-stage fourth-order method:𝑏₁=568,𝑏₂=4742,𝑏₃=−512,𝑏₄=80357.(3.2) The embedded 6(4) pair of hybrid methods is denoted as EHM6(4) method and represented by the tableau:−01000001000005471251121500010119−12000120070010002−11204−1744−1744145301681412284535507275684472−15238507(3.3)

3.2. Variable Step-Size Exponentially Fitted Explicit Hybrid Method

We consider the four-stage explicit hybrid method given by the following tableau: −010000010000057𝑎₃₁𝑎₃₂000101192000𝑎⁴²−1𝑎₄₃002−11204−1744𝑎⁵³𝑎₅₄0𝑏₁𝑏₂𝑏₃𝑏₄𝑏₅(3.4) Note that the 𝐜-values and some of the 𝐀-values are taken from the sixth-order hybrid method of Franco as described briefly in Section 3.1. The formula of a four-stage explicit hybrid method is𝑌𝑌₁=𝑦_𝑛−1,𝑌₂=𝑦_𝑛,(3.5)₃𝑌=1+𝑐₃𝑦_𝑛−𝑐₃𝑦_𝑛−1+ℎ²𝑎₃₁𝑓_𝑛−1+𝑎₃₂𝑓_𝑛,(3.6)₄𝑌=1+𝑐₄𝑦_𝑛−𝑐₄𝑦_𝑛−1+ℎ²𝑎₄₁𝑓_𝑛−1+𝑎₄₂𝑓_𝑛+𝑎₄₃𝑓𝑥_𝑛+𝑐₃ℎ,𝑌₃,(3.7)₅𝑦=1+𝑐₅𝑦_𝑛−𝑐₅𝑦_𝑛−1+ℎ²𝑎₅₁𝑓_𝑛−1+𝑎₅₂𝑓_𝑛+𝑎₅₃𝑓𝑥_𝑛+𝑐₃ℎ,𝑌₃+𝑎₅₄𝑓𝑥_𝑛+𝑐₄ℎ,𝑌₄,(3.8)_𝑛+1=2𝑦_𝑛−𝑦_𝑛−1+ℎ²𝑏₁𝑓_𝑛−1+𝑏₂𝑓_𝑛+𝑏₃𝑓𝑥_𝑛+𝑐₃ℎ,𝑌₃+𝑏₄𝑓𝑥_𝑛+𝑐₄ℎ,𝑌₄+𝑏₅𝑓𝑥_𝑛+𝑐₅ℎ,𝑌₅(.3.9) Associate with (3.6), (3.7), and (3.8) the linear operator 𝐿, we have𝐿[ℎ,𝐚]𝑦(𝑥)=𝑦𝑥+𝑐_𝑖ℎ−1+𝑐_𝑖𝑦(𝑥)+𝑐_𝑖𝑦(𝑥−ℎ)−ℎ^2𝑖−1_𝑗=1𝑎_𝑖𝑗𝑦^𝑥+𝑐_𝑖ℎ,𝑖=3,4,5.(3.10) Note that the operator 𝐿 integrates exactly 1 and 𝑥 for any 𝑥 and ℎ. The coefficients of the variable step-size exponentially fitted hybrid method are obtained by solving equations arising from the choice of integers 𝑀, 𝐾, and 𝑃 described as follows.

Step 1. Consider the linear operator 𝐿 that associated with (3.6). Solving equations obtained by assuming that 𝑐₃=1/5 and imposing 𝐿 to integrate exactly 𝑥² and 𝑥³, we obtain 𝑎₃₁=4/125 and 𝑎₃₂=11/125 which is the same value as the coefficients of EHM6(4) method. This means that 𝑀=4. Choosing 𝐾=1,  𝑃=0 gives us {1,𝑥,𝑒^±𝜇𝑥}. Setting 𝑐₃=1/5 and imposing 𝐿 to integrate exactly 𝑒^±𝜇𝑥 leads the following equations: exp15𝜇ℎ−56+5exp1(𝜇ℎ)−ℎ²𝑎₃₁𝜇²exp(𝜇ℎ)−ℎ²exp−1𝑎₃₂𝜇²=0,5𝜇ℎ−65+15exp(𝜇ℎ)−ℎ²𝑎₃₁𝜇²exp(𝜇ℎ)−ℎ²𝑎₃₂𝜇²=0.(3.11) Solving these equations, we get coefficients 𝑎₃₁ and 𝑎₃₂.

Step 2. Consider the linear operator 𝐿 that associated with (3.7). Solving equations obtained by assuming that 𝑐₃=1/5,  𝑐₄=7/10,  𝑎₄₁=119/2000 and imposing 𝐿 to integrate exactly 𝑥² and 𝑥³, we obtain 𝑎₄₂=1071/2000 and 𝑎₄₃=0 which are of the same value as the coefficients of EHM6(4) method. This means that 𝑀=4. Choosing 𝐾=1, 𝑃=0 gives us {1,𝑥,𝑒^±𝜇𝑥}. Setting 𝑐₃=1/5,  𝑐₄=7/10,  𝑎₄₁=119/2000, and 𝐿[ℎ,𝐚]𝑒^±𝜇𝑥=0, we get the following equations: exp710𝜇ℎ−1170+10exp7(𝜇ℎ)−119ℎ²𝜇²2000exp(𝜇ℎ)−ℎ²𝑎₄₂𝜇²−ℎ²𝑎₄₃𝜇²exp115𝜇ℎ=0,exp((7/10)𝜇ℎ)−1170+170exp(𝜇ℎ)−2101090ℎ²𝜇²exp(𝜇ℎ)−ℎ²𝑎₄₂𝜇²−ℎ²𝑎₄₃𝜇²exp((1/5)𝜇ℎ)=0.(3.12) Solving these equations, we get 𝑎₄₂ and 𝑎₄₃.

Step 3. Consider the linear operator 𝐿 that associated with (3.8). Solving equations obtained by assuming that 𝑐₃=1/5,  𝑐₄=7/10,𝑐₅=−1/2,  𝑎₅₁=−11/204,  𝑎₅₂=−7/144 and imposing 𝐿 to integrate exactly 𝑥² and 𝑥³, we obtain 𝑎₅₃=−7/144 and 𝑎₅₄=4/153 which are of the same value as the coefficients of EHM6(4) method. This means that 𝑀=4. Choosing 𝐾=1 and 𝑃=0 gives us {1,𝑥,𝑒^±𝜇𝑥}. Setting 𝑐₃=1/5,  𝑐₄=7/10,  𝑐₅=−1/2,  𝑎₅₁=−11/204,  𝑎₅₂=−7/144, and 𝐿[ℎ,𝐚]𝑒^±𝜇𝑥=0 leads to the following equations: exp−12𝜇ℎ−21−2exp1(𝜇ℎ)+11ℎ²𝜇²204exp(𝜇ℎ)+1744ℎ²𝜇²−ℎ²𝑎₅₃𝜇²exp1−ℎ₂5𝜇ℎ𝑎₅₄𝜇²exp7110𝜇ℎ=0,exp((−1/2)𝜇ℎ)−12−12exp(𝜇ℎ)+21014ℎ²𝜇²exp(𝜇ℎ)+7144ℎ²𝜇²−ℎ²𝑎₅₃𝜇²−ℎexp((1/5)𝜇ℎ)₂𝑎₅₄𝜇²exp((7/10)𝜇ℎ)=0.(3.13) Solving these equations, we get 𝑎₅₃ and 𝑎₅₄.

Step 4. The linear operator associating with (3.9) is given as 𝐿[ℎ,𝐚]𝑦(𝑥)=−𝑦(ℎ𝑥₂+ℎ)−2𝑦(𝑥)+𝑦(𝑥−ℎ)𝑏₁𝑦^(𝑥−ℎ)+𝑏₂𝑦^(𝑥)+𝑏₃𝑦^𝑥+𝑐₃ℎ+𝑏₄𝑦^𝑥+𝑐₄ℎ+𝑏₅𝑦^𝑥+𝑐₅ℎ,(3.14) where 𝐚=[𝑏₁,𝑏₂,𝑏₃,𝑏₄,𝑏₅]. Solving equations obtained by assuming that 𝑐₃=1/5,  𝑐₄=7/10,  𝑐₅=−1/2 and imposing 𝐿 to integrate exactly 𝑥², 𝑥³, 𝑥⁴, 𝑥⁵ and 𝑥⁶, we obtain 𝑏₁=168,𝑏²=1142,𝑏³=2584,𝑏⁴=50357,𝑏⁵=27(3.15) which are of the same value as the coefficients of EHM6(4) method. This means that 𝑀=7. Choosing 𝐾=4 and 𝑃=0 gives us {1,𝑥,𝑥²,𝑥³,𝑥⁴,𝑒^±𝜇𝑥}. Setting 𝑐₃=1/5,𝑐₄=7/10,  𝑐₅=−1/2,  𝐿[ℎ,𝐚]𝑥²=0,𝐿[ℎ,𝐚]𝑥³=0,𝐿[ℎ,𝐚]𝑥⁴=0 and 𝐿[ℎ,𝐚]𝑒^±𝜇𝑥=0, the following equations are obtained: (𝑥(𝑥++ℎℎ)₃)²−2𝑥²+(𝑥−ℎ)²−ℎ²2𝑏₁+2𝑏₂+2𝑏₃+2𝑏₄+2𝑏₅=0,−2𝑥³+(𝑥−ℎ)³−+ℎ𝑏²₃𝑏¹(6𝑥−6ℎ)+6𝑏₂𝑥(𝑥+ℎ)₄6𝑥+56ℎ+𝑏⁴6𝑥+251ℎ+𝑏⁵(6𝑥−3ℎ)=0,−2𝑥⁴+(𝑥−ℎ)⁴−ℎ²12𝑏¹(𝑥−ℎ)²+12𝑏₂𝑥²+12𝑏₃𝑥+15ℎ²+12𝑏₄𝑥+710ℎ²+12𝑏₅𝑥−1exp(𝜇ℎ)−2+12ℎ²=0,exp(𝜇ℎ)−ℎ²𝑏₁𝜇²exp(𝜇ℎ)−ℎ²𝑏₂𝜇²−ℎ²𝑏₃𝜇²exp15𝜇ℎ−ℎ²𝑏₄𝜇²exp7−ℎ₂10𝜇ℎ𝑏₅𝜇²exp−112𝜇ℎ=0,exp(𝜇ℎ)−2+exp(𝜇ℎ)−ℎ²𝑏₁𝜇²exp(𝜇ℎ)−ℎ²𝑏₂𝜇²−ℎ²𝑏₃𝜇²exp((1/5)𝜇ℎ)−ℎ²𝑏₄𝜇²exp((7/10)𝜇ℎ)−ℎ²𝑏₅𝜇²exp((−1/2)𝜇ℎ)=0.(3.16) Solving these equations, we get 𝑏₁, 𝑏₂, 𝑏₃, 𝑏₄, and 𝑏₅.

Step 5. Consider 𝑦_𝑛+1=2𝑦_𝑛−𝑦_𝑛−1+ℎ²𝑏₁𝑓_𝑛−1+𝑏₂𝑓_𝑛+𝑏₃𝑓𝑥_𝑛+𝑐₃ℎ,𝑌₃+𝑏₄𝑓𝑥_𝑛+𝑐₄ℎ,𝑌₄.(3.17) Associate with (3.17) is the linear operator 𝐿: 𝐿[ℎ,𝐚]𝑦(𝑥)=−𝑦(ℎ𝑥₂+ℎ)−2𝑦(𝑥)+𝑦(𝑥−ℎ)𝑏₁𝑦^(𝑥−ℎ)+𝑏₂𝑦^(𝑥)+𝑏₃𝑦^𝑥+𝑐₃ℎ+𝑏₄𝑦^𝑥+𝑐₄ℎ,(3.18) where 𝐚=[𝑏₁,𝑏₂,𝑏₃,𝑏₄]. Solving equations obtained by assuming that 𝑐₃=1/5,  𝑐₄=7/10 and imposing 𝐿 to integrate exactly 𝑥²,𝑥³,𝑥⁴, and 𝑥⁵, we obtain 𝑏₁=5/68,𝑏₂=47/42,𝑏₃=−5/12,𝑏₄=80/357 which are of the same value as the coefficients of EHM6(4) method. This means that 𝑀=6. Choosing 𝐾=3 and 𝑃=0 gives us {𝑥²,𝑥³,𝑒^±𝜇𝑥}. Setting 𝑐₃=1/5,𝑐₄=7/10 and 𝐿[ℎ,𝐚]𝑥²=0,𝐿[ℎ,𝐚]𝑥³=0,𝐿[ℎ,𝐚]𝑒^±𝜇𝑥=0, the following equations are obtained: ((𝑥𝑥++ℎℎ))²₃−2𝑥²+(𝑥−ℎ)²−ℎ²2𝑏₁+2𝑏₂+2𝑏₃+2𝑏₄=0,−2𝑥³+(𝑥−ℎ)³−ℎ²𝑏₁(6𝑥−6ℎ)+6𝑏₂𝑥+𝑏₃6𝑥+65ℎ+𝑏₄6𝑥+21exp(𝜇ℎ)−2+15ℎ=0,exp(𝜇ℎ)−ℎ²𝑏₁𝜇²exp(𝜇ℎ)−ℎ²𝑏₂𝜇²−ℎ²𝑏₃𝜇²exp15𝜇ℎ−ℎ²𝑏₄𝜇²exp7110𝜇ℎ=0,exp(𝜇ℎ)−2+exp(𝜇ℎ)−ℎ²𝑏₁𝜇²exp(𝜇ℎ)−ℎ²𝑏₂𝜇²−ℎ²𝑏₃𝜇²exp−1−ℎ₂5𝜇ℎ𝑏₄𝜇²exp−710𝜇ℎ=0.(3.19) Solving these equations we get 𝑏₁, 𝑏₂, 𝑏₃ and 𝑏₄.

Let 𝑧=𝜇ℎ. The Taylor series expansions for coefficients of the variable step-size exponentially fitted explicit hybrid method:𝑎₃₁=4125−41678275𝑧²[phone]𝑧⁴−24188615234375𝑧⁶+2763464692138671875𝑧⁸𝑎−⋯,₃₂=11125−187735700𝑧²[phone]𝑧⁴−155785339375000000𝑧⁶+35428292388593750000000𝑧⁸𝑎−⋯,₄₂=10712000−81040309090𝑧²−1716269480000000𝑧⁴−7931374315760000000000𝑧⁶−90299601593456000000000000𝑧⁸𝑎−⋯,₄₃=658072400000𝑧²+24244211440000000𝑧⁴+2208900615760000000000𝑧⁶+472955352710368000000000000𝑧⁸𝑎+⋯,₅₃=−7144+7122070𝑧²[phone]𝑧⁴+125191559362880000000𝑧⁶+14183860314515200000000𝑧⁸𝑎+⋯,₅₄=4153−24818300𝑧²−4209143200000𝑧⁴−3046063207360000000𝑧⁶−127567691200000000𝑧⁸𝑏+⋯,₁=168−47160𝑧²[phone]𝑧⁴−47349392000000𝑧⁶+13599174231906500000000𝑧⁸𝑏−⋯,₂=1142−9380𝑧²[phone]𝑧⁴−804157624000000𝑧⁶+13599172904249600000000𝑧⁸𝑏+⋯,₃=2584+3192𝑧²−411975680𝑧⁴[phone]𝑧⁶−13599173485099520000000𝑧⁸𝑏−⋯,₄=50357−33132𝑧²[phone]𝑧⁴−47334574400000𝑧⁶+135991729623345920000000𝑧⁸𝑏+⋯,₅=27+9180𝑧²−414939200𝑧⁴+8041172872000000𝑧⁶−13599178712748800000000𝑧⁸−⋯,𝑏₁=568−628100𝑧²[phone]𝑧⁴−45671912500000𝑧⁶+46981667298960400000000000𝑧⁸−⋯,𝑏₂=4742+114300𝑧²−40051132300000𝑧⁴[phone]𝑧⁶−25871517139122245200000000000𝑧⁸+⋯,𝑏₃=−512−4100𝑧²[phone]𝑧⁴−124332700000000𝑧⁶+11777016978731800000000000𝑧⁸−⋯,𝑏₄=80357−219175𝑧²+1229331124550000𝑧⁴−4809722950000000𝑧⁶+304306431471039084200000000000𝑧⁸−⋯.(3.20) It is noted that, if 𝑧→0 then the coefficients of EHM6(4) method will be recovered. This is the important feature of the new method which will be denoted by EEHM6(4) method. In our program code, we use 𝜇 as an imaginary number and thus, 𝑧 can be written as 𝑧=𝑖𝑤ℎ. We first convert all coefficients (which are in term of 𝑧=𝜇ℎ) to trigonometric form by using computer algebra package Maple. Then, we substitute 𝜇=𝑖𝑤 into the resulting expressions. For small 𝑧, the coefficients are subject to heavy cancellations. Since 𝑧 varies throughout the integration, therefore it is convenient to use Taylor series expansions for the coefficients of the method.

4. Stability

Now our interest lies in the stability analysis for an exponentially fitted explicit hybrid method of the form (2.1) in which its coefficients are functions of one frequency 𝜇 and the step-length ℎ. Applying this method to the test problem 𝑦^=−𝜆²𝑦,𝜆>0,(4.1) we obtain 𝐘𝑦=(𝐞+𝐜)𝑦_𝑛−𝐜𝑦_𝑛−1−𝐻²𝐀𝐘,(4.2)_𝑛+1=2𝑦_𝑛−𝑦_𝑛−1−𝐻²𝐛^𝑇𝐘(4.3) which is written in vector form where 𝐻=𝜆ℎ, 𝐘=(𝑌₁,…,𝑌_𝑠)^𝑇, and 𝐞=(1,…,1)^𝑇. Solving for 𝐘 in (4.2) and then substitute it into (4.3) leads to the following recursion: 𝑦_𝑛+1−𝑆𝐻²,𝑧𝑦_𝑛+𝑃𝐻²,𝑧𝑦_𝑛−1=0,(4.4) where 𝑆(𝐻²,𝑧)=2−𝐻²𝐛^𝑇(𝐈+𝐻²𝐀)⁻¹(𝐞+𝐜) and 𝑃(𝐻²,𝑧)=1−𝐻²𝐛^𝑇(𝐈+𝐻²𝐀)⁻¹𝐜.

The characteristic equation associated with (4.4) is 𝜉²−𝑆𝐻²,𝑧𝜉+𝑃𝐻²,𝑧=0.(4.5) The subsequent definitions are introduced to provide the stability concept for exponentially fitted hybrid methods corresponding to the characteristic equation (4.5). Here, we follow the main ideas given by Coleman and Ixaru [3].

Definition 4.1. For exponentially fitted hybrid methods corresponding to the characteristic equation (4.5) and satisfying 𝑃(𝐻²,𝑧)=1, a region of stability Ω is a region of the 𝐻-𝑧 plane throughout which |𝑆(𝐻²,𝑧)|<2. Any closed curve defined by |𝑆(𝐻²,𝑧)|=2 is a stability boundary.

Definition 4.2. For exponentially fitted hybrid methods corresponding to the characteristic equation (4.5) and satisfying |𝑃(𝐻²,𝑧)|<1, a region of stability Ω is a region of the 𝐻-𝑧 plane such that, for every (𝐻,𝑧)∈Ω,||𝑆𝐻²,𝑧||<𝑃𝐻²,𝑧+1.(4.6) Any closed curve defined by 𝑃(𝐻²,𝑧)=1 and |𝑆(𝐻²,𝑧)|=𝑃(𝐻²,𝑧)+1 is a stability boundary.
In Figure 1, we illustrate the stability region of the exponentially fitted method which is based on the higher-order formula of EHM6(4) method whereas, in Figure 2, we show the stability region of the exponentially fitted method which is based on the lower-order formula of EHM6(4) method. From Figures 1 and 2, it is observed that the stability regions get smaller as 𝐻 increases.

(a)

(b)

(c)

Figure 1 

Stability region of the exponentially fitted method which is based on the higher-order formula of EHM6(4) method. (a) (𝐻=0.1, 𝑧∈𝐂), (b) (𝐻=0.5, 𝑧∈𝐂), (c) (𝐻=1, 𝑧∈𝐂).

(a)

(b)

(c)

Figure 2 

Stability region of the exponentially fitted method which is based on the lower-order formula of EHM6(4) method. (a) (𝐻=0.5, 𝑧∈𝐂), (b) (𝐻=1, 𝑧∈𝐂), (c) (𝐻=2, 𝑧∈𝐂).

5. Numerical Experiments

The following are abbreviations representing the codes that have been used in the comparisons.(i)EEHM6(4): Variable step-size exponentially fitted explicit hybrid method derived in this paper. This method requires four function evaluations at each step.(ii)EHM6(4): Embedded explicit hybrid method 6(4) pair derived in this paper. This method requires four function evaluations at each step.(iii)FRKN4(3): Modified embedded explicit RKN 4(3) pair proposed by Van de Vyver [10]. This method requires three function evaluations at each step.

We have applied all the above codes to four problems to evaluate their efficiency. For controlling the step-size in FRKN4(3) code, we use the step-size control procedure discussed in [10]. In Figures 3, 4, 5, and 6, the efficiency curves of natural logarithm of maximum global error (MAXGE) versus the computational effort measured by the natural logarithm of number of function evaluations (NFE) are depicted. The maximum global error (MAXGE) is given by this formula MAXGE=max‖‖𝑦𝑥_𝑛−𝑦_𝑛‖‖,(5.1) where 𝑦(𝑥_𝑛) is the true solution and 𝑦_𝑛 is the computed solution. In the following test problems, the tolerances used are tol=10^−2𝑖, 𝑖=1,…,6. The exact solutions are used to obtain all the starting values.

Figure 3 

Efficiency curves for Problem 1.

Figure 4 

Efficiency curves for Problem 2.

Figure 5 

Efficiency curves for Problem 3.

Figure 6 

Efficiency curves for Problem 4.

Problem 1 (perturbed system). Source: Franco [11]: 𝑦^1+100𝑦₁+2𝑦¹𝑦₂𝑦²¹+𝑦²²𝑦=𝑓₁(𝑥),𝑦₁(0)=1,𝑦^1(0)=𝜀,_2+25𝑦₂+𝑦²¹−𝑦²²𝑦²¹+𝑦²²=𝑓₂(𝑥),𝑦₂(0)=−𝜀,𝑦^2(0)=5(5.2) with 𝜀=10⁻³ and 𝑓₁(𝑥)=2cos(10𝑥)sin(5𝑥)+2𝜀(sin(5𝑥)sin(𝑥)−cos(10𝑥)cos(𝑥))−𝜀²sin(2𝑥)cos²(10𝑥)+sin²(5𝑥)+2𝜀(sin(𝑥)cos(10𝑥)−cos(𝑥)sin(5𝑥))+𝜀²𝑓+99𝜀sin(𝑥),₂(𝑥)=cos²(10𝑥)−sin²(5𝑥)+2𝜀(sin(𝑥)cos(10𝑥)+cos(𝑥)sin(5𝑥))−𝜀²cos(2𝑥)cos²(10𝑥)+sin²(5𝑥)+2𝜀(sin(𝑥)cos(10𝑥)−cos(𝑥)sin(5𝑥))+𝜀²−24𝜀cos(𝑥).(5.3) Solution: 𝑦₁(𝑥)=cos(10𝑥)+𝜀sin(𝑥),𝑦₂(𝑥)=sin(5𝑥)−𝜀cos(𝑥),0≤𝑥≤10.(5.4) For EEHM6(4) and FRKN4(3) codes, we choose 𝑧=10𝑖 for the first component whereas 𝑧=5𝑖 for the second component.

Problem 2 (linear oscillatory problem). Source: Franco [6]: 𝑦^1𝑦=−13𝑦₁+12𝑦₂+9cos(2𝑥)−12sin(2𝑥),𝑦₁(0)=1,𝑦^1(0)=−4,_2=12𝑦₁−13𝑦₂−12cos(2𝑥)+9sin(2𝑥),𝑦₂(0)=0,𝑦^2(0)=8(5.5) with the theoretical solution: 𝑦₁(𝑥)=sin(𝑥)−sin(5𝑥)+cos(2𝑥),𝑦₂(𝑥)=sin(𝑥)+sin(5𝑥)+sin(2𝑥),0≤𝑥≤10(5.6) For EEHM6(4) and FRKN4(3) codes, we choose 𝑧=5𝑖.

Problem 3 (perturbed system). Source: Hans Van de Vyver [10]: 𝑦𝑦^1+25𝑦₁+𝜀𝑦²¹+𝑦²²=𝜀𝑓₁(𝑥),𝑦₁(0)=1,𝑦^1(0)=0,_1+25𝑦₂+𝜀𝑦²¹+𝑦²²=𝜀𝑓₂(𝑥),𝑦₂(0)=𝜀,𝑦^1(0)=5(5.7) with 𝜀=10⁻³ and 𝑓𝑓₁₂(𝑥)=1+𝜀²+2𝜀sin5𝑥+𝑥²+2cos𝑥²+25−4𝑥²sin𝑥²,(𝑥)=1+𝜀²+2𝜀sin5𝑥+𝑥²−2sin𝑥²+25−4𝑥²cos𝑥².(5.8) Solution: 𝑦₁(𝑥)=cos(5𝑥)+𝜀sin𝑥²,𝑦₂(𝑥)=sin(5𝑥)+𝜀cos𝑥²,0≤𝑥≤5(5.9) For EEHM6(4) and FRKN4(3) codes, we choose 𝑧=5𝑖.

Problem 4 (the undamped Duffing equation). Source: Van de Vyver [10]: 𝑦^=−𝑦−𝑦³+𝐵cos(𝑣𝑥),𝑦(0)=0.200426728067,𝑦^(0)=0,(5.10) where 𝐵=1/500 and 𝑣=1.01. The exact solution computed by the Galerkin method with a precision 10⁻¹² of the coefficients is given by 𝑦(𝑥)=𝐴₁cos(𝑣𝑥)+𝐴₃cos(3𝑣𝑥)+𝐴₅cos(5𝑣𝑥)+𝐴₇cos(7𝑣𝑥)+𝐴₉cos(9𝑣𝑥),0≤𝑥≤20,(5.11) where 𝐴₁=0.200179477536, 𝐴₃=2.46946143⋅10⁻⁴, 𝐴₅=3.04014⋅10⁻⁷, 𝐴₇=3.74⋅10⁻¹⁰, and 𝐴₉=0.000000000000. For EEHM6(4) and FRKN4(3) codes, we choose 𝑧=𝑖.

From Figures 3–5, it is obvious that EEHM6(4) performs more efficiently than EHM6(4) and FRKN4(3) codes for solving Problems 1–3. For Problem 4, EEHM6(4) is as efficient as EHM6(4) code.

6. Conclusions

The embedding approach for hybrid methods has been devised by adopting the embedding concept in RKN methods. By using this embedding approach, we present Franco’s sixth-order explicit hybrid method for variable step-size code, denoted as EHM6(4) method. Based on EHM6(4) method, the new variable step-size exponentially fitted hybrid method has been developed and it is denoted as EEHM6(4) method. The numerical results show that EEHM6(4) method is efficient for the numerical solution of oscillatory problems. Furthermore, EEHM6(4) method is preferable to handle perturbed oscillators.