# A Family ofA-Stable Optimized Hybrid Block Methods for Integrating Stiff Differential Systems

Authors: Rajat Singla; Gurjinder Singh; Higinio Ramos; V. Kanwar
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this article, a family of one-step hybrid block methods having two intrastep points is developed for solving first-order initial value stiff differential systems that occur frequently in science and engineering. In each method of the family, an intrastep point controls the order of the main method and a second one has a control over the stability features of the method. The approach used to develop the class ofA-stable methods is based on interpolation and collocation procedures. The methods exhibit hybrid nature and produce numerical solutions at several points simultaneously. These methods can also be formulated as Runge-Kutta (RK) methods. Comparisons between the RK and block formulations of the proposed methods reveal a better performance of the block formulation in terms of computational efficiency. Furthermore, the efficiency of the methods is improved when they are formulated as adaptive step-size solvers using an error-control approach. Some methods of the proposed class have been tested to solve some well-known stiff differential systems. The numerical experiments show that the proposed family of methods performs well in comparison with some of the existing methods in the scientific literature.


---

## Body

## 1. Introduction

Differential equations are used to formulate most physical phenomena in mathematical terms. Due to the lack of availability of analytical methods, in many cases, they must be solved approximately using numerical methods [1–46]. Many of these applications often result into stiff ordinary differential equations (ODEs) that model atmospheric phenomena, chemical kinetics (for instance, explosions), engineering control systems, lasers, mechanics, electronic circuits, or molecular dynamics [1, 2]. It is well known that only a restricted number of differential equations can be solved analytically. In the absence of explicit solutions to these problems, numerous discrete numerical methods have been developed that accomplish the required task by providing numerical solutions. This article deals with the development of an implicit family of numerical integrators for solving the following initial value problem numerically (see the monographs [40, 41])(1)Z′t=Ft,Zt;Zt0=Z0,where t∈t0,tN,Z:t0,tN⟶ℝm,F:t0,tN×ℝm⟶ℝm, assuming that the conditions to have a unique solution hold.In the scientific literature, numerous numerical integrators are available to solve the differential system (1). For instance, the Runge-Kutta (RK) and multistep methods are two well-defined classes of integrators that have been widely used for this purpose. Among the last ones, the backward differentiation formulas (BDFs) are used to solve numerically stiff problems. Apart from them, some new approaches like hybrid and block methods were introduced by researches, which were recently used for solving these problems. Initially, Milne [3] proposed the block methods to overcome some of the drawbacks of predictor-corrector methods. Later on, Rosser [4] extended Milne’s concept to general purpose codes. Block methods have been shown to be more efficient in terms of accuracy and reduced number of function evaluations, thus requiring lower computational cost. Another advantage of using block methods is that they are self-starting. The block methods produce the numerical solution of problem (1) at several points simultaneously, since their structure consists of main and additional methods that are applied successively in block intervals. Many researchers have elongated the existing block methods by adding off-step points into it which are also referred to as hybrid points. Hybrid methods were constructed as a bridge between extrapolation and substitution methods. These formulas were simultaneously proposed by Gragg and Stetter [17], Butcher [18], and Gear [19]. These methods also have the characteristic of Runge-Kutta methods of using appropriate data at points other than the end points of the interval, along with certain properties of linear multistep methods. The main characteristic of hybrid methods is that they overcome the zero-stability barrier which restricts the order in the class of zero-stable linear multistep methods. Due to having this feature, hybrid methods are sometimes called modified linear multistep methods. For the development and application of hybrid and block methods on different types of ODEs one can see the pioneering works [6–16, 45, 46] and the references therein.Many authors use the hybrid block methods to solve ODEs numerically with arbitrary values of the intermediate points; then Ramos et al. proposed for the first time an optimization strategy to get the best choice of values of intrastep points by achieving a higher order of accuracy in the main formulas of the method at some points of the block (one can see the procedure adopted in [14] concerning second-order problems).This paper consists of the hybrid formulation of a class of single-step methods, involving two transitional points related to parametersr,s,0<r,s<1. The appropriate value of r is obtained by following an optimization strategy based on the minimization of the local truncation error of the main formula at the final point of the block. Going a little further, a novel approach is used to get appropriate values for the intermediate point s by imposing the conditions of A-stability to the resulting method obtained earlier. Therefore, we get a family of A-stable methods. The proposed family of methods was also formulated as Runge-Kutta methods. Comparisons have been made between the block formulation and the RK formulation for solving well-known differential systems, showing that the proposed family exhibits a better performance. In the RK formulation, one of the methods belongs to the LobattoIIIA family as a special case. In many practical applications, it is difficult to solve stiff ODEs, since the choice of many integrators for these problems is limited by the requirements of numerical stability. In deriving an efficient numerical integrator by preserving or improving the existing efficiency, the challenge is to enhance the stability region. Therefore, the aim is to develop a new class of implicit hybrid methods with two intermediate points in which one point controls the order of the main formula, while the other point controls the size of the absolute stability region of the proposed numerical method. Further, in order to improve the performance of the methods, an adaptive step-size formulation has been discussed using an embedded-type procedure. In the embedded-type approach, two methods of different orders are executed simultaneously. The lower order scheme is used to estimate the local error at each step of integration and the higher order method is used to advance the integration process. In this article, we have developed lower order methods in order to formulate some of the proposed methods in adaptive step-size mode. The lower order methods have been designed in such a way that they share the same function evaluations used by the main methods.The paper is summarized as follows: Section2 presents the development of one-step hybrid block methods having at least fourth orders of convergence with two intrastep points. In Section 3, the formulation of a class of A-stable methods and their main characteristics have been discussed. Section 4 deals with some particular cases of the developed family. In Section 5 the formulation in adaptive step-size mode is discussed. Section 6 deals with the Runge-Kutta formulation of the proposed family of A-stable methods. Comparisons have been made between the block and the RK formulations by performing some numerical experiments. Some numerical experiments were also performed to assess its efficiency in comparison with some other existing methods in the scientific literature. Finally, some conclusions were presented in Section 8.


## 2. Development of a Family of Hybrid Block Methods

The aim is to construct a family of one-step hybrid block methods for solving the differential system (1) numerically. Before further proceeding, system (1) is assumed to have a unique differentiable solution denoted byZt. Consider the interval of interest t0,tN discretized as follows:(2)tn=t0+nh,n=0,1,2,…,N.With constant step-sizeh=tn+1−tn,n=0,1,2,…,N−1. The numerical solution of (1) at tn will be denoted as Zn≈Ztn. For the sake of simplicity, the derivation of the method has been done for m=1, and thus the differential equation results in z′t=ft,zt. In the construction, a purely interpolation and collocation approach is considered. Let us approximate the true solution of problem (1) in tn,tn+1 by an interpolating polynomial of the form(3)zt≈Θt=∑j=04bjΨjt,where bj are the coefficients to be determined and Ψjt=t−tnj represents the basis polynomial functions. To make the proposed method of hybrid nature, we consider two points tn+r=tn+rh,tn+s=tn+sh such that 0<r,s<1. The unknown coefficients can be evaluated by imposing the following interpolating and collocation conditions:(i)
ztn≃zn=Θtn,(ii)
z′tn+k≃fn+k=Θ′tn+k,wherefn+k=ftn+k,zn+k and zn+k≃ztn+k for k=0,r,s,1. A system of five equations with five unknowns bj,j=0,1,2,3,4, is obtained, which can be solved by any computer algebra system (CAS). We have solved the system with the help of Mathematica and got the values of the coefficients:(4)b0=zn,b1=fn,b2=−−1+rr−s−1+sr+s+rsfn+s2r2−r+sfn+1−−1+sfn+r2h−1+rrr−s−1+ss+−1+rr2fn+s2h−1+rrr−s−1+ss,b3=−1+rr−s−1+s1+r+sfn+rs−r2+s2fn+1−s−1+s2fn+r3h2−1+rrr−s−1+ss+r−1+r2fn+s3h2−1+rrr−s−1+ss,b4=−−1+rr−s−1+sfn+rr−ssfn+1+−1+ssfn+r−−1+rrfn+s4h3−1+rrr−s−1+ss.Then, by inserting these values and changing the variablet⟶tn+mh, (3) can be written in a convenient form as(5)ztn+mh≈Θtn+mh=α0zn+hβ0fn+βrfn+r+βsfn+s+β1fn+1.And(6)α0=1,β0=m−3m3+12rs+4m21+r+s−6mr+s+rs12rs,βr=m23m2+6s−4m1+s12−1+rrr−s,βs=m2−6r+m4−3m+4r12r−s−1+ss,β1=m23m2+6rs−4mr+s12−1+r−1+s.Evaluating expression (3) at m=1, the approximation to the theoretical solution of problem (1) at t=tn+1 is obtained as(7)zn+1=zn+hAr,s12−1+rrr−s−1+ss.And(8)Ar,s=−1+rr−s−1+s1−2s+r−2+6sfn+rr−ss3+6rs−4r+sfn+1+−1+ss−1+2sfn+r+1−2r−1+rrfn+s.The appropriate values ofr are obtained from the local truncation error of the formula in (4), which is obtained after expanding in Taylor series around tn, which results in(9)ℒztn+1,h=3−5s+5r−1+2sy5tnh51440+Oh6.Equating the principal term of this error to zero and solving forr, we get the value of r in terms of s, that is,(10)r=−3+5s5−1+2s.Now, in order to get approximations ofztn+r,ztn+s,ztn+1, we use the value in (4) and take m=r,s,1 in (3), thus obtaining the one-parameter family given by(11)zn+r=zn+hA1s−1+ss−2+5s3+10−1+ss1500−1+2s3,zn+s=zn+hA2s12−1+s−3+5s−2+5s3+10−1+ss,zn+1=zn+hA3s12−1+ss−3+5s−2+5s3+10−1+ss.And(12)A1s=−3+5sfn+s3−5s2−2+5s−7+15s+fn−1+s−2+5s3+10−1+ss−21+5s40+s−113+100s+s−fn+13−5s23+10−1+ss3+5s−3+4s+25fn+r1−2s2−1+s−33+5s44+s−101+80s,A2s=sfn−1+s−2+5s3+10−1+ss−18+s46+5s−7+2s+fn+s−3+5s−2+5s−18+s62+15s−5+2s+s2−125fn+r−2+s−1+s−1+2s3+fn+1−3+5s3+10−1+ss6+5s−3+2s,A3s=−6fn+s+fn−1+s−2+5s3+10−1+ss1+2s−4+5s+s25−fn+s+5fn+r1−2s4−1+s+fn+1−3+5s3+10−1+ss3+2s−6+5s.In the next section the characteristics of the methods and restrictions on the intermediate values will be derived to make the above family of methods a family of A-stable methods.


