# Generalized Hybrid One-Step Block Method Involving Fifth Derivative for Solving Fourth-Order Ordinary Differential Equation Directly

Authors: Mohammad Alkasassbeh; Zurni Omar
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

A general one-step three-hybrid (off-step) points block method is proposed for solving fourth-order initial value problems of ordinary differential equations directly. A power series approximate function is employed for deriving this method. The approximate function is interpolated atxn,xn+r,xn+s,xn+t while its fourth and fifth derivatives are collocated at all points xi, i=0,r,s,t,1, in the interval of approximation. Several fourth-order initial value problems of ordinary differential equations are then solved to compare the performance of the proposed method with the derived methods. The analysis of the method reveals that the method is consistent and zero stable concluding that the method is also convergent. The numerical results demonstrate the superiority of the new method over the existing ones in terms of error.


---

## Body

## 1. Introduction

In this article, we consider the following general fourth-order initial value problems (IVPs) as shown:(1)yiv=fx,y,y′,y′′,y′′′,ya=α1,y0′a=α2,y0′′a=α3,y0′′′a=α4,x∈a,b, with the assumption that f∈C5a,b. Several phenomena in physical fields such as neural networks, electric circuits, and ship dynamics can be expressed in differential equations (DEs) forms (see [1–3]).Block method is one of the efficient methods proposed in 1953 by [4] to enhance the performance of the numerical methods. In 1967, [5] employed block method to provide the essential starting values needed for the predictor schemes. Subsequently, hybrid methods were initiated by [6], which involve the evaluation of functions off-step (nonstep) points. The introduction of hybrid (off-step) points in block methods has many advantages such as the ability to change step size, utilizing data off-step points, and the most important feature according to [7] is the capability to circumvent zero stability barrier condition ([8]).To increase the accuracy further and to solve the stiffness problem in ordinary differential equations (ODEs), [9] derived second derivative multistep methods for stiff ODEs. Recently, authors like in [10, 11] proposed high derivative methods for the same reason. The former developed block hybrid-second derivative method for stiff systems, while the latter introduced a Simpson’s type second derivative method for the solution of a first-order stiff system of IVP. In addition, [12] proposed a continuous fourth derivative method for third-order boundary value problems. Following these scholars’ footsteps a new generalized three-hybrid one-step fifth derivative method for solving fourth-order ODEs directly using the approach of interpolation and collocation will be proposed.This article consists of five sections: Section1 is for the introduction, while Section 2 illustrates the method derivation, where we consider three-off-step points and employing the collocation approach. The analysis of the approach is discussed in Section 3 which includes zero stability, order, consistency, and convergence. Section 4 covers the solution of some mathematical problems to show the performance of the developed method. Finally, a brief conclusion is performed in Section 5.


## 2. Development of the Method

