Table of Contents Author Guidelines Submit a Manuscript
International Journal of Mathematics and Mathematical Sciences
Volume 2011 (2011), Article ID 328197, 17 pages
http://dx.doi.org/10.1155/2011/328197
Research Article

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

Department of Mathematics, Faculty of Science, Universiti Putra Malaysia (UPM), Selangor, 43400 Serdang, Malaysia

Received 24 May 2011; Revised 9 September 2011; Accepted 9 September 2011

Academic Editor: Naseer Shahzad

Copyright © 2011 F. Samat et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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,𝑥,𝑥2,𝑥3,𝑥4,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=𝐲0,𝐲𝑥0=𝐲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 [13]) 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, 𝛼0, and 𝛽0 do not both vanish. Associate with the method is the linear operator𝐿[],𝐚𝑦(𝑥)=𝑘𝑗=0𝛼𝑗𝑦(𝑥+𝑗)2𝑘𝑗=0𝛽𝑗𝑦(𝑥+𝑗),(1.3) where 𝐚=[𝛼0,𝛼1,,𝛼𝑘,𝛽0,𝛽1,,𝛽𝑘]. 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,𝑥,𝑥2,,𝑥𝐾,𝑒±𝜇𝑥,𝑥𝑒±𝜇𝑥,𝑥2𝑒±𝜇𝑥,,𝑥𝑃𝑒±𝜇𝑥(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,𝑥,𝑥2,,𝑥𝑀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,𝑥,𝑥2,,𝑥𝐾,𝑒±𝜇𝑥,𝑥𝑒±𝜇𝑥,𝑥2𝑒±𝜇𝑥,,𝑥𝑃𝑒±𝜇𝑥(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,𝑌2=𝑦𝑛,𝑌𝑖=1+𝑐𝑖𝑦𝑛𝑐𝑖𝑦𝑛1+2𝑖1𝑗=1𝑎𝑖𝑗𝑓𝑥𝑛+𝑐𝑗,𝑌𝑗𝑦,𝑖=3,,𝑠,𝑛+1=2𝑦𝑛𝑦𝑛1+2𝑏1𝑓𝑛1+𝑏2𝑓𝑛+𝑠𝑖=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:100000𝑐00003𝑎31𝑎3200𝑐𝑠𝑎𝑠1𝑎𝑠2𝑎𝑠,𝑠10𝑏1𝑏2𝑏𝑠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<divtol, then 𝑛+1=𝑛,

(ii) if LTEtol/div, then 𝑛+1=2𝑛,

(iii) if LTEdivtol, then 𝑛+1=(1/2)𝑛 and repeat the step,

where, from numerical experiments, div is chosen to be 217. 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:4𝑖=1𝑏𝑖=1,4𝑖=1𝑏𝑖𝑐𝑖=0,4𝑖=1𝑏𝑖𝑐2𝑖=16,44𝑖=1𝑗=1𝑏𝑖𝑎𝑖𝑗=1,124𝑖=1𝑏𝑖𝑐3𝑖=0,44𝑖=1𝑗=1𝑏𝑖𝑐𝑖𝑎𝑖𝑗=1,1244𝑖=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=1,𝑐2=0,  𝑎21=0, and 𝑎𝑖𝑗=0  (𝑗𝑖). Substituting 𝑐1=1, 𝑐2=0,  𝑎21=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:𝑏1=5,68𝑏2=47,42𝑏35=,12𝑏4=80.357(3.2) The embedded 6(4) pair of hybrid methods is denoted as EHM6(4) method and represented by the tableau:1000000100000541251171250001011920001071120000002117204714441440153168114225845023577568475421280357(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: 10000001000005𝑎31𝑎32700010119𝑎200042𝑎431002117204𝑎14453𝑎540𝑏1𝑏2𝑏3𝑏4𝑏5(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=𝑦𝑛1,𝑌2=𝑦𝑛,𝑌(3.5)3=1+𝑐3𝑦𝑛𝑐3𝑦𝑛1+2𝑎31𝑓𝑛1+𝑎32𝑓𝑛𝑌,(3.6)4=1+𝑐4𝑦𝑛𝑐4𝑦𝑛1+2𝑎41𝑓𝑛1+𝑎42𝑓𝑛+𝑎43𝑓𝑥𝑛+𝑐3,𝑌3𝑌,(3.7)5=1+𝑐5𝑦𝑛𝑐5𝑦𝑛1+2𝑎51𝑓𝑛1+𝑎52𝑓𝑛+𝑎53𝑓𝑥𝑛+𝑐3,𝑌3+𝑎54𝑓𝑥𝑛+𝑐4,𝑌4,𝑦(3.8)𝑛+1=2𝑦𝑛𝑦𝑛1+2𝑏1𝑓𝑛1+𝑏2𝑓𝑛+𝑏3𝑓𝑥𝑛+𝑐3,𝑌3+𝑏4𝑓𝑥𝑛+𝑐4,𝑌4+𝑏5𝑓𝑥𝑛+𝑐5,𝑌5.(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 𝑐3=1/5 and imposing 𝐿 to integrate exactly 𝑥2 and 𝑥3, we obtain 𝑎31=4/125 and 𝑎32=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 𝑐3=1/5 and imposing 𝐿 to integrate exactly 𝑒±𝜇𝑥 leads the following equations: 1exp56𝜇5+15exp(𝜇)2𝑎31𝜇2exp(𝜇)2𝑎32𝜇21=0,exp56𝜇5+15exp(𝜇)2𝑎31𝜇2exp(𝜇)2𝑎32𝜇2=0.(3.11) Solving these equations, we get coefficients 𝑎31 and 𝑎32.

Step 2. Consider the linear operator 𝐿 that associated with (3.7). Solving equations obtained by assuming that 𝑐3=1/5,  𝑐4=7/10,  𝑎41=119/2000 and imposing 𝐿 to integrate exactly 𝑥2 and 𝑥3, we obtain 𝑎42=1071/2000 and 𝑎43=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 𝑐3=1/5,  𝑐4=7/10,  𝑎41=119/2000, and 𝐿[,𝐚]𝑒±𝜇𝑥=0, we get the following equations: 7exp10𝜇17+71010exp(𝜇)1192𝜇22000exp(𝜇)2𝑎42𝜇22𝑎43𝜇21exp51𝜇=0,exp((7/10)𝜇)17+71010exp(𝜇)11920002𝜇2exp(𝜇)2𝑎42𝜇22𝑎43𝜇2exp((1/5)𝜇)=0.(3.12) Solving these equations, we get 𝑎42 and 𝑎43.

Step 3. Consider the linear operator 𝐿 that associated with (3.8). Solving equations obtained by assuming that 𝑐3=1/5,  𝑐4=7/10,𝑐5=1/2,  𝑎51=11/204,  𝑎52=7/144 and imposing 𝐿 to integrate exactly 𝑥2 and 𝑥3, we obtain 𝑎53=7/144 and 𝑎54=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 𝑐3=1/5,  𝑐4=7/10,  𝑐5=1/2,  𝑎51=11/204,  𝑎52=7/144, and 𝐿[,𝐚]𝑒±𝜇𝑥=0 leads to the following equations: 1exp21𝜇21+2exp(𝜇)112𝜇2+7204exp(𝜇)1442𝜇22𝑎53𝜇21exp5𝜇2𝑎54𝜇27exp110𝜇=0,1exp((1/2)𝜇)212exp(𝜇)+112042𝜇27exp(𝜇)+1442𝜇22𝑎53𝜇2exp((1/5)𝜇)2𝑎54𝜇2exp((7/10)𝜇)=0.(3.13) Solving these equations, we get 𝑎53 and 𝑎54.

Step 4. The linear operator associating with (3.9) is given as 𝐿[],𝐚𝑦(𝑥)=𝑦(𝑥+)2𝑦(𝑥)+𝑦(𝑥)2𝑏1𝑦(𝑥)+𝑏2𝑦(𝑥)+𝑏3𝑦𝑥+𝑐3+𝑏4𝑦𝑥+𝑐4+𝑏5𝑦𝑥+𝑐5,(3.14) where 𝐚=[𝑏1,𝑏2,𝑏3,𝑏4,𝑏5]. Solving equations obtained by assuming that 𝑐3=1/5,  𝑐4=7/10,  𝑐5=1/2 and imposing 𝐿 to integrate exactly 𝑥2, 𝑥3, 𝑥4, 𝑥5 and 𝑥6, we obtain 𝑏1=168,𝑏2=1142,𝑏3=2584,𝑏4=50357,𝑏5=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,𝑥,𝑥2,𝑥3,𝑥4,𝑒±𝜇𝑥}. Setting 𝑐3=1/5,𝑐4=7/10,  𝑐5=1/2,  𝐿[,𝐚]𝑥2=0,𝐿[,𝐚]𝑥3=0,𝐿[,𝐚]𝑥4=0 and 𝐿[,𝐚]𝑒±𝜇𝑥=0, the following equations are obtained: (𝑥+)22𝑥2+(𝑥)222𝑏1+2𝑏2+2𝑏3+2𝑏4+2𝑏5=0,(𝑥+)32𝑥3+(𝑥)32𝑏1(6𝑥6)+6𝑏2𝑥+𝑏366𝑥+5+𝑏46𝑥+215+𝑏5(6𝑥3)=0,(𝑥+)42𝑥4+(𝑥)4212𝑏1(𝑥)2+12𝑏2𝑥2+12𝑏31𝑥+52+12𝑏47𝑥+102+12𝑏51𝑥221=0,exp(𝜇)2+exp(𝜇)2𝑏1𝜇2exp(𝜇)2𝑏2𝜇22𝑏3𝜇21exp5𝜇2𝑏4𝜇27exp10𝜇2𝑏5𝜇21exp21𝜇=0,exp(𝜇)2+exp(𝜇)2𝑏1𝜇2exp(𝜇)2𝑏2𝜇22𝑏3𝜇2exp((1/5)𝜇)2𝑏4𝜇2exp((7/10)𝜇)2𝑏5𝜇2exp((1/2)𝜇)=0.(3.16) Solving these equations, we get 𝑏1, 𝑏2, 𝑏3, 𝑏4, and 𝑏5.

Step 5. Consider 𝑦𝑛+1=2𝑦𝑛𝑦𝑛1+2𝑏1𝑓𝑛1+𝑏2𝑓𝑛+𝑏3𝑓𝑥𝑛+𝑐3,𝑌3+𝑏4𝑓𝑥𝑛+𝑐4,𝑌4.(3.17) Associate with (3.17) is the linear operator 𝐿: 𝐿[],𝐚𝑦(𝑥)=𝑦(𝑥+)2𝑦(𝑥)+𝑦(𝑥)2𝑏1𝑦(𝑥)+𝑏2𝑦(𝑥)+𝑏3𝑦𝑥+𝑐3+𝑏4𝑦𝑥+𝑐4,(3.18) where 𝐚=[𝑏1,𝑏2,𝑏3,𝑏4]. Solving equations obtained by assuming that 𝑐3=1/5,  𝑐4=7/10 and imposing 𝐿 to integrate exactly 𝑥2,𝑥3,𝑥4, and 𝑥5, we obtain 𝑏1=5/68,𝑏2=47/42,𝑏3=5/12,𝑏4=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 {𝑥2,𝑥3,𝑒±𝜇𝑥}. Setting 𝑐3=1/5,𝑐4=7/10 and 𝐿[,𝐚]𝑥2=0,𝐿[,𝐚]𝑥3=0,𝐿[,𝐚]𝑒±𝜇𝑥=0, the following equations are obtained: (𝑥+)22𝑥2+(𝑥)222𝑏1+2𝑏2+2𝑏3+2𝑏4=0,(𝑥+)32𝑥3+(𝑥)32𝑏1(6𝑥6)+6𝑏2𝑥+𝑏366𝑥+5+𝑏46𝑥+215e1=0,xp(𝜇)2+exp(𝜇)2𝑏1𝜇2exp(𝜇)2𝑏2𝜇22𝑏3𝜇21exp5𝜇2𝑏4𝜇27exp110𝜇=0,exp(𝜇)2+exp(𝜇)2𝑏1𝜇2exp(𝜇)2𝑏2𝜇22𝑏3𝜇21exp5𝜇2𝑏4𝜇27exp10𝜇=0.(3.19) Solving these equations we get 𝑏1, 𝑏2, 𝑏3 and 𝑏4.

Let 𝑧=𝜇. The Taylor series expansions for coefficients of the variable step-size exponentially fitted explicit hybrid method:𝑎31=4125172𝑧468752+1352𝑧3515625424188𝑧6152343756+2763464𝑧6921386718758𝑎,32=11125737𝑧1875002+11099𝑧2812500041557853𝑧393750000006+354282923𝑧885937500000008𝑎,42=1071200014399𝑧80000021716269𝑧4800000004793137431𝑧576000000000069029960159𝑧34560000000000008𝑎,43=65807𝑧24000002+2424421𝑧14400000004+220890061𝑧57600000000006+4729553527𝑧103680000000000008𝑎+,537=+144127𝑧72002+117733𝑧216000004+125191559𝑧3628800000006+141838603𝑧145152000000008𝑎+,54=4153413𝑧28800242091𝑧4320000043046063𝑧2073600000006127567𝑧6912000000008𝑏+,1=1168𝑧47602+41𝑧239904004473𝑧493920000006+1359917𝑧42319065000000008𝑏,2=11342𝑧9802+41𝑧164640048041𝑧576240000006+1359917𝑧29042496000000008𝑏+,3=25+184𝑧392241𝑧19756804+8041𝑧6914880000061359917𝑧34850995200000008𝑏,4=501357𝑧33322+41𝑧167932804473𝑧345744000006+1359917𝑧296233459200000008𝑏+,5=27+1𝑧980241𝑧49392004+8041𝑧17287200000061359917𝑧87127488000000008,𝑏1=56821𝑧68002+5141𝑧5355000044567𝑧19125000006+469816672𝑧989604000000000008,𝑏2=47+4213𝑧1400240051𝑧1323000004+24539𝑧2700000000625871517139𝑧1222452000000000008+,𝑏35=112𝑧4002+263𝑧2700000412433𝑧27000000006+1177701697𝑧87318000000000008,𝑏4=8035711𝑧29752+122933𝑧1124550000448097𝑧229500000006+30430643147𝑧10390842000000000008.(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 𝑦=𝜆2𝑦,𝜆>0,(4.1) we obtain 𝐘=(𝐞+𝐜)𝑦𝑛𝐜𝑦𝑛1𝐻2𝑦𝐀𝐘,(4.2)𝑛+1=2𝑦𝑛𝑦𝑛1𝐻2𝐛𝑇𝐘(4.3) which is written in vector form where 𝐻=𝜆, 𝐘=(𝑌1,,𝑌𝑠)𝑇, and 𝐞=(1,,1)𝑇. Solving for 𝐘 in (4.2) and then substitute it into (4.3) leads to the following recursion: 𝑦𝑛+1𝐻𝑆2𝑦,𝑧𝑛𝐻+𝑃2𝑦,𝑧𝑛1=0,(4.4) where 𝑆(𝐻2,𝑧)=2𝐻2𝐛𝑇(𝐈+𝐻2𝐀)1(𝐞+𝐜) and 𝑃(𝐻2,𝑧)=1𝐻2𝐛𝑇(𝐈+𝐻2𝐀)1𝐜.

The characteristic equation associated with (4.4) is 𝜉2𝐻𝑆2𝐻,𝑧𝜉+𝑃2,𝑧=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 𝑃(𝐻2,𝑧)=1, a region of stability Ω is a region of the 𝐻-𝑧 plane throughout which |𝑆(𝐻2,𝑧)|<2. Any closed curve defined by |𝑆(𝐻2,𝑧)|=2 is a stability boundary.

Definition 4.2. For exponentially fitted hybrid methods corresponding to the characteristic equation (4.5) and satisfying |𝑃(𝐻2,𝑧)|<1, a region of stability Ω is a region of the 𝐻-𝑧 plane such that, for every (𝐻,𝑧)Ω,||𝑆𝐻2||𝐻,𝑧<𝑃2,𝑧+1.(4.6) Any closed curve defined by 𝑃(𝐻2,𝑧)=1 and |𝑆(𝐻2,𝑧)|=𝑃(𝐻2,𝑧)+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.

fig1
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, 𝑧𝐂).
fig2
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=102𝑖, 𝑖=1,,6. The exact solutions are used to obtain all the starting values.

328197.fig.003
Figure 3: Efficiency curves for Problem 1.
328197.fig.004
Figure 4: Efficiency curves for Problem 2.
328197.fig.005
Figure 5: Efficiency curves for Problem 3.
328197.fig.006
Figure 6: Efficiency curves for Problem 4.

Problem 1 (perturbed system). Source: Franco [11]: 𝑦1+100𝑦1+2𝑦1𝑦2𝑦21+𝑦22=𝑓1(𝑥),𝑦1(0)=1,𝑦1𝑦(0)=𝜀,2+25𝑦2+𝑦21𝑦22𝑦21+𝑦22=𝑓2(𝑥),𝑦2(0)=𝜀,𝑦2(0)=5(5.2) with 𝜀=103 and 𝑓1(𝑥)=2cos(10𝑥)sin(5𝑥)+2𝜀(sin(5𝑥)sin(𝑥)cos(10𝑥)cos(𝑥))𝜀2sin(2𝑥)cos2(10𝑥)+sin2(5𝑥)+2𝜀(sin(𝑥)cos(10𝑥)cos(𝑥)sin(5𝑥))+𝜀2𝑓+99𝜀sin(𝑥),2(𝑥)=cos2(10𝑥)sin2(5𝑥)+2𝜀(sin(𝑥)cos(10𝑥)+cos(𝑥)sin(5𝑥))𝜀2cos(2𝑥)cos2(10𝑥)+sin2(5𝑥)+2𝜀(sin(𝑥)cos(10𝑥)cos(𝑥)sin(5𝑥))+𝜀224𝜀cos(𝑥).(5.3) Solution: 𝑦1(𝑥)=cos(10𝑥)+𝜀sin(𝑥),𝑦2(𝑥)=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𝑦1+12𝑦2+9cos(2𝑥)12sin(2𝑥),𝑦1(0)=1,𝑦1𝑦(0)=4,2=12𝑦113𝑦212cos(2𝑥)+9sin(2𝑥),𝑦2(0)=0,𝑦2(0)=8(5.5) with the theoretical solution: 𝑦1(𝑥)=sin(𝑥)sin(5𝑥)+cos(2𝑥),𝑦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𝑦1𝑦+𝜀21+𝑦22=𝜀𝑓1(𝑥),𝑦1(0)=1,𝑦1𝑦(0)=0,1+25𝑦2𝑦+𝜀21+𝑦22=𝜀𝑓2(𝑥),𝑦2(0)=𝜀,𝑦1(0)=5(5.7) with 𝜀=103 and 𝑓1(𝑥)=1+𝜀2+2𝜀sin5𝑥+𝑥2𝑥+2cos2+254𝑥2𝑥sin2,𝑓2(𝑥)=1+𝜀2+2𝜀sin5𝑥+𝑥2𝑥2sin2+254𝑥2𝑥cos2.(5.8) Solution: 𝑦1𝑥(𝑥)=cos(5𝑥)+𝜀sin2,𝑦2𝑥(𝑥)=sin(5𝑥)+𝜀cos2,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]: 𝑦=𝑦𝑦3+𝐵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 1012 of the coefficients is given by 𝑦(𝑥)=𝐴1cos(𝑣𝑥)+𝐴3cos(3𝑣𝑥)+𝐴5cos(5𝑣𝑥)+𝐴7cos(7𝑣𝑥)+𝐴9cos(9𝑣𝑥),0𝑥20,(5.11) where 𝐴1=0.200179477536, 𝐴3=2.46946143104, 𝐴5=3.04014107, 𝐴7=3.741010, and 𝐴9=0.000000000000. For EEHM6(4) and FRKN4(3) codes, we choose 𝑧=𝑖.

From Figures 35, it is obvious that EEHM6(4) performs more efficiently than EHM6(4) and FRKN4(3) codes for solving Problems 13. 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.

References

  1. T. E. Simos and J. Vigo-Aguiar, “A dissipative exponentially-fitted method for the numerical solution of the Schrödinger equation and related problems,” Computer Physics Communications, vol. 152, no. 3, pp. 274–294, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  2. T. E. Simos and P. S. Williams, “P-stable hybrid exponentially-fitted method for the numerical integration of the Schrodinger equation,” Computer Physics Communications, vol. 131, no. 1, pp. 109–119, 2000. View at Publisher · View at Google Scholar
  3. J. P. Coleman and L. Gr. Ixaru, “P-stability and exponential-fitting methods for y=f(x,y),” IMA Journal of Numerical Analysis, vol. 16, no. 2, pp. 179–199, 1996. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  4. G. Vanden Berghe and M. Van Daele, “Trigonometric polynomial or exponential fitting approach?” Journal of Computational and Applied Mathematics, vol. 233, no. 4, pp. 969–979, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  5. R. D'Ambrosio, E. Esposito, and B. Paternoster, “Exponentially fitted two-step hybrid methods for y=f(x,y),” Journal of Computational and Applied Mathematics, vol. 235, no. 16, pp. 4888–4897, 2011. View at Publisher · View at Google Scholar
  6. J. M. Franco, “A class of explicit two-step hybrid methods for second-order IVPs,” Journal of Computational and Applied Mathematics, vol. 187, no. 1, pp. 41–57, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  7. J. R. Dormand, M. E. A. El-Mikkawy, and P. J. Prince, “Families of Runge-Kutta-Nyström formulae,” IMA Journal of Numerical Analysis, vol. 7, no. 2, pp. 235–250, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  8. S. Filippi and J. Graf, “New Runge-Kutta-Nyström formula-pairs of order 8(7), 9(8), 10(9) and 11(10) for differential equations of the form y=f(x,y),” Journal of Computational and Applied Mathematics, vol. 14, no. 3, pp. 361–370, 1986. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  9. J. P. Coleman, “Order conditions for a class of two-step methods for y=f(x,y),” IMA Journal of Numerical Analysis, vol. 23, no. 2, pp. 197–220, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. H. Van de Vyver, “A Runge-Kutta-Nyström pair for the numerical integration of perturbed oscillators,” Computer Physics Communications, vol. 167, no. 2, pp. 129–142, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  11. J. M. Franco, “Exponentially fitted explicit Runge-Kutta-Nyström methods,” Journal of Computational and Applied Mathematics, vol. 167, no. 1, pp. 1–19, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH