# A Trigonometrically Fitted Block Method for Solving Oscillatory Second-Order Initial Value Problems and Hamiltonian Systems

Authors: F. F. Ngwane; S. N. Jator
Journal: International Journal of Differential Equations
 (2017)

---

## Abstract

In this paper, we present a block hybrid trigonometrically fitted Runge-Kutta-Nyström method (BHTRKNM), whose coefficients are functions of the frequency and the step-size for directly solving general second-order initial value problems (IVPs), including Hamiltonian systems such as the energy conserving equations and systems arising from the semidiscretization of partial differential equations (PDEs). Four discrete hybrid formulas used to formulate the BHTRKNM are provided by a continuous one-step hybrid trigonometrically fitted method with an off-grid point. We implement BHTRKNM in a block-by-block fashion; in this way, the method does not suffer from the disadvantages of requiring starting values and predictors which are inherent in predictor-corrector methods. The stability property of the BHTRKNM is discussed and the performance of the method is demonstrated on some numerical examples to show accuracy and efficiency advantages.


---

## Body

## 1. Introduction

In what follows, we consider the numerical solution of the general second-order IVPs of the form(1)y′′=fx,y,y′,yx0=y0,y′x0=y0′,x∈x0,xN,where f:R×R2m→R2m, N>0 is an integer, and m is the dimension of the system. Problems of form (1) frequently arise in several areas of science and engineering such as classical mechanics, celestial mechanics, quantum mechanics, control theory, circuit theory, astrophysics, and biological sciences. Equation (1) is traditionally solved by reducing it into a system of first-order IVPs of double dimension and then solved using the various methods that are available for solving systems of first-order IVPs (see Lambert [1, 2], Hairer and Wanner in [3], Hairer [4], and Brugnano and Trigiante [5, 6]).Nevertheless, there are numerous methods for directly solving the special second-order IVPs in which the first derivative does not appear explicitly and it has been shown that these methods have the advantages of requiring less storage space and fewer number of function evaluations (see Hairer [4], Hairer et al. [7], Simos [8], Lambert and Watson, and [9], Twizell and Khaliq [10]). Fewer methods have been proposed for directly solving second-order IVPs in which the first derivative appears explicitly (see Vigo-Aguiar and Ramos [11], Awoyemi [12], Chawla and Sharma [13], Mahmoud and Osman [14], Franco [15], and Jator [16, 17]). It is also the case that some of these IVPs possess solutions with special properties that may be known in advance and take advantage of when designing numerical methods. In this light, several methods have been presented in the literature which take advantage of the special properties of the solution that may be known in advance (see Coleman and Duxbury [18], Coleman and Ixaru [19], Simos [20], Vanden Berghe et al. [21], Vigo-Aguiar and Ramos [11], Fang et al. [22], Nguyen et al. [23], Ramos and Vigo-Aguiar [24], Franco and Gómez [25], and Ozawa [26]). However, most of these methods are restricted to solving special second-order IVPs in a predictor-corrector mode.Our objective is to present a BHTRKNM that is implemented in a block-by-block fashion; in this way, the method does not suffer from the disadvantages of requiring starting values and predictors which are inherent in predictor-corrector methods (see Jator et al. [27], Jator [16], and Ngwane and Jator [28]). We note that multiderivative trigonometrically fitted block methods for y′′=f(x,y,y′) have been proposed in Jator [29] and Jator [16]. However, the BHTRKNM proposed in this paper avoids the computation of higher order derivatives which have the potential to increase computational cost, especially, when applied to nonlinear systems. In this paper, we propose a BHTRKNM which is of order 3 and its application is extended to solving oscillatory systems, PDEs, and Hamiltonian systems including the energy conserving equation.The organization of this article is as follows. In Section2, we derive the BHTRKNM for solving (1). The analysis and implementation of the BHTRKNM are discussed in Section 3. Numerical examples are given in Section 4 to show the accuracy and efficiency of the BHTRKNM. Finally, the conclusion of the paper is given in Section 5.


## 2. Development of the BHTRKNM

In order to numerical integrate (1) we define the BHTRKNM as consisting of the following four discrete formulas:(2)yn+1=yn+hyn′+h2∑j=01βjfn+j+βn+vfn+v,yn+v=yn+hyn′+h2∑j=01γjfn+j+γn+vfn+v,hyn+1′=hyn′+h2∑j=01βj′fn+j+βn+v′fn+v,hyn+v′=hyn′+h2∑j=01γj′fn+j+γn+v′fn+v,where βj, βj′, γj, and γj′ are coefficients that depend on the step-length h and frequency w. In general, the frequency w is chosen near the exact frequency of the true solution (see [30]). The coefficients of the method are chosen so that the method integrates the IVP (1) exactly where the solutions are members of the linear space 1,x,x2,sin⁡wx,cos⁡wx.The main method has the form(3)yn+1=α0yn+δ0hyn′+h2∑j=01βjfn+j+βn+vfn+v,where α0, δ0, and β0, βv, and β1 are to be determined coefficient functions of the frequency and step-size. In order to derive the main method and additional methods we initially seek a continuous local approximation Π(x) on the interval [xn,xn+1] of the form(4)Πx=α0xyn+δ0xhyn′+h2∑j=01βjxfn+j+βn+vxfn+v,where α0(x), δ0(x), and βj(x), j=0,v,1, are continuous coefficients. The first derivative of (4) is given by(5)Π′x=ddxΠx.We assume that yn+j=Π(xn+j) is the numerical approximation to the analytical solution y(xn+j), yn+j′=Π′(xn+j) is the numerical approximation to y′(xn+j), and fn+j=Π′′(xn+j) is an approximation to y′(xn+j), j=0,v,1.The following theorem shows how the continuous method (4) is constructed. This is done by requiring that on the interval from xn to xn+1=xn+h the exact solution is locally approximated by function (4) with (5) obtained as a consequence.Theorem 1.
LetFi(x)=xi, i=0,1,2, F3(x)=sin⁡wx, and F4x=cos⁡wx be basis functions and let V=(yn,yn′,fn,fn+v,fn+1)T be a vector, where T is the transpose. Define the matrix G by(6)G=F0xn⋯F4xnF0′xn⋯F4′xnF0′′xn⋯F4′′xnF0′′xn+v⋯F4′′xn+vF0′′xn+1⋯F4′′xn+1and Gi is obtained by replacing the ith column of G by the vector V. Let the following conditions be satisfied:(7)Πxn=yn,Π′xn=yn′,Π′′xn+j=fn+j,j=0,v,1;then the continuous representations (4) and (5) are equivalent to the following:(8)Πx=∑i=04det⁡Gidet⁡GFix,(9)Π′x=ddx∑i=04det⁡Gidet⁡GFix.Proof.
To prove this theorem, we use the approach given in Jator [17] with appropriate notational modification. We start by requiring that the method (4) be defined by the assumed basis functions(10)α0x=∑i=04αi+1,0Fix,hδ0x=∑i=04hδi+1,0Fix,h2βjx=∑i=04h2βi+1,jFix,where αi+1,0, hδi+1,0, and h2βi+1,j are coefficients to be determined. Substituting (10) into (4) we get (11)Πx=∑i=04αi+1,0Fixyn+∑i=04hδi+1,0Fixyn′+∑j=01∑i=04h2βi+1,jFixfn+jwhich is simplified to (12)Πx=∑i=04αi+1,0Fixyn+hδi+1,0Fixyn′+∑j=01h2βi+1,jFixfn+jand expressed as(13)Πx=∑i=04liFix,where (14)li=αi+1,0yn+hδi+1,0yn′+∑j=01h2βi+1,jfn+j.By imposing conditions (7) on (13), we obtain a system of five equations which can be expressed as (15)GL=V,where L=(l0,l1,…,l4)T is a vector whose coefficients are determined via Cramer’s rule as (16)li=det⁡Gidet⁡G,i=0,1,…,4,where Gi is obtained by replacing the ith column of G by V. In order to obtain the continuous approximation, we use the elements of L to rewrite (13) as(17)Πx=∑i=04det⁡Gidet⁡GFix,whose first derivative is given by (18)Π′x=ddx∑i=04det⁡Gidet⁡GFix.Remark 2.
We note that, in the derivation of the BHTRKNM, the basis functionsFi(x)=xi, i=0,1,2, F3(x)=sin⁡wx, and F4(x)=cos⁡wx are chosen because they are simple to analyze. Nevertheless, other possible bases are possible (see Nguyen et al. [23]).


