Variable stepsize construction of a two-step optimized hybrid block method with relative stability

1 Introduction

The initial value problems of the following form are the most frequently used problems in several fields of studies:

where x ∈ [ a , b ] ⊂ R and w ( x ) , g ( x , w ( x ) ) ∈ R n . The assumption on the existence of a unique continuously differentiable solution w ( x ) under suitable conditions is made. The assumption confirms that problem (1) is well-posed. In every discipline of science, the aforementioned type of equations often arise. Several models relying on ordinary differential equations have been suggested to understand the dynamic behavior of the coronavirus pandemic, shown in refs [1,2,3], reaction-rate equations [4], exponential growth/decay [5], van der Pol oscillator [6], nonlinear corneal shape model [7], electrical, hydraulic and mechanical systems with variable mass [8], double pendulum [9], logistic growth [10], and many more.

The non-linearity 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, multi-derivative 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 time-efficiency. The block techniques have been prevalent among the scientific community due to their self-starting characteristics and ability to avoid overlapping piece-wise solutions. These methods, which include both primary and supplemental procedures, can be used to obtain an approximate solution at multiple intra-step 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 seven-step 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, resistor-inductor-capacitor 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 two-step 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 two-step hybrid block method with two intra-step points, where the optimization strategy is explained in-depth, 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 two-step 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 two-step optimized block technique with two intra-step points, where both intra-step points are optimized from the local truncation error (LTE) of the main formula. Let us consider the partition a = x 0 < x 1 < x 2 , ⋯ < x N = b on the integration interval [ a , b ] with stepsize Δ x = x k + 1 − x k , k = 0 , 1 , 2 , … , N − 1 , and assume that on a generic sub-interval [ x n , x n + 2 ] , the true solution w ( x ) of (1) can be approximated by a polynomial Q ( x ) in the following form:

where δ j ∈ R represent real unknown parameters. When Eq. (2) is differentiated, the following result is obtained:

Consider two intra-step points, x n + r = x n + r Δ x , x n + s = x n + s Δ x with 0 < r < 1 < s < 2 , to compute the approximate solution of the initial value problem (IVP) (1) at the point x n + 2 , assuming that w n = w ( x n ) . To start the procedure, consider the approximation in (2) determined at x n , and its first-order derivative determined at the points x n , x n + r , x n + 1 , x n + s , x n + 2 . By so doing, we obtain the following linear system of six equations in six real unknown parameters δ j , j = 0 , 1 , … 5 :

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 Q ( x ) in terms of w n , g n , g n + r , g n + 1 , g n + s , and g n + 2 . Putting these values in (2) while using the change of variable x = x n + t Δ x , we reach the following:

where

To obtain the two-step block method, we evaluate Q ( t ) at the collocation points x n + r , x n + 1 , x n + s , and x n + 2 , that is, we take t = r , 1 , s , 2 . This results in the following formulas that can be considered as a hybrid block method with two parameters:

where w n + i ≃ w ( x n + i Δ x ) are approximations of the true solution, and g n + i = g ( x n + i , w n + i ) , for i = r , 1 , s , 2 . Two unknown parameters r , s are related to the intra-step points x r , x s in the above-obtained approximations. We will equate the first leading term of the LTE of w n + 1 and w n + 2 to zero to achieve appropriate values for these parameters. As a result, optimal values for the parameters will be obtained, and at the end of the sub-interval [ x n , x n + 2 ] , the value w n + 2 will be the only value required to forward the integration to the next sub-interval. As a result, we consider the LTE of the formula given in (9) as:

At this stage, there are two parameters ( r and s ). We find the values of these two unknown parameters from leading terms of the LTEs as given in the aforementioned formulas for w n + 1 and w n + 2 . Therefrom, we obtain the following optimal values of the required parameters:

Substituting these values in (10) and (11), we obtain:

Thus, the two optimized parameters (intra-step points: r and s ) yielded the following two-step optimized hybrid block method:

