# Solving the Telegraph and Oscillatory Differential Equations by a Block Hybrid Trigonometrically Fitted Algorithm

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

---

## Abstract

We propose a block hybrid trigonometrically fitted (BHT) method, whose coefficients are functions of the frequency and the step-size for directly solving general second-order initial value problems (IVPs), including systems arising from the semidiscretization of hyperbolic Partial Differential Equations (PDEs), such as the Telegraph equation. The BHT is formulated from eight discrete hybrid formulas which are provided by a continuous two-step hybrid trigonometrically fitted method with two off-grid points. The BHT 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. The stability property of the BHT 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 the form of (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 [9], and 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]). It is also the case that some of these IVPs possess solutions with special properties that may be known in advance and taken 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 [17], Coleman and Ixaru [18], Simos [19], Vanden Berghe et al. [20], Vigo-Aguiar and Ramos [11], Fang et al. [21], Nguyen et al. [22], Ramos and Vigo-Aguiar [23], Franco and Gómez [24], and Ozawa [25]). However, most of these methods are restricted to solving special second-order IVPs in a predictor-corrector mode.Our objective is to present a BHT 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. [26], Jator [27], 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 [27]. However, the BHT 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. We note that Ramos et al. [30] recently proposed a trigonometrically fitted optimized two-step hybrid block method for solving the general second-order IVPs with oscillatory solutions. However, the method given in [30] is of an order 2; hence, in this paper, we propose a BHT which is of order 5 and its application is extended to solving PDEs such as the Telegraph equation.The organization of this paper is as follows. In Section2, we derive the BHT for solving (1). The analysis and implementation of the BHT are discussed in Section 3. Numerical examples are given in Section 4 to show the accuracy and efficiency of the BHT. Finally, the conclusion of the paper is given in Section 5.


## 2. Development of Method

Consider(2)yn+2=α2,0yn+α2,1yn+1+h2∑j=02β2,jfn+j+β2,n+vfn+v+β2,n+ϵfn+ϵ,yn+v=αv,0yn+αv,1yn+1+h2∑j=02βv,jfn+j+βv,n+vfn+v+βv,n+ϵfn+ϵ,yn+ϵ=αϵ,0yn+αϵ,1yn+1+h2∑j=02βϵ,jfn+j+βϵ,n+vfn+v+βϵ,n+ϵfn+ϵwhich are used together with additional methods given as(3)hyn+2′=α2,0′yn+α2,1′yn+1+h2∑j=02β2,j′fn+j+β2,n+v′fn+v+β2,n+ϵ′fn+ϵ,hyn+v′=αv,0′yn+αv,1′yn+1+h2∑j=02βv,j′fn+j+βv,n+v′fn+v+βv,n+ϵ′fn+ϵ,hyn+1′=α1,0′yn+α1,1′yn+1+h2∑j=02β1,j′fn+j+β1,n+v′fn+v+β1,n+ϵ′fn+ϵ,hyn+ϵ′=αϵ,0′yn+αϵ,1′yn+1+h2∑j=02βϵ,j′fn+j+βϵ,n+v′fn+v+βϵ,n+ϵ′fn+ϵ,hyn′=α0,0′yn+α0,1′yn+1+h2∑j=02β0,j′fn+j+β0,n+v′fn+v+β0,n+ϵ′fn+ϵ,where αj,0, αj,1, and βj,j, j=0,v,1,ϵ,2, are coefficients that depend on the step-length h and frequency w. The coefficients of the method are chosen so that the method integrates IVP (1) exactly where the solutions are members of the linear space 〈1,x,x2,x3,x4,sin(wx),cos(wx)〉.The main method has the form(4)yn+2=∑i=01αiyn+i+h2∑j=02βjfn+j+βn+vfn+v+βn+ϵfn+ϵ,where αi, i=0,1, and βj, j=0,v,1,ϵ,2, are to be determined coefficient functions of the frequency and step-size. To derive the main method and additional methods, we initially seek a continuous local approximation Π(x) on the interval [xn,xn+2] of the form(5)Πx=α0xyn+α1xyn+1+h2∑jβjxfn+j,j=0,v,1,ϵ,2which represents our CHT and where α0(x), α1(x), and βj(x), j=0,v,1,ϵ,2, are continuous coefficients. The first derivative of (5) is given by(6)Π′x=ddxΠx.We assume that yn+j=Π(xn+j) is the numerical approximation to the analytical solution yxn+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,ϵ,2.The following theorem shows how continuous method (5) is constructed. This is done by requiring that on the interval from xn to xn+2=xn+2h the exact solution is locally approximated by function (5) with (6) obtained as a consequence.Theorem 1.
LetFi(x)=xi, i=0,1,2,3,4, F5(x)=sinwx, and F6(x)=coswx be basis functions and K=(yn,yn+1,fn,fn+v,fn+1,fn+ϵ,fn+2)T a vector, where T is the transpose. Define the matrix G by (7)G=F0xn⋯F6xnF0xn+1⋯F6xn+1F0′′xn⋯F6′′xnF0′′xn+v⋯F6′′xn+vF0′′xn+1⋯F6′′xn+1F0′′xn+ϵ⋯F6′′xn+ϵF0′′xn+2⋯F6′′xn+2and Gi is obtained by replacing the ith column of G by the vector K. Let the following conditions be satisfied:(8)Πxn=yn,Π′xn=yn′,Π′′xn+j=fn+j,j=0,v,1,ϵ,2;then continuous representations (5) and (6) are equivalent to the following:(9)Πx=∑i=06det⁡Gidet⁡GFix,(10)Π′x=ddx∑i=06det⁡Gidet⁡GFix.Proof.
We use the approach given in Jator [16] with appropriate notational modification. Let method (5) be defined by the assumed basis functions:(11)α0x=∑i=06αi+1,0Fix,α1x=∑i=06αi+1,1Fix,h2βjx=∑i=06h2βi+1,jFix,where αi+1,0,  αi+1,1,  and h2βi+1,j, are coefficients to be determined. Substituting (11) into (5), we get(12)Πx=∑i=06αi+1,0Fixyn+∑i=06αi+1,1Fixyn+1+∑j∑i=06h2βi+1,jFixfn+jwhich is simplified to(13)Πx=∑i=06αi+1,0Fixyn+αi+1,1Fixyn+1+∑jh2βi+1,jFixfn+jand expressed as(14)Πx=∑i=06liFix,where(15)li=αi+1,0yn+αi+1,1yn+1+∑jh2βi+1,jfn+j.By imposing conditions (8) on (14), we obtain a system of six equations which can be expressed as (16)GL=K,where L=(l0,l1,…,l6)T is a vector whose coefficients are determined via Cramer’s rule as (17)li=det⁡Gidet⁡G,i=0,1,…,6,where Gi is obtained by replacing the ith column of G by K. In order to obtain the continuous approximation, we use the elements of L to rewrite (14) as(18)Πx=∑i=06det⁡Gidet⁡GFix,whose first derivative is given by (19)Π′x=ddx∑i=06det⁡Gidet⁡GFix.Remark 2.
In the derivation of the BHT, the basis functionsFi(x)=xi, i=0,1,2,3,4, F5(x)=sinwx, and F6(x)=coswx are chosen because they are simple to analyze. Nevertheless, other possible bases are possible (see Nguyen et al. [22]).


### 2.1. Specification of the Method