### 2.1. Specification of the Method

The continuous methods (8) and (9) which are equivalent to forms (4) and (5) are used to generate two discrete methods and two additional methods. The discrete and additional methods are then applied as a BHTRKNM for solving (1). We choose v=1/2 and evaluating (8) at x=xn+v and x=xn+1, respectively, gives the two discrete methods yn+v=Π(xn+vh) and yn+1=Π(xn+h) which takes the form of the main method. Evaluating (9) at x=xn+v and x=xn+1, respectively, gives the additional methods yn+v′=Π′(xn+vh) and yn+1′=Π′xn+h. The coefficients and their corresponding Taylor series equivalence of yn+v, yn+1, hyn+v′, and hyn+1′ are, respectively, given as follows:(19)αv,0=1,δv,0=12,βv,0=csc⁡u/4csc⁡u/28+u2cos⁡u/4-42cos⁡3u/4+usin⁡3u/416u2=796+7u27680+71u43870720+53u6123863040+23u82179989504+⋯,βv,v=-csc⁡u/428+-8+u2cos⁡u/2-4usin⁡u/216u2=116-u22304-u4276480-u634406400-u84459069440+⋯,βv,1=csc⁡u/42u+4cot⁡u/2-4csc⁡u/232u=-196-11u223040-19u41290240-247u6619315200-1013u898099527680+⋯,(20)α1,0=1,δ1,0=1,β1,0=csc⁡u/4csc⁡u/22+u2cos⁡u/4-2cos⁡3u/4+usin⁡3u/44u2=16+u2480+19u4483840+17u619353600+29u81362493440+⋯,β1,v=-csc⁡u/42csc⁡u/2-2+2cos⁡u+usin⁡u8u=13-u2720-u480640-u69676800-u81226244096+⋯,β1,1=csc⁡u/42-4+u2+4cos⁡u/2+2ucot⁡u/2-2ucsc⁡u/28u2=-u21440-13u4483840-u61290240-251u812262440960+⋯,(21)αv,0=0,δv,0=1,γv,0′=-cot⁡u/2u+18csc⁡u42=524+19u25760+23u4322560+263u6154828800+1033u824524881920+⋯,γv,v′=-csc⁡u/42csc⁡u/2-2+2cos⁡u+usin⁡u8u=13-u2720-u480640-u69676800-u81226244096+⋯,γv,1′=csc⁡u/42u+4cot⁡u/2-4csc⁡u/28u=-124-11u25760-19u4322560-247u6154828800-1013u824524881920+⋯,(22)α1,0=0,δ1,0=1,β1,0′=csc⁡u/42u-2sin⁡u/24u=16+u2720+u480640+u69676800+u81226244096+⋯,β1,v′=csc⁡u/42u-2sin⁡u/24u=23-u2360-u440320-u64838400-u8613122048+⋯,β1,1′=csc⁡u/42-4+u2+4cos⁡u/2+2ucot⁡u/2-2ucsc⁡u/28u2=16+u2720+u480640+u69676800+u81226244096+⋯.Remark 3.
We note that the Taylor series expansions in (19) through (22) must be used when u→0 because the corresponding trigonometric coefficients given in these equations are vulnerable to heavy cancelations (see [8]).


### 2.2. Block Form

BHTRKNM is formulated from the four discrete hybrid formulas stated in (2) which are provided by the continuous one-step hybrid trigonometrically fitted method with one off-grid point given by (4) and its first derivative (5). We define the following vectors: (23)Yμ+1=yn+v,yn+1,hyn+v′,hyn+1′T,Yμ=yn-v,yn,hyn-v′,hyn′T,Fμ+1=fn+v,fn+1,hfn+v′,hfn+1′T,Fμ=fn-v,fn,hfn-v′,hfn′T,where μ=0,…,N, n=0,…,N. The methods in (2) specified by the coefficients (19)–(22) are combined to give the BHTRKNM, which is expressed as(24)A1Yμ+1=A0Yμ+h2B0Fμ+B1Fμ+1,where A0, A1, B0, and B1 are matrices of dimension four whose elements characterize the method and are given by the coefficients of (2).


## 2.1. Specification of the Method