## 3. Characteristics of the Methods

In this section, some basic characteristics of the family in (5) will be discussed, such as order of accuracy, zero-stability analysis, conditions on the parameter s to reach A-stability, and linear stability analysis.


### 3.1. Consistency and Order of Accuracy

Firstly, consider the difference operatorℒ¯ related to the proposed family of methods (7) (see [41])(13)ℒ¯zt,h=∑kρ¯kzt+kh−h∑kδ¯kz′t+kh,k=0,r,s,1,where ρ¯k and δ¯k are precisely the corresponding coefficients of the method. Assuming that zt has continuous derivatives of sufficiently high order and expanding zt+kh and z′t+kh in Taylor series about t we get(14)ℒ¯zt,h=Φ¯0zt+Φ¯1hz′t+Φ¯2h2z″t+…+Φ¯phpzpt+….Definition 1.
The linear difference operator (8) and the corresponding block method are said to be of order p if, in (9),(15)Φ¯0=Φ¯1=Φ¯2=…=Φ¯p=0andΦ¯p+1≠0.
Here,Φ¯i’s are vectors of coefficients and Φ¯p+1 is called the vector of error constants, and thus,(16)ℒ¯zt,h=Φ¯p+1hp+1zp+1t+Ohp+2.
For the proposed methods, we getΦ¯0=Φ¯1=Φ¯2=Φ¯3=Φ¯4=0.
And(17)Φ¯5=6−25s+25s23750000−1+2s5,−1+s3s3240−1+2s,0T.
This indicates that the proposed methods have at least orderp=4, which makes them consistent, as p≥1 is the sufficient feature in this concern (see [41]).


### 3.2. Zero-Stability

In order to address the zero-stability analysis of the proposed methods (7), let us consider h⟶0, which reduces the proposed methods to zn+r=zn,zn+s=zn,zn+1=zn. Hence, the methods can be written as(18)N∧ζ=D¯N∧ζ−1,where(19)D¯=001001001.And(20)N∧ζ=zn+r,zn+s,zn+1T,N∧ζ−1=zn+r−1,zn+s−1,znT.The first characteristic polynomial for the proposed methods (7) is(21)ρη=detI3η−D¯=η2η−1.As the roots ofρη=0 satisfy ηj≤1 and the root of modulus one is simple then the proposed methods are zero-stable (see [20]).


### 3.3. Convergence

As per establishment of Lax equivalence theorem, a method is said to be convergent if and only if the method is consistent and zero-stable, which implies that the methods in (5) are convergent.


### 3.4. Linear Stability

The concept of linear stability is quite different from that of zero-stability, since it is not practical to deal withh⟶0. In practice we always deal with some finite value of h>0. In what follows, by applying the theory of linear stability analysis, the region of parameter s is carried out which makes the family of A-stable methods, among which some of them can be considered as almost ℒ-stable methods.Definition 2.
The region of absolute stability,S, is defined as [21](22)S=H∈ℂ:ϕH<1,where ϕH is the stability function. If ℂ−⊆S, that is, all values of the left-half complex plane are included in the region of absolute stability, then the method is said to be A-stable.Theorem 1.
For eachs∈I=5−5/10,2/5∪5+5/10,1, the corresponding method in (7) is A-stable.Proof.
Consider the well-known Dahlquist’s test equation, given as(23)z′t=αzt,Reα<0.
The theoretical solution of problem (15) is zt=keαt, which tends to zero as t⟶∞. If we solve problem (15) numerically with any of the proposed methods, it is expected that the solution obtained will behave similarly to the analytical solution. Our aim is to obtain the values of s for which the numerical method reproduces the behaviour of the true solution of the test problem. In that regard, the proposed methods are applied to problem (15), thus obtaining(24)zn+r=zn+hB1s−1+ss−2+5s3+10−1+ss1500−1+2s3,zn+s=zn+hB2s12−1+s−3+5s−2+5s3+10−1+ss,zn+1=zn+hB3s12−1+ss−3+5s−2+5s3+10−1+ss,where(25)B1s=α−3+5szn+s3−5s2−2+5s−7+15s+zn−1+s−2+5s3+10−1+ss−21+5s40+s−113+100s+s−zn+13−5s23+10−1+ss3+5s−3+4s+25zn+r1−2s2−1+s−33+5s44+s−101+80s,B2s=αszn−1+s−2+5s3+10−1+ss−18+s46+5s−7+2s+zn+s−3+5s−2+5s−18+s62+15s−5+2s+s2−125zn+r−2+s−1+s−1+2s3+zn+1−3+5s3+10−1+ss6+5s−3+2s,B3s=α−6zn+s+zn−1+s−2+5s3+10−1+ss1+2s−4+5s+s25−zn+s+5zn+r1−2s4−1+s+zn+1−3+5s3+10−1+ss3+2s−6+5s.
After eliminatingzn+r,zn+s from (16) and changing the variables hα⟶H and zn+1⟶κ,zn⟶1, we get the characteristic equation given by(26)36H−18H2+3H3+120−48H+6H23−5s−120−48H+6H2sκ=36H+18H2+3H3+120+72H+18H2+2H33−5s−120+72H+18H2+2H3s,
provided that(27)−1+s≠0,s≠0,−1+2s≠0,−3+5s≠0,−2+5s≠0,3−10s+10s2≠0.
In this case the stability functionϕH is readily obtained from the characteristic equation as(28)ϕH=120−240s+12H6−15s+5s2+6H23−9s+5s2+H32−7s+5s2120−240s+H3s−3+5s+6H21+s−5s2+12H−4+5s+5s2.
In order to make the methodsA-stable, we apply the limit when H⟶±∞ to (16) and get(29)Gs=limH⟶±∞ϕH=5s2−7s+2s5s−3.
Now, in order to get the values ofs for which the resulting method is A-stable, we imposed on Gs that(30)Gs≤1.
The values ofs∈0,1 for which inequality (18) is satisfied are s∈I1=5−5/10,1/2∪5+5/10,1. For these values of s, the corresponding method is A-stable.
Since we are interested in intermediate points, we have to discard those values ofs for which the value of r as given in (4) becomes greater than 1 or less than 0. Hence, we solve the system of inequalities(31)−3+5s5−1+2s>1,−3+5s5−1+2s<0.
We get that fors∈I2=2/5,3/5−1/2 it is r<0 or r>1, which implies that I=I1−I2 is the set for which the corresponding method becomes A-stable, with 0<r,s<1. Figure 1 shows this set.Figure 1
Plot of the values ofs satisfying the inequality given in (20). (a) Plot of function Gs in (19). (b) Plot of values of s on real line.
(a)(b)Definition 3.
If the method isA-stable and the function ϕH tends to 0 as H⟶±∞, where ϕ is the stability function of the method, the method is said to be ℒ-stable.Theorem 2.
There are no values ofr and s with 0<r<s<1 for which the method in (14) is ℒ-stable.Proof.
According to (15), the method is not defined for s=2/5,1. For other values of s we have that limH⟶±∞ϕH≠0, and thus the method is not ℒ-stable.
From the above analysis we have seen that fors=5−5/10≈0.276393 the resulting method is A-stable and the stability region is exactly the left half plane (one can see the plot of the stability region in Figure 2(b)). As s increases from 5−5/10 to 2/5 the stability region also increases (see Figures 2(b)–2(e)). Similarly, for s=5+5/10≈0.723607 the stability region is exactly the left half plane (see Figure 2(g)) and this region increases as s approaches towards 1 (see Figures 2(g)–2(k)). The following plots in Figure 2 show the absolute stability regions of some methods of the proposed family for different values of s.Figure 2
Absolute stability regions of some methods of the proposed family for different values ofs. (a) For s=2/10∉I. (b) For s=5−5/10∈I. (c) For s=3/10∈I. (d) For s=35/100∈I. (e) For s=39/100∈I. (f) For s=7/10∉I. (g) For s=5+5/10∈I. (h) For s=8/10∈I. (i) For s=9/10∈I. (j) For s=21/22∈I. (k) For s=99/100∈I.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)


