# Functionally Fitted Block Method for Solving the General Oscillatory Second-Order Initial Value Problems and Hyperbolic Partial Differential Equations

Authors: S. N. Jator; F. F. Ngwane; N. O. Kirby
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

We present a block hybrid functionally fitted Runge–Kutta–Nyström method (BHFNM) which is dependent on the stepsize and a fixed frequency. Since the method is implemented in a block-by-block fashion, the method does not require starting values and predictors inherent to other predictor-corrector methods. Upon deriving our method, stability is illustrated, and it is used to numerically solve the general second-order initial value problems as well as hyperbolic partial differential equations. In doing so, we demonstrate the method’s relative accuracy and efficiency.


---

## Body

## 1. Introduction

We consider the numerical solution of the second-order IVP(1)y′′=fx,y,y′,yx0=y0,y′x0=y0′,x∈a,b,where f:R×Rm×Rm→Rm is a smooth function and m is the dimension of the system and its application to hyperbolic partial differential equations. Frequently, whether analyzing or numerically solving systems like (1), the system is reduced to an equivalent system of first-order IVP of dimension 2m to leverage methods developed for first-order systems (cf. Lambert [1, 2], Hairer et al. [3], Hairer [4], and Brugnano et al. [5, 6] in the numerical setting).There are a wide range of methods that avoid order-reduction provided they have the special formy′′=f(x,y). Importantly, these methods require less storage space and fewer function evaluations (cf. Hairer [4], Hairer et al. [7], Simos [8], Lambert et al., and [9], Twizell et al. [10]). Fewer methods have been developed which both avoid the order-reduction and solve the general second-order IVP (1) (see Vigo-Aguiar et al. [11], Awoyemi [12], Chawla and Sharma [13], Mahmoud and Osman [14], Franco [15], and Jator [16, 17]). Still other methods have been developed for classes of second-order IVP which use qualitative features of solutions (see Coleman and Duxbury [18], Coleman and Ixaru [19], Simos [20], Vanden et al. [21], Vigo-Aguiar et al. [11], Fang et al. [22], Nguyen et al. [23], Ramos and Vigo-Aguiar [24], Franco and Gomez [25], and Ozawa [26]).Our objective is to present a block hybrid functionally fitted Runge–Kutta–Nyström method (BHFNM) that is implemented in a block-by-block fashion, which does not suffer from the disadvantages of requiring starting values and predictors inherent to 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 [16, 29]. However, the BHFNM 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. It is also the case that the proposed method performs better than those given in Ngwane and Jator [30, 31]. In this paper, we propose a BHFNM which is of order 5 and its application is extended to solving oscillatory systems and hyperbolic partial differential equations. It is in the same spirit as those presented by Ngwane and Jator [30, 31].The organization of this article is as follows. In Section2, we derive the BHFNM for solving (1). The analysis and implementation of the BHFNM are discussed in Section 3. Numerical examples are given in Section 4 to show the accuracy and efficiency of the BHFNM. Finally, the conclusion of the paper is given in Section 5.


## 2. Development of the BHFNM

To simplify the presentation we introduce the following notation: leth>0 denote the stepsize in the method, xi=x0+ih for any real number i, and (2)J=r<s<v<1be a set of positive real numbers that are to be chosen, which will be used to specify intermediate points used in the method described below. Further, let J0=J∪{0}.In order to numerically integrate (1), we define the BHFNM which is applied on the partition QN={xn:n=0,1,…,N}, in which the step [xn,yn,yn′]↦[xn+1,yn+1,yn+1′] is given by combining the main methods: for each j∈J,(3)yn+j=αj,0yn+hγj,0yn′+h2∑k∈J0βj,kfn+k,and(4)hyn+j′=hγj,0′yn′+h2∑k∈J0βj,k′fn+k,where αj,0, γj,0, γj,0′, βj,k, βj,k′ for j∈J and k∈J0 are coefficients to be determined and fn+j=f(xn+j,yn+j,yn+j′). We note that yn+j is the numerical approximation to the analytical solution evaluated at xn+j.The coefficients of the method are chosen so that the method integrates the IVP (1) exactly whenever the solutions are members of the linear space 〈1,x,x2,x3,x4,sin⁡wx,cos⁡(wx)〉, and they will depend on the frequency w as well as the step-length h.In order to derive the main methods and additional methods, we initially seek a continuous local approximationΠ(x) on the interval [xn,xn+1], expressed in vector form as(5)Πx=Ω0x⋯Ω6xl0l1⋮l6,where L=[l0,l1,…,l6]⊤ is a vector coefficients to be uniquely determined and Ωi(x)=xi, i=0,1,2,3,4, Ω5(x)=sin⁡wx, and Ω6(x)=cos⁡wx are basis functions. In order to determine the coefficients given by L, we impose the following conditions on (5):(6)Πxn=yn,Π′xn=yn′,Π′′xn+j=fn+j,j∈J0,and we obtain a system of equations which is expressed as(7)Ω0xn⋯Ω6xnΩ0′xn⋯Ω6′xnΩ0′′xn⋯Ω6′′xnΩ0′′xn+r⋯Ω6′′xn+rΩ0′′xn+s⋯Ω6′′6xn+sΩ0′′xn+v⋯Ω6′′xn+vΩ0′′xn+1⋯Ω6′′xn+1l0l1l2l3l4l5l6=ynyn′fnfn+rfn+sfn+vfn+1We note thatyn+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), for j∈J0.System (7) is solved with the aid of Mathematica to obtain the coefficients given by L, which are substituted into (5) to give the continuous scheme given by(8)Πx=Ω0x⋯Ω6xΩ0xn⋯Ω6xnΩ0′xn⋯Ω6′xnΩ0′′xn⋯Ω6′′xnΩ0′′xn+r⋯Ω6′′xn+rΩ0′′xn+s⋯Ω6′′6xn+sΩ0′′xn+v⋯Ω6′′xn+vΩ0′′xn+1⋯Ω6′′xn+1-1ynyn′fnfn+rfn+sfn+vfn+1which after simplification takes the form(9)Πx=α0xyn+δ0xhyn′+h2∑j∈J0βjxfn+j,where α0(x), δ0(x), and βj(x) for j∈J0 are continuous coefficients. The first derivative of (9) is given by(10)Π′x=ddxΠx.Remark 1.
We note that in the derivation of the BHFNM, the basis functionsΩi(x)=xi, i=0,1,2,3,4,Ω5(x)=sin⁡wx, and Ω6(x)=cos⁡wx are chosen because they are simple to analyze. Nevertheless, other possible bases are possible (see Nguyen et al. [23] and Hoang et al. [23]). We note that the zeros of the Chebyshev’s polynomial of the first kind or the zeros of the Legendre’s polynomial can also be chosen. However, we elect not to use them in this case because upon testing them the results were found to be the same as for the chosen points.


### 2.1. Specification of the Method

