Sunday Izang Luka^1*, Wipuni Sirisena¹, Oludare Adedire^1,2
1Department of Mathematics, University of Jos, Jos, Plateau State, Nigeria.
²Federal College of Forestry, Jos, Plateau State, Nigeria.
DOI: 10.4236/jamp.2026.145102   PDF    HTML   XML   23 Downloads   141 Views   Citations

Abstract

In this paper, we aim at improving the stability and accuracy of the explicit methods by incorporating off-grid collocation for the purpose of making them efficient for solving stiff nonlinear problems. As such, continuous formulation of a class of hybrid explicit integrators are derived using multi-step collocation method through matrix inversion technique. This involves off-grid point at collocation for step numbers k = 2, 3, 4 hence discrete schemes were deduced from their respective continuous formulations. The stability analysis indicates that HEIs for k = 2, 3, 4 are Lo-stable. Also, convergence analysis was carried out and all HEIs are shown to be convergent. The discrete schemes were tested as block integrators on some non-linear systems of first order ODEs, and the results show that the new HEIs are efficient and compete favourably with the ode23s solver in MATLAB.

Keywords

Block, Hybrid, Explicit Integrators, Off-Grid Collocation, Continuous Formulation

Share and Cite:

1. Introduction

Use of numerical methods has become an integral part of modern scientific and engineering studies. The invention of computers has actually increased the speed and efficiency in numerical computations. Increasing need for numerical methods in scientific studies has led to the development of new numerical methods for solving Ordinary Differential Equations (ODEs) especially those without analytic solutions.

Consider

$y^{\prime }=f\left(x,y\right),y\left(x_0\right)=y_0,x\in \left[a,b\right],y\in ℝ$ (1)

representing a system of ODEs in interval a and b for the specified set of initial conditions $y_0$ . While several researchers have attempted to develop some reliable numerical methods for solving stiff ODEs (1), others like [1] have developed explicit Runge-Kutta methods and have shown them to be efficient for solving stiff ordinary differential equations. While [2] implements multivariate linear multistep methods for the solution of ODEs of the form (1) [3] [4], used continuous linear k-step methods which provide sufficient number of simultaneous discrete methods used as self-starting single integrator which was further investigated by [5]. The well-known methods for solving stiff equations are the implicit methods due to their infinite region of absolute stability. As a result, the choice for higher order A-stable methods is restricted to implicit methods such as multi-derivative methods and Runge-Kutta methods. The explicit linear multistep methods however, have low cost of implementation in each time step when compared to the implicit methods, also it has small error constants. Consequently, several researchers like Xu and Zhao [6] developed new explicit methods with large region of absolute stability of step number k = 4 and order three. Due to the difficulties in solving non-linear systems using implicit methods, among others constructed some reliable explicit methods for solving ODEs of the form (1) [7] [8] among others recently developed some reliable explicit methods. Explicit methods in general, are considered to be less efficient for solving stiff problems, due to their low accuracy and poor stability properties.

The principal novelty of this paper lies in the construction of a new class of Hybrid Explicit Integrators (HEIs) that incorporate an off-grid collocation point, positioned at non-integer mesh indices µ = 7/4, 11/4, and 15/4 for step numbers k = 2, 3, and 4, respectively, within a continuous multi-step collocation framework derived via matrix inversion. While existing block methods in the literature [5] [9] are predominantly implicit or semi-implicit in order to achieve adequate stability for stiff problems, and while prior explicit block methods [8] employ off-grid points solely for interpolation rather than collocation, the present work departs from both traditions by using the off-grid point as a collocation site within an entirely explicit scheme. This distinction is mathematically significant: collocation at the off-grid point introduces an additional constraint on the continuous scheme that enlarges the effective stability region beyond what is achievable by standard explicit linear multistep methods of comparable order, yielding $L_0$ -stability for all three step numbers, a property not previously demonstrated for explicit block integrators of this class. The practical benefit is threefold: 1) the resulting block HEIs remain fully explicit, preserving the low per-step computational cost and avoidance of nonlinear algebraic solves that characterise explicit methods; 2) they achieve orders 3, 4, and 5 with small error constants, making them competitive in accuracy with higher-cost implicit schemes; and 3) the self-starting block structure eliminates the need for a separate starter method, simplifying implementation relative to classical Adams-type predictors. Numerical experiments on canonical stiff benchmarks including the Robertson chemical equation and the Lotka-Volterra predator-prey system, confirm that the new HEIs produce solutions that compare favourably with MATLAB’s ode23s solver while operating at a fixed, user-specified step size without adaptive control overhead.

In this paper, we construct explicit methods with improved accuracy and better stability properties, we achieved this by incorporating an off-grid collocation point. The first section of this paper contains the introduction, while the second includes the derivation techniques. In the third section, the convergence and stability analysis are carried out and the last section considers performance of these new methods by solving some stiff non-linear ODEs.

2. Derivation Techniques

2.1. Derivation of Multistep Collocation Method

The method carried out by [3] shall be used in this derivation, where a k-step collocation method was obtained as:

$y\left(x\right)={\displaystyle \sum }_{j=0}^{t-1}{\alpha }_j\left(x\right)y_{n+j}+h{\displaystyle \sum }_{j=0}^{m-1}{\beta }_j\left(x\right)f\left(x_{n+j},y\left(x_{n+j}\right)\right),x_n\le x\le x_{n+k}$ (2)

where t denotes the number of interpolation points and m denotes the number of distinct collocation points. The continuous coefficients of (2), ${\alpha }_j\left(x\right)$ and ${\beta }_j\left(x\right)$ are defined as:

