Next Article in Journal
Finite Difference–Collocation Method for the Generalized Fractional Diffusion Equation
Next Article in Special Issue
Local Fractional Homotopy Perturbation Method for Solving Coupled Sine-Gordon Equations in Fractal Domain
Previous Article in Journal
Exponential Stability of Highly Nonlinear Hybrid Differently Structured Neutral Stochastic Differential Equations with Unbounded Delays
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Second Derivative Block Hybrid Methods for the Numerical Integration of Differential Systems

by
Dauda Gulibur Yakubu
1,
Ali Shokri
2,
Geoffrey Micah Kumleng
3 and
Daniela Marian
4,*
1
Department of Mathematical Sciences, Abubakar Tafawa Balewa University, Bauchi PMB 0248, Nigeria
2
Department of Mathematics, Faculty of Science, University of Maragheh, Maragheh 83111-55181, Iran
3
Department of Mathematics, University of Jos, Jos PMB 2084, Nigeria
4
Department of Mathematics, Technical University of Cluj-Napoca, 28 Memorandumului Street, 400114 Cluj-Napoca, Romania
*
Author to whom correspondence should be addressed.
Fractal Fract. 2022, 6(7), 386; https://doi.org/10.3390/fractalfract6070386
Submission received: 3 June 2022 / Revised: 27 June 2022 / Accepted: 7 July 2022 / Published: 10 July 2022
(This article belongs to the Special Issue Advances in Nonlinear Differential Equations with Applications)

Abstract

:
The second derivative block hybrid method for the continuous integration of differential systems within the interval of integration was derived. The second derivative block hybrid method maintained the stability properties of the Runge–Kutta methods suitable for solving stiff differential systems. The lack of such stability properties makes the continuous solution not reliable, especially in solving large stiff differential systems. We derive these methods by using one intermediate off-grid point in between the familiar grid points for continuous solution within the interval of integration. The new family had a high accuracy, non-overlapping piecewise continuous solution with very low error constants and converged under the suitable conditions of stability and consistency. The results of computational experiments are presented to demonstrate the efficiency and usefulness of the methods, which also indicate that the block hybrid methods are competitive with some strong stability stiff integrators.

1. Introduction