The continuous form (8) and its first derivative, which are equivalent to the forms (9) and (10), are used to generate four discrete methods and four additional methods. The discrete and additional methods are then applied as a BHFNM for solving (1). We choose J={r,s,v,1}={1/4,1/2,3/4,1}, and evaluating (9) at x=xn+j, j∈J, gives the four discrete methods yn+j=Π(xn+jh), which takes the form of the main methods (3). Evaluating (10) at x=xn+j gives the additional methods yn+j′=Π′(xn+jh), with j∈J, which takes the form of the additional methods (4). The coefficients and the corresponding Taylor series equivalence of yn+j, hyn+j′ where j∈J are, respectively, given in terms of u=wh as follows:(11)αr,0=1,γr,0=12βr,0=csc4u/8sec⁡u/85u2cos⁡u/8-9u2+128cos⁡3u/8+32usin⁡5u/8+4cos⁡5u/82048u2=36723040+199u26193152+1543u414863564800+281u6597939978240+7261781u82742391916199936000+31117u101914688174219591680+…βr,r=39u2-4utan⁡u/8-384+672-22u2csc2u/8+276ucot⁡u/8+3uu-18sin⁡u/4csc4u/8384u2=3128-337u23870720-457u41857945600-289u6290665267200-166697u831163544502272000-2138501u1065817405988798464000+…βr,s=csc4u/8576-11u2cos⁡u/2-7u2+576cos⁡u/4+144usin⁡u/21536u2=-473840+353u25160960+19u4165150720+569u63487983206400+72887u8457065319366656000-151u104178882919923712000+…βr,v=7u2+480csc2u/8+3ucsc4u/8sec⁡u/8ucos⁡3u/8-23sin⁡3u/8+sin⁡5u/8-384384u2=295760-u2241920+43u4464486400+127u6163499212800+12577u82448564210892800+20593u10632859672969216000+…βr,1=csc4u/8u2+384-7u2cos⁡u/4+48usin⁡3u/8sec⁡u/8-3843072u2=-77680-289u230965760-139u42123366400-2899u66975966412800-7116007u82742391916199936000-8563517u10526539247910387712000+…(12)αs,0=1,γs,0=12βs,0=csc4u/87u2-10u2cos⁡u/4+12sec⁡u/8usin⁡5u/8+2cos⁡5u/8-24384u2=531440+19u2241920+37u4154828800+163u6163499212800+461287u885699747381248000+12773u10391770273742848000+…βs,r=20u2+3u2csc4u/8-16u2-15csc2u/8-6utan⁡u/8+6ucot⁡u/823-9csc2u/8-28896u2=110-13u260480-67u4116121600-29u613624934400-46559u84284987369062400-26753u10411358787429990400+…βs,s=-csc4u/82u2cos⁡u/4+u2-24cos⁡u/2-12usin⁡u/2+2464u2=-148+u25760+23u477414400+u62477260800+331u81098714710016000-703u101371195958099968000+…βs,v=3csc2u/8-2-2u2+u2csc2u/8+48+6utan⁡u/8+6ucot⁡u/89-7csc2u/896u2=190-u260480+u45529600+13u68174960640+32027u83060705263616000+15019u10228532659683328000+…βs,1=-csc4u/8u2+2u2-24cos⁡u/4-12usin⁡3u/8sec⁡u/8+48384u2=-1480-u248384-13u492897280-47u654499737600-452681u885699747381248000-269639u108227175748599808000+…(13)αv,0=1,γv,0=34βv,0=csc4u/827u2-45u2cos⁡u/4+16sec⁡u/83usin⁡5u/8+4cos⁡5u/8-641024u2=1472560+141u21146880+67u4183500800+1157u6775107379200+8233u81025960312832000+316627u106500484542103552000+…βv,r=315u2+96csc2u/8+ucsc4u/8sec⁡u/83ucos⁡3u/8-2sin⁡3u/8+3sin⁡5u/8-128128u2=117640-3u28960-u41146880-19u66055526400-5647u8352673857536000-9811u10101570070970368000+…βv,s=-3csc4u/83u2-64cos⁡u/2+15u2+64cos⁡u/4-48usin⁡u/2512u2=271280+153u2573440+39u491750400+59u6129184563200-211u81880927240192000-107u1043336563614023680+…βv,v=9u2+9u2csc4u/8+96-18u2csc2u/8+12utan⁡u/8+12ucot⁡u/89-7csc2u/8-128128u2=3128-3u2143360+u43276800+7u62768240640+3257u8201527918592000+81163u10812560567762944000+…βv,1=3csc4u/8sec⁡u/88sin⁡3u/8-3ucos3u/8512u=-92560-39u21146880-41u4183500800-1039u6775107379200-18197u82257112688230400-321977u106500484542103552000+…(14)α1,0=1,γ1,0=1β1,0=csc4u/8sec⁡u/82usin⁡5u/8+cos⁡5u/8-u2+2cos⁡3u/832u2=790+u26048+u42073600+u6510935040+28319u82678117105664000+83u101285496210718720+…β1,r=csc4u/8sec⁡u/88u2cos⁡u/8+u2+12cos⁡3u/8+3u2-4cos⁡5u/8-6usin⁡3u/8-18usin⁡5u/896u2=415-u22160-17u414515200-u6243302400-56153u82678117105664000-32783u10257099242143744000+…β1,s=-csc4u/8-18sin⁡u/2+ucos⁡u/2+8ucos⁡u/448u=115+u22520+u41612800+u61703116800-97u8178541140377600-139u1028566582460416000+…β1,v=4u2+3u2csc4u/8-6u2+4csc2u/8+3utan⁡u/8+3ucot⁡u/89-7csc2u/8+4824u2=445-u215120+u42903040+17u65109350400+8299u8382588157952000+34451u10257099242143744000+…β1,1=-csc4u/8sec⁡u/82u2cos⁡u/8+u2-6cos⁡3u/8-6usin⁡3u/8+6cos⁡5u/896u2=-u230240-u43628800-u6567705600-7201u8669529276416000-17u10256842399744000+…(15)γr,0′=1βr,0′=csc4u/8sec⁡u/848sin⁡5u/8-sin⁡3u/8+11ucos⁡u/8-23ucos⁡3u/8768u=2512880+863u23870720+71u488473600+1409u6373712486400+1461541u868559797904998400+408887u103134162189942784000+…βr,r′=-24tan⁡u/8+u6csc4u/8-37csc2u/8+46+cot⁡u/8264-48csc2u/896u=3231440-1159u21935360-589u4309657600-10519u61307993702400-135067u83116354450227200-2869897u1010969567664799744000+…βr,s′=-csc4u/85ucos⁡u/4+4-9sin⁡u/2+9sin⁡u/4+ucos⁡u/296u=-11120+37u280640+23u425804800+41u627249868800+6049u82856658246041600+961u10457065319366656000+…βr,v′=10u+24tan⁡u/8-19ucsc2u/8+32u+sin⁡u/4-5sin⁡u/2csc4u/896u=531440-5u2387072+221u4309657600+1579u6261598740480+1388953u834279898952499200+567829u102193913532959948800+…βr,1′=-csc4u/8sec⁡u/848sin⁡u/8-sin⁡3u/8+7ucos⁡u/8+5ucos⁡3u/8768u=-192880-271u23870720-313u4619315200-8551u62615987404800-1413149u868559797904998400-258803u101994466848145408000+…(16)γs,0′=1βs,0′=csc4u/84u+3sin⁡u/2-3sin⁡u/4-7ucos⁡u/496u=29360+37u2241920+23u477414400+41u681749606400+6049u88569974738124800+961u101371195958099968000+…βs,r′=csc4u/811u-36sin⁡u/2+24sin⁡u/4+7ucos⁡u/2-6ucos⁡u/496u=3190-53u2120960-31u438707200-7u65839257600-5273u84284987369062400+151u10685597979049984000+…βs,s′=-csc4u/8-18sin⁡u/2+ucos⁡u/2+8ucos⁡u/448u=115+u22520+u41612800+u61703116800-97u8178541140377600-139u1028566582460416000+…βs,v′=cot2u/8csc2u/8ucos⁡u/4+2-12sin⁡u/424u=190-11u2120960-u438707200+17u640874803200+8377u84284987369062400+4297u10685597979049984000+…βs,1′=-csc4u/8ucos⁡u/4+2-12sin⁡u/496u=-1360-u248384-u411059200-u63269984256-691u8779088612556800-u10430516784332800+…(17)γv,0′=1βv,0′=csc4u/8sec⁡u/816sin⁡u/8+sin⁡5u/8-3ucos⁡u/8+3cos⁡3u/8256u=27320+29u2143360+7u49830400+1007u6290665267200+17293u8846417258086400+44623u10348240243326976000+…βv,r′=-8tan⁡u/8+3u2csc4u/8-5csc2u/8+6+cot⁡u/888-48csc2u/832u=51160-37u271680-53u434406400-991u6145332633600-1531u838473511731200-307553u101218840851644416000+…βv,s′=3csc4u/84sin⁡u/2+sin⁡u/4-3ucos⁡u/432u=940+3u28960+u42867200-u63027763200-113u835267385753600-601u1050785035485184000+…βv,v′=6u+8tan⁡u/8+6ucsc4u/8-9ucsc2u/8+8cot⁡u/85-6csc2u/832u=21160+u214336+37u434406400+211u629066526720+18649u8423208629043200+65357u10243768170328883200+…βv,1′=csc4u/832sin⁡3u/8sec⁡u/8-6ucos⁡u/4+1256u=-3320-13u2143360-41u468812800-1039u6290665267200-18197u8846417258086400-321977u102437681703288832000+…(18)γ1,0′=1β1,0′=csc4u/8u+6sin⁡u/2-4ucos⁡u/448u=790+u27560+u44838400+u65109350400-97u8535623421132800-139u1085699747381248000+…β1,r′=csc4u/8ucos⁡u/2+2-6sin⁡u/212u=1645-u21890-u41209600-u61277337600+97u8133905855283200+139u1021424936845312000+…β1,s′=-csc4u/8-18sin⁡u/2+ucos⁡u/2+8ucos⁡u/424u=215+u21260+u4806400+u6851558400-97u889270570188800-139u1014283291230208000+…β1,v′=csc4u/8ucos⁡u/2+2-6sin⁡u/212u=1645-u21890-u41209600-u61277337600+97u8133905855283200+139u1021424936845312000+…β1,1′=csc4u/8u+6sin⁡u/2-4ucos⁡u/448u=790+u27560+u44838400+u65109350400-97u8535623421132800-139u1085699747381248000+…Remark 2.
We note that the Taylor series expansions in (11) through (18) must be used when u→0 because the corresponding trigonometric coefficients given in these equations are vulnerable to heavy cancelations (see [8]).


