# Implicit Two-Point Block Method for Solving Fourth-Order Initial Value Problem Directly with Application

Authors: Reem Allogmany; Fudziah Ismail; Zanariah Abdul Majid; Zarina Bibi Ibrahim
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper proposes an implicit block method with two-point to directly solve the fourth-order Initial Value Problems (IVPs). The implicit block method is derived by adopting Hermite interpolating polynomial as the basis function, incorporating the first derivative offt,y,y′,y′′,y′′′ to enhance the solution’s accuracy. A block formulation is presented to acquire the numerical approximation at two points simultaneously. The introduced method’s basic properties, including order, zero stability, and convergence, are presented. Numerical experiments are carried out to verify the accuracy and efficiency of the proposed method compared with those of the several existing methods. Application in ship dynamics is also presented which yield impressive results for the proposed two-point block method.


---

## Body

## 1. Introduction

Higher order differential equations (ODEs) can be used to model problems arising from the field of applied sciences and engineering [1]. The problem of static deflection of a uniform beam which can be modelled as a fourth-order initial value problem (IVP) is a good example of a real problem in engineering [2, 3]. The classical way of solving them is by reducing the equation into the system of first-order ODEs, but this process is too rigorous compared with the direct methods [4, 5]. In addition, many researchers have presented direct methods to avoid the reduction effort [6–15]. To enhance the efficacy of numerical methods, the block method is introduced with the idea of producing simultaneously r-point of the approximate solutions at one time step. Kuboye and Omar [7] and Jacob [8] introduced efficient zero-stable numerical block methodologies for fourth-order initial value problems (IVPs). For block multistep method, approximations are used to generate the approximate solutions at the new block. Omar and Kuboye [9] proposed a block method for directly solving general fourth-order IVPs by increasing the step number; Adeyeye and Omar [14] also presented a direct six-step block method. Recently, researchers have been exploring numerical methods with more functional evaluations with the aim of obtaining numerical solutions with very high precision. This is how hybrid methods in nature arises. Jator [16] and Yap and Ismail [10] developed four-step hybrid block methods, and Abdelrahim and Omar [11] generalized three off-step points. Kayode et al. [12] modified the implicit hybrid block method to directly solve the IVPs associated with fourth-order ODEs. Generally, collocation and interpolation techniques are utilized as direct methods. The points need to be collocated and interpolated after which a system of linear equations must be resolved in order to obtain the method’s coefficients. Therefore, we developed a two-point block the implicit method by using a strategy in our proposed method which can be implemented in a straightforward manner. The method replaces the function by interpolating and integrating the polynomial. Our derivation only requires us to do interpolation and integration to obtain the coefficients. The method also involves the fifth derivative of the solution to equation (1), which aims at acquiring better accuracy. In general, the accuracy of a method increases with increase in the order of the method. But, with the idea of incorporation of higher derivative of the solution in the process, higher and better accuracy is achieved without corresponding increase in the order of the method. We are concerned with the development of the numerical solution of IVPs for fourth-order ODEs of the form:(1)y4=ft,y,y′,y″,y″′ya=y0,y′a=y0′,y″a=y0″,y″′a=y0″,a≤t≤b,where the fifth derivative of the solution to equation (1) is(2)y5=f′t,y,y′,y″,y″′=ft+fyy′+fy′y″+fy″y″′+ffy″′=gt,y,y′,y″,y″′.We assume thatf in equation (1) is differentiable to a desired order in region R and ft,y,y′,y′′,y′′′ satisfies Lipchitz condition in its second, third, fourth, and fifth terms as follows:(3)ft,y1,y′,y″,y‴−ft,y2,y′,y″,y‴≤Ly1−y2,ft,y,y1′,y″,y‴−ft,y,y2′,y″,y‴≤Ly1′−y2′,ft,y,y′,y1″,y‴−ft,y,y′,y2″,y‴≤Ly1″−y2″,ft,y,y′,y″,y1‴−ft,y,y′,y″,y2‴≤Ly1‴−y2‴,for all points t,yi,y′,y″,y‴, t,y,yi′,y″,y‴, t,y,y′,yi″,y‴, and t,y,y′,y″,yi‴, i=1,2 in the region R. Then, the IVPs in equation (1) have a unique solution in R ([17, 18]).The paper is structured in the following manner: in Section2, a discussion is undertaken about the derivation of the block method of order six (I2PBDO6) along with the basic idea on how the block works. The implementation of this method to solve general fourth-order ODEs is proposed in Section 3. Numerical experiments and the application of the ship dynamics problem are given in Section 4 to show the accuracy and efficiency of the proposed method. The study’s conclusion is finally provided in Section 5.


## 2. Derivation of Method