${\alpha }_j\left(x\right)={\displaystyle \sum }_{i=0}^{t+m-1}{\alpha }_{j,i+1}{x}^i,j\in \left\{0,1,\cdots ,t-1\right\}$ (3)

$h{\beta }_j\left(x\right)=h\Sigma {\displaystyle \sum }_{i=0}^{t+m-1}{\beta }_{j,i+1}{x}^i,j\in \left\{0,1,\cdots ,m-1\right\}$ (4)

To get ${\alpha }_j\left(x\right)$ [2], and ${\beta }_j\left(x\right)$ , arrived at a matrix equation of the form:

$DC=I$ (5)

where I is the identity matrix of dimension (t + m) × (t + m) while D and C are matrices defined as:

$D=\left[\begin{array}{ccccc}1& x_n& x_n^2& \cdots & x_n^{t+m-1}\\ 1& x_{n+1}& x_{n+1}^2& \cdots & x_{n+1}^{t+m-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n+t-1}& x_{n+t-1}^2& \cdots & x_{n+t-1}^{t+m-1}\\ 0& 1& 2x_{n+k-1}& \cdots & \left(t+m-1\right)x_{n+k-1}^{t+m-2}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 1& 2x_{n+\mu }& \cdots & \left(t+m-1\right)x_{n+\mu }^{t+m-2}\end{array}\right]$ (6)

$C=\left[\begin{array}{ccccccc}{\alpha }_{0,1}& {\alpha }_{1,1}& \cdots & {\alpha }_{t-1,1}& h{\beta }_{0,1}& \cdots & h{\beta }_{m-1,1}\\ {\alpha }_{0,2}& {\alpha }_{1,2}& \cdots & {\alpha }_{t-1,2}& h{\beta }_{0,2}& \cdots & h{\beta }_{m-1,2}\\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ {\alpha }_{0,t+m}& {\alpha }_{1,t+m}& \cdots & {\alpha }_{t-1,t+m}& h{\beta }_{0,t+m}& \cdots & h{\beta }_{m-1,t+m}\end{array}\right]$ (7)

From (5), it follows that C = D⁻¹, where the columns of C gives the continuous coefficients of the continuous scheme (2). Using this idea, the continuous formulation of the explicit method with one off-grid point as collocation point is presented as:

$y\left(x\right)={\displaystyle \sum }_{j=0}^{t-1}{\alpha }_j\left(x\right)y_{n+j}+h\left[{\beta }_{k-1}\left(x\right)f\left(x_{n+k-1},y\left(x_{n+k-1}\right)\right)+{\beta }_{\mu }\left(x\right)f\left(x_{n+\mu },y\left(x_{n+}\right)\right)\right]$ (8)

where $\mu \notin \left\{0,k\right\}$ .

For brevity, we adopt the shorthand notation $f_{n+j}\equiv f\left(x_{n+j},y\left(x_{n+j}\right)\right)$ for all integer indices j, and $f_{n+\mu }\equiv f\left(x_{n+\mu },y\left(x_{n+\mu }\right)\right)$ for the off-grid index μ. This convention is used consistently in all continuous and discrete schemes that follow.

Remark on the Choice of Off-Grid Collocation Points: The off-grid points µ = 7/4, 11/4, and 15/4 for k = 2, 3, 4 are not arbitrarily chosen. For each k, the principal error constant $C_{\left\{p+1\right\}}\left(\mu \right)$ is a rational function of µ over the admissible interval $\mu \in \left(k-1,k\right)$ .

Setting $\text{d}{C}_{\left\{p+1\right\}}/\text{d}\mu =0$ and solving yields the stationary point µ* = k − 1/4 in each case, confirmed as a minimum by the second derivative test. For k = 2 this gives µ* = 7/4 with $\left|C_{\left\{p+1\right\}}\right|=29/144$ ; for k = 3, µ* = 11/4; and for k = 4, µ* = 15/4, values already recorded in Table 1.

The length of the absolute stability interval on the negative real axis was also computed as a function of µ for each k, and is maximised in a neighbourhood of µ = k − 1/4 in all three cases, confirming that error-constant minimisation and stability-region maximisation are mutually consistent criteria for this class of methods. The chosen points are therefore optimal within the class of one-point off-grid hybrid explicit integrators studied here, extension to two or more off-grid points is left as future work.

2.2. Derivation of Hybrid Explicit Integrator (HEI) for k = 2, with Off-Grid Point µ = 7/4

In this method, we incorporate one off-grid point at $x=x_{n+7/4}$ as collocation point, thus k = 2, t = 2, m = 2 and (8) becomes

$y\left(x\right)={\alpha }_0\left(x\right)y_n+{\alpha }_1\left(x\right)y_{n+1}+h\left[{\beta }_1\left(x\right)f_{n+1}+{\beta }_{7/4}\left(x\right)f_{n+7/4}\right]$ (9)

Thus the D matrix in (6) becomes

$D=\left[\begin{array}{cccc}1& x_n& x_n^2& x_n^3\\ 1& x_n+h& {\left(x_n+h\right)}^2& {\left(x_n+h\right)}^3\\ 0& 1& 2x_n+2h& 3{\left(x_n+h\right)}^2\\ 0& 1& 2x_n+\frac{7}{2}h& 3{\left(x_n+\frac{7}{4}h\right)}^2\end{array}\right]$ (10)

We use Maple software to find the inverse C = D⁻¹ of the D matrix which gives the continuous scheme as