## 2.1. Specification of the Method

The continuous form (8) and its first derivative, which are equivalent to the forms (9) and (10), are used to generate four discrete methods and four additional methods. The discrete and additional methods are then applied as a BHFNM for solving (1). We choose J={r,s,v,1}={1/4,1/2,3/4,1}, and evaluating (9) at x=xn+j, j∈J, gives the four discrete methods yn+j=Π(xn+jh), which takes the form of the main methods (3). Evaluating (10) at x=xn+j gives the additional methods yn+j′=Π′(xn+jh), with j∈J, which takes the form of the additional methods (4). The coefficients and the corresponding Taylor series equivalence of yn+j, hyn+j′ where j∈J are, respectively, given in terms of u=wh as follows:(11)αr,0=1,γr,0=12βr,0=csc4u/8sec⁡u/85u2cos⁡u/8-9u2+128cos⁡3u/8+32usin⁡5u/8+4cos⁡5u/82048u2=36723040+199u26193152+1543u414863564800+281u6597939978240+7261781u82742391916199936000+31117u101914688174219591680+…βr,r=39u2-4utan⁡u/8-384+672-22u2csc2u/8+276ucot⁡u/8+3uu-18sin⁡u/4csc4u/8384u2=3128-337u23870720-457u41857945600-289u6290665267200-166697u831163544502272000-2138501u1065817405988798464000+…βr,s=csc4u/8576-11u2cos⁡u/2-7u2+576cos⁡u/4+144usin⁡u/21536u2=-473840+353u25160960+19u4165150720+569u63487983206400+72887u8457065319366656000-151u104178882919923712000+…βr,v=7u2+480csc2u/8+3ucsc4u/8sec⁡u/8ucos⁡3u/8-23sin⁡3u/8+sin⁡5u/8-384384u2=295760-u2241920+43u4464486400+127u6163499212800+12577u82448564210892800+20593u10632859672969216000+…βr,1=csc4u/8u2+384-7u2cos⁡u/4+48usin⁡3u/8sec⁡u/8-3843072u2=-77680-289u230965760-139u42123366400-2899u66975966412800-7116007u82742391916199936000-8563517u10526539247910387712000+…(12)αs,0=1,γs,0=12βs,0=csc4u/87u2-10u2cos⁡u/4+12sec⁡u/8usin⁡5u/8+2cos⁡5u/8-24384u2=531440+19u2241920+37u4154828800+163u6163499212800+461287u885699747381248000+12773u10391770273742848000+…βs,r=20u2+3u2csc4u/8-16u2-15csc2u/8-6utan⁡u/8+6ucot⁡u/823-9csc2u/8-28896u2=110-13u260480-67u4116121600-29u613624934400-46559u84284987369062400-26753u10411358787429990400+…βs,s=-csc4u/82u2cos⁡u/4+u2-24cos⁡u/2-12usin⁡u/2+2464u2=-148+u25760+23u477414400+u62477260800+331u81098714710016000-703u101371195958099968000+…βs,v=3csc2u/8-2-2u2+u2csc2u/8+48+6utan⁡u/8+6ucot⁡u/89-7csc2u/896u2=190-u260480+u45529600+13u68174960640+32027u83060705263616000+15019u10228532659683328000+…βs,1=-csc4u/8u2+2u2-24cos⁡u/4-12usin⁡3u/8sec⁡u/8+48384u2=-1480-u248384-13u492897280-47u654499737600-452681u885699747381248000-269639u108227175748599808000+…(13)αv,0=1,γv,0=34βv,0=csc4u/827u2-45u2cos⁡u/4+16sec⁡u/83usin⁡5u/8+4cos⁡5u/8-641024u2=1472560+141u21146880+67u4183500800+1157u6775107379200+8233u81025960312832000+316627u106500484542103552000+…βv,r=315u2+96csc2u/8+ucsc4u/8sec⁡u/83ucos⁡3u/8-2sin⁡3u/8+3sin⁡5u/8-128128u2=117640-3u28960-u41146880-19u66055526400-5647u8352673857536000-9811u10101570070970368000+…βv,s=-3csc4u/83u2-64cos⁡u/2+15u2+64cos⁡u/4-48usin⁡u/2512u2=271280+153u2573440+39u491750400+59u6129184563200-211u81880927240192000-107u1043336563614023680+…βv,v=9u2+9u2csc4u/8+96-18u2csc2u/8+12utan⁡u/8+12ucot⁡u/89-7csc2u/8-128128u2=3128-3u2143360+u43276800+7u62768240640+3257u8201527918592000+81163u10812560567762944000+…βv,1=3csc4u/8sec⁡u/88sin⁡3u/8-3ucos3u/8512u=-92560-39u21146880-41u4183500800-1039u6775107379200-18197u82257112688230400-321977u106500484542103552000+…(14)α1,0=1,γ1,0=1β1,0=csc4u/8sec⁡u/82usin⁡5u/8+cos⁡5u/8-u2+2cos⁡3u/832u2=790+u26048+u42073600+u6510935040+28319u82678117105664000+83u101285496210718720+…β1,r=csc4u/8sec⁡u/88u2cos⁡u/8+u2+12cos⁡3u/8+3u2-4cos⁡5u/8-6usin⁡3u/8-18usin⁡5u/896u2=415-u22160-17u414515200-u6243302400-56153u82678117105664000-32783u10257099242143744000+…β1,s=-csc4u/8-18sin⁡u/2+ucos⁡u/2+8ucos⁡u/448u=115+u22520+u41612800+u61703116800-97u8178541140377600-139u1028566582460416000+…β1,v=4u2+3u2csc4u/8-6u2+4csc2u/8+3utan⁡u/8+3ucot⁡u/89-7csc2u/8+4824u2=445-u215120+u42903040+17u65109350400+8299u8382588157952000+34451u10257099242143744000+…β1,1=-csc4u/8sec⁡u/82u2cos⁡u/8+u2-6cos⁡3u/8-6usin⁡3u/8+6cos⁡5u/896u2=-u230240-u43628800-u6567705600-7201u8669529276416000-17u10256842399744000+…(15)γr,0′=1βr,0′=csc4u/8sec⁡u/848sin⁡5u/8-sin⁡3u/8+11ucos⁡u/8-23ucos⁡3u/8768u=2512880+863u23870720+71u488473600+1409u6373712486400+1461541u868559797904998400+408887u103134162189942784000+…βr,r′=-24tan⁡u/8+u6csc4u/8-37csc2u/8+46+cot⁡u/8264-48csc2u/896u=3231440-1159u21935360-589u4309657600-10519u61307993702400-135067u83116354450227200-2869897u1010969567664799744000+…βr,s′=-csc4u/85ucos⁡u/4+4-9sin⁡u/2+9sin⁡u/4+ucos⁡u/296u=-11120+37u280640+23u425804800+41u627249868800+6049u82856658246041600+961u10457065319366656000+…βr,v′=10u+24tan⁡u/8-19ucsc2u/8+32u+sin⁡u/4-5sin⁡u/2csc4u/896u=531440-5u2387072+221u4309657600+1579u6261598740480+1388953u834279898952499200+567829u102193913532959948800+…βr,1′=-csc4u/8sec⁡u/848sin⁡u/8-sin⁡3u/8+7ucos⁡u/8+5ucos⁡3u/8768u=-192880-271u23870720-313u4619315200-8551u62615987404800-1413149u868559797904998400-258803u101994466848145408000+…(16)γs,0′=1βs,0′=csc4u/84u+3sin⁡u/2-3sin⁡u/4-7ucos⁡u/496u=29360+37u2241920+23u477414400+41u681749606400+6049u88569974738124800+961u101371195958099968000+…βs,r′=csc4u/811u-36sin⁡u/2+24sin⁡u/4+7ucos⁡u/2-6ucos⁡u/496u=3190-53u2120960-31u438707200-7u65839257600-5273u84284987369062400+151u10685597979049984000+…βs,s′=-csc4u/8-18sin⁡u/2+ucos⁡u/2+8ucos⁡u/448u=115+u22520+u41612800+u61703116800-97u8178541140377600-139u1028566582460416000+…βs,v′=cot2u/8csc2u/8ucos⁡u/4+2-12sin⁡u/424u=190-11u2120960-u438707200+17u640874803200+8377u84284987369062400+4297u10685597979049984000+…βs,1′=-csc4u/8ucos⁡u/4+2-12sin⁡u/496u=-1360-u248384-u411059200-u63269984256-691u8779088612556800-u10430516784332800+…(17)γv,0′=1βv,0′=csc4u/8sec⁡u/816sin⁡u/8+sin⁡5u/8-3ucos⁡u/8+3cos⁡3u/8256u=27320+29u2143360+7u49830400+1007u6290665267200+17293u8846417258086400+44623u10348240243326976000+…βv,r′=-8tan⁡u/8+3u2csc4u/8-5csc2u/8+6+cot⁡u/888-48csc2u/832u=51160-37u271680-53u434406400-991u6145332633600-1531u838473511731200-307553u101218840851644416000+…βv,s′=3csc4u/84sin⁡u/2+sin⁡u/4-3ucos⁡u/432u=940+3u28960+u42867200-u63027763200-113u835267385753600-601u1050785035485184000+…βv,v′=6u+8tan⁡u/8+6ucsc4u/8-9ucsc2u/8+8cot⁡u/85-6csc2u/832u=21160+u214336+37u434406400+211u629066526720+18649u8423208629043200+65357u10243768170328883200+…βv,1′=csc4u/832sin⁡3u/8sec⁡u/8-6ucos⁡u/4+1256u=-3320-13u2143360-41u468812800-1039u6290665267200-18197u8846417258086400-321977u102437681703288832000+…(18)γ1,0′=1β1,0′=csc4u/8u+6sin⁡u/2-4ucos⁡u/448u=790+u27560+u44838400+u65109350400-97u8535623421132800-139u1085699747381248000+…β1,r′=csc4u/8ucos⁡u/2+2-6sin⁡u/212u=1645-u21890-u41209600-u61277337600+97u8133905855283200+139u1021424936845312000+…β1,s′=-csc4u/8-18sin⁡u/2+ucos⁡u/2+8ucos⁡u/424u=215+u21260+u4806400+u6851558400-97u889270570188800-139u1014283291230208000+…β1,v′=csc4u/8ucos⁡u/2+2-6sin⁡u/212u=1645-u21890-u41209600-u61277337600+97u8133905855283200+139u1021424936845312000+…β1,1′=csc4u/8u+6sin⁡u/2-4ucos⁡u/448u=790+u27560+u44838400+u65109350400-97u8535623421132800-139u1085699747381248000+…Remark 2.
We note that the Taylor series expansions in (11) through (18) must be used when u→0 because the corresponding trigonometric coefficients given in these equations are vulnerable to heavy cancelations (see [8]).