Let the following power series polynomial be the approximate solution of (1)(2)ψx=∑j=02m+l-1ajx-xnhj,where l=4 and m=5, the number of interpolation and collocation points. Now, differentiating (2) four and five times yields(3)ψivx=∑j=42m+l-1ajj!h4j-4!x-xnhj-4=fx,y,y′,y′′,y′′′,ψvx=∑j=52m+l-1ajj!h5j-5!x-xnhj-5=gx,y,y′,y′′,y′′′.Interpolate (2) at xn+i=xn+ih, i=0,r,s,t, and collocate (3) at all points xn+i=xn+ih, i=0,r,s,t,1, where r,s,t∈0,1 produces a system of equations in matrix form as below(4)1000000⋯01rr2r3r4r5r6⋯rv^1ss2s3s4s5s6⋯sv^1tt2t3t4t5t6⋯tv^00004!0!h400⋯000004!0!h45!r1!h46!r22!h4⋯v^!rv^-4v^-4!h400004!0!h45!s1!h46!s22!h4⋯v^!sv^-4v^-4!h400004!0!h45!t1!h46!t22!h4⋯v^!tv^-4v^-4!h400004!0!h45!1!h46!2!h4⋯v^!v^-4!h4000005!0!h50⋯0000005!0!h56!r1!h5⋯v^!rv^-5v^-5!h5000005!0!h56!s1!h5⋯v^!sv^-5v^-5!h5000005!0!h56!t1!h5⋯v^!tv^-5v^-5!h5000005!0!h56!1!h5⋯v^!v^-5!h5a0a1a2a3a4a5a6a7a8a9a10a11a12a13=ynyn+ryn+syn+tfnfn+rfn+sfn+tfn+1gngn+rgn+sgn+tgn+1, where v^=2m+l-1. Matrix manipulation is then employed to solve the resulting system (4) for the unknown coefficients aj′s, j=0,1,…,13. Substituting the obtained values of the coefficients back into (2) yields(5)ψx=∑i=0,r,s,tαiyx+ih+h4∑i=r,s,tβifx+ih+∑i=01βifx+ih+h5∑i=r,s,tγigx+ih+∑i=01γigx+ih,where n=0,1,2,…,N-1, h=xn-xn-1 is the constant step size for the partition πN of the interval a,b which is given by πN=a=x0<x1<⋯<xN-1<xN=b, αi,βi, and γi are undetermined constants listed in Appendix I in the Supplementary Material (see Supplementary Material available online at https://doi.org/10.1155/2017/7637651). For simplicity we shall use yx+ih=yn+i, fx+ih=fn+i, gx+ih=gn+i, and gn+i=fn+i′.Calculating the first, second, and third derivatives of (5) produces(6)ψ′x=ddx∑i=0,r,s,tαiyn+i+h4ddx∑i=r,s,tβifn+i+∑i=01βifn+i+h5ddx∑i=r,s,tγign+i+∑i=01γign+i,(7)ψ′′x=d2dx2∑i=0,r,s,tαiyn+i+h4d2dx2∑i=r,s,tβifn+i+∑i=01βifn+i+h5d2dx2∑i=r,s,tγign+i+∑i=01γign+i,(8)ψ′′′x=d3dx3∑i=0,r,s,tαiyn+i+h4d3dx3∑i=r,s,tβifn+i+∑i=01βifn+i+h5d3dx3∑i=r,s,tγign+i+∑i=01γign+i. Evaluating (5) at the noninterpolating points xn+t,xn+1 with (6) and (7) at all points xn+i, i=0,r,s,t,1, gives the following general equations in block form:(9)A0Ym1=A1Ym0+∑i=01BiFmi+∑i=01DiGmi,where A0 is an identity matrix of order 12 and (10)A1=0001000rh000r2h22000r3h360001000sh000s2h22000s3h360001000th000t2h22000t3h360001000h000h22000h3600000001000rh000r2h2200000001000sh000s2h2200000001000th000t2h2200000001000h000h22000000000001000rh000000000001000sh000000000001000th000000000001000h0000000000000001000000000000000100000000000000010000000000000001,B0=000000000000000B1160000000000000000B2160000000000000000B3160000000000000000B4160000000000000000B5160000000000000000B6160000000000000000B7160000000000000000B8160000000000000000B9160000000000000000B10160000000000000000B11160000000000000000B12160000000000000000B13160000000000000000B14160000000000000000B15160000000000000000B16160,B1=B111B121B131B141B211B221B231B241B311B321B331B341B411B421B431B441B511B521B531B541B611B621B631B641B711B721B731B741B811B821B831B841B911B921B931B941B1011B1021B1031B1041B1111B1121B1131B1141B1211B1221B1231B1241B1311B1321B1331B1341B1411B1421B1431B1441B1511B1521B1531B1541B1611B1621B1631B1641,D0=000000000000000D1160000000000000000D2160000000000000000D3160000000000000000D4160000000000000000D5160000000000000000D6160000000000000000D7160000000000000000D8160000000000000000D9160000000000000000D10160000000000000000D11160000000000000000D12160000000000000000D13160000000000000000D14160000000000000000D15160000000000000000D16160,D1=D111D121D131D141D211D221D231D241D311D321D331D341D411D421D431D441D511D521D531D541D611D621D631D641D711D721D731D741D811D821D831D841D911D921D931D941D1011D1021D1031D1041D1111D1121D1131D1141D1211D1221D1231D1241D1311D1321D1331D1341D1411D1421D1431D1441D1511D1521D1531D1541D1611D1621D1631D1641.The entries of B0, B1, D0, and D1 are listed in Appendix II in the Supplementary Material section, while the vectors Ym0, Ym1, Fm0, Fm1, Gm0, Gm1 are defined as follows: (11)Ym0=yn-tyn-syn-rynyn-t′yn-s′yn-r′yn′yn-t′′yn-s′′yn-r′′yn′′yn-t′′′yn-s′′′yn-r′′′yn′′′,Ym1=yn+ryn+syn+tyn+1yn+r′yn+s′yn+t′yn+1′yn+r′′yn+s′′yn+t′′yn+1′′yn+r′′′yn+s′′′yn+t′′′yn+1′′′,Fm0=fn-15fn-14fn-13fn-12fn-11fn-10fn-9fn-8fn-7fn-6fn-5fn-4fn-3fn-2fn-1fn,Fm1=fn+rfn+sfn+tfn+1,Gm0=gn-15gn-14gn-13gn-12gn-11gn-10gn-9gn-8gn-7gn-6gn-5gn-4gn-3gn-2gn-1gn,Gm1=gn+rgn+sgn+tgn+1.


## 3. Analysis of the Method

### 3.1. Zero Stability

Definition 1.
The hybrid block method formula (9) is said to be zero stable if no root Rm of the first characteristic equation ρR has modulus greater than one; that is, Rm⩽1, and if Rm=1 then the multiplicity of Rm must not exceed four.To prove that the roots of the first characteristic equation satisfy the previous definition(12)ρR=detRA0-A1=R00-1000-rh000-r2h22000-r3h360R0-1000-sh000-s2h22000-s3h3600R-1000-th000-t2h22000-t3h36000R-1000-h000-h22000-h360000R00-1000-rh000-r2h2200000R0-1000-sh000-s2h22000000R-1000-th000-t2h220000000R-1000-h000-h2200000000R00-1000-rh000000000R0-1000-sh0000000000R-1000-th00000000000R-1000-h000000000000R00-10000000000000R0-100000000000000R-1000000000000000R-1=0,which imply that (13)R12R-14=0⟹Rm=0,ifm=1112,1,ifm=13116.As a result, the developed method is zero stable.


### 3.2. Order of the Method

The linear operatorΘ associated with the hybrid block methods formula (9) is defined as (14)Θyx;h=A0Ym1-A1Ym0-∑i=01BiFmi-∑i=01DiGmi. Expanding the above equation in Taylor series and combining like terms imply(15)Θyx;h=C^0h0yx+C^1h1y′x+C^2h2y′′x+⋯+C^p+4hp+4yp+4x+⋯. According to [7, 13] method (9) is said to be of order p if (16)C^0=C^1=⋯=C^p+3=0,C^p+4≠0.The term C^p+4 is called the error constant and the local truncation error is given by (17)tn+k=C^p+4yp+4hp+4+Ohp+5. Comparing like terms of yi and hi in (15) produces the coefficients C^0=C^1=⋯=C^13=0 with vector of error constants (18)C^14=C^141C^142C^143C^144C^145C^146C^147C^148C^149C^1410C^1411C^1412C^1413C^1414C^1415C^1416T, where(19)C^141=r8784604620800020r6-70r5s-70r5t-70r5+65r4s2+260r4st+260r4s+65r4t2+260r4t+65r4-260r3s2t-260r3s2-260r3st2-1040r3st-260r3s-260r3t2-260r3t+286r2s2t2+1144r2s2t+286r2s2+1144r2st2+1144r2st+286r2t2-1430rs2t2-1430rs2t-1430rst2+2145s2t2,C^142=s8784604620800065r2s4-260r2s3t-260r2s3+286r2s2t2+1144r2s2t+286r2s2-1430r2st2-1430r2st+2145r2t2-70rs5+260rs4t+260rs4-260rs3t2-1040rs3t-260rs3+1144rs2t2+1144rs2t-1430rst2+20s6-70s5t-70s5+65s4t2+260s4t+65s4-260s3t2-260s3t+286s2t2,C^143=t87846046208000286r2s2t2-1430r2s2t+2145r2s2-260r2st3+1144r2st2-1430r2st+65r2t4-260r2t3+286r2t2-260rs2t3+1144rs2t2-1430rs2t+260rst4-1040rst3+1144rst2-70rt5+260rt4-260rt3+65s2t4-260s2t3+286s2t2-70st5+260st4-260st3+20t6-70t5+65t4,C^144=178460462080002145r2s2t2-1430r2s2t+286r2s2-1430r2st2+1144r2st-260r2s+286r2t2-260r2t+65r2-1430rs2t2+1144rs2t-260rs2+1144rst2-1040rst+260rs-260rt2+260rt-70r+286s2t2-260s2t+65s2-260st2+260st-70s+65t2-70t+20,C^145=r7130767436800028r6-91r5s-91r5t-91r5+78r4s2+312r4st+312r4s+78r4t2+312r4t+78r4-286r3s2t-286r3s2-286r3st2-1144r3st-286r3s-286r3t2-286r3t+286r2s2t2+1144r2s2t+286r2s2+1144r2st2+1144r2st+286r2t2-1287rs2t2-1287rs2t-1287rst2+1716s2t2,C^146=s7130767436800078r2s4-286r2s3t-286r2s3+286r2s2t2+1144r2s2t+286r2s2-1287r2st2-1287r2st+1716r2t2-91rs5+312rs4t+312rs4-286rs3t2-1144rs3t-286rs3+1144rs2t2+1144rs2t-1287rst2+28s6-91s5t-91s5+78s4t2+312s4t+78s4-286s3t2-286s3t+286s2t2,C^147=t71307674368000286r2s2t2-1287r2s2t+1716r2s2-286r2st3+1144r2st2-1287r2st+78r2t4-286r2t3+286r2t2-286rs2t3+1144rs2t2-1287rs2t+312rst4-1144rst3+1144rst2-91rt5+312rt4-286rt3+78s2t4-286s2t3+286s2t2-91st5+312st4-286st3+28t6-91t5+78t4,C^148=113076743680001716r2s2t2-1287r2s2t+286r2s2-1287r2st2+1144r2st-286r2s+286r2t2-286r2t+78r2-1287rs2t2+1144rs2t-286rs2+1144rst2-1144rst+312rs-286rt2+312rt-91r+286s2t2-286s2t+78s2-286st2+312st-91s+78t2-91t+28,C^149=r610059033600014r6-42r5s-42r5t-42r5+33r4s2+132r4st+132r4s+33r4t2+132r4t+33r4-110r3s2t-110r3s2-110r3st2-440r3st-110r3s-110r3t2-110r3t+99r2s2t2+396r2s2t+99r2s2+396r2st2+396r2st+99r2t2-396rs2t2-396rs2t-396rst2+462s2t2,C^1410=s610059033600033r2s4-110r2s3t-110r2s3+99r2s2t2+396r2s2t+99r2s2-396r2st2-396r2st+462r2t2-42rs5+132rs4t+132rs4-110rs3t2-440rs3t-110rs3+396rs2t2+396rs2t-396rst2+14s6-42s5t-42s5+33s4t2+132s4t+33s4-110s3t2-110s3t+99s2t2,C^1411=t610059033600099r2s2t2-396r2s2t+462r2s2-110r2st3+396r2st2-396r2st+33r2t4-110r2t3+99r2t2-110rs2t3+396rs2t2-396rs2t+132rst4-440rst3+396rst2-42rt5+132rt4-110rt3+33s2t4-110s2t3+99s2t2-42st5+132st4-110st3+14t6-42t5+33t4,C^1412=1100590336000462r2s2t2-396r2s2t+99r2s2-396r2st2+396r2st-110r2s+99r2t2-110r2t+33r2-396rs2t2+396rs2t-110rs2+396rst2-440rst+132rs-110rt2+132rt-42r+99s2t2-110s2t+33s2-110st2+132st-42s+33t2-42t+14,C^1413=r55029516800028r6-77r5s-77r5t-77r5+55r4s2+220r4st+220r4s+55r4t2+220r4t+55r4-165r3s2t-165r3s2-165r3st2-660r3st-165r3s-165r3t2-165r3t+132r2s2t2+528r2s2t+132r2s2+528r2st2+528r2st+132r2t2-462rs2t2-462rs2t-462rst2+462s2t2,C^1414=s55029516800055r2s4-165r2s3t-165r2s3+132r2s2t2+528r2s2t+132r2s2-462r2st2-462r2st+462r2t2-77rs5+220rs4t+220rs4-165rs3t2-660rs3t-165rs3+528rs2t2+528rs2t-462rst2+28s6-77s5t-77s5+55s4t2+220s4t+55s4-165s3t2-165s3t+132s2t2,C^1415=t550295168000132r2s2t2-462r2s2t+462r2s2-165r2st3+528r2st2-462r2st+55r2t4-165r2t3+132r2t2-165rs2t3+528rs2t2-462rs2t+220rst4-660rst3+528rst2-77rt5+220rt4-165rt3+55s2t4-165s2t3+132s2t2-77st5+220st4-165st3+28t6-77t5+55t4,C^1416=150295168000462r2s2t2-462r2s2t+132r2s2-462r2st2+528r2st-165r2s+132r2t2-165r2t+55r2-462rs2t2+528rs2t-165rs2+528rst2-660rst+220rs-165rt2+220rt-77r+132s2t2-165s2t+55s2-165st2+220st-77s+55t2-77t+28, which concludes that the order of the developed method is p=10.


### 3.3. Consistency

Definition 2.
A block method is said to be consistent if its orderp is greater than one.Consistency property is achieved for the hybrid block method from the above analysis since the orderp=10≥1.


### 3.4. Convergence

Theorem 3 (see [14]).
Consistency and zero stability are sufficient conditions for a linear multistep method to be convergent.The hybrid block method (9) is convergent since it fulfils both the consistency and zero stability conditions.


### 3.5. Region of Absolute Stability

Stability region of the hybrid block (9) is discussed in spirit of [7]. The test problem of the forms y′=λy, y′′=λ2y, y′′′=λ3y, yiv=λ4y, and yv=λ5y, where λ<0 are substituted into the main methods of block (9), is shown below(20)A¯0Y¯m+1=∑i=14A¯iY¯mi+∑i=01B¯iF¯mi+∑i=01D¯iG¯mi,where (21)A¯0=1000010000100001,A¯1=0001000100010001,A¯2=000r000s000t0001,A¯3=000r22000s22000t2200012,A¯4=000r36000s36000t3600016,B¯0=1h4000B1160000B2160000B3160000B4160,B¯1=1h4B111B121B131B141B211B221B231B241B311B321B331B341B411B421B431B441,D¯0=1h5000D1160000D2160000D3160000D4160,D¯1=1h5D111D121D131D141D211D221D231D241D311D321D331D341D411D421D431D441,while the vectors Y¯m+1, Y¯m1, Y¯m2, Y¯m3, F¯m0, F¯m1, G¯m0, and G¯m1 are defined as follows: (22)Y¯m+1=yn+ryn+syn+tyn+1,Y¯m1=yn-tyn-syn-ryn,Y¯m2=yn-t′yn-s′yn-r′yn′,Y¯m3=yn-t′′yn-s′′yn-r′′yn′′,Y¯m4=yn-t′′′yn-s′′′yn-r′′′yn′′′,F¯m0=fn-3fn-2fn-1fn,F¯m1=fn+rfn+sfn+tfn+1,G¯m0=gn-3gn-2gn-1gn,G¯m1=gn+rgn+sgn+tgn+1, substituting q=λh in (20) yields(23)Y¯m+1=MqY¯m1,where the amplification matrix Mq is given by(24)Mq=A¯0-q4B¯1-q5D¯1-1A¯1+qA¯2+q2A¯3+q3A¯4+q4B¯0+q5D¯0.The analysis shows that the matrix Mq has eigenvalues η1,η2,η3,η4=0,0,0,η4, where the dominant eigenvalue η4 is a function of q subject to the values r,s,t=1/4,1/2,3/4 given by(25)η4=∑i=020ciqiK∑j=021cjqj, where K=10376293541461622784 and the values ci and cj are listed in Table 1.Table 1
c-value
q
i coefficients
q
j coefficients
c
0
−19991259433786613699768237287038031044280320
0
c
1
−19991259433786613699768237287038031044280320
192662816967402505371648000
c
2
−99956297168933068498841186435190155221401600
0
c
3
−33318765722977689499613728811730051740467200
0
c
4
−78083522390790487043756727396404132472422400
0
c
5
−12270936643443632883203625005043197725900800
502432963738164264960000
c
6
−994813550529796720162894578620563036569600
−79502926099016318976000
c
7
59773830578760679822424316692208078028800
0
c
8
31534047193242772464338416945297159618560
0
c
9
4794949827573668397142419401567455150080
133171866570211983360
c
10
359049264395827172053492731789485015040
−14458011284982988800
c
11
−3356754007384253296917217569574748160
268452685173657600
c
12
−4527291309580589178991255608586403840
0
c
13
−632549152682975031987854772642562048
5775969118252032
c
14
−50983394698341254755895386446446592
−659482112163840
c
15
−2347009723724983385107341899462400
26223671235440
c
16
3581359970953878842628978864128
−529590387200
c
17
8897056249188627790299370721984
67083520896
c
18
741826115832694210392858425664
−8638140672
c
19
30585824514148316274163930311
506275974
c
20
582915589719914249658541800
−15909210
c
21
0
189567The region of absolute stability is depicted in the dark area in Figure1.Figure 1
Region of absolute stability.


### 3.6. Pseudocode of the Developed Method

The following pseudocode illustrates how the developed method is implemented for solving fourth-order initial value problems.Step 1.
The valuesyx0,y′x0,y′′x0,y′′′x0 on the interval a,b from the problem are given.Step 2.
Set the values of the hybrid pointsr,s,t and step size h.Step 3.
Forj=a:h:b, set xn+j=xn+jh, j=r,s,t,1.Step 4.
Evaluate the approximate valuesyxn+j, y′xn+j, y′′xn+j, y′′′xn+j using the direct hybrid block method.Step 5.
Find the solution for the resultant system using the built in function fsolve in Matlab.Step 6.
Calculate the maximum error =Exactsolution-approximatesolution.


## 3.1. Zero Stability

Definition 1.
The hybrid block method formula (9) is said to be zero stable if no root Rm of the first characteristic equation ρR has modulus greater than one; that is, Rm⩽1, and if Rm=1 then the multiplicity of Rm must not exceed four.To prove that the roots of the first characteristic equation satisfy the previous definition(12)ρR=detRA0-A1=R00-1000-rh000-r2h22000-r3h360R0-1000-sh000-s2h22000-s3h3600R-1000-th000-t2h22000-t3h36000R-1000-h000-h22000-h360000R00-1000-rh000-r2h2200000R0-1000-sh000-s2h22000000R-1000-th000-t2h220000000R-1000-h000-h2200000000R00-1000-rh000000000R0-1000-sh0000000000R-1000-th00000000000R-1000-h000000000000R00-10000000000000R0-100000000000000R-1000000000000000R-1=0,which imply that (13)R12R-14=0⟹Rm=0,ifm=1112,1,ifm=13116.As a result, the developed method is zero stable.


## 3.2. Order of the Method

The linear operatorΘ associated with the hybrid block methods formula (9) is defined as (14)Θyx;h=A0Ym1-A1Ym0-∑i=01BiFmi-∑i=01DiGmi. Expanding the above equation in Taylor series and combining like terms imply(15)Θyx;h=C^0h0yx+C^1h1y′x+C^2h2y′′x+⋯+C^p+4hp+4yp+4x+⋯. According to [7, 13] method (9) is said to be of order p if (16)C^0=C^1=⋯=C^p+3=0,C^p+4≠0.The term C^p+4 is called the error constant and the local truncation error is given by (17)tn+k=C^p+4yp+4hp+4+Ohp+5. Comparing like terms of yi and hi in (15) produces the coefficients C^0=C^1=⋯=C^13=0 with vector of error constants (18)C^14=C^141C^142C^143C^144C^145C^146C^147C^148C^149C^1410C^1411C^1412C^1413C^1414C^1415C^1416T, where(19)C^141=r8784604620800020r6-70r5s-70r5t-70r5+65r4s2+260r4st+260r4s+65r4t2+260r4t+65r4-260r3s2t-260r3s2-260r3st2-1040r3st-260r3s-260r3t2-260r3t+286r2s2t2+1144r2s2t+286r2s2+1144r2st2+1144r2st+286r2t2-1430rs2t2-1430rs2t-1430rst2+2145s2t2,C^142=s8784604620800065r2s4-260r2s3t-260r2s3+286r2s2t2+1144r2s2t+286r2s2-1430r2st2-1430r2st+2145r2t2-70rs5+260rs4t+260rs4-260rs3t2-1040rs3t-260rs3+1144rs2t2+1144rs2t-1430rst2+20s6-70s5t-70s5+65s4t2+260s4t+65s4-260s3t2-260s3t+286s2t2,C^143=t87846046208000286r2s2t2-1430r2s2t+2145r2s2-260r2st3+1144r2st2-1430r2st+65r2t4-260r2t3+286r2t2-260rs2t3+1144rs2t2-1430rs2t+260rst4-1040rst3+1144rst2-70rt5+260rt4-260rt3+65s2t4-260s2t3+286s2t2-70st5+260st4-260st3+20t6-70t5+65t4,C^144=178460462080002145r2s2t2-1430r2s2t+286r2s2-1430r2st2+1144r2st-260r2s+286r2t2-260r2t+65r2-1430rs2t2+1144rs2t-260rs2+1144rst2-1040rst+260rs-260rt2+260rt-70r+286s2t2-260s2t+65s2-260st2+260st-70s+65t2-70t+20,C^145=r7130767436800028r6-91r5s-91r5t-91r5+78r4s2+312r4st+312r4s+78r4t2+312r4t+78r4-286r3s2t-286r3s2-286r3st2-1144r3st-286r3s-286r3t2-286r3t+286r2s2t2+1144r2s2t+286r2s2+1144r2st2+1144r2st+286r2t2-1287rs2t2-1287rs2t-1287rst2+1716s2t2,C^146=s7130767436800078r2s4-286r2s3t-286r2s3+286r2s2t2+1144r2s2t+286r2s2-1287r2st2-1287r2st+1716r2t2-91rs5+312rs4t+312rs4-286rs3t2-1144rs3t-286rs3+1144rs2t2+1144rs2t-1287rst2+28s6-91s5t-91s5+78s4t2+312s4t+78s4-286s3t2-286s3t+286s2t2,C^147=t71307674368000286r2s2t2-1287r2s2t+1716r2s2-286r2st3+1144r2st2-1287r2st+78r2t4-286r2t3+286r2t2-286rs2t3+1144rs2t2-1287rs2t+312rst4-1144rst3+1144rst2-91rt5+312rt4-286rt3+78s2t4-286s2t3+286s2t2-91st5+312st4-286st3+28t6-91t5+78t4,C^148=113076743680001716r2s2t2-1287r2s2t+286r2s2-1287r2st2+1144r2st-286r2s+286r2t2-286r2t+78r2-1287rs2t2+1144rs2t-286rs2+1144rst2-1144rst+312rs-286rt2+312rt-91r+286s2t2-286s2t+78s2-286st2+312st-91s+78t2-91t+28,C^149=r610059033600014r6-42r5s-42r5t-42r5+33r4s2+132r4st+132r4s+33r4t2+132r4t+33r4-110r3s2t-110r3s2-110r3st2-440r3st-110r3s-110r3t2-110r3t+99r2s2t2+396r2s2t+99r2s2+396r2st2+396r2st+99r2t2-396rs2t2-396rs2t-396rst2+462s2t2,C^1410=s610059033600033r2s4-110r2s3t-110r2s3+99r2s2t2+396r2s2t+99r2s2-396r2st2-396r2st+462r2t2-42rs5+132rs4t+132rs4-110rs3t2-440rs3t-110rs3+396rs2t2+396rs2t-396rst2+14s6-42s5t-42s5+33s4t2+132s4t+33s4-110s3t2-110s3t+99s2t2,C^1411=t610059033600099r2s2t2-396r2s2t+462r2s2-110r2st3+396r2st2-396r2st+33r2t4-110r2t3+99r2t2-110rs2t3+396rs2t2-396rs2t+132rst4-440rst3+396rst2-42rt5+132rt4-110rt3+33s2t4-110s2t3+99s2t2-42st5+132st4-110st3+14t6-42t5+33t4,C^1412=1100590336000462r2s2t2-396r2s2t+99r2s2-396r2st2+396r2st-110r2s+99r2t2-110r2t+33r2-396rs2t2+396rs2t-110rs2+396rst2-440rst+132rs-110rt2+132rt-42r+99s2t2-110s2t+33s2-110st2+132st-42s+33t2-42t+14,C^1413=r55029516800028r6-77r5s-77r5t-77r5+55r4s2+220r4st+220r4s+55r4t2+220r4t+55r4-165r3s2t-165r3s2-165r3st2-660r3st-165r3s-165r3t2-165r3t+132r2s2t2+528r2s2t+132r2s2+528r2st2+528r2st+132r2t2-462rs2t2-462rs2t-462rst2+462s2t2,C^1414=s55029516800055r2s4-165r2s3t-165r2s3+132r2s2t2+528r2s2t+132r2s2-462r2st2-462r2st+462r2t2-77rs5+220rs4t+220rs4-165rs3t2-660rs3t-165rs3+528rs2t2+528rs2t-462rst2+28s6-77s5t-77s5+55s4t2+220s4t+55s4-165s3t2-165s3t+132s2t2,C^1415=t550295168000132r2s2t2-462r2s2t+462r2s2-165r2st3+528r2st2-462r2st+55r2t4-165r2t3+132r2t2-165rs2t3+528rs2t2-462rs2t+220rst4-660rst3+528rst2-77rt5+220rt4-165rt3+55s2t4-165s2t3+132s2t2-77st5+220st4-165st3+28t6-77t5+55t4,C^1416=150295168000462r2s2t2-462r2s2t+132r2s2-462r2st2+528r2st-165r2s+132r2t2-165r2t+55r2-462rs2t2+528rs2t-165rs2+528rst2-660rst+220rs-165rt2+220rt-77r+132s2t2-165s2t+55s2-165st2+220st-77s+55t2-77t+28, which concludes that the order of the developed method is p=10.


## 3.3. Consistency

Definition 2.
A block method is said to be consistent if its orderp is greater than one.Consistency property is achieved for the hybrid block method from the above analysis since the orderp=10≥1.


## 3.4. Convergence

Theorem 3 (see [14]).
Consistency and zero stability are sufficient conditions for a linear multistep method to be convergent.The hybrid block method (9) is convergent since it fulfils both the consistency and zero stability conditions.


## 3.5. Region of Absolute Stability

Stability region of the hybrid block (9) is discussed in spirit of [7]. The test problem of the forms y′=λy, y′′=λ2y, y′′′=λ3y, yiv=λ4y, and yv=λ5y, where λ<0 are substituted into the main methods of block (9), is shown below(20)A¯0Y¯m+1=∑i=14A¯iY¯mi+∑i=01B¯iF¯mi+∑i=01D¯iG¯mi,where (21)A¯0=1000010000100001,A¯1=0001000100010001,A¯2=000r000s000t0001,A¯3=000r22000s22000t2200012,A¯4=000r36000s36000t3600016,B¯0=1h4000B1160000B2160000B3160000B4160,B¯1=1h4B111B121B131B141B211B221B231B241B311B321B331B341B411B421B431B441,D¯0=1h5000D1160000D2160000D3160000D4160,D¯1=1h5D111D121D131D141D211D221D231D241D311D321D331D341D411D421D431D441,while the vectors Y¯m+1, Y¯m1, Y¯m2, Y¯m3, F¯m0, F¯m1, G¯m0, and G¯m1 are defined as follows: (22)Y¯m+1=yn+ryn+syn+tyn+1,Y¯m1=yn-tyn-syn-ryn,Y¯m2=yn-t′yn-s′yn-r′yn′,Y¯m3=yn-t′′yn-s′′yn-r′′yn′′,Y¯m4=yn-t′′′yn-s′′′yn-r′′′yn′′′,F¯m0=fn-3fn-2fn-1fn,F¯m1=fn+rfn+sfn+tfn+1,G¯m0=gn-3gn-2gn-1gn,G¯m1=gn+rgn+sgn+tgn+1, substituting q=λh in (20) yields(23)Y¯m+1=MqY¯m1,where the amplification matrix Mq is given by(24)Mq=A¯0-q4B¯1-q5D¯1-1A¯1+qA¯2+q2A¯3+q3A¯4+q4B¯0+q5D¯0.The analysis shows that the matrix Mq has eigenvalues η1,η2,η3,η4=0,0,0,η4, where the dominant eigenvalue η4 is a function of q subject to the values r,s,t=1/4,1/2,3/4 given by(25)η4=∑i=020ciqiK∑j=021cjqj, where K=10376293541461622784 and the values ci and cj are listed in Table 1.Table 1
c-value
q
i coefficients
q
j coefficients
c
0
−19991259433786613699768237287038031044280320
0
c
1
−19991259433786613699768237287038031044280320
192662816967402505371648000
c
2
−99956297168933068498841186435190155221401600
0
c
3
−33318765722977689499613728811730051740467200
0
c
4
−78083522390790487043756727396404132472422400
0
c
5
−12270936643443632883203625005043197725900800
502432963738164264960000
c
6
−994813550529796720162894578620563036569600
−79502926099016318976000
c
7
59773830578760679822424316692208078028800
0
c
8
31534047193242772464338416945297159618560
0
c
9
4794949827573668397142419401567455150080
133171866570211983360
c
10
359049264395827172053492731789485015040
−14458011284982988800
c
11
−3356754007384253296917217569574748160
268452685173657600
c
12
−4527291309580589178991255608586403840
0
c
13
−632549152682975031987854772642562048
5775969118252032
c
14
−50983394698341254755895386446446592
−659482112163840
c
15
−2347009723724983385107341899462400
26223671235440
c
16
3581359970953878842628978864128
−529590387200
c
17
8897056249188627790299370721984
67083520896
c
18
741826115832694210392858425664
−8638140672
c
19
30585824514148316274163930311
506275974
c
20
582915589719914249658541800
−15909210
c
21
0
189567The region of absolute stability is depicted in the dark area in Figure1.Figure 1
Region of absolute stability.


## 3.6. Pseudocode of the Developed Method

The following pseudocode illustrates how the developed method is implemented for solving fourth-order initial value problems.Step 1.
The valuesyx0,y′x0,y′′x0,y′′′x0 on the interval a,b from the problem are given.Step 2.
Set the values of the hybrid pointsr,s,t and step size h.Step 3.
Forj=a:h:b, set xn+j=xn+jh, j=r,s,t,1.Step 4.
Evaluate the approximate valuesyxn+j, y′xn+j, y′′xn+j, y′′′xn+j using the direct hybrid block method.Step 5.
Find the solution for the resultant system using the built in function fsolve in Matlab.Step 6.
Calculate the maximum error =Exactsolution-approximatesolution.


## 4. Numerical Examples

The general three-hybrid one-step hybrid block method (9) with order p=10 in this section was employed to solve five problems from the literature. In order to show the performance of the proposed method, the problems are solved with different values of step size h and hybrid points r,s,t.Problem 1 (linear fourth-order problem).
(26)
y
i
v
=
y
′
′
′
+
y
′
′
+
y
′
+
2
y
,
0
≤
x
≤
2
,
y
0
=
0
,
y
′
0
=
0
,
y
′
′
0
=
0
,
y
′
′
′
0
=
30
.
Exact Solution. y=2e2x-5e-x+3cos⁡x-9sin⁡x.
Source. See [15].Problem 2 (nonlinear fourth-order problem).
(27)
y
i
v
=
y
′
2
-
y
y
′
′
′
-
4
x
2
+
e
x
1
+
x
2
-
4
x
,
0
≤
x
≤
1
,
y
0
=
1
,
y
′
0
=
1
,
y
′
′
0
=
3
,
y
′
′
′
0
=
1
.
Exact Solution. y=x2+ex with h=1/10.
Source. See [15].Problem 3 (application to problem from ship dynamics).
In this example, the derived method is implemented to solve a physical problem from ship dynamics. As stated by [3], when a sinusoidal wave of frequency passes along a ship or off-shore structure, the resultant fluid actions vary with time. In a particular case study by [3], the fourth-order problem is defined as (28)yiv+3y′′+y2+ϵcos⁡ωx=0,x>0,y0=1,y′0=y′′ρ=y′′′0=0.Exact Solution. yx=2cos⁡x+cos⁡2x.Problem 4 (system of linear fourth-order IVPs).
(29)
y
1
i
v
=
y
1
+
y
4
-
e
2
x
,
y
0
=
1
,
y
′
0
=
1
,
y
′
′
0
=
1
,
y
′
′
′
0
=
1
,
0
≤
x
≤
2
,
y
2
i
v
=
y
2
+
y
3
-
x
e
x
,
y
0
=
1
,
y
′
0
=
-
1
,
y
′
′
0
=
1
,
y
′
′
′
0
=
-
1
,
0
≤
x
≤
2
,
y
3
i
v
=
4
y
1
+
y
3
,
y
0
=
0
,
y
′
0
=
1
,
y
′
′
0
=
2
,
y
′
′
′
0
=
3
,
0
≤
x
≤
2
,
y
4
i
v
=
y
1
+
y
2
+
16
y
4
-
e
x
-
e
-
x
,
y
0
=
1
,
y
′
0
=
2
,
y
′
′
0
=
4
,
y
′
′
′
0
=
8
,
0
≤
x
≤
2
,
Exact Solution. y1=ex, y2=e-x, y3=xex, and y4=e2x with h=1/5.Problem 5 (system of nonlinear fourth-order IVPs).
(30)
y
1
i
v
=
y
1
+
y
2
-
cos
⁡
x
,
y
0
=
0
,
y
′
0
=
1
,
y
′
′
0
=
0
,
y
′
′
′
0
=
-
1
,
0
≤
x
≤
π
,
y
2
i
v
=
y
1
+
y
2
-
sin
⁡
x
,
y
0
=
1
,
y
′
0
=
0
,
y
′
′
0
=
-
1
,
y
′
′
′
0
=
0
,
0
≤
x
≤
π
,
y
3
i
v
=
32
y
1
y
2
,
y
0
=
0
,
y
′
0
=
2
,
y
′
′
0
=
0
,
y
′
′
′
0
=
-
8
,
0
≤
x
≤
π
,
y
4
i
v
=
16
y
2
2
-
y
1
2
,
y
0
=
1
,
y
′
0
=
0
,
y
′
′
0
=
-
2
,
y
′
′
′
0
=
0
,
0
≤
x
≤
π
,
Exact Solution. y1=sin⁡x, y2=cos⁡x, y3=sin⁡2x, and y4=x2+cos⁡2x, with h=π/5.


## 5. Conclusion

A hybrid one-step block method of orderp=10 with generalized three-off-step points has been derived successfully. The developed method was employed to solve general fourth-order IVPs of ODEs directly. The numerical analysis performed shows that the developed method is consistent and zero stable which leads to the convergence of the method. The numerical results were then compared with the results obtained by the existing methods in terms of error. The new method is found to have better performance over the other methods (refer to Tables 2–12).Table 2
Comparison of the derived method with [15, 16] and Adams Bashforth-Adams Moulton method.
Step size
Method employed
Absolute error atx=2
0.1
Derived method
1.34
E
(
-
11
)
BHCM4
1.74
E
(
-
8
)
Jator
1.26
E
(
-
4
)
Adams
2.11
E
(
-
3
)
0.05
Derived method
1.95
E
(
-
11
)
BHCM4
8.45
E
(
-
11
)
Jator
1.91
E
(
-
6
)
Adams
5.37
E
(
-
4
)Table 3
Comparison of the derived method with [15, 17] and Adams Bashforth-Adams Moulton method.
Step size
Method employed
Absolute error atx=1
0.2
Derived method
2.22
E
(
-
15
)
BHCM4
2.38
E
(
-
12
)
Adams
5.01
E
(
-
7
)
Awoyemi
5.84
E
(
-
4
)
0.1
Derived method
1.77
E
(
-
15
)
BHCM4
1.95
E
(
-
14
)
Adams
2.44
E
(
-
6
)
Awoyemi
9.26
E
(
-
5
)Table 4
Comparison of the proposed method with [15] for the case ϵ=0.
Step size
Method employed
Absolute error atx=15
0.25
Derived method
1.1
E
(
-
16
)
BHCM4
5.2
E
(
-
7
)
Twizell
1.9
E
(
-
4
)
Adams
4.9
E
(
-
3
)
0.1
Derived method
1.1
E
(
-
14
)
BHCM4
2.8
E
(
-
10
)
Cortell
3.7
E
(
-
5
)
Adams
8.4
E
(
-
5
)Table 5
Exact and approximate solutions for solvingy1 using the developed method.
X-value
Exact solution
Approximate solution
Error in solvingy1
0.20
1.221402758160169900
1.221402758160169900
0.000000
E
(
+
00
)
0.40
1.491824697641270300
1.491824697641270300
0.000000
E
(
+
00
)
0.60
1.822118800390509100
1.822118800390508900
2.220446
E
(
-
16
)
0.80
2.225540928492467900
2.225540928492467400
4.440892
E
(
-
16
)
1.00
2.718281828459045500
2.718281828459045100
4.440892
E
(
-
16
)
1.20
3.320116922736547200
3.320116922736547700
4.440892
E
(
-
16
)
1.40
4.055199966844674500
4.055199966844674500
0.000000
E
(
+
00
)
1.60
4.953032424395114000
4.953032424395114900
8.881784
E
(
-
16
)
1.80
6.049647464412944800
6.049647464412945700
8.881784
E
(
-
16
)
2.00
7.389056098930648600
7.389056098930650400
1.776357
E
(
-
15
)Table 6
Exact and approximate solutions for solvingy2 using the developed method.
X-value
Exact solution
Approximate solution
Error in solvingy2
0.20
0.818730753077981820
0.818730753077981820
0.000000
E
(
+
00
)
0.40
0.670320046035639330
0.670320046035639330
0.000000
E
(
+
00
)
0.60
0.548811636094026390
0.548811636094026390
0.000000
E
(
+
00
)
0.80
0.449328964117221560
0.449328964117221620
5.551115
E
(
-
17
)
1.00
0.367879441171442330
0.367879441171442330
0.000000
E
(
+
00
)
1.20
0.301194211912202080
0.301194211912202080
0.000000
E
(
+
00
)
1.40
0.246596963941606490
0.246596963941606460
2.775558
E
(
-
17
)
1.60
0.201896517994655440
0.201896517994655380
5.551115
E
(
-
17
)
1.80
0.165298888221586560
0.165298888221586500
5.551115
E
(
-
17
)
2.00
0.135335283236612730
0.135335283236612650
8.326673
E
(
-
17
)Table 7
Exact and approximate solutions for solvingy3 using the developed method.
X-value
Exact solution
Approximate solution
Error in solvingy3
0.20
0.244280551632033990
0.244280551632033970
2.775558
E
(
-
17
)
0.40
0.596729879056508210
0.596729879056508090
1.110223
E
(
-
16
)
0.60
1.093271280234305600
1.093271280234305400
2.220446
E
(
-
16
)
0.80
1.780432742793974300
1.780432742793974100
2.220446
E
(
-
16
)
1.00
2.718281828459045500
2.718281828459045100
4.440892
E
(
-
16
)
1.20
3.984140307283856600
3.984140307283857000
4.440892
E
(
-
16
)
1.40
5.677279953582544000
5.677279953582544000
0.000000
E
(
+
00
)
1.60
7.924851879032181600
7.924851879032183400
1.776357
E
(
-
15
)
1.80
10.889365435943299000
10.889365435943303000
3.552714
E
(
-
15
)
2.00
14.778112197861295000
14.778112197861301000
5.329071
E
(
-
15
)Table 8
Exact and approximate solutions for solvingy4 using the developed method.
X-value
Exact solution
Approximate solution
Error in solvingy4
0.20
1.491824697641270300
1.491824697641270300
0.000000
E
(
+
00
)
0.40
2.225540928492467900
2.225540928492467400
4.440892
E
(
-
16
)
0.60
3.320116922736548100
3.320116922736547700
4.440892
E
(
-
16
)
0.80
4.953032424395114900
4.953032424395114900
0.000000
E
(
+
00
)
1.00
7.389056098930650400
7.389056098930650400
0.000000
E
(
+
00
)
1.20
11.023176380641601000
11.023176380641601000
0.000000
E
(
+
00
)
1.40
16.444646771097048000
16.444646771097048000
0.000000
E
(
+
00
)
1.60
24.532530197109342000
24.532530197109349000
7.105427
E
(
-
15
)
1.80
36.598234443677974000
36.598234443677988000
1.421085
E
(
-
14
)
2.00
54.598150033144215000
54.598150033144236000
2.131628
E
(
-
14
)Table 9
Exact and approximate solutions for solvingy1 using the developed method.
X-value
Exact solution
Approximate solution
Error in solvingy1
0.6283185
0.587785252292473140
0.587785252292475800
2.664535
E
(
-
15
)
1.2566371
0.951056516295153530
0.951056516295303300
1.497691
E
(
-
13
)
1.8849556
0.951056516295153640
0.951056516296290070
1.136424
E
(
-
12
)
2.5132741
0.587785252292473250
0.587785252297075340
4.602096
E
(
-
12
)
3.1415927
0.000000000000000122
0.000000000013537650
1.353753
E
(
-
11
)Table 10
Exact and approximate solutions for solvingy2 using the developed method.
X-value
Exact solution
Approximate solution
Error in solvingy2
0.6283185
0.809016994374947450
0.809016994374955000
7.549517
E
(
-
15
)
1.2566371
0.309016994374947450
0.309016994375326810
3.793632
E
(
-
13
)
1.8849556
−0.309016994374947340
−0.309016994372823590
2.123746
E
(
-
12
)
2.5132741
−0.809016994374947340
−0.809016994368452200
6.495138
E
(
-
12
)
3.1415927
−1.000000000000000000
−0.999999999984998780
1.500122
E
(
-
11
)Table 11
Exact and approximate solutions for solvingy3 using the developed method.
X-value
Exact solution
Approximate solution
Error in solvingy3
0.6283185
0.951056516295153530
0.951056516302223540
7.070011
E
(
-
12
)
1.2566371
0.587785252292473250
0.587785252378536070
8.606282
E
(
-
11
)
1.8849556
−0.587785252292473030
−0.587785251895384220
3.970888
E
(
-
10
)
2.5132741
−0.951056516295153640
−0.951056515173170910
1.121983
E
(
-
09
)
3.1415927
−0.000000000000000245
0.000000002364014064
2.364014
E
(
-
09
)Table 12
Exact and approximate solutions for solvingy4 using the developed method.
X-value
Exact solution
Approximate solution
Error in solvingy4
0.6283185
0.703801170418521750
0.703801170427732160
9.210410
E
(
-
12
)
1.2566371
0.770119709799349850
0.770119709892383430
9.303358
E
(
-
11
)
1.8849556
2.744040590017221300
2.744040590323378700
3.061573
E
(
-
10
)
2.5132741
6.625563811072136200
6.625563811678808300
6.066720
E
(
-
10
)
3.1415927
10.869604401089358000
10.869604401982164000
8.928058
E
(
-
10
)



---

*Source: 7637651-2017-11-12.xml*