This section depicts the proposed method’s derivation based on Hermite Interpolating PolynomialP2t, which interpolates ft,y,y′,y″,y″′ at two points. The polynomial has the form(4)Pt=∑i=0n∑k=0mi−1fikLi,kt,where(5)fi=fti,yi,yi′,,yi″,yi″′,ti=a+ih,i=0,1,…,n,and h=b−a/n; Li,kt is the generalized Lagrange polynomial, i=0,1,…,n,k=0,1,…,mi, and n is a positive integer.The method computes the approximate solutionsyn+1 and yn+2 concurrently at tn+1 and tn+2, where tn becomes the starting point and tn+2 is the last point in the block with step size 2h. The evaluation solution of yn+2 at the point tn+2 will be used as the next iteration’s initial value:(6)P2t=L0,0tf0+L1,0tf1+L2,0tf2+L0,1tg0+L1,1tg1+L2,1tg2,where g is the fifth derivative of the solution and(7)L0,0t=t−tn+1tn−tn+12t−tn+2tn−tn+22−2tn−tn+2+2tn−tn+1t−tnt−tn+1tn−tn+12t−tn+2tn−tn+22,L1,0t=t−tntn+1−tn2t−tn+2tn+1−tn+22−2tn+1−tn+2+2tn+1−tnt−tn+1t−tntn+1−tn2t−tn+2tn+1−tn+22,L2,0t=t−tntn+2−tn2t−tn+1tn+2−tn+12−2tn+2−tn+1+2tn+2−tnt−tn+2t−tntn+2−tn2t−tn+1tn+2−tn+12,L0,1t=t−tnt−tn+1tn−tn+12t−tn+2tn−tn+22,L1,1t=t−tn+1t−tntn+1−tn2t−tn+2tn+1−tn+22,L2,1t=t−tn+2t−tntn+2−tn2t−tn+1tn+2−tn+12.The approximate solutionyn+1 at the point tn+1 can be retrieved by integrating equation (1) once, twice, thrice, and four times, in relation to t over the interval tn,tn+1. The following formula can be reached:(8)∫tntn+1y4dt=∫tntn+1ft,y,y′,y″,y″′dt,∫tntn+1∫tnty4dtdt=∫tntn+1∫tntft,y,y′,y″,y″′dtdt,∫tntn+1∫tnt∫tnty4dtdtdt=∫tntn+1∫tnt∫tntft,y,y′,y″,y″′dtdtdt,∫tntn+1∫tnt∫tnt∫tnty4dtdtdtdt=∫tntn+1∫tnt∫tnt∫tntft,y,y′,y″,y″′dtdtdtdt.Lettn+1=tn+h which gives(9)y″′tn+1=y″′tn+∫tntn+1ft,y,y′,y″,y″′dt,(10)y″tn+1=y″tn+hy″′tn+∫tntn+1tn+1−tft,y,y′,y″,y″′dt,(11)y′tn+1=y′tn+hy″tn+h22y″′tn+∫tntn+1tn+1−t22ft,y,y′,y″,y″′dt,(12)ytn+1=ytn+hy′tn+h22y″tn+h36y″′tn+∫tntn+1tn+1−t36ft,y,y′,y″,y″′dt.Then,ft,y,y′,y″,y″′ in equations (9)–(12) will be replaced by Hermite interpolating polynomial P2t.Takingt=tn+2+sh leads to s=t−tn+2/h and dt=hds and then making changes in the integration limit from −2 to −1, we have(13)y″′tn+1=y″′tn+∫−2−1f0L0,0s+f1L1,0s+f2L2,0s+g0L0,1s+g1L1,1s+g2L2,1shds,(14)y″tn+1=y″tn+hy″′tn+∫−2−1−hs+1f0L0,0s+f1L1,0s+f2L2,0s+g0L0,1s+g1L1,1s+g2L2,1shds,(15)y′tn+1=y′tn+hy″tn+h22y″′tn+∫−2−1hs+122f0L0,0s+f1L1,0s+f2L2,0s+g0L0,1s+s+g2L2,1shds,(16)ytn+1=ytn+hy′tn+h22y″tn+h36y″′tn+∫−2−1−hs+136f0L0,0s+f1L1,0s+f2L2,0s+g0L0,1s+s+g2L2,1shds,where(17)L0,0=14s+12s21+32+s,L1,0=2+s2s2,L2,0=142+s2s+121−3s,L0,1=h42+ss+12s2,L1,1=hs+12+s2s2,L2,1=h4ss+122+s2.Evaluating the integrals in equations (13)–(16) gives(18)yn+1‴=yn‴+h101240fn+128240fn+1+11240fn+2+h213240gn−40240gn+1−3240gn+2,(19)yn+1″=yn″+hyn″′+h21342fn+742fn+1+142fn+2+h3591680gn−1281680gn+1−111680gn+2,(20)yn+1′=yn′+hyn″+h22yn‴+h38176720fn+2566720fn+1+476720fn+2+h4836720gn−1406720gn+1−136720gn+2,(21)yn+1=yn+hyn′+h22yn″+h33!yn‴+h4672016fn+142016fn+1+32016fn+2+h518560480gn−25660480gn+1−2560480gn+2.Now, we have four formulas, one for the approximate solution and the others for the approximation of the first, second, and third derivatives of the proposed solution at the block’s end pointtn,tn+1. Therefore, we need to consider the evaluation of p2t at the point yn+2 over tn+1,tn+2 to have a two-point implicit block method. By applying the same technique as for the formula for yn+1, we have the following formula at tn+2:(22)yn+2‴=yn+1‴+h11240fn+128240fn+1+101240fn+2+h23240gn+40240gn+1−13240gn+2,(23)yn+2″=yn+1″+hyn+1‴+h2371680fn+6161680fn+1+1871680fn+2+h35840gn+76840gn+1+16840gn+2,(24)yn+2′=yn+1′+hyn+1″+h22yn+1‴+h3416720fn+9286720fn+1+1516720fn+2+h4116720gn+1886720gn+1−296720gn+2,(25)yn+2=yn+1+hyn+1′+h22yn+1″+h33!yn+1‴+h454032fn+1484032fn+1+154032fn+2+h51030240gn+19030240gn+1−2330240gn+2.


### 2.1. Order of the Method