It is worth to be noted that the above strategy concerning the optimization of the intra-step points was first introduced in ref. [20]. The pseudo-code 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:

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 A -stable algorithm, the theory of order stars is a relatively new concept in the field of stability analysis. One of the fundamental publications in ref. [23] should be read in order to obtain a better understanding of the theoretical aspects of order stars that arise from the structure of numerical techniques such as those presented here. Consider two complex-valued polynomial type functions denoted by W 1 and W 2 having degree i and j , respectively, with the quotient denoted by Ψ ( ρ ) = ρ 1 ρ 2 . As an example, consider the following stability function given in ref. [18] for the two-step optimized block method:

The A -stability characteristics have been depicted in Figure 1 for the method given in (13). The solution for the denominator ρ 2 in Ψ ( ρ ) is called a pole. Another assumption is to consider a complex-valued function Ω ( ρ ) that plays a crucial role in rest of the analysis. An order star Φ ( ρ ) divides the complex plane C into three distinct areas which are shown as the triplet written as { T + , T 0 , T − } . Members of the triplet can be defined in the following way:

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 R + is obvious and in the meantime contractivity region is shown by the triplet T − . The triplet R 0 evidently represents the boundary for the remaining two triplets in C . With the complex exponential function Ω ( ρ ) = exp ( ρ ) , our primary focus is upon the order stars of first kind. In this pursuit, the order star for the stability function Φ ( ρ ) given in (17) gives rise to the following three types of regions:

When the above discussed regions are plotted in C , the obvious star-like fingers became the reason for the shape to have name of order stars. An important concept related to the A -acceptability has to be understood before we proceed with the said regions.

Figure 1

Region of absolute stability of the two-step optimized hybrid block method given in (13).

The plot of order stars for the optimized block method is shown in Figure 2, where the shaded regions depict the triplet T + . It is worth noting that there is nothing common between T + and the imaginary axis. Additionally, no poles are seen in the left side of the complex plane, that is, where Re ( ρ ) < 0 . These essential characteristics show that the stability function of the optimized block method under investigation in the present research analysis has the properties of being A -acceptable with Re ( ρ ) < 0 .

Figure 2

The plot of order stars for the two-step optimized hybrid block method given in (13).

4 A variable stepsize formulation

With the help of an embedded-type 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 α and order β with β < α and performs their simultaneous execution. The method (13) is taken as a higher order method with α = 5 , and we choose a lower order method with β = 2 . The lower order method is used to obtain an approximation of the IVP (1) at the endpoint of the block x = x n + 1 , denoted by w n + 1 ∗ . This value is used to estimate the local error in w n + 1 ∗ in comparison to the more accurate approximation w n + 2 . The value w n + 2 is used for advancing the integration process, called local extrapolation. The lower order formula is constructed to share the FEs with the method to control the computational cost. To obtain a variable stepsize formulation, we used the following steps: Once the system in (13) has been solved, a new approximation of the solution is obtained with the lower order method at x = x n + 1 , which is named as w n + 1 ∗ . The corresponding local error l e n + 2 is

where w ( x ) shows the exact solution of the IVP in (1). For the solution w n + 2 provided by the higher order method, we have w ( x n + 2 Δ x ) − w n + 2 = O ( Δ x α + 1 ) , and thus the local error estimation, EST, is obtained as

This estimation corresponds to the local error of the lower order formula, since l e n + 2 is O ( Δ x β + 1 ) and hence dominates in (20) for enough small values of Δ x as explained in ref. [24]. In the present study, we consider the following lower order formula called the implicit trapezoidal rule [25] with second-order convergence:

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 second-order 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 ( Δ x old ) to the new one ( Δ x new ), using the following formula:

where β represents the order of the lower order method and 0 < κ < 1 is a safety factor whose purpose is to avoid failure steps. We assumed κ = 0.95 throughout the simulations for each method under consideration. In this process, the solver will keep trying until it obtains a stepsize for which the magnitude of EST is no more significant than the tolerance. Sometimes the solver might take much work to predict the new stepsize. To avoid this situation and prevent large fluctuations in the stepsizes, we impose some restrictions on the new stepsize ( Δ x new ), it is required to lie between Δ x min (minimum stepsize) and Δ x max (maximum stepsize), that is,