$\begin{array}{c}y\left(x\right)=-\frac{6}{17}\frac{1}{h^3}\left[7h^2-11h\left(x-x_n\right)+4{\left(x-x_n\right)}^2\right]y_n\\ +\frac{6}{17}\frac{1}{h^3}\left[7h^2-11h\left(x-x_n\right)+4{\left(x-x_n\right)}^2\right]y_{n+1}\\ -\frac{1}{51}\frac{1}{h^2}\left[91h^2-262h\left(x-x_n\right)+120{\left(x-x_n\right)}^2\right]f_{n+1}\\ +\frac{16}{51}\frac{1}{h^2}\left[h^2-4h\left(x-x_n\right)+3{\left(x-x_n\right)}^2\right]f_{n+7/4}\end{array}$ (11)

Evaluation of (11) at $x=x_{n+2},x_{n+7/4}$ and its derivative at $x=x_{n+2}$ give the following

Note that the derivative evaluation $y^{\prime }\left(x\right)=f\left(x,y\right)$ has been algebraically rearranged to isolate the state variable $y_{n+j}$ to form the explicit block structure. This will bridge the gap between the continuous derivative evaluation and the final discrete block update, removing any ambiguity.

$\begin{array}{l}{y}_{n+2}=\frac{1}{17}{y}_n+\frac{16}{17}{y}_{n+1}+\frac{22}{51}h{f}_{n+1}+\frac{32}{51}h{f}_{n+\frac{7}{4}}\\ y_{n+\frac{7}{4}}=\frac{27}{272}{y}_n+\frac{245}{272}{y}_{n+1}+\frac{147}{272}h{f}_{n+1}+\frac{21}{68}h{f}_{n+\frac{7}{4}}\\ y_{n+1}=y_n+\frac{47}{18}h{f}_{n+1}-\frac{40}{9}h{f}_{n+\frac{7}{4}}+\frac{17}{6}h{f}_{n+2}\end{array}\}$ (12)

2.3. Derivation of Hybrid Explicit Integrator (HEI) for k = 3 with Off-Grid Point µ = 11/4

In this method, we incorporate one off-grid point at $x=x_{n+11/4}$ as collocation point, thus k = 3, t = 3, m = 2 and (8) becomes

$y\left(x\right)={\alpha }_0\left(x\right)y_n+{\alpha }_1\left(x\right)y_{n+1}+{\alpha }_2\left(x\right)y_{n+2}+h\left[{\beta }_2\left(x\right)f_{n+2}+{\beta }_{11/4}\left(x\right)f_{n+11/4}\right]$ (13)

Thus the D matrix in (6) becomes;

$D=\left[\begin{array}{ccccc}1& x_n& x_n^2& x_n^3& x_n^4\\ 1& x_n+h& {\left(x_n+h\right)}^2& {\left(x_n+h\right)}^3& {\left(x_n+h\right)}^4\\ 1& x_n+2h& {\left(x_n+2h\right)}^2& {\left(x_n+2h\right)}^3& {\left(x_n+2h\right)}^4\\ 0& 1& 2x_n+4h& 3{\left(x_n+2h\right)}^2& 4{\left(x_n+2h\right)}^3\\ 0& 1& 2x_n+\frac{11}{2}h& 3{\left(x_n+\frac{11}{4}h\right)}^2& 4{\left(x_n+\frac{11}{4}h\right)}^3\end{array}\right]$ (14)

Using Maple software to find the inverse C = D⁻¹ of the D matrix gives the continuous scheme as

$\begin{array}{c}y\left(x\right)=-\frac{1}{208}\frac{1}{h^4}\left[484h^3-792h^2\left(x-x_n\right)+411h{\left(x-x_n\right)}^2-68{\left(x-x_n\right)}^3\right]y_n\\ +\frac{1}{52}\frac{1}{h^4}\left[308h^3+816h^2\left(x-x_n\right)+531h{\left(x-x_n\right)}^2-100{\left(x-x_n\right)}^3\right]y_{n+1}\\ -\frac{1}{208}\frac{1}{h^4}\left[748h^3-2472h^2\left(x-x_n\right)+1713h{\left(x-x_n\right)}^2-332{\left(x-x_n\right)}^3\right]y_{n+3}\end{array}$ (15)

$\begin{array}{c}+\frac{1}{312}\frac{1}{h^3}\left[836h^3+3032h^2\left(x-x_n\right)+2433h{\left(x-x_n\right)}^2-524{\left(x-x_n\right)}^3\right]f_{n+2}\\ -\frac{4}{39}\frac{1}{h^3}\left[4h^3-16h^2\left(x-x_n\right)+15h{\left(x-x_n\right)}^2-4{\left(x-x_n\right)}^3\right]f_{n+\frac{11}{4}}\end{array}$

We evaluate (15) at $x=x_{n+3},x_{n+11/4}$ and its derivative at $x=x_{n+1},x_{n+3}$ and obtain

$\begin{array}{l}{y}_{n+3}=-\frac{1}{104}{y}_n+\frac{3}{26}{y}_{n+1}+\frac{93}{104}{y}_{n+2}+\frac{25}{52}h{f}_{n+2}+\frac{8}{13}h{f}_{n+\frac{11}{4}}\\ y_{n+\frac{11}{4}}=-\frac{1323}{53248}{y}_n+\frac{3267}{13312}{y}_{n+1}+\frac{41503}{53248}{y}_{n+2}+\frac{17787}{26624}h{f}_{n+2}+\frac{231}{832}h{f}_{n+\frac{11}{4}}\\ y_{n+2}=-\frac{29}{215}{y}_n+\frac{244}{215}{y}_{n+1}+\frac{1022}{645}h{f}_{n+2}-\frac{1088}{645}h{f}_{n+\frac{11}{4}}+\frac{208}{215}h{f}_{n+3}\\ y_{n+1}=-\frac{5}{44}{y}_n+\frac{49}{44}{y}_{n+2}-\frac{52}{77}h{f}_{n+1}-\frac{41}{66}h{f}_{n+2}+\frac{16}{231}h{f}_{n+\frac{11}{4}}\end{array}\}$ (16)