The equations (15)–(22) can be written in a matrix difference equation to obtain the order of the new method as follows:(26)αYm=hβYm′+h2γYm″+h3δYm″′+h4ζFm+h5λGm,where β,γ, and δ are defined as(27)α=00000000000000000000000000000−110000000000000000000000000000000−11,β=0000000000000000000001−100000010000000000000000000000001−100000010,γ=00000000000001−1000000100000001200000000000000001−100000010000000120,δ=000001−10000001000000012000000016000000001−100000010000000120000000160,ζ=00000101240128240112400000013427421420000081767202566720476720000006720161420163201600000112401282401012400000037168011301871680000004167209286720151672000000540321484032154032,λ=0000013240−40240−324000000591680−1281680−11168000000836720−1406720−1367200000018560480−25660480−256048000000324040240−1324000000584076840−16840000001167201886720−29672000000103024019030240−2330240,Ym=yn−5yn−4yn−3yn−2yn−1ynyn+1yn+2,Ym′=yn−5′yn−4′yn−3′yn−2′yn−1′yn′yn+1′yn+2′,Ym″=yn−5″yn−4″yn−3″yn−2″yn−1″yn″yn+1″yn+2″,Ym″′=yn−5″′yn−4″′yn−3″′yn−2″′yn−1″′yn″′yn+1″′yn+2″′,Fm=fn−5fn−4fn−3fn−2fn−1fnfn+1fn+2,Gm=gn−5gn−4gn−3gn−2gn−1gngn+1gn+2.It is possible to define the linear operator that is linked to equation (26) as(28)ℓyt:h=αYm−hβYm′−h2γYm″−h3δYm″′−h4ζFm−h5λGm.Expanding equation (28) in the Taylor series yields(29)ℓyt:h=C0yt+C1hy′t+C2h2y″t+C3h3y″′t+⋯+Cphpypt+Cp+1hp+1yp+1t+⋯.The linear operator and the new method have orderp if C0=C1=⋯=Cp=⋯=Cp+3=0, Cp+4≠0, which is the error constant. Hence, in our method, C0=C1=⋯=C9=0 and C10=1.058199×10−4,5.78707×10−5,1.763674×10−5,3.858×10−6,1.0581959×10−4,4.795011×10−5,1.267633×10−5,2.479443×10−6T. Thus, it can be concluded that the order of the proposed method is 6.


### 2.2. Zero Stability

The method (18)–(25) is deemed zero-stable in case the roots zi=1,2,…,N of the first characteristic polynomial ρz=detzA0−A1 is found to satisfy zi|≤1.Now, in our proposed method, we will use the following technique to find the matrix of the first characteristic polynomial.By replacing equation (18) into equation (22), we have(30)yn+2‴=yn‴+h112240fn+256240fn+1+112240fn+2+h216240gn−16240gn+2.Also, substituting equations (18) and (19) into equation (23), we get(31)yn+2″=yn″+2hyn‴+h212641680fn+17921680fn+1+3041680fn+2+h310105gn−16105gn+1−4105gn+2.Substituting equations (18)–(21) into equation (24), we have(32)yn+2′=yn′+2hyn″+2h2yn‴+h368105fn+64105fn+1+8105fn+2+h45126720gn−11656720gn+1−1286720gn+2.By substituting equations (18)–(22) into equation (25), we have(33)yn+2=yn+2hyn′+2h2yn″+4h3yn‴+h42463fn+1663fn+1+263fn+2+h525606720gn−51206720gn+1−5126720gn+2.By using the equations (18)–(22) and (30)–(33), the general form of the matrix A0andA1 can be written as A0=8×8, identity matrix,(34)A1=0000100000000100000000100000000100001000000001000000001000000001.Then,ρz=detzA0−A1=0 implies that ρz=detz4z−14. Hence, z=0,0,0,0,1,1,1,1Therefore, it can be inferred that the proposed method has zero stability.


## 2.1. Order of the Method

The equations (15)–(22) can be written in a matrix difference equation to obtain the order of the new method as follows:(26)αYm=hβYm′+h2γYm″+h3δYm″′+h4ζFm+h5λGm,where β,γ, and δ are defined as(27)α=00000000000000000000000000000−110000000000000000000000000000000−11,β=0000000000000000000001−100000010000000000000000000000001−100000010,γ=00000000000001−1000000100000001200000000000000001−100000010000000120,δ=000001−10000001000000012000000016000000001−100000010000000120000000160,ζ=00000101240128240112400000013427421420000081767202566720476720000006720161420163201600000112401282401012400000037168011301871680000004167209286720151672000000540321484032154032,λ=0000013240−40240−324000000591680−1281680−11168000000836720−1406720−1367200000018560480−25660480−256048000000324040240−1324000000584076840−16840000001167201886720−29672000000103024019030240−2330240,Ym=yn−5yn−4yn−3yn−2yn−1ynyn+1yn+2,Ym′=yn−5′yn−4′yn−3′yn−2′yn−1′yn′yn+1′yn+2′,Ym″=yn−5″yn−4″yn−3″yn−2″yn−1″yn″yn+1″yn+2″,Ym″′=yn−5″′yn−4″′yn−3″′yn−2″′yn−1″′yn″′yn+1″′yn+2″′,Fm=fn−5fn−4fn−3fn−2fn−1fnfn+1fn+2,Gm=gn−5gn−4gn−3gn−2gn−1gngn+1gn+2.It is possible to define the linear operator that is linked to equation (26) as(28)ℓyt:h=αYm−hβYm′−h2γYm″−h3δYm″′−h4ζFm−h5λGm.Expanding equation (28) in the Taylor series yields(29)ℓyt:h=C0yt+C1hy′t+C2h2y″t+C3h3y″′t+⋯+Cphpypt+Cp+1hp+1yp+1t+⋯.The linear operator and the new method have orderp if C0=C1=⋯=Cp=⋯=Cp+3=0, Cp+4≠0, which is the error constant. Hence, in our method, C0=C1=⋯=C9=0 and C10=1.058199×10−4,5.78707×10−5,1.763674×10−5,3.858×10−6,1.0581959×10−4,4.795011×10−5,1.267633×10−5,2.479443×10−6T. Thus, it can be concluded that the order of the proposed method is 6.