## 3.1. Consistency and Order of Accuracy

Firstly, consider the difference operatorℒ¯ related to the proposed family of methods (7) (see [41])(13)ℒ¯zt,h=∑kρ¯kzt+kh−h∑kδ¯kz′t+kh,k=0,r,s,1,where ρ¯k and δ¯k are precisely the corresponding coefficients of the method. Assuming that zt has continuous derivatives of sufficiently high order and expanding zt+kh and z′t+kh in Taylor series about t we get(14)ℒ¯zt,h=Φ¯0zt+Φ¯1hz′t+Φ¯2h2z″t+…+Φ¯phpzpt+….Definition 1.
The linear difference operator (8) and the corresponding block method are said to be of order p if, in (9),(15)Φ¯0=Φ¯1=Φ¯2=…=Φ¯p=0andΦ¯p+1≠0.
Here,Φ¯i’s are vectors of coefficients and Φ¯p+1 is called the vector of error constants, and thus,(16)ℒ¯zt,h=Φ¯p+1hp+1zp+1t+Ohp+2.
For the proposed methods, we getΦ¯0=Φ¯1=Φ¯2=Φ¯3=Φ¯4=0.
And(17)Φ¯5=6−25s+25s23750000−1+2s5,−1+s3s3240−1+2s,0T.
This indicates that the proposed methods have at least orderp=4, which makes them consistent, as p≥1 is the sufficient feature in this concern (see [41]).


## 3.2. Zero-Stability

In order to address the zero-stability analysis of the proposed methods (7), let us consider h⟶0, which reduces the proposed methods to zn+r=zn,zn+s=zn,zn+1=zn. Hence, the methods can be written as(18)N∧ζ=D¯N∧ζ−1,where(19)D¯=001001001.And(20)N∧ζ=zn+r,zn+s,zn+1T,N∧ζ−1=zn+r−1,zn+s−1,znT.The first characteristic polynomial for the proposed methods (7) is(21)ρη=detI3η−D¯=η2η−1.As the roots ofρη=0 satisfy ηj≤1 and the root of modulus one is simple then the proposed methods are zero-stable (see [20]).


## 3.3. Convergence

As per establishment of Lax equivalence theorem, a method is said to be convergent if and only if the method is consistent and zero-stable, which implies that the methods in (5) are convergent.


## 3.4. Linear Stability

The concept of linear stability is quite different from that of zero-stability, since it is not practical to deal withh⟶0. In practice we always deal with some finite value of h>0. In what follows, by applying the theory of linear stability analysis, the region of parameter s is carried out which makes the family of A-stable methods, among which some of them can be considered as almost ℒ-stable methods.Definition 2.
The region of absolute stability,S, is defined as [21](22)S=H∈ℂ:ϕH<1,where ϕH is the stability function. If ℂ−⊆S, that is, all values of the left-half complex plane are included in the region of absolute stability, then the method is said to be A-stable.Theorem 1.
For eachs∈I=5−5/10,2/5∪5+5/10,1, the corresponding method in (7) is A-stable.Proof.
Consider the well-known Dahlquist’s test equation, given as(23)z′t=αzt,Reα<0.
The theoretical solution of problem (15) is zt=keαt, which tends to zero as t⟶∞. If we solve problem (15) numerically with any of the proposed methods, it is expected that the solution obtained will behave similarly to the analytical solution. Our aim is to obtain the values of s for which the numerical method reproduces the behaviour of the true solution of the test problem. In that regard, the proposed methods are applied to problem (15), thus obtaining(24)zn+r=zn+hB1s−1+ss−2+5s3+10−1+ss1500−1+2s3,zn+s=zn+hB2s12−1+s−3+5s−2+5s3+10−1+ss,zn+1=zn+hB3s12−1+ss−3+5s−2+5s3+10−1+ss,where(25)B1s=α−3+5szn+s3−5s2−2+5s−7+15s+zn−1+s−2+5s3+10−1+ss−21+5s40+s−113+100s+s−zn+13−5s23+10−1+ss3+5s−3+4s+25zn+r1−2s2−1+s−33+5s44+s−101+80s,B2s=αszn−1+s−2+5s3+10−1+ss−18+s46+5s−7+2s+zn+s−3+5s−2+5s−18+s62+15s−5+2s+s2−125zn+r−2+s−1+s−1+2s3+zn+1−3+5s3+10−1+ss6+5s−3+2s,B3s=α−6zn+s+zn−1+s−2+5s3+10−1+ss1+2s−4+5s+s25−zn+s+5zn+r1−2s4−1+s+zn+1−3+5s3+10−1+ss3+2s−6+5s.
After eliminatingzn+r,zn+s from (16) and changing the variables hα⟶H and zn+1⟶κ,zn⟶1, we get the characteristic equation given by(26)36H−18H2+3H3+120−48H+6H23−5s−120−48H+6H2sκ=36H+18H2+3H3+120+72H+18H2+2H33−5s−120+72H+18H2+2H3s,
provided that(27)−1+s≠0,s≠0,−1+2s≠0,−3+5s≠0,−2+5s≠0,3−10s+10s2≠0.
In this case the stability functionϕH is readily obtained from the characteristic equation as(28)ϕH=120−240s+12H6−15s+5s2+6H23−9s+5s2+H32−7s+5s2120−240s+H3s−3+5s+6H21+s−5s2+12H−4+5s+5s2.
In order to make the methodsA-stable, we apply the limit when H⟶±∞ to (16) and get(29)Gs=limH⟶±∞ϕH=5s2−7s+2s5s−3.
Now, in order to get the values ofs for which the resulting method is A-stable, we imposed on Gs that(30)Gs≤1.
The values ofs∈0,1 for which inequality (18) is satisfied are s∈I1=5−5/10,1/2∪5+5/10,1. For these values of s, the corresponding method is A-stable.
Since we are interested in intermediate points, we have to discard those values ofs for which the value of r as given in (4) becomes greater than 1 or less than 0. Hence, we solve the system of inequalities(31)−3+5s5−1+2s>1,−3+5s5−1+2s<0.
We get that fors∈I2=2/5,3/5−1/2 it is r<0 or r>1, which implies that I=I1−I2 is the set for which the corresponding method becomes A-stable, with 0<r,s<1. Figure 1 shows this set.Figure 1
Plot of the values ofs satisfying the inequality given in (20). (a) Plot of function Gs in (19). (b) Plot of values of s on real line.
(a)(b)Definition 3.
If the method isA-stable and the function ϕH tends to 0 as H⟶±∞, where ϕ is the stability function of the method, the method is said to be ℒ-stable.Theorem 2.
There are no values ofr and s with 0<r<s<1 for which the method in (14) is ℒ-stable.Proof.
According to (15), the method is not defined for s=2/5,1. For other values of s we have that limH⟶±∞ϕH≠0, and thus the method is not ℒ-stable.
From the above analysis we have seen that fors=5−5/10≈0.276393 the resulting method is A-stable and the stability region is exactly the left half plane (one can see the plot of the stability region in Figure 2(b)). As s increases from 5−5/10 to 2/5 the stability region also increases (see Figures 2(b)–2(e)). Similarly, for s=5+5/10≈0.723607 the stability region is exactly the left half plane (see Figure 2(g)) and this region increases as s approaches towards 1 (see Figures 2(g)–2(k)). The following plots in Figure 2 show the absolute stability regions of some methods of the proposed family for different values of s.Figure 2
Absolute stability regions of some methods of the proposed family for different values ofs. (a) For s=2/10∉I. (b) For s=5−5/10∈I. (c) For s=3/10∈I. (d) For s=35/100∈I. (e) For s=39/100∈I. (f) For s=7/10∉I. (g) For s=5+5/10∈I. (h) For s=8/10∈I. (i) For s=9/10∈I. (j) For s=21/22∈I. (k) For s=99/100∈I.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)


## 4. Some Particular Cases

In this section we consider some specific values ofs∈I. In order to satisfy the condition 0<r<s<1, we take the values of s in 5+5/10,1.


### 4.1. Method fors=5+5/10