We note that continuous methods (9) and (10) which are equivalent to forms (5) and (6) are used to generate three discrete methods and five additional methods. The discrete and additional methods are then applied as a BHT for solving (1). We choose v=1/2, ϵ=3/2 and evaluating (9) at x=xn+2, x=xn+v, and x=xn+ϵ, respectively, gives the three discrete methods yn+2=Π(xn+2h), yn+v=Π(xn+vh), and yn+ϵ=Π(xn+ϵh) which take the form of the main method. Evaluating (10) at x=xn, x=xn+v, x=xn+1, x=xn+ϵ, and x=xn+2, respectively, gives the additional methods yn′=Π′(xn), yn+v′=Π′(xn+vh), yn+1′=Π′(xn+h), yn+ϵ′=Π′(xn+ϵh), and yn+2′=Π′(xn+2h). The coefficients and their corresponding Taylor series equivalence of yn+v,  yn+ϵ,  yn+2,  hyn′,  hyn+v′,  hyn+1′,  hyn+ϵ′, and hyn+2′ are, respectively, given as follows:(20)αv,0=12,αv,1=12,βv,0=-7u2+13u2-144cos⁡u/2+48cos⁡u+96csc4⁡u/4768u2=-191920-221u21935360-233u4232243200-199u627249868800-1543u834780741632000-26281u10128549621071872000+⋯,βv,v=-h223u2-12u2+32cos⁡u/2+13u2+144cos⁡u+240csc4⁡u/4768u2=-17h2160+79h2u2241920+79h2u429030400+61h2u63406233600+4099h2u847823519744000+1739h2u1016068702633984000+⋯,βv,1=17u2-144cos⁡u/2+u2+72cos⁡u+72csc4⁡u/4384u2=-7960-19u264512-83u438707200-u6100925440+577u863764692992000+773u10856997473812480+⋯,βv,ϵ=-5u2+u2+48cos⁡u/2-48cot2⁡u/4csc2⁡u/4192u2=-1480+u215120+u47257600-u6212889600-167u81707982848000-2633u102008587829248000+⋯,βv,2=5u2+u2+48cos⁡u/2-48csc4⁡u/4768u2=11920+31u21935360+67u4232243200+109u627249868800+18127u8382588157952000+64931u10128549621071872000+⋯,(21)αϵ,0=-12,αϵ,1=32,βϵ,0=-u2+5u2-16cos⁡u/2+16cos⁡ucsc4⁡u/4256u2=171920+53u2645120+11u425804800-19u627249868800-6427u8127529385984000-11509u1014283291230208000+⋯,βϵ,v=23u2-4u2cos⁡u/2+5u2+48cos⁡u-48csc4⁡u/4256u2=21160-u25040-u42419200+u670963200+167u8569327616000+2633u10669529276416000+⋯,βϵ,1=3-7u2+16cos⁡u/2+u2-8cos⁡u+24csc4⁡u/4128u2=67320+11u2107520-17u412902400-173u64541644800-12277u821254897664000-49729u107141645615104000+⋯,βϵ,ϵ=19u2+4u2+32cos⁡u/2+u2+16cos⁡u-144csc4⁡u/4256u2=13480+u216128+u4460800+5u6136249344+691u81449197568000+u10186856243200+⋯,βϵ,2=-5u2+u2+48cos⁡u/2-48csc4⁡u/4256u2=-1640-31u2645120-67u477414400-109u69083289600-18127u8127529385984000-64931u1042849873690624000+⋯,(22)α2,0=1,α2,1=2,β2,0=-u2+2u2cos⁡u/2+6cos⁡u-6csc4⁡u/448u2=160+u27560+u43628800-u6106444800-167u8853991424000-2633u101004293914624000+⋯,β2,v=5u2+u2+12cos⁡u-12csc4⁡u/424u2=415-u21890-u4907200+u626611200+167u8213497856000+2633u10251073478656000+⋯,β2,1=-10u2cos⁡u/2+u2-18cos⁡u+18csc4⁡u/424u2=1330+u21260+u4604800-u617740800-167u8142331904000-2633u10167382319104000+⋯,β2,ϵ=5u2+u2+12cos⁡u-12csc4⁡u/424u2=415-u21890-u4907200+u626611200+167u8213497856000+2633u10251073478656000+⋯,β2,2=-u2+2u2cos⁡u/2+6cos⁡u-6csc4⁡u/448u2=160+u27560+u43628800-u6106444800-167u8853991424000-2633u101004293914624000+⋯,(23)α0,0′=-1,α0,1′=1,β0,0′=csc4⁡u/4-7u2+10u2cos⁡u/2-6cos⁡5u/4sec⁡u/4-6usin⁡5u/4sec⁡u/4+696u2=-53360-19u215120-37u42419200-163u6638668800-461287u883691159552000-12773u1095647039488000+⋯,β0,v′=-csc4⁡u/4sec⁡u/45u2cos⁡5u/4+6u2+4cos⁡u/4+u2-6cos⁡3u/4-6usin⁡3u/4-18usin⁡5u/4-18cos⁡5u/496u2=-25+13u23780+29u653222400+67u41814400+46559u84184557977600+26753u10100429391462400+⋯,β0,1′=csc4⁡u/42u2cos⁡u/2+u2-6cos⁡u-6usin⁡u+616u2=112-u2360-23u41209600-331u81072963584000-u69676800+703u10334764638208000+⋯,β0,ϵ′=-csc4⁡u/4sec⁡u/4u2cos⁡5u/4+6u2+4cos⁡u/4+5u2-18cos⁡3u/4-18usin⁡3u/4-6usin⁡5u/4-6cos⁡5u/496u2=-245-u486400+u23780-13u631933440-32027u82988969984000-15019u1055794106368000+⋯,β0,2′=csc4⁡u/4u2+2u2-6cos⁡u/2-6usin⁡3u/4sec⁡u/4+1296u2=1120+u23024+13u41451520+47u6212889600+452681u883691159552000+269639u102008587829248000+⋯,(24)αv,0′=-1,αv,1′=1,βv,0′=csc3⁡u/4sec⁡u/4u2+8sin⁡u/2-4u+8sin⁡u-8ucos⁡u/264u2=13480+17u232256+67u46451200+1549u66812467200+1205803u8223176425472000+9739u1072873934848000+⋯,βv,v′=-u2+24sin⁡u/2-12u3cos⁡u/2+2csc3⁡u/4sec⁡u/496u2=7144-109u280640-139u45806080-4951u610218700800-246931u822317642547200-4329211u1016068702633984000+⋯,βv,1′=csc4⁡u/4u2cos⁡u/2-u2+18cos⁡u-18usin⁡u/2+1848u2=-110+u21120+23u42419200+19u6212889600+7213u89299017728000+4289u10669529276416000+⋯,βv,ϵ′=-u2+24sin⁡u/2-12ucos⁡u/2+2csc3⁡u/4sec⁡u/496u2=7240+13u280640+109u49676800+83u6227082240+1119247u8111588212736000+1397321u105356234211328000+⋯,βv,2′=csc2⁡u/4u2-12ucsc⁡u/2+2496u2=-71440-37u2161280-419u458060800-577u62919628800-577u62919628800-1148099u8223176425472000-169051u101285496210718720+⋯,(25)α1,0′=-1,α1,1′=1,β1,0′=-csc4⁡u/4sec⁡u/42u2cos⁡3u/4+u2-6cos⁡u/4+6usin⁡u/4+cos⁡5u/496u2=172-u230240-u4172800-61u6319334400-215521u841845579776000-1769u1013390585528320+⋯,β1,v′=csc4⁡u/4sec⁡u/4u2cos⁡5u/4+5u2-12cos⁡u/4+6usin⁡u/4+3cos⁡3u/4+9cos⁡5u/448u2=1345-u215120+41u43628800+5u612773376+17569u81673823191040+267983u101004293914624000+⋯,β1,1′=-10u2cos⁡u/2+u2-18cos⁡u+18csc4⁡u/448u2=1360+u22520+u41209600-u635481600-167u8284663808000-2633u10334764638208000+⋯,β1,ϵ′=csc4⁡u/4sec⁡u/45u2-12cos⁡u/4+u2+9cos⁡3u/4+3cos⁡5u/4-2usin⁡u/448u2=-145-u22160-u480640-113u6319334400-406493u841845579776000-85817u10334764638208000+⋯,β1,2′=-csc4⁡u/4sec⁡u/4u2-2cos⁡u/4+2cos⁡3u/4-usin⁡u/432u2=1360+u26048+11u41814400+29u6159667200+103669u820922789888000+5911u1045649723392000+⋯,(26)αϵ,0′=-1,αϵ,1′=1,βϵ,0′=-csc4⁡u/4sec⁡u/47u2cos⁡3u/4+17u2-24cos⁡u/4+24cos⁡5u/4-usin⁡u/4384u2=311440+5u213824+29u43870720+3847u620437401600+3313369u8669529276416000+28177u10218621804544000+⋯,βϵ,v′=csc4⁡u/4sec⁡u/47u2cos⁡5u/4+127u2-8cos⁡u/4+5u2+24cos⁡3u/4-72usin⁡u/4-24usin⁡3u/4+72cos⁡5u/4384u2=1980-167u2241920-359u429030400-1117u63406233600-47629u85150225203200-4023451u1016068702633984000+⋯,βϵ,1′=csc4⁡u/4-7u2cos⁡u/2+u2-6cos⁡u+6usin⁡u/2+616u2=815-u210080-19u42419200-31u6212889600-54371u827897053184000-14821u10669529276416000+⋯,βϵ,ϵ′=csc4⁡u/4sec⁡u/45u2cos⁡5u/4+127u2-8cos⁡u/4+7u2+72cos⁡3u/4-24usin⁡u/4-72usin⁡3u/4+24cos⁡5u/4384u2=157720+199u2241920+221u49676800+97u6185794560+3965821u8334764638208000+499747u101785411403776000+⋯,βϵ,2′=-csc4⁡u/47u2+5u2+24cos⁡u/2-12usin⁡3u/4sec⁡u/4-24192u2=-196-191u2483840-587u458060800-1613u66812467200-3748337u8669529276416000-79621u10584316459417600+⋯,(27)α2,0′=-1,α2,1′=1,β2,0′=-2cos⁡u/2+1csc4⁡u/4sec⁡u/4u2-2cos⁡u/4+2cos⁡3u/4-usin⁡u/432u2=1120-u25040-u4115200-7u630412800-156349u827897053184000-6109u1044635285094400+⋯,β2,v′=csc4⁡u/4sec⁡u/426u2-24cos⁡u/4+7u2+6cos⁡3u/4+3u2+6cos⁡5u/4-2u3sin⁡3u/4+sin⁡5u/496u2=1445-u21260+19u41814400+71u6159667200+2291u8199264665600+12767u1045649723392000+⋯,β2,1′=-csc4⁡u/426u2cos⁡u/2+u2+18cos⁡u-18usin⁡u+148u2=720+u2280+u448384+u621288960-4021u84649508864000-5969u10334764638208000+⋯,β2,ϵ′=csc4⁡u/4sec⁡u/47u2cos⁡5u/4+26u2-24cos⁡u/4+3u2+6cos⁡3u/4-6usin⁡3u/4-18usin⁡5u/4+6cos⁡5u/496u2=23-u2252-23u4604800-u61971200-72143u86974263296000-6119u1023911759872000+⋯,β2,2′=-csc4⁡u/4-5u2+27u2+6cos⁡u/2-6usin⁡5u/4sec⁡u/4-1296u2=59360+u2720+113u47257600+157u6638668800+148307u827897053184000+262967u102008587829248000+⋯.Remark 3.
We note that the Taylor series expansions in (20) through (27) 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