## 3. Properties of the Method

### 3.1. Local Truncation Error

We define the local truncation errors (LTEs) of (3) and (4) specified by the coefficients in (11) through (18) as(19)LJyxn;h=yxn+jh-αj,0yxn+hγj,0y′xn+h2∑k∈J0βj,ky′′xn+kh,LJ′yxn;h=hy′xn+jh-hγj,0′y′xn+h2∑k∈J0βj,k′y′′xn+kh.Assuming that y(x) is sufficiently differentiable, we can expand the terms in LJ as a Taylor series about the point xn to obtain the expressions for the LTEs LJ and LJ′ for J={r,s,v,1}={1/4,1/2,3/4,1} as(20)LJyxn;h=107h7165150720w2y5xn+y7xn+Oh8,h7645120w2y5xn+y7xn+Oh8,9h73670016w2y5xn+y7xn+Oh8,h7322560w2y5xn+y7xn+Oh8,(21)LJ′yxn;h=3h7655360w2y5xn+y7xn+Oh8,h7368640w2y5xn+y7xn+Oh8,3h7655360w2y5xn+y7xn+Oh8,-h81935360w2y6xn+y8xn+Oh9.Remark 3.
The local truncation error constants of the BHFNM formulated from(2) and (3) and specified by the coefficients in (11) through (18) are given, respectively, by (107/165150720, 1/645120, 9/3670016, 1/322560, 3/655360, 1/368640, 3/655360,0)⊤ and orders (5,5,5,5,5,5,5,6)T; hence the BHFNM has order at least 5.


### 3.2. Block Form

BHFNM is formulated from the four discrete hybrid formulas stated in (3) which are provided by the continuous one-step hybrid trigonometrically fitted method with three off-grid points given by (9) and its first derivative (10). We define the following vectors:(22)Yn+1=yn+r,yn+s,yn+v,yn+1,hyn+r′,hyn+s′,hyn+v′,hyn+1′⊤,Yn=yn-v,yn-s,yn-r,yn,hyn-v′,hyn-s′,hyn-r′,hyn′⊤,Fn+1=fn+r,fn+s,fn+v,fn+1,hfn+r′,hfn+s′,hfn+v′,hfn+1′⊤,Fn=fn-v,fn-s,fn-r,fn,hfn-v′,hfn-s′,hfn-r′,hfn-′⊤,where n=0,…,N. The methods in (3) specified by coefficients (11)-(18) are combined to give the BHFNM, which is expressed as(23)A1Yn+1=A0Yn+h2B0Fn+B1Fn+1where A0, A1, B0, and B1 are square matrices of dimension eight whose elements characterize the method and are given by the coefficients of (3).