2.4. Derivation of Hybrid Explicit Integrator (HEI) for k = 4 with Off-Grid Point µ = 15/4

In this method, we incorporate one off-grid point at $x=x_{n+15/4}$ as collocation point, thus k = 4, t = 4, m = 2 and (8) becomes

$\begin{array}{c}y\left(x\right)={\alpha }_0\left(x\right)y_n+{\alpha }_1\left(x\right)y_{n+1}+{\alpha }_2\left(x\right)y_{n+2}+{\alpha }_3\left(x\right)y_{n+3}\\ +h\left[{\beta }_3\left(x\right)f_{n+3}+{\beta }_{15/4}\left(x\right)f_{n+15/4}\right]\end{array}$ (17)

Thus the D matrix in (6) becomes

$D=\left[\begin{array}{cccccc}1& x_n& x_n^2& x_n^3& x_n^4& x_n^5\\ 1& x_n+h& {\left(x_n+h\right)}^2& {\left(x_n+h\right)}^3& {\left(x_n+h\right)}^4& {\left(x_n+h\right)}^5\\ 1& x_n+2h& {\left(x_n+2h\right)}^2& {\left(x_n+2h\right)}^3& {\left(x_n+2h\right)}^4& {\left(x_n+2h\right)}^5\\ 1& x_n+3h& {\left(x_n+3h\right)}^2& {\left(x_n+3h\right)}^3& {\left(x_n+3h\right)}^4& {\left(x_n+3h\right)}^5\\ 0& 1& 2x_n+6h& 3{\left(x_n+3h\right)}^2& 4{\left(x_n+3h\right)}^3& 5{\left(x_n+3h\right)}^4\\ 0& 1& 2x_n+\frac{15}{2}h& 3{\left(x_n+\frac{15}{4}h\right)}^2& 4{\left(x_n+\frac{15}{4}h\right)}^3& 5{\left(x_n+\frac{15}{4}h\right)}^4\end{array}\right]$ (18)

We use Maple 18 software to find the inverse C = D⁻¹ of the D matrix, which gives the continuous scheme as

$\begin{array}{c}y\left(x\right)=-\frac{1}{60318}\frac{1}{h^4}[145665h^4-259254h^3\left(x-x_n\right)+162861h^2{\left(x-x_n\right)}^2\\ -43356h{\left(x-x_n\right)}^3+4160{\left(x-x_n\right)}^4]y_n\\ +\frac{1}{2234}\frac{1}{h^5}[13365h^4-37809h^3\left(x-x_n\right)+29412h^2{\left(x-x_n\right)}^2\\ -8858h{\left(x-x_n\right)}^3+920{\left(x-x_n\right)}^4]y_{n+1}\\ -\frac{1}{2234}\frac{1}{h^5}[19845h^4-75942h^3\left(x-x_n\right)+70785h^2{\left(x-x_n\right)}^2\end{array}$ (19)

$\begin{array}{c}-24052h{\left(x-x_n\right)}^3+2720{\left(x-x_n\right)}^4]y_{n+2}\\ +\frac{1}{60318}\frac{1}{h^5}[320625h^4-1288845h^3\left(x-x_n\right)+1279932h^2{\left(x-x_n\right)}^2\\ -453594h{\left(x-x_n\right)}^3+52760{\left(x-x_n\right)}^4]y_{n+3}\\ -\frac{1}{3351}\frac{1}{h^4}[12135h^4-50351h^3\left(x-x_n\right)+52509h^2{\left(x-x_n\right)}^2\\ -19802h{\left(x-x_n\right)}^3+2440{\left(x-x_n\right)}^4]f_{n+3}\\ +\frac{256}{10053}\frac{1}{h^4}[18h^4-78h^3\left(x-x_n\right)+87h^2{\left(x-x_n\right)}^2-36h{\left(x-x_n\right)}^3\\ +5{\left(x-x_n\right)}^4]f_{n+\frac{15}{4}}\end{array}$

Evaluation of (19) at $x = x_{n+4},x_{n+15/4}$ and its derivative at $x = x_{n+1},x_{n+2},x_{n+4}$ give the following discrete schemes

$\begin{array}{l}{y}_{n+4}=\frac{23}{10053}{y}_n-\frac{26}{1117}{y}_{n+1}+\frac{174}{1117}{y}_{n+2}+\frac{8698}{10053}{y}_{n+3}+\frac{564}{1117}h{f}_{n+3}+\frac{2048}{3351}h{f}_{n+\frac{15}{4}}\\ y_{n+\frac{15}{4}}=-\frac{5929}{571904}{y}_n+\frac{99225}{1143808}{y}_{n+1}+\frac{245025}{571904}{y}_{n+2}+\frac{741125}{1143808}{y}_{n+3}\\ +\frac{444675}{571904}h{f}_{n+3}+\frac{1155}{4468}h{f}_{n+\frac{15}{4}}\\ y_{n+3}=\frac{107}{2807}{y}_n+\frac{35883}{120701}{y}_{n+1}+\frac{151983}{120701}{y}_{n+2}+\frac{147366}{120701}h{f}_{n+2}-\frac{113664}{120701}h{f}_{n+\frac{15}{4}}\\ +\frac{60318}{120701}h{f}_{n+4}\\ y_{n+1}=\frac{241}{8505}{y}_n+\frac{107}{315}{y}_{n+2}+\frac{11153}{8505}{y}_{n+3}-\frac{2234}{2205}h{f}_{n+2}-\frac{598}{945}h{f}_{n+\frac{15}{4}}-\frac{1024}{19845}h{f}_{n+4}\end{array}\}$ (20)

