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

Authors: F. Samat; F. Ismail; M. B. Suleiman
Journal: International Journal of Mathematics and Mathematical Sciences
 (2011)

---

## 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,x,exp(μx),exp(−μx)},μ∈C, while the last stage of this method integrates exactly systems whose solutions are linear combinations of  {1,x,x2,x3,x4,exp(μx),exp(−μx)}. 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.


---

## Body

## 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(1.1)y′′(x)=f(x,y(x)),y(x0)=y0,y′(x0)=y0′
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 lineark-step methods of the form(1.2)∑j=1kαjyn+j=h2∑j=1kβjyn+j′′,
where αk=1, α0, and β0 do not both vanish. Associate with the method is the linear operator(1.3)L[h,a]y(x)=∑j=0kαjy(x+jh)-h2∑j=0kβjy′′(x+jh),
where a=[α0,α1,…,αk,β0,β1,…,βk]. The linear operator L is said to integrate exactly the function y(x) if L[h,a]y(x)=0. We now impose the requirement that the operator L integrates exactly the following functions:(1.4)1,x,x2,…,xK,e±μx,xe±μx,x2e±μx,…,xPe±μx
with K+2P=M-3. Then, we solve the algebraic equations for αj and βj, j=0,1,…,k. In particular, the step-by-step exponential-fitting procedure is described as follows.Step 1.
Find the maximalM such that L integrates exactly the set of power functions {1,x,x2,…,xM-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.5)1,x,x2,…,xK,e±μx,xe±μx,x2e±μx,…,xPe±μx
with K+2P=M-3. Using M obtained from Step 1, we find suitable P and K which correspond to a finite set consisting of power functions and exponentials. Then, solve equations obtained by assuming that L 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]:
(2.1)Y1=yn-1,Y2=yn,Yi=(1+ci)yn-ciyn-1+h2∑j=1i-1aijf(xn+cjh,Yj),i=3,…,s,yn+1=2yn-yn-1+h2[b1fn-1+b2fn+∑i=3sbif(xn+cih,Yi)],
where fn-1 and fn represent f(xn-1,yn-1) and f(xn,yn), respectively. The method requires s-1 function evaluations or stages at each step of integration. The explicit hybrid method can be represented by the Butcher tableau:(2.2)-1000⋯00000⋯0c3a31a320⋯0⋮⋮⋮⋱⋮⋮csas1as2⋯as,s-10b1b2⋯bs-1bs=cAaabT
Adopting the concept of embedded p(q) pair of Runge-Kutta-Nystrom (RKN) methods (see [7, 8]), we define an embedded p(q) pair of hybrid methods to be based on the hybrid method (c,A,b) of order p and another hybrid method (c,A,b¯) of order q<p, represented by the following tableau:
(2.3)cAaabTbT¯
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
(2.4)LTE=‖yn+1-y¯n+1‖,
where yn+1 and y¯n+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) iftol/div<LTE<div·tol, then hn+1=hn,(ii) ifLTE≤tol/div, then hn+1=2hn,(iii) ifLTE≥div·tol, then hn+1=(1/2)hn 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 s=4 as given in [9] are listed as follows.Order conditions:(3.1)∑i=14b¯i=1,∑i=14b¯ici=0,∑i=14b¯ici2=16,∑i=14∑j=14b¯iaij=112,∑i=14b¯ici3=0,∑i=14∑j=14b¯iciaij=112,∑i=14∑j=14b¯iaijcj=0.
The two-step explicit hybrid methods (2.1) are a subclass of the two-step hybrid methods considered by Coleman [9] by taking c1=-1,c2=0,  a21=0, and aij=0  (j≥i). Substituting c1=-1, c2=0,  a21=0, and aij=0  (j≥i) together with coefficients ci and aij of the sixth-order formula into (3.1), we get the following values of coefficients b¯ for the three-stage fourth-order method:(3.2)b¯1=568,b¯2=4742,b¯3=-512,b¯4=80357.
The embedded 6(4) pair of hybrid methods is denoted as EHM6(4) method and represented by the tableau:(3.3)-10000000000015412511125000710119200010712000000-12-11204-7144-714441530∑avavaaa1681142258450357275684742-51280357