On the other hand, when ∣ ∣ EST ∣ ∣ < TOL , we can increase the step size for the next step, taking Δ x new = c Δ x old with c > 1 . Here, in the numerical experiments, we used Δ x new = 2 Δ x old . A suitable initial stepsize ( Δ x ini ) is required to start the process, for which different approaches can be used (see ref. [26]). One of the most commonly used approaches is the one given in ref. [27], which is about selecting a small starting stepsize, and later the algorithm will fine-tune it according to the changing step size approach defined by the user. The pseudo-code for the variable stepsize version of the two-step optimized block method (13) is provided in Algorithm 2 under Appendix A.

5 Numerical dynamics with results and discussion

This section describes how the two-step 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 n = 0 , the method (13) is solved as a system with w ( 0 ) = w 0 known to be the initial state of the system. After solving the problem at the initial state, we used the widely known second-order convergent Newton–Raphson method to obtain w ( x 1 ) = w 1 . The value w ( x 2 ) = w 2 is then calculated using w 1 from the preceding block as the initial value, and the procedure is repeated until the last value at x M is reached. We choose the duration of the integration interval as a multiple of 2 Δ x (i.e., x M − x 0 = k ( 2 Δ x ) , k ∈ N ) because the approach under consideration is two-step. The Newton–Raphson approach has been implemented using the Find Root command from Mathematica 12.1. It may be noted that all numerical computations are performed in Mathematica 12.1 on a personal computer running on Windows OS with Intel(R) Core(TM) i7-1065G7 CPU @ 1.30 GHz 1.50 GHz processor having 24.0 GB installed RAM. For comparison of the numerical simulations, we used three at least fifth-order methods, and one is the LobattoIIIA method of sixth-order. Five numerical experiments are presented, with two being scalar, while the other three are two- and four-dimensional systems. Considering the potentiality of block methods to deal with stiff differential models, we have taken some well-known stiff applied problems that have appeared several times in recent literature. These problems are solved with the following methods:

• Two-step at least fifth-order OBM (Block) with two intra-step points as given in (13).

• Two-step at least fifth-order OBM in reformulated version (RBlock) with two intra-step points.

• Two-step sixth-order OBM in RK form (RKBlock) with two intra-step points.

• One-step at least fifth-order OBM appeared in ref. [17].

• Laguerre polynomial hybrid block method (LPHBM) of at least fifth-order appeared in ref. [28].

• Fully-implicit RK-type fifth-order method called (Radau-I) appeared in [29, p. 226].

• Implicit RK-type sixth-order Lobatto method (Lob-III A) appeared in ref. [29, p. 228].

In Problem 1, a two-dimensional system of ODEs is simulated with both fixed and variable stepsize approaches of the two-step optimized block method (13). In Table 1, the tolerance ( ε ) is set to 1 0 − 3 while the initial stepsize ( Δ x ini ) is assumed to be 1 0 − 1 . The number of steps ( N ) and the number of FEs are identical in both approaches for a fair comparison. The numerical errors, including norm, root mean square, and the arithmetic mean, are considerably smaller in the variable stepsize approach than the approach taken with a fixed step size. This comparison confirms the usefulness of the variable stepsize approach for the block method. It may also be noted that the superiority of the latter approach is further confirmed in Table 2 when the tolerance is set to 1 0 − 6 . Moreover, since the variable stepsize approach is the major focus of this research investigation, numerical results are presented in the forthcoming simulations only with the variable stepsize approach for the two-step optimized block method given in (13), its reformulated version given in (14), and for its RK form as given by the Butcher tableau.

The numerical simulations are performed for Problem 2 under the variable stepsize approach setting the tolerance to 1 0 − 3 and 1 0 − 4 . Several methods, as mentioned before, are used to compute different types of errors and the number of FEs. For example, it is observed in Table 3 that the block method and its reformulated version used 65 evaluations to yield the errors with 1 0 − 9 , including OBM with a slightly more significant error. On the other hand, the RK version of the block method under consideration used 108 evaluations but reduced the amount of error to 1 0 − 10 . The rest of the methods, though performed well enough excluding Radau-I, took more evaluations to yield the errors with 1 0 − 9 and 1 0 − 10 . A similar sort of behavior is analyzed in Table 4 for Problem 2 under the variable stepsize approach when the tolerance was set to 1 0 − 4 .