In this subsection, the BHT method is formulated from the eight discrete hybrid formulas stated in (2) and (3). We emphasize that these eight formulas are provided by the continuous two-step hybrid trigonometrically fitted method with two off-grid points given by (5) and its first derivative (6). First, we define the following vectors:(28)Yμ+1=yn+v,yn+1,yn+ϵ,yn+2,hyn+v′,hyn+1′,hyn+ϵ′,hyn+2′T,Yμ=yn-ϵ,yn-1,yn-v,yn,hyn-v′,hyn-1′,hyn-ϵ′,hyn′T,Fμ+1=fn+v,fn+1,fn+ϵ,fn+2,hfn+v′,hfn+1′,hfn+ϵ,hfn+2′T,Fμ=fn-ϵ,fn-1,fn-v,fn,hfn-v′,hfn-1′,hfn-ϵ′,hfn′T,where μ=0,…,N;  n=0,…,N. The three discrete methods whose coefficients are specified by (2) and the five additional methods in (3) whose coefficients are specified by (20) to (27) are combined to give the BHT method, which is expressed as(29)A1Yμ+1=A0Yμ+h2B0Fμ+B1Fμ+1,where A0, A1, B0, and B1 are matrices of dimension eight whose elements characterize the method and are given by the coefficients of (2) and (3).


## 2.1. Specification of the Method

