J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 Journal of the Nigerian Society of Physical Sciences Original Research A Class of Block Multi-derivative Numerical Techniques for Singular Advection Equations M. O. Ogunniran∗ Department of Mathematical Sciences, Osun State University Osogbo, Nigeria Abstract The integration of some differential equations is hard to acquire because of the presence of singular point(s) in these equations. These equations are best solved by some unique technique. Multi-derivative techniques have a long history of powerful integration of such equations yet till date, a couple of class of this technique has been explored for integrating partial differential equations. This work centers around the development, analysis and implementation of a class of multi-derivative technique on partial differential equations. The approaches were effectively analyzed and were turned out to be consistent, stable and convergent. Numerical outcomes got likewise demonstrated the approximation quality of the technique over existing techniques in literature. Keywords: Advection equations, Singular point, Multi-derivative, Conservation law, off-step point Article History : Received: 21 April 2019 Received in revised form: 03 June 2019 Accepted for publication: 04 June 2019 Published: 13 July 2019 c©2019 Journal of the Nigerian Society of Physical Sciences. All Rights Reserved. Communicated by: T. Latunde 1. Introduction Considering an equation which is known as the linear ad- vection equation for a quantity u(x, t): ∂u ∂t + v(x, t) ∂u ∂x = 0 (1) where v(x, t) is known, This type of equation in (1) is classified as first order hyperbolic partial differential equation. Equation (1) is termed as a conservation law since it expresses conser- vation of mass, energy or momentum under the conditions for which it is derived, i.e. the assumptions on which the equation is based. ∗Corresponding Author Tel. No: +2347039104606 Email address: muideen.ogunniran@uniosun.edu.ng (M. O. Ogunniran ) In physical phenomenon, v(x, t) is linear or flow velocity. Although equation (1) is possibly the simplest partial differen- tial equation, this simplicity is deceptive in the sense that it can be very difficult to integrate numerically due to the presence of singularity, a distinctive feature of first order hyperbolic partial differential equations. Equation (1) will have to be supplemented with an initial condition: u(x, 0) = u0(x) (2) and if necessary boundary condition(s): u(a, t) = u1(t), u(b, t) = u2(t). x ∈ [a, b] (3) One of the most popular techniques applied to the solution of Advection equation is finite difference method [1,2]. In recent times, the finite element and the finite volume methods were 62 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 63 also introduced. Many special methods have also been devel- oped to deal with Singularly Perturbed Problems. Many Au- thors [3,4,5,6] are devoted to Singularly Perturbed Problems of ordinary differential equation and the authors discussed the sit- uation and width of boundary layer(s) and give some effective numerical algorithms. In the work of Xianyl [7], a linear hybrid variable method was developed for the solution of Advection equations and an efficient numerical treatment of non-linearities in the advection-diffusion-reaction equations was presented by Utku et al. [8]. An extensive work was done by Gaetamo and Silvia [9] on the fractional diffusion-advection equations for fluids and plasmas. Adrian [10] presented a generalization of prefactored compact schemes for the solution of advection equations. 2. Method and Analysis The methods were derived by assuming a continuous ap- proximant for un(x) of a multi-step multi-derivative method of the form: un(x) = αn(x)yn + l∑ i=1 hi(βi,0(x) f (i−1) n + βi,v(x) f (i−1) n+v + k∑ j=1 β j,k(x) f (i−1) n+ j ) ≈ y(x) (4) for solving IVP: y′ = f (x, y), y(0) = y0, where f (x, y) is con- tinuous and jth differentiable, yn is an approximation to y(xn), xn = nh, h > 0 and f ( j) m = f ( j)(xm, ym) such that: f (0)(x, y) = f (x, y), f ( j)(x, y) = ∂ f ( j−1)(x, y) ∂x + f (x, y) ∂ f ( j−1)(x, y) ∂y , where k is the step number, l is derivative order, v ∈ (0, k) is an off-step point and j is the derivative order of f (x, y). However, the continuous coefficients βi, j(x), i = 0, v, k, j = 1(1)k were determined. To this end, approximation of the exact solution y(x) was sought by evaluating the function: u(x) = r+ls−1∑ j=0 a j x j (5) where a j, ( j = 0, 1, ..., r + ls−1) are coefficients determined, x j are the basis functions of degree r + ls − 1, l, r and s being the derivative order, interpolating and collocations points respec- tively. While ensuring that the function (4) corresponds with the an- alytical solution at the end point xn, the following conditions were imposed on u(x) and its derivatives u(k)(x) to get the coef- ficients of the desired methods: u(xn+ j) = yn+ j, j = 0 u′(xn+ j) = fn+ j, j ∈ [0, · · · , k] u′′(xn+ j) = f ′n+ j = gn+ j, j ∈ [0, · · · , k] u′′′(xn+ j) = f ′′n+ j = hn+ j, j ∈ [0, · · · , k] ... u(k)(xn+ j) = f (k−1) n+ j , j ∈ [0, · · · , k]  (6) The coefficients obtained were substituted into (5) to obtain the continuous coefficients in (4) which were then evaluated at xn+k to obtain the desired methods. 2.1. One-step, First Derivative and Two off-step Method Here, the first derivative evaluation was considered taking into consideration, one step and two off-step points. The proce- dure is as follows: For k = 1, l = 1 and v = 17 , 2 7 (4) becomes: u(x) = αn(x)yn + h(β1,0(x) fn + β1, 17 (x) fn+ 17 +β1, 27 (x) fn+ 27 + β1,1(x) fn+1) (7) while taking r = 1 and s = 4, (5) and its derivative become: u(x) = 4∑ j=0 a j x j (8) u′(x) = 4∑ j=1 ja j x j−1 (9) Imposing the following conditions obtained from (6): u(xn+ j) = yn+ j, j = 0. (10) u′(xn+ j) = fn+ j, j = 0, 1 7 , 2 7 , 1 (11) a system of five algebraic equations with five unknowns were obtained. The 5 × 5 system of equations obtained are given be- low: a0 = yn (12) a1 = fn (13) a1 + 2 7 a2 + 3 49 a3 + 4 343 a4 = fn+ 17 (14) a1 + 4 7 a2 + 12 49 a3 + 32 343 a4 = fn+ 27 (15) a1 + 2a2 + 3a3 + 4a4 = fn+1. (16) Solving the resulting equations gives: a0 = yn (17) a1 = fn+1 (18) a2 = − 23 4 fn + 49 6 fn+ 17 − 49 20 fn+ 27 + 1 30 fn+1 (19) a3 = 35 3 fn − 49 2 fn+ 17 + 196 15 fn+ 27 − 7 30 fn+1 (20) a4 = − 343 40 fn+ 27 + 49 120 fn+1 − 49 8 fn + 343 24 fn+ 17 (21) a4 = − 49 8 fn + 343 24 fn+ 17 − 343 40 fn+ 27 + 49 120 fn+1 (22) Substituting these values in (7) and collecting like terms in yn, fn, fn+ 17 , fn+ 27 , fn+1 yields the continuous method: y(x) = αn(x)yn + h(β0,1(x) fn + β1, 17 (x) fn+ 17 +β1, 27 (x) fn+ 27 + β1,1(x) fn+1) (23) 63 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 64 where: αn(x) = 1 β1,0(x) = − 49 (2 x−1)4 8 + 35 (2 x−1)3 3 − 23 (2 x−1)2 4 + 2 x − 1 β1, 17 (x) = 343 (2 x−1) 4 24 − 49 (2 x−1)3 2 + 49 (2 x−1)2 6 β1, 27 (x) = −343 (2 x−1) 4 40 + 196 (2 x−1)3 15 − 49 (2 x−1)2 20 β1,1(x) = 49 (2 x−1)4 120 − 7 (2 x−1)3 30 + 1/30 (2 x − 1) 2  (24) Evaluating (23) at xn+1 yields the one-step, first derivative main method possessing two off-grid points: yn+1 = yn + h 24 (19 fn − 49 fn+ 17 + 49 fn+ 27 + 5 fn+1) (25) Also, evaluating (23) at xn+ 17 and xn+ 27 yields the additional method as: yn+ 17 = yn + h 5880 (335 fn + 595 fn+ 17 − 91 fn+ 27 + fn+1) yn+ 27 = yn + h 21 ( fn + 4 fn+ 17 + fn+ 27 ) (26) The one-step, first derivative and two off-grid points method is presented in block form as follows: yn+1 = yn + h 24 (19 fn − 49 fn+ 17 + 49 fn+ 27 + 5 fn+1) yn+ 17 = yn + h 5880 (335 fn + 595 fn+ 17 − 91 fn+ 27 + fn+1) yn+ 27 = yn + h 21 ( fn + 4 fn+ 17 + fn+ 27 ) (27) 2.2. Two-step, First Derivative and Two off-step Method Here, a new member of the methods is obtained. The pro- cedure is as follows: For k = 2, l = 1, and v = 17 and 9 7 , (4) becomes: u(x) = αn(x)yn + h(β1,0(x) fn + β1, 17 (x) fn+ 17 +β1, 97 (x) fn+ 97 + β1,1(x) fn+1 + β1,2(x) fn+2) (28) taking r = 1 and s = 5, (5) and its derivative become: u(x) = 5∑ j=0 a j x j (29) u′(x) = 5∑ j=1 ja j x j−1 (30) Imposing the following conditions obtained from (6): u(xn+ j) = yn+ j, j = 0. (31) u′(xn+ j) = fn+ j, j = 0, 1 7 , 1, 9 7 , 2. (32) The two-step, first derivative and two off-grid point method is presented in block form as follows: yn+2 = yn + h(− 53 135 fn + 2401 2340 fn+ 17 + 11 4 fn+1 + 24012700 fn+ 97 + 227 975 fn+2) yn+ 17 = yn + h( 3391 52920 fn + 10579 131040 fn+ 17 − 41 8820 fn+1 + 439 151200 fn+ 97 − 79 382200 fn+2) yn+1 = yn + h(− 257 1080 fn + 14749 18720 fn+ 17 + 127 180 fn+1 − 5831 21600 fn+ 97 + 113 7800 fn+2) yn+ 97 = yn + h(− 3141 13720 fn + 11259 14560 fn+ 17 + 59136860 fn+1 − 747 5600 fn+ 97 + 11421 891800 fn+2)  (33) 2.3. Two-step, Second Derivative and One off-step Method Here k = 2, l = 2, and v = 17 is the chosen off-step point. The new approximation now becomes: u(x) = αn(x)yn + h ( β1,0(x) fn + β1, 17 (x) fn+ 17 +β1,1(x) fn+1 + β1,2(x) fn+2 ) +h2 ( β2,0(x) fn + β2, 17 (x)gn+ 17 +β2,1(x)gn+1 + β2,2(x)gn+2 ) (34) where g(xn+ j) = f ′(xn+ j) = u′′(xn+ j), j = 0, · · · , k. Setting r = 1 and s = 4, the interpolating function now be- comes: u(x) = 8∑ j=0 a j x j. (35) Obtaining the first and second derivatives, we have: u′(x) = ∑8 j=1 ja j x j−1 u′′(x) = ∑8 j=2 j( j − 1)a j x j−2. (36) Imposing the following conditions as obtained from (6): u(xn+ j) = yn+ j, j = 0 (37) u′(xn+ j) = fn+ j, j = 0, 1 7 , 1, 2. (38) u′′(xn+ j) = gn+ j, j = 0, 1 7 , 1, 2. (39) Thus, the required method is given in block method as: yn+2 = − 31h2 gn+2 845 + 16h2 gn+1 45 + 19208h2 gn+ 17 7605 + 7 5 h2 gn + 781h fn+2 2197 + 4 9 h fn+1 − 470596h fn+ 17 19773 + 25h fn + yn yn+ 17 = − 15997h2 gn+2 66805808160 − 33203h2 gn+1 3557705760 − 516419h2 gn+ 17 250417440 + 565619h2 gn 395300640 + 191117h fn+2 173695101216 + 9703h fn+1 304946208 + 3016667h fn+ 17 39862368 + 5308663h fn 79060128 + yn yn+1 = − 43h2 gn+2 81120 − 197h2 gn+1 4320 + 333739h2 gn+ 17 730080 + 101h2 gn 480 + 539h fn+2 210912 + 919h fn+1 2592 − 19176787h fn+ 17 5694624 + 385h fn 96 + yn (40) 64 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 65 3. Analysis of Methods 3.1. Order and Error Constant In this section, analysis of the methods using Taylor’s se- ries method for the derivation of order and error constant of the methods were carried out. Following Lambert (1991) and Fatunla (1991), the linear difference operator L associated with the method is defined by: L[y(x); h] = k∑ i=0 αi(x)y(x + ih) − l∑ i=1 hi(βi,0(x)y (i)(x + ih) +βi,v(x)y (i)(x + vh) + k∑ j=1 β j,k(x)y ( j)(x + jh)) (41) where y(x) is an arbitrary function, continuously differentiable on [a, b]. Recall that; while expanding y(x + ih) and its derivatives y( j)(x + ih) as Taylor’s series about x, and collecting terms gives: L[y(x); h] = c0y(x) + c1hy (1)(x) +· · ·+ cqh qy(q)(x) +· · ·(42) where cq, q = 0, 1, 2, · · · are the constant coefficients are given as follows: c0 = 1 −αn c1 = k −βi,0 −βi,v − ∑k j=1 β j,k ... cq = kq q! − vq−1 (q − 1)! βi,v − kq−1 + 1 (q − 1)! ∑k j=1 β j,k  (43) According to Henrici (1962), the main method (4) has order p if: L[y(x); h] = O(hp+1), (44) where: c0 = c1 = · · · = cp = 0, cp+1 , 0. Therefore, cp+1 is the error constant and cp+1hp+1(xn) is the principal local truncation error at point xn. 3.2. Consistency and Zero-stability 3.2.1. Consistency According to Lambert (1991), a numerical scheme is said to be consistent if the order p of the method is such that p > 1. 3.2.2. Zero Stability A numerical method is said to be zero-stable if no root of the first characteristic polynomial has modulus greater than one that is if every root with modulus one is simple (Lambert, 1991). A Linear Multi-step Method is said to be zero-stable if no root of the first characteristic polynomial defined by: ρ(r) = det[rA(0) − A(1)] (45) where A(0) is an identity matrix of dimension k, A(1) is a matrix of dimension k and satisfy |rs| ≤ 1 and every root of |rs| = 1, s = 1(1)k, has multiplicity not exceeding as h → 0. In what follows, the k-step, m-derivative block method is written as a matrix finite difference equation of the form: A(0)Yµ+1 = A (1)Yµ + h[B (0) Fµ+1 + B (1) Fµ] +h2[C(0)Gµ+1 + C (1)Gµ] + · · · + hm[Z(0)Zµ+1 + Z (1)Zµ] (46) where m is the derivative order. Yµ+1, Yµ, Fµ+1, Fµ, Gµ+1, Gµ, · · · , Zµ+1, Zµ are column matrix determined from the derived methods. A(0), A(1), B(0), B(1), C(0), C(1), · · · , Z(0), Z(1) are coefficients matrices of the derived methods. To this end, the zero stability of the derived methods are discussed below. For method (27), the matrix difference equation is: A(0)Yµ+1 = A (1)Yµ + h[B (0) Fµ+1 + B (1) Fµ] (47) where: Yµ+1 = (yn+ 17 , yn+ 27 , yn+1) T Yµ = (yn−17 , yn− 27 , yn) T Fµ+1 = ( fn+ 17 , fn+ 27 , fn+1) T Fµ = ( fn− 17 , fn− 27 , fn) T A(0) =  1 0 0 0 1 0 0 0 1  , A(1) =  0 0 1 0 0 1 0 0 1  B(0) =  595 5880 − 91 5880 1 5880 4 21 1 21 0 − 49 24 49 24 5 24  , B(1) =  0 0 3355880 0 0 121 0 0 1924   (48) Zero-stability is a concept concerned with the limiting factor of the difference equation as h → 0. Therefore, the difference equation is reduced to: A(0)Yµ+1 = A (1)Yµ (49) in which the first characteristic polynomial: ρ(r) = 0. (50) The solution of equation (50) corresponds to equation (45), thus solving (50) with parameters obtained in (48) gives the roots of r. ρ(r) = det r  1 0 0 0 1 0 0 0 1  −  0 0 1 0 0 1 0 0 1   = 0 (51) ρ(r) = det   r 0 −1 0 r −1 0 0 r − 1   = 0 (52) from (52) we have: r2(r − 1) = 0 r = 0, 0, 1. (53) 65 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 66 Table 1: Table of Orders and Error Constants of the Derived Multi-Derivative Methods Method Evaluating points Order Error Constant xn+1 4 −31 35280 = −8.79E − 04 27 xn+ 17 4 −1 246960 = −4.05E − 06 xn+ 27 4 −1 1512630 = −6.61E − 07 xn+2 5 - 2 1575 = 1.27E − 03 33 xn+ 17 5 2491 282357600 = 8.82E − 06 xn+1 5 - 53 117600 = 4.51E − 04 xn+ 97 5 − 39447 94119200 = 4.09E − 04 xn+2 8 157 38896200 = 4.04E − 06 40 xn+ 17 8 162133 1025046183571200 = 1.58E − 10 xn+1 8 307 1244678400 = 2.47E − 07 Since no root of the characteristic equation violates the condi- tion: |rs| ≤ 1, s = 1(1)m and m is the highest power of the characteristics equation ρ(r) = 0. It is concluded that method (27) is zero-stable. Table 2: Table of roots of the characteristic equations for the Derived Multi- Derivative Methods Method Absolute value of Roots, |rs| 27 0,0,1. 33 0,0,0,1. 40 0,0,1. 3.3. Convergence Stating without proof and following the fundamental the- orem of Dahliquist as reported by Henrici (1962) for a linear multi-step method. 3.3.1. Theorem 1 The necessary and sufficient conditions for a numerical method to be convergent are that it is consistent and zero-stable. Here, it is sufficient to say that the derived methods are convergent since they are all consistent and zero-stable. 3.4. Linear Stability The linear stability properties of the derived methods are determined by expressing them in a form applicable to the test problem: y′ = λy, for which y′n = λyn, y ′′ n = λ 2 yn, · · · , y (n) n = λ nyn., λ < 0 (54) to yield: Yµ+1 = M(z)Yµ, z = hλ (55) where the amplification matrix M(z) is given by: M(z) = (A(0) − zB(0) − z2C(0) − ·· ·− zlZ(0))−1(A(1) +zB(1) + z2C(1) + · · · + zlZ(1)). (56) where z = hλ. The matrix M(z) has eigenvalues ξ1,ξ2, · · · ,ξm = 0, 0, · · · ,ξm, where the dominant eigenvalue ξm is the stability function R(z) which is a rational function with real coefficients, m is the order of R(z). The stability function and plot for each of the derived methods are as given below: Figure 1: Stability Region of Method (27) Figure 2: Stability Region of Method (33) 66 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 67 Figure 3: Stability Region of Method (40) 4. Numerical Examples and Results 5. Computational Details In the implementation of the derived methods, a system of nonlinear equations must be solved in order to obtain the desired approximation. To solve these nonlinear systems, a Newton-Krylov solver, nsoli.m or a modified Newton solver, nsold was used and compared numerical results with standard and existing methods in literature. It is important to point that the numerical methods were programmed via MATLAB 9.2 version on a personal computer with the following specifica- tions: • System name- HP Pavilion x360 Convertible • Processor- Intel(R) Core(TM) i3-7100U CPU @ 2.40GHz • Installed memory (RAM)- 8.00GB • System Type- 64-bits Operating System, x64-based pro- cessor • Operating system- Windows 10 Moreso, computational experiments were done with software optimization and only the points of emphases were shown. 5.0.1. Example 1 Advection equation where v(x, t) is constant: ∂u ∂t = −v ∂u ∂x u(x, 0) = exp(− x 2 L2 ) Exact Solution:u(x, t) = exp(− (x−vt) 2 L2 ) v = 0.5, L = 10, x ∈ [0, 1], t ∈ [0, 1] (57) is considered. 5.0.2. Example 2 An advection equation where v(x, t) = 1 + x2 1 + 2xt + 2x2 + x4 ∂u ∂t = −v(x, t) ∂u ∂x u(x, 0) = exp(− x 2 L2 ) Exact Solution:u(x, t) = exp(− 1L2 (x − t 1+x2 )) L = 10, x ∈ [0, 1], t ∈ [0, 1] (58) is also considered. 5.0.3. Example 3 The inviscid Burgers equation where v(x, t) = u(x, t) ∂u ∂t = −u(x, t) ∂u ∂x u(x, 0) = x Exact Solution: u(x, t) = x t + 1 x ∈ [0, 1], t ∈ [0, 1]. (59) 5.0.4. Example 4 A nonlinear advection equation with perturbation term: � ∂u ∂t + u ∂u ∂x = f (x, t), (x, t) ∈ D = [0, 1],×[0, 1] u(0, t) = �t, u(x, 0) = x, (60) where f (x, t) = �2 +�t + x, u(x, t) = �t + x is the true solution, � = 10−4 is also considered. 6. Tables of Results and Comparison In this section, the standards methods of Upwind, Lax method, the method of Wang et al. [11] and the Method were employed successfully for solving the advection equations. The numeri- cal results of the proposed method show some superiority over other methods in-terms of accuracy and reliability. Moreover, the Method is also effective for solving the equation with sin- gular terms and other nonlinear singular perturbation problems. Wang et al. [11] used a perturbation method and reproducing kernel method in solving a class of singularly perturbed par- tial differential equation. From Table 3 above, it could be ob- served that the Method shows high computational strength as it compares favourably well with some known standard methods reported. 1 From the table above, it could also be observed that the Method shows high computational strength as it compares favourably well with some known standard methods reported. 2 The superiority of the proposed methods could be estab- lished in the table above. From the table above, it could also be 1(x, t) is the coordinate of x and t, u(x, t) is the exact solution, UM is the Upwind method, LM is the Lax method, Method represents the solution using the Method and Relative error=| u(x, t) - Method |/u(x, t). Percentage error = Relative error × 100%. CPU time is the method’s computational time. 2(x, t) is the coordinate of x and t, u(x, t) is the exact solution, UM is the Upwind method, LM is the Lax method, Method represents the solution using the Method and Relative error=| u(x, t) - Method |/u(x, t). Percentage error = Relative error × 100%. CPU time is the method’s computational time. 67 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 68 Table 3: Numerical Comparison for Example 1 for v(x, t) = 0.5, ∆x = 0.1, ∆t = 0.05, L = 10 (x, t) u(x, t) UM LM Method (27) Relative error Percentage error (%) (0.1000,0.1000) 0.9999750003 0.9999453153 0.9998539643 0.9999937500 1.87502 E -05 0.0019 (0.2000,0.2000) 0.9999000050 0.9998334259 0.9996180921 0.9999000050 0.00000 0.00 (0.3000,0.3000) 0.9997750253 0.9996699630 0.9993277560 0.9996000800 1.74985 E -04 0.017 (0.4000,0.4000) 0.9996000800 0.9994562666 0.9989757869 0.9991004049 4.99875 E -04 0.050 (0.5000,0.5000) 0.9993751953 0.9991927826 0.9985937720 0.9984012793 9.74525 E -04 0.097 (0.6000,0.6000) 0.9991004049 0.9988796763 0.9981806744 0.9975031224 1.59872 E -03 0.16 (0.7000,0.7000) 0.9987757500 0.9985170174 0.9977971450 0.9964064722 2.37218 E -03 0.24 (0.8000,0.8000) 0.9984012793 0.9981048482 0.9974581631 0.9951119854 3.29456 E -03 0.33 (0.9000,0.9000) 0.9979770489 0.9976432089 0.9972837348 0.9936204364 4.36544 E -03 0.44 (1.0000,1.0000) 0.9975031224 0.9971321487 0.9973810576 0.9919327166 5.58435 E -03 0.56 CPU time NA 0.460000s 0.430000s 0.399600s Table 4: Numerical Comparison of Example 2 for ∆x = 0.1, ∆t = 0.05, L = 10 (x, t) u(x, t) UM LM Method (33) Relative error Percentage error (%) (0.1000,0.1000) 0.9999900991 1.000150539 1.000294258 1.0005001250 5.10031 E -04 0.051 (0.2000,0.2000) 0.9999230799 1.000304064 1.000350801 0.9999000050 2.30767 E -05 0.0023 (0.3000,0.3000) 0.9997523243 1.000452525 1.000336322 0.9996000800 1.52282 E -04 0.015 (0.4000,0.4000) 0.9994484280 1.000569182 1.000136052 0.9991004049 3.48215 E -04 0.035 (0.5000,0.5000) 0.9990004998 1.001112699 0.9998403420 0.9984012793 5.99820 E -04 0.060 (0.6000,0.6000) 0.9984130253 0.9930245386 0.9993584603 0.9975031224 9.11349 E -04 0.091 (0.7000,0.7000) 0.9977006342 1.041218330 0.9987797082 0.9964064722 1.29714 E -03 0.13 (0.8000,0.8000) 0.9968829170 0.9006435880 0.9980268384 0.9951119854 1.77647 E -03 0.18 (0.9000,0.9000) 0.9959804757 1.111721077 0.9972063456 0.9936204364 2.36956 E -03 0.24 (1.0000,1.0000) 0.9950124792 0.9510926234 0.9963208653 0.9919327166 3.09520 E -03 0.31 CPU time NA 0.510000s 0.370000s 0.250000s Table 5: Numerical Comparison of Example 3 for, ∆x = 0.1, ∆t = 0.05 (x, t) u(x, t) UM LM Method (33) Relative error Percentage error (%) (0.1000,0.1000) 0.09090909091 0.0904875000 0.0904875000 0.09090909091 1.000000 E -16 0.00 (0.2000,0.2000) 0.1666666667 0.1653231744 0.1653231744 0.1666666667 1.000000 E -16 0.00 (0.3000,0.3000) 0.2307692308 0.2283156714 0.2283156715 0.2307692308 1.000000 E -16 0.00 (0.4000,0.4000) 0.2857142857 0.2821173696 0.2821192653 0.2857142857 1.000000 E -16 0.00 (0.5000,0.5000) 0.3333333333 0.3289352242 0.3286585068 0.3333333333 1.000000 E -16 0.00 (0.6000,0.6000) 0.3750000000 0.3543939265 0.3694124359 0.3750000000 1.000000 E -16 0.00 (0.7000,0.7000) 0.4117647059 0.6111136606 0.4055658801 0.4117647059 1.000000 E -16 0.00 (0.8000,0.8000) 0.4444444444 0.00861594345 0.4381714651 0.4444444444 1.000000 E -16 0.00 (0.9000,0.9000) 0.4736842105 38.46493835 0.4681540744 0.4736842105 1.000000 E -16 0.00 (1.0000,1.0000) 0.5000000000 0.1220675888 0.4964220517 0.5000000000 1.000000 E -16 0.00 CPU time NA 0.350000s 0.320000s 0.281250s observed that the Method shows high computational strength as it compares favourably well with some known standard meth- ods reported. 3 From table 6 above, it could also be observed that the Method shows some superiority in-terms of accuracy when compared with the method of Wang et al. [11] and also compares favourably 3(x, t) is the coordinate of x and t, u(x, t) is the exact solution, UM is the Upwind method, LM is the Lax method, Method represents the solution using the Method and Relative error=| u(x, t) - Method |/u(x, t). Percentage error = Relative error × 100%. CPU time is the method’s computational time. well with some known standard methods reported. 4 7. Summary and Discussion of Results This work contains a new class of multi-derivative method. The multi-derivative method was considered up to the second 4(x, t) is the coordinate of x and t, u(x, t) is the exact solution, UM is the Upwind method, LM is the Lax method, Method represents the solution using the Method and Absolute error=| u(x, t) - Method |. CPU time is the method’s computational time. 68 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 69 Table 6: Numerical Comparison for Example 4 for v(x, t) = u(x, t), ∆x = 0.01, ∆t = 0.01, � = 10−4 (x, t) u(x, t) UM LM Wang et al. [11] Method (40) Absolute error (0.0000,0.0000) 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 (0.0100,0.0100) 0.01000100000 0.01009999010 0.01009999010 0.01000000000 0.01000010000 9.000000 E -07 (0.0300,0.0300) 0.03000300000 0.03089993849 0.03089993849 0.03000000000 0.03000030000 2.700000 E -06 (0.0500,0.0500) 0.05000500000 0.05249984036 0.08002383910 0.05000000000 0.05000050000 4.500000 E -06 (0.0600,0.0600) 0.06000600000 0.06359977227 0.1648091770 0.06000000000 0.06000060000 5.400000 E -06 (0.0800,0.0800) 0.08000800000 0.08639959684 0.4698987551 0.08000000000 0.08000080000 7.200000 E -06 (0.1000,0.1000) 0.1000100000 1.009904485 0.8868283796 0.1000000000 0.1000010000 4.500000 E -06 CPU time NA 0.400000s 0.171875s NA 0.078125s Table 7: Error Comparison for Example 1 for v(x, t) = 1.0. N Adrain (2019) PC4 Scheme Adrain (2019) PC4 Scheme Proposed Method 40 6.43092 E -03 1.81277 E -04 5.58435 E -04 60 9.73768 E -04 2.78696 E -07 3.36544 E -09 80 2.71484 E -04 5.66109 E -09 3.24956 E -10 100 1.06781 E -04 2.87180 E -09 1.33524 E -10 Figure 4: Graph of Exact Solution of Example 1 derivative (k = 2). The method was developed using some selected off-step points in its formation (this is to improve its accuracy), evaluation such as collocation and interpolation of these points were carried out and a system of algebraic equa- tions was obtained for each value of k together with number of off-step points considered. These equations were solved to ob- tain the unknown parameters. The interpolating function is then evaluated at the selected collocation points to obtain a block form of this class of methods. These methods were then ana- lyzed and were found to be stable and convergent. Table 3-7 show the numerical computational results on prob- lems advection equations. Extensive comparison was done with existing and standard methods as shown in the literature. These methods were demonstrated on some Examples and the supe- riority of the derived methods over these existing methods was established. It is worthy to say that the derived methods exhibit stronger computational strength. Also, Figures 4-11, in pairs, Figure 5: Graph of Numerical Solution of Example 1 show the graphical comparison of derived methods with exact for the advection problems considered. 8. Conclusion Multi-derivative techniques are numerical tools designed for approximating singular problems. These methodologies use the evaluations of several derivatives at some prior points. Due to the characterized multi-step nature of some class of these meth- ods, additional methods of evaluation called the predictor are always required in computation of previous values. This study produces a block form of these methodologies which need no predictions and thereby taking care of this drawback. Multi- derivative techniques have a long history of improvement in integrating Ordinary Differential Equations, but to date, a few class of these techniques have been explored in the integration of Partial Differential Equations. This research work explored 69 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 70 Figure 6: Graph of Exact Solution of Example 2 Figure 7: Graph of Numerical Solution of Example 2 Figure 8: Graph of Exact Solution of Example 3 Figure 9: Graph of Numerical Solution of Example 3 Figure 10: Graph of Exact Solution of Example 4 Figure 11: Graph of Numerical Solution of Example 4 the possibility of multi-derivative methods in integrating singu- lar advection equations. From the computational results of this study, it could be observed that the derived methodologies solve satisfactorily the advection problems. 70 M. O. Ogunniran / J. Nig. Soc. Phys. Sci. 1 (2019) 62–71 71 Acknowledgment I thank the referees for the positive enlightening comments and suggestions, which have greatly helped us in making im- provements to this paper. References [1] A. Balzano, “Mosquito: An Efficient Finite Difference Scheme for Nu- merical Simulation of 2D Advection Equation”, International Journal for Numerical Methods in Fluids 32 (1999) 2. [2] B. J. Noye and M. Dehghan, ”New explicit finite difference schemes for two-dimensional diffusion subject to specification of mass”, Numerical Methods For Partial Differential Equations 15 (1999) 4. [3] Y. Han, X. Zhang, L. Liu, Y. and Wu, ”Multiple positive solutions of sin- gular nonlinear Sturm-Liouville problems with Caratheodory perturbed term”, Journal of Applied Mathematics, (2012) 160891. [4] F. Z. Geng and S. P. Qian, ”Solving singularly perturbed multipantograph delay equations based on the reproducing kernel method”, Abstract and Applied Analysis, (2014) 794716. [5] M. Stojanovic, ”Splines difference methods for a singular perturbation problem”, Applied Numerical Mathematics, 21 (1996) 3. [6] S. Bushnaq, B. Maayah, S. Momani, and A. Alsaedi, “A Reproducing Kernel Hilbert Space Method for Solving Systems of Fractional Integro differential Equations”, Abstract and Applied Analysis (2014) 103016. [7] Z. Xianyl, ”Linear Hybrid-variable Methods for Advection Equations”, Advances in Computational Mathematics, Journal of Differential Equa- tions, Conf 23 (2016) https://doi.org/10.1007/s10444-018-9647-z. [8] E. Utku, S. Murat and K. Huseyin, “Efficient Numerical Treatment of Nonlinearities in the Advection-Diffusion-Reaction Equation”, In- ternational Journal of Numerical Methods for Heat and Fluid Flow, http://doi.org/10.1108/HFF-02-2017-0198. [9] Z. Gaetamo and P. Silvia, “On the Fractional Diffusion-Advection Equa- tion for Fluids and Plasmas”, Fluid 4 (2019) 62. [10] S. Adrian, “A Generalization of Prefactored Compact Schemes for the So- lution of Adevection Equations”, arXiv:1902.04429v1 math.NA (2019). [11] Y. Wang, H. Yu, G. Tan, and S. Qu, “Solving a Class of Singularly Per- turbed Partial Differential Equations by Using the Perturbation Method and Reproducing Kernel Method”, Hindawi Abstract and Applied Anal- ysis (2014) 5. 71