The continuous methods (8) and (9) which are equivalent to forms (4) and (5) are used to generate two discrete methods and two additional methods. The discrete and additional methods are then applied as a BHTRKNM for solving (1). We choose v=1/2 and evaluating (8) at x=xn+v and x=xn+1, respectively, gives the two discrete methods yn+v=Π(xn+vh) and yn+1=Π(xn+h) which takes the form of the main method. Evaluating (9) at x=xn+v and x=xn+1, respectively, gives the additional methods yn+v′=Π′(xn+vh) and yn+1′=Π′xn+h. The coefficients and their corresponding Taylor series equivalence of yn+v, yn+1, hyn+v′, and hyn+1′ are, respectively, given as follows:(19)αv,0=1,δv,0=12,βv,0=csc⁡u/4csc⁡u/28+u2cos⁡u/4-42cos⁡3u/4+usin⁡3u/416u2=796+7u27680+71u43870720+53u6123863040+23u82179989504+⋯,βv,v=-csc⁡u/428+-8+u2cos⁡u/2-4usin⁡u/216u2=116-u22304-u4276480-u634406400-u84459069440+⋯,βv,1=csc⁡u/42u+4cot⁡u/2-4csc⁡u/232u=-196-11u223040-19u41290240-247u6619315200-1013u898099527680+⋯,(20)α1,0=1,δ1,0=1,β1,0=csc⁡u/4csc⁡u/22+u2cos⁡u/4-2cos⁡3u/4+usin⁡3u/44u2=16+u2480+19u4483840+17u619353600+29u81362493440+⋯,β1,v=-csc⁡u/42csc⁡u/2-2+2cos⁡u+usin⁡u8u=13-u2720-u480640-u69676800-u81226244096+⋯,β1,1=csc⁡u/42-4+u2+4cos⁡u/2+2ucot⁡u/2-2ucsc⁡u/28u2=-u21440-13u4483840-u61290240-251u812262440960+⋯,(21)αv,0=0,δv,0=1,γv,0′=-cot⁡u/2u+18csc⁡u42=524+19u25760+23u4322560+263u6154828800+1033u824524881920+⋯,γv,v′=-csc⁡u/42csc⁡u/2-2+2cos⁡u+usin⁡u8u=13-u2720-u480640-u69676800-u81226244096+⋯,γv,1′=csc⁡u/42u+4cot⁡u/2-4csc⁡u/28u=-124-11u25760-19u4322560-247u6154828800-1013u824524881920+⋯,(22)α1,0=0,δ1,0=1,β1,0′=csc⁡u/42u-2sin⁡u/24u=16+u2720+u480640+u69676800+u81226244096+⋯,β1,v′=csc⁡u/42u-2sin⁡u/24u=23-u2360-u440320-u64838400-u8613122048+⋯,β1,1′=csc⁡u/42-4+u2+4cos⁡u/2+2ucot⁡u/2-2ucsc⁡u/28u2=16+u2720+u480640+u69676800+u81226244096+⋯.Remark 3.
We note that the Taylor series expansions in (19) through (22) must be used when u→0 because the corresponding trigonometric coefficients given in these equations are vulnerable to heavy cancelations (see [8]).


## 2.2. Block Form

BHTRKNM is formulated from the four discrete hybrid formulas stated in (2) which are provided by the continuous one-step hybrid trigonometrically fitted method with one off-grid point given by (4) and its first derivative (5). We define the following vectors: (23)Yμ+1=yn+v,yn+1,hyn+v′,hyn+1′T,Yμ=yn-v,yn,hyn-v′,hyn′T,Fμ+1=fn+v,fn+1,hfn+v′,hfn+1′T,Fμ=fn-v,fn,hfn-v′,hfn′T,where μ=0,…,N, n=0,…,N. The methods in (2) specified by the coefficients (19)–(22) are combined to give the BHTRKNM, which is expressed as(24)A1Yμ+1=A0Yμ+h2B0Fμ+B1Fμ+1,where A0, A1, B0, and B1 are matrices of dimension four whose elements characterize the method and are given by the coefficients of (2).


## 3. Error Analysis and Stability

### 3.1. Local Truncation Error (LTE)

We define the local truncation error of (24) as(25)LZx;h=Zμ+1-AZμ+h2BF¯μ+h2CF¯μ+1,where (26)Zμ+1=yxn+v,yxn+1,hy′xn+v,hy′xn+1T,Zμ=yxn-v,yxn,hy′xn-v,hy′xnT,F¯μ+1=fxn+v,yn+v,fxn+1,yn+1,hf′xn+v,yn+v,hf′xn+1,yn+1T,F¯μ=fxn-v,yn-v,fxn,yn,hf′xn-v,yn-v,hf′xn,ynT,and L[Z(x);h]=[L1[z(x);h],L2[z(x);h],…,L4[z(x);h]]T is linear different operator.Suppose thatZ(x) is sufficiently differentiable. Then, a Taylor series expansion of the terms in (25) about the point x gives the following expression for local truncation error:(27)LZx;h=C0Zx+C1hZ′x+⋯+CqhqZqx+⋯,where Ci, i=0,1,…, are constant coefficients (see [17]).Definition 4.
The block method (24) has algebraic order at least p≥1 provided there exists a constant Cp+2≠0 such that the local truncation error Eμ satisfies Eμ=Cp+2hp+2+O(hp+3), where · is the maximum norm.Remark 5.
(i) The local truncation error constants(C¯p+2) of (yn+v,yn+1,hyn+v′,hyn+1′)T of the block method (24) are given, respectively, by C¯5=1/1440,1/720,1/384,0T, where C¯0=C¯1=C¯2=C¯3=C¯4=0.
(ii) From the local truncation error constant computation, it follows that the method (24) has order p at least three.


### 3.2. Stability

The linear-stability of the BHTRKNM is discussed by applying the method to the test equationy′′=-λ2y, where λ is a real constant (see [18]). Letting Υ=λh, it is easily shown as in [19] that the application of (24) to the test equation yields(28)Yμ+1=MΥ2;uYμ,MΥ2;u≔A1-Υ2B1-1A0+Υ2B0,where the matrix M(Υ2;u) is the amplification matrix which determines the stability of the method. In the spirit of [22], the spectral radius of ρ(M(Υ2;u)) can be obtained from the characteristics equation(29)ρ2-2ΓΥ2;uρ+ΘΥ2;u=0,where u=ωh, ΓΥ2;u=traceM(Υ2;u), and Θ(Υ2;u)=det⁡M(Υ2;u) are rational functions. We let q=λh in the following definition.Definition 6.
A region of stability is a region in theq-u plane, throughout which ρ(M(Υ2;u))≤1 and any closed curve given by ρ(M(Υ2;u))=1 defines the stability boundary of the method (see [22]). We note that the plot for the stability region of the BHTRKNM method is given in Figure 1.Figure 1
The stability region plotted in the (q,u)-plane.Remark 7.
It is observed that, in theq-u plane, the BHTRKNM is stable for q∈[0,47.96] and u∈[-π,π] (see Figure 1).