We note that continuous methods (9) and (10) which are equivalent to forms (5) and (6) are used to generate three discrete methods and five additional methods. The discrete and additional methods are then applied as a BHT for solving (1). We choose v=1/2, ϵ=3/2 and evaluating (9) at x=xn+2, x=xn+v, and x=xn+ϵ, respectively, gives the three discrete methods yn+2=Π(xn+2h), yn+v=Π(xn+vh), and yn+ϵ=Π(xn+ϵh) which take the form of the main method. Evaluating (10) at x=xn, x=xn+v, x=xn+1, x=xn+ϵ, and x=xn+2, respectively, gives the additional methods yn′=Π′(xn), yn+v′=Π′(xn+vh), yn+1′=Π′(xn+h), yn+ϵ′=Π′(xn+ϵh), and yn+2′=Π′(xn+2h). The coefficients and their corresponding Taylor series equivalence of yn+v,  yn+ϵ,  yn+2,  hyn′,  hyn+v′,  hyn+1′,  hyn+ϵ′, and hyn+2′ are, respectively, given as follows:(20)αv,0=12,αv,1=12,βv,0=-7u2+13u2-144cos⁡u/2+48cos⁡u+96csc4⁡u/4768u2=-191920-221u21935360-233u4232243200-199u627249868800-1543u834780741632000-26281u10128549621071872000+⋯,βv,v=-h223u2-12u2+32cos⁡u/2+13u2+144cos⁡u+240csc4⁡u/4768u2=-17h2160+79h2u2241920+79h2u429030400+61h2u63406233600+4099h2u847823519744000+1739h2u1016068702633984000+⋯,βv,1=17u2-144cos⁡u/2+u2+72cos⁡u+72csc4⁡u/4384u2=-7960-19u264512-83u438707200-u6100925440+577u863764692992000+773u10856997473812480+⋯,βv,ϵ=-5u2+u2+48cos⁡u/2-48cot2⁡u/4csc2⁡u/4192u2=-1480+u215120+u47257600-u6212889600-167u81707982848000-2633u102008587829248000+⋯,βv,2=5u2+u2+48cos⁡u/2-48csc4⁡u/4768u2=11920+31u21935360+67u4232243200+109u627249868800+18127u8382588157952000+64931u10128549621071872000+⋯,(21)αϵ,0=-12,αϵ,1=32,βϵ,0=-u2+5u2-16cos⁡u/2+16cos⁡ucsc4⁡u/4256u2=171920+53u2645120+11u425804800-19u627249868800-6427u8127529385984000-11509u1014283291230208000+⋯,βϵ,v=23u2-4u2cos⁡u/2+5u2+48cos⁡u-48csc4⁡u/4256u2=21160-u25040-u42419200+u670963200+167u8569327616000+2633u10669529276416000+⋯,βϵ,1=3-7u2+16cos⁡u/2+u2-8cos⁡u+24csc4⁡u/4128u2=67320+11u2107520-17u412902400-173u64541644800-12277u821254897664000-49729u107141645615104000+⋯,βϵ,ϵ=19u2+4u2+32cos⁡u/2+u2+16cos⁡u-144csc4⁡u/4256u2=13480+u216128+u4460800+5u6136249344+691u81449197568000+u10186856243200+⋯,βϵ,2=-5u2+u2+48cos⁡u/2-48csc4⁡u/4256u2=-1640-31u2645120-67u477414400-109u69083289600-18127u8127529385984000-64931u1042849873690624000+⋯,(22)α2,0=1,α2,1=2,β2,0=-u2+2u2cos⁡u/2+6cos⁡u-6csc4⁡u/448u2=160+u27560+u43628800-u6106444800-167u8853991424000-2633u101004293914624000+⋯,β2,v=5u2+u2+12cos⁡u-12csc4⁡u/424u2=415-u21890-u4907200+u626611200+167u8213497856000+2633u10251073478656000+⋯,β2,1=-10u2cos⁡u/2+u2-18cos⁡u+18csc4⁡u/424u2=1330+u21260+u4604800-u617740800-167u8142331904000-2633u10167382319104000+⋯,β2,ϵ=5u2+u2+12cos⁡u-12csc4⁡u/424u2=415-u21890-u4907200+u626611200+167u8213497856000+2633u10251073478656000+⋯,β2,2=-u2+2u2cos⁡u/2+6cos⁡u-6csc4⁡u/448u2=160+u27560+u43628800-u6106444800-167u8853991424000-2633u101004293914624000+⋯,(23)α0,0′=-1,α0,1′=1,β0,0′=csc4⁡u/4-7u2+10u2cos⁡u/2-6cos⁡5u/4sec⁡u/4-6usin⁡5u/4sec⁡u/4+696u2=-53360-19u215120-37u42419200-163u6638668800-461287u883691159552000-12773u1095647039488000+⋯,β0,v′=-csc4⁡u/4sec⁡u/45u2cos⁡5u/4+6u2+4cos⁡u/4+u2-6cos⁡3u/4-6usin⁡3u/4-18usin⁡5u/4-18cos⁡5u/496u2=-25+13u23780+29u653222400+67u41814400+46559u84184557977600+26753u10100429391462400+⋯,β0,1′=csc4⁡u/42u2cos⁡u/2+u2-6cos⁡u-6usin⁡u+616u2=112-u2360-23u41209600-331u81072963584000-u69676800+703u10334764638208000+⋯,β0,ϵ′=-csc4⁡u/4sec⁡u/4u2cos⁡5u/4+6u2+4cos⁡u/4+5u2-18cos⁡3u/4-18usin⁡3u/4-6usin⁡5u/4-6cos⁡5u/496u2=-245-u486400+u23780-13u631933440-32027u82988969984000-15019u1055794106368000+⋯,β0,2′=csc4⁡u/4u2+2u2-6cos⁡u/2-6usin⁡3u/4sec⁡u/4+1296u2=1120+u23024+13u41451520+47u6212889600+452681u883691159552000+269639u102008587829248000+⋯,(24)αv,0′=-1,αv,1′=1,βv,0′=csc3⁡u/4sec⁡u/4u2+8sin⁡u/2-4u+8sin⁡u-8ucos⁡u/264u2=13480+17u232256+67u46451200+1549u66812467200+1205803u8223176425472000+9739u1072873934848000+⋯,βv,v′=-u2+24sin⁡u/2-12u3cos⁡u/2+2csc3⁡u/4sec⁡u/496u2=7144-109u280640-139u45806080-4951u610218700800-246931u822317642547200-4329211u1016068702633984000+⋯,βv,1′=csc4⁡u/4u2cos⁡u/2-u2+18cos⁡u-18usin⁡u/2+1848u2=-110+u21120+23u42419200+19u6212889600+7213u89299017728000+4289u10669529276416000+⋯,βv,ϵ′=-u2+24sin⁡u/2-12ucos⁡u/2+2csc3⁡u/4sec⁡u/496u2=7240+13u280640+109u49676800+83u6227082240+1119247u8111588212736000+1397321u105356234211328000+⋯,βv,2′=csc2⁡u/4u2-12ucsc⁡u/2+2496u2=-71440-37u2161280-419u458060800-577u62919628800-577u62919628800-1148099u8223176425472000-169051u101285496210718720+⋯,(25)α1,0′=-1,α1,1′=1,β1,0′=-csc4⁡u/4sec⁡u/42u2cos⁡3u/4+u2-6cos⁡u/4+6usin⁡u/4+cos⁡5u/496u2=172-u230240-u4172800-61u6319334400-215521u841845579776000-1769u1013390585528320+⋯,β1,v′=csc4⁡u/4sec⁡u/4u2cos⁡5u/4+5u2-12cos⁡u/4+6usin⁡u/4+3cos⁡3u/4+9cos⁡5u/448u2=1345-u215120+41u43628800+5u612773376+17569u81673823191040+267983u101004293914624000+⋯,β1,1′=-10u2cos⁡u/2+u2-18cos⁡u+18csc4⁡u/448u2=1360+u22520+u41209600-u635481600-167u8284663808000-2633u10334764638208000+⋯,β1,ϵ′=csc4⁡u/4sec⁡u/45u2-12cos⁡u/4+u2+9cos⁡3u/4+3cos⁡5u/4-2usin⁡u/448u2=-145-u22160-u480640-113u6319334400-406493u841845579776000-85817u10334764638208000+⋯,β1,2′=-csc4⁡u/4sec⁡u/4u2-2cos⁡u/4+2cos⁡3u/4-usin⁡u/432u2=1360+u26048+11u41814400+29u6159667200+103669u820922789888000+5911u1045649723392000+⋯,(26)αϵ,0′=-1,αϵ,1′=1,βϵ,0′=-csc4⁡u/4sec⁡u/47u2cos⁡3u/4+17u2-24cos⁡u/4+24cos⁡5u/4-usin⁡u/4384u2=311440+5u213824+29u43870720+3847u620437401600+3313369u8669529276416000+28177u10218621804544000+⋯,βϵ,v′=csc4⁡u/4sec⁡u/47u2cos⁡5u/4+127u2-8cos⁡u/4+5u2+24cos⁡3u/4-72usin⁡u/4-24usin⁡3u/4+72cos⁡5u/4384u2=1980-167u2241920-359u429030400-1117u63406233600-47629u85150225203200-4023451u1016068702633984000+⋯,βϵ,1′=csc4⁡u/4-7u2cos⁡u/2+u2-6cos⁡u+6usin⁡u/2+616u2=815-u210080-19u42419200-31u6212889600-54371u827897053184000-14821u10669529276416000+⋯,βϵ,ϵ′=csc4⁡u/4sec⁡u/45u2cos⁡5u/4+127u2-8cos⁡u/4+7u2+72cos⁡3u/4-24usin⁡u/4-72usin⁡3u/4+24cos⁡5u/4384u2=157720+199u2241920+221u49676800+97u6185794560+3965821u8334764638208000+499747u101785411403776000+⋯,βϵ,2′=-csc4⁡u/47u2+5u2+24cos⁡u/2-12usin⁡3u/4sec⁡u/4-24192u2=-196-191u2483840-587u458060800-1613u66812467200-3748337u8669529276416000-79621u10584316459417600+⋯,(27)α2,0′=-1,α2,1′=1,β2,0′=-2cos⁡u/2+1csc4⁡u/4sec⁡u/4u2-2cos⁡u/4+2cos⁡3u/4-usin⁡u/432u2=1120-u25040-u4115200-7u630412800-156349u827897053184000-6109u1044635285094400+⋯,β2,v′=csc4⁡u/4sec⁡u/426u2-24cos⁡u/4+7u2+6cos⁡3u/4+3u2+6cos⁡5u/4-2u3sin⁡3u/4+sin⁡5u/496u2=1445-u21260+19u41814400+71u6159667200+2291u8199264665600+12767u1045649723392000+⋯,β2,1′=-csc4⁡u/426u2cos⁡u/2+u2+18cos⁡u-18usin⁡u+148u2=720+u2280+u448384+u621288960-4021u84649508864000-5969u10334764638208000+⋯,β2,ϵ′=csc4⁡u/4sec⁡u/47u2cos⁡5u/4+26u2-24cos⁡u/4+3u2+6cos⁡3u/4-6usin⁡3u/4-18usin⁡5u/4+6cos⁡5u/496u2=23-u2252-23u4604800-u61971200-72143u86974263296000-6119u1023911759872000+⋯,β2,2′=-csc4⁡u/4-5u2+27u2+6cos⁡u/2-6usin⁡5u/4sec⁡u/4-1296u2=59360+u2720+113u47257600+157u6638668800+148307u827897053184000+262967u102008587829248000+⋯.Remark 3.
We note that the Taylor series expansions in (20) through (27) 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

