Exponentially-Fitted Fourth-Derivative Single-Step Obrechkoff Method for Oscillatory/Periodic Problems

Abstract

The quest for accurate and more efficient methods for solving periodic/oscillatory problems is gaining more attention in recent time. This paper presents the construction and implementation of a family of exponentially-fitted Obrechkoff methods using a six-step flowchart discussed in the literature. A single-step Obrechkoff method involving terms up to the fourth derivative was used as the base method. We also present the stability and convergence properties of the constructed family of methods. Two numerical examples were used to illustrate the performance of the constructed methods.

1. Introduction

Ordinary Differential Equations (ODEs) that exhibit pronounced oscillatory or periodic behaviour in their solutions are often encountered in fields such as chemistry, engineering, electronic, mechanics and astrophysics [1,2,3]. Many of the classical methods for solving prominent classes of problems in ODE have been developed using only monomials as a basis [1,4,5,6,7,8]. In [9], the authors used a hybrid method to examine the direct solution of higher order (second, third and fourth order) initial value problem (IVP) of ordinary differential equations. However, in practice, many classical methods usually perform poorly when applied to problems with pronounced periodic or oscillatory behaviour in their solution [3,10,11]. This is due to the fact that for better accuracy to be achieved, a very small step size would be required with corresponding decrease in performance, especially in terms of efficiency [1]. One way to overcome this barrier is to adapt classical methods for such problems. The adaptation which is called “exponential/trigonometric fitting” involves the replacement of some of the highest order monomials of the basis by exponentials or trigonometric [3,11].Detailed analysis of the oscillation–preserving behaviour of some existing RKN–type methods were analysed from the point of view of geometric integration in [12]. Authors in [13] presented surveys on recent advances in the allied challenges of discretizing highly oscillatory ordinary differential equations and computing numerical quadrature of highly oscillatory integrals. They also attempted to sketch the mathematical foundations of a general approach to these issues [13] .

A pioneering work in the use of exponentially–fitted formulae for differential equations was conducted by [14]. The authors in [14] constructed integration formulae which contains free parameters - chosen so that a given function $exp\left(q\right)$ where q is real, satisfies the integration formulae exactly. The proposed method in [14] was on a one-step formulae, however, in [15], A-stable fourth order exponentially-fitted formulae based on a linear two-step formula was derived. Using the concept proposed in [15], the author in [16] proposed a Multiderivative Linear Multistep Method (MLMM) with k = 1 in the second derivative formulae. Many authors have proposed specially adapted Runge–-Kutta (RK) algorithms to solve this class of problems [17,18,19,20]. In this direction, exponentially–fitted RK (EFRK) methods which integrate exactly first–order systems whose solutions can be expressed as linear combinations of functions of the form $\{exp(\lambda t),exp(-\lambda t\left)\right\}$ or $\{cos(\omega t),sin(\omega t\left)\right\}$ were introduced in [21,22]. The construction of an implicit trigonometrically-fitted single-step method with a second derivative using trigonometric basis function was proposed in [23].

In this work, we used the six–step flowchart described in [3] to construct a class of exponentially–fitted single-step fourth-derivative Obrechkoff methods suitable for solving

$$\begin{array}{c}\hfill y^{\prime }=f(x,y),x\in \left[x_0,X\right],y\left(x_0\right)=y_0.\end{array}$$

(1)

2. Construction of Method

A classical fourth-derivative single-step Obrechkoff method for solving the first order initial value problem (1) can generally be written as

$$\begin{array}{cc}\hfill y_{j+1}& =a_0{y}_j+h\left(b_1{f}_{j+1}+b_0{f}_j\right)+h^2\left(c_1{f}_{j+1}^{\prime }+c_0{f}_j^{\prime }\right)+\hfill \\ \hfill & h^3\left(d_1{f}_{j+1}^{″}+d_0{f}_j^{″}\right)+h^4\left(e_1{f}_{j+1}^{\prime \prime \prime }+e_0{f}_j^{\prime \prime \prime }\right)\hfill \end{array}$$

(2)

where $a_0,b_0,b_1,c_0,c_1,d_0,d_1,e_0$ and $e_1$ are coefficients to be determined. Here, we present the construction of the exponentially-fitted variants of (2) using the six-step flowchart described in [3]. Following the six-step flowchart, the corresponding linear difference operator $\mathcal{L}[h,\mathbf{a}]$ is obtained as