The system whose numerical approximation is sought is written in the form
{ d y d x = f ( x , y ) , ( a x b ) , y ( x 0 ) = y 0 .
In Equation (1), y : [ a , b ] R m and f : [ a , b ] × R m R m are differentiable. To obtain accurate integration methods, which combine, to some extent, the advantages of the Runge–Kutta methods (RKMs) and linear multistep methods (LMMs), the use of the multistep collocation technique has been proposed by many authors, for example, Onumanyi et al. [1], Chollom and Jackiewicz [2], Chollom and Onumanyi [3], and Jator [4]. In this work, methods were designed for finding a continuous approximate solution of the system in Equation (1), where y(x) belongs to C1([a, b], Rm) and the set of points defined as
Ω : a = x 0 < x 1 < < x n + 1 = b ,  
so that
x n : x n = x 0 + n h ,   n = 0 , 1 , , N 1 , h = x n + 1 x n   and   N = ( b a ) / h .
The h in this paper, for simplicity, is a constant and N is a positive integer. Some of the methods derived for Equation (1) were, in fact, to evaluate the solution only at the first derivative of Equation (1). Long before our consideration of introducing the off-step points, Gragg and Stetter [5], Butcher [6,7], and Gear [8] had already considered introducing some off-step points and referred to them as generalized multistep predictor–corrector methods, a modified multistep method, and hybrid methods, respectively. Similarly, the introduction of the second derivative terms had already been considered by many authors. For example, earlier, Urabe [9] worked on a second derivative method with y ( x ) = g ( x , y ) to obtain a starting method for the single-step integrator in the paper. In [10], Mitsui changed slightly Urabe’s PC pair to improve the performance of the method. To unify and extend this result after some years, Cash [11] generalized the PC pair of Urabe’s type of method. Gupta, in [12], derived and implemented second-derivative methods. Shaintani [13,14] suggested some integration algorithms very similar to RKMs, with y ( x ) = f ( x , y ) and y ( x ) = g ( x , y ) . In [15], the author constructed ( p , q ) -stage RKMs, which exhibit y ( x ) = g ( x , y ) evaluation. Chan and Tsai [16] considered explicit two-derivative RKMs, which are cheaper to calculate with fewer function evaluations than the standard RKMs. Recently, many authors have worked on methods to obtain better approximate solutions to differential equations or on stability properties to improve the accuracy and efficiency of solution of differential equations (see, for example, [17,18,19,20,21,22,23,24]).
In this article, we extend the work of Yakubu et al. [25] to derive block-hybrid methods that show a high order of accuracy with very low error constants and large regions of absolute stability and converge rapidly to the required solution. We should also point out that the effectiveness of this class of methods for the treatment of stiff systems is shown on the basis of their attractive properties and the efficient technique to deal with a large system of a stiff initial value problem of ordinary differential equations.
Definition 1
([26]). Let  Y m and  F m be vectors given by
Y m = ( y n , y n + 1 , , y n + r 1 ) T ,
F m = ( f n , f n + 1 , , f n + r 1 ) T .
Then the k-block method is of the form
Y m = i = 1 k A i Y m i + i = 0 k B i F m i .
If r = 1, then the above equation in Equation (2) is just the classical k-step method. When B0 = 0, Equation (2) is explicit; otherwise, it is implicit.
The below diagram depicts the idea of the new methods.
In Figure 1, [a, b] is divided into a series of equal lengths of a block of six points with size or length h. The approximate solutions { y n + u , y n + 1 , y n + v , y n + 2 , y n + w , y n + 3 } are computed simultaneously in the block at the points { x n + u , x n + 1 , x n + v , x n + 2 , x n + w , x n + 3 } in the kth block. Since the methods are self-starting, we do not need predictors to start the block methods.

2. The Block Hybrid Methods

The block hybrid methods in this segment are based on the polynomial of the form
y ( x ) = α 0 + α 1 x + α 2 x 2 + + α p 1 x p 1 = i = 0 p 1 α i x i
and are referred to as interpolation polynomials, which is twice continuously differentiable. The y ( x ) is interpolated at { x n + j } , and y ( x ) and y ( x ) are collocated at { c n + j } to obtain the system,
y ( x n + j ) = y ( x n + j h ) , j { 0 , 1 , 2 , , r 1 } ,
y ( c n + j ) = f n + j = f ( x n + j h , y ( x n + j h ) ) , j = 0 , 1 , 2 , , s 1 ,
y ( c n + j ) = g n + j = f x + f y y = f x + f f y , j = 0 , 1 , 2 , , t 1 .
Following Yakubu et al. [25], we put Equations (4)–(6) as:
V α = y
to have
V = ( 1 x n x n 2 x n 3 x n 4 x n p 1 1 x n + r 1 x n + r 1 2 x n + r 1 3 x n + r 1 4 x n + r 1 p 1 0 1 2 x n 3 x n 2 4 x n 3 D x n p 2 0 1 2 x n + s 1 3 x n + s 1 2 4 x n + s 1 3 D x n + s 1 p 2 0 0 2 6 x n 12 x n 2 D x n p 3 0 0 2 6 x n + t 1 12 x n + t 1 2 D x n + t 1 p 3 )
α = ( α 0 , α 1 , α 2 , , α p 1 ) T ,   y = ( y n , , y n + r 1 , f n , , f n + s 1 , g n , , g n + t 1 ) T .
The D = ( p 1 ) and D = ( p 1 ) ( p 2 ) are first and second derivatives. From Equation (7), we have
α = U y ,   where   U = V 1
which is rearranged to have
y ( x ) = j = 0 r 1 α j ( x ) y n + j + h j = 0 s 1 β j ( x ) f n + j + h 2 j = 0 t 1 γ j ( x ) g n + j ,
where
α j ( x ) = i = 0 p 1 α j , i + 1 x i , j = 0 , 1 , 2 , , r 1 ,
h β j ( x ) = h i = 0 p 1 β j , i + 1 x i , j = 0 , 1 , 2 , , s 1 ,
h 2 γ j ( x ) = h 2 i = 0 p 1 γ j , i + 1 x i , j = 0 , 1 , 2 , , t 1 .
In fact, the coefficients in Equation (11) can be calculated from the inverse of the matrix U, as in Equation (9), or written as
U = V 1 .
Insert Equation (11) into Equation (10) (see Yakubu et al. [25,27]) to have
y ( x ) = j = 0 r 1 i = 0 r + s + t 1 α i + 1 , j y n + j P i ( x ) + h j = 0 s 1 i = 0 r + s + t 1 β i + 1 , j f n + j P i ( x ) + h 2 j = 0 t 1 i = 0 r + s + t 1 γ i + 1 , j g n + j P i ( x )
which is factorized to obtain
y ( x ) = i = 0 r + s + t 1 { j = 0 r 1 α i + 1 , j y n + j + h j = 0 s 1 β i + 1 , j f n + j + h 2 j = 0 t 1 γ i + 1 , j g n + j } P i ( x ) = i = 0 r + s + t 1 φ i P i ( x )
where
ϕ i = j = 0 r 1 α i + 1 , j y n + j + h j = 0 s 1 β i + 1 , j f n + j + h 2 j = 0 t 1 γ i + 1 , j g n + j .
Then Equation (13) becomes
y ( x ) = { j = 0 r 1 α j , 1 y n + j + h j = 0 s 1 β j , 1 f n + j + h 2 j = 0 t 1 γ j , 1 g n + j ,
j = 0 r 1 α j , 2 y n + j + h j = 0 s 1 β j , 2 f n + j + h 2 j = 0 t 1 γ j , 2 g n + j ,
j = 0 r 1 α j , r + s + t 1 y n + j + h j = 0 s 1 β j , r + s + t 1 f n + j + h 2 j = 0 t 1 γ j , r + s + t 1 g n + j } ( 1 , x , x 2 , , x r + s + t 1 ) T .
Expanding Equation (14) fully, we obtain
y ( x ) = ( y n , , y n + r 1 , f n , , f n + s 1 , g n , , g n + t 1 ) T U T ( 1 , x , , x r + s + t 1 ) T .
The T in Equation (15) denotes the transpose of.

3. Specification of the Multistep Block Hybrid Methods

3.1. Block Hybrid Method of Seventh Order

In this segment, we use the multistep approach for the construction of the new block hybrid method with symmetric points of order seven. We introduce three off-step points, u = 1 2 , v = 3 2 , and w = 5 2 , and η = ( x x n ) for the construction of the continuous scheme. These points are carefully chosen to guarantee the convergence of the method, as pointed out by [28,29,30]. From Equation (10), putting r = 1 and s = 7 gives the block hybrid method of the form
y ( x ) = α 0 ( x ) y n + h j = 0 6 β j ( x ) f n + j .
Simplifying Equation (16), the interpolation and collocation polynomial in Equation (10) reduces to the proposed continuous scheme of the form in Equation (15), as follows:
y ( x ) = α 0 ( x ) y n + h [ β 0 ( x ) f n + β 1 ( x ) f n + u + β 2 ( x ) f n + 1 + β 3 ( x ) f n + v + β 4 ( x ) f n + 2 + β 5 ( x ) f n + w + β 6 ( x ) f n + 3 ]
where
α 0 ( x ) = 1 ,
β 0 ( x ) = [ 96 η 7 1176 h η 6 + 5880 h 2 η 5 15,435 h 3 η 4 + 22,736 h 4 η 3 18,522 h 5 η 2 + 7560 h 6 η 7560 h 6 ] ,
β 1 ( x ) = [ 24 η 7 + 280 h η 6 1302 h 2 η 5 + 3045 h 3 η 4 3654 h 4 η 3 + 1890 h 5 η 2 315 h 6 ] ,
β 2 ( x ) = [ 480 η 7 5320 h η 6 + 23,016 h 2 η 5 48,405 h 3 η 4 + 49,140 h 4 η 3 18,900 h 5 η 2 2520 h 6 ] ,
β 3 ( x ) = [ 240 η 7 + 2520 h η 6 10,164 h 2 η 5 + 19,530 h 3 η 4 17,780 h 4 η 3 + 6300 h 5 η 2 945 h 6 ] ,
β 4 ( x ) = [ 480 η 7 4760 h η 6 + 17,976 h 2 η 5 32,235 h 3 η 4 + 27,720 h 4 η 3 9450 h 5 η 2 2520 h 6 ] ,
β 5 ( x ) = [ 24 η 7 + 224 h η 6 798 h 2 η 5 + 1365 h 3 η 4 1134 h 4 η 3 + 378 h 5 η 2 315 h 6 ] ,
β 6 ( x ) = [ 96 η 7 840 h η 6 + 2856 h 2 η 5 4725 h 3 η 4 + 3836 h 4 η 3 1260 h 5 η 2 7560 h 6 ] ,
Evaluate Equation (17) at x n + u , x n + 1 x n + v , x n + 2 , x n + w , and x n + 3 to obtain the method:
y n + u = y n + h 120,960 [ 19,087 f n + 65,112 f n + u 46,461 f n + 1 + 37,504 f n + v 20,211 f n + 2 + 6312 f n + w 863 f n + 3 ]
y n + 1 = y n + h 7560 [ 1139 f n + 5640 f n + u + 33 f n + 1 + 1328 f n + v 807 f n + 2 + 264 f n + w 37 f n + 3 ]
y n + v = y n + h 4480 [ 685 f n + 3240 f n + u + 1161 f n + 1 + 2176 f n + v 729 f n + 2 + 216 f n + w 29 f n + 3 ]
y n + 2 = y n + h 945 [ 143 f n + 696 f n + u + 192 f n + 1 + 752 f n + v + 87 f n + 2 + 24 f n + w 4 f n + 3 ]
y n + w = y n + h 24,192 [ 3715 f n + 17,400 f n + u + 6375 f n + 1 + 16,000 f n + v + 11,625 f n + 2 + 5640 f n + w 275 f n + 3 ]
y n + 3 = y n + h 280 [ 41 f n + 216 f n + u + 27 f n + 1 + 272 f n + v + 27 f n + 2 + 216 f n + w + 41 f n + 3 ]

3.2. Second-Derivative Block Hybrid Method of Order 14

Here, we introduce the second-derivative term to have the block hybrid method of order 14 whereby we have the following interpolation and collocation polynomial of the form
y ( x ) = α 0 ( x ) y n + h j = 0 6 β j ( x ) f n + j + h 2 j = 0 6 γ j ( x ) g n + j
Simplify Equation (19) to obtain the proposed continuous scheme of the form in Equation (15) as:
y ( x ) = α 0 ( x ) y n + h [ β 0 ( x ) f n + β 1 ( x ) f n + u + β 2 ( x ) f n + 1 + β 3 ( x ) f n + v + β 4 ( x ) f n + 2 + β 5 ( x ) f n + w + β 6 ( x ) f n + 3 ] + h 2 [ γ 0 ( x ) g n + γ 1 ( x ) g n + u + γ 2 ( x ) g n + 1 + γ 3 ( x ) g n + v + γ 4 ( x ) g n + 2 + γ 5 ( x ) g n + w + γ 6 ( x ) g n + 3 ]
where
α 0 ( x ) = 1 ,
β 0 ( x ) = [ 56 η 14 10,125 h 13 16,384 η 13 131,625 h 12 + 38,339 η 12 30,375 h 11 15,428 η 11 2025 h 10 + 1,028,069 η 10 33,750 h 9 2,578,814 η 9 30,375 h 8 + 54,733,637 η 8 324,000 h 7 2,737,391 η 7 11,340 h 6 + 29,786,393 η 6 121,500 h 5 173,613,232 η 5 101,250 h 4 + 1,383,221 η 4 18,000 h 3 48,587 η 3 2700 h 2 + η ] ,
β 1 ( x ) = [ 352 η 14 3375 h 13 100,064 η 13 43,875 h 12 + 3016 η 12 135 h 11 320,048 η 11 2475 h 10 + 2,762,014 η 10 5625 h 9 12,992,474 η 9 10,125 h 8 + 1,588,042 η 8 675 h 7 14,281,196 η 7 4725 h 6 + 2,995,066 η 6 1125 h 5 2,864,446 η 5 1875 h 4 + 12,798 η 4 25 h 3 375 η 3 5 h 2 ] ,
β 2 ( x ) = [ 8 η 14 27 h 13 2192 η 13 351 h 12 + 1583 η 12 27 h 11 31984 η 11 99 h 10 + 104,287 η 10 90 h 9 229,583 η 9 81 h 8 + 4,159,523 η 8 864 h 7 2,142,395 η 7 378 h 6 + 36,235 η 6 8 h 5 35,194 η 5 15 h 4 + 2865 η 4 4 h 3 100 η 3 h 2 ] ,
β 3 ( x ) = [ 256 η 13 1053 h 12 128 η 12 27 h 11 + 36,224 η 11 891 h 10 5440 η 10 27 h 9 + 154,928 η 9 243 h 8 12,088 η 8 9 h 7 + 1,079,840 η 7 567 h 6 48,448 η 6 27 h 5 + 436,624 η 5 405 h 4 10,160 η 4 27 h 3 + 1600 η 3 27 h 2 ] ,
β 4 ( x ) = [ 8 η 14 27 h 13 + 2176 η 13 351 h 12 1559 η 12 27 h 11 + 31,244 η 11 99 h 10 101,107 η 10 90 h 9 + 221,290 η 9 81 h 8 3,999,539 η 8 864 h 7 + 4,132,745 η 7 756 h 6 158,999 η 6 36 h 5 + 70,051 η 5 30 h 4 11,745 η 4 16 h 3 + 425 η 3 4 h 2 ] ,
β 5 ( x ) = [ 352 η 14 3375 h 13 + 92,128 η 13 43,875 h 12 63,496 η 12 3375 h 11 + 244,912 η 11 2475 h 10 1,908,574 η 10 5625 h 9 + 8,060,938 η 9 10,125 h 8 4,402,922 η 8 3375 h 7 + 1,411,868 η 7 945 h 6 146,794 η 6 125 h 5 + 1,136,462 η 5 1875 h 4 23,334 η 4 125 h 3 + 664 η 3 25 h 2 ] ,
β 6 ( x ) = [ 56 η 14 10,125 h 13 + 14,192 η 13 131,625 h 12 1139 η 12 1215 h 11 + 9712 η 11 2025 h 10 180,403 η 10 11,250 h 9 + 1,117,297 η 9 30,375 h 8 3,827,857 η 8 64,800 h 7 + 1,884,301 η 7 28,350 h 6 12,502,381 η 6 243,000 h 5 + 1,326,917 η 5 50,625 h 4 21,559 η 4 2700 h 3 + 152 η 3 135 h 2 ] ,
γ 0 ( x ) = [ 8 η 14 14,175 h 12 112 η 13 8775 h 11 + 791 η 12 6075 h 10 392 η 11 495 h 9 + 21,581 η 10 6750 h 8 18,277 η 9 2025 h 7 + 1,184,153 η 8 64,800 h 6 3613 η 7 135 h 5 + 685,307 η 6 24,300 h 4 23,569 η 5 1125 h 3 + 37,849 η 4 3600 h 2 49 η 3 15 h + η 2 2 ] ,
γ 1 ( x ) = [ 32 η 14 1575 h 12 1312 η 13 2925 h 11 + 40 η 12 95 h 10 4304 η 11 165 h 9 + 12,594 η 10 125 h 8 60,514 η 9 225 h 7 + 22,796 η 8 45 h 6 212,308 η 7 315 h 5 + 46,658 η 6 75 h 4 47618 η 5 125 h 3 + 702 η 4 5 h 2 24 η 3 h ] ,
γ 2 ( x ) = [ 8 η 14 63 h 12 320 η 13 117 h 11 + 79 η 12 3 h 10 1644 η 11 11 h 9 + 16,649 η 10 30 h 8 38,182 η 9 27 h 7 + 725,969 η 8 288 h 6 791,491 η 7 252 h 5 + 21,449 η 6 8 h 4 29,929 η 5 20 h 3 + 495 η 4 h 2 75 η 3 h ] ,
γ 3 ( x ) = [ 128 η 14 567 h 12 128 η 13 27 h 11 + 10,784 η 12 243 h 10 6592 η 11 27 h 9 + 118,264 η 10 135 h 8 19,352 η 9 9 h 7 + 298,168 η 8 81 h 6 830,608 η 7 189 h 5 + 872,360 η 6 243 h 4 258,952 η 5 135 h 3 + 5480 η 4 9 h 2 800 η 3 9 h ] ,
γ 4 ( x ) = [ 8 η 14 63 h 12 304 η 13 117 h 11 + 71 η 12 3 h 10 1392 η 11 11 h 9 + 13,229 η 10 30 h 8 28,373 η 9 27 h 7 + 503,201 η 8 288 h 6 255,529 η 7 126 h 5 + 19,361 η 6 12 h 4 4208 η 5 5 h 3 + 4185 η 4 16 h 2 75 η 3 2 h ] ,
γ 5 ( x ) = [ 32 η 14 1575 h 12 1184 η 13 2925 h 11 + 808 η 12 225 h 10 3088 η 11 165 h 9 + 7954 η 10 125 h 8 33,338 η 9 225 h 7 + 54,256 η 8 225 h 6 86,468 η 7 315 h 5 + 5366 η 6 25 h 4 13,786 η 5 125 h 3 + 846 η 4 25 h 2 24 η 3 5 h ] ,
γ 6 ( x ) = [ 8 η 14 14,175 h 12 32 η 13 2925 h 11 + 23 η 12 243 h 10 716 η 11 1485 h 9 + 10,841 η 10 6750 h 8 7436 η 9 2025 h 7 + 76,213 η 8 12,960 h 6 8317 η 7 1260 h 5 + 247,819 η 6 48,600 h 4 35,009 η 5 13,500 h 3 + 71 η 4 90 h 2 η 3 9 h ] ,
Evaluating the continuous scheme in Equation (20) as usual at x n + u , x n + 1 x n + v , x n + 2 , x n + w , and x n + 3 , we obtain the method
y n + u = y n + h 1,245,404,160,000 [ 199,368,819,177 f n 68,951,829,552 f n + u 380,416,470,375 f n + 1 + 300,642,304,000 f n + v + 457,138,998,375 f n + 2 + 110,327,270,448 f n + w + 4,592,987,927 f n + 3 ] h 2 249,080,832,000 [ 1,784,098,013 g n 33,488,665,488 g n + u 71,514,207,675 g n + 1 77,935,000,000 g n + v 3,164,886,075 g n + 2 3,963,034,512 g n + w 90,441,763 g n + 3 ]
y n + 1 = y n + h 4,864,860,000 [ 783,720,817 f n + 706,775,424 f n + u 457,058,625 f n + 1 + 1,387,808,000 f n + v + 1,957,353,375 f n + 2 + 466,919,808 f n + w + 19,341,201 f n + 3 ] h 2 972,972,000 [ 7,057,013 g n 117,681,984 g n + u 337,970,925 g n + 1 336,793,600 g n + v 134,615,475 g n + 2 16,742,592 g n + w 380,629 g n + 3 ]
y n + v = y n + h 5,125,120,000 [ 826,473,395 f n + 775,497,456 f n + u + 688,759,875 f n + 1 + 2,699,264,000 f n + v + 2,168,488,125 f n + 2 + 508,254,480 f n + w + 20,942,669 f n + 3 ] h 2 1,025,024,000 [ 7,450,095 g n 123,030,576 g n + u 333,689,625 g n + 1 390,561,600 g n + v 147,584,025 g n + 2 18,189,360 g n + w 411,921 g n + 3 ]
y n + 2 = y n + h 152,026,875 [ 24,532,563 f n + 23,488,800 f n + u + 23,587,500 f n + 1 + 116,768,000 f n + v + 99,037,875 f n + 2 + 15,993,312 f n + w + 645,700 f n + 3 ] h 2 30,405,375 [ 221,317 g n 3,633,120 g n + u 9,727,200 g n + 1 10,524,800 g n + v 5,041,125 g n + 2 567,648 g n + w 12,680 g n + 3 ]
y n + w = y n + h 1,992,646,656 [ 322,126,585 f n + 322,599,120 f n + u + 379,475,625 f n + 1 + 1,617,920,000 f n + v + 1,719,564,375 f n + 2 + 609,445,680 f n + w + 10,485,255 f n + 3 ] h 2 1,992,646,656 [ 14,560,225 g n 235,515,600 g n + u 614,964,375 g n + 1 623,480,000 g n + v 210,324,375 g n + 2 64,098,000 g n + w 1,010,975 g n + 3 ]
y n + 3 = y n + h 20,020,000 [ 3,310,219 f n + 5,014,656 f n + u + 11,161,125 f n + 1 + 21,088,000 f n + v + 11,161,125 f n + 2 + 5,014,656 f n + w + 3,310,219 f n + 3 ] h 2 4,004,000 [ 30,711 g n 409,536 g n + u 726,975 g n + 1 + 726,975 g n + 2 + 409,536 g n + w 30,711 g n + 3 ]
The order and error constants for the constructed block hybrid methods are presented in Table 1. It is clear from the table that the members of the block hybrid method without a second-derivative evaluation are all of order seven except the last member in the block, which has an order higher than the remaining members in the block (order eight). The members of the block hybrid method with a second derivative are of uniform accuracy of order 14 with smaller error constants and, hence, are more accurate than those without a second derivative.

4. Regions of Absolute Stability (RAS) of the Block Hybrid Methods

Generally, in designing a new numerical method, it is very important to consider the stability properties of the method. Therefore, in this paper, we reformulate the block hybrid methods, as in [31,32], by the partitioning ( s + r ) × ( s + r ) of the form
[ Y [ n ] y [ n 1 ] ] = [ A | U B | V ] [ h f ( Y [ n ] ) y [ n ] ] ,   n = 1 ,   2 ,   ,   N ,
where
Y [ n ] = [ Y 1 [ n ] Y 2 [ n ] Y s [ n ] ] ,   y [ n 1 ] = [ y 1 [ n 1 ] y 2 [ n 1 ] y r [ n 1 ] ] ,   f ( Y [ n ] ) = [ f ( Y 1 [ n ] ) f ( Y 2 [ n ] ) f ( Y s [ n ] ) ] ,   y [ n ] = [ y 1 [ n ] y 2 [ n ] y r [ n ] ] ,
A = [ 0 0 A B ] ,   U = [ I 0 0 0 μ e μ ] ,   B = [ A B 0 0 ν T ω T ] ,   V = [ I μ e μ 0 0 I 0 0 I θ ] ,
and e = [ 1 , , 1 ] T R m .
Thus, Equation (22a) is
[ Y 1 [ n ] Y 2 [ n ] Y s [ n ] y 1 [ n ] y r [ n ] ] = [ A | U B | V ] [ h f ( Y 1 [ n ] ) h f ( Y 2 [ n ] ) h f ( Y s [ n ] ) y 1 [ n 1 ] y r [ n 1 ] ] .
The values r and s denote output and stage values, respectively. Applying Equation (22) to the linear test equation y = λ y , x 0 and λ C , we have M ( z ) as
M ( z ) = V + z B ( 1 z A ) 1 U
and the stability polynomial ρ ( η , z ) of the method can easily be obtained as
ρ ( η , z ) = d e t ( η I M ( z ) ) .
The region of absolute stability of the method is defined as
= x C : ρ ( η , z ) = 1 | η | 1 .
Computing the stability function gives the stability polynomial of the method, which is plotted to produce the required graph of the region of absolute stability of the method, as shown in Figure 2.
Remark 1.
In the stable block hybrid second-derivative implicit method, we added the matrix D1 obtained from the coefficients of h 2 to the matrices A, C, B, and D, which enabled us to plot the region of absolute stability of the new method. The region of absolute stability of method (18) is A( α )-stable while the region of absolute stability of the second-derivative implicit method (21) is A-stable since the region contains the complex plane outside the enclosed figure.

5. Numerical Illustrations

For the illustration of the performance of the derived methods, we consider both linear and nonlinear challenging systems. To provide a direct comparison, Matlab software codes were written for the preliminary test experiments using a fixed step length. We present the calculated results in tables and depict the curves in figures. Here, nfe and Ext are the function evaluations and exact values, respectively.
Example 1.
The Kaps problem [30].
We consider the nonlinear Kaps stiff system,
[ y 1 ( x ) y 2 ( x ) ] = [ 1002   y 1 ( x ) + 1000   y 2 2 ( x ) y 1 ( x ) y 2 ( x ) ( 1 + y 2 ( x ) ) ] , [ y 1 ( 0 ) y 2 ( 0 ) ] = [ 1 1 ] .
The exact value of the system is
[ y 1 ( x ) y 2 ( x ) ] = [ exp ( 2 x ) exp ( x ) ] .
The solutions of this example are shown in Table 2 and the solution curves are depicted in Figure 3.
Example 2.
Consider the linear stiff system.
[ y 1 ( x ) y 2 ( x ) y 3 ( x ) ] = [ 0 1 0 1 0 0 25 1 25 ] [ y 1 ( x ) y 2 ( x ) y 3 ( x ) ] , [ y 1 ( 0 ) y 2 ( 0 ) y 3 ( 0 ) ] = [ 0 1 2 ] .
The exact value is
[ y 1 ( x ) y 2 ( x ) y 3 ( x ) ] = [ s i n ( x ) c o s ( x ) s i n ( x ) + 2   e x p ( L x ) ] .
The results of using the newly constructed methods are shown in Table 3 and the solution curves are in Figure 4.
Example 3.
The linear problem by Enright [33] is given by:
[ y 1 ( x ) y 2 ( x ) y 3 ( x ) y 4 ( x ) ] = [ 1 0 0 0 0 10 0 0 0 0 100 0 0 0 0 1000 ] [ y 1 ( x ) y 2 ( x ) y 3 ( x ) y 4 ( x ) ] , [ y 1 ( 0 ) y 2 ( 0 ) y 3 ( 0 ) y 4 ( 0 ) ] = [ 1 1 1 1 ] .
The results of the integration are largely self-explanatory. If we examine the accuracy obtained, however, we see that the newly constructed methods are considerably accurate (see Table 4). The plotted curves are displayed in Figure 5.
Example 4.
This is given by Gear [34]:
y 1 = = 0.013   y 1 1000   y 1 y 3 , y 1 ( 0 ) = 1 , y 2 = 2500   y 2 y 3 ,     y 2 ( 0 ) = 1 , y 3 = 0.013   y 1 1000   y 1 y 3 2500   y 2 y 3 ,     y 3 ( 0 ) = 0 .
We solve this problem, and the solution curves are presented in Figure 6.
Example 5.
Here, the present problem was solved by [35]. Therefore, for comparison, we present the graphical plots of this example in Figure 7, comparing with the exact solution curves. The application of the newly derived methods to this problem is to demonstrate their performance. However, we considered only the first four components { y 1 , y 2 , y 3 , y 4 } , as shown in Table 5.
[ y 1 ( x ) y 2 ( x ) y 3 ( x ) y 4 ( x ) y 5 ( x ) y 6 ( x ) ] = [ 10 100 0 0 0 0 100 10 0 0 0 0 0 0 4 0 0 0 0 0 0 1 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0.1 ] [ y 1 ( x ) y 2 ( x ) y 3 ( x ) y 4 ( x ) y 5 ( x ) y 6 ( x ) ] , [ y 1 ( 0 ) y 2 ( 0 ) y 3 ( 0 ) y 4 ( 0 ) y 5 ( 0 ) y 6 ( 0 ) ] = [ 1 1 1 1 1 1 ]

6. Concluding Remarks

The presented second-derivative block hybrid method for a stiff system of ordinary differential equations is suitable for large systems. The second-derivative block hybrid time integrator provides good performance. Numerical results for the new second-derivative block hybrid method are promising and are demonstrably comparable to those obtained from popular high-order stiff time integrators found in the literature. Their stability properties, based on Remark 1, indicate that they are good candidates for large stiff systems. The next step of our research is to further apply some new methods to modeled differential equations that arise in other areas of scientific fields, such as chemical reaction, enzyme kinetics, cardiac electrophysiology, models of drug magnetic nanoparticle transport, and a model of tumor immune interaction, to mention just a few.

Author Contributions

Conceptualization and methodology D.G.Y.; formal analysis and software G.M.K., A.S.; investigation: D.M. and A.S.; writing—original draft preparation, D.G.Y.; writing—review and editing, D.M.; validation and visualization A.S. All authors planned the scheme, developed the mathematical modeling, and examined the theory validation. The manuscript was written through the contributions of all authors. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Tertiary Education Trust Fund (TETFund) Ref. No. TETF/ DAST&D. D/6.13/NOM-CA/BAS&BNAS. The first and third authors gratefully acknowledged the financial support of the TETFUND. The authors gratefully acknowledged the reviewers for their thorough and very fair comments and observations.

Conflicts of Interest

The authors declare that there is no conflict of interest for the study.

References

  1. Onumanyi, P.; Awoyemi, D.O.; Jator, S.N.; Sirisena, U.W. New linear multi-step methods with continuous coefficients for first order initial value problems. J. Nig. Math. Soc. 1994, 13, 37–51. [Google Scholar]
  2. Chollom, J.P.; Jackiewicz, Z. Construction of two step Runge-Kutta (TSRK) methods with large regions of absolute stability. J. Compt. Appl. Math. 2003, 157, 125–137. [Google Scholar] [CrossRef] [Green Version]
  3. Chollom, J.P.; Onumanyi, P. Variable order A-stable Adams Moulton type block hybrid methods for solution of stiff first order ODEs. J. Math. Assoc. Niger. Abacus (31) B 2004, 2, 177–192. [Google Scholar]
  4. Jator, S.N. Leaping type of algorithm for parabolic partial differential equations. JMS Natl. Math. Cent. Abuja Niger. 2013, 2, 149–172. [Google Scholar]
  5. Gragg, W.B.; Stetter, H.J. Generalized multistep predictor-corrector methods. J. Assoc. Comput. Mach. 1964, 11, 188–209. [Google Scholar] [CrossRef]
  6. Butcher, J.C. Modified multistep method for the numerical integration of ordinary differential equations. J. Assoc. Comput. Mach. 1965, 12, 124–135. [Google Scholar] [CrossRef]
  7. Butcher, J.C. A multistep generalization of Runge-Kutta method with four or five stages. J. Ass. Comput. Mach. 1967, 14, 84–99. [Google Scholar] [CrossRef]
  8. Gear, W.C. Hybrid multistep method for initial value in ordinary differential equations. J. SIAM Numer. Anal. 1965, 2, 69–86. [Google Scholar]
  9. Urabe, M. An implicit one-step method of high-order accuracy for the numerical integration of ordinary differential equations. J. Numer. Math. 1970, 15, 151–164. [Google Scholar] [CrossRef]
  10. Mitsui, T. A modified version of Urabe’s implicit single-step method. J. Comp. Appl. Math. 1987, 20, 325–332. [Google Scholar] [CrossRef] [Green Version]
  11. Cash, J.R. High order methods for the numerical integration of ordinary differential equations. Numer. Math. 1978, 30, 385–409. [Google Scholar] [CrossRef]
  12. Gupta, G.K. Implementing second-derivative multistep methods using the Nordsieck polynomial representation. J. Math. Comput. 1978, 32, 13–18. [Google Scholar] [CrossRef]
  13. Shintani, H. On one-step methods utilizing the second derivative. Hiroshima Math. J. 1971, 1, 349–372. [Google Scholar] [CrossRef]
  14. Shintani, H. On explicit one-step methods utilizing the second derivative. Hiroshima Math. J. 1972, 2, 353–368. [Google Scholar] [CrossRef]
  15. Mitsui, T. Runge-Kutta type integration formulas including the evaluation of the second-derivative, Part I. Publ. Res. Inst. Math. Sci. 1982, 18, 325–364. [Google Scholar] [CrossRef] [Green Version]
  16. Chan, R.P.K.; Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. Numer. Algorithm 2010, 53, 171–194. [Google Scholar] [CrossRef]
  17. Marian, D.; Ciplea, S.A.; Lungu, N. On the Ulam-Hyers Stability of Biharmonic Equation. Univ. Politeh. Buchar. Sci. Bull.-Ser. A-Appl. Math. Phys. 2020, 8, 141–148. [Google Scholar]
  18. Shokri, A. The Symmetric P-Stable Hybrid Obrenchkoff Methods for the numerical solution of second Order IVPS. TWMS J. Pure Appl. Math. 2012, 5, 28–35. [Google Scholar]
  19. Shokri, A. An explicit trigonometrically fitted ten-step method with phase-lag of order infinity for the numerical solution of the radial Schrödinger equation. J. Appl. Comput. Math. 2015, 14, 63–74. [Google Scholar]
  20. Shokri, A.; Saadat, H. P-stability, TF and VSDPL technique in Obrechkoff methods for the numerical solution of the Schrödinger equation. Bull. Iran. Math. Soc. 2016, 42, 687–706. [Google Scholar]
  21. Marian, D.; Ciplea, S.A.; Lungu, N. Ulam-Hyers stability of Darboux-Ionescu problem. Carpathian J. Math. 2021, 37, 211–216. [Google Scholar] [CrossRef]
  22. Sunday, J.; Shokri, A.; Marian, D. Variable step hybrid block method for the approximation of Kepler problem. Fractal Fract. 2022, 6, 343. [Google Scholar] [CrossRef]
  23. Kwami, A.M.; Kumleng, G.M.; Kolo, A.M.; Yakubu, D.G. Block hybrid multistep methods for the numerical integration of stiff systems of ordinary differential equations arising from chemical reactions. Abacus J. Math. Assoc. Nig. 2015, 42, 134–164. [Google Scholar]
  24. Singh, G.; Garg, A.; Kanwar, V.; Ramos, H. An efficient optimized adaptive step-size hybrid block method for integrating differential systems. Appl. Math. Comp. 2019, 362, 124567. [Google Scholar] [CrossRef]
  25. Yakubu, D.G.; Aminu, M.; Tumba, P.; Abdulhameed, M. An efficient family of second-derivative Runge-Kutta collocation methods for oscillatory systems. J. Nig. Math. Soc. 2018, 37, 111–138. [Google Scholar]
  26. Chu, M.T.; Hamilton, H. Parallel solution of ordinary differential equations by multi-block methods. SIAM J. Sci. Stat. Comp. 1987, 8, 342–353. [Google Scholar] [CrossRef]
  27. Yakubu, D.G.; Aminu, M.; Aminu, A. The numerical integration of stiff systems using stable multistep multi-derivative methods. Mod. Meth. Numer. Math. 2017, 8, 99–117. [Google Scholar] [CrossRef] [Green Version]
  28. Fatunla, S.O. Block methods for second order ODEs. Intern. J.Comput. Math. 1991, 41, 55–63. [Google Scholar] [CrossRef]
  29. Lambert, J.D. Computational Methods in Ordinary Differential Equations; John Willey and Sons: New York, NY, USA, 1973. [Google Scholar]
  30. Lambert, J.D. Numerical Methods for Ordinary Differential Systems; John Wiley: New York, NY, USA, 1991. [Google Scholar]
  31. Burrage, K.; Butcher, J.C. Non-linear stability for a general class of differential equation method. BIT Numer. Math. 1980, 20, 185–203. [Google Scholar] [CrossRef]
  32. Butcher, J.C. Numerical Methods for Ordinary Differential Equations, 2nd ed.; John Wiley & Sons, Ltd.: New York, NY, USA, 2008. [Google Scholar]
  33. Enright, W.H. Second derivative multistep methods for stiff ordinary differential equations. SIAM J. Numer. Anal. 1974, 11, 321–341. [Google Scholar] [CrossRef]
  34. Gear, W.C. DIFSUB for Solution of ordinary differential equations. Comm. ACM 1971, 14, 185–190. [Google Scholar] [CrossRef]
  35. Fatunla, S.O. Numerical integrators for stiff and highly oscillatory differential equations. J. Math. Comput. 1980, 34, 373–390. [Google Scholar] [CrossRef]
Figure 1. Schematic representation of the block hybrid methods for stiff system of initial value problems.
Figure 1. Schematic representation of the block hybrid methods for stiff system of initial value problems.
Fractalfract 06 00386 g001
Figure 2. Regions of absolute stability of the block hybrid methods. (a) Method (18) is A( α )-stable. (b) Method (21) is A-stable.
Figure 2. Regions of absolute stability of the block hybrid methods. (a) Method (18) is A( α )-stable. (b) Method (21) is A-stable.
Fractalfract 06 00386 g002
Figure 3. Graphical plots of Example 1 using block hybrid methods with nfe = 500. (a) Solution curve of Example 1 using (18); (b) solution curve of Example 1 using (21).
Figure 3. Graphical plots of Example 1 using block hybrid methods with nfe = 500. (a) Solution curve of Example 1 using (18); (b) solution curve of Example 1 using (21).
Fractalfract 06 00386 g003
Figure 4. Graphical plots of Example 2 using block hybrid methods with nfe = 500. (a) Solution curve of Example 2 using (18); (b) solution curve of Example 2 using (21).
Figure 4. Graphical plots of Example 2 using block hybrid methods with nfe = 500. (a) Solution curve of Example 2 using (18); (b) solution curve of Example 2 using (21).
Fractalfract 06 00386 g004
Figure 5. Graphical plots of Example 3 using block hybrid methods with nfe = 500. (a) Solution curve of Example 3 using (18); (b) solution curve of Example 3 using (21).
Figure 5. Graphical plots of Example 3 using block hybrid methods with nfe = 500. (a) Solution curve of Example 3 using (18); (b) solution curve of Example 3 using (21).
Fractalfract 06 00386 g005
Figure 6. Graphical plots of Example 4 using the block hybrid methods with nfe = 500. (a) Solution curve of Example 4 using (18); (b) solution curve of Example 4 using (21).
Figure 6. Graphical plots of Example 4 using the block hybrid methods with nfe = 500. (a) Solution curve of Example 4 using (18); (b) solution curve of Example 4 using (21).
Fractalfract 06 00386 g006
Figure 7. Graphical plots of Example 5 using block hybrid methods with nfe = 500. (a) Solution curve of Example 5 using method (18); (b) solution curve of Example 5 using method (21).
Figure 7. Graphical plots of Example 5 using block hybrid methods with nfe = 500. (a) Solution curve of Example 5 using method (18); (b) solution curve of Example 5 using method (21).
Fractalfract 06 00386 g007
Table 1. Order and error constants for the block hybrid methods.
Table 1. Order and error constants for the block hybrid methods.
MethodOrderError Constants
Block method (18)(i) yn + u, P = 7C8 = 4.4403 × 10−5
(ii) yn + 1, P = 7C8 = 3.3068 × 10−5
(iii) yn + v, P = 7C8 = 3.9236 × 10−5
(iv) yn + 2, P = 7C8 = 3.3068 × 10−5
(v) yn + w, P = 7C8 = 4.4403 × 10−5
(vi) yn + 3, P = 8C9 = 1.2555 × 10−5
Uniform order block method (21)(i) yn + u, P = 14C15 = 1.4789 × 10−12
(ii) yn + 1, P = 14C15 = 1.5718 × 10−12
(iii) yn + v, P = 14C15 = 1.5989 × 10−12
(iv) yn + 2, P = 14C15 = 1.6261 × 10−12
(v) yn + w, P = 14C15 = 1.7190 × 10−12
(vi) yn + 3, P = 14C15 = 3.1979 × 10−12
Table 2. Absolute errors in the numerical integration of example 1.
Table 2. Absolute errors in the numerical integration of example 1.
x y i Method (18)Method (21)
y 1 1.223052805026881 × 10−31.228938367083599 × 10−3
5 y 2 1.290570363021715 × 10−61.800318343625484 × 10−6
y 1 3.320709446422848 × 10−53.325679258575631 × 10−5
50 y 2 9.887815172193726 × 10−85.804723043345561 × 10−7
y 1 3.619658989642897 × 10−123.622719245691676 × 10−12
250 y 2 2.523305960607913 × 10−112.101212666995355 × 10−10
y 1 7.167561881971770 × 10−217.173620185942641 × 10−21
500 y 2 1.122741130992365 × 10−159.350493168888896 × 10−15
Table 3. Absolute errors in the numerical integration of example 2.
Table 3. Absolute errors in the numerical integration of example 2.
x y i Method (18)Method (21)
y 1 00
5 y 2 1.110223024625157 × 10−16 1.110223024625157 × 10−16
y 3 7.993605777301127 × 10−150
y 1 5.551115123125783 × 10−174.163336342344337 × 10−17
50 y 2 3.330669073875470 × 10−163.330669073875470 × 10−16
y 3 1.038058528024521 × 10−140
y 1 2.220446049250313 × 10−161.110223024625157 × 10−16
250 y 2 1.110223024625157 × 10−161.110223024625157 × 10−16
y 3 3.330669073875470 × 10−161.665334536937735 × 10−16
y 1 4.440892098500626 × 10−163.330669073875470 × 10−16
500 y 2 2.220446049250313 × 10−161.110223024625157 × 10−16
y 3 3.330669073875470 × 10−162.220446049250313 × 10−16
Table 4. Absolute errors in the numerical integration of example 3.
Table 4. Absolute errors in the numerical integration of example 3.
x y i Method (18)Method (21)
y 1 00
5 y 2 00
y 3 00
y 4 1.110223024625157 × 10−161.110223024625157 × 10−16
y 1 2.220446049250313 × 10−162.220446049250313 × 10−16
50 y 2 4.440892098500626 × 10−164.440892098500626 × 10−16
y 3 1.110223024625157 × 10−161.110223024625157 × 10−16
y 4 2.220446049250313 × 10−161.110223024625157 × 10−16
y 1 5.551115123125783 × 10−165.551115123125783 × 10−16
250 y 2 7.771561172376096 × 10−167.771561172376096 × 10−16
y 3 7.771561172376096 × 10−161.110223024625157 × 10−16
y 4 5.204170427930421 × 10−185.204170427930421 × 10−18
y 1 2.220446049250313 × 10−162.220446049250313 × 10−16
500 y 2 5.551115123125783 × 10−165.551115123125783 × 10−16
y 3 3.885780586188048 × 10−162.775557561562891 × 10−16
y 4 1.626303258728257 × 10−193.388131789017201 × 10−20
Table 5. Absolute errors of numerical integration of example 5.
Table 5. Absolute errors of numerical integration of example 5.
xyiMethod (18)Method (21)
y 1 2.024105327791403 × 10−102.220446049250313 × 10−16
y 2 4.056337835067758 × 10−101.318389841742373 × 10−16
5 y 3 00
y 4 00
y 1 1.721994824510631 × 10−93.330669073875470 × 10−16
y 2 1.453979242560521 × 10−97.771561172376096 × 10−16
50 y 3 4.440892098500626 × 10−164.440892098500626 × 10−16
y 4 01.110223024625157 × 10−16
y 1 2.077217382476237 × 10−106.591949208711867 × 10−17
y 2 1.233960850166582 × 10−111.734723475976807 × 10−18
250 y 3 2.775557561562891 × 10−178.326672684688674 × 10−17
y 4 6.661338147750939 × 10−166.661338147750939 × 10−16
y 1 2.711908290569135 × 10−126.810144895924575 × 10−19
y 2 6.195955749334348 × 10−132.710505431213761 × 10−20
500 y 3 6.938893903907228 × 10−184.857225732735060 × 10−17
y 4 3.885780586188048 × 10−163.330669073875470 × 10−16
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Yakubu, D.G.; Shokri, A.; Kumleng, G.M.; Marian, D. Second Derivative Block Hybrid Methods for the Numerical Integration of Differential Systems. Fractal Fract. 2022, 6, 386. https://doi.org/10.3390/fractalfract6070386

AMA Style

Yakubu DG, Shokri A, Kumleng GM, Marian D. Second Derivative Block Hybrid Methods for the Numerical Integration of Differential Systems. Fractal and Fractional. 2022; 6(7):386. https://doi.org/10.3390/fractalfract6070386

Chicago/Turabian Style

Yakubu, Dauda Gulibur, Ali Shokri, Geoffrey Micah Kumleng, and Daniela Marian. 2022. "Second Derivative Block Hybrid Methods for the Numerical Integration of Differential Systems" Fractal and Fractional 6, no. 7: 386. https://doi.org/10.3390/fractalfract6070386

APA Style

Yakubu, D. G., Shokri, A., Kumleng, G. M., & Marian, D. (2022). Second Derivative Block Hybrid Methods for the Numerical Integration of Differential Systems. Fractal and Fractional, 6(7), 386. https://doi.org/10.3390/fractalfract6070386

Article Metrics

Back to TopTop