### 3.3. Stability

The linear-stability of the BHFNM 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 (23) to the test equation yields(24)Yn+1=MΥ2;uYn,MΥ2;u≔A1-Υ2B1-1A0+Υ2B0,where the matrix M(Υ2;u) is the amplification matrix which determines the stability of the method.Definition 4.
A region of stability is a region in the(q,u)-plane, throughout which the spectral radius ρ(M(Υ2;u))≤1 and any closed curve given by ρ(M(Υ2;u))=1 define the stability boundary of the method (see [22]). We note that the plot for the stability region of the BHFNM method is given in Figure 1.Figure 1
The stability region plotted in the(Υ,u)-plane.Remark 5.
It is observed that in the(Υ,u)-plane the BHFNM is stable for Υ∈[0,75.09] and u∈[-π,π] (see Figure 1).


### 3.4. Implementation

The methods given by (11)-(18) are combined to form the block method BHFNM (23), which is used to solve (1) without requiring starting values and predictors. BHFNM is implemented in a block-by-block fashion using a Mathematica 10.0 code, enhanced by the features NSolve[] and FindRoot[] (see Keiper and Gear [32]) for linear and nonlinear problems, respectively. In particular, the BHFNM is applied over x0,x1, [x1,x2],…,[xN-1,xN] by choosing N,h=(b-a)/N, and using (23) to simultaneously obtain the values of (25)yn+1/4,yn+1/2,yn+3/4,yn+1,yn+1/4′,yn+1/2′,yn+3/4′,yn+1′⊤,n=0,1…,N-1.


## 3.1. Local Truncation Error

We define the local truncation errors (LTEs) of (3) and (4) specified by the coefficients in (11) through (18) as(19)LJyxn;h=yxn+jh-αj,0yxn+hγj,0y′xn+h2∑k∈J0βj,ky′′xn+kh,LJ′yxn;h=hy′xn+jh-hγj,0′y′xn+h2∑k∈J0βj,k′y′′xn+kh.Assuming that y(x) is sufficiently differentiable, we can expand the terms in LJ as a Taylor series about the point xn to obtain the expressions for the LTEs LJ and LJ′ for J={r,s,v,1}={1/4,1/2,3/4,1} as(20)LJyxn;h=107h7165150720w2y5xn+y7xn+Oh8,h7645120w2y5xn+y7xn+Oh8,9h73670016w2y5xn+y7xn+Oh8,h7322560w2y5xn+y7xn+Oh8,(21)LJ′yxn;h=3h7655360w2y5xn+y7xn+Oh8,h7368640w2y5xn+y7xn+Oh8,3h7655360w2y5xn+y7xn+Oh8,-h81935360w2y6xn+y8xn+Oh9.Remark 3.
The local truncation error constants of the BHFNM formulated from(2) and (3) and specified by the coefficients in (11) through (18) are given, respectively, by (107/165150720, 1/645120, 9/3670016, 1/322560, 3/655360, 1/368640, 3/655360,0)⊤ and orders (5,5,5,5,5,5,5,6)T; hence the BHFNM has order at least 5.


## 3.2. Block Form

BHFNM is formulated from the four discrete hybrid formulas stated in (3) which are provided by the continuous one-step hybrid trigonometrically fitted method with three off-grid points given by (9) and its first derivative (10). We define the following vectors:(22)Yn+1=yn+r,yn+s,yn+v,yn+1,hyn+r′,hyn+s′,hyn+v′,hyn+1′⊤,Yn=yn-v,yn-s,yn-r,yn,hyn-v′,hyn-s′,hyn-r′,hyn′⊤,Fn+1=fn+r,fn+s,fn+v,fn+1,hfn+r′,hfn+s′,hfn+v′,hfn+1′⊤,Fn=fn-v,fn-s,fn-r,fn,hfn-v′,hfn-s′,hfn-r′,hfn-′⊤,where n=0,…,N. The methods in (3) specified by coefficients (11)-(18) are combined to give the BHFNM, which is expressed as(23)A1Yn+1=A0Yn+h2B0Fn+B1Fn+1where A0, A1, B0, and B1 are square matrices of dimension eight whose elements characterize the method and are given by the coefficients of (3).


## 3.3. Stability

The linear-stability of the BHFNM 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 (23) to the test equation yields(24)Yn+1=MΥ2;uYn,MΥ2;u≔A1-Υ2B1-1A0+Υ2B0,where the matrix M(Υ2;u) is the amplification matrix which determines the stability of the method.Definition 4.
A region of stability is a region in the(q,u)-plane, throughout which the spectral radius ρ(M(Υ2;u))≤1 and any closed curve given by ρ(M(Υ2;u))=1 define the stability boundary of the method (see [22]). We note that the plot for the stability region of the BHFNM method is given in Figure 1.Figure 1
The stability region plotted in the(Υ,u)-plane.Remark 5.
It is observed that in the(Υ,u)-plane the BHFNM is stable for Υ∈[0,75.09] and u∈[-π,π] (see Figure 1).


## 3.4. Implementation

The methods given by (11)-(18) are combined to form the block method BHFNM (23), which is used to solve (1) without requiring starting values and predictors. BHFNM is implemented in a block-by-block fashion using a Mathematica 10.0 code, enhanced by the features NSolve[] and FindRoot[] (see Keiper and Gear [32]) for linear and nonlinear problems, respectively. In particular, the BHFNM is applied over x0,x1, [x1,x2],…,[xN-1,xN] by choosing N,h=(b-a)/N, and using (23) to simultaneously obtain the values of (25)yn+1/4,yn+1/2,yn+3/4,yn+1,yn+1/4′,yn+1/2′,yn+3/4′,yn+1′⊤,n=0,1…,N-1.


## 4. Numerical Examples