From (4), we have that r=5−5/10. Using these values in (5) we get the proposed hybrid method, whose coefficients are listed in Table 1.Table 1
Coefficients of the methodOHBM1.
α0β0βrβsβ1zn+r111+5/12025−5/12025−135/120−1+5/120zn+s111−5/12025+135/12025+5/120−1−5/120zn+111/125/125/121/12After expanding in Taylor series aroundtn the associated linear operator in (6), we get values Φ¯i=0,i=0,…,4, and the vector of error constants is Φ¯5=−1/60005,1/60005,0T with(32)ℒ¯zt,h=−160005h5z5t+Oh6160005h5z5t+Oh611512000h7z7t+Oh8.Thus, the method has at least orderp=4. Considering the test problem (15), we get that the stability function is given by(33)ϕH=120+60H+12H2+H3120−60H+12H2−H3.AndH=hα. The stability region is plotted in Figure 2(g).This particular case can also be derived by setting the first two terms of the local truncation error in (5) equal to zero. We get the system of nonlinear equations(34)3−5s+5r−1+2s=0,1−s−s2+r2−1+2s+r−1+s+2s2=0.This system corresponds to two curves in thers-plane which are symmetric with respect to the line r=s, having a unique solution with the condition 0<r<s<1, given by r=5−5/10,s=5+5/10. The method will be named as OHBM1.


### 4.2. Method fors=21/22

Now, the corresponding value ofr is r=39/100. The resulting coefficients are listed in Table 2.Table 2
Coefficients of the methodOHBM2.
α0β0βrβsβ1zn+r1296543/200000064961/219600−2924207/900000033040683/122000000zn+s18435/692125359375/9499347287/396−37044/81191zn+11599/49142500000/443207729282/39123−53/122After the expansion in Taylor series aroundtn, we get that Φ¯i=0,i=0,…,4, with the vector of error constants Φ¯5=4488095313/17600000000000,−3087/8245811200,0T and(35)ℒ¯zt,h=448809531317600000000000h5z5t+Oh6−30878245811200h5z5t+Oh6−37915840000h6z6t+Oh7.Hence, the method has at least orderp=4. Applying the method to the test problem (15), it results in the stability function(36)ϕH=52800+21852H+3006H2+61H352800−30948H+7554H2−819H3.AndH=hα. The stability region is plotted in Figure 2(j). The method will be named as OHBM2.


### 4.3. Method fors=99/100

The corresponding value ofr is r=39/98. The coefficients of the resulting method are listed in Table 3.Table 3
Coefficients of the methodOHBM3.
α0β0βrβsβ1zn+r12348515/155296683378661/11182388−5389515625/375429723919160037/13882582zn+s13233439/26000000427023163413/7416890000001120383/386800−307584783/118000000zn+112881/2316611529602/200256032500000/861597−307/118After the expansion in Taylor series aroundtn we get that Φ¯i=0,i=0,…,4, with the vector of error constants Φ¯5=4060958967/14462732748800,−323433/78400000000000,0T and(37)ℒ¯zt,h=406095896714462732748800h5z5t+Oh6−32343378400000000000h5z5t+Oh6−190170560000h6z6t+Oh7.Hence, the method has at least orderp=4. Applying the method to the test problem (15), it results in the stability function(38)ϕH=235200+94788H+12114H2+59H3235200−140412H+34926H2−3861H3.AndH=hα. The stability region is plotted in Figure 2(k). The method will be named as OHBM3.


## 4.1. Method fors=5+5/10

From (4), we have that r=5−5/10. Using these values in (5) we get the proposed hybrid method, whose coefficients are listed in Table 1.Table 1
Coefficients of the methodOHBM1.
α0β0βrβsβ1zn+r111+5/12025−5/12025−135/120−1+5/120zn+s111−5/12025+135/12025+5/120−1−5/120zn+111/125/125/121/12After expanding in Taylor series aroundtn the associated linear operator in (6), we get values Φ¯i=0,i=0,…,4, and the vector of error constants is Φ¯5=−1/60005,1/60005,0T with(32)ℒ¯zt,h=−160005h5z5t+Oh6160005h5z5t+Oh611512000h7z7t+Oh8.Thus, the method has at least orderp=4. Considering the test problem (15), we get that the stability function is given by(33)ϕH=120+60H+12H2+H3120−60H+12H2−H3.AndH=hα. The stability region is plotted in Figure 2(g).This particular case can also be derived by setting the first two terms of the local truncation error in (5) equal to zero. We get the system of nonlinear equations(34)3−5s+5r−1+2s=0,1−s−s2+r2−1+2s+r−1+s+2s2=0.This system corresponds to two curves in thers-plane which are symmetric with respect to the line r=s, having a unique solution with the condition 0<r<s<1, given by r=5−5/10,s=5+5/10. The method will be named as OHBM1.


## 4.2. Method fors=21/22

Now, the corresponding value ofr is r=39/100. The resulting coefficients are listed in Table 2.Table 2
Coefficients of the methodOHBM2.
α0β0βrβsβ1zn+r1296543/200000064961/219600−2924207/900000033040683/122000000zn+s18435/692125359375/9499347287/396−37044/81191zn+11599/49142500000/443207729282/39123−53/122After the expansion in Taylor series aroundtn, we get that Φ¯i=0,i=0,…,4, with the vector of error constants Φ¯5=4488095313/17600000000000,−3087/8245811200,0T and(35)ℒ¯zt,h=448809531317600000000000h5z5t+Oh6−30878245811200h5z5t+Oh6−37915840000h6z6t+Oh7.Hence, the method has at least orderp=4. Applying the method to the test problem (15), it results in the stability function(36)ϕH=52800+21852H+3006H2+61H352800−30948H+7554H2−819H3.AndH=hα. The stability region is plotted in Figure 2(j). The method will be named as OHBM2.


## 4.3. Method fors=99/100

The corresponding value ofr is r=39/98. The coefficients of the resulting method are listed in Table 3.Table 3
Coefficients of the methodOHBM3.
α0β0βrβsβ1zn+r12348515/155296683378661/11182388−5389515625/375429723919160037/13882582zn+s13233439/26000000427023163413/7416890000001120383/386800−307584783/118000000zn+112881/2316611529602/200256032500000/861597−307/118After the expansion in Taylor series aroundtn we get that Φ¯i=0,i=0,…,4, with the vector of error constants Φ¯5=4060958967/14462732748800,−323433/78400000000000,0T and(37)ℒ¯zt,h=406095896714462732748800h5z5t+Oh6−32343378400000000000h5z5t+Oh6−190170560000h6z6t+Oh7.Hence, the method has at least orderp=4. Applying the method to the test problem (15), it results in the stability function(38)ϕH=235200+94788H+12114H2+59H3235200−140412H+34926H2−3861H3.AndH=hα. The stability region is plotted in Figure 2(k). The method will be named as OHBM3.


## 5. Formulation as an Adaptive Step-Size Algorithm

To formulate efficiently any of the methods whose coefficients are listed in Tables1–3, an embedded-type procedure is adopted following the guidelines in Shampine et al. [22]. The strategy is similar to the one used with embedded Runge-Kutta methods. In the formulation, there is a simultaneous execution of two methods; one is the proposed method and the other one is a lower order method which is useful in estimating the local error at each step. The higher order method is used for advancing the integration process. The lower order method is considered in such a way that it uses the same evaluations of the function in the proposed method, in order not to increase the computational cost. The procedure is elaborated as follows:Firstly, a second approximation of the solution at the final point of the intervaltn,tn+1 is obtained using the lower order method (say of order q). Denoting this approximation by zn+1∗ an estimate of the local error εn+1 is(39)εn+1=ztn+h−zn+1∗,where zt represents the true solution. Now, the proposed method of higher order (say p) is applied to get the approximation zn+1 in the current step. Then, we have(40)ztn+h−zn+1=Ohp+1,And thus(41)Esn=zn+1−zn+1∗=ztn+h−zn+1∗−ztn+h−zn+1=εn+1+Ohp+1.The estimationEsn is easy to compute since the local error of the lower order method, εn+1, is an Ohq+1, which dominates in (31) for small enough values of h. By comparing it to the more precise solution zn+1, the error at zn+1∗ can be estimated. In this procedure, we are not able to calculate the value of the local error at each step of integration exactly, but an assumption has been made that it is rather smaller than the estimated local error [24].Below we present the lower order formulas that will be considered according to the strategy described, for each of the particular methods. For the method in Table1, we consider the third-order formula(42)zn+1∗=zn+h16fn+5+512fn+r+5−512fn+s.And local truncation errorLTE=z4tnh4/360+Oh5.For the method in Table2, we consider the third-order formula(43)zn+1∗=zn+h2202457fn+1004339123fn+r+4750072657fn+s.And local truncation errorLTE=−53z4tnh4/26400+Oh5.For the method in Table3, we consider the third-order formula(44)zn+1∗=zn+h98011583fn+232897339417fn+r+197500861597fn+s.And local truncation errorLTE=−307z4tnh4/117600+Oh5. Note that the formulas in (20)–(22) use the same function evaluations as the proposed methods. To advance the integration, the required step-size for the next step will be determined with the help of the error estimate, Esn. According to the usual strategy with predefined tolerance, TOL, the algorithm automatically provides the new step-size given by(45)hn=ξhoTOL‖Esn‖1/q+1,where q is the order of the lower order method, and hn and ho represent the new step-size and the previous step-size. In expression (23), a safety factor 0<ξ<1 is introduced to avoid the failure of steps.Some conditions were imposed on the new step-sizehn in the form of an If statement. This is used to avoid large fluctuations in the step-sizes: the step-size is not allowed to decrease more than hmin (minimum step-size allowed) or to increase more than hmax (maximum step-size allowed). This may be included in the implementation using an If statement:(46)Ifhmin≤hn≤hmaxthenho=hn.If‖Esn‖<TOL, then a successful step is considered and the strategy will be used successively in determining the step-size to be used for the next step. In the literature, different strategies are available for selecting the appropriate starting step hi (one can see Shampine et al. [24] and Watts [25]). There is also a simple way to select the starting step-size as discussed in Sedgwick [26], one selects a very small step-size, and then the algorithm will change the value of the chosen step-size according to the requirements of the problem.