In this subsection, the BHT method is formulated from the eight discrete hybrid formulas stated in (2) and (3). We emphasize that these eight formulas are provided by the continuous two-step hybrid trigonometrically fitted method with two off-grid points given by (5) and its first derivative (6). First, we define the following vectors:(28)Yμ+1=yn+v,yn+1,yn+ϵ,yn+2,hyn+v′,hyn+1′,hyn+ϵ′,hyn+2′T,Yμ=yn-ϵ,yn-1,yn-v,yn,hyn-v′,hyn-1′,hyn-ϵ′,hyn′T,Fμ+1=fn+v,fn+1,fn+ϵ,fn+2,hfn+v′,hfn+1′,hfn+ϵ,hfn+2′T,Fμ=fn-ϵ,fn-1,fn-v,fn,hfn-v′,hfn-1′,hfn-ϵ′,hfn′T,where μ=0,…,N;  n=0,…,N. The three discrete methods whose coefficients are specified by (2) and the five additional methods in (3) whose coefficients are specified by (20) to (27) are combined to give the BHT method, which is expressed as(29)A1Yμ+1=A0Yμ+h2B0Fμ+B1Fμ+1,where A0, A1, B0, and B1 are matrices of dimension eight whose elements characterize the method and are given by the coefficients of (2) and (3).


## 3. Error Analysis and Stability

### 3.1. Local Truncation Error (LTE)

Define the local truncation error of (29) as(30)LZx;h=Zμ+1-AZμ+h2BF¯μ+h2CF¯μ+1,where (31)Zμ+1=yxn+v,yxn+1,yxn+ϵ,yxn+2,hy′xn+v,hy′xn+1,hy′xn+ϵ,hy′xn+2T,Zμ=yxn-ϵ,yxn-1,yxn-v,yxn,hy′xn-ϵ,hy′xn-1,hy′xn-v,hy′xnT,F¯μ+1=fxn+ϵ,yn+ϵ,fxn+1,yn+1,fxn+v,yn+v,fxn+2,yn+2,hf′xn+v,yn+v,hf′xn+1,yn+1,hf′xn+ϵ,yn+ϵ,hf′xn+2,yn+2T,F¯μ=fxn-ϵ,yn-ϵ,fxn-1,yn-1,fxn-v,yn-v,fxn,yn,hf′xn-ϵ,yn-ϵ,hf′xn-1,yn-1,hf′xn-v,yn-v,hf′xn,ynT,LZx;h=L1zx;h,L2zx;h,…,L8zx;hTis linear different operator.Suppose thatZ(x) is sufficiently differentiable. Then, a Taylor series expansion of the terms in (30) about the point x gives the following expression for local truncation error:(32)LZx;h=C0Zx+C1hZ′x+⋯+CqhqZqx+⋯,where Ci, i=0,1,…, are constant coefficients (see [16]).Definition 4.
Block method (29) has algebraic order of 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,yn+ϵ,yn+2,hyn+v′,hyn+1′,hyn+ϵ′,hyn+2′)T of block method (29) are given, respectively, by C¯7=(-1/61440,-1/5040,1/61440,0,-1/5040, 61/645120,-1/40320,61/645120,-1/5040)T, where C¯0=C¯1=C¯2=C¯3=C¯4=C¯5=C¯6=0.
(ii) From the local truncation error constant computation, it follows that method (29) has order p of at least five.


### 3.2. Stability

The linear stability of the BHT is discussed by applying the method to the test equationy′′=-λ2y, where λ is a real constant (see [17]). Letting Υ=λh, it is easily shown as in [18] that the application of (29) to the test equation yields(33)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 [21], the spectral radius of ρ(M(Υ2;u)) can be obtained from the characteristics equation(34)ρ2-2ΓΥ2;uρ+ΘΥ2;u=0,where Γ(Υ2;u)=traceM(Υ2;u) and Θ(Υ2;u)=detM(Υ2;u) are rational functions.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 [21]). We note that the plot for the stability region of the BHT method is given in Figure 1.Figure 1
The stability region for the BHT plotted in the(q,u)-plane.


### 3.3. Implementation

We emphasize that the main method and the additional methods specified by (20)–(27) are combined to form block method BHT (29), which is used to solve (1) without requiring starting values and predictors. BHT 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 [35]). 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 BHT is applied.Step 1.
ChooseN, h=(b-a)/N, and the number of blocks Γ=N/2. Using (29), n=0, μ=0, the values of (y1/2,y1,y3/2,y2)T and (y1/2′,y1′,y3/2′,y2′)T are simultaneously obtained over the subinterval [x0,x2], as y0 and y0′ are known from IVP (1).Step 2.
Forn=2, μ=1, the values of (y5/2,y3,y7/2,y4)T and (y5/2′,y3′,y7/2′,y4′)T are simultaneously obtained over the subinterval [x2,x4], as y2 and y2′ are known from the previous block.Step 3.
The process is continued forn=4,…,N-2 and μ=2,…,Γ to obtain the numerical solution to (1) on the subintervals [x0,x2],[x2,x4],…,[xN-2,xN].


## 3.1. Local Truncation Error (LTE)