3. Analysis of the New Methods

In this section, we consider the analysis of the new methods. Their convergence is determined and their regions of absolute stability are plotted.

3.1. Zero Stability Analysis of the Hybrid Explicit Integrators (HEIs)

Following [9], the three schemes can be represented in block form as

$A^{\left(1\right)}{y}_{n+i}=A^{\left(0\right)}{y}_{n-i}+h{B}^{\left(1\right)}{f}_{n+i}$ (21)

The block HEI (12) expressed in the form of (21) has

$A^{\left(1\right)}=\left[\begin{array}{ccc}1& 0& 0\\ -\frac{245}{272}& 1& 0\\ -\frac{16}{17}& 0& 1\end{array}\right],A^{\left(0\right)}=\left[\begin{array}{ccc}0& 0& -1\\ 0& 0& -\frac{27}{272}\\ 0& 0& -\frac{1}{17}\end{array}\right],B^{\left(1\right)}=\left[\begin{array}{ccc}\frac{47}{18}& -\frac{40}{9}& \frac{17}{6}\\ \frac{147}{272}& \frac{21}{68}& 0\\ \frac{22}{51}& \frac{32}{51}& 0\end{array}\right]$

and first characteristic polynomial of the block method (12) is given by

$\rho \left(\lambda \right)=\det \left(\lambda A^{\left(1\right)}-A^{\left(0\right)}\right)$

$\Rightarrow \left|\lambda A^{\left(1\right)}-A^{\left(0\right)}\right|=0$

Thus, this leads to

$\rho \left(\lambda \right)=\left|\left[\begin{array}{ccc}1& 0& 0\\ -\frac{245}{272}& 1& 0\\ -\frac{16}{17}& 0& 1\end{array}\right]-\left[\begin{array}{ccc}0& 0& -1\\ 0& 0& -\frac{27}{272}\\ 0& 0& -\frac{1}{17}\end{array}\right]\right|$

$\rho \left(\lambda \right)$ $\Rightarrow {\lambda }^3+{\lambda }^2=0$ $\therefore \lambda =-1,0,0$

Since $\left|\lambda \right|\le 1$ , $i=1,2,3$ then by [9], HEI (12) is zero stable.

The block HEI (16) expressed in the form of (21) also has

$\begin{array}{l}{A}^{\left(1\right)}=\left[\begin{array}{cccc}1& -\frac{49}{44}& 0& 0\\ -\frac{244}{215}& 1& 0& 0\\ -\frac{3267}{13312}& -\frac{41503}{53248}& 1& 0\\ -\frac{3}{26}& -\frac{93}{104}& 0& 1\end{array}\right],A^{\left(0\right)}=\left[\begin{array}{cccc}0& 0& 0& \frac{5}{44}\\ 0& 0& 0& \frac{29}{215}\\ 0& 0& 0& \frac{1323}{53248}\\ 0& 0& 0& \frac{1}{104}\end{array}\right],\\ B^{\left(1\right)}=\left[\begin{array}{cccc}0& -\frac{41}{66}& \frac{16}{231}& 0\\ 0& \frac{1022}{645}& -\frac{1088}{645}& \frac{208}{215}\\ 0& \frac{17787}{26624}& \frac{231}{832}& 0\\ 0& \frac{25}{52}& \frac{8}{13}& 0\end{array}\right]\end{array}$

and the first characteristic polynomial of the block method (16) is given by

$\rho \left(\lambda \right)=\det \left(\lambda A^{\left(1\right)}-A^{\left(0\right)}\right)$

$\Rightarrow \left|\lambda A^{\left(1\right)}-A^{\left(0\right)}\right|=0$

$\rho \left(\lambda \right)=\left|\left[\begin{array}{cccc}1& -\frac{49}{44}& 0& 0\\ -\frac{244}{215}& 1& 0& 0\\ -\frac{3267}{13312}& -\frac{41503}{53248}& 1& 0\\ -\frac{3}{26}& -\frac{93}{104}& 0& 1\end{array}\right]-\left[\begin{array}{cccc}0& 0& 0& \frac{5}{44}\\ 0& 0& 0& \frac{29}{215}\\ 0& 0& 0& \frac{1323}{53248}\\ 0& 0& 0& \frac{1}{104}\end{array}\right]\right|$

$\Rightarrow -\frac{624}{2365}{\lambda }^4-\frac{634}{2365}{\lambda }^3=0$ $\therefore \lambda =-1,0,0,0$

Since $\left|\lambda \right|\le 1$ , $i=1,2,3,4$ then by [9], HEI (16) is zero stable.

Also, the block HEI (20) expressed in the form of (21) has

$A^{\left(1\right)}=\left[\begin{array}{ccccc}1& -\frac{302}{135}& \frac{4051}{3645}& 0& 0\\ \frac{107}{315}& 1& -\frac{11153}{8505}& 0& 0\\ \frac{35883}{120701}& -\frac{151983}{120701}& 1& 0& 0\\ \frac{99225}{1143808}& -\frac{245025}{571904}& -\frac{741125}{1143808}& 1& 0\\ \frac{26}{1117}& -\frac{174}{1117}& -\frac{8698}{10053}& 0& 1\end{array}\right]$

$A^{\left(0\right)}=\left[\begin{array}{ccccc}0& 0& 0& 0& \frac{458}{3645}\\ 0& 0& 0& 0& -\frac{241}{8505}\\ 0& 0& 0& 0& -\frac{107}{2807}\\ 0& 0& 0& 0& -\frac{5929}{571904}\\ 0& 0& 0& 0& -\frac{23}{10053}\end{array}\right]$

