Abstract
Several numerical techniques for solving initial value problems arise in physical and natural sciences. In many cases, these problems require numerical treatment to achieve the required solution. However, in today’s modern era, numerical algorithms must be costeffective with suitable convergence and stability features. At least the fifthorder convergent twostep optimized hybrid block method recently proposed in the literature is formulated in this research work with its variable stepsize approach for numerically solving first and higherorder initialvalue problems in ordinary differential equations. It has been constructed using a continuous approximation achieved through interpolation and collocation techniques at two intrastep points chosen by optimizing the local truncation errors of the main formulae. The theoretical analysis, including order stars for the relative stability, is considered. Both fixed and variable stepsize approaches are presented to observe the superiority of the latter approach. When tested on challenging differential systems, the method gives better accuracy, as revealed by the efficiency plots and the error distribution tables, including the machine time measured in seconds.
1 Introduction
The initial value problems of the following form are the most frequently used problems in several fields of studies:
where
The nonlinearity and stiffness of some models make it more difficult for applied mathematicians to devise efficient ways of obtaining approximate solutions with sufficient accuracy in a short amount of time. Various numerical methods, such as explicit and implicit Runge–Kutta (RK)type methods, diagonally implicit RK methods, singly diagonally implicit RK methods, linear multistep methods, including Adams–Bashforth/Moulton methods, backward and forward differential methods, multiderivative methods, rational/nonlinear methods, trigonometrically fitted methods, and hybrid block methods, have been created in the past and recent literature in search of better accuracy and timeefficiency. The block techniques have been prevalent among the scientific community due to their selfstarting characteristics and ability to avoid overlapping piecewise solutions. These methods, which include both primary and supplemental procedures, can be used to obtain an approximate solution at multiple intrastep points at the same time, as shown in refs [11,12, 13,14]. Several optimal block techniques have recently been developed in the literature to solve first and higher order ordinary differential equations, including one, two, four, five, six, and sevenstep block techniques. Nonetheless, only a small fraction of these approaches are considered for practical purposes owing to their limitations on the accuracy and convergence features. Therefore, we have been motivated to devise some strategies for the block methods to reduce functional evaluations (FEs) and computational costs. This is achieved in the present work by the formulation of an adaptive version of a block method available in the literature. As long as the block methods are concerned, there are many applications wherein these methods are helpful. For example, the nonlinear logistic growth model is often used in population dynamics, the mass spring damper system in physics, resistorinductorcapacitor series circuits in electronics, the Prothero–Robinson problem in chemistry, the periodic orbit system in quantum mechanics, and many other fields.
We try to formulate a variable stepsize version of the twostep optimal block method (OBM) with reasonably regular use and growing preference for efficient and robust block methods. It may be noted that the constant stepsize version is also reformulated here for a quick reference. Although several researchers have recently developed some of the optimized block techniques [15,16,17], they are either computationally expensive, have a lower order of convergence, or are only suitable for a specific class of initial value problems. In addition, most of them have not been represented in their RK form and reformulation versions. The unique feature related to the reformulated version of a block method was first proposed by Ramos in ref. [18] whose algorithm is being employed herein for the variable stepsize formulation and discussing its relative stability. The improvement of the method considering the variable step size implementation in the present research work is an advance in the performance of the strategy as adopted in ref. [19].
The present study is designed as follows: Section 2 contains the derivation of the optimized twostep hybrid block method with two intrastep points, where the optimization strategy is explained indepth, as well as the reformulation of the method, including its implicit RK structure. Sections 3 and 4 include the order stars (relative stability) and variable stepsize approach, respectively. In Section 5, various challenging nonlinear stiff models from different fields of engineering and science are considered, where numerical results are achieved using both constant and variable stepsize approaches of the twostep optimized hybrid block method and some existing methods with similar properties. Finally, Section 6 concludes the research findings with some recommendations for future study.
2 Derivation of the optimized block method
In this section, we derive the twostep optimized block technique with two intrastep points, where both intrastep points are optimized from the local truncation error (LTE) of the main formula. Let us consider the partition
where
Consider two intrastep points,
Solution of the above linear system gives values of the six unknown parameters that are not shown here for brevity. These parameters will determine the coefficients of the polynomial
where
To obtain the twostep block method, we evaluate
where
At this stage, there are two parameters (
Substituting these values in (10) and (11), we obtain:
Thus, the two optimized parameters (intrastep points:
It is worth to be noted that the above strategy concerning the optimization of the intrastep points was first introduced in ref. [20]. The pseudocode for the above method is provided in Algorithm 1 under Appendix A. It is worth noting that the reformulation of the optimized block method produces substantial savings in the computational cost. The idea of reformulation is based on the strategy proposed by Ramos [18]. The reduction in the computational time comes from the fact that the number of FEs is reduced as follows:
The above reformulation is later abbreviated as RBlock in the forthcoming sections. Moreover, the implicit RK structure of the optimized block method (13) is also presented through the usual Butcher tableau as follows:
0  0  0  0  0  0 


















1 





1 





In the section of numerical simulations, the abbreviation RKBlock is used for the method given in the above Butcher tableau.
3 Order stars
It has been commonly observed in the literature that the stability properties of various block methods generated through various techniques and steps are rarely investigated for their application to the theory of order stars (see refs [21,22]). Because there has been little discussion of this problem, we were motivated to perform a thorough analysis of the stability theory via order stars of the development of a new RK method and its economical implementation, which were taken into account in the current research work. These order stars play a crucial role in determining the relative stability of numerical techniques, as indicated by the order stars. Despite the fact that it may be used to examine regions of linear absolute stability of a
The
In available literature for the theory of order stars, two of its types are commonly used and referred to as the first and second type. Mathematically,
There are always some distinguished features observed among the three triplets mentioned above. However, these features can be visualized graphically with the help of different shades given to their plotted regions. In this way, the region of growth of relative stability for the set
When the above discussed regions are plotted in
Definition 1
(
The plot of order stars for the optimized block method is shown in Figure 2, where the shaded regions depict the triplet
4 A variable stepsize formulation
With the help of an embeddedtype procedure, this section discusses the formulation of the optimized block method (13) as a variable stepsize integrator. The procedure considers the combination of two methods of order
where
This estimation corresponds to the local error of the lower order formula, since
The aforementioned method is used because of the number of FEs to avoid any additional cost incurred by the lower order method. We have also used this secondorder method for each method chosen in the present research work for the numerical simulations carried out under the variable stepsize approach. A local tolerance TOL is predefined during the implementation process by the user. If the estimated error EST is greater than tolerance (TOL), the algorithm automatically adjusts the step size from the previous value (
where
On the other hand, when
5 Numerical dynamics with results and discussion
This section describes how the twostep optimized block method given in (13) performs while solving various types of mild and highly stiff application problems we come across in various fields from physical sciences. The method does not require several initial conditions or a predictor. While
Twostep at least fifthorder OBM (Block) with two intrastep points as given in (13).
Twostep at least fifthorder OBM in reformulated version (RBlock) with two intrastep points.
Twostep sixthorder OBM in RK form (RKBlock) with two intrastep points.
Onestep at least fifthorder OBM appeared in ref. [17].
Laguerre polynomial hybrid block method (LPHBM) of at least fifthorder appeared in ref. [28].
Fullyimplicit RKtype fifthorder method called (RadauI) appeared in [29, p. 226].
Implicit RKtype sixthorder Lobatto method (LobIII A) appeared in ref. [29, p. 228].
Problem 1
We consider the following stiff system of firstorder ODEs taken from ref. [30]:
where
In Problem 1, a twodimensional system of ODEs is simulated with both fixed and variable stepsize approaches of the twostep optimized block method (13). In Table 1, the tolerance

FE  Approach 

RMS  Mean 

25  125  Fixed 



25  125  Variable 



N  FE  Approach 

RMS  Mean 

112  1,060  Fixed 



112  1,060  Variable 



Problem 2
Consider another twodimensional stiff problem [31]:
where the exact solution is
The numerical simulations are performed for Problem 2 under the variable stepsize approach setting the tolerance to
Method 

RMSE  Mean  FE 

Block 



65 
RBlock 



65 
RKBlock 



108 
OBM 



65 
LPHBM 



135 
RadauI 



252 
LobIIIA 



90 
Method 

RMSE  Mean  FE 

Block 



125 
RBlock 



125 
RKBlock 



204 
OBM 



125 
LPHBM 



360 
RadauI 



260 
LobIIIA 



170 
Problem 3
The fourdimensional periodic orbit system taken from ref. [32] is considered:
where the exact solution over the interval
The numerical simulations are performed for Problem 3 under the variable stepsize approach setting the tolerance to
Method 

RMSE  Mean  FE 

Block 



55 
RBlock 



55 
RKBlock 



84 
OBM 



55 
Cheby 



295 
RadauI 



56 
LobIIIA 



70 
Method 

RMSE  Mean  FE 

Block 



230 
RBlock 



230 
RKBlock 



342 
OBM 



230 
LPHBM 



26320 
RadauI 



228 
LobIIIA 



285 
Problem 4
The Prothero–Robinson problem [33]:
where
The numerical simulations are performed for the Prothero–Robinson problem given in Problem 4 under the variable stepsize approach setting the tolerance to

Block  RBlock  RKBlock  OBM  LPHBM  RadauI  LobIIIA 

10 























Method  MaxErr  LastErr  Norm  RMSE  FE 

Block 




105 
RBlock 




105 
RKBlock 




126 
OBM 




105 
LPHBM 




225 
RadauI 




84 
LobIIIA 




105 
Method  MaxErr  LastErr  Norm  RMSE  FE 

Block 




210 
RBlock 




210 
RKBlock 




252 
OBM 




210 
LPHBM 




660 
RadauI 




168 
LobIIIA 




210 
Method  MaxErr  LastErr  Norm  RMSE  FE 

Block 




435 
RBlock 




435 
RKBlock 




522 
OBM 




435 
LPHBM 




2040 
RadauI 




348 
LobIIIA 




435 
Problem 5
The highly oscillatory problem [34]:
where the exact solution is
The numerical simulations are performed for the highly oscillatory problem given in Problem 5 under the variable stepsize approach setting the tolerance to
N  Block  RBlock  RKBlock  OBM  LPHBM  RadauI  LobIIIA 

10 























Method  MaxErr  LastErr  Norm  RMSE  FE 

Block 




65 
RBlock 




65 
RKBlock 




168 
OBM 




65 
LPHBM 




110 
RadauI 




112 
LobIIIA 




140 
Method  MaxErr  LastErr  Norm  RMSE  FE 

Block 




110 
RBlock 




110 
RKBlock 