### 3.3. Implementation

The main method and the additional methods specified by (19)–(22) are combined to form the block method BHTRKNM (24), which is used to solve (1) without requiring starting values and predictors. BHTRKNM is implemented in a block-by-block fashion using a Mathematica 10.0 code, enhanced by the feature NSolve[] for linear problems while nonlinear problems were solved by Newton’s method enhanced by the feature FindRoot[] (see Keiper and Gear [33]). Mathematica can symbolically compute derivatives and so the entries of the Jacobian matrix which involve partial derivatives are automatically generated. In what follows, we summarize how BHTRKNM is applied.Step 1.
ChooseN, h=(b-a)/N, and the number of blocks Γ=N. Using (24), n=0, and μ=0, the values of (y1/2,y1,y1/2′,y1′)T are simultaneously obtained over the subinterval [x0,x1], as y0 and y0′ are known from the IVP (1).Step 2.
Forn=1 and μ=1, the values of (y3/2,y2,y3/2′,y2′)T are simultaneously obtained over the subinterval [x1,x2], as y1 and y1′ are known from the previous block.Step 3.
The process is continued forn=2,…,N-1 and μ=2,…,Γ-1 to obtain the numerical solution to (1) on the subintervals x0,x1,[x1,x2],…,[xN-1,xN].In order to illustrate the efficiency of our method, we solved a variety of problems including oscillatory systems, PDEs such as the Telegraph equation, and Hamiltonian systems. The following methods are selected for comparison:(i)
BHTRKNM given in this paper.(ii)
ARKN: adapted Runge-Kutta-Nyström method in [34] which has order five.(iii)
(DS3.12): difference scheme (3.12) in [32].(iv)
ESDIRK: explicit singly diagonally implicit Runge-Kutta method in [26].(v)
FESDIRK: functionally fitted ESDIRK in [26].(vi)
EFRK: exponentially fitted Runge-Kutta method (Method (b)) in Simos [8].(vii)
N4: fourth-order standard Runge-Kutta-Nyström method in [35].


## 3.1. Local Truncation Error (LTE)

We define the local truncation error of (24) as(25)LZx;h=Zμ+1-AZμ+h2BF¯μ+h2CF¯μ+1,where (26)Zμ+1=yxn+v,yxn+1,hy′xn+v,hy′xn+1T,Zμ=yxn-v,yxn,hy′xn-v,hy′xnT,F¯μ+1=fxn+v,yn+v,fxn+1,yn+1,hf′xn+v,yn+v,hf′xn+1,yn+1T,F¯μ=fxn-v,yn-v,fxn,yn,hf′xn-v,yn-v,hf′xn,ynT,and L[Z(x);h]=[L1[z(x);h],L2[z(x);h],…,L4[z(x);h]]T is linear different operator.Suppose thatZ(x) is sufficiently differentiable. Then, a Taylor series expansion of the terms in (25) about the point x gives the following expression for local truncation error:(27)LZx;h=C0Zx+C1hZ′x+⋯+CqhqZqx+⋯,where Ci, i=0,1,…, are constant coefficients (see [17]).Definition 4.
The block method (24) has algebraic order at least p≥1 provided there exists a constant Cp+2≠0 such that the local truncation error Eμ satisfies Eμ=Cp+2hp+2+O(hp+3), where · is the maximum norm.Remark 5.
(i) The local truncation error constants(C¯p+2) of (yn+v,yn+1,hyn+v′,hyn+1′)T of the block method (24) are given, respectively, by C¯5=1/1440,1/720,1/384,0T, where C¯0=C¯1=C¯2=C¯3=C¯4=0.
(ii) From the local truncation error constant computation, it follows that the method (24) has order p at least three.


## 3.2. Stability

The linear-stability of the BHTRKNM is discussed by applying the method to the test equationy′′=-λ2y, where λ is a real constant (see [18]). Letting Υ=λh, it is easily shown as in [19] that the application of (24) to the test equation yields(28)Yμ+1=MΥ2;uYμ,MΥ2;u≔A1-Υ2B1-1A0+Υ2B0,where the matrix M(Υ2;u) is the amplification matrix which determines the stability of the method. In the spirit of [22], the spectral radius of ρ(M(Υ2;u)) can be obtained from the characteristics equation(29)ρ2-2ΓΥ2;uρ+ΘΥ2;u=0,where u=ωh, ΓΥ2;u=traceM(Υ2;u), and Θ(Υ2;u)=det⁡M(Υ2;u) are rational functions. We let q=λh in the following definition.Definition 6.
A region of stability is a region in theq-u plane, throughout which ρ(M(Υ2;u))≤1 and any closed curve given by ρ(M(Υ2;u))=1 defines the stability boundary of the method (see [22]). We note that the plot for the stability region of the BHTRKNM method is given in Figure 1.Figure 1
The stability region plotted in the (q,u)-plane.Remark 7.
It is observed that, in theq-u plane, the BHTRKNM is stable for q∈[0,47.96] and u∈[-π,π] (see Figure 1).


## 3.3. Implementation