### 3.2. Variable Step-Size Exponentially Fitted Explicit Hybrid Method

We consider the four-stage explicit hybrid method given by the following tableau:(3.4)-10000000000015a31a320007101192000a42a4300-12-11204-7144a53a540b1b2b3b4b5
Note that the c-values and some of the A-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(3.5)Y1=yn-1,Y2=yn,(3.6)Y3=(1+c3)yn-c3yn-1+h2(a31fn-1+a32fn),(3.7)Y4=(1+c4)yn-c4yn-1+h2[a41fn-1+a42fn+a43f(xn+c3h,Y3)],(3.8)Y5=(1+c5)yn-c5yn-1+h2[a51fn-1+a52fn+a53f(xn+c3h,Y3)+a54f(xn+c4h,Y4)],(3.9)yn+1=2yn-yn-1+h2[b1fn-1+b2fn+b3f(xn+c3h,Y3)+b4f(xn+c4h,Y4)+b5f(xn+c5h,Y5)].
Associate with (3.6), (3.7), and (3.8) the linear operator L, we have(3.10)L[h,a]y(x)=y(x+cih)-(1+ci)y(x)+ciy(x-h)-h2∑j=1i-1aijy′′(x+cih),i=3,4,5.
Note that the operator L integrates exactly 1 and x for any x and h. The coefficients of the variable step-size exponentially fitted hybrid method are obtained by solving equations arising from the choice of integers M, K, and P described as follows.Step 1.
Consider the linear operatorL that associated with (3.6). Solving equations obtained by assuming that c3=1/5 and imposing L to integrate exactly x2 and x3, we obtain a31=4/125 and a32=11/125 which is the same value as the coefficients of EHM6(4) method. This means that M=4. Choosing K=1,  P=0 gives us {1,x,e±μx}. Setting c3=1/5 and imposing L to integrate exactly e±μx leads the following equations:
(3.11)exp⁡⁡(15μh)-65+15exp⁡(μh)-h2a31μ2exp⁡(μh)-h2a32μ2=0,exp⁡⁡(-15μh)-65+15exp⁡(μh)-h2a31μ2exp⁡(μh)-h2a32μ2=0.
Solving these equations, we get coefficients a31 and a32.Step 2.
Consider the linear operatorL that associated with (3.7). Solving equations obtained by assuming that c3=1/5,  c4=7/10,  a41=119/2000 and imposing L to integrate exactly x2 and x3, we obtain a42=1071/2000 and a43=0 which are of the same value as the coefficients of EHM6(4) method. This means that M=4. Choosing K=1, P=0 gives us {1,x,e±μx}. Setting c3=1/5,  c4=7/10,  a41=119/2000, and L[h,a]e±μx=0, we get the following equations:
(3.12)exp⁡(710μh)-1710+710exp⁡(μh)-119h2μ22000exp⁡(μh)-h2a42μ2-h2a43μ2exp⁡(15μh)=0,1exp⁡((7/10)μh)-1710+710exp⁡(μh)-1192000h2μ2exp⁡(μh)-h2a42μ2-h2a43μ2exp⁡((1/5)μh)=0.
Solving these equations, we get a42 and a43.Step 3.
Consider the linear operatorL that associated with (3.8). Solving equations obtained by assuming that c3=1/5,  c4=7/10,c5=-1/2,  a51=-11/204,  a52=-7/144 and imposing L to integrate exactly x2 and x3, we obtain a53=-7/144 and a54=4/153 which are of the same value as the coefficients of EHM6(4) method. This means that M=4. Choosing K=1 and P=0 gives us {1,x,e±μx}. Setting c3=1/5,  c4=7/10,  c5=-1/2,  a51=-11/204,  a52=-7/144, and L[h,a]e±μx=0 leads to the following equations:
(3.13)exp⁡(-12μh)-12-12exp⁡(μh)+11h2μ2204exp⁡(μh)+7144h2μ2-h2a53μ2exp⁡(15μh)-h2a54μ2exp⁡(710μh)=0,1exp⁡((-1/2)μh)-12-12exp⁡(μh)+11204h2μ2exp⁡(μh)+7144h2μ2-h2a53μ2exp⁡((1/5)μh)-h2a54μ2exp⁡((7/10)μh)=0.
Solving these equations, we get a53 and a54.Step 4.
The linear operator associating with (3.9) is given as
(3.14)L[h,a]y(x)=y(x+h)-2y(x)+y(x-h)-h2[b1y′′(x-h)+b2y′′(x)+b3y′′(x+c3h)+b4y′′(x+c4h)+b5y′′(x+c5h)],
where a=[b1,b2,b3,b4,b5]. Solving equations obtained by assuming that c3=1/5,  c4=7/10,  c5=-1/2 and imposing L to integrate exactly x2, x3, x4, x5 and x6, we obtain
(3.15)b1=168,b2=1142,b3=2584,b4=50357,b5=27
which are of the same value as the coefficients of EHM6(4) method. This means that M=7. Choosing K=4 and P=0 gives us {1,x,x2,x3,x4,e±μx}. Setting c3=1/5,c4=7/10,  c5=-1/2,  L[h,a]x2=0,L[h,a]x3=0,L[h,a]x4=0 and L[h,a]e±μx=0, the following equations are obtained:
(3.16)(x+h)2-2x2+(x-h)2-h2(2b1+2b2+2b3+2b4+2b5)=0,(x+h)3-2x3+(x-h)3-h2[(6x+215h)b1(6x-6h)+6b2x+b3(6x+65h)+b4(6x+215h)+b5(6x-3h)]=0,(x+h)4-2x4+(x-h)4-h2[(x-12h)212b1(x-h)2+12b2x2+12b3(x+15h)2+12b4(x+710h)2+12b5(x-12h)2]=0,exp⁡(μh)-2+1exp⁡(μh)-h2b1μ2exp⁡(μh)-h2b2μ2-h2b3μ2exp⁡(15μh)-h2b4μ2exp⁡(710μh)-h2b5μ2exp⁡(-12μh)=0,1exp⁡(μh)-2+exp⁡(μh)-h2b1μ2exp⁡(μh)-h2b2μ2-h2b3μ2exp⁡((1/5)μh)-h2b4μ2exp⁡((7/10)μh)-h2b5μ2exp⁡((-1/2)μh)=0.
Solving these equations, we get b1, b2, b3, b4, and b5.Step 5.
Consider(3.17)y¯n+1=2yn-yn-1+h2(b¯1fn-1+b¯2fn+b¯3f(xn+c3h,Y3)+b¯4f(xn+c4h,Y4)).
Associate with (3.17) is the linear operator L:
(3.18)L[h,a]y(x)=y(x+h)-2y(x)+y(x-h)-h2[b¯1y′′(x-h)+b¯2y′′(x)+b¯3y′′(x+c3h)+b¯4y′′(x+c4h)],
where a=[b¯1,b¯2,b¯3,b¯4]. Solving equations obtained by assuming that c3=1/5,  c4=7/10 and imposing L to integrate exactly x2,x3,x4, and x5, we obtain b¯1=5/68,b¯2=47/42,b¯3=-5/12,b¯4=80/357 which are of the same value as the coefficients of EHM6(4) method. This means that M=6. Choosing K=3 and P=0 gives us {x2,x3,e±μx}. Setting c3=1/5,c4=7/10 and L[h,a]x2=0,L[h,a]x3=0,L[h,a]e±μx=0, the following equations are obtained:
(3.19)(x+h)2-2x2+(x-h)2-h2(2b¯1+2b¯2+2b¯3+2b¯4)=0,(x+h)3-2x3+(x-h)3-h2(b¯1(6x-6h)+6b¯2x+b¯3(6x+65h)+b¯4(6x+215h))=0,exp(μh)-2+1exp⁡(μh)-h2b¯1μ2exp⁡(μh)-h2b¯2μ2-h2b¯3μ2exp⁡(15μh)-h2b¯4μ2exp⁡(710μh)=0,1exp⁡(μh)-2+exp⁡(μh)-h2b¯1μ2exp⁡(μh)-h2b¯2μ2-h2b¯3μ2exp⁡(-15μh)-h2b¯4μ2exp⁡(-710μh)=0.
Solving these equations we get b¯1, b¯2, b¯3 and b¯4.Letz=μh. The Taylor series expansions for coefficients of the variable step-size exponentially fitted explicit hybrid method:(3.20)a31=4125-17246875z2+13523515625z4-24188615234375z6+2763464692138671875z8-⋯,a32=11125-737187500z2+1109928125000z4-155785339375000000z6+35428292388593750000000z8-⋯,a42=10712000-14399800000z2-1716269480000000z4-7931374315760000000000z6-90299601593456000000000000z8-⋯,a43=658072400000z2+24244211440000000z4+2208900615760000000000z6+472955352710368000000000000z8+⋯,a53=-7144+1277200z2+11773321600000z4+125191559362880000000z6+14183860314515200000000z8+⋯,a54=4153-41328800z2-4209143200000z4-3046063207360000000z6-127567691200000000z8+⋯,b1=168-14760z2+4123990400z4-47349392000000z6+13599174231906500000000z8-⋯,b2=1142-3980z2+411646400z4-804157624000000z6+13599172904249600000000z8+⋯,b3=2584+1392z2-411975680z4+804169148800000z6-13599173485099520000000z8-⋯,b4=50357-13332z2+4116793280z4-47334574400000z6+135991729623345920000000z8+⋯,b5=27+1980z2-414939200z4+8041172872000000z6-13599178712748800000000z8-⋯,b¯1=568-216800z2+514153550000z4-45671912500000z6+46981667298960400000000000z8-⋯,b¯2=4742+131400z2-40051132300000z4+245392700000000z6-25871517139122245200000000000z8+⋯,b¯3=-512-1400z2+2632700000z4-124332700000000z6+11777016978731800000000000z8-⋯,b¯4=80357-112975z2+1229331124550000z4-4809722950000000z6+304306431471039084200000000000z8-⋯.
It is noted that, if z→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, z can be written as z=iwh. We first convert all coefficients (which are in term of z=μh) to trigonometric form by using computer algebra package Maple. Then, we substitute μ=iw into the resulting expressions. For small z, the coefficients are subject to heavy cancellations. Since z varies throughout the integration, therefore it is convenient to use Taylor series expansions for the coefficients of the method.