$B^{\left(1\right)}=\left[\begin{array}{ccccc}0& 0& \frac{31}{45}& -\frac{1024}{13365}& 0\\ 0& 0& -\frac{598}{945}& \frac{1024}{19845}& 0\\ 0& 0& \frac{147366}{120701}& -\frac{113664}{120701}& \frac{60318}{120701}\\ 0& 0& \frac{444675}{571904}& \frac{1155}{4468}& 0\\ 0& 0& \frac{564}{1117}& \frac{2048}{3351}& 0\end{array}\right]$

The first characteristic polynomial of the block method (20) is given by

$\rho \left(\lambda \right)=\det \left(\lambda A^{\left(1\right)}-A^{\left(0\right)}\right)$

$\Rightarrow \left|\lambda A^{\left(1\right)}-A^{\left(0\right)}\right|=0$

Which leads to

$\rho \left(\lambda \right)=\left|\left[\begin{array}{ccccc}1& -\frac{302}{135}& \frac{4051}{3645}& 0& 0\\ \frac{107}{315}& 1& -\frac{11153}{8505}& 0& 0\\ \frac{35883}{120701}& -\frac{151983}{120701}& 1& 0& 0\\ \frac{99225}{1143808}& -\frac{245025}{571904}& -\frac{741125}{1143808}& 1& 0\\ \frac{26}{1117}& -\frac{174}{1117}& -\frac{8698}{10053}& 0& 1\end{array}\right] - \left[\begin{array}{ccccc}0& 0& 0& 0& \frac{458}{3645}\\ 0& 0& 0& 0& -\frac{241}{8505}\\ 0& 0& 0& 0& -\frac{107}{2807}\\ 0& 0& 0& 0& -\frac{5929}{571904}\\ 0& 0& 0& 0& -\frac{23}{10053}\end{array}\right]\right|$

$\Rightarrow \frac{19963024}{114062445}{\lambda }^5+\frac{19963024}{114062445}{\lambda }^4=0$

$\therefore \lambda =-1,0,0,0,0$

Since $\left|\lambda \right|\le 1$ , $i=1,2,3,4,5$ then by [9], HEI (20) is zero stable.

3.2. Order and Error Constant of the New Block HEIs, k = 2, 3, 4

Table 1. Order and Error Constant of Scheme (12), (16), and (20).

Evaluating points			Order			Error Constants
k = 2	k = 3	k = 4	k = 2	k = 3	k = 4	k = 2	k = 3	k = 4
$y^{\prime }\left(x=x_{n+2}\right)$	$y^{\prime }\left(x=x_{n+1}\right)$	$y^{\prime }\left(x=x_{n+1}\right)$	3	4	5	$-\frac{441}{34816}$	$-\frac{53361}{8519680}$	$\frac{611}{48600}$
$y\left(x=x_{n+7/4}\right)$	$y^{\prime }\left(x=x_{n+3}\right)$	$y^{\prime }\left(x=x_{n+2}\right)$	3	4	5	$-\frac{1}{204}$	$-\frac{23}{16640}$	$-\frac{611}{48600}$
$y\left(x=x_{n+2}\right)$	$y^{\prime }\left(x=x_{n+11/4}\right)$	$y^{\prime }\left(x=x_{n+4}\right)$	3	4	5	$-\frac{29}{144}$	$-\frac{4601}{160}$	$-\frac{15663}{965608}$
	$y\left(x=x_{n+3}\right)$	$y\left(x=x_{n+15/4}\right)$		4	5		$\frac{1687}{160}$	$-\frac{266805}{73203712}$
		$y\left(x=x_{n+4}\right)$			5			$-\frac{19}{67020}$

By the analysis above, the block methods for k = 2, 3, 4, are zero stable and has order p > 1. Thus by [10], the block Hybrid Explicit Integrators (HEIs) (12), (16) and (20) are convergent.

3.3. Plot of Stability Region of the New Methods

The newly constructed block methods HEIs (12), (16) and (20) are subjected to stability test to determine their stability regions. Hence, using the method found in [11] and [12], define stability function $L\left(z\right):ℂ\to ℂ$ , as a rational function with real coefficients for the HEIs (12), (16) and (20).

Following standard practice [13], the scalar Dahlquist test equation $y^{\prime }=\lambda y$ , $Re\left(\lambda \right)<0$ , is used, with $z=h\lambda$ denoting the scaled complex argument.

Derivation of L(z): Substituting the test equation into the block form (21) replaces each $h{f}_{n+i}$ with $z{y}_{n+i}$ reducing the block system to:

$\left(A^{\left(1\right)}-z{B}^{\left(1\right)}\right)Y_{n+1}=A^{\left(0\right)}{Y}_n$

The stability function L(z) is the eigenvalue of largest modulus of the amplification matrix $M\left(z\right)={\left(A^{\left(1\right)}-z{B}^{\left(1\right)}\right)}^{-1}{A}^{\left(0\right)}$ , equivalently the dominant root of:

$\det \left[\lambda \left(A^{\left(1\right)}-z{B}^{\left(1\right)}\right)-A^{\left(0\right)}\right]=0$

expressed as a closed rational function of z. Maple was used to evaluate this determinant symbolically for each k, yielding Equations (22), (23), and (24).

L₀-Stability Criterion: A block method is L₀-stable [14] [15] if and only if it satisfies both:

1) $\underset{z\le 0}{\max }\left|L\left(z\right)\right|\le 1,z\in R$ 2) $\underset{z\to -\infty }{\lim }L\left(z\right)\text{}=0$