$$\begin{array}{cc}\hfill \mathcal{L}[h,\mathbf{a}]y\left(x\right)& =y(x+h)-a_{0}y\left(x\right)-h\left(b_1{y}^{\prime }(x+h)+b_0{y}^{\prime }\left(x\right)\right)-\hfill \\ \hfill & h^2\left(c_1{y}^{″}(x+h)+c_0{y}^{″}\left(x\right)\right)-\hfill \\ \hfill & h^3\left(d_1{y}^{\prime \prime \prime }(x+h)+d_0{y}^{\prime \prime \prime }\left(x\right)\right)-\hfill \\ \hfill & h^4\left(e_1{y}^{\prime \prime \prime \prime }(x+h)+e_0{y}^{\prime \prime \prime \prime }\left(x\right)\right)\hfill \end{array}$$

(3)

where $\mathbf{a}:=(a_0,b_0,b_1,c_0,c_1,d_0,d_1,e_0,e_1)$. Step II of the procedure requires that we get the maximum value of M such that the algebraic system

$$\left\{L_m^{*}\left(\mathbf{a}\right)=h^{-m}\mathcal{L}[h,\mathbf{a}]x^m{|}_{x=0}=0|m=0,1,2,\cdots ,M-1\right\}$$

can be solved. The above results in

$$\begin{array}{ccc}\hfill L_0^{*}\left(\mathbf{a}\right)& =& 1-a_0=0\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill L_1^{*}\left(\mathbf{a}\right)& =& -b_0-b_1+1=0\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill L_2^{*}\left(\mathbf{a}\right)& =& -2b_1-2c_0-2c_1+1=0\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill L_3^{*}\left(\mathbf{a}\right)& =& -3b_1-6c_1-6d_0-6d_1+1=0\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill L_4^{*}\left(\mathbf{a}\right)& =& -4b_1-12c_1-24d_1-24e_0-24e_1+1=0\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill L_5^{*}\left(\mathbf{a}\right)& =& -5b_1-20c_1-60d_1-120e_1+1=0\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill L_6^{*}\left(\mathbf{a}\right)& =& -6b_1-30c_1-120d_1-360e_1+1=0\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill L_7^{*}\left(\mathbf{a}\right)& =& -7b_1-42c_1-210d_1-840e_1+1=0\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill L_8^{*}\left(\mathbf{a}\right)& =& -8b_1-56c_1-336d_1-1680e_1+1=0\hfill \end{array}$$

(12)

and the algebraic system is compatible when $M=9$. In addition, the solution only results in the coefficients of the associated classical method to be adapted. To exponentially fit the associated classical method, we proceed to step III of the six–step flowchart and obtain expressions for $G^{+}(Z,\mathbf{a})$ and $G^{-}(Z,\mathbf{a})$ which are respectively defined as

$$\begin{array}{ccc}\hfill G^{+}(Z,\mathbf{a})& =& \frac{1}{2}\left(E_0^{*}(z,\mathbf{a})+E_0^{*}(-z,\mathbf{a})\right)\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill G^{-}(Z,\mathbf{a})& =& \frac{1}{2z}\left(E_0^{*}(z,\mathbf{a})-E_0^{*}(-z,\mathbf{a})\right)\hfill \end{array}$$

(14)

where $E_0^{*}(\pm z,\mathbf{a})=e^{\mp \omega x}\mathcal{L}[h,\mathbf{a}]e^{\pm \omega x}$ and $Z=z^2$. The expressions for $G^{+}(Z,\mathbf{a})$ and $G^{-}(Z,\mathbf{a})$ are respectively obtained as

$$\begin{array}{ccc}\hfill G^{+}(Z,\mathbf{a})& =& -a_0+sinh\left(\sqrt{Z}\right)\left(-b_1\sqrt{Z}-d_1{Z}^{3/2}\right)+\hfill \\ & & cosh\left(\sqrt{Z}\right)\left(-c_{1}Z-e_1{Z}^2+1\right)-c_{0}Z-e_0{Z}^2\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill G^{-}(Z,\mathbf{a})& =& cosh\left(\sqrt{Z}\right)\left(-b_1-d_{1}Z\right)-b_0+\hfill \\ & & sinh\left(\sqrt{Z}\right)\left(-c_1\sqrt{Z}-e_1{Z}^{3/2}+\frac{1}{\sqrt{Z}}\right)-d_{0}Z\hfill \end{array}$$