The comparison of the following methods are given in this section:(i)
Block Nyström method of order 5 based on polynomial basis and given in [33].(ii)
Block Hybrid Trigonometrically Fitted Algorithm of order 5 given in [30].(iii)
Block hybrid trigonometrically fitted Runge–Kutta–Nyström method of order 3 given in [31].(iv)
Block hybrid functionally fitted Runge–Kutta–Nyström method (BHFNM) of order 5 given in this paper.Example 6.
We consider the following general second-order IVP given in [34] (26)y′′+ω2y=-δy′,y0=1,y′0=-δ2,and the analytical solution is(27)yt=e-δ/2tcostω2-δ24,ω=1,δ=10-3.This example is solved using the Block Nyström method based on polynomial basis given in [33], Block Hybrid Trigonometrically Fitted Algorithm of order 5 given in [30], Block hybrid trigonometrically fitted Runge–Kutta–Nyström method of order 3 given in [31], and the BHFNM given in this paper. The errors obtained which are compared for different step sizes and displayed in Table 1 show that BHFNM is more accurate.Table 1
Comparison of global errors for Example6 on [0,100], w=1.
h BHFNM Method in [33] Method in [30] Method in [31] 2 6.30 × 1 0 - 6 6.46 × 1 0 - 3 2.74 × 1 0 - 4 6.48 × 1 0 - 4 1 1.09 × 1 0 - 7 1.17 × 1 0 - 4 6.31 × 1 0 - 6 4.39 × 1 0 - 5 1/2 1.53 × 1 0 - 9 1.88 × 1 0 - 6 1.16 × 1 0 - 7 2.99 × 1 0 - 6 1/4 2.90 × 1 0 - 11 2.96 × 1 0 - 8 1.85 × 1 0 - 9 1.88 × 1 0 - 7 1/8 5.16 × 1 0 - 13 4.64 × 1 0 - 10 2.92 × 1 0 - 11 1.18 × 1 0 - 8Example 7.
We consider the oscillatory system that is solved in [34](28)y′′+13-12-1213y=12ε2532-23y′+ε2f1tf2t,y0=ε2ε2,y′0=-46,with f1(t)=36/5sin⁡t+24sin⁡5t and f2(t)=-24/5sin⁡t-36sin⁡5t.
The exact solution of this system is(29)y=sin⁡t-sin⁡5t+εcos⁡tsin⁡t+sin⁡5t+εcos⁡t.This example is solved using the Block Nyström method based on polynomial basis given in [33], Block Hybrid Trigonometrically Fitted Algorithm of order 5 given in [30], Block hybrid trigonometrically fitted Runge–Kutta–Nyström method of order 3 given in [31], and the BHFNM given in this paper. The errors obtained are compared for different step sizes and displayed in Table 2 show that BHFNM is more accurate.Table 2
Comparison of global error for Example7 on the interval [0,100], w=5.
h BHFNM-Error(y1) Method in [33]-Error(y1) BHFNM-Error(y2) Method in [33]-Error(y2) Method in [30]-Error(y1)  Method in [30]-Error(y2) Method in [31]-Error(y1)  Method in [31]-Error(y2) 1 3.05 × 1 0 - 3 1.97 × 1 0 0 3.04 × 1 0 - 3 1.97 × 1 0 0 2.92 × 1 0 - 1  2.92 × 1 0 - 1 2.05 × 1 0 0  2.05 × 1 0 0  1/2 4.79 × 1 0 - 5 1.24 × 1 0 - 1 4.79 × 1 0 - 5 1.24 × 1 0 - 1 3.26 × 1 0 - 3  3.26 × 1 0 - 3 9.07 × 1 0 - 2  9.07 × 1 0 - 2  1/4 7.38 × 1 0 - 7 2.30 × 1 0 - 3 7.38 × 1 0 - 7 2.30 × 1 0 - 3 4.79 × 1 0 - 5  4.79 × 1 0 - 5 4.95 × 1 0 - 3  4.95 × 1 0 - 3  1/8 1.17 × 1 0 - 8 3.76 × 1 0 - 5 1.18 × 1 0 - 8 3.76 × 1 0 - 5 7.40 × 1 0 - 7  7.40 × 1 0 - 7 3.00 × 1 0 - 4  3.00 × 1 0 - 4  1/16 1.80 × 1 0 - 10 5.93 × 1 0 - 7 1.80 × 1 0 - 10 5.93 × 1 0 - 7 1.19 × 1 0 - 8  1.20 × 1 0 - 8 1.86 × 1 0 - 5  1.86 × 1 0 - 5  1/32 2.92 × 1 0 - 12 9.38 × 1 0 - 9 3.06 × 1 0 - 12 9.38 × 1 0 - 9 1.80 × 1 0 - 10  1.80 × 1 0 - 10 1.16 × 1 0 - 6  1.16 × 1 0 - 6


### 4.1. Application to a PDE

L ∞ and L2 are calculated as follows: L2=(∑iM‍ei,j2)1/2 and L∞=max⁡ei,j, 1≤i≤M, where ei,j=ui,j-uxi,tj, j=1,…,N.Example 8.
We consider the telegraph equation which was also solved in [33](30)∂2u∂t2+2π∂u∂t+π2u=∂2u∂x2+π2sin⁡πxsin⁡πt+2cos⁡πt,0≤x≤1,0≤t≤1.The analytical solution is given by(31)ux,t=sin⁡πxsin⁡πtand the initial conditions are defined to match with the analytical solution.This PDE is solved by first discretization the spatial variablex via the finite difference method to obtain(32)∂2uit∂t2=ui+1t-2uit+ui-1tΔx2-2π∂uit∂t-π2uit+git,0<t<1,where Δx=(x0-xM)/M, xi=x0+iΔx, i=0,1,…,M, ui(t)≈u(xi,t), U(t)=[u1(t),…,uM(t)]T, U′(t)=[u1′(t),…,uM′(t)]T, G=[g1(t)-2πu1′(t)-π2u1(t),…,gM(t)-2πuM′(t)-π2uM(t)]T, and gi(t)≈g(xi,t)=π2sin⁡πxi(sin⁡πt+cos⁡(πt)), which can be written in the form(33)U′′=ft,U,U′,subject to the initial conditions U(0)=0, i=1,…,M, U′(0)=πsin⁡(πxi), i=1,…,M, w=π, where f(t,U,U′)=SU+G, andS is a M×M, matrix arising from the semidiscretized system.This problem is solved using the BHFNM and the numerical results are displayed in Table3 and Figure 2. The table gives L∞ and L2 errors, while a comparison of the exact and numerical solutions is given in the figure plotted in 3 dimensions.Table 3
L ∞ and L2 errors for Example 8, for M=20 and N=10 at different times.
T L ∞ L 2 1 7.86979 × 1 0 - 4 2.48865 × 1 0 - 3 2 8.20935 × 1 0 - 4 2.59602 × 1 0 - 3 3 8.22400 × 1 0 - 4 2.60066 × 1 0 - 3 4 8.22463 × 1 0 - 4 2.60086 × 1 0 - 3 5 8.22467 × 1 0 - 4 2.60087 × 1 0 - 3Figure 2
3D plots of exact solution (a) and numerical solution (b) for Example8 up to t=2, M=20, and N=10.
(a)
Exact solution (b)
Numerical solutionExample 9.
We consider the dissipative nonlinear wave equation given in [35](34)∂2u∂t2=∂2u∂x2-2u∂u∂t+π2-1-2sin⁡πxsin⁡tsin⁡πxcos⁡πt,0≤x≤1,0≤t≤1,with initial conditions(35)ux,0=sin⁡πx,utx,0=0.
The analytical solution is given by(36)ux,t=sin⁡πxcos⁡tand the initial conditions are defined to match with the analytical solution.This PDE is solved by first discretization the spatial variablex via the finite difference method to obtain(37)∂2uit∂t2=ui+1t-2uit+ui-1tΔx2-2π∂uit∂t-π2uit+git,0<t<1,where Δx=(x0-xM)/M, xi=x0+iΔx, i=0,1,…,M, ui(t)≈u(xi,t), U(t)=[u1(t),…,uM(t)]T, U′(t)=[u1′(t),…,uM′(t)]T, G=[g1(t)-2u1(t)u1′(t),…,gM(t)-2uM(t)uM′(t)]T, and gi(t)≈g(xi,t)=(π2-1-2sin⁡πxisin⁡t)(sin⁡πxicos⁡(πt)), which can be written in the form(38)U′′=ft,U,U′,subject to the initial conditions U(0)=sin⁡(πxi), i=1,…,M, U′(0)=0, i=1,…,M, where f(t,U,U′)=SU+G, andS is a M×M,  w=π, matrix arising from the semidiscretized system.This problem is solved using the BHFNM and the numerical results are displayed in Table4 and Figure 3. The table gives L∞ and L2 errors, while a comparison of the exact and numerical solutions is given in the figure plotted in 3 dimensions.Table 4
L ∞ and L2 errors for Example 9, for M=N=10 at different times.
T L ∞ L 2 1 1.12220 × 1 0 - 2 1.74162 × 1 0 - 2 2 8.62553 × 1 0 - 3 1.46484 × 1 0 - 2 3 6.58683 × 1 0 - 3 1.52150 × 1 0 - 2 4 6.97169 × 1 0 - 3 1.46843 × 1 0 - 2 5 7.99197 × 1 0 - 3 1.74162 × 1 0 - 2Figure 3
3D plots of exact solution (a) and numerical solution (b) for Example9 up to t=2, M=20, and N=10.
(a)
Exact solution (b)
Numerical solutionExample 10.
We consider the dissipative nonlinear wave equation given in [35](39)∂2u∂t2=∂2u∂x2+Υu2-1∂u∂t+π2+Υ2e-2Υtsin2πxe-Υtsin⁡πx,0≤x≤1,0≤t≤1,with initial conditions(40)ux,0=sin⁡πx,utx,0=-Υsin⁡πx.
The analytical solution is given by(41)ux,t=e-Υtsin⁡πxand the initial conditions are defined to match with the analytical solution.This PDE is solved by first discretization the spatial variablex via the finite difference method to obtain(42)∂2uit∂t2=ui+1t-2uit+ui-1tΔx2+Υu2-1∂uit∂t+git,0<t<1,where Δx=(x0-xM)/M, xi=x0+iΔx, i=0,1,…,M, ui(t)≈u(xi,t), U(t)=[u1(t),…,uM(t)]T, U′(t)=[u1′(t),…,uM′(t)]T, G=[g1(t)-Υ(u1(t)2-1)u1′(t),…,gM(t)-Υ(uM(t)2-1)uM′(t)]T, and gi(t)≈g(xi,t)=π2+Υ2e-2Υtsin2πxie-Υtsin⁡(πxi), which can be written in the form(43)U′′=ft,U,U′,subject to the initial conditions U(0)=sin⁡(πxi), i=1,…,M, U′(0)=-Υsin⁡(πxi), i=1,…,M, where f(t,U,U′)=SU+G, andS is a M×M,w=π, matrix arising from the semidiscretized system.This problem is solved using the BHFNM and the numerical results are displayed in Table5 and Figure 4. The table gives L∞ and L2 errors, while a comparison of the exact and numerical solutions are given in the figure plotted in 3 dimensions.Table 5
L ∞ and L2 errors for Example 10, for M=N=10 and Υ=1 at different times.
T L ∞ L 2 1 8.26147 × 1 0 - 3 1.59350 × 1 0 - 2 2 7.51186 × 1 0 - 3 1.54297 × 1 0 - 2 3 3.78627 × 1 0 - 3 7.50972 × 1 0 - 2 4 2.87725 × 1 0 - 3 5.80621 × 1 0 - 2 5 3.78271 × 1 0 - 3 6.09587 × 1 0 - 2Figure 4
3D plots of exact solution (a) and numerical solution (b) for Example10 up to t=2, M=20, N=10, and Υ=2.
(a)
Exact solution (b)
Numerical solution