## 2.2. Zero Stability

The method (18)–(25) is deemed zero-stable in case the roots zi=1,2,…,N of the first characteristic polynomial ρz=detzA0−A1 is found to satisfy zi|≤1.Now, in our proposed method, we will use the following technique to find the matrix of the first characteristic polynomial.By replacing equation (18) into equation (22), we have(30)yn+2‴=yn‴+h112240fn+256240fn+1+112240fn+2+h216240gn−16240gn+2.Also, substituting equations (18) and (19) into equation (23), we get(31)yn+2″=yn″+2hyn‴+h212641680fn+17921680fn+1+3041680fn+2+h310105gn−16105gn+1−4105gn+2.Substituting equations (18)–(21) into equation (24), we have(32)yn+2′=yn′+2hyn″+2h2yn‴+h368105fn+64105fn+1+8105fn+2+h45126720gn−11656720gn+1−1286720gn+2.By substituting equations (18)–(22) into equation (25), we have(33)yn+2=yn+2hyn′+2h2yn″+4h3yn‴+h42463fn+1663fn+1+263fn+2+h525606720gn−51206720gn+1−5126720gn+2.By using the equations (18)–(22) and (30)–(33), the general form of the matrix A0andA1 can be written as A0=8×8, identity matrix,(34)A1=0000100000000100000000100000000100001000000001000000001000000001.Then,ρz=detzA0−A1=0 implies that ρz=detz4z−14. Hence, z=0,0,0,0,1,1,1,1Therefore, it can be inferred that the proposed method has zero stability.


## 3. Implementation

In this section, we briefly mention the developed method’s implementation. The approximate solutions in equations (15)–(22) will be estimated using the predictor-corrector schemes. The predictor equations using Taylor method are(35)yn+i″′p=yn+i−1″′c+hfn+i−1c,yn+i″p=yn+i−1″c+hyn+i−1″′c+h22!fn+i−1c,yn+i′p=yn+i−1′c+hyn+i−1″c+h22!y″′cn+i−1+h33!fn+i−1c,yn+ip=yn+i−1c+hyn+i−1′c+h22!yn+i−1″c+h33!yn+i−1″′c+h44!fn+i−1c,fn+ipt,y,y′,y″,y″′=ftn+i,yn+ip,yn+i′p,yn+i″p,yn+i″′p,gn+ipt,y,y′,y″,y″′=f′tn+i,yn+ip,yn+i′p,yn+i″p,yn+i″′p,i=1,2.Define equation (35) as the initial approximation and apply the derived method in equations (18)–(25) as the corrector. In order to obtain the corrector iteration, the following equations will be used(36)yn+1″′c=yn″′c+h240101fnc+128fn+1p+11fn+2p+h224013gnc−40gn+1p−3gn+2p,yn+1″c=yn″c+hyn″′c+h4213fnc+7fn+1p+1fn+2p+h3168059gnc−128gn+1p−11gn+2p,yn+1′c=yn′c+hyn″c+h22yn″c+h36720817fnc+256fn+1p+47fn+2p+h4672083gnc−140gn+1p−13gn+2p,yn+1c=ync+hyn′c+h22y′′nc+h33!yn″′c+h4201667fnc+14fn+1p+3fn+2p+h560480185gnc−256gn+1p−25gn+2p.Then, we evaluate the functionsfn+1c and gn+1c which we use to compute the following approximate solutions at the second point:(37)yn+2″′c=yn+1″′c+h24011fnc+128fn+1c+101fn+2p+h22403gnc+40gn+1c−13gn+2p,yn+2″c=yn+1″c+hyn+1″′c+h2168037fnc+616fn+1c+187fn+2p+h38405gnc+76gn+1c16gn+2p,yn+2′c=yn+1′c+hyn+1″c+h22yn+1″′c+h3672041fnc+928fn+1c+151fn+2p+h4672011gnc+188gn+1c−29gn+2p,yn+2c=yn+1c+hyn+1′c+h22yn+1″c+h33!yn+1″′c+h440325fnc+148fn+1c+15fn+2p+h53024010gnc+190gn+1c−23gn+2p.Next, we evaluate the functionsfn+2c and gn+2c that we will use in the next corrector iteration. Hence, we repeat the procedure given in equations (36) and (37) to perform the next corrector iterations until the end of the interval.


## 4. Results and Discussion

Within this section, we tested the I2PBDO6 method to solve single as well as a system of initial value challenges in the form of equation (1). We will see how the method of two-point block (order 6) produces comparable results to the existing direct methods of similar characteristics. Most of the existing methods are of higher order in comparison to the proposed method. The following notations will be used in the tables:(i)
I2PBDO6: the direct implicit two-point block method proposed in this paper, order 6(ii)
MCM: direct maximal order multiderivative collocation method [6](iii)
ZSB7: the direct zero-stable block method of order seven [7](iv)
ZSDM6: the direct zero-stable method of order six [8](v)
IHB8: the direct implicit hybrid block method of order eight [12](vi)
FDM6: the direct linear multistep method of order six [16](vii)
BHC8: the direct block hybrid collocation method of order eight [10](viii)
FSHB8: the direct four-step hybrid block method of order eight [11](ix)
SSB8: the direct seven-step block method of order eight [9](x)
LIBO7: the direct implicit six-step block method of order seven [14](xi)
RKTF5: the direct RKTF method of order five with four-stage [13](xii)
RKF5: the RK method with six-stage given in Lambert [19](xiii)
h: step size(xiv)
FCN: number of function evaluations(xv)
NS: number of steps(xvi)
MAXE: maximum absolute error