(16)

where $\omega$, the frequency of oscillation is real or imaginary, $z=\omega h={\omega }_h$. (For the trigonometric case, i.e., $\omega$ is imaginary, we choose $z=\omega h=i\mu h$, i.e., $z^2=-{\mu }^2{h}^2=Z$.)

To implement step IV, consider the reference set of M functions:

$$\left\{1,x,\cdots ,x^K,exp(\pm \omega x),xexp(\pm \omega x),\cdots ,x^{P}exp(\pm \omega x)\right\}$$

with $K+2P=M-3$. Since for our method $M=9$, we have five possibilities, which we shall respectively refer to as S1, S2, S3, S4 and S5:

• S1: $K=8,P=-1$, the classical case with the set

  $$1,x,x^2,x^3,x^4,x^5,x^6,x^7,x^8$$
• S2: $K=6,P=0$, the mixed case with the set

  $$1,x,x^2,x^3,x^4,x^5,x^6,exp(\pm \omega x)$$
• S3: $K=4,P=1$, the mixed case with the set

  $$1,x,x^2,x^3,x^4,exp(\pm \omega x),xexp(\pm \omega x)$$
• S4: $K=2,P=2$, the mixed case with the set

  $$1,x,x^2,exp(\pm \omega x),xexp(\pm \omega x),x^{2}exp(\pm \omega x)$$
• S5: $K=0,P=3$, the mixed case with the set

  $$1,exp(\pm \omega x),xexp(\pm \omega x),x^{2}exp(\pm \omega x),x^{3}exp(\pm \omega x)$$

In order to get the corresponding coefficients of the method associated with each case, we implement step V of the algorithm by solving the algebraic system

$$L_k^{*}=0,0\le k\le K,G^{\left(p\right)\pm }(Z,\mathbf{a})=0,0\le p\le P$$

and the coefficients of the methods associated with each case are respectively obtained as follows:

S1 :: (K,P) = (8,−1)

$$\begin{array}{c}\hfill \left.\begin{array}{c}{a}_0=1\hfill \\ b_0=b_1=\frac{1}{2}\hfill \\ c_0=-c_1=\frac{3}{28}\hfill \\ d_0=d_1=\frac{1}{84}\hfill \\ e_0=-e_1=\frac{1}{1680}\hfill \end{array}\right\}\end{array}$$

(17)

Equation (17) gives the coefficients of the classical method associated with Equation (2).

S2 :: (K,P) = (6,0)

$$\begin{array}{c}\hfill \left.\begin{array}{c}{a}_0=1\hfill \\ b_0=b_1=\frac{1}{2}\hfill \\ c_0=-c_1=\frac{\left(z^4-120\right)sinh\left(\frac{z}{2}\right)-5z\left(z^2-12\right)cosh\left(\frac{z}{2}\right)}{10z^2\left(\left(z^2+12\right)sinh\left(\frac{z}{2}\right)-6zcosh\left(\frac{z}{2}\right)\right)}\hfill \\ d_0=d_1=\frac{\left(z^4-60z^2-720\right)sinh\left(\frac{z}{2}\right)+360zcosh\left(\frac{z}{2}\right)}{120z^2\left(\left(z^2+12\right)sinh\left(\frac{z}{2}\right)-6zcosh\left(\frac{z}{2}\right)\right)}\hfill \\ e_0=-e_1=\frac{z\left(z^2+60\right)cosh\left(\frac{z}{2}\right)-12\left(z^2+10\right)sinh\left(\frac{z}{2}\right)}{120z^2\left(\left(z^2+12\right)sinh\left(\frac{z}{2}\right)-6zcosh\left(\frac{z}{2}\right)\right)}\hfill \end{array}\right\}\end{array}$$

(18)

S3 :: (K,P) = (4,1)