### 4.2. Estimating the Frequency

Some techniques for estimating the frequency are given in [24, 36–39]. In the spirit of [36], where the technique for estimating the frequency was to determine the roots of the polynomial associated with the LTE, we rewrite the principal terms given in (20) and (21) as follows:(44)ϖj=ϱjwn2y5xn+y7xn=ϱjD5wn2+D2yxn,ϱ8wn2y6xn+y8xn=ϱ8D6wn2+D2yxn.where ϱj, j=1,…,8, are constants, y(j)=djy/dxj, the positive j represents the jth derivative, and D=d/dx is a linear differential operator and therefore invertible, which facilitates the simplification of ϖj. We estimate wn by setting any member of ϖj=0 and solving for wn to obtain two possible solutions of wn as follows:(45)wn=±-D7yxnD5yxn=±-D2yxnyxn.This estimate has been tested on Example 6 on [0,100], for w=wn, where wn=ω2+δy′(xn)/y(xn).We observe that whenw=wn was used, an appropriate estimate for the parameter was provided at each step, in which case the computational time was increased and results were slightly less accurate than those given by the exact frequency, w=1. We then noticed that, for this example, the results produced by setting w=w0 were similar to those given by w=1. Therefore, for this example, the appropriate choice for the frequency is w=w0. We do not claim that this approach will work in all cases; nevertheless, we will rigorously conduct research on this topic in our future work. Details of the numerical results are given in Table 6.Table 6
Comparison of global errors for Example6 on [0,100], w=1, w=wn, and w=w0.
h BHFNM (w=1) BHFNM (w=wn) BHFNM (w=w0) 2 6.30 × 1 0 - 6 6.28 × 1 0 - 6 6.30 × 1 0 - 6 1 1.09 × 1 0 - 7 5.35 × 1 0 - 7 1.09 × 1 0 - 7 1/2 1.53 × 1 0 - 9 4.83 × 1 0 - 9 1.45 × 1 0 - 9 1/4 2.90 × 1 0 - 11 5.74 × 1 0 - 11 2.90 × 1 0 - 11 1/8 5.16 × 1 0 - 13 8.29 × 1 0 - 13 4.41 × 1 0 - 13


## 4.1. Application to a PDE