The main method and the additional methods specified by (19)–(22) are combined to form the block method BHTRKNM (24), which is used to solve (1) without requiring starting values and predictors. BHTRKNM is implemented in a block-by-block fashion using a Mathematica 10.0 code, enhanced by the feature NSolve[] for linear problems while nonlinear problems were solved by Newton’s method enhanced by the feature FindRoot[] (see Keiper and Gear [33]). Mathematica can symbolically compute derivatives and so the entries of the Jacobian matrix which involve partial derivatives are automatically generated. In what follows, we summarize how BHTRKNM is applied.Step 1.
ChooseN, h=(b-a)/N, and the number of blocks Γ=N. Using (24), n=0, and μ=0, the values of (y1/2,y1,y1/2′,y1′)T are simultaneously obtained over the subinterval [x0,x1], as y0 and y0′ are known from the IVP (1).Step 2.
Forn=1 and μ=1, the values of (y3/2,y2,y3/2′,y2′)T are simultaneously obtained over the subinterval [x1,x2], as y1 and y1′ are known from the previous block.Step 3.
The process is continued forn=2,…,N-1 and μ=2,…,Γ-1 to obtain the numerical solution to (1) on the subintervals x0,x1,[x1,x2],…,[xN-1,xN].In order to illustrate the efficiency of our method, we solved a variety of problems including oscillatory systems, PDEs such as the Telegraph equation, and Hamiltonian systems. The following methods are selected for comparison:(i)
BHTRKNM given in this paper.(ii)
ARKN: adapted Runge-Kutta-Nyström method in [34] which has order five.(iii)
(DS3.12): difference scheme (3.12) in [32].(iv)
ESDIRK: explicit singly diagonally implicit Runge-Kutta method in [26].(v)
FESDIRK: functionally fitted ESDIRK in [26].(vi)
EFRK: exponentially fitted Runge-Kutta method (Method (b)) in Simos [8].(vii)
N4: fourth-order standard Runge-Kutta-Nyström method in [35].


## 4. Numerical Examples

In this section, numerical experiments are performed using a code in Mathematica 10.0 to illustrate the accuracy and efficiency of the method.Example 1.
We consider the following inhomogeneous IVP by Simos [8].(30)y′′=-100y+99sin⁡x,y0=1,y′0=11,x∈0,1000,where the analytical solution is given by(31)Exact:yx=cos⁡10x+sin⁡10x+sin⁡x.This example was solved using the order 3 BHTRKNM and the endpoint errors (Err=|y(xN)-yN|) obtained were compared to the order 4 exponentially fitted method given in Simos [8]. In Table 1 it is shown that BHTRKNM is more efficient than the method in Simos [8]. We also compare the computational efficiency of the two methods in Figure 2 by considering the FNEs (number of function evaluations) over N integration steps for each method. This example illustrates that the BHTRKNM performs better.Table 1
Results, withω=10, for Example 1.
Our method Simos [8] N Err N Err 1000 2.14 × 1 0 - 3 1000 1.4×10-1 2000 5.98 × 1 0 - 5 2000 3.5×10-2 4000 2.06 × 1 0 - 5 4000 1.1×10-3 8000 1.26 × 1 0 - 6 8000 8.4×10-5 16000 7.79 × 1 0 - 8 16000 5.5×10-6 32000 4.67 × 1 0 - 9 32000 3.5×10-7Figure 2
Efficiency curve for Example1.Example 2.
We consider the nonlinear Duffing equation which was also solved by Simos [8] and Ixaru and Vanden Berghe [31]:(32)y′′+y+y3=Bcos⁡Ωx,y0=C0,y′0=0.The analytical solution is given by(33)Exact:yx=C1cos⁡Ωx+C2cos⁡3Ωx+C3cos⁡5Ωx+C4cos⁡7Ωx,where Ω=1.01, B=0.002, C0=0.200426728069, C1=0.200179477536, C2=0.246946143×10-3, C3=0.304016×10-6, and C4=0.374×10-9. We choose ω=1.01.We compare the endpoint global errors for our method with those of Simos [8] and Ixaru and Vanden Berghe [31]. We see from Table 2 that the results produced by our method are competitive to those given in Simos [8] and Ixaru and Vanden Berghe [31]. Hence our method is more accurate and efficient as demonstrated in Figure 3.Table 2
Results, withω=1.01, for Example 2.
Our method Simos [8] Ixaru and VandenBerghe [31] N Err N Err N Err 300 7.52×10-5 300 1.7×10-3 300 1.1 × 1 0 - 3 600 2.47×10-6 600 1.9×10-4 600 5.4 × 1 0 - 5 1200 1.34×10-7 1200 1.4×10-5 1200 1.9 × 1 0 - 6 2400 8.11×10-9 2400 8.7×10-7 2400 6.2 × 1 0 - 8Figure 3
Efficiency curves for Example2.Example 3.
We consider the following two-body problem which was solved by Ozawa [26] on [0,50π]:(34)y1′′=-y1r3,y2′′=-y2r3,r=y12+y22,y10=1-e,y1′0=0,y20=0,y2′0=1+e1-e,where e(0≤e<1) is an eccentricity. The exact solution of this problem is(35)Exact:y1x=cos⁡k-e,y2x=1-e2sin⁡k,where k is the solution of Kepler’s equation k=x+esin⁡(k). We choose ω=1.We show in Table3 that the results obtained using the BHTRKNM method are more accurate than the explicit singly diagonally implicit Runge-Kutta (ESDIRK) and the functionally fitted ESDIRK (FESDIRK) methods given in Ozawa [26]. In Figure 4, we also illustrate the efficiency advantage of the BHTRKNM method over those in Ozawa [26].Table 3
Steps and absolute errors, withω=1, for Example 3  [0,50π].
Our method FESDIRK4(3) [26] ESDIRK4(3) [26] Steps Err Steps Err Steps Err 200 4.42 × 1 0 - 4 381 1.40 × 1 0 - 3 884 9.36×10-3 300 3.2 × 1 0 - 5 680 1.69 × 1 0 - 4 1573 6.20×10-4 400 5.39 × 1 0 - 8 1207 1.85 × 1 0 - 5 2796 4.42×10-5 600 4.25 × 1 0 - 7 2144 1.94 × 1 0 - 6 4970 3.41×10-6 1000 1.06 × 1 0 - 8 3806 1.99 × 1 0 - 7 8833 2.85×10-7 1200 1.76 × 1 0 - 9 6762 2.02 × 1 0 - 8 15706 2.53×10-8Figure 4
Efficiency curves for Example3.Example 4.
We consider the stiff second-order IVP (see [16] and references herein) (36)y1′′=ϵ-2y1+2ϵ-2y2,y2′′=1-ϵy1+1-2ϵy2,y10=2,y1′0=0,y20=-1,y2′0=0,ϵ=2500,w=1,x∈0,100.y1(x)=2cos⁡x; y2(x)=-cos⁡x where ϵ is an arbitrary parameter.This problem was chosen to demonstrate the stability of the BHTRKNM (Figure5). As mentioned in Remark 7, the method is stable when q∈[0,47.06] and u∈[-π,π].Figure 5
These figures illustrate the stability of the BHTRKNM applied to Example4. In (a) the method is stable with N=722, q∈[0,47.96], and the global error is 1.7×10-10, whereas in (b) the method is unstable with N=721, q>47.96, and the global error is 7005.78.
(a)
0 ≤ q ≤ 47 . 96 (b)
q > 47 . 96