Define the local truncation error of (29) as(30)LZx;h=Zμ+1-AZμ+h2BF¯μ+h2CF¯μ+1,where (31)Zμ+1=yxn+v,yxn+1,yxn+ϵ,yxn+2,hy′xn+v,hy′xn+1,hy′xn+ϵ,hy′xn+2T,Zμ=yxn-ϵ,yxn-1,yxn-v,yxn,hy′xn-ϵ,hy′xn-1,hy′xn-v,hy′xnT,F¯μ+1=fxn+ϵ,yn+ϵ,fxn+1,yn+1,fxn+v,yn+v,fxn+2,yn+2,hf′xn+v,yn+v,hf′xn+1,yn+1,hf′xn+ϵ,yn+ϵ,hf′xn+2,yn+2T,F¯μ=fxn-ϵ,yn-ϵ,fxn-1,yn-1,fxn-v,yn-v,fxn,yn,hf′xn-ϵ,yn-ϵ,hf′xn-1,yn-1,hf′xn-v,yn-v,hf′xn,ynT,LZx;h=L1zx;h,L2zx;h,…,L8zx;hTis linear different operator.Suppose thatZ(x) is sufficiently differentiable. Then, a Taylor series expansion of the terms in (30) about the point x gives the following expression for local truncation error:(32)LZx;h=C0Zx+C1hZ′x+⋯+CqhqZqx+⋯,where Ci, i=0,1,…, are constant coefficients (see [16]).Definition 4.
Block method (29) has algebraic order of 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,yn+ϵ,yn+2,hyn+v′,hyn+1′,hyn+ϵ′,hyn+2′)T of block method (29) are given, respectively, by C¯7=(-1/61440,-1/5040,1/61440,0,-1/5040, 61/645120,-1/40320,61/645120,-1/5040)T, where C¯0=C¯1=C¯2=C¯3=C¯4=C¯5=C¯6=0.
(ii) From the local truncation error constant computation, it follows that method (29) has order p of at least five.


## 3.2. Stability

The linear stability of the BHT is discussed by applying the method to the test equationy′′=-λ2y, where λ is a real constant (see [17]). Letting Υ=λh, it is easily shown as in [18] that the application of (29) to the test equation yields(33)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 [21], the spectral radius of ρ(M(Υ2;u)) can be obtained from the characteristics equation(34)ρ2-2ΓΥ2;uρ+ΘΥ2;u=0,where Γ(Υ2;u)=traceM(Υ2;u) and Θ(Υ2;u)=detM(Υ2;u) are rational functions.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 [21]). We note that the plot for the stability region of the BHT method is given in Figure 1.Figure 1
The stability region for the BHT plotted in the(q,u)-plane.


## 3.3. Implementation

We emphasize that the main method and the additional methods specified by (20)–(27) are combined to form block method BHT (29), which is used to solve (1) without requiring starting values and predictors. BHT 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 [35]). 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 BHT is applied.Step 1.
ChooseN, h=(b-a)/N, and the number of blocks Γ=N/2. Using (29), n=0, μ=0, the values of (y1/2,y1,y3/2,y2)T and (y1/2′,y1′,y3/2′,y2′)T are simultaneously obtained over the subinterval [x0,x2], as y0 and y0′ are known from IVP (1).Step 2.
Forn=2, μ=1, the values of (y5/2,y3,y7/2,y4)T and (y5/2′,y3′,y7/2′,y4′)T are simultaneously obtained over the subinterval [x2,x4], as y2 and y2′ are known from the previous block.Step 3.
The process is continued forn=4,…,N-2 and μ=2,…,Γ to obtain the numerical solution to (1) on the subintervals [x0,x2],[x2,x4],…,[xN-2,xN].


## 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]:(35)y′′=-100y+99sinx,y0=1,y′0=11,x∈0,1000,where the analytical solution is given by(36)Exact:yx=cos⁡10x+sin⁡10x+sin⁡x.This example was solved using the order 5 BHT and the end-point 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 BHT 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. Our method (BHT) requires fewer number of function evaluations. Hence, for this example, BHT performs better.Table 1
Results, withω=10, for Example 1.
Our method
Simos [8]
N
Err
N
Err
1000
1.9
×
1
0
-
3
1000
1.4
×
1
0
-
1
2000
8.9
×
1
0
-
6
2000
3.5
×
1
0
-
2
4000
4.2
×
1
0
-
8
4000
1.1
×
1
0
-
3
8000
9.7
×
1
0
-
11
8000
8.4
×
1
0
-
5
16000
6.7
×
1
0
-
11
16000
5.5
×
1
0
-
6
32000
4.3
×
1
0
-
13
32000
3.5
×
1
0
-
7Figure 2
Efficiency curve for Example1.Example 2.
We consider the nonlinear Duffing equation which was also solved by Simos [8] and Ixaru and Berghe [31]:(37)y′′+y+y3=Bcos⁡Ωx,y0=C0,y′0=0.The analytical solution is given by(38)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 end-point global errors for our method with those of Simos [8] and Ixaru and 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 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 Berghe [31]
N
Err
N
Err
N
Err
300
7.7
×
1
0
-
5
300
1.7
×
1
0
-
3
300
1.1
×
1
0
-
3
600
1.7
×
1
0
-
6
600
1.9
×
1
0
-
4
600
5.4
×
1
0
-
5
1200
1.4
×
1
0
-
8
1200
1.4
×
1
0
-
5
1200
1.9
×
1
0
-
6
2400
1.9
×
1
0
-
10
2400
8.7
×
1
0
-
7
2400
6.2
×
1
0
-
8Figure 3
Efficiency curves for Example2.Example 3.
We consider the nonlinear perturbed system on the range[0,10], with ε=10-3 that was also solved in Fang et al. [21]:(39)y1′′+25y1+εy12+y22=εφ1x,y10=1,y1′0=0,y2′′+25y2+εy12+y22=εφ2x,y20=ε,y2′0=5,where(40)φ1x=1+ε2+2εsin⁡5x+x2+2cos⁡x2+25-4x2sin⁡x2,φ2x=1+ε2+2εsin⁡5x+x2-2sin⁡x2+25-4x2cos⁡x2,and the exact solution is given by y1(x)=cos(5x)+εsin(x2), y2(x)=sin(5x)+εcos(x2), representing a periodic motion of constant frequency with small perturbation of variable frequency.This problem was solved using the BHT and the maximum global errors (Err=Max|y(x)-y|) obtained were compared to the variable step-size trigonometrically fitted Runge-Kutta-Nyström method (TFARKN5(3)) given in Fang et al. [21] and a Runge-Kutta-Nyström method (ARKN5(3)) which was constructed by Franco [15]. In Table 3, the maximum global errors for the three methods are compared. In general, the TFARKN5(3) and ARKN5(3) are expected to perform better because of their variable-step implementation advantage. Nevertheless, the BHT which is implemented in fixed step-size mode is highly competitive to these methods.Table 3
A comparison of methods for Example3.
ARKN5(3)
TFARKN5(3)
BHT
N (rejected)
-
l
o
g
10(Err)
N (rejected)
-
l
o
g
10(Err)
N
-
l
o
g
10(Err)
42 (15)
2.82
29 (6)
2.78
50
3.42
86 (7)
4.96
88 (9)
5.33
100
4.61
260 (5)
7.16
262 (8)
7.85
260
7.52
812 (3)
9.37
811 (4)
10.38
810
10.43Example 4.
We consider the following two-body problem which was solved by Ozawa [25] on [0,50π]: (41)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(42)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 Table4 that the results obtained using the BHT method are more accurate than the explicit singly diagonally implicit Runge-Kutta (ESDIRK) and the functionally fitted ESDIRK (FESDIRK) methods given in Ozawa [25]. In Figure 4, we also illustrate the efficiency advantage of the BHT method over those in Ozawa [25].Table 4
Steps and absolute errors, withω=1, for Example 4  [0,50π].
Our method
FESDIRK4(3) [25]
ESDIRK4(3) [25]
Steps
Err
Steps
Err
Steps
Err
220
3.52
×
1
0
-
5
170
2.866
×
1
0
-
1
277
2.153
×
1
0
0
300
1.13
×
1
0
-
6
225
7.846
×
1
0
-
3
496
1.494
×
1
0
-
1
400
1.03
×
1
0
-
7
381
1.399
×
1
0
-
3
884
9.359
×
1
0
-
3
600
3.49
×
1
0
-
10
680
1.690
×
1
0
-
4
1573
6.200
×
1
0
-
4
800
1.14
×
1
0
-
11
1207
1.846
×
1
0
-
5
2796
4.416
×
1
0
-
5
1000
7.68
×
1
0
-
13
2144
1.938
×
1
0
-
6
4970
3.412
×
1
0
-
6
1200
2.8
×
1
0
-
14
3806
1.993
×
1
0
-
7
8833
2.848
×
1
0
-
7
2400
1.02
×
1
0
-
13
6762
2.021
×
1
0
-
8
15706
2.530
×
1
0
-
8Figure 4
Efficiency curve for Example4.