## 6. Formulation as a Runge-Kutta Method

Vigo-Aguiar and Ramos describe in [33] the parallelism between some class of RK collocation methods and special family of multistep formulas (block formulation) with an explanation about the construction of these formulas. Any of the block methods in (5) can be formulated as an RK method. The method in Table 1 is formulated as an RK method taking the vector of stage values c=0,5−5/10,5+5/10,1T, which results in the following Butcher tableau:(47)The block methods in Tables2 and 3 can be also formulated as RK methods with stage vectors c=0,39/100,21/22,1T and c=0,39/98,99/100,1T and Butcher tableaux(48)When the method in Table1 is formulated into a RK structure it results in the LobattoIIIA method of order p=6 (see [40]). However, when other members of the family of block methods are formulated in the RK framework, they are novel ones since they have not been found in the present literature neither can be written as a linear combination of the four fundamental LobattoIIIA, IIIB, IIIC, or IIIC RK methods. As is known, RK methods perform well on stiff problems if they are stiffly accurate. The above formulated RK methods also exhibit the same property, consisting in that the numerical solution of the current step coincides with the last stage (see Hairer and Wanner [21]). The methods in Tables 2 and 3 formulated as RK methods are named as RK15 and RK25 and have order p=5. As we anticipated, using the theory of RK methods, it can be verified that all the methods of the family are A-stable and have order 5, except for the LobattoIIIA, which has order six. In order to verify that the block structure is more efficient than the RK structure, in this section the well-known problem named as “The Brusselator system” is solved using both approaches, as a block method (OHBM1, OHBM2, and OHBM3) and as a RK one (LobattoIIIA, RK15, and RK25). It should be noted that the RK ones required greater computational effort in all cases.


### 6.1. The Brusselator System

In this example the diffusion-free “Brusselator system” [20] is solved numerically, which was also solved in [30, 33, 34]. The problem consists in a system of two equations(49)z1′t=L+z12tz2t−M+1z1t,z10=z10,z2′t=Mz1t−z12tz2t,z20=z20,where L and M are positive real constants [35]. The critical point of the system is z1∗,z2∗=L,M/L. The values M=3,L=1 are assumed for the numerical computations, with the initial values z10=1.5,z20=3. As suggested by Butcher and Podhaisky in [34], the interval of integration considered is 0,20. In the absence of a true solution, the following reference solution (see [30]) at the end point of the interval, tN=20, will be used to compare errors(50)Z1=0.4986370712683478483331816235,Z2=4.5967803494520111826429803773.The problem has been solved with step-sizesh=0.1,0.05,0.025,0.0125. The data in Table 4 shows that the method in block formulation performs better in comparison with the RK formulation in terms of computational cost, while the errors are similar, as expected. In Figure 3, we have plotted the efficiency curves concerning the absolute global errors on logarithmic scale versus the computational times. It is clear that the formulations as block methods (OHBM1, OHBM2, and OHBM3) perform better than the corresponding RK ones (LobattoIIIA (p=6), RK15p=5, RK25p=5). The abbreviation Ctime is used for the CPU time.Table 4
Comparison of block structure and RK structure for “the Brusselator system.”
FevalOHBM1LobattoIIIAOHBM2RK15OHBM2RK258004.10035×10−94.10035×10−91.14453×10−61.14453×10−61.28289×10−61.28289×10−6Ctime0.6406250.843750.406250.6406250.468750.76562516001.68977×10−101.68977×10−103.50600×10−83.50600×10−83.94564×10−83.94564×10−8Ctime1.218751.687500.8281251.109380.812501.2656332002.6068×10−122.61391×10−121.10406×10−91.10406×10−91.24294×10−91.24295×10−9Ctime2.453134.015631.718752.31251.703132.3750064004.52971×10−144.44089×10−143.46034×10−113.46061×10−113.89626×10−113.89644×10−11Ctime4.859386.656253.453134.81253.390634.73438Figure 3
Efficiency curves of the absolute global errors versus CPU times for “the Brusselator system.”


## 6.1. The Brusselator System

In this example the diffusion-free “Brusselator system” [20] is solved numerically, which was also solved in [30, 33, 34]. The problem consists in a system of two equations(49)z1′t=L+z12tz2t−M+1z1t,z10=z10,z2′t=Mz1t−z12tz2t,z20=z20,where L and M are positive real constants [35]. The critical point of the system is z1∗,z2∗=L,M/L. The values M=3,L=1 are assumed for the numerical computations, with the initial values z10=1.5,z20=3. As suggested by Butcher and Podhaisky in [34], the interval of integration considered is 0,20. In the absence of a true solution, the following reference solution (see [30]) at the end point of the interval, tN=20, will be used to compare errors(50)Z1=0.4986370712683478483331816235,Z2=4.5967803494520111826429803773.The problem has been solved with step-sizesh=0.1,0.05,0.025,0.0125. The data in Table 4 shows that the method in block formulation performs better in comparison with the RK formulation in terms of computational cost, while the errors are similar, as expected. In Figure 3, we have plotted the efficiency curves concerning the absolute global errors on logarithmic scale versus the computational times. It is clear that the formulations as block methods (OHBM1, OHBM2, and OHBM3) perform better than the corresponding RK ones (LobattoIIIA (p=6), RK15p=5, RK25p=5). The abbreviation Ctime is used for the CPU time.Table 4
Comparison of block structure and RK structure for “the Brusselator system.”
FevalOHBM1LobattoIIIAOHBM2RK15OHBM2RK258004.10035×10−94.10035×10−91.14453×10−61.14453×10−61.28289×10−61.28289×10−6Ctime0.6406250.843750.406250.6406250.468750.76562516001.68977×10−101.68977×10−103.50600×10−83.50600×10−83.94564×10−83.94564×10−8Ctime1.218751.687500.8281251.109380.812501.2656332002.6068×10−122.61391×10−121.10406×10−91.10406×10−91.24294×10−91.24295×10−9Ctime2.453134.015631.718752.31251.703132.3750064004.52971×10−144.44089×10−143.46034×10−113.46061×10−113.89626×10−113.89644×10−11Ctime4.859386.656253.453134.81253.390634.73438Figure 3
Efficiency curves of the absolute global errors versus CPU times for “the Brusselator system.”


## 7. Results and Discussion

In this section, the proposed methods using fixed step-size and adaptive step-size formulations are implemented on a set of well-known differential systems existing in the literature. Some of the considered differential systems model different chemical reactions where some variables change rapidly whereas others vary slowly, indicating stiffness. In the first subsection an example is solved numerically by the proposed methods. The results presented in Table5 show that in some cases the methods with large stability regions performed better even though the main formulas have less order. The main purpose of this subsection is to establish that the stability region also plays a vital role in solving the differential equations numerically, along with the order of the method. In the second subsection, some well-known problems are solved numerically with the proposed methods. The results presented in Tables 6–12 using the strategy discussed in Section 5 and some well-known numerical integrators reveal a better performance of the proposed methods. Note that the data corresponding to the adaptive step-size version of the method in Table 2 are not presented here as they are similar to those obtained with the adaptive step-size version of the method in Table 3, since they have similar stability characteristics and order.


### 7.1. Nomenclature

The following nomenclature has been used in the article:(i)
N: Number of integration steps(ii)
FEVAL’s: Number of function evaluations, including derivatives(iii)
MaxErr: Maximum absolute errors along the integration interval given by(51)MaxErr=max1≤k≤mmax0≤n≤Nzktn−zkn.(iv)
h: Fixed step-size(v)
hi: Initial step-size(vi)
TOL: Predefined user tolerancewherezktn and zkn represent the analytical and computed kth-component of the solution of a m-system differential problem at tn.


### 7.2. Implementation with Fixed Step-Size

This experiment is performed using the proposed methods in Tables1–3, which clearly shows that the method with larger absolute stability regions performs better in case of large step-sizes.


### 7.3. An Example