### 4.1. Problems Wherey′ Appears Explicitly

Example 5.
We consider the harmonic oscillator with frequencyΩ and small perturbation δ that was solved in Franco [15] and Guo and Yan [34]:(37)y′′+δy′+Ω2y=0,y0=0,y′0=-δ2,x∈0,1000,where the analytical solution is given by(38)Exact:yx=eδ/2xcos⁡Ω2-δ24,where Ω=1, δ=10-6, and δ=10-10. Guo and Yan [34] solved this problem using ARKN method. The results in Table 4 show that the BHTRKNM is competitive with the order 5 Runge-Kutta-Nyström method.Table 4
Results for Example5.
Our method ARKN N Error(δ=10-6) E r r o r   ( δ = 1 0 - 10 ) N Error(δ=10-6) E r r o r   ( δ = 1 0 - 10 ) 2000 1.82 × 1 0 - 8 2.00 × 1 0 - 12 2000 9.05×10-8 9.00×10-12 4000 1.14 × 1 0 - 9 8.32 × 1 0 - 14 4000 5.43×10-9 7.06×10-13 8000 7.13 × 1 0 - 11 3.50 × 1 0 - 13 8000 2.03×10-10 2.87×10-13 16000 4.33 × 1 0 - 12 1.17 × 1 0 - 13 16000 7.25×10-12 3.56×10-13 32000 5.22 × 1 0 - 14 1.59 × 1 0 - 12 32000 3.45×10-13 5.91×10-13


### 4.2. Hyperbolic PDE

Example 6.
We consider the given Telegraph equation (see Ding et al. [32]). (39)∂2u∂t2+2π∂u∂t+π2u=∂2u∂x2+π2sin⁡πxsin⁡πt+cos⁡πt0≤x≤1,0≤t≤1.The exact solution is given byu(x,y)=sin⁡(πx)sin⁡(πt).In order to solve this PDE using the BHTRKNM, we carry out the semidiscretization of the spatial variablex using the second-order finite difference method to obtain the following second-order system in the second variable t.(40)∂2um∂t2+2π∂um∂t+π2um-um+1-2um+um-1Δx2=gm,0<t<1,m=1,…,M-1,uxm,0=um,utxm,0=um′,where Δx=(b-a)/M, xm=a+mΔx, m=0,1,…,M, u=[u1(t),…,uM(t)]T, g=[g1(t),…,gm(t)]T, um(t)≈u(xm,t), and gm(t)≈g(xm,t)=π2sin⁡πxm(sin⁡πt+cos⁡(πt)), which can be written in the form(41)u′′=ft,u,u′,subject to the boundary conditions u(t0)=u0 and u′(t0)=u0′ where ft,u,u′=Au+g, and A is (M-1)×(M-1), matrix arising from the semidiscretized system, and g is a vector of constants.The boundary conditions are chosen accordingly. This example was chosen to demonstrate that the BHTRKNM can be used to solve the Telegraph equation. In Table5, the results produced by the BHTRKNM using Δt=1/100 and space step Δx=1/100 are compared to scheme (3.12) (λ1=1/12 and λ2=5/6), time step Δt=1/200, and space step Δx=1/100, given in Ding et al. [32]. It is obvious from Table 5 that the BHTRKNM is more accurate than the method given in [32]. Moreover, the errors produced by BHTRKNM method using Δt=1/100 and space step Δx=1/100 are given in Figure 6.Table 5
Results, withω=π, for Example 6.
Our method Ding et al. [32] x Err Err 0.2 2.46 × 1 0 - 10 9.62×10-10 0.4 3.96 × 1 0 - 10 1.56×10-9 0.6 3.98 × 1 0 - 10 1.56×10-9 0.8 2.46 × 1 0 - 10 9.62×10-10Figure 6
Absolute errors for Example6.
(a)
N = M = 10 (b)
N = M = 100


### 4.3. Hamiltonian Systems and Energy Conservation