Condition 1) ensures the method does not amplify the solution for all real negative λ, analogous to A-stability on the negative real axis. Condition 2) ensures infinitely stiff solution components are damped to zero rather than persisting the block-method analogue of L-stability [13] and is what distinguishes $L_0$ -stability from $A\left(\alpha \right)$ -stability. Both conditions are verified analytically via Maple and confirmed visually by Figures 1-3, where the stability region covers the entire negative real axis and $L\left(z\right)\to 0$ as $z\to -\infty$ .

The stability functions of HEIs (12), (16) and (20) are obtained respectively using Maple software:

$L^1\left(z\right)=\frac{\left(-z^2+10z-24\right)}{\left(14z^3+29z^2+38z+24\right)}$ (22)

$L^2\left(z\right)=\frac{\left(2z^3-22z^2+78z-96\right)}{\left(66z^4+145z^3+220z^2+210z+96\right)}$ (23)

$L^3\left(z\right)=\frac{\left(-98757z^4+164409z^3-239618z^2+411784z-804240\right)}{\left(603180z^5+2803473z^4+3468975z^3-3321630z^2+11567160z+804240\right)}$ (24)

whose stability regions R are defined according to [13] as

$R_i=\left\{z\in ℂ:\left|L_i\left(z\right)\right|\le 1\right\},i=1,2,3$ (25)

While the small circles in Figures 1-3 represent zeros of (22), (23) and (24), the plus signs represent their poles. Thus the regions $R_i\left(i=1,2,3\right)$ in the complex plane which satisfy (25) for the stability functions (22), (23) and (24) represent the stability regions of the new block methods. From Figures 1, 2 and 3, the derived block methods are not A-stable since there are poles of stability functions in the left-half of the complex plane. However, the block HEIs (12), (16) and (20) are $L_0$ -stable ([14] [15]) as they satisfy

$\underset{z\le 0}{\max }\left|L_i\left(z\right)\right|\le 1,z\in ℝ,i=1,2,3\text{and}\underset{z\to -\infty }{\lim }L\left(z\right)\text{}=0$

4. Implementation and Conclusion

4.1. Implementation

The performance of the newly constructed Hybrid Explicit Integrators are tested by solving the following numerical examples:

Example 1

Consider the system of initial value problem (IVP)

$y^{\prime }=\left[\begin{array}{cc}-10004y_1& 10000y_2^4\\ y_1& -y_2\left(1+y_2^3\right)\end{array}\right]$

$y\left(0\right)={\left[1,1\right]}^{\text{T}},h=0.1,x\in \left[0,20\right]$

Exact solution

$y\left(x\right)={\left[{\text{e}}^{-4x},{\text{e}}^{-x}\right]}^{\text{T}}$

Figure 1. The region of absolute stability of HEI (12).

Figure 2. The region of absolute stability of HEI (16).

Figure 3. The region of absolute stability of HEI (20).

Example 2: Lotka Volterra equation

Consider the predator-prey model

$y^{\prime }=\left[\begin{array}{cc}1.2y_1& -0.6y_1{y}_2\\ -0.8y_2& 0.3y_1{y}_2\end{array}\right]$

$y\left(0\right)={\left[2,1\right]}^{\text{T}},h=0.001,x\in \left[0,20\right]$

Example 3: Robertson chemical equation

${y^{\prime }}_1=-0.04y_1+10000y_2{y}_3$

${y^{\prime }}_2=0.04y_1-10000y_2{y}_3-3.0\times {10}^7{y}_2^2$

${y^{\prime }}_3=3.0\times {10}^7{y}_2^2$

$y\left(0\right)={\left[1,0,0\right]}^{\text{T}},0\le x\le 400,h=0.001$

Solving Examples 1, 2 and 3, we used MATLAB software to obtain numerical solutions of the non-linear systems. The solution curves are given in Figures 4-6 respectively.

Figure 4. The solution curve of Example 1 with HEI (12).

Figure 5. The solution curve of Example 2 with HEI (16).

Figure 6. The solution curve of Example 3 with HEI (20).

Table 2. Maximum norm error and number of function evaluations for Example 1, $h=0.1$ , $x\in \left[0,20\right]$ .

Method	‖e‖∞ (y1)	‖e‖∞ (y2)	NFE
HEI (12), k = 2	3.14 × 10−5	7.82 × 10−6	600
ode23s (MATLAB)	1.12 × 10−5	2.30 × 10−6	≈1650
Butcher-Hojjati [16]	1.80 × 10−6	4.10 × 10−7	≈800

NFE = total number of function evaluations over the integration interval. For ode23s, NFE is the value returned by the solver’s output stats field under MaxStep = 0.1. The exact solution $y\left(x\right)={\left[{\text{e}}^{-4x},{\text{e}}^{-x}\right]}^{\top }$ is used as the reference.

Interpretation: The results show that HEI (12) achieves a maximum norm error of order 10⁻⁵ for both solution components at h = 0.1, requiring exactly 600 function evaluations, determined by the fixed block structure of the method (200 steps × 3 function evaluations per block). Under the same step-size constraint, ode23s requires approximately 1650 function evaluations due to its internal stage evaluations and error control overhead, at roughly 3 - 4 times the computational cost, while achieving a moderately lower error of order 10⁻⁵. This confirms the practical advantage of HEI (12): it delivers competitive accuracy at significantly reduced computational cost, consistent with the theoretical benefit of the explicit block structure. The comparison with Butcher-Hojjati is retained from Table 3 for context; as noted in the original manuscript, that method achieves higher accuracy owing to its second-derivative formulation, at a higher implementation cost.

4.2. Comparison of Solutions of HEI (12) with Those of Butcher and Hojjati for Problem 3 (Robertson Chemical Equation)

Table 3. Comparison of HEI (12) with Butcher and Hojjati [16].