We consider the following inhomogeneous stiff problem whose eigenvalues of the coefficient matrix areλ1=−1,λ2=−1000(52)z1′t=−2z1t+z2t+2sint,z2′t=998z1t−999z2t+999cost−sint.To activate the transient component the initial data arez10=0,z20=0, which yield the exact solution as(53)z1t=−1999e−t+1999e−1000t+sint,z2t=−1999e−t−998999e−1000t+cost.The problem is solved numerically using the proposed integrators OHBM1, OHBM2, and OHBM3 over the intervalt∈0,10 with various fixed step-sizes h=1,0.4,0.2. The numerical data are presented in Table 5. Figure 4 shows the exact solution and the numerical solutions with OHBM1, OHBM2, and OHBM3 for different step-sizes. One can observe clearly that the proposed method OHBM1 produces oscillations in the discrete solution for step-sizes h≥0.005, whereas the methods OHBM2 and OHBM3 match the exact solution even with large step-sizes.Table 5
Numerical results for problem.
hNMethodMaxErrFEVALs110OHBM19.75312×10−140OHBM27.00562×10−2OHBM31.19781×10−20.425OHBM19.40822×10−1100OHBM27.39160×10−2OHBM37.41226×10−30.250OHBM18.86038×10−1200OHBM25.41812×10−2OHBM31.90508×10−4Figure 4
Exact and discrete solutions ofz2t for problem using the proposed methods with different values of h. (a) For h=0.4. (b) For h=0.2. (c) For h=0.1. (d) For h=0.05.
(a)(b)(c)(d)


### 7.4. Implementation with Adaptive Step-Size