In this section we present additional examples to show that the BHTRKNM preserves energy. To do so we consider Hamiltonian systems of the form(42)p′=-∇qHp,q,q′=-∇pHp,q,where H(p,q) is an arbitrary scalar function of the variables (p,q). Let M be a positive definite matrix and let U(q) be a potential and the total energy H expressed as the sum of the kinetic and potential energy namely in the form(43)Hp,q=12pTM-1p+Uq;then systems (42) can be written as a system of first-order differential equations(44)q′=v,v′=fq,where the momentra p=Mv is in terms of the velocities and f(q)=-M-1∇U(q) is in terms of the negative gradient of a potential. See [36–39] and references therein for further details. The Hamiltonian function, H(y), defined by H(y)=H(p,q) is a polynomial in the variables p and q. The Hamiltonian function conserves energy if(45)Hyn+1=Hyn,∀n,h>0.Example 7.
We consider the perturbed Kepler’s problem in [40] given by(46)q1′′=-q1q12+q223/2-2ϵ+ϵ2q1q12+q225/2,q10=1,q1′0=0,q2′′=-q2q12+q223/2-2ϵ+ϵ2q2q12+q225/2,q20=0,q2′0=1+ϵ.The exact solution of this problem is(47)q1t=cos⁡t+ϵt,q2t=sin⁡t+ϵt.The Hamiltonian is(48)H=12q1′2+q2′2-1q12+q22-2ϵ+ϵ23q12+q223/2.The system also has the angular momentum L=q1q2′-q2q1′ as a first integral. We take the parameter value ϵ=10-3.This problem is solved in the interval[0,1000] using the BHTRKNM for various values of h=0.1/2i-1, i=0,1,2,3,4. The BHTRKNM preserves the Hamiltonian energy and to demonstrate this, we plot the logarithm of the global error of the Hamiltonian EH=|Hn-H0| and the momentum EL=|Ln-L0| as given in Figures 7(a) and 7(b), respectively. The problem was also solved using N4 given in Sommeijer [35] and in Figure 7(c), the efficiency curves for the BHTRKNM and N4 are compared showing that the BHTRKNM is superior.Figure 7
Perturbed Kepler problem: the logarithm of the global error of the HamiltonianEH=|Hn-H0| against t for h=0.2 and the momentum EL=|Ln-L0| are presented in (a) and (b), respectively. In (c) we compare the efficiency curves for the BHTRKNM and N4. Timing comparison is provided in (d). It is clear from the timing curves that BHTRKNM is very efficient.
(a)
N = 1000; h=0.2 (b)
N = 1000; h=0.2 (c)
Efficiency curves for BHTRKNM and N4 (d)
Timing curves for BHTRKNM and N4 (e)
BHTRKNM solution,h=0.2 (f)
N4 solution,h=0.2Example 8.
We consider the pendulum oscillator in [36] (and references herein) given by(49)q′′=-sin⁡qwith initial conditions(50)q0=0,q′0=1.5and Hamiltonian(51)H=12q′2-cos⁡q.This problem is solved using the BHTRKNM on the interval[0,50] for h=1 and w=1,1.4,2 and the results for the phase diagrams produced by the BHTRKNM in the q-q′ plane are presented in Figures 8(a), 8(b), and 8(c), respectively. We observe that the BHTRKNM gives good results for all the values of w, since all the diagrams follow the exact flow of the pendulum problem as given in Figure 8(e). As illustrated in these Figures, the numerical solutions are periodic and in accordance with the fact that the pendulum equation has a periodic solution. We note that Van Daele and Vanden Berghe in [36] obtained similar results for a smaller step-size h=0.5 and w=1 using S/V method including other versions of the S/V method and it was observed that the SE1EF method in [36] produced better numerical results for w=1 than for w=1.4. The problem was also solved using the fourth-order Runge-Kutta method (RK4) and the results presented in Figure 8(d) show a distortion in the flow diagram for the RK4; hence the BHTRKNM is superior. The pendulum problem was also presented in Figure 9.Figure 8
The pendulum problem: phase diagrams for BHTRKNM withw=1,1.4,2 are presented in (a), (b), and (c), respectively. (d) illustrates a distortion in the flow for the RK4, while (e) shows exact flow of the pendulum problem. In the diagrams, y=q; y′=q′.
(a)
BHTRKNM,w=1; h=1 (b)
BHTRKNM,w=1.4; h=1 (c)
BHTRKNM,w=2; h=1 (d)
RK4,h=1 (e)
Exact flowFigure 9
The pendulum problem: Hamiltonian errorEH=|Hn-H0| using RK4 with h=0.5 is presented in (a), while (b) shows EH for BHTRKNM. We also solve the problem on a larger interval of integration [0,500], using the BHTRNM in (c) and (d) and the RK4 in (e) and (f), respectively.
(a)
RK4-EH,h=0.5 (b)
BHTRKNM-EH,w=1.0; h=0.5 (c)
BHTRKNM,w=1.0; h=0.1 (d)
BHTRKNM-EH,w=1.0; h=0.1 (e)
RK4,w=1.0; h=0.1 (f)
RK4-EH,w=1.0; h=0.1


## 4.1. Problems Wherey′ Appears Explicitly

Example 5.
We consider the harmonic oscillator with frequencyΩ and small perturbation δ that was solved in Franco [15] and Guo and Yan [34]:(37)y′′+δy′+Ω2y=0,y0=0,y′0=-δ2,x∈0,1000,where the analytical solution is given by(38)Exact:yx=eδ/2xcos⁡Ω2-δ24,where Ω=1, δ=10-6, and δ=10-10. Guo and Yan [34] solved this problem using ARKN method. The results in Table 4 show that the BHTRKNM is competitive with the order 5 Runge-Kutta-Nyström method.Table 4
Results for Example5.
Our method ARKN N Error(δ=10-6) E r r o r   ( δ = 1 0 - 10 ) N Error(δ=10-6) E r r o r   ( δ = 1 0 - 10 ) 2000 1.82 × 1 0 - 8 2.00 × 1 0 - 12 2000 9.05×10-8 9.00×10-12 4000 1.14 × 1 0 - 9 8.32 × 1 0 - 14 4000 5.43×10-9 7.06×10-13 8000 7.13 × 1 0 - 11 3.50 × 1 0 - 13 8000 2.03×10-10 2.87×10-13 16000 4.33 × 1 0 - 12 1.17 × 1 0 - 13 16000 7.25×10-12 3.56×10-13 32000 5.22 × 1 0 - 14 1.59 × 1 0 - 12 32000 3.45×10-13 5.91×10-13


## 4.2. Hyperbolic PDE

Example 6.
We consider the given Telegraph equation (see Ding et al. [32]). (39)∂2u∂t2+2π∂u∂t+π2u=∂2u∂x2+π2sin⁡πxsin⁡πt+cos⁡πt0≤x≤1,0≤t≤1.The exact solution is given byu(x,y)=sin⁡(πx)sin⁡(πt).In order to solve this PDE using the BHTRKNM, we carry out the semidiscretization of the spatial variablex using the second-order finite difference method to obtain the following second-order system in the second variable t.(40)∂2um∂t2+2π∂um∂t+π2um-um+1-2um+um-1Δx2=gm,0<t<1,m=1,…,M-1,uxm,0=um,utxm,0=um′,where Δx=(b-a)/M, xm=a+mΔx, m=0,1,…,M, u=[u1(t),…,uM(t)]T, g=[g1(t),…,gm(t)]T, um(t)≈u(xm,t), and gm(t)≈g(xm,t)=π2sin⁡πxm(sin⁡πt+cos⁡(πt)), which can be written in the form(41)u′′=ft,u,u′,subject to the boundary conditions u(t0)=u0 and u′(t0)=u0′ where ft,u,u′=Au+g, and A is (M-1)×(M-1), matrix arising from the semidiscretized system, and g is a vector of constants.The boundary conditions are chosen accordingly. This example was chosen to demonstrate that the BHTRKNM can be used to solve the Telegraph equation. In Table5, the results produced by the BHTRKNM using Δt=1/100 and space step Δx=1/100 are compared to scheme (3.12) (λ1=1/12 and λ2=5/6), time step Δt=1/200, and space step Δx=1/100, given in Ding et al. [32]. It is obvious from Table 5 that the BHTRKNM is more accurate than the method given in [32]. Moreover, the errors produced by BHTRKNM method using Δt=1/100 and space step Δx=1/100 are given in Figure 6.Table 5
Results, withω=π, for Example 6.
Our method Ding et al. [32] x Err Err 0.2 2.46 × 1 0 - 10 9.62×10-10 0.4 3.96 × 1 0 - 10 1.56×10-9 0.6 3.98 × 1 0 - 10 1.56×10-9 0.8 2.46 × 1 0 - 10 9.62×10-10Figure 6
Absolute errors for Example6.
(a)
N = M = 10 (b)
N = M = 100