### 4.1. Problems

Problem 1.
(38)y4+y′′=0,0≤t≤π2,y0=0,y′0=−1.172−50π,y′′0=1144−100π,y′′′0=1.2144−100π.
Exact solution:(39)yt=1−t−cost−1.2sint144−100π.Problem 2.
(40)y4=y′2−yy′′′−4t2+et1−4t+t2,0≤t≤1,y0=1,y′0=1,y″0=3,y″′0=1.
Exact solution:(41)yt=t2+et.Problem 3.
(42)y4=y″′+y″+y′+2y,t∈0,2,y0=0,y′0=0,y″0=0,y″′0=30.
Exact solution:(43)yt=2e2t−5e−t+3cost−9sint.Problem 4.
(44)y4−t4+14t3+49t2+32t−12et=0,y0=0,y′0=0,y″0=2,y″′0=−6.
Exact solution:(45)yt=t21−t2et.Problem 5.
(46)y4=−ty−8+7t+t3et,t∈0,2,y0=0,y′0=1,y″0=0,y″′0=−3.
Exact solution:(47)yt=t1−tet.Problem 6.
(48)y14=−y2″,y10=1,y1′0=1,y1′′0=1,y1″′0=1,y24=−y1″,y20=−1,y2′0=−1,y2″0=−1,y2″′0=−1,y34=−y3″−y3−cost,y30=−1,y3′0=0,y3″0=1,y3″′0=0,y44=−y4″−y4−2cost,y40=−2,y4′0=0,y4″0=2,y4‴0=0.
Exact solution:(49)y1t=et,y2t=−et,y3t=−cost,y4t=−2cost.Problem 7.
(50)y14=e3ty4,y10=1,y1′0=−1,y1″0=1,y1″′0=−1,y24=16e−ty1,y20=1,y2′0=−2,y2″0=4,y2″′0=−8,y34=81e−ty2,y30=1,y3′0=−3,y3″0=9,y3″′0=−27,y44=256e−ty3,y40=1,y4′0=−4,y4″0=16,y4″′0=−64.
Exact solution:(51)y1t=e−t,y2t=e−2t,y3t=e−3t,y4t=e−4t,0≤t≤3.
The aim of this study was to demonstrate the accuracy, efficiency, and applicability of the proposed I2PBDO6 method. Tables1–4 show the absolute errors at different points as well as different step size. It is remarkable that the errors are zero at some points of the interval which means that the computed solutions are identical to the exact solutions with 20 digits at least. The two-point block method can be seen to significantly outperform MCM [6], ZSB7 [7], LIBO7 [14], ZSDM6 [8], and IHB8 [12] with regard to accuracy. For Problems 3 to 5, the results of the proposed I2PBDO6 method are tabulated alongside others in Tables 5–7. We have considered the absolute errors along the interval using different number of steps in terms of accuracy as well as steps taken. It can be observed that increasing the number of steps leads to decrease in the error except that in Table 6 at the number of Step 5004 and 50004, since Problem 4 involves exponential function which requires an infinite number of digits for numerical representation. However, these numbers must be represented by a finite number of digits in computer for numerical computation that causes round-off errors. Thus, at a large number of steps such as 5004 and 50004, a sequence of many operations and each subject to a small rounding error is obtained, where the error accumulates and affects the accuracy of the final solutions. The proposed method shows very good performance compared with FDM6 [16], BHC8 [10], FSHB8 [11], SSB8 [9], and LIBO7 [14]. The system in Problem 6 and Problem 7 has integrated in the interval t∈0,10 and t∈0,3, respectively. A comparison of numerical outcomes is made with the solution of the same group of systems when they are solved by using RKTF5 [13] and RKF5 [19] methods. We have solved these systems to consider the maximum errors and the number of function evaluations at different step sizes, as shown in Tables 8 and 9. As we can see, the efficiency of the I2PBDO6 method is improving clearly when smaller step sizes are used. Based on the efficiency curves in Figures 1 and 2 of log (max error) versus log (function evaluations), the proposed I2PBDO6 method achieved smaller error for the same number of function evaluations compared with RKTF5 and RKF5. Hence, I2PBDO6 accumulated more total function calls; nevertheless, Figures 1 and 2 display the best performances and efficiency of the proposed I2PBDO6 method.
Overall, the proposed method I2PBDO6 of order 6 has shown remarkable convergence since the approximate solution is very close to the exact solution. Moreover, the method is more efficient compared with the all present methods whose order is nearly equal to or greater as compared with that of the proposed method. The data given in Tables1–9 and Figures 1 and 2 show the superiority of I2PBDO6 with regard to efficiency, accuracy, and total number of steps taken at different points of t or different step sizes.Table 1
Numerical findings to solve Problem1 with step size h=1/320.
tMCM [6]I2PBDO60.1031250.215807406667 − 122.165587 − 190.206250.235634428636 − 122.162928 − 180.306250.250504044591 − 126.913642 − 180.406250.258852835000 − 121.122300 − 170.506250.257955115601 − 121.466220 − 170.6031250.245233520990 − 124.308768 − 180.7031250.216478744652 − 128.610773 − 180.8031250.168292463298 − 121.893054 − 170.9031250.969693075836 − 132.923863 − 17Table 2
Numerical findings to solve Problem1 with step size h=0.01.
tZSB7 [7]I2PBDO60.015.421011 − 205.420312 − 200.025.421011 − 205.419616 − 200.032.710505 − 192.709462 − 190.041.084202 − 190.000000 + 000.053.252607 − 193.250530 − 190.063.252607 − 193.250120 − 190.081.734723 − 180.000000 + 000.094.336809 − 182.165937 − 190.108.456777 − 180.000000 + 00Table 3
Numerical findings to solve Problem1 with step size h=0.1.
tLIBO7 [14]I2PBDO60.16.51 − 190.000000 + 000.28.67 − 198.652339 − 190.34.34 − 194.321318 − 190.41.73 − 181.590085 − 180.58.67 − 193.450126 − 190.60.00 + 001.723552 − 190.70.00 + 000.000000 + 000.80.00 + 000.000000 + 000.90.00 + 000.000000 + 001.00.00 + 001.719064 − 19Table 4
Numerical findings to solve Problem2 with step size h=1/320.
tZSDM6 [8]IHB8 [12]I2PBDO60.1031250.666133814775 − 144.261834 − 140.0000000 + 000.2062500.974775815621 − 137.314820 − 133.9099252 − 160.3062500.461852778244 − 121.625321 − 127.2441945 − 160.4062500.142241773915 − 116.568972 − 129.9936929 − 160.5062500.342548212018 − 112.849147 − 121.2947886 − 150.6031250.687361279006 − 114.231457 − 120.0000000 + 000.7031250.126365584663 − 101.246852 − 111.5163355 − 150.8031250.213882245248 − 105.146142 − 112.8632323 − 150.9031250.339896999435 − 103.765423 − 114.5622733 − 15Table 5
Numerical findings to solve Problem3.
NSFDM6BHC8FSHB8I2PBDO6201.25612 − 041.74 − 088.07 − 101.851538 − 10401.90752 − 068.45 − 113.22 − 122.328543 − 12802.96411 − 083.69 − 13—2.115566 − 141008.65117 − 097.11 − 14—1.079081 − 14Table 6
Numerical findings to solve Problem4, t∈0,1.
NSSSB8I2PBDO6543.547029 − 111.455049 − 145041.141416 − 106.226973 − 1550041.439730 − 094.668800 − 14500042.466095 − 094.156048 − 13Table 7
Numerical findings to solve Problem5.
NSLIBO7 [14]I2PBDO6832.61 − 132.727680 − 141452.84 − 144.776630 − 142343.02 − 143.818367 − 15Table 8
Numerical findings to solve Problem6.
hFCNMethodMAXE0.1700I2PBDO63.924317 − 10404RKTF55.327998 − 52424RKF52.064967 − 30.0252800I2PBDO61.691679 − 131600RKTF54.254471 − 79600RKF51.917218 − 60.01255600I2PBDO68.014250 − 143204RKTF51.786611 − 819224RKF56.053233 − 80.0062511200I2PBDO83.008073 − 146404RKTF56.002665 − 938424RKF54.878530 − 9Table 9
Numerical findings to solve Problem7.
hFCNMethodMAXE0.1210I2PBDO61.375600 − 3120RKTF51.534759 − 1720RKF56.153184 − 10.025840I2PBDO64.992115 − 7484RKTF53.920592 − 52904RKF56.983629 − 40.01251680I2PBDO68.294011 − 9960RKTF52.560210 − 65760RKF51.960240 − 50.006253360I2PBDO81.330643 − 101920RKTF51.468912 − 711520RKF54.878530 − 9Figure 1
Efficiency curves for I2PBDO6, RKTF5, and RKF5 for solving Problem6.Figure 2
Efficiency curves for I2PBDO6, RKTF5, and RKF5 for solving Problem7.


