J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 Journal of the Nigerian Society of Physical Sciences Efficient Hybrid Block Method For The Numerical Solution Of Second-order Partial Differential Problems via the Method of Lines Olumide O. Olaiyaa,, Razaq A. Azeezb, Mark I. Modebeia aDepartment of Mathematics Programme, National Mathematical Centre, Abuja, Nigeria bDepartment of Mathematics, University of Abuja, Abuja, Nigeria Abstract This study is therefore aimed at developing classes of efficient numerical integration schemes, for direct solution of second-order Partial Differ- ential Equations (PDEs) with the aid of the method of lines. The power series polynomials were used as basis functions for trial solutions in the derivation of the proposed schemes via collocation and interpolation techniques at some appropriately chosen grid and off-grid points the derived schemes are consistent, zero-stable and convergent. the proposed methods perform better in terms of accuracy than some existing methods in the literature. DOI:10.46481/jnsps.2021.141 Keywords: Initial Value Problem, Boundary Value Problem; Block method, Linear Multistep Method, Hybrid method, method od lines Article History : Received: 02 October 2020 Received in revised form: 10 December 2020 Accepted for publication: 23 January 2021 Published: 27 February 2021 c©2021 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: T. Latunde 1. Introduction In this work, we consider the second-order PDE of the form A ∂y2 ∂x2 + B ∂y2 ∂t2 +C ∂y ∂t = 0, x∈ [a, b] , t > 0 With any of the following initial-boundary conditions y(x,α1) = ξ1, y(x,α2) = ξ2, y(β1t) = ζ1, y(β2, t) = ζ2, y(x,α1) = ξ1, y(β2, t) = ζ2 Email address: olaiyaolumide.o@gmail.com (Olumide O. Olaiya ) , where αi, βi, ζi, ξi, i = 1, 2, a, b, A, B, C are all real constants and A , 0. y ∈ C2 [a, b] × [c, d] , x ∈ (a, b) , t ∈ (c, d). Second-order PDEs can either be of Laplacian, Poisson forms which could either be heat or waveform of equations. They find their applications in numerous areas of human endeavours, especially in mathematical sciences and engineering. The developed differential equations require solutions, ei- ther in closed form (analytic form) or in numerical forms. In most cases, closed-form of solutions is rare to come by as there exist limited methods for solving such models in the form of differential equations. This brings to light the use of numerical methods/techniques to solve the modelled problem. Numerical techniques are numerous and the types know no bounds. They include; The Euler method, Runge-Kutta meth- ods, linear Multistep method, shooting method, Finite differ- ence method, finite element methods, e.g Galerkin method, 26 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 27 Spectral element method. Other methods are; Spectral method base on Fourier transformation; Method of lines reduces the PDE to a large system of the ordinary differential equation (ODE) Boundary Element Method (BEM) based on transforming a PDE to an integral equation on the boundary of the domains and it is popular in computational fluid dynamics, the list is endless. Authors who have worked extensively on numerical methods for approximation of solution of a differential equation include but no limited to Brugnano and Trigiante [1], Onumanyi et al. [2, 3], Jator [4, 5], Fatunla [6], Yusuph and Onumanyi [7], Siraj-ul-Islam, et al.[8], Adewale et al. [9]. We adopt the method of lines approach which is commonly used for solving time-dependent partial differential equations (PDEs), whereby the spatial derivatives are replaced by finite difference approximations see Ngwane and Jator [10]. The Method of Lines (MOLs) allows the conversion of PDEs into ODEs by complete or partial discretisation of the independent variables resulting in algebraic equations. If partial discretization is carried out and with only one remaining independent variable, then this results in the system of ODEs which is an approximation of the original PDE. Thus, one of the underscored features of the MOL is the use of existing, and generally well established, numerical methods for ODEs, for more literature on this approach see Brugnano and Trigiante [1], Ramos and Vigo-Aguiar [11], Cash [12] and Jator and Li [13]. 2. Derivation of the Method Second order ordinary differential equation of the form is considered y′′ = f (x, y) (1) subject certain conditions, where a, b are real numbers, f is a continuous function on (a, b) and y ∈ C2[a, b]. A 2-step block methods for the problem of the form (1) is considered. The grid points given by xn, xn+1 = xn + h, xn+2 = xn + 2h, are considered for solving the problem in (1) on the interval [xn, xn+2]. We assume a trial solution y(x) of (1) by a polynomial p(x) given by y(x) ' p(x) = m−1∑ i=0 ai x i (2) which on differentiating yields y′′(x) ' p′′(x) = m−1∑ i=2 i(i − 1)ai x i−2 (3) with the ai ∈ R real unknown parameters to be determined. and m = r + s; r is the number of interpolation points and s is the number of collocation points. 2.1. Specification of the method In this work the interval of integration considered is [xn, xn+2], we thus consider two different categories of off-set points viz-a- viz the points x i 3 , for i = 1, 2, 4, 5 and x i 4 , for i = 1, 2, 3, 5, 6, 7. 2.1.1. Case 1 Here, we consider the specification where the off-set points are x i 3 , for i = 1, 2, 4, 5. interpolating (2) at the points x i 3 , for i = 1, 2 implies r = 2 and collocating (3) at points x i 3 , for i = 0(1)6 implies s = 7 so that (2) and (3) becomes y(x) ' p(x) = 8∑ i=0 ai x i = a0 + a1 x + a2 x 2 + · · · + a8 x 8 (4) which on differentiating twice yields y′′(x) ' p′′(x) = 8∑ i=2 i(i − 1)ai x i−2 = 2a2 + 6a3 x + 12a4 x 2 + · · · + 56a8 x 6 (5) 27 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 28 From the imposed collocation condition, the following system of algebraic equation is obtained A =  1 xn+ 13 x2 n+ 13 x3 n+ 13 x4 n+ 13 x5 n+ 13 x6 n+ 13 x7 n+ 13 x8 n 13 1 xn+ 23 x2 n+ 23 x3 n+ 23 x4 n+ 23 x5 n+ 23 x6 n+ 23 x7 n+ 23 x8 n 23 0 0 2 6xn 12x2n 20x 3 n 30x 4 n 42x 5 n 56x 6 n 0 0 2 6xn+ 13 12x2 n+ 13 20x3 n+ 13 30x4 n+ 13 42x5 n+ 13 56x6 n+ 13 0 0 2 6xn+ 23 12x2 n+ 23 20x3 n+ 23 30x4 n+ 23 42x5 n+ 23 56x6 n+ 23 0 0 2 6xn+1 12x2n+1 20x 3 n+1 30x 4 n+1 42x 5 n+1 56x 6 n+1 0 0 2 6xn+ 43 12x2 n+ 43 20x3 n+ 43 30x4 n+ 43 42x5 n+ 43 56x6 n+ 43 0 0 2 6xn+ 53 12x2 n+ 53 20x3 n+ 53 30x4 n+ 53 42x5 n+ 53 56x6 n+ 53 0 0 2 6xn+2 12x2n+2 20x 3 n+2 30x 4 n+2 42x 5 n+2 56x 6 n+2  , x=  a0 a1 a2 a3 a4 a5 a6 a7 a8  ;b=  yn+ 13 yn+ 23 h2 fn h2 fn+ 13 h2 fn+ 23 h2 fn+1 h2 fn+ 43 h2 fn+ 53 h2 fn+2  e= (1,x, x2, x3, x4, x5, x6, x7, x8)T Where yn ≈ y (xn) , fn ≈ f (xn, yn, y′n). Hence, we state the following theorem without proof. Theorem 2.1. [14] Let (4) and (5) be satisfied, then the 2-step continuous linear hybrid multistep method is equivalent to the equation y(x) =bT (A−1k ) T e (6) where b, A and e are as defined above. Applying the above theorem, the following continuous hybrid method is derived y(x) = 2∑ i=1 α i 3 yn+ i3 +h 2 2∑ j=0 β i 3 fn+ i3 (7) where α and β are function of t given as α 1 3 = −3t+2, α 2 3 = 3t−1 β0=h2 [ 21 32 t 6− 459 4480 t− 27 160 t 7+ 12 t 2+ 814480 t 8+ 863108894− 49 40 t 3− 441 320 t 5+ 203120 t 4 ] β 1 3 =h2 [ −243 2240 t 8+ 2728 t 7− 279 80 t 6+ 26140 t 5− 261 40 t 4+3t3−2762360480 t+ 8999 90720 ] β 2 3 =h2 [ 243 896 t 8− 513 224 t 7+ 1233160 t 6− 4149 320 t 5+ 35132 t 4− 15 4 t 3+ 18689120960 t− 769 181440 ] β1=h2 [ −81 224 t 8+ 8128 t 7− 363 40 t 6+ 27920 t 5− 127 12 t 4+ 103 t 3− 139 864 t+ 1987 136080 ] β 4 3 =h2 [ 243 896 t 8− 459 224 t 7+ 963160 t 6− 2763 320 t 5+ 9916 t 4− 15 8 t 3+ 10921120960 t− 1609 181440 ] β 5 3 =h2 [ −243 2240 t 8+ 2735 t 7− 171 80 t 6+ 11740 t 5− 81 40 t 4+ 35 t 3− 347 12096 t+ 263 90720 ] β2=h2 [ 81 4480 t 8− 27 224 t 7+ 51160 t 6− 27 64 t 5+ 137480 t 4− 1 12 t 3+ 479120960 t− 221 544320 ] (8) where t= x−xnh , Evaluating (7) at non-interpolating points i.e at the points x=xn+ i3 for i= 0, 3, 4, 5, 6 which is equivalent to t = i 3 yn= 2yn+ 13 −yn+ 23 +h 2 ( 863 108864 fn+ 8999 90720 fn+ 13 − 769 181440 fn+ 23 + 1987 136080 fn+1 − 1609 181440 fn+ 43 + 263 90720 fn+ 53 − 221 544320 fn+2 ) yn+1= −yn+ 13 +2yn+ 23 +h 2 ( −221 544320 fn+ 977 90720 fn+ 13 + 16451 181440 fn+ 23 + 1357 136080 fn+1 + 71181440 fn+ 43 − 31 90720 fn+ 53 + 31 544320 fn+2 ) yn+ 43 = −2yn+ 13 +3yn+ 23 +h 2 ( −137 181440 fn+ 209 10080 fn+ 13 + 433 2240 fn+ 23 + 4927 45360 fn+1 + 25720160 fn+ 43 − 1 672 fn+ 53 + 31 181440 fn+2 ) yn+ 53 = −3yn+ 13 +4yn+ 23 +h 2 ( −19 181440 fn+ 17 560 fn+ 13 + 2987 10080 fn+ 23 + 4927 22680 fn+1 + 3893360 fn+ 43 + 41 5040 fn+ 53 − 11 90720 fn+2 ) yn+2= −4yn+ 13 +5yn+ 23 +h 2 ( −95 54432 fn+ 389 9072 fn+ 13 + 7085 18144 fn+ 23 + 4633 13608 fn+1 − 3893 18144 fn+ 43 + 1061 9072 fn+ 53 + 409 54432 fn+2 ) (9) 28 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 29 Differentiating α 1 3 , α 2 3 and all the β i 3 , i = 0(1)6 and evaluating the derivative of (8) at the points x=xn+ i3 ,i = 0(1)6, equivalent to t= i3 . hy′n+3yn+ 13 −3yn+ 23 = h 2 [ −459 4480 fn− 27623 60480 fn+ 13 + 18689 120960 fn+23− 139 864 fn+1+ 10921 120960 fn+ 43 − 347 12096 fn+ 53 + 479 120960 fn+2 ] hy ′ n+ 13 +3yn+ 13 −3yn+ 23 = h 2 [ 199 72576 fn− 1973 20160 fn+ 13 − 18 128 fn+23 + 4157 90720 fn+1− 851 40320 fn+ 43 + 416720 fn+ 53 + 289 120960 fn+2 ] y ′ n+ 23 +3yn+ 13 −3yn+ 23 = h 2 [ −731 362880 fn− 13 320 fn+ 13 + 6347 40320 fn+23− 3971 90720 fn+1+ 257 13440 fn+ 43 − 109 20160 fn+ 53 + 253 262880 fn+2 ] hy ′ n+1+3yn+ 13 −3yn+ 23 = h 2 [ −1 1920 fn+ 1537 60480 fn+ 13 + 39587 120960 fn+23 + 4927 30240 fn+1− 2201 120960 fn+ 43 + 2096720 fn+ 53 − 43 120960 fn+2 ] hy ′ n+ 43 +3yn+ 13 −3yn+ 23 = h 2 [ −571 362880 fn+ 691 20160 fn+ 13 + 1299 4480 fn+23 + 33533 90720 fn+1+ 1223 8064 fn+ 43 + 796720 fn+ 53 + 59 51840 fn+2 ] hy ′ n+ 53 +3yn+ 13 −3yn+ 23 = h 2 [ −29 362880 fn+ 51 2240 fn+ 13 + 13313 40320 fn+23 + 5081 18144 fn+1− 5519 13440 fn+ 43 + 73576 fn+ 53 + 1313 362880 fn+2 ] (10) The schemes in (9) and (10) form the requited method for solving (1) numerically. 2.1.2. Case 2 Here, we consider the specification where the off-set points are also x i 3 , for i= 1, 2, 4, 5. Interpolating (2) at the points xi+n, for i= 0, 1 implies r= 2 and collocating (3) at points x i 3 , for i = 0(1)6 implies s= 7 so that (2) and (3) becomes y(x)'p(x) = 8∑ i=0 ai x i =a0+a1 x+a2 x 2 +· · ·+a8 x 8 (11) which on differentiating yields y′′(x)'p′′(x) = 8∑ i=2 i(i−1)ai x i−2 = 2a2+6a3 x+12a4 x 2 +· · ·+56a8 x 6 (12) From the imposed collocation condition, the following system of algebraic equation is obtained A =  1 xn x2n x 3 n x 4 n x 5 n x 6 n x 7 n x 8 n 1 xn+1 x2n+1 x 3 n+1 x 4 n+1 x 5 n+1 x 6 n+1 x 7 n+1 x 8 n+1 0 0 2 6xn 12x2n 20x 3 n 30x 4 n 42x 5 n 56x 6 n 0 0 2 6xn+ 13 12x2 n+ 13 20x3 n+ 13 30x4 n+ 13 42x5 n+ 13 56x6 n+ 13 0 0 2 6xn+ 23 12x2 n+ 23 20x3 n+ 23 30x4 n+ 23 42x5 n+ 23 56x6 n+ 23 0 0 2 6xn+1 12x2n+1 20x 3 n+1 30x 4 n+1 42x 5 n+1 56x 6 n+1 0 0 2 6xn+ 43 12x2 n+ 43 20x3 n+ 43 30x4 n+ 43 42x5 n+ 43 56x6 n+ 43 0 0 2 6xn+ 53 12x2 n+ 53 20x3 n+ 53 30x4 n+ 53 42x5 n+ 53 56x6 n+ 53 0 0 2 6xn+2 12x2n+2 20x 3 n+2 30x 4 n+2 42x 5 n+2 56x 6 n+2  , x=  a0 a1 a2 a3 a4 a5 a6 a7 a8  ;b=  yn yn+1 h2 fn h2 fn+ 13 h2 fn+ 23 h2 fn+1 h2 fn+ 43 h2 fn+ 53 h2 fn+2  e= (1,x, x2, x3, x4, x5, x6, x7, x8)T Applying the above theorem, the following continuous hybrid method is derived y(x) = 1∑ i=0 αiyn+i+h 2 2∑ i=0 β i 3 fn+ i3 (13) 29 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 30 where α and β are function of t given as α0= −t α1= 1+t β0= 105h2 t−112h2 t3 +56h2 t4 +378h2 t5−252h2 t6−324h2 t7 +243h2 t8 13440 β 1 3 = − 3(−85h2 t−56h2 t3 +42h2 t4 +168h2 t5−168h2 t6−72h2 t7 +81h2 t8) 2240 β 2 3 = 3(347h2 t−560h2 t3 +840h2 t4 +546h2 t5−1092h2 t6−180h2 t7 +405h2 t8) 4480 β1= 563h2 t+1680h2 t2−3430h2 t4 +3528h2 t6−1215h2 t8 3360 β 4 3 = 3(−41h2 t+560h2 t3 +840h2 t4−546h2 t5−1092h2 t6 +180h2 t7 +405h2 t8) 4480 β 5 3 = − 3(−5h2 t+56h2 t3 +42h2 t4−168h2 t5−168h2 t6 +72h2 t7 +81h2 t8) 2240 β2= −11h2 t+112h2 t3 +56h2 t4−378h2 t5−252h2 t6 +324h2 t7 +243h2 t8 13440  (14) where t= x−xn−1h . Evaluating (13) at the points x=xn+ i3 for i= 1, 2, 4, 5, 6 which is equivalent to t= −2/3,−1/3, 1/3, 2/3, 1 the following main methods are obtained yn+ 13 = 2yn 3 + yn+1 3 −h 2 ( 2803 fn 544320− 1265 f n+ 13 18144 − 1657 f n+ 23 60480 − 1777 fn+1 136080 + 1049 f n+ 43 181440 − 11 f n+ 53 6048 + 137 fn+2 544320 ) yn+ 23 = yn 3 + 2yn+1 3 −h 2 ( 1291 fn 544320− 1217 f n+ 13 30240 − 10711 f n+ 23 181440 − 1567 fn+1 136080 + 163 f n+ 43 60480 − 67 f n+ 53 90720 + 53 fn+2 544320 ) yn+ 43 = − yn 3 + 4yn+1 3 +h 2 ( 661 fn 272160 + 1789 f n+ 13 45360 + 2147 f n+ 23 30240 + 6817 fn+1 68040 + 841 f n+ 43 90720 − f n+ 53 15120− 11 fn+2 272160 ) yn+ 53 = − 2yn 3 + 5yn+1 3 +h 2 ( 535 fn 108864 + 475 f n+ 13 6048 + 5167 f n+ 23 36288 + 5725 fn+1 27216 + 1321 f n+ 43 12096 + 193 f n+ 53 18144 − 53 fn+2 108864 ) yn+2= −yn+2yn+1+h2 ( 47 fn 6720 + 27 224 fn+ 13 + 459 f n+ 23 2240 + 563 fn+1 1680 + 459 f n+ 43 2240 + 27 224 fn+ 53 + 47 fn+2 6720 ) (15) differentiating (14) and evaluating the derivative of (13) at the points x=xn+ 3i for i= 0(1)6 which is equivalent to t= −1,−2/3,−1/3, 0, 1/3, 2/3, 1 the following additional methods are obtained y′n= − y h n + yn h +1 −h ( 253 fn 2688 − 165 448 fn+ 13 + 267 f n+ 23 4480 − 5 32 fn+1+ 363 f n+ 43 4480 − 57 f n+ 53 2240 + 47 fn+2 13440 ) , y′n+ 13 = − yn h + yn+1 h +h ( 4019h fn 362880 − 571h f n+ 13 60480 − 679h f n+ 23 3456 + 4577h fn+1 90720 − 3673h f n+ 43 120960 + 113h f n+ 53 12096 − 457h fn+2 362880 ) , y′n+ 23 = − yn h + yn+1 h +h ( 2293h fn 362880 + 223h f n+ 13 1728 + 7561h f n+ 23 120960 − 3551h fn+1 90720 + 1193h f n+ 43 120960 − 131h f n+ 53 60480 + 17h fn+2 72576 ) , y′n+1= − yn h + yn+1 h +h ( h fn 128 + 51 448 h fn+ 13 + 1041h f n+ 23 4480 + 563h fn+1 3360 − 123h f n+ 43 4480 + 3 448 h fn+ 53 − 11h fn+2 13440 ) , y′n+ 43 = − yn h + yn+1 h +h ( 2453h fn 362880 + 7421h f n+ 13 60480 + 23593h f n+ 23 120960 + 33953h fn+1 90720 + 3445h f n+ 43 24192 − 103h f n+ 53 12096 + 7h fn+2 10368 ) , y′n+ 53 = − yn h + yn+1 h +h ( 599h fn 72576 + 1345h f n+ 13 12096 + 28459h f n+ 23 120960 + 5165h fn+1 18144 + 48551h f n+ 43 120960 + 1123h f n+ 53 8640 − 1481h fn+2 362880 ) , y′n+2= − yn h + yn+1 h +h ( 47h fn 13440 + 327h f n+ 13 2240 + 111 896 h fn+ 23 + 1651h fn+1 3360 + 93 640 h fn+ 43 + 219 448 h fn+ 53 + 453h fn+2 4480 ) (16) Here, it is noteworthy that (9) and (10) are combined to form a block for case 1 while (15) and (16) form another block for case 11. For each case, (1) is solved which is a system of second-order ordinary differential equations resulting from the semi- discretization of a second-order PDE. 3. Analysis of the Method 3.1. Order and Local Truncation Error (LTE) The LMMs (8), and (13) are said to be of order p if C0=C1=C2= · · · =C p+µ−1 = 0, C p+µ,0. 30 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 31 Here C p+µ is the error constant and C p+µh p+µy( p+µ)(xn) is the principal Local Truncation Error (LTE) at the point xn. The C′s are given by C0=α0+α1+α2+· · ·+αk C1= (α0+α1+α2+· · ·+αk) − (β0+β1+· · ·+βk) Cq= 1 q! (α1+2 qα2+· · ·+k qαk)− 1 (q−3)! (β1+2 q−1β2+· · ·+k q−3βk),q= 2, 3, . . . The LTE associated with any of (8) and (13) is given by the difference operator L[y(x) :h] = 2∑ i=1 α i 3 y(xn+ i 3 )−h2 2∑ j=0 β i 3 y′′(xn+ i 3 ) (17) where y∈C2[a, b] is an arbitrary function. Expanding (17) in Taylor’s series about the point xn, the expression is obtained as: L[y(xn) :h] =C0y(xn)+C1hy ′(xn)+C2h 2y′′(xn)+· · ·+Cρ+2h ρ+2yρ+2(xn) (18) Expanding each scheme in (9) and (10), the following principal truncation errors are obtained: C0p+2=− 19y(9)[x]h9 119042784 +O[h]10, C1p+2= 31y(9)[x]h9 1190427840 +O[h]10, C2p+2= y(9)[x]h9 4723920 +O[h]10 C 4 3 p+2= 31y(9)[x]h9 595213920 +O[h]10, C 5 3 p+2= 31y(9)[x]h9 595213920 +O[h]10, C′0 p+2= 8881y(9)[x]h9 5952139200 +O[h]10 C′ 1 3 p+2= − 409y(9)[x]h9 1700611200 +O[h]10, C′ 2 3 p+2= 1201y(9)[x]h9 5952139200 +O[h]10, C′1 p+2= − 463y(9)[x]h9 11904278400 +O[h]10 C′ 4 3 p+2= 1201y(9)[x]h9 5952139200 +O[h]10, C′ 5 3 p+2= − 409y(9)[x]h9 1700611200 +O[h]10, C′2 p+2= 8881y(9)[x]h9 5952139200 +O[h]10 The above blocked method (9) and (10) is of uniform order p= 7 Expanding each scheme in (15) and (16), the following principal truncation errors are obtained: C 1 3 p+2= 349y(9)(x)h9 3571283520 +O(h)10, C 2 3 p+2= 2y(9)(x)h9 55801305 +O(h)10, C 4 3 p+2= − 2y(9)(x)h9 55801305 +O(h)10 C 5 3 p+2= − 349y(9)(x)h9 3571283520 +O(h)10, C2p+2= y(9)(x)h9 32659200 +O(h)10, C′0p+2= y(9)(x)h9 765450 +O(h)10 C′ 1 3 p+2= − 1691y(9)(x)h9 3968092800 +O(h)10, C′ 2 3 p+2= y(9)(x)h9 62001450 +O(h)10, C′1p+2− 11y(9)(x)h9 48988800 +O(h)10 C′ 4 3 p+2 y(9)(x)h9 62001450 +O(h)10, C′ 5 3 p+2− 1691y(9)(x)h9 3968092800 +O(h)10, C′2p+2 y(9)(x)h9 765450 ++O(h)10 The above blocked method (15) and (16) is of uniform order p= 7 The LMM (8) (same for (13)) is said to be consistent if it has order p≥1 and the first and second characteristic polynomials which are defined respectively, as ρ(r) = k∑ j=0 α jz j (19) and σ(r) = k∑ j=0 β jz j (20) 31 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 32 where r is the principal root, satisfy the following conditions: k∑ j=0 α j= 0 (21) ρ(1) =ρ′(1) = 0 (22) and ρ′′(1) = 2!σ(1) (23) Henrici [15], Lambert[16]. Consider the main method in (9) given as yn+2= −4yn+ 13 +5yn+ 23 +h 2 ( −95 54432 fn+ 389 9072 fn+ 13 + 7085 18144 fn+ 23 + 4633 13608 fn+1 − 3893 18144 fn+ 43 + 1061 9072 fn+ 53 + 409 54432 fn+2 ) (24) The condition (21) is satisfied. the first characteristic equation for (9) is given as: ρ(r) =r2+4r 1 3 −5r 2 3 (25) ρ′(r) = 2r+ 4 3r2/3 − 10 3r1/3 (26) Here ρ(r) = 0, ρ′(r) = 0. Therefore, (22) is satisfied. The second characteristic polynomial for (9) is given as σ(r) = −95 54432 + 389 9072 r 1 3 + 7085 18144 r 2 3 + 4633 13608 r+ 3893 18144 r 4 3 + 1061 9072 r 5 3 + 409 54432 r2 (27) σ(1) = 10 9 (28) ρ′′(r) = 2− 8 9r5/3 + 10 9r4/3 . ρ′′(r) = 20 9 (29) Hence condition (23) is satisfied. Conclusively, the hybrid method is consistent. Consider the main method in (15) given as yn+2= −yn+2yn+1+h2 ( 47 fn 6720 + 27 224 fn+ 13 + 459 f n+ 23 2240 + 563 fn+1 1680 + 459 f n+ 43 2240 + 27 224 fn+ 53 + 47 fn+2 6720 ) (30) The condition (21) is satisfied. the first characteristic equation for (15) is given as: ρ(r) =r2−2r+1 (31) ρ′(r) = 2r−2 (32) Here ρ(1) = 0, ρ′(1) = 0. Therefore, (22) is satisfied. The second characteristic polynomial for (15) is given as σ(r) = 47 6720 + 27 224 r 1 3 + 459 2240 r 2 3 + 563 1680 r+ 459 2240 r 4 3 + 27 224 r 5 3 + 47 6720 r2 (33) σ(1) = 1 (34) ρ′′(r) = 2. ρ′′(1) = 2 (35) Hence condition (23) is satisfied. Conclusively, the hybrid method is consistent. 32 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 33 3.2. Zero Stability To establish hat the methods are zero stable, each of the method in block form are solved simultaneously to obtain all the yi and y′i’s for appropriate index i, (see Modebei et al.[17]). For the method (9) and its additional methods in (10), they are taken in block form and solved simultaneously to obtain y i 3 for i = 1(1)6 to obtain the following block method. For block method (9)-(10) yn+ 13 = yn+ hy′n 3 +h 2 ( 28549 fn 1088640 + 275 f n+ 13 5184 − 5717 f n+ 23 120960 + 10621 fn+1 272160 − 7703 f n+ 43 362880 + 403 f n+ 53 60480 − 199 fn+2 217728 ) yn+ 23 = yn+ 2hy′n 3 +h 2 ( 1027 fn 17010 + 194 945 fn+ 13 − 8 81 fn+ 23 + 788 fn+1 8505 − 97 f n+ 43 1890 + 46 f n+ 53 2835 − 19 fn+2 8505 ) yn+1= yn+hy′n+h 2 ( 253 fn 2688 + 165 448 fn+ 13 − 267 f n+ 23 4480 + 5 fn+1 32 − 363 f n+ 43 4480 + 57 f n+ 53 2240 − 47 fn+2 13440 ) yn+ 43 = yn+ 4hy′n 3 +h 2 ( 1088 fn 8505 + 1504 f n+ 13 2835 − 8 945 fn+ 23 + 2624 fn+1 8505 − 8 81 fn+ 43 + 32 945 fn+ 53 − 8 fn+2 1701 ) yn+ 53 = yn+ 5hy′n 3 +h 2 ( 35225 fn 217728 + 8375 f n+ 13 12096 + 3125 f n+ 23 72576 + 25625 fn+1 54432 − 625 f n+ 43 24192 + 275 f n+ 53 5184 − 1375 fn+2 217728 ) yn+2= yn+2hy′n+h 2 ( 41 fn 210 + 6 7 fn+ 13 + 3 35 fn+ 23 + 68 fn+1 105 + 3 70 fn+ 43 + 6 35 fn+ 53 ) y′n+ 13 = y ′ n+h ( 19087 fn 181440 + 2713 f n+ 13 7560 − 15487 f n+ 23 60480 + 586 fn+1 2835 − 6737 f n+ 43 60480 + 263 f n+ 53 7560 − 863 fn+2 181440 ) y′n+ 23 = y ′ n+h ( 1139 fn 11340 + 94 189 fn+ 13 + 11 f n+ 23 3780 + 332 fn+1 2835 − 269 f n+ 43 3780 + 22 945 fn+ 53 − 37 fn+2 11340 ) y′n+1= y ′ n+h ( 137 fn 1344 + 27 56 fn+ 13 + 387 f n+ 23 2240 + 34 fn+1 105 − 243 f n+ 43 2240 + 9 280 fn+ 53 − 29 fn+2 6720 ) y′n+ 43 = y ′ n+h ( 286 fn 2835 + 464 945 fn+ 13 + 128 945 fn+ 23 + 1504 fn+1 2835 + 58 945 fn+ 43 + 16 945 fn+ 53 − 8 fn+2 2835 ) y′n+ 53 = y ′ n+h ( 3715 fn 36288 + 725 f n+ 13 1512 + 2125 f n+ 23 12096 + 250 fn+1 567 + 3875 f n+ 43 12096 + 235 f n+ 53 1512 − 275 fn+2 36288 ) y′n+2= y ′ n+h ( 41 fn 420 + 18 35 fn+ 13 + 9 140 fn+ 23 + 68 fn+1 105 + 9 140 fn+ 43 + 18 35 fn+ 53 + 41 fn+2 420 ) (36) For block method (15)-(16) similar operation is carried out; A numerical method is zero-stable if the solutions remain bounded as h→0, which means that the method does not provide solutions that grow unbounded as the number of steps increases, Modebei et al.[17]. To show the zero-stability of the block method (36), we take h→0 the method may be rewritten in matrix form as A0Yn=A1Yn−1 (37) Yn= (Y 0 n , Y 1 n ) T Y 0n = (yn+ 13 , yn+ 23 , yn+1, yn+ 43 , yn+ 53 , yn+2) T Y 1n = (y ′ n+ 13 , y′n+ 23 , y ′ n+1, y ′ n+ 43 , y′n+ 53 , y ′ n+2) T For method (36) A0=I12×12 identity matrix and A1=I12×12 matrix given by A1=  A11 0 0 A22  , A11=  1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0  , A22=  0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0  33 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 34 The characteristic polynomial of the matrix A11 is given as |A11−λI|, that is λ5(λ−1)= 0 with root λ j= 0 for j= 1, . . . , 5 and λ6= 1. The characteristic polynomial of the matrix A22 is given as |A22−λI|, that is λ5(λ−1) = 0 with root λ j= 0 for j= 1, . . . , 5 and λ6= 1. For method (15)-(16) A0=I12×12 identity matrix and A1=I12×12 matrix given by A1=  A11 0 0 A22  , A11=A22=  1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0  The characteristic polynomial of the matrix Aii is given as |Aii−λI|, i= 1, 2, that is λ5(λ−1) = 0 with root λ j= 0 for j= 1, . . . , 5 and λ6= 1. Definition 3.1. The two step hybrid block method (9)-(10) (or (15)-(16)) is said to be zero stabile if the number of root of the first characteristic equation |ρ(r)| < 1 and if |ρ(r)| = 1, then the multiplicity of ρ(r) must not exceed 2. Hence, the are zero stable. 3.3. Convergence of the Methods Definition 3.2. Convergence: An LMM is said to be convergent if and only if it is consistent and zero-stable. By the above definition, the derived hybrid methods are convergent. 4. Numerical Examples In this section, the performance of the developed two-step hybrid block scheme is examined. the exact and approximate solution are tabulated. The tables below show the numerical re- sults of the newly developed scheme with the exact solution for solving the problem and the result of the developed scheme are more accurate than existing methods. For simplicity, method in (9)-(10) would be termed Hybrid 2-step Block Method 1 (H2BM1), and method in (15)-(16) would be termed Hybrid 2-step Block Method 2 (H2BM2) Example 1. Consider the PDE (Ngwane and Jator [10]). κ ∂y2 ∂x2 − ∂y ∂t = 0 y(0,t) =y(1,t) = 0, y(x, 0) = sinπx+sinωπx, κ>1 (38) The analytic solution is y(x, t) =e−π 2κtsinπx+e−ω 2π2κtsinπx Following Ngwane and Jator [10], (38) becomes dy2m dx2 = 1 κ y(x , tm+1)−y(x , tm−1)] (∆t) (39) ym(0,tm) =ym(1,tm) = 0, ym(x, 0) = sinπx+sinωπx, κ>1 Table 1: Exact and Numerical solution for Example 1 x Exact H2BM1 Error 0.0 1.65341E-9 1.65341E-9 0 0.1 6.16242E-10 6.16242E-10 2.78E-12 0.2 2.29678E-10 2.29678E-10 4.45E-12 0.3 8.56029E-11 8.56029E-11 4.20E-12 0.4 3.19048E-11 3.19048E-11 5.00E-12 0.5 1.18911E-11 1.18911E-11 6.45E-12 0.6 4.43194E-12 4.43194E-12 7.31E-12 0.7 1.65181E-12 1.65181E-12 3.28E-12 0.8 6.15646E-13 6.15646E-13 4.11E-13 0.9 2.29456E-13 2.29456E-13 5.99E-13 Table 2: Exact and Numerical solution for Example 1 9 Exact H2BM2 Error 0.0 1.65341E-9 1.65341E-9 0 0.1 6.16242E-10 6.16242E-10 3.24E-12 0.2 2.29678E-10 2.29678E-10 7.21E-12 0.3 8.56029E-11 8.56029E-11 2.22E-12 0.4 3.19048E-11 3.19048E-11 8.45E-12 0.5 1.18911E-11 1.18911E-11 1.15E-12 0.6 4.43194E-12 4.43194E-12 1.18E-12 0.7 1.65181E-12 1.65181E-12 5.08E-12 0.8 6.15646E-13 6.15646E-13 2.48E-13 0.9 2.29456E-13 2.29456E-13 4.79E-13 where tm=m∆t, m= 0, 1, . . . ,; M ym(x)≈y(x, tm),y(x) = [y0(x), y1(x), . . . ,yM−1(x)] T , hence (39) becomes the system d 2 ym (x) dx2 = f (x , ym) which is in the form of (1), where f (x, tm) =Ay+G and A is an M−1 square matrix, G is a vector of constants. BHSDA is L -Stable Block Hybrid Second Derivative Al- gorithm in Ngwane and Jator [10]. Tables 1 and 2 shows the comparison of exact solution and the mothers H2BM1 and H2BM2 respectively. For example 1. Table 3 shows the comparison of maximum errors obtained for example 1 using the derived methods and the method in Ngwane and Jator[10]. This shows the superiorly of the derived methods over existing methods. Figure 1 show the surface plots for the exact solution and Numerical solutions for example 1. Example 2. Consider the PDE ([14]): ∂y2 ∂x2 + ∂y2 ∂t2 = −32π 2sin(4πx), x∈[0, 1] y(±1,t) =y(x,±1) = 0, t>0 (40) Table 3: Comparison of maximum errors obtained in different methods for Ex- ample 1 at t= 1. κ H2BM1 H2BM2 BHSDA 1 1.076 × 10−11 1.022 × 10−12 2.64 × 10−6 2 1.024 × 10−12 1.085 × 10−12 1.32 × 10−6 3 1.045 × 10−12 1.045 × 10−12 1.32 × 10−6 5 1.035 × 10−12 1.055 × 10−12 1.32 × 10−6 10 1.019 × 10−12 1.041 × 10−12 1.32 × 10−6 34 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 35 Figure 1: Surface plots for the Exact and numerical solution for Example 2 Table 4: Exact and Numerical solution using H2BM1 for Example 2. κ H2BM1 H2BM2 BHSDA x Exact H2BM1 Error 0.0 6.634320126E-16 6.634320126E-16 0. 0.2 0.5590169122545 0.5590169122520 2.51E-12 0.4 -0.9045084378512 -0.9045084378511 1.00E-12 0.6 0.9045084354472 0.9045084354462 1.00E-12 0.8 -0.5590169311454 -0.5590169311421 3.37E-12 1.0 -4.658833273E-16 -4.658833273E-16 0. for N= 10, x∈[−1, 1] Table 5: Exact and Numerical solution using H2BM2 for Example 2. x Exact H2BM2 Error 0.0 6.634320126E-16 6.634320126E-16 0. 0.2 0.5590169122545 0.5590169122589 4.49E-12 0.4 -0.9045084378512 -0.9045084378577 6.52E-12 0.6 0.9045084354472 0.9045084354428 5.53E-12 0.8 -0.5590169311454 -0.5590169311433 2.12E-12 1.0 -4.658833273E-16 -4.658833273E-16 0. for N= 10, x∈[−1, 1] The analytic solution is y(x, t) = sin4πxsin4πt. Following Ngwane and Jator [10], (38) becomes dy2m dx2 = − y (x , tm+1)−2y (x , tm) +y(x , tm−1)] (2∆t) −32π2sin(4πx) (41) ym(±1,tm) =ym(x,±1) = 0, where tm=m∆t, m= 0, 1, . . . ,M; ym(x)≈y(x, tm),y(x) = [y0(x), y1(x), . . . ,yM−1(x)] T , hence (41) becomes the system d 2 ym (x) dx2 = f (x , ym) which is in the form of (1), where f (x, tm) =Ay+G and A is an M−1 square matrix, G is a vector of constants. BVM and BUM are Boundary Value Methods and the Block Unification Methods in Biala [14]. Tables 4-5 shows the comparison of the exact solution and the methods H2BM1 and H2BM2 respectively for example 2. Table 6-7 shows the maximum error and CPU time obtained for different methods. Table 9 shows the maximum error and CPU time obtained in Biala [14]. Table 6: Comparison of maximum errors obtained in different methods for Ex- ample 2 at t= 1 N H2BM1 l∞ error CPU Time H2BM2 l∞ error CPU Time 16 2.257E-7 0.112 5.547E-7 0.121 32 2.787E-7 0.898 1.712E-7 0.871 64 8.234E-7 2.785 2.337E-7 2.662 128 3.114E-7 11.211 4.785E-7 12.009 256 2.779E-7 31.812 1.112E-7 30.101 Table 7: Comparison of maximum errors obtained in different methods for Ex- ample 2 at t= 1 N BVM l∞ error CPU Time BUM l∞ error CPU Time 16 9.662E-0 0.483 1.251E-1 0.531 32 2.582E-2 1.235 2.578E-2 1.031 64 6.433E-3 5.358 6.459E-3 5.516 128 1.607E-3 43.641 1.607E-3 46.923 256 2.000E-0 512.843 4.016E-4 532.657 This show that the derived methods performs accurately, su- periorly and affluently in terms of the computer time, and errors obtained in examples 2. Figure 2 shows the surface plots for the exact and Numerical solution for examples 2. Example 3. We consider the PDE (Jator [15]). ∂y2 ∂t2 + ∂y2 ∂x2 = sin(y), x∈[−3, 3] y(x, 0) = 4arctan(e x√ 1−c2 ), yt(x, 0) = − 4ce x√ 1−c2 √ 1−c2 (1+e 2 x√ 1−c2 ) , 0 < t<1 (42) The analytic solution is y(x, t) = 4arctan(sech(x)t), c is ve- locity of the wave. The problem is solved for c= 0.5, ∆t= 0.125 Following Ngwane and Jator [10], (38) becomes dy2m dx2 = − y (x , tm+1)−2y (x , tm) +y(x , tm−1)] (2∆t) +sin(ym)(43) ym (x, 0) = 4arctan ( e x√ 1−c2 ) , ymt (x, 0) = − 4ce x√ 1−c2 √ 1−c2 1+e2 x√1−c2 , 0 < t<1 35 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 36 Figure 2: Surface plots for the Exact and numerical solution for Example 3 Figure 3: Surface plots for the Exact and numerical solution for Example 3 Table 8: Exact and Numerical solution using different methods for Example 3 x Exact H2BM1 H2BM2 SBVM 0.125 0.12195641127 0.12195641122 0.12195641178 0.121956 0.25 0.11346954831 0.11346954852 0.11346954877 0.113469 0.375 0.10557254177 0.10557254144 0.10557254129 0.105573 0.5 0.09822454418 0.09822454411 0.09822454412 0.0982256 0.625 0.09138842135 0.09138842122 0.09138842129 0.0913892 0.75 0.08502738529 0.08502738541 0.08502738557 0.0850284 0.875 0.07910885264 0.07910885215 0.07910885213 0.0791101 1.00 0.07360212510 0.07360212548 0.07360212544 0.0736035 Table 9: Approximate and Numerical solution for Example 3 x H2BM1 Error H2BM2 Error SBVM Error 0.125 7.2E-10 5.1E-10 1.30E-7 0.25 2.1E-11 4.6E-10 2.94E-7 0.375 3.3E-11 4.8E-10 4.51E-7 0.5 7.2E-11 8.0E-10 6.49E-7 0.625 1.4E-10 5.5E-10 8.41E-7 0.75 1.2E-10 3.7E-10 1.03E-6 0.875 4.9E-10 3.3E-10 1.25E-6 1.00 3.8E-10 7.7E-10 1.44E-6 where tm=m∆t, m= 0, 1, . . . ,M; ym(x)≈y(x, tm),y(x) = [ y0(x), y1(x), . . . ,yM−1(x)] T , hence (43) becomes the system d 2 ym (x) dx2 = f (x , ym) which is in the form of (1). where f (x, tm) =Ay+G and A is an M−1 square matrix, G is a vector of constants. SBVM is symmetric boundary value method in Jator [15]. Table 8 shows the exact and numerical solution using the differ- ent methods for example 3 while Table 9 shows error obtained for example 3. Figure 3 shows the surface plots for the Exact and numerical solutions of H2BM1 and H2BM2 for Example 3 5. Conclusion The development of some numerical schemes has been pro- posed in this work. This was developed via the interpolation and collocation techniques using power series function as trial solutions. The methods were effectively illustrated some partial differential equations (PDE) and the results obtained were ac- curate. The analysis of the new methods showed that all satisfy the properties of numerical methods for the solution of differen- tial equations. Namely, Consistency, Zero- Stability, Continuity and convergence. References [1] L. Brugnano, L. & D. Trigiante, “Solving Differential Problems by Mul- tistep Initial and Boundary Value Methods”, Gordon and Breach Science Publishers, Amsterdam, (1998). [2] P. Onumanyi, U.W. Sirisena. & S.N Jator, “Continuous finite difference approximations for solving differential equations”, International. Journal of. Computer Mathematics. 72 (1999) 15. [3] P. Onumanyi, D. O. Awoyemi, S. N. Jator & U. W. Sirisena, “New linear multistep methods with continuous coefficients for first order initial value problems”, 19 Journal of Nigeria Mathematics Society 13, (1994) 37. [4] S. N. Jator “On the Hybrid Method with Three-Off Step Points for Ini- tial Value Problems”, International Journal of Mathematical Education in Science and Technology, 41 (2010) 110 36 Olaiya et al. / J. Nig. Soc. Phys. Sci. 3 (2022) 26–37 37 [5] S. N. Jator “Trigonometric symmetric boundary value method for oscillat- ing solutions including the sine-Gordon and Poisson equations”, Cogent Mathematics 3 (2016) 1271. [6] S. O. Fatunla, “Block Methods for second-order IVP”, International Jour- nal of Mathematics 55 (1999). [7] Y. Yusuph & P. Onumanyi; “New multiple FDM through multistep col- location” Proceedings of the conference, National Mathematical Centre, Abuja (2005). [8] Siraj-ul-Islam, Imran Aziz & Bozidar Sarler, “The numerical solution of second-order boundary-value problems by collocation method with the Haar wavelets”, Mathematical and Computer Modelling 52 (2010) 1577. [9] J. Adewale, A. Olaide. & J. Sunday, “Continuous Block Method for the Solution of second rder Initial Value Problems of Ordinary Diff. Equa- tion” International Journal of Pure and Applied mathematics 82 (2013) 405. [10] F. F. Ngwane & S. N. Jator, “L-Stable Block Hybrid Second Derivative Algorithm for Parabolic Partial Differential Equations”, American Jour- nal of Computational Mathematics 4 (2014) 87. [11] J. Vigo-Aguiar & H. Ramos, “A family of A-Stable Collocation Methods of Higher Order for Initial-Value Problems”, IMA Journal of Numerical Analysis 27 (2007) 798. [12] J. R, Cash “Two New Finite Difference Schemes for Parabolic Equa- tions”, SIAM Journal of Numerical Analysis 21 (1984) 433. [13] S. N. Jator & J. Li, “An algorithm for Second Order Initial and Boundary Value Problems” Numerical Algorithm (2012) 59. [14] T. A. Biala, “Computational Study of the Boundary Value Methods and the Block Unification Methods ”, Abstract and Applied Analysis (2012) Article ID 8465103, 14 pages [15] P. Henrici, “Discrete Variable Methods in Ordinary Differential Equa- tions”, John Wiley, New York, 1962. [16] J. D. Lambert, “Computational Methods in Ordinary Differential Equa- tions”, John Wiley and Sons, New York, NY, USA, 1973. [17] M. I. Modebei, R. B. Adeniyi, S. N. Jator & H. C. Ramos “A block hybrid integrator for numerically solving fourth-order Initial Value Problems”, Applied Mathematics and Computation 346 (2019) 680. 37