L ∞ and L2 are calculated as follows: L2=(∑iM‍ei,j2)1/2 and L∞=max⁡ei,j, 1≤i≤M, where ei,j=ui,j-uxi,tj, j=1,…,N.Example 8.
We consider the telegraph equation which was also solved in [33](30)∂2u∂t2+2π∂u∂t+π2u=∂2u∂x2+π2sin⁡πxsin⁡πt+2cos⁡πt,0≤x≤1,0≤t≤1.The analytical solution is given by(31)ux,t=sin⁡πxsin⁡πtand the initial conditions are defined to match with the analytical solution.This PDE is solved by first discretization the spatial variablex via the finite difference method to obtain(32)∂2uit∂t2=ui+1t-2uit+ui-1tΔx2-2π∂uit∂t-π2uit+git,0<t<1,where Δx=(x0-xM)/M, xi=x0+iΔx, i=0,1,…,M, ui(t)≈u(xi,t), U(t)=[u1(t),…,uM(t)]T, U′(t)=[u1′(t),…,uM′(t)]T, G=[g1(t)-2πu1′(t)-π2u1(t),…,gM(t)-2πuM′(t)-π2uM(t)]T, and gi(t)≈g(xi,t)=π2sin⁡πxi(sin⁡πt+cos⁡(πt)), which can be written in the form(33)U′′=ft,U,U′,subject to the initial conditions U(0)=0, i=1,…,M, U′(0)=πsin⁡(πxi), i=1,…,M, w=π, where f(t,U,U′)=SU+G, andS is a M×M, matrix arising from the semidiscretized system.This problem is solved using the BHFNM and the numerical results are displayed in Table3 and Figure 2. The table gives L∞ and L2 errors, while a comparison of the exact and numerical solutions is given in the figure plotted in 3 dimensions.Table 3
L ∞ and L2 errors for Example 8, for M=20 and N=10 at different times.
T L ∞ L 2 1 7.86979 × 1 0 - 4 2.48865 × 1 0 - 3 2 8.20935 × 1 0 - 4 2.59602 × 1 0 - 3 3 8.22400 × 1 0 - 4 2.60066 × 1 0 - 3 4 8.22463 × 1 0 - 4 2.60086 × 1 0 - 3 5 8.22467 × 1 0 - 4 2.60087 × 1 0 - 3Figure 2
3D plots of exact solution (a) and numerical solution (b) for Example8 up to t=2, M=20, and N=10.
(a)
Exact solution (b)
Numerical solutionExample 9.
We consider the dissipative nonlinear wave equation given in [35](34)∂2u∂t2=∂2u∂x2-2u∂u∂t+π2-1-2sin⁡πxsin⁡tsin⁡πxcos⁡πt,0≤x≤1,0≤t≤1,with initial conditions(35)ux,0=sin⁡πx,utx,0=0.
The analytical solution is given by(36)ux,t=sin⁡πxcos⁡tand the initial conditions are defined to match with the analytical solution.This PDE is solved by first discretization the spatial variablex via the finite difference method to obtain(37)∂2uit∂t2=ui+1t-2uit+ui-1tΔx2-2π∂uit∂t-π2uit+git,0<t<1,where Δx=(x0-xM)/M, xi=x0+iΔx, i=0,1,…,M, ui(t)≈u(xi,t), U(t)=[u1(t),…,uM(t)]T, U′(t)=[u1′(t),…,uM′(t)]T, G=[g1(t)-2u1(t)u1′(t),…,gM(t)-2uM(t)uM′(t)]T, and gi(t)≈g(xi,t)=(π2-1-2sin⁡πxisin⁡t)(sin⁡πxicos⁡(πt)), which can be written in the form(38)U′′=ft,U,U′,subject to the initial conditions U(0)=sin⁡(πxi), i=1,…,M, U′(0)=0, i=1,…,M, where f(t,U,U′)=SU+G, andS is a M×M,  w=π, matrix arising from the semidiscretized system.This problem is solved using the BHFNM and the numerical results are displayed in Table4 and Figure 3. The table gives L∞ and L2 errors, while a comparison of the exact and numerical solutions is given in the figure plotted in 3 dimensions.Table 4
L ∞ and L2 errors for Example 9, for M=N=10 at different times.
T L ∞ L 2 1 1.12220 × 1 0 - 2 1.74162 × 1 0 - 2 2 8.62553 × 1 0 - 3 1.46484 × 1 0 - 2 3 6.58683 × 1 0 - 3 1.52150 × 1 0 - 2 4 6.97169 × 1 0 - 3 1.46843 × 1 0 - 2 5 7.99197 × 1 0 - 3 1.74162 × 1 0 - 2Figure 3
3D plots of exact solution (a) and numerical solution (b) for Example9 up to t=2, M=20, and N=10.
(a)
Exact solution (b)
Numerical solutionExample 10.
We consider the dissipative nonlinear wave equation given in [35](39)∂2u∂t2=∂2u∂x2+Υu2-1∂u∂t+π2+Υ2e-2Υtsin2πxe-Υtsin⁡πx,0≤x≤1,0≤t≤1,with initial conditions(40)ux,0=sin⁡πx,utx,0=-Υsin⁡πx.
The analytical solution is given by(41)ux,t=e-Υtsin⁡πxand the initial conditions are defined to match with the analytical solution.This PDE is solved by first discretization the spatial variablex via the finite difference method to obtain(42)∂2uit∂t2=ui+1t-2uit+ui-1tΔx2+Υu2-1∂uit∂t+git,0<t<1,where Δx=(x0-xM)/M, xi=x0+iΔx, i=0,1,…,M, ui(t)≈u(xi,t), U(t)=[u1(t),…,uM(t)]T, U′(t)=[u1′(t),…,uM′(t)]T, G=[g1(t)-Υ(u1(t)2-1)u1′(t),…,gM(t)-Υ(uM(t)2-1)uM′(t)]T, and gi(t)≈g(xi,t)=π2+Υ2e-2Υtsin2πxie-Υtsin⁡(πxi), which can be written in the form(43)U′′=ft,U,U′,subject to the initial conditions U(0)=sin⁡(πxi), i=1,…,M, U′(0)=-Υsin⁡(πxi), i=1,…,M, where f(t,U,U′)=SU+G, andS is a M×M,w=π, matrix arising from the semidiscretized system.This problem is solved using the BHFNM and the numerical results are displayed in Table5 and Figure 4. The table gives L∞ and L2 errors, while a comparison of the exact and numerical solutions are given in the figure plotted in 3 dimensions.Table 5
L ∞ and L2 errors for Example 10, for M=N=10 and Υ=1 at different times.
T L ∞ L 2 1 8.26147 × 1 0 - 3 1.59350 × 1 0 - 2 2 7.51186 × 1 0 - 3 1.54297 × 1 0 - 2 3 3.78627 × 1 0 - 3 7.50972 × 1 0 - 2 4 2.87725 × 1 0 - 3 5.80621 × 1 0 - 2 5 3.78271 × 1 0 - 3 6.09587 × 1 0 - 2Figure 4
3D plots of exact solution (a) and numerical solution (b) for Example10 up to t=2, M=20, N=10, and Υ=2.
(a)
Exact solution (b)
Numerical solution


## 4.2. Estimating the Frequency

Some techniques for estimating the frequency are given in [24, 36–39]. In the spirit of [36], where the technique for estimating the frequency was to determine the roots of the polynomial associated with the LTE, we rewrite the principal terms given in (20) and (21) as follows:(44)ϖj=ϱjwn2y5xn+y7xn=ϱjD5wn2+D2yxn,ϱ8wn2y6xn+y8xn=ϱ8D6wn2+D2yxn.where ϱj, j=1,…,8, are constants, y(j)=djy/dxj, the positive j represents the jth derivative, and D=d/dx is a linear differential operator and therefore invertible, which facilitates the simplification of ϖj. We estimate wn by setting any member of ϖj=0 and solving for wn to obtain two possible solutions of wn as follows:(45)wn=±-D7yxnD5yxn=±-D2yxnyxn.This estimate has been tested on Example 6 on [0,100], for w=wn, where wn=ω2+δy′(xn)/y(xn).We observe that whenw=wn was used, an appropriate estimate for the parameter was provided at each step, in which case the computational time was increased and results were slightly less accurate than those given by the exact frequency, w=1. We then noticed that, for this example, the results produced by setting w=w0 were similar to those given by w=1. Therefore, for this example, the appropriate choice for the frequency is w=w0. We do not claim that this approach will work in all cases; nevertheless, we will rigorously conduct research on this topic in our future work. Details of the numerical results are given in Table 6.Table 6
Comparison of global errors for Example6 on [0,100], w=1, w=wn, and w=w0.
h BHFNM (w=1) BHFNM (w=wn) BHFNM (w=w0) 2 6.30 × 1 0 - 6 6.28 × 1 0 - 6 6.30 × 1 0 - 6 1 1.09 × 1 0 - 7 5.35 × 1 0 - 7 1.09 × 1 0 - 7 1/2 1.53 × 1 0 - 9 4.83 × 1 0 - 9 1.45 × 1 0 - 9 1/4 2.90 × 1 0 - 11 5.74 × 1 0 - 11 2.90 × 1 0 - 11 1/8 5.16 × 1 0 - 13 8.29 × 1 0 - 13 4.41 × 1 0 - 13


## 5. Conclusion

We have presented a BHFNM with coefficients that depend on a fixed frequency and the stepsize and that is implemented in a block-by-block fashion, which does not suffer from the disadvantages of requiring starting values and predictors inherent to predictor-corrector methods. It can be used to solve general second-order IVP, in particular, those that arise from the semidiscretization of hyperbolic PDE. Numerical studies have been presented to establish accuracy of the method. Since the exact value of the parameterw is not always available, we have attempted to estimate an appropriate value of the parameter w, in the spirit of [36], where the technique for estimating the frequency was to determine the roots of the polynomial associated with the LTE. We do not claim that this approach will work for all cases; hence we are determined to conduct rigorous research on this topic in our future work. We also note that the methods presented were not designed to incorporate a variable stepsize implementation, since we did not consider additional methods which are required to provide the error estimates as we proceed on the interval of interest. Nevertheless, our future research will include a variable stepsize implementation of the methods.



---

*Source: 1535430-2019-02-26.xml*