This section deals with the implementation of the proposed methods in variable step-size. Comparisons were made with some well-known numerical integrators as described as follows:(1)
ODE15s: ODE15s is an integrated solver coded in MATLAB, specially created for solving stiff-systems. The concept used in the construction of this solver is based on the numerical differentiation formulas (NDFs) of orders 1 to 5. It is a variable step, variable order (VSVO) IVP solver. For a detailed survey on this code, one can see [22, 23].(2)
RADAU: This solver consists of variable order (1, 5, 9, 13) with an adapting step-size control feature. The method is based on the implicit Runge-Kutta methods Radau-IIa. For the experiment purposes, we used the MATLAB code developed by Hairer (see the MatlabStiff package inhttp://www.unige.ch/hairer/software.html).(3)
VSOHBM1 and VSOHBM3: These are the proposed hybrid schemes whose coefficients are listed in Tables1 and 3, which were implemented with the variable step-size strategy described in Section 5. The codes of these new schemes were developed in Mathematica. At each step of integration, a system of linear/nonlinear equations has to be solved due to their implicit nature. In order to solve the resulting nonlinear systems the command FindRoot is used, while for linear systems the command NSolve was considered.Note: The RADAU and ODE15s integrators use different strategies for changing the step-size than the proposed schemes. RADAU and ODE15s consider variable step as well as variable order methods, while the proposed schemes are implemented using the embedded procedure discussed in Section5. In the implementation of RADAU and ODE15s codes, AbsTOL and RelTOL have been used, which are abbreviations for absolute and relative tolerances, while for the proposed schemes the TOL will be used. For numerically solving the test problems, we consider TOL=AbsTOL=RelTOL. The software MATLABR2009b was used for the implementation of RADAU and ODE15s on test problems whereas the proposed methods were implemented in Wolfram Mathematica 11, version 11.0.1.0. The runs have been done using a laptop with processor i3-4030U CPU @ 1.90 GHz, Win10.


## 7.1. Nomenclature

The following nomenclature has been used in the article:(i)
N: Number of integration steps(ii)
FEVAL’s: Number of function evaluations, including derivatives(iii)
MaxErr: Maximum absolute errors along the integration interval given by(51)MaxErr=max1≤k≤mmax0≤n≤Nzktn−zkn.(iv)
h: Fixed step-size(v)
hi: Initial step-size(vi)
TOL: Predefined user tolerancewherezktn and zkn represent the analytical and computed kth-component of the solution of a m-system differential problem at tn.


## 7.2. Implementation with Fixed Step-Size

This experiment is performed using the proposed methods in Tables1–3, which clearly shows that the method with larger absolute stability regions performs better in case of large step-sizes.


## 7.3. An Example

We consider the following inhomogeneous stiff problem whose eigenvalues of the coefficient matrix areλ1=−1,λ2=−1000(52)z1′t=−2z1t+z2t+2sint,z2′t=998z1t−999z2t+999cost−sint.To activate the transient component the initial data arez10=0,z20=0, which yield the exact solution as(53)z1t=−1999e−t+1999e−1000t+sint,z2t=−1999e−t−998999e−1000t+cost.The problem is solved numerically using the proposed integrators OHBM1, OHBM2, and OHBM3 over the intervalt∈0,10 with various fixed step-sizes h=1,0.4,0.2. The numerical data are presented in Table 5. Figure 4 shows the exact solution and the numerical solutions with OHBM1, OHBM2, and OHBM3 for different step-sizes. One can observe clearly that the proposed method OHBM1 produces oscillations in the discrete solution for step-sizes h≥0.005, whereas the methods OHBM2 and OHBM3 match the exact solution even with large step-sizes.Table 5
Numerical results for problem.
hNMethodMaxErrFEVALs110OHBM19.75312×10−140OHBM27.00562×10−2OHBM31.19781×10−20.425OHBM19.40822×10−1100OHBM27.39160×10−2OHBM37.41226×10−30.250OHBM18.86038×10−1200OHBM25.41812×10−2OHBM31.90508×10−4Figure 4
Exact and discrete solutions ofz2t for problem using the proposed methods with different values of h. (a) For h=0.4. (b) For h=0.2. (c) For h=0.1. (d) For h=0.05.
(a)(b)(c)(d)


## 7.4. Implementation with Adaptive Step-Size

This section deals with the implementation of the proposed methods in variable step-size. Comparisons were made with some well-known numerical integrators as described as follows:(1)
ODE15s: ODE15s is an integrated solver coded in MATLAB, specially created for solving stiff-systems. The concept used in the construction of this solver is based on the numerical differentiation formulas (NDFs) of orders 1 to 5. It is a variable step, variable order (VSVO) IVP solver. For a detailed survey on this code, one can see [22, 23].(2)
RADAU: This solver consists of variable order (1, 5, 9, 13) with an adapting step-size control feature. The method is based on the implicit Runge-Kutta methods Radau-IIa. For the experiment purposes, we used the MATLAB code developed by Hairer (see the MatlabStiff package inhttp://www.unige.ch/hairer/software.html).(3)
VSOHBM1 and VSOHBM3: These are the proposed hybrid schemes whose coefficients are listed in Tables1 and 3, which were implemented with the variable step-size strategy described in Section 5. The codes of these new schemes were developed in Mathematica. At each step of integration, a system of linear/nonlinear equations has to be solved due to their implicit nature. In order to solve the resulting nonlinear systems the command FindRoot is used, while for linear systems the command NSolve was considered.Note: The RADAU and ODE15s integrators use different strategies for changing the step-size than the proposed schemes. RADAU and ODE15s consider variable step as well as variable order methods, while the proposed schemes are implemented using the embedded procedure discussed in Section5. In the implementation of RADAU and ODE15s codes, AbsTOL and RelTOL have been used, which are abbreviations for absolute and relative tolerances, while for the proposed schemes the TOL will be used. For numerically solving the test problems, we consider TOL=AbsTOL=RelTOL. The software MATLABR2009b was used for the implementation of RADAU and ODE15s on test problems whereas the proposed methods were implemented in Wolfram Mathematica 11, version 11.0.1.0. The runs have been done using a laptop with processor i3-4030U CPU @ 1.90 GHz, Win10.


## 8. Numerical Experiments Performed with VSOHBM3

### 8.1. Robertson’s Chemistry Problem

Firstly, a well-known classical problem of chemical sciences which is formulated from the kinetics of an autocatalytic reaction (Hairer & Wanner, 1996) is considered. For having detailing on the origin of the problem, one can see ‘Test Set for Initial Value Problem Solvers’ [27]. The problem is also studied by Cash [28, 29]. Problem was also solved with numerical integrators in [30, 32, 33]. It consists of a stiff system of three nonlinear ordinary differential equations(54)z1′t=−0.04z1t+104z2tz3t,z10=1,z2′t=0.04z1t−104z2tz3t−3×107z22t,z20=0,z3′t=3×107z2t2,z30=0.The problem is numerically solved by discretizing the interval0,40 as discussed in [28]. In the absence of analytical solution of this problem in the present literature, the errors were compared by considering the reference solution at the end point tN=40 which is as follows:(55)z1=0.7158270687194135,z2=9.185534764558135×10−6,z3=0.28416374574582.As described in [30], the numerical computations have been done by considering hi,TOL=10−ξ,10−ξ+4,ξ=2,3. The data presented in Table 6 show that the proposed algorithm VSOHBM3 performs better than the other numerical integrators.Table 6
Data for “Robertson’s chemistry problem.”
hiTOLMethodMaxErrFEVALs10−210−6RADAU7.9719×10−8300ODE15s1.0229×10−6158VSOHBM31.62210×10−910410−310−7RADAU1.4542×10−8321ODE15s3.6332×10−7212VSOHBM31.23246×10−10176


### 8.2. Van der Pol System

Another well-known system considered was the Van der Pol system in [37]. The system comprises two ordinary differential equations(56)z1′t=z2t,z10=2,z2′t=1−z12tz2t−z1tε,z20=−23+1081ε−2922187ε2−181419683ε3.The problem was solved numerically withε=10−1 over the interval of integration 0,0.55139 as in [37]. The analytical solution of the system is unknown due to which the errors were compared by the reference solution [30].(57)z1=1.5633739442300918,z2=−1.0000208318542727.The numerical results were calculated considering the parameters(58)hi,TOL=10−ξ,10−ξ+3,ξ=3,4,5.The results in Table7 reveal the outperformance of the proposed adaptive step-size scheme compared with the other integrators.Table 7
Data for problem.
hiTOLMethodMaxErrz1MaxErrz2NFEVALs10−310−6RADAU1.1000×10−71.5400×10−613119ODE15s2.3320×10−63.2750×10−62951VSOHBM38.78801×10−118.78801×10−883210−410−7RADAU1.2200×10−81.832×10−716148ODE15s3.5330×10−75.8800×10−74982VSOHBM31.81322×10−129.66674×10−10208010−510−8RADAU1.2900×10−92.0750×10−1018154ODE15s4.2330×10−86.9190×10−865108VSOHBM32.93099×10−145.30709×10−1135140


### 8.3. Gear’s Problem

In this example the problem suggested by Gear [38] is considered. The system consists of three ordinary differential equations:(59)z1′t=−0.013z1t−1000z1tz3t,z10=1,z2′t=−2500z2tz3t,z20=1,z3′t=−0.013z1t−1000z1tz3t−2500z2tz3t,z30=0.Similarly as in [38] the problem was integrated over the fixed interval 0,50. The numerical results in Table 8 have been obtained by considering(60)hi,TOL=10−ξ,10−ξ+3,ξ=1,2.Table 8
Data for “gear’s problem.”
hiTOLMethodMaxErrNFEVALs10−110−4RADAU2.8810×10−514103ODE15s2.1249×10−41832VSOHBM34.10906×10−762410−210−5RADAU1.4970×10−61185ODE15s1.5589×10−62339VSOHBM33.01563×10−8936The following reference solution at the end pointtN=50 was considered for computing the errors in the numerical solution [30].(61)z1=0.59765469806558128638,z2=1.40234340854787827842,z3=−1.8933865404351958485×10−6.The proposed algorithm again gives the better outcomes for this problem.


### 8.4. Stiff System of Three Equations

Consider the following first-order ODEs stiff system which was also solved in [39](62)z1′t=−1000z1t,z10=1,z2′t=−2000z2t,z20=1,z3′t=−3000z1tz2t,z30=1.The problem was solved numerically over the interval0,20. The errors in the numerical solution have been compared to the exact solution of the system(63)z1t=e−1000t,z2t=e−2000t,z3t=e−3000t.The numerical results were computed by considering(64)hi,TOL=10−ξ,10−ξ+1,ξ=1,2,3.The data given in Table9 show a superior performance of the proposed adaptive step-size method.Table 9
Data for the stiff system.
hiTOLMethodMaxErrNFEVALs10−110−2RADAU1.4820×10−124143ODE15s1.1900×10−24292VSOHBM33.20551×10−4114410−210−3RADAU1.2100×10−229170ODE15s2.6000×10−352104VSOHBM32.03312×10−5156010−310−4RADAU2.1000×10−337210ODE15s3.6169×10−5110180VSOHBM31.81507×10−82288


### 8.5. A Nonlinear System

Here, the following system of first-order ODEs is considered(65)z1′t=−200z2t2,z10=1,z2′t=−100z2t,z20=1.Fort∈0,20, while solving the problem numerically the parameters chosen with varying step-size and tolerances were(66)hi,TOL=10−ξ,10−ξ+1,ξ=2,3.The analytical solution of the problem is(67)z1t=e−200t,z2t=e−100t.The numerical data presented in Table10 show the good performance of the proposed variable step-size integrator.Table 10
Data for problem.
hiTOLMethodMaxErrNFEVALs10−210−3RADAU5.3000×10−325166ODE15s2.4000×10−35097VSOHBM32.58958×10−5135210−310−4RADAU1.8000×10−324159ODE15s2.9898×10−478126VSOHBM32.20958×10−62080


## 8.1. Robertson’s Chemistry Problem

Firstly, a well-known classical problem of chemical sciences which is formulated from the kinetics of an autocatalytic reaction (Hairer & Wanner, 1996) is considered. For having detailing on the origin of the problem, one can see ‘Test Set for Initial Value Problem Solvers’ [27]. The problem is also studied by Cash [28, 29]. Problem was also solved with numerical integrators in [30, 32, 33]. It consists of a stiff system of three nonlinear ordinary differential equations(54)z1′t=−0.04z1t+104z2tz3t,z10=1,z2′t=0.04z1t−104z2tz3t−3×107z22t,z20=0,z3′t=3×107z2t2,z30=0.The problem is numerically solved by discretizing the interval0,40 as discussed in [28]. In the absence of analytical solution of this problem in the present literature, the errors were compared by considering the reference solution at the end point tN=40 which is as follows:(55)z1=0.7158270687194135,z2=9.185534764558135×10−6,z3=0.28416374574582.As described in [30], the numerical computations have been done by considering hi,TOL=10−ξ,10−ξ+4,ξ=2,3. The data presented in Table 6 show that the proposed algorithm VSOHBM3 performs better than the other numerical integrators.Table 6
Data for “Robertson’s chemistry problem.”
hiTOLMethodMaxErrFEVALs10−210−6RADAU7.9719×10−8300ODE15s1.0229×10−6158VSOHBM31.62210×10−910410−310−7RADAU1.4542×10−8321ODE15s3.6332×10−7212VSOHBM31.23246×10−10176


## 8.2. Van der Pol System

Another well-known system considered was the Van der Pol system in [37]. The system comprises two ordinary differential equations(56)z1′t=z2t,z10=2,z2′t=1−z12tz2t−z1tε,z20=−23+1081ε−2922187ε2−181419683ε3.The problem was solved numerically withε=10−1 over the interval of integration 0,0.55139 as in [37]. The analytical solution of the system is unknown due to which the errors were compared by the reference solution [30].(57)z1=1.5633739442300918,z2=−1.0000208318542727.The numerical results were calculated considering the parameters(58)hi,TOL=10−ξ,10−ξ+3,ξ=3,4,5.The results in Table7 reveal the outperformance of the proposed adaptive step-size scheme compared with the other integrators.Table 7
Data for problem.
hiTOLMethodMaxErrz1MaxErrz2NFEVALs10−310−6RADAU1.1000×10−71.5400×10−613119ODE15s2.3320×10−63.2750×10−62951VSOHBM38.78801×10−118.78801×10−883210−410−7RADAU1.2200×10−81.832×10−716148ODE15s3.5330×10−75.8800×10−74982VSOHBM31.81322×10−129.66674×10−10208010−510−8RADAU1.2900×10−92.0750×10−1018154ODE15s4.2330×10−86.9190×10−865108VSOHBM32.93099×10−145.30709×10−1135140


## 8.3. Gear’s Problem

In this example the problem suggested by Gear [38] is considered. The system consists of three ordinary differential equations:(59)z1′t=−0.013z1t−1000z1tz3t,z10=1,z2′t=−2500z2tz3t,z20=1,z3′t=−0.013z1t−1000z1tz3t−2500z2tz3t,z30=0.Similarly as in [38] the problem was integrated over the fixed interval 0,50. The numerical results in Table 8 have been obtained by considering(60)hi,TOL=10−ξ,10−ξ+3,ξ=1,2.Table 8
Data for “gear’s problem.”
hiTOLMethodMaxErrNFEVALs10−110−4RADAU2.8810×10−514103ODE15s2.1249×10−41832VSOHBM34.10906×10−762410−210−5RADAU1.4970×10−61185ODE15s1.5589×10−62339VSOHBM33.01563×10−8936The following reference solution at the end pointtN=50 was considered for computing the errors in the numerical solution [30].(61)z1=0.59765469806558128638,z2=1.40234340854787827842,z3=−1.8933865404351958485×10−6.The proposed algorithm again gives the better outcomes for this problem.


## 8.4. Stiff System of Three Equations

Consider the following first-order ODEs stiff system which was also solved in [39](62)z1′t=−1000z1t,z10=1,z2′t=−2000z2t,z20=1,z3′t=−3000z1tz2t,z30=1.The problem was solved numerically over the interval0,20. The errors in the numerical solution have been compared to the exact solution of the system(63)z1t=e−1000t,z2t=e−2000t,z3t=e−3000t.The numerical results were computed by considering(64)hi,TOL=10−ξ,10−ξ+1,ξ=1,2,3.The data given in Table9 show a superior performance of the proposed adaptive step-size method.Table 9
Data for the stiff system.
hiTOLMethodMaxErrNFEVALs10−110−2RADAU1.4820×10−124143ODE15s1.1900×10−24292VSOHBM33.20551×10−4114410−210−3RADAU1.2100×10−229170ODE15s2.6000×10−352104VSOHBM32.03312×10−5156010−310−4RADAU2.1000×10−337210ODE15s3.6169×10−5110180VSOHBM31.81507×10−82288


## 8.5. A Nonlinear System

Here, the following system of first-order ODEs is considered(65)z1′t=−200z2t2,z10=1,z2′t=−100z2t,z20=1.Fort∈0,20, while solving the problem numerically the parameters chosen with varying step-size and tolerances were(66)hi,TOL=10−ξ,10−ξ+1,ξ=2,3.The analytical solution of the problem is(67)z1t=e−200t,z2t=e−100t.The numerical data presented in Table10 show the good performance of the proposed variable step-size integrator.Table 10
Data for problem.
hiTOLMethodMaxErrNFEVALs10−210−3RADAU5.3000×10−325166ODE15s2.4000×10−35097VSOHBM32.58958×10−5135210−310−4RADAU1.8000×10−324159ODE15s2.9898×10−478126VSOHBM32.20958×10−62080


## 9. Numerical Experiments Performed with VSOHBM1

### 9.1. Jacobi Elliptic Functions

In this example, the well-known system involving Jacobi elliptic functionssn,cn, and dn is discussed which also appeared in [36](68)sn′t=cntdnt,sn0=0,cn′t=−sntdnt,cn0=1,dn′t=−ksntcnt,dn0=1,where k=1/2. The analytical solution of the problem can be described in the form of series as(69)snt=2π1/2M∑p=0∞bp+1/21−b2p+1sin2p+1u,cnt=2π1/2M∑p=0∞bp+1/21+b2p+1cos2p+1u,dnt=π2M+2πM∑p=1∞bp1+b2pcos2pu.Hereb=e−π,u=πt/2M, and M=∫0π/2dϕ/1−1/2sin2ϕ≈1.85. The problem was solved numerically over the integration interval 0,50. For this problem, the parameters considered were(70)hi,TOL=10−ξ,10−ξ+2,ξ=1,2.The numerical results showed in Table11 elaborate the superiority of the new scheme.Table 11
Data for problem.
hiTOLMethodMaxErrNFEVALs10−110−3RADAU3.7×10−31301106ODE15s3.53×10−2210441VSOHBM11.38182×10−48132410−210−4RADAU9.8342×10−41521114ODE15s2.89×10−2359735VSOHBM15.37351×10−6139556


### 9.2. A Nonlinear System of Two Equations

Consider the nonlinear system of first-order ODEs taken from Rene [42] and also solved in [39](71)z1′t=z2t−z1t2−1+t,z10=1,z2′t=1−20z2t2−1+t2,z20=1.where t∈0,10. The analytical solution of the system is(72)z1t=11+t,z2t=1+t.And there are the eigenvalues of the Jacobian matrix−2/1+t and −401+t. The numerical computations have been obtained by varying the initial step-sizes and tolerances as(73)hi,TOL=10−ξ,10−ξ+4,ξ=1,2.The data computed in Table12 support that the proposed block method gives the best results in terms of accuracy and number of function evaluations.Table 12
Data for problem.
hiTOLMethodMaxErrFEVALs10−110−5RADAU8.9642×10−8151ODE15s6.2839×10−6134VSOHBM11.67505×10−148010−210−6RADAU3.8021×10−8157ODE15s1.4991×10−6172VSOHBM11.77636×10−15140


### 9.3. The Oregonator Problem

Lastly, the system consisting of three stiff nonlinear equations is considered which is formulated from the celebrated Belousov-Zhabotinsky reaction. The origin of the problem is discussed in the test set [27]. The differential system is(74)z1′t=αz2t+z1t1−βz1t−z2t,z2′t=z3t−1+z1tz2tα,z3′t=γz1t−z3t.Andα=77.27,β=8.375×10−6,γ=0.161 and initial values z0=1,2,3T. The numerical computations were done by considering the integration interval 0,360 and the numerical errors were compared by considering the reference solution at tN=360 as in the test set [27].(75)z1tN=0.1000814870318523×101,z2tN=0.1228178521549917×104,z3tN=0.1320554942846706×103.The problem was solved numerically by consideringhi,TOL=10−ξ,10−ξ+1,ξ=1,2,3,4. In Figure 5 we have plotted the efficiency curves concerning maximum absolute error versus the number of function evaluations for the methods considered. From Figure 5, one can easily judge that the proposed algorithm executes better than the other chosen integrators in terms of accuracy.Figure 5
Efficiency curves for “the Oregonator problem.”


## 9.1. Jacobi Elliptic Functions

In this example, the well-known system involving Jacobi elliptic functionssn,cn, and dn is discussed which also appeared in [36](68)sn′t=cntdnt,sn0=0,cn′t=−sntdnt,cn0=1,dn′t=−ksntcnt,dn0=1,where k=1/2. The analytical solution of the problem can be described in the form of series as(69)snt=2π1/2M∑p=0∞bp+1/21−b2p+1sin2p+1u,cnt=2π1/2M∑p=0∞bp+1/21+b2p+1cos2p+1u,dnt=π2M+2πM∑p=1∞bp1+b2pcos2pu.Hereb=e−π,u=πt/2M, and M=∫0π/2dϕ/1−1/2sin2ϕ≈1.85. The problem was solved numerically over the integration interval 0,50. For this problem, the parameters considered were(70)hi,TOL=10−ξ,10−ξ+2,ξ=1,2.The numerical results showed in Table11 elaborate the superiority of the new scheme.Table 11
Data for problem.
hiTOLMethodMaxErrNFEVALs10−110−3RADAU3.7×10−31301106ODE15s3.53×10−2210441VSOHBM11.38182×10−48132410−210−4RADAU9.8342×10−41521114ODE15s2.89×10−2359735VSOHBM15.37351×10−6139556


## 9.2. A Nonlinear System of Two Equations

Consider the nonlinear system of first-order ODEs taken from Rene [42] and also solved in [39](71)z1′t=z2t−z1t2−1+t,z10=1,z2′t=1−20z2t2−1+t2,z20=1.where t∈0,10. The analytical solution of the system is(72)z1t=11+t,z2t=1+t.And there are the eigenvalues of the Jacobian matrix−2/1+t and −401+t. The numerical computations have been obtained by varying the initial step-sizes and tolerances as(73)hi,TOL=10−ξ,10−ξ+4,ξ=1,2.The data computed in Table12 support that the proposed block method gives the best results in terms of accuracy and number of function evaluations.Table 12
Data for problem.
hiTOLMethodMaxErrFEVALs10−110−5RADAU8.9642×10−8151ODE15s6.2839×10−6134VSOHBM11.67505×10−148010−210−6RADAU3.8021×10−8157ODE15s1.4991×10−6172VSOHBM11.77636×10−15140


## 9.3. The Oregonator Problem

Lastly, the system consisting of three stiff nonlinear equations is considered which is formulated from the celebrated Belousov-Zhabotinsky reaction. The origin of the problem is discussed in the test set [27]. The differential system is(74)z1′t=αz2t+z1t1−βz1t−z2t,z2′t=z3t−1+z1tz2tα,z3′t=γz1t−z3t.Andα=77.27,β=8.375×10−6,γ=0.161 and initial values z0=1,2,3T. The numerical computations were done by considering the integration interval 0,360 and the numerical errors were compared by considering the reference solution at tN=360 as in the test set [27].(75)z1tN=0.1000814870318523×101,z2tN=0.1228178521549917×104,z3tN=0.1320554942846706×103.The problem was solved numerically by consideringhi,TOL=10−ξ,10−ξ+1,ξ=1,2,3,4. In Figure 5 we have plotted the efficiency curves concerning maximum absolute error versus the number of function evaluations for the methods considered. From Figure 5, one can easily judge that the proposed algorithm executes better than the other chosen integrators in terms of accuracy.Figure 5
Efficiency curves for “the Oregonator problem.”


## 10. Conclusions

In this article, a family of one-stepA-stable hybrid block methods using two intrastep points, namely, r and s, were developed. The values of intrastep point r have been specified in terms of s after optimizing the local truncation error of the approximate value zn+1. A suitable interval for the intrastep point s has been found over which the proposed family satisfies the condition of A-stability, making it worthy for solving stiff problems. Hence, in the family of methods, one intrastep point controls the order of the method using an optimization strategy and the second intrastep point controls the stability region. The development of the method is based on two well-known approaches viz. hybrid and block, coupled with an optimization strategy. The family of methods may be formulated as RK methods of four stages in which one of the methods from the family results to be the LobattoIIIA method. We have conducted numerical experiments with both formulations, that is, block and RK formulations of the proposed family. It is worth to mention here that this straightforward block formulation could be considered for other RK methods, claiming almost certainly that there will be a gain in computational cost. Some numerical results have been presented to show that the efficiency of the proposed methods with larger stability regions (OHBM2 and OHBM3) is acceptable in spite of having main formulas with less order. An improved version of the new proposed family considering an adaptive step-size formulation is also presented. In this regard, new methods of lower order are designed in order to formulate the proposed methods in adaptive step-size mode. Numerical experiments reveal that the new proposed schemes are a good alternative for solving the types of problems considered. In particular, the proposed VSOHBM1 and VSOHBM3 may be competitive with the RadauIIA of variable order (1, 5, 9, 13) and with ODE15s based on numerical differentiation formulas of orders 1 to 5.



---

*Source: 5576891-2022-05-16.xml*