$$\begin{array}{c}\hfill \left.\begin{array}{c}{a}_0=1\hfill \\ b_0=b_1=\frac{1}{2}\hfill \\ c_0=-c_1=\frac{z^4-12z^2+\left(z^2-36\right)zsinh\left(z\right)+96cosh\left(z\right)-96}{12z^2\left(z^2+zsinh\left(z\right)-4cosh\left(z\right)+4\right)}\hfill \\ d_0=d_1=\frac{-4\left(z^2+6\right)+\left(z^2+24\right)cosh\left(z\right)-9zsinh\left(z\right)}{6z^2\left(z^2+zsinh\left(z\right)-4cosh\left(z\right)+4\right)}\hfill \\ e_0=-e_1=-\frac{z^4+12z^2-\left(z^2+48\right)zsinh\left(z\right)+12\left(z^2+4\right)cosh\left(z\right)-48}{12z^4\left(z^2+zsinh\left(z\right)-4cosh\left(z\right)+4\right)}\hfill \end{array}\right\}\end{array}$$

(19)

S4 :: (K,P) = (2,2)

$$\begin{array}{c}\hfill \left.\begin{array}{c}{a}_0=1\hfill \\ b_0=b_1=\frac{1}{2}\hfill \\ c_0=-c_1=\frac{sinh\left(\frac{z}{2}\right)\left(2z^3+21z+24sinh\left(z\right)\right)-18z^{2}cosh\left(\frac{z}{2}\right)-3zsinh\left(\frac{3z}{2}\right)}{2z^2\left(\left(2z^2+1\right)cosh\left(\frac{z}{2}\right)-2zsinh\left(\frac{z}{2}\right)-cosh\left(\frac{3z}{2}\right)\right)}\hfill \\ d_0=d_1=-\frac{z^{3}cosh\left(\frac{z}{2}\right)+sinh\left(\frac{z}{2}\right)\left(3z^2+3zsinh\left(z\right)-16cosh\left(z\right)+16\right)}{z^3\left(\left(2z^2+1\right)cosh\left(\frac{z}{2}\right)-2zsinh\left(\frac{z}{2}\right)-cosh\left(\frac{3z}{2}\right)\right)}\hfill \\ e_0=-e_1=\frac{sinh\left(\frac{z}{2}\right)\left(-2z^3-21z+12sinh\left(z\right)\right)+6z^{2}cosh\left(\frac{z}{2}\right)-zsinh\left(\frac{3z}{2}\right)}{2z^4\left(\left(2z^2+1\right)cosh\left(\frac{z}{2}\right)-2zsinh\left(\frac{z}{2}\right)-cosh\left(\frac{3z}{2}\right)\right)}\hfill \end{array}\right\}\end{array}$$

(20)

S5 :: (K,P) = (0,3)

$$\begin{array}{c}\hfill \left.\begin{array}{c}{a}_0=1\hfill \\ b_0=b_1=-\frac{2\left(4z^3-6\left(2z^2+1\right)sinh\left(z\right)+2\left(z^2+3\right)zcosh\left(z\right)-6z+3sinh\left(2z\right)\right)}{z\left(\left(z^2+3\right)z^2-2\left(z^2+3\right)zsinh\left(z\right)+6z^{2}cosh\left(z\right)-3{sinh}^2\left(z\right)\right)}\hfill \\ c_0=-c_1=-\frac{2\left(z^4-3z^2(4cosh\left(z\right)+3)+9{sinh}^2\left(z\right)+12zsinh\left(z\right)\right)}{z^2\left(\left(z^2+3\right)z^2-2\left(z^2+3\right)zsinh\left(z\right)+6z^{2}cosh\left(z\right)-3{sinh}^2\left(z\right)\right)}\hfill \\ d_0=d_1=\frac{2\left(2z^3(cosh\left(z\right)+2)+6sinh\left(z\right)-3sinh\left(2z\right)-6z(cosh\left(z\right)-1)\right)}{z^3\left(\left(z^2+3\right)z^2-2\left(z^2+3\right)zsinh\left(z\right)+6z^{2}cosh\left(z\right)-3{sinh}^2\left(z\right)\right)}\hfill \\ e_0=-e_1=\frac{z^4+2z^{3}sinh\left(z\right)+z^2(3-6cosh\left(z\right))-3{sinh}^2\left(z\right)+6zsinh\left(z\right)}{z^4\left(\left(z^2+3\right)z^2-2\left(z^2+3\right)zsinh\left(z\right)+6z^{2}cosh\left(z\right)-3{sinh}^2\left(z\right)\right)}\hfill \end{array}\right\}\end{array}$$