### 4.1. Problems Wherey′ Appears Explicitly

Example 5 (Bessel’s IVP).
We consider the Bessel differential equation that was also solved by Vigo-Aguiar and Ramos [11]:(43)x2y′′+xy′+x2-0.25y=0,y1=2πsin⁡1≈0.671397071418031,y′1=2cos⁡1-sin⁡12π≈0.0954005144474746,x∈1,8,where the exact (analytical) solution is given by(44)Exact:yx=J1/2x=2πxsin⁡x.This problem was chosen to demonstrate the performance of our method on the general second-order IVP with variable coefficients. We compare our results with the variable-step Falkner method of order eight that was implemented in predictor-corrector mode by Vigo-Aguiar and Ramos [32]. The results displayed in Table 5 show that the BHT method performs better.Table 5
Results, withω=1, for Example 5.
Our method
Vigo-Aguiar and Ramos [32]
N
Err
N
Err
67
1.14
×
1
0
-
9
67
7.11
×
1
0
-
7
82
3.5
×
1
0
-
10
82
9.26
×
1
0
-
8
97
1.3
×
1
0
-
10
97
87.8
×
1
0
-
9
112
5.5
×
1
0
-
11
112
1.12
×
1
0
-
10
125
2.9
×
1
0
-
11
125
2.71
×
1
0
-
11Example 6.
We consider the harmonic oscillator with frequencyΩ and small perturbation δ that was solved in Franco [15] and Guo and Yan [36]:(45)y′′+δy′+Ω2y=0,y0=0,y′0=-δ2,x∈0,1000,where the analytical solution is given by(46)Exact:yx=eδ/2xcos⁡Ω2-δ24,where Ω=1, δ=10-6, and δ=10-10. The problem was solved in Guo and Yan [36] using ARKN method. In Table 6, the errors are compared at x=1000. We observed that the BHT is competitive with the order 5 Runge-Kutta-Nyström method.Table 6
Errors atx=1000, for Example 6.
BHT
ARKN
h
Error(δ=10-6)
Error(δ=10-10)
h
Error(δ=10-6)
Error(δ=10-10)
1
4.12
×
1
0
-
8
1.11
×
1
0
-
11
1/2
9.05
×
1
0
-
8
9.00
×
1
0
-
12
1/2
7.06
×
1
0
-
10
2.19
×
1
0
-
13
1/4
5.43
×
1
0
-
9
7.06
×
1
0
-
13
1/4
1.23
×
1
0
-
11
3.12
×
1
0
-
13
1/8
2.03
×
1
0
-
10
2.87
×
1
0
-
13
1/8
5.23
×
1
0
-
12
5.44
×
1
0
-
12
1/16
7.25
×
1
0
-
12
3.56
×
1
0
-
13
1/16
5.62
×
1
0
-
12
2.94
×
1
0
-
12
1/32
3.45
×
1
0
-
13
5.91
×
1
0
-
13


### 4.2. Hyperbolic PDEs

Example 7.
We consider the given Telegraph equation (see Ding et al. [33]): (47)∂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 BHT, 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:(48)∂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(49)u′′=ft,u,u′,subject to the boundary conditions u(t0)=u0, u′t0=u0′, where f(t,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 BHT can be used to solve the Telegraph equation. In Table7, the results produced by the BHT 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. [33]. It is obvious from Table 5 that the BHT is more accurate than the method given in [33]. Moreover, the errors produced by BHT method using Δt=1/100 and space step Δx=1/100 are given in Figure 5.Table 7
Results, withω=π, for Example 7.
Our method
Ding et al. [33]
x
Err
Err
0.2
2.06
×
1
0
-
10
9.62
×
1
0
-
10
0.4
3.33
×
1
0
-
10
1.56
×
1
0
-
9
0.6
3.33
×
1
0
-
10
1.56
×
1
0
-
9
0.8
2.06
×
1
0
-
10
9.62
×
1
0
-
10Figure 5
Absolute errors for Example7.
(a)
N
=
M
=
10
(b)
N
=
M
=
100Example 8.
We consider the wave equation given in Franco [15]. A problem representing a vibrating string with speed ω is given by(50)∂2u∂t2-x1-x∂2u∂x2+ω2-2u=00<x<1,0<t≤5,u0,t=0,u1,t=0,ux,0=x1-x,utx,0=0,where the initial and Dirichlet boundary conditions have been chosen such that the solution is given by u(x,t)=x(1-x)cosωt. In order to solve this PDE using the BHT, we carry out the semidiscretization of the spatial variable x using the second-order finite difference method to obtain the following second-order system in the second variable t:(51)∂2um∂t2-xm1-xm+um+1-2um+um-1Δx2+ω2-2um=gm,m=1,…,M-1,uxm,0=xm1-xm,utxm,0=0,0<t≤5,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)=0, which can be written in the form(52)u′′=ft,u,subject to the boundary conditions u(t0)=u0, u′(t0)=u0′, where f(t,u)=Au+g, and A is M-1×M-1, matrix arising from the semidiscretized system, and g is a vector of constants.In Figure6, we give the errors produced by the BHT which show that the method performs very well on this problem.Figure 6
Absolute errors for Example8.
(a)
N
=
M
=
2
(b)
N
=
M
=
4Example 9.
We consider the following mildly stiff IVP which was also solved in [34]:(53)y′′=-1001y′-1000y,y0=1,y′0=-1,x∈0,10Exact:yx=e-x.This example is given to show that the method still performs well on problems with nontrigonometric solutions. The problem was solved using the BHT and the results obtained were compared with the polynomial based method given in [34] and the standard fourth-order Runge-Kutta method (RK4). The results given in Table 8 show that the BHT is competitive with the method in [34] and is superior to RK4 which are designed for problems with nontrigonometric solutions.Table 8
Comparison of methods for Example9.
N
Method in [34]
Our method
N
RK4
10
1.12
×
1
0
-
4
2.23
×
1
0
-
4
250
Overflow
20
1.69
×
1
0
-
6
3.36
×
1
0
-
6
500
Overflow
40
1.22
×
1
0
-
8
2.44
×
1
0
-
8
1000
Overflow
80
9.79
×
1
0
-
11
1.96
×
1
0
-
10
2000
Overflow
160
1.06
×
1
0
-
12
2.13
×
1
0
-
12
4000
1.20
×
1
0
-
13
320
1.54
×
1
0
-
14
2.48
×
1
0
-
14
8000
7.94
×
1
0
-
15