## 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 s=4 as given in [9] are listed as follows.Order conditions:(3.1)∑i=14b¯i=1,∑i=14b¯ici=0,∑i=14b¯ici2=16,∑i=14∑j=14b¯iaij=112,∑i=14b¯ici3=0,∑i=14∑j=14b¯iciaij=112,∑i=14∑j=14b¯iaijcj=0.
The two-step explicit hybrid methods (2.1) are a subclass of the two-step hybrid methods considered by Coleman [9] by taking c1=-1,c2=0,  a21=0, and aij=0  (j≥i). Substituting c1=-1, c2=0,  a21=0, and aij=0  (j≥i) together with coefficients ci and aij of the sixth-order formula into (3.1), we get the following values of coefficients b¯ for the three-stage fourth-order method:(3.2)b¯1=568,b¯2=4742,b¯3=-512,b¯4=80357.
The embedded 6(4) pair of hybrid methods is denoted as EHM6(4) method and represented by the tableau:(3.3)-10000000000015412511125000710119200010712000000-12-11204-7144-714441530∑avavaaa1681142258450357275684742-51280357


## 3.2. Variable Step-Size Exponentially Fitted Explicit Hybrid Method

We consider the four-stage explicit hybrid method given by the following tableau:(3.4)-10000000000015a31a320007101192000a42a4300-12-11204-7144a53a540b1b2b3b4b5
Note that the c-values and some of the A-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(3.5)Y1=yn-1,Y2=yn,(3.6)Y3=(1+c3)yn-c3yn-1+h2(a31fn-1+a32fn),(3.7)Y4=(1+c4)yn-c4yn-1+h2[a41fn-1+a42fn+a43f(xn+c3h,Y3)],(3.8)Y5=(1+c5)yn-c5yn-1+h2[a51fn-1+a52fn+a53f(xn+c3h,Y3)+a54f(xn+c4h,Y4)],(3.9)yn+1=2yn-yn-1+h2[b1fn-1+b2fn+b3f(xn+c3h,Y3)+b4f(xn+c4h,Y4)+b5f(xn+c5h,Y5)].
Associate with (3.6), (3.7), and (3.8) the linear operator L, we have(3.10)L[h,a]y(x)=y(x+cih)-(1+ci)y(x)+ciy(x-h)-h2∑j=1i-1aijy′′(x+cih),i=3,4,5.
Note that the operator L integrates exactly 1 and x for any x and h. The coefficients of the variable step-size exponentially fitted hybrid method are obtained by solving equations arising from the choice of integers M, K, and P described as follows.Step 1.
Consider the linear operatorL that associated with (3.6). Solving equations obtained by assuming that c3=1/5 and imposing L to integrate exactly x2 and x3, we obtain a31=4/125 and a32=11/125 which is the same value as the coefficients of EHM6(4) method. This means that M=4. Choosing K=1,  P=0 gives us {1,x,e±μx}. Setting c3=1/5 and imposing L to integrate exactly e±μx leads the following equations:
(3.11)exp⁡⁡(15μh)-65+15exp⁡(μh)-h2a31μ2exp⁡(μh)-h2a32μ2=0,exp⁡⁡(-15μh)-65+15exp⁡(μh)-h2a31μ2exp⁡(μh)-h2a32μ2=0.
Solving these equations, we get coefficients a31 and a32.Step 2.
Consider the linear operatorL that associated with (3.7). Solving equations obtained by assuming that c3=1/5,  c4=7/10,  a41=119/2000 and imposing L to integrate exactly x2 and x3, we obtain a42=1071/2000 and a43=0 which are of the same value as the coefficients of EHM6(4) method. This means that M=4. Choosing K=1, P=0 gives us {1,x,e±μx}. Setting c3=1/5,  c4=7/10,  a41=119/2000, and L[h,a]e±μx=0, we get the following equations:
(3.12)exp⁡(710μh)-1710+710exp⁡(μh)-119h2μ22000exp⁡(μh)-h2a42μ2-h2a43μ2exp⁡(15μh)=0,1exp⁡((7/10)μh)-1710+710exp⁡(μh)-1192000h2μ2exp⁡(μh)-h2a42μ2-h2a43μ2exp⁡((1/5)μh)=0.
Solving these equations, we get a42 and a43.Step 3.
Consider the linear operatorL that associated with (3.8). Solving equations obtained by assuming that c3=1/5,  c4=7/10,c5=-1/2,  a51=-11/204,  a52=-7/144 and imposing L to integrate exactly x2 and x3, we obtain a53=-7/144 and a54=4/153 which are of the same value as the coefficients of EHM6(4) method. This means that M=4. Choosing K=1 and P=0 gives us {1,x,e±μx}. Setting c3=1/5,  c4=7/10,  c5=-1/2,  a51=-11/204,  a52=-7/144, and L[h,a]e±μx=0 leads to the following equations:
(3.13)exp⁡(-12μh)-12-12exp⁡(μh)+11h2μ2204exp⁡(μh)+7144h2μ2-h2a53μ2exp⁡(15μh)-h2a54μ2exp⁡(710μh)=0,1exp⁡((-1/2)μh)-12-12exp⁡(μh)+11204h2μ2exp⁡(μh)+7144h2μ2-h2a53μ2exp⁡((1/5)μh)-h2a54μ2exp⁡((7/10)μh)=0.
Solving these equations, we get a53 and a54.Step 4.
The linear operator associating with (3.9) is given as
(3.14)L[h,a]y(x)=y(x+h)-2y(x)+y(x-h)-h2[b1y′′(x-h)+b2y′′(x)+b3y′′(x+c3h)+b4y′′(x+c4h)+b5y′′(x+c5h)],
where a=[b1,b2,b3,b4,b5]. Solving equations obtained by assuming that c3=1/5,  c4=7/10,  c5=-1/2 and imposing L to integrate exactly x2, x3, x4, x5 and x6, we obtain
(3.15)b1=168,b2=1142,b3=2584,b4=50357,b5=27
which are of the same value as the coefficients of EHM6(4) method. This means that M=7. Choosing K=4 and P=0 gives us {1,x,x2,x3,x4,e±μx}. Setting c3=1/5,c4=7/10,  c5=-1/2,  L[h,a]x2=0,L[h,a]x3=0,L[h,a]x4=0 and L[h,a]e±μx=0, the following equations are obtained:
(3.16)(x+h)2-2x2+(x-h)2-h2(2b1+2b2+2b3+2b4+2b5)=0,(x+h)3-2x3+(x-h)3-h2[(6x+215h)b1(6x-6h)+6b2x+b3(6x+65h)+b4(6x+215h)+b5(6x-3h)]=0,(x+h)4-2x4+(x-h)4-h2[(x-12h)212b1(x-h)2+12b2x2+12b3(x+15h)2+12b4(x+710h)2+12b5(x-12h)2]=0,exp⁡(μh)-2+1exp⁡(μh)-h2b1μ2exp⁡(μh)-h2b2μ2-h2b3μ2exp⁡(15μh)-h2b4μ2exp⁡(710μh)-h2b5μ2exp⁡(-12μh)=0,1exp⁡(μh)-2+exp⁡(μh)-h2b1μ2exp⁡(μh)-h2b2μ2-h2b3μ2exp⁡((1/5)μh)-h2b4μ2exp⁡((7/10)μh)-h2b5μ2exp⁡((-1/2)μh)=0.
Solving these equations, we get b1, b2, b3, b4, and b5.Step 5.
Consider(3.17)y¯n+1=2yn-yn-1+h2(b¯1fn-1+b¯2fn+b¯3f(xn+c3h,Y3)+b¯4f(xn+c4h,Y4)).
Associate with (3.17) is the linear operator L:
(3.18)L[h,a]y(x)=y(x+h)-2y(x)+y(x-h)-h2[b¯1y′′(x-h)+b¯2y′′(x)+b¯3y′′(x+c3h)+b¯4y′′(x+c4h)],
where a=[b¯1,b¯2,b¯3,b¯4]. Solving equations obtained by assuming that c3=1/5,  c4=7/10 and imposing L to integrate exactly x2,x3,x4, and x5, we obtain b¯1=5/68,b¯2=47/42,b¯3=-5/12,b¯4=80/357 which are of the same value as the coefficients of EHM6(4) method. This means that M=6. Choosing K=3 and P=0 gives us {x2,x3,e±μx}. Setting c3=1/5,c4=7/10 and L[h,a]x2=0,L[h,a]x3=0,L[h,a]e±μx=0, the following equations are obtained:
(3.19)(x+h)2-2x2+(x-h)2-h2(2b¯1+2b¯2+2b¯3+2b¯4)=0,(x+h)3-2x3+(x-h)3-h2(b¯1(6x-6h)+6b¯2x+b¯3(6x+65h)+b¯4(6x+215h))=0,exp(μh)-2+1exp⁡(μh)-h2b¯1μ2exp⁡(μh)-h2b¯2μ2-h2b¯3μ2exp⁡(15μh)-h2b¯4μ2exp⁡(710μh)=0,1exp⁡(μh)-2+exp⁡(μh)-h2b¯1μ2exp⁡(μh)-h2b¯2μ2-h2b¯3μ2exp⁡(-15μh)-h2b¯4μ2exp⁡(-710μh)=0.
Solving these equations we get b¯1, b¯2, b¯3 and b¯4.Letz=μh. The Taylor series expansions for coefficients of the variable step-size exponentially fitted explicit hybrid method:(3.20)a31=4125-17246875z2+13523515625z4-24188615234375z6+2763464692138671875z8-⋯,a32=11125-737187500z2+1109928125000z4-155785339375000000z6+35428292388593750000000z8-⋯,a42=10712000-14399800000z2-1716269480000000z4-7931374315760000000000z6-90299601593456000000000000z8-⋯,a43=658072400000z2+24244211440000000z4+2208900615760000000000z6+472955352710368000000000000z8+⋯,a53=-7144+1277200z2+11773321600000z4+125191559362880000000z6+14183860314515200000000z8+⋯,a54=4153-41328800z2-4209143200000z4-3046063207360000000z6-127567691200000000z8+⋯,b1=168-14760z2+4123990400z4-47349392000000z6+13599174231906500000000z8-⋯,b2=1142-3980z2+411646400z4-804157624000000z6+13599172904249600000000z8+⋯,b3=2584+1392z2-411975680z4+804169148800000z6-13599173485099520000000z8-⋯,b4=50357-13332z2+4116793280z4-47334574400000z6+135991729623345920000000z8+⋯,b5=27+1980z2-414939200z4+8041172872000000z6-13599178712748800000000z8-⋯,b¯1=568-216800z2+514153550000z4-45671912500000z6+46981667298960400000000000z8-⋯,b¯2=4742+131400z2-40051132300000z4+245392700000000z6-25871517139122245200000000000z8+⋯,b¯3=-512-1400z2+2632700000z4-124332700000000z6+11777016978731800000000000z8-⋯,b¯4=80357-112975z2+1229331124550000z4-4809722950000000z6+304306431471039084200000000000z8-⋯.
It is noted that, if z→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, z can be written as z=iwh. We first convert all coefficients (which are in term of z=μh) to trigonometric form by using computer algebra package Maple. Then, we substitute μ=iw into the resulting expressions. For small z, the coefficients are subject to heavy cancellations. Since z 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 h. Applying this method to the test problem
(4.1)y′′=-λ2y,λ>0,
we obtain
(4.2)Y=(e+c)yn-cyn-1-H2AY,(4.3)yn+1=2yn-yn-1-H2bTY
which is written in vector form where H=λh, Y=(Y1,…,Ys)T, and e=(1,…,1)T. Solving for Y in (4.2) and then substitute it into (4.3) leads to the following recursion:
(4.4)yn+1-S(H2,z)yn+P(H2,z)yn-1=0,
where S(H2,z)=2-H2bT(I+H2A)-1(e+c) and P(H2,z)=1-H2bT(I+H2A)-1c.The characteristic equation associated with (4.4) is
(4.5)ξ2-S(H2,z)ξ+P(H2,z)=0.
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 P(H2,z)=1, a region of stability Ω is a region of the H-z plane throughout which |S(H2,z)|<2. Any closed curve defined by |S(H2,z)|=2 is a stability boundary.Definition 4.2.
For exponentially fitted hybrid methods corresponding to the characteristic equation (4.5) and satisfying |P(H2,z)|<1, a region of stability Ω is a region of the H-z plane such that, for every (H,z)∈Ω,(4.6)|S(H2,z)|<P(H2,z)+1.
Any closed curve defined by P(H2,z)=1 and |S(H2,z)|=P(H2,z)+1 is a stability boundary.
In Figure1, 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 H increases.Stability region of the exponentially fitted method which is based on the higher-order formula of EHM6(4) method. (a) (H=0.1, z∈C), (b) (H=0.5, z∈C), (c) (H=1, z∈C).
(a)(b)(c)Stability region of the exponentially fitted method which is based on the lower-order formula of EHM6(4) method. (a) (H=0.5, z∈C), (b) (H=1, z∈C), (c) (H=2, z∈C).
(a)(b)(c)