### 4.2. Application to Ship Dynamics Problem

We take into consideration a physical challenge from the ship dynamics as shown in Wu et al. [20], where a sinusoidal wave passes via a ship or offshore structure, with the accompanying fluid actions varying with time t. In the aforementioned study [20], the differential equation is defined as(52)y4+3y″+y2+εcosΩt=0,t>0,with(53)y0=1,y′0=y″0=y″′=0,where ε=0 with regard to the theoretical solution, yt=2cost−cos2t. On the other hand, the theoretical solution can be seen to be indeterminate when ε≠0 as can be seen in [21]. Numerical investigations to resolve the ordinary differential equations of the fourth-order have been expanded for resolving a problem from ship dynamics. Numerical consideration was offered when ε=0 and ε=1 with Ω=0.252−1 as can be seen in Twizell [21] and Cortell [22]. However, Twizell [21] and Cortell [22] derived the ODE method by lowering the equations to system of first-order (ODEs) rather than directly resolving it. Moreover, Cortell [22] introduced the extension of the classical RK method. Figure 3 depicts the exact and approximate solutions by using I2PBDO6 for the application of ship dynamics problem with ε=0. However, the results of equation (52), ε=1, and Ω=0.252−1 in 0,15 are shown in 4. It is obtained by I2PBDO6 with the solution obtained by Mathematica in-built package NDSolve, where h=0.1. From Figures 3 and 4, it is clear that the solutions that are obtained by I2PBDO6 strongly agree with not only the exact solutions but also with those of Mathematica in-built package NDSolve, respectively.Figure 3
Response curve concerning equation (52) with ε=0 and h=0.1 in 0,15.Figure 4
Response curve concerning equation (52) with ε=1 and Ω=0.252−1 and h=0.1 in 0,15.