(21)

As expected, the exponentially fitted variants reduce to the classical method as $z\to 0$.

3. Error Analysis :: Local Truncation Error (lte)

The leading term of the local truncation error (lte) for the exponentially-fitted method with respect to the basis functions

$$\begin{array}{c}\hfill \left\{1,x,\cdots ,x^K,exp(\pm \omega x),xexp(\pm \omega x),\cdots ,x^{P}exp(\pm \omega x)\right\}\end{array}$$

(22)

is of the form

$$\begin{array}{c}\hfill lt{e}^{EF}\left(t\right)={(-1)}^{P+1}{h}^M\frac{{\mathcal{L}}_{K+1}^{*}\left(\mathbf{a}\left(Z\right)\right)}{(K+1)!Z^{P+1}}{D}^{K+1}{(D^2-{\omega }^2)}^{P+1}y\left(x\right)\end{array}$$

(23)

with K, P and M satisfying the condition $K+2P=M-3$, [3].

For the five methods constructed in this work, the leading terms of the local truncation error are obtained as follows:

• S1 :: $(K,P)=(8,-1)$

  $$\begin{array}{ccc}\hfill lt{e}_{EF}\left(t\right)& =& \frac{h^9{y}^{\left(9\right)}\left(x\right)}{\mathrm{25,401,600}}\hfill \end{array}$$

  (24)
• S2 :: $(K,P)=(6,0)$

  $$\begin{array}{ccc}\hfill lt{e}_{EF}\left(t\right)& =& -h^9\frac{\left(20z\left(z^2+42\right)cosh\left(\frac{z}{2}\right)-\left(z^4+180z^2+1680\right)sinh\left(\frac{z}{2}\right)\right)}{\mathrm{100,800}{z}^4\left(\left(z^2+12\right)sinh\left(\frac{z}{2}\right)-6zcosh\left(\frac{z}{2}\right)\right)}\hfill \\ & & \times \left(u^{\left(9\right)}\left(t\right)-{\omega }^2{u}^{\left(7\right)}\left(t\right)\right)\hfill \end{array}$$

  (25)
• S3 :: $(K,P)=(4,1)$

  $$\begin{array}{ccc}\hfill lt{e}_{EF}\left(t\right)& =& \frac{h^9}{720z^8\left(z^2+zsinh\left(z\right)-4cosh\left(z\right)+4\right)}\hfill \\ \hfill & & \times \left(2880-240z^2+24z^4+z^6-\right.\hfill \\ \hfill & & 24\left(z^4+50z^2+120\right)cosh\left(z\right)+\hfill \\ \hfill & & \left.z\left(z^4+240z^2+2880\right)sinh\left(z\right)\right)\hfill \\ & & \times \left(y^{\left(9\right)}\left(x\right)-2{\omega }^2{y}^{\left(7\right)}\left(x\right)+{\omega }^4{y}^{\left(5\right)}\left(x\right)\right)\hfill \end{array}$$

  (26)
• S4 :: $(K,P)=(2,2)$

  $$\begin{array}{ccc}\hfill lt{e}_{EF}\left(t\right)& =& -\frac{h^9}{12z^9\left(\left(2z^2+1\right)cosh\left(\frac{z}{2}\right)-2zsinh\left(\frac{z}{2}\right)-cosh\left(\frac{3z}{2}\right)\right)}\hfill \\ \hfill & & \times \left(-z^3\left(2z^2+85\right)cosh\left(\frac{z}{2}\right)+z^{3}cosh\left(\frac{3z}{2}\right)+\right.\hfill \\ \hfill & & 2sinh\left(\frac{z}{2}\right)\left(192+99z^2+7z^4-192cosh\left(z\right)+\right.\hfill \\ \hfill & & \left.108zsinh\left(z\right))-18z^{2}sinh\left(\frac{3z}{2}\right)\right)\hfill \\ & & \times \left(y^{\left(9\right)}\left(x\right)-3{\omega }^2{y}^{\left(7\right)}\left(x\right)+3{\omega }^4{y}^{\left(5\right)}\left(x\right)-{\omega }^6{y}^{\left(3\right)}\left(x\right)\right)\hfill \end{array}$$

  (27)