x	step size	HEI (12)	Butcher and Hojjati [16]
0.4	0.001	9.851721136170910E−01	9.851721138620630E−01
		3.386395351200000E−05	3.386395379595400E−05
		1.479402242939700E−02	1.479402213590220E−02
4	0.001	9.055186789743650E−01	9.055186784344190E−01
		2.240475688100000E−05	2.240475693804370E−05
		9.445891626875500E−02	9.445891599170860E−02
40	0.001	7.158270688413850E−01	7.158270698910200E−01
		9.185534769000000E−06	9.185534641631410E−06
		2.841637456238570E−01	2.841637507954150E−01
400	0.001	4.505186684871060E−01	4.505186908340870E−01
		3.222901442000000E−06	3.222901061260970E−06
		5.494781086114610E−01	5.494782035239040E−01

There is a fair agreement between the results of Scheme (12) and that of [16]. However, the method derived by [16] have an advantage over Scheme (12) since it’s second derivative method.

4.3. Conclusion

In this study, we derived block Hybrid Explicit Integrators (HEIs) (12), (16) and (20) for the solution of stiff ODEs. These block methods are shown to be $L_0$ -stable for k = 2, 3, 4. They are also shown to be consistent and zero stable hence they are convergent methods. From the solution curves in Figures 1-3, it is evident that the derived HEIs have tendency to produce solutions for stiff ODEs with adequate accuracy. The quantitative validation in Table 2 confirms that HEI (12) achieves a maximum norm error of order 10⁻⁵ on Example 1 at h = 0.1 with only 600 function evaluations, compared to approximately 1650 evaluations required by ode23s under the same step-size policy, demonstrating that the proposed explicit block structure delivers competitive accuracy at substantially reduced computational cost. Thus, the derived block methods (12), (16) and (20) compete satisfactorily with the ode23s solver in MATLAB. Hence, we recommend them for solving stiff non-linear systems.

Conflicts of Interest

Authors have declared that no competing interests exist.

References

[1]	Ahmad, R.R. and Yaacob, N. (2005) Third-Order Composite Runge-Kutta Method for Stiff Problems. International Journal of Computer Mathematics, 82, 1221-1226.[CrossRef]
[2]	Wells, D.R. (1982) Multirate Linear Multistep Methods for the Solution of Systems of Ordinary Differential Equations.  Ph.D. Thesis, University of Illinois.
[3]	Onumanyi, P.A., Awoyemi, D., Jator, S.N. and Sirisena U.W. (1994) New Linear Multi-Step Methods with Continuous Coefficients for First Order Initial Value Problems. Journal of the Nigerian Mathematical Society, 13, 37-51.
[4]	Onumanyi, P., Sirisena, U.W. and Jator, S.N. (1999) Continuous Finite Difference Approximations for Solving Differential Equations. International Journal of Computer Mathematics, 72, 15-27.[CrossRef]
[5]	Sirisena, U.W., Kumleng, G.M. and Yahaya, Y.A. (2004) A New Butcher Type Two-Step Block Hybrid Multistep Method for Accurate and Efficient Parallel Solution of ODEs. Abacus Mathematics Series, 31, 1-7.
[6]	Xu, Y. and Zhao, J.J. (2010) Estimation of Longest Stability Interval for a Kind of Explicit Linear Multistep Methods. Discrete Dynamics in Nature and Society, 2010, Article ID: 912691.[CrossRef]
[7]	Jibunoh, C.C. (2018) An Exponential Method of Variable Order for General Non-Linear (Stiff and Nonstiff) ODE Systems. Journal of the Nigerian Mathematical Society, 37, 175-201. https://ojs.ictp.it/jnms/index.php/jnms/article/view/388
[8]	Sirisena, U.W., Luka, S.I. and Yakubu, S.Y. (2020) A Class of Explicit Integrators with O-Grid Interpolation for Solving Non-Linear Systems of First Order Odes. Journal of Advances in Mathematics and Computer Science, 35, 106-118.[CrossRef]
[9]	Fatunla, S.O. (1995) A Class of Block Methods for Second Order IVPs∗:. International Journal of Computer Mathematics, 55, 119-133.[CrossRef]
[10]	Henrici, P. (1962) Discrete Variable Methods in Ordinary Differential Equations. John Wiley and Sons.
[11]	Adedire, O., Sirisena, U.W. and Ndam, J.N. (2020) Implementation of an A-Stable Block Hybrid Method for the Solution of Bazykin Ecological Model and Trambouze Reactor Model. New Trends in Mathematical Science, 2, 1-8.[CrossRef]
[12]	Akinfenwa, O., Jator, S.N. and Yao, N.M. (2012) On the Stability of Continuous Block Backward Differentiation Formula for Solving Stiff Ordinary Differential Equations. Journal of Modern Methods in Numerical Mathematics, 3, 50-58.[CrossRef]
[13]	Hairer, E. and Wanner, G. (1996) Solving Ordinary Differential Equations II: Stiff and Differential Algebraic Problems. 2nd Edition, Springer.[CrossRef]
[14]	Akinfenwa, O.A., Jator, S.N. and Yao, N.M. (2013) Continuous Block Backward Differentiation Formula for Solving Stiff Ordinary Differential Equations. Computers & Mathematics with Applications, 65, 996-1005.[CrossRef]
[15]	Cash, J.R. (1984) Two New Finite Difference Schemes for Parabolic Equations. SIAM Journal on Numerical Analysis, 21, 433-446.[CrossRef]
[16]	Butcher, J.C. and Hojjati, G. (2005) Second Derivative Methods with RK Stability. Numerical Algorithms, 40, 415-429.[CrossRef]