## 4.3. Hamiltonian Systems and Energy Conservation

In this section we present additional examples to show that the BHTRKNM preserves energy. To do so we consider Hamiltonian systems of the form(42)p′=-∇qHp,q,q′=-∇pHp,q,where H(p,q) is an arbitrary scalar function of the variables (p,q). Let M be a positive definite matrix and let U(q) be a potential and the total energy H expressed as the sum of the kinetic and potential energy namely in the form(43)Hp,q=12pTM-1p+Uq;then systems (42) can be written as a system of first-order differential equations(44)q′=v,v′=fq,where the momentra p=Mv is in terms of the velocities and f(q)=-M-1∇U(q) is in terms of the negative gradient of a potential. See [36–39] and references therein for further details. The Hamiltonian function, H(y), defined by H(y)=H(p,q) is a polynomial in the variables p and q. The Hamiltonian function conserves energy if(45)Hyn+1=Hyn,∀n,h>0.Example 7.
We consider the perturbed Kepler’s problem in [40] given by(46)q1′′=-q1q12+q223/2-2ϵ+ϵ2q1q12+q225/2,q10=1,q1′0=0,q2′′=-q2q12+q223/2-2ϵ+ϵ2q2q12+q225/2,q20=0,q2′0=1+ϵ.The exact solution of this problem is(47)q1t=cos⁡t+ϵt,q2t=sin⁡t+ϵt.The Hamiltonian is(48)H=12q1′2+q2′2-1q12+q22-2ϵ+ϵ23q12+q223/2.The system also has the angular momentum L=q1q2′-q2q1′ as a first integral. We take the parameter value ϵ=10-3.This problem is solved in the interval[0,1000] using the BHTRKNM for various values of h=0.1/2i-1, i=0,1,2,3,4. The BHTRKNM preserves the Hamiltonian energy and to demonstrate this, we plot the logarithm of the global error of the Hamiltonian EH=|Hn-H0| and the momentum EL=|Ln-L0| as given in Figures 7(a) and 7(b), respectively. The problem was also solved using N4 given in Sommeijer [35] and in Figure 7(c), the efficiency curves for the BHTRKNM and N4 are compared showing that the BHTRKNM is superior.Figure 7
Perturbed Kepler problem: the logarithm of the global error of the HamiltonianEH=|Hn-H0| against t for h=0.2 and the momentum EL=|Ln-L0| are presented in (a) and (b), respectively. In (c) we compare the efficiency curves for the BHTRKNM and N4. Timing comparison is provided in (d). It is clear from the timing curves that BHTRKNM is very efficient.
(a)
N = 1000; h=0.2 (b)
N = 1000; h=0.2 (c)
Efficiency curves for BHTRKNM and N4 (d)
Timing curves for BHTRKNM and N4 (e)
BHTRKNM solution,h=0.2 (f)
N4 solution,h=0.2Example 8.
We consider the pendulum oscillator in [36] (and references herein) given by(49)q′′=-sin⁡qwith initial conditions(50)q0=0,q′0=1.5and Hamiltonian(51)H=12q′2-cos⁡q.This problem is solved using the BHTRKNM on the interval[0,50] for h=1 and w=1,1.4,2 and the results for the phase diagrams produced by the BHTRKNM in the q-q′ plane are presented in Figures 8(a), 8(b), and 8(c), respectively. We observe that the BHTRKNM gives good results for all the values of w, since all the diagrams follow the exact flow of the pendulum problem as given in Figure 8(e). As illustrated in these Figures, the numerical solutions are periodic and in accordance with the fact that the pendulum equation has a periodic solution. We note that Van Daele and Vanden Berghe in [36] obtained similar results for a smaller step-size h=0.5 and w=1 using S/V method including other versions of the S/V method and it was observed that the SE1EF method in [36] produced better numerical results for w=1 than for w=1.4. The problem was also solved using the fourth-order Runge-Kutta method (RK4) and the results presented in Figure 8(d) show a distortion in the flow diagram for the RK4; hence the BHTRKNM is superior. The pendulum problem was also presented in Figure 9.Figure 8
The pendulum problem: phase diagrams for BHTRKNM withw=1,1.4,2 are presented in (a), (b), and (c), respectively. (d) illustrates a distortion in the flow for the RK4, while (e) shows exact flow of the pendulum problem. In the diagrams, y=q; y′=q′.
(a)
BHTRKNM,w=1; h=1 (b)
BHTRKNM,w=1.4; h=1 (c)
BHTRKNM,w=2; h=1 (d)
RK4,h=1 (e)
Exact flowFigure 9
The pendulum problem: Hamiltonian errorEH=|Hn-H0| using RK4 with h=0.5 is presented in (a), while (b) shows EH for BHTRKNM. We also solve the problem on a larger interval of integration [0,500], using the BHTRNM in (c) and (d) and the RK4 in (e) and (f), respectively.
(a)
RK4-EH,h=0.5 (b)
BHTRKNM-EH,w=1.0; h=0.5 (c)
BHTRKNM,w=1.0; h=0.1 (d)
BHTRKNM-EH,w=1.0; h=0.1 (e)
RK4,w=1.0; h=0.1 (f)
RK4-EH,w=1.0; h=0.1


## 5. Conclusion

This paper presents a BHTRKNM whose coefficients are functions of the frequency and the step-size for directly solving general second-order initial value problems (IVPs), oscillatory systems, and Hamiltonian systems, as well as systems arising from the semidiscretization of hyperbolic PDEs, such as the Telegraph equation. We implement the BHTRKNM in a block-by-block fashion; thus the method does not need starting values and predictors which are inherent in predictor-corrector methods. Numerical experiments presented in this paper clearly demonstrate that our method has a reasonably wide stability region and enjoys accuracy and efficiency advantages when compared to existing methods in the literature. Technique for accurately estimating the frequency as suggested in [30, 41] as well as implementing the method in a variable step mode will be considered in future.



---

*Source: 9293530-2017-01-22.xml*