• S5 :: $(K,P)=(0,3)$

  $$\begin{array}{ccc}\hfill lt{e}_{EF}\left(t\right)& =& \frac{h^9\left(\frac{4\left(4z^3-6\left(2z^2+1\right)sinh\left(z\right)+2\left(z^2+3\right)zcosh\left(z\right)-6z+3sinh\left(2z\right)\right)}{z\left(\left(z^2+3\right)z^2-2\left(z^2+3\right)zsinh\left(z\right)+6z^{2}cosh\left(z\right)-3{sinh}^2\left(z\right)\right)}+1\right)}{z^8}\hfill \\ & & \times \left(y^{\left(9\right)}\left(x\right)-4{\omega }^2{y}^{\left(7\right)}\left(x\right)+6{\omega }^4{y}^{\left(5\right)}\left(x\right)-4{\omega }^6{y}^{\left(3\right)}\left(x\right)+{\omega }^8{y}^{\prime }\left(x\right)\right)\hfill \end{array}$$

  (28)

4. Convergence and Stability Analysis

Theorem 1

(Dahlquist Theorem). The necessary and sufficient conditions for a linear multistep method to be convergent are that it be consistent and zero-stable [4].

Dahlquist Theorem (1) also holds true for exponetially-fitted-based algorithms but the concepts of consistency and stability have to be adapted since their coefficients are no longer constants.

Definition 1.

An exponentially-fitted method associated with the fitting space (22) is said to be of exponential order q, relative to the frequency ω if q is the maximum value of M such that the algebraic system $\left\{{\mathcal{L}}_m^{*}\left(\mathbf{a}\right)=0|m=0,\cdots ,M-1\right\}$ is compatible [3].

Definition 2.

A linear multistep method is said to be consistent if it has order $\mathcal{P}\ge 1$ [1,4].

Since the order of the constructed method, $M=9\ge 1$ for all the constructed schemes, the consistency requirement is satisfied. Hence, the constructed schemes are all consistent.

Definition 3.

The method Equation (2) is zero stable if no root of the first characteristic polynomial has modulus greater than one and if every root with modulus one is simple [1,2].

In order to establish the stability of (2), we apply Equation (2) to the test problems $y^{\prime }=\lambda y$ and obtain the stability function $R\left(q\right)$, of the class of methods as

$$\begin{array}{c}\hfill \frac{y_{n+1}}{y_n}=R\left(q\right)=\frac{1+b_{0}q+c_0{q}^2+d_0{q}^3+e_0{q}^4}{1-b_{1}q-c_1{q}^2-d_1{q}^3-e_1{q}^4},\mathrm{with}q=\lambda h.\end{array}$$

(29)

Definition 4.

A region of absolute stability is a region in the complex plane, throughout which $\left|R\left(q\right)\right|<1$. Any closed curve defined by $\left|R\left(q\right)\right|=1$ is an absolute stability boundary. In addition, any interval $(\alpha ,\beta )$ of the real line is said to be the interval of absolute stability if the method is stable for all $q\in (\alpha ,\beta )$ [1,4].

The absolute stability regions for all the methods constructed in this work are the same and given in Figure 1.

From Figure 1, it can be seen that stability region of the methods contains the entire left half plane, hence they are all A-stable and have their absolute stability interval as $\left(-\infty ,0\right]$.

5. Numerical Results

In this section, we considered two test problems. The constructed methods are implemented on these test problems and the obtained results were compared with those of the classical eighth-order Runge-Kutta (RK-8) method.

5.1. Problem 1

The first test problem considered in this work is the initial value problem given as

$$y^{″}-y=0.001cosx,y\left(0\right)=1,y^{\prime }=0$$

with exact solution

$$y\left(x\right)=cosx+0.0005xsinx$$

This problem was studied in [24,25]. Using different stepsizes, we implement the constructed methods on this problem and present the maximum absolute errors in

Table 1. Maximum absolute error for constructed methods on Section 5.1 with step-size $h=2^{-i}\pi ,i=0,1,2,3$.