## 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
(5.1)MAXGE=max⁡(‖y(xn)-yn‖),
where y(xn) is the true solution and yn is the computed solution. In the following test problems, the tolerances used are tol=10-2i, i=1,…,6. The exact solutions are used to obtain all the starting values.Figure 3
Efficiency curves for Problem1.Figure 4
Efficiency curves for Problem2.Figure 5
Efficiency curves for Problem3.Figure 6
Efficiency curves for Problem4.Problem 1 (perturbed system).
Source: Franco [11]:
(5.2)y1′′+100y1+2y1y2y12+y22=f1(x),y1(0)=1,y1′(0)=ε,y2′′+25y2+y12-y22y12+y22=f2(x),y2(0)=-ε,y2′(0)=5
with ε=10-3 and
(5.3)f1(x)=2cos⁡(10x)sin⁡(5x)+2ε(sin⁡(5x)sin⁡(x)-cos⁡(10x)cos⁡(x))-ε2sin⁡(2x)cos⁡2(10x)+sin⁡2(5x)+2ε(sin⁡(x)cos⁡(10x)-cos⁡(x)sin⁡(5x))+ε2+99εsin⁡(x),f2(x)=cos⁡2(10x)-sin⁡2(5x)+2ε(sin⁡(x)cos⁡(10x)+cos⁡(x)sin⁡(5x))-ε2cos⁡(2x)cos⁡2(10x)+sin⁡2(5x)+2ε(sin⁡(x)cos⁡(10x)-cos⁡(x)sin⁡(5x))+ε2-24εcos⁡(x).
Solution:
(5.4)y1(x)=cos⁡⁡(10x)+εsin⁡(x),y2(x)=sin⁡(5x)-εcos⁡(x),0≤x≤10.
For EEHM6(4) and FRKN4(3) codes, we choose z=10i for the first component whereas z=5i for the second component.Problem 2 (linear oscillatory problem).
Source: Franco [6]:
(5.5)y1′′=-13y1+12y2+9cos⁡(2x)-12sin⁡(2x),y1(0)=1,y1′(0)=-4,y2′′=12y1-13y2-12cos⁡(2x)+9sin⁡(2x),y2(0)=0,y2′(0)=8
with the theoretical solution:
(5.6)y1(x)=sin⁡(x)-sin⁡(5x)+cos⁡(2x),y2(x)=sin⁡(x)+sin⁡(5x)+sin⁡(2x),0≤x≤10
For EEHM6(4) and FRKN4(3) codes, we choose z=5i.Problem 3 (perturbed system).
Source: Hans Van de Vyver [10]:
(5.7)y1′′+25y1+ε(y12+y22)=εf1(x),y1(0)=1,y1′(0)=0,y1′′+25y2+ε(y12+y22)=εf2(x),y2(0)=ε,y1′(0)=5
with ε=10-3 and
(5.8)f1(x)=1+ε2+2εsin⁡(5x+x2)+2cos⁡(x2)+(25-4x2)sin⁡(x2),f2(x)=1+ε2+2εsin⁡(5x+x2)-2sin⁡(x2)+(25-4x2)cos⁡(x2).
Solution:
(5.9)y1(x)=cos⁡(5x)+εsin⁡(x2),y2(x)=sin⁡(5x)+εcos⁡(x2),0≤x≤5
For EEHM6(4) and FRKN4(3) codes, we choose z=5i.Problem 4 (the undamped Duffing equation).
Source: Van de Vyver [10]:
(5.10)y′′=-y-y3+Bcos⁡(vx),y(0)=0.200426728067,y′(0)=0,
where B=1/500 and v=1.01. The exact solution computed by the Galerkin method with a precision 10-12 of the coefficients is given by
(5.11)y(x)=A1cos⁡(vx)+A3cos⁡(3vx)+A5cos⁡(5vx)+A7cos⁡(7vx)+A9cos⁡(9vx),0≤x≤20,
where A1=0.200179477536, A3=2.46946143·10-4, A5=3.04014·10-7, A7=3.74·10-10, and A9=0.000000000000. For EEHM6(4) and FRKN4(3) codes, we choose z=i.From Figures3–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.



---

*Source: 328197-2011-11-30.xml*