## 4.1. Problems Wherey′ Appears Explicitly

Example 5 (Bessel’s IVP).
We consider the Bessel differential equation that was also solved by Vigo-Aguiar and Ramos [11]:(43)x2y′′+xy′+x2-0.25y=0,y1=2πsin⁡1≈0.671397071418031,y′1=2cos⁡1-sin⁡12π≈0.0954005144474746,x∈1,8,where the exact (analytical) solution is given by(44)Exact:yx=J1/2x=2πxsin⁡x.This problem was chosen to demonstrate the performance of our method on the general second-order IVP with variable coefficients. We compare our results with the variable-step Falkner method of order eight that was implemented in predictor-corrector mode by Vigo-Aguiar and Ramos [32]. The results displayed in Table 5 show that the BHT method performs better.Table 5
Results, withω=1, for Example 5.
Our method
Vigo-Aguiar and Ramos [32]
N
Err
N
Err
67
1.14
×
1
0
-
9
67
7.11
×
1
0
-
7
82
3.5
×
1
0
-
10
82
9.26
×
1
0
-
8
97
1.3
×
1
0
-
10
97
87.8
×
1
0
-
9
112
5.5
×
1
0
-
11
112
1.12
×
1
0
-
10
125
2.9
×
1
0
-
11
125
2.71
×
1
0
-
11Example 6.
We consider the harmonic oscillator with frequencyΩ and small perturbation δ that was solved in Franco [15] and Guo and Yan [36]:(45)y′′+δy′+Ω2y=0,y0=0,y′0=-δ2,x∈0,1000,where the analytical solution is given by(46)Exact:yx=eδ/2xcos⁡Ω2-δ24,where Ω=1, δ=10-6, and δ=10-10. The problem was solved in Guo and Yan [36] using ARKN method. In Table 6, the errors are compared at x=1000. We observed that the BHT is competitive with the order 5 Runge-Kutta-Nyström method.Table 6
Errors atx=1000, for Example 6.
BHT
ARKN
h
Error(δ=10-6)
Error(δ=10-10)
h
Error(δ=10-6)
Error(δ=10-10)
1
4.12
×
1
0
-
8
1.11
×
1
0
-
11
1/2
9.05
×
1
0
-
8
9.00
×
1
0
-
12
1/2
7.06
×
1
0
-
10
2.19
×
1
0
-
13
1/4
5.43
×
1
0
-
9
7.06
×
1
0
-
13
1/4
1.23
×
1
0
-
11
3.12
×
1
0
-
13
1/8
2.03
×
1
0
-
10
2.87
×
1
0
-
13
1/8
5.23
×
1
0
-
12
5.44
×
1
0
-
12
1/16
7.25
×
1
0
-
12
3.56
×
1
0
-
13
1/16
5.62
×
1
0
-
12
2.94
×
1
0
-
12
1/32
3.45
×
1
0
-
13
5.91
×
1
0
-
13


## 4.2. Hyperbolic PDEs

Example 7.
We consider the given Telegraph equation (see Ding et al. [33]): (47)∂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 BHT, 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:(48)∂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(49)u′′=ft,u,u′,subject to the boundary conditions u(t0)=u0, u′t0=u0′, where f(t,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 BHT can be used to solve the Telegraph equation. In Table7, the results produced by the BHT 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. [33]. It is obvious from Table 5 that the BHT is more accurate than the method given in [33]. Moreover, the errors produced by BHT method using Δt=1/100 and space step Δx=1/100 are given in Figure 5.Table 7
Results, withω=π, for Example 7.
Our method
Ding et al. [33]
x
Err
Err
0.2
2.06
×
1
0
-
10
9.62
×
1
0
-
10
0.4
3.33
×
1
0
-
10
1.56
×
1
0
-
9
0.6
3.33
×
1
0
-
10
1.56
×
1
0
-
9
0.8
2.06
×
1
0
-
10
9.62
×
1
0
-
10Figure 5
Absolute errors for Example7.
(a)
N
=
M
=
10
(b)
N
=
M
=
100Example 8.
We consider the wave equation given in Franco [15]. A problem representing a vibrating string with speed ω is given by(50)∂2u∂t2-x1-x∂2u∂x2+ω2-2u=00<x<1,0<t≤5,u0,t=0,u1,t=0,ux,0=x1-x,utx,0=0,where the initial and Dirichlet boundary conditions have been chosen such that the solution is given by u(x,t)=x(1-x)cosωt. In order to solve this PDE using the BHT, we carry out the semidiscretization of the spatial variable x using the second-order finite difference method to obtain the following second-order system in the second variable t:(51)∂2um∂t2-xm1-xm+um+1-2um+um-1Δx2+ω2-2um=gm,m=1,…,M-1,uxm,0=xm1-xm,utxm,0=0,0<t≤5,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)=0, which can be written in the form(52)u′′=ft,u,subject to the boundary conditions u(t0)=u0, u′(t0)=u0′, where f(t,u)=Au+g, and A is M-1×M-1, matrix arising from the semidiscretized system, and g is a vector of constants.In Figure6, we give the errors produced by the BHT which show that the method performs very well on this problem.Figure 6
Absolute errors for Example8.
(a)
N
=
M
=
2
(b)
N
=
M
=
4Example 9.
We consider the following mildly stiff IVP which was also solved in [34]:(53)y′′=-1001y′-1000y,y0=1,y′0=-1,x∈0,10Exact:yx=e-x.This example is given to show that the method still performs well on problems with nontrigonometric solutions. The problem was solved using the BHT and the results obtained were compared with the polynomial based method given in [34] and the standard fourth-order Runge-Kutta method (RK4). The results given in Table 8 show that the BHT is competitive with the method in [34] and is superior to RK4 which are designed for problems with nontrigonometric solutions.Table 8
Comparison of methods for Example9.
N
Method in [34]
Our method
N
RK4
10
1.12
×
1
0
-
4
2.23
×
1
0
-
4
250
Overflow
20
1.69
×
1
0
-
6
3.36
×
1
0
-
6
500
Overflow
40
1.22
×
1
0
-
8
2.44
×
1
0
-
8
1000
Overflow
80
9.79
×
1
0
-
11
1.96
×
1
0
-
10
2000
Overflow
160
1.06
×
1
0
-
12
2.13
×
1
0
-
12
4000
1.20
×
1
0
-
13
320
1.54
×
1
0
-
14
2.48
×
1
0
-
14
8000
7.94
×
1
0
-
15


## 5. Conclusion

We have presented a BHT method whose coefficients are functions of the frequency and the step-size for directly solving general second-order initial value problems (IVPs), including systems arising from the semidiscretization of hyperbolic PDEs, such as the Telegraph equation. The BHT 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. We have also shown that the BHT method has a reasonably wide stability region and enjoys accuracy and efficiency advantages when compared to existing methods in the literature. Our future research will be to incorporate a technique for accurately estimating the frequency as suggested in [37, 38] as well as implementing the method in a variable-step mode.



---

*Source: 347864-2015-11-24.xml*