## 4.1. Problems

Problem 1.
(38)y4+y′′=0,0≤t≤π2,y0=0,y′0=−1.172−50π,y′′0=1144−100π,y′′′0=1.2144−100π.
Exact solution:(39)yt=1−t−cost−1.2sint144−100π.Problem 2.
(40)y4=y′2−yy′′′−4t2+et1−4t+t2,0≤t≤1,y0=1,y′0=1,y″0=3,y″′0=1.
Exact solution:(41)yt=t2+et.Problem 3.
(42)y4=y″′+y″+y′+2y,t∈0,2,y0=0,y′0=0,y″0=0,y″′0=30.
Exact solution:(43)yt=2e2t−5e−t+3cost−9sint.Problem 4.
(44)y4−t4+14t3+49t2+32t−12et=0,y0=0,y′0=0,y″0=2,y″′0=−6.
Exact solution:(45)yt=t21−t2et.Problem 5.
(46)y4=−ty−8+7t+t3et,t∈0,2,y0=0,y′0=1,y″0=0,y″′0=−3.
Exact solution:(47)yt=t1−tet.Problem 6.
(48)y14=−y2″,y10=1,y1′0=1,y1′′0=1,y1″′0=1,y24=−y1″,y20=−1,y2′0=−1,y2″0=−1,y2″′0=−1,y34=−y3″−y3−cost,y30=−1,y3′0=0,y3″0=1,y3″′0=0,y44=−y4″−y4−2cost,y40=−2,y4′0=0,y4″0=2,y4‴0=0.
Exact solution:(49)y1t=et,y2t=−et,y3t=−cost,y4t=−2cost.Problem 7.
(50)y14=e3ty4,y10=1,y1′0=−1,y1″0=1,y1″′0=−1,y24=16e−ty1,y20=1,y2′0=−2,y2″0=4,y2″′0=−8,y34=81e−ty2,y30=1,y3′0=−3,y3″0=9,y3″′0=−27,y44=256e−ty3,y40=1,y4′0=−4,y4″0=16,y4″′0=−64.
Exact solution:(51)y1t=e−t,y2t=e−2t,y3t=e−3t,y4t=e−4t,0≤t≤3.
The aim of this study was to demonstrate the accuracy, efficiency, and applicability of the proposed I2PBDO6 method. Tables1–4 show the absolute errors at different points as well as different step size. It is remarkable that the errors are zero at some points of the interval which means that the computed solutions are identical to the exact solutions with 20 digits at least. The two-point block method can be seen to significantly outperform MCM [6], ZSB7 [7], LIBO7 [14], ZSDM6 [8], and IHB8 [12] with regard to accuracy. For Problems 3 to 5, the results of the proposed I2PBDO6 method are tabulated alongside others in Tables 5–7. We have considered the absolute errors along the interval using different number of steps in terms of accuracy as well as steps taken. It can be observed that increasing the number of steps leads to decrease in the error except that in Table 6 at the number of Step 5004 and 50004, since Problem 4 involves exponential function which requires an infinite number of digits for numerical representation. However, these numbers must be represented by a finite number of digits in computer for numerical computation that causes round-off errors. Thus, at a large number of steps such as 5004 and 50004, a sequence of many operations and each subject to a small rounding error is obtained, where the error accumulates and affects the accuracy of the final solutions. The proposed method shows very good performance compared with FDM6 [16], BHC8 [10], FSHB8 [11], SSB8 [9], and LIBO7 [14]. The system in Problem 6 and Problem 7 has integrated in the interval t∈0,10 and t∈0,3, respectively. A comparison of numerical outcomes is made with the solution of the same group of systems when they are solved by using RKTF5 [13] and RKF5 [19] methods. We have solved these systems to consider the maximum errors and the number of function evaluations at different step sizes, as shown in Tables 8 and 9. As we can see, the efficiency of the I2PBDO6 method is improving clearly when smaller step sizes are used. Based on the efficiency curves in Figures 1 and 2 of log (max error) versus log (function evaluations), the proposed I2PBDO6 method achieved smaller error for the same number of function evaluations compared with RKTF5 and RKF5. Hence, I2PBDO6 accumulated more total function calls; nevertheless, Figures 1 and 2 display the best performances and efficiency of the proposed I2PBDO6 method.
Overall, the proposed method I2PBDO6 of order 6 has shown remarkable convergence since the approximate solution is very close to the exact solution. Moreover, the method is more efficient compared with the all present methods whose order is nearly equal to or greater as compared with that of the proposed method. The data given in Tables1–9 and Figures 1 and 2 show the superiority of I2PBDO6 with regard to efficiency, accuracy, and total number of steps taken at different points of t or different step sizes.Table 1
Numerical findings to solve Problem1 with step size h=1/320.
tMCM [6]I2PBDO60.1031250.215807406667 − 122.165587 − 190.206250.235634428636 − 122.162928 − 180.306250.250504044591 − 126.913642 − 180.406250.258852835000 − 121.122300 − 170.506250.257955115601 − 121.466220 − 170.6031250.245233520990 − 124.308768 − 180.7031250.216478744652 − 128.610773 − 180.8031250.168292463298 − 121.893054 − 170.9031250.969693075836 − 132.923863 − 17Table 2
Numerical findings to solve Problem1 with step size h=0.01.
tZSB7 [7]I2PBDO60.015.421011 − 205.420312 − 200.025.421011 − 205.419616 − 200.032.710505 − 192.709462 − 190.041.084202 − 190.000000 + 000.053.252607 − 193.250530 − 190.063.252607 − 193.250120 − 190.081.734723 − 180.000000 + 000.094.336809 − 182.165937 − 190.108.456777 − 180.000000 + 00Table 3
Numerical findings to solve Problem1 with step size h=0.1.
tLIBO7 [14]I2PBDO60.16.51 − 190.000000 + 000.28.67 − 198.652339 − 190.34.34 − 194.321318 − 190.41.73 − 181.590085 − 180.58.67 − 193.450126 − 190.60.00 + 001.723552 − 190.70.00 + 000.000000 + 000.80.00 + 000.000000 + 000.90.00 + 000.000000 + 001.00.00 + 001.719064 − 19Table 4
Numerical findings to solve Problem2 with step size h=1/320.
tZSDM6 [8]IHB8 [12]I2PBDO60.1031250.666133814775 − 144.261834 − 140.0000000 + 000.2062500.974775815621 − 137.314820 − 133.9099252 − 160.3062500.461852778244 − 121.625321 − 127.2441945 − 160.4062500.142241773915 − 116.568972 − 129.9936929 − 160.5062500.342548212018 − 112.849147 − 121.2947886 − 150.6031250.687361279006 − 114.231457 − 120.0000000 + 000.7031250.126365584663 − 101.246852 − 111.5163355 − 150.8031250.213882245248 − 105.146142 − 112.8632323 − 150.9031250.339896999435 − 103.765423 − 114.5622733 − 15Table 5
Numerical findings to solve Problem3.
NSFDM6BHC8FSHB8I2PBDO6201.25612 − 041.74 − 088.07 − 101.851538 − 10401.90752 − 068.45 − 113.22 − 122.328543 − 12802.96411 − 083.69 − 13—2.115566 − 141008.65117 − 097.11 − 14—1.079081 − 14Table 6
Numerical findings to solve Problem4, t∈0,1.
NSSSB8I2PBDO6543.547029 − 111.455049 − 145041.141416 − 106.226973 − 1550041.439730 − 094.668800 − 14500042.466095 − 094.156048 − 13Table 7
Numerical findings to solve Problem5.
NSLIBO7 [14]I2PBDO6832.61 − 132.727680 − 141452.84 − 144.776630 − 142343.02 − 143.818367 − 15Table 8
Numerical findings to solve Problem6.
hFCNMethodMAXE0.1700I2PBDO63.924317 − 10404RKTF55.327998 − 52424RKF52.064967 − 30.0252800I2PBDO61.691679 − 131600RKTF54.254471 − 79600RKF51.917218 − 60.01255600I2PBDO68.014250 − 143204RKTF51.786611 − 819224RKF56.053233 − 80.0062511200I2PBDO83.008073 − 146404RKTF56.002665 − 938424RKF54.878530 − 9Table 9
Numerical findings to solve Problem7.
hFCNMethodMAXE0.1210I2PBDO61.375600 − 3120RKTF51.534759 − 1720RKF56.153184 − 10.025840I2PBDO64.992115 − 7484RKTF53.920592 − 52904RKF56.983629 − 40.01251680I2PBDO68.294011 − 9960RKTF52.560210 − 65760RKF51.960240 − 50.006253360I2PBDO81.330643 − 101920RKTF51.468912 − 711520RKF54.878530 − 9Figure 1
Efficiency curves for I2PBDO6, RKTF5, and RKF5 for solving Problem6.Figure 2
Efficiency curves for I2PBDO6, RKTF5, and RKF5 for solving Problem7.