i	RK-8	$(\mathbf{K},\mathbf{P})=(8,-1)$	$(\mathbf{K},\mathbf{P})=(6,0)$	$(\mathbf{K},\mathbf{P})=(4,2)$	$(\mathbf{K},\mathbf{P})=(2,1)$	$(\mathbf{K},\mathbf{P})=(0,3)$
0	8.256847$\times {10}^{-1}$	1.701678$\times {10}^{-5}$	0.000000	0.000000	0.000000	2.220446$\times {10}^{-16}$
1	5.921881$\times {10}^{-3}$	1.485065$\times {10}^{-5}$	1.516774$\times {10}^{-8}$	1.160079$\times {10}^{-13}$	4.170189$\times {10}^{-14}$	9.992007$\times {10}^{-16}$
2	9.221257$\times {10}^{-5}$	6.125912$\times {10}^{-8}$	6.400131$\times {10}^{-11}$	1.488434$\times {10}^{-11}$	6.878617$\times {10}^{-12}$	7.895672$\times {10}^{-13}$
3	1.154806$\times {10}^{-6}$	2.425554$\times {10}^{-10}$	3.548604$\times {10}^{-11}$	1.3173$\times {10}^{-9}$	3.317455$\times {10}^{-10}$	7.396323$\times {10}^{-11}$

.

For this problem, the exponentially fitted methods gave better results compared with their classical counterpart and the Runge–Kutta method as seen from Table 1.

5.2. Problem 2

The inhomogeneous equation

$$y^{″}+100y=99sinx,y\left(0\right)=1,y^{\prime }=11.$$

with exact solution

$$y\left(x\right)=sin\left(x\right)+sin\left(10x\right)+cos\left(10x\right)$$

is considered as the second test case. This problem has also been studied by [24,26]. The constructed methods were implemented on it with different stepsizes and the results obtained were also compared with those of the Runge-Kutta method. The table of maximum absolute errors is given in

Table 2. Maximum absolute error for constructed methods on Section 5.2 with step-size $h=2^{-i}\pi ,i=0,1,2,3,4,5$.

i	RK-8	$(\mathbf{K},\mathbf{P})=(8,-1)$	$(\mathbf{K},\mathbf{P})=(6,0)$	$(\mathbf{K},\mathbf{P})=(4,2)$	$(\mathbf{K},\mathbf{P})=(2,1)$	$(\mathbf{K},\mathbf{P})=(0,3)$
0	2.385689$\times {10}^{27}$	2.385402	6.994405$\times {10}^{-15}$	5.77316$\times {10}^{-15}$	1.210143$\times {10}^{-14}$	3.719247$\times {10}^{-14}$
1	2.779917$\times {10}^{37}$	2.410011	1.507892$\times {10}^{-7}$	1.467305$\times {10}^{-4}$	4.454081$\times {10}^{-3}$	7.805366
2	1.036922$\times {10}^{39}$	2.413778	2.714209$\times {10}^{-7}$	4.764676$\times {10}^{-5}$	6.398771$\times {10}^{-3}$	1.120104
3	1.466513$\times {10}^{12}$	2.377447$\times {10}^{-1}$	4.035951$\times {10}^{-10}$	4.450952$\times {10}^{-8}$	4.873151$\times {10}^{-6}$	5.310473$\times {10}^{-4}$
4	4.728532$\times {10}^{-1}$	1.336235$\times {10}^{-3}$	1.392109$\times {10}^{-12}$	1.123486$\times {10}^{-10}$	1.137188$\times {10}^{-8}$	1.153992$\times {10}^{-6}$
5	5.249947$\times {10}^{-3}$	5.711746$\times {10}^{-6}$	1.162043$\times {10}^{-11}$	6.347695$\times {10}^{-11}$	3.808664$\times {10}^{-11}$	4.012365$\times {10}^{-9}$

.

Again, the exponentially fitted methods gave better results compared with their classical counterpart and the Runge-Kutta method as seen from Table 2.

6. Conclusions

An exponentially-fitted one-step fourth-derivative Obrechkoff method for oscillatory problems was constructed. The new methods are self-starting and of algebraic order eight. The stability and convergence properties of the constructed method were analysed and we showed that the new methods are A-stable. The results obtained from the numerical examples show that the new methods are suitable for solving periodic/oscillatory problems.