## 4.2. Application to Ship Dynamics Problem

We take into consideration a physical challenge from the ship dynamics as shown in Wu et al. [20], where a sinusoidal wave passes via a ship or offshore structure, with the accompanying fluid actions varying with time t. In the aforementioned study [20], the differential equation is defined as(52)y4+3y″+y2+εcosΩt=0,t>0,with(53)y0=1,y′0=y″0=y″′=0,where ε=0 with regard to the theoretical solution, yt=2cost−cos2t. On the other hand, the theoretical solution can be seen to be indeterminate when ε≠0 as can be seen in [21]. Numerical investigations to resolve the ordinary differential equations of the fourth-order have been expanded for resolving a problem from ship dynamics. Numerical consideration was offered when ε=0 and ε=1 with Ω=0.252−1 as can be seen in Twizell [21] and Cortell [22]. However, Twizell [21] and Cortell [22] derived the ODE method by lowering the equations to system of first-order (ODEs) rather than directly resolving it. Moreover, Cortell [22] introduced the extension of the classical RK method. Figure 3 depicts the exact and approximate solutions by using I2PBDO6 for the application of ship dynamics problem with ε=0. However, the results of equation (52), ε=1, and Ω=0.252−1 in 0,15 are shown in 4. It is obtained by I2PBDO6 with the solution obtained by Mathematica in-built package NDSolve, where h=0.1. From Figures 3 and 4, it is clear that the solutions that are obtained by I2PBDO6 strongly agree with not only the exact solutions but also with those of Mathematica in-built package NDSolve, respectively.Figure 3
Response curve concerning equation (52) with ε=0 and h=0.1 in 0,15.Figure 4
Response curve concerning equation (52) with ε=1 and Ω=0.252−1 and h=0.1 in 0,15.


## 5. Conclusions

In this article, we proposed a two-point block method with extra derivative to directly solve initial value problems of the fourth order. The derivation of the new method can be easily implemented. This implicit block method is found to satisfy the property of convergence as suggested by the significant enhancement with regard to accuracy in the numerical results. Therefore, the proposed method I2PBDO6 is suitable for attaining high precision when solving fourth-order ODEs.



---

*Source: 6351279-2020-09-15.xml*