Implicit Four-Point Hybrid Block Integrator for the Simulations of Stiff Models

Over the years, the systematic search for stiff model solvers that are near-optimal has attracted the attention of many researchers. An attempt has been made in this research to formulate an implicit Four-Point Hybrid Block Integrator (FPHBI) for the simulations of some renowned rigid stiff models. The integrator is formulated by using the Lagrange polynomial as basis function. The properties of the integrator which include order, consistency, and convergence were analyzed. Further analysis showed that the proposed integrator has an A-stability region. The A-stability nature of the integrator makes it more robust and fitted for the simulation of stiff models. To test the computational reliability of the new integrator, few well-known technical stiff models such as the pharmacokinetics, Robertson and Van der Pol models were solved. The results generated were then compared with those of some existing methods including the MATLAB solid solvent, ode 15s. From the results generated, the new implicit FPHBI performed better than the ones with which we compared our results with. DOI:10.46481/jnsps.2022.777


Introduction
In this research article, an A-stable implicit FPHBI is formulated for the simulations of stiff models. The concept of Astability was introduced by [1] to address the effect of stiffness in differential equations. The desire for this special property (A-stability) implies that suitable methods have to be derived for the simulations of stiff models, [2]. This has, therefore, motivated us to formulate an implicit FPHBI for the simulations of stiff models of the form, where y T = (y 1 , y 2 , ..., y m ), µ T = (µ 1 , µ 2 , ..., µ m ), A is an m × m matrix and ϕ is an m-dimensional vector. Equation (1) is assumed to satisfy the Lipschitz conditions, [3]. Equations of the form (1) often arise in many applications in engineering and sciences. Suffice to say that most of these equations often results to stiff differential equations which in some cases do not have closed form solutions. Stiff models are found in description of pharmacokinetics, chemical reactions, molecular dynamics, control systems, mechanics, electronic circuits, lasers, [4][5][6][7][8].
In general, equation (1) may be stiff or non-stiff. A nonstiff differential equation is one in which there is a simultaneous evolution of solution components and also has time-scales that are comparable. They are often solved using explicit methods with some error control [9]. On the contrary, stiff differential equations do not have a universally accepted definition. The term 'stiff' first appeared in [10]. Equation (1) is called stiff if it satisfies any of the following definitions: 1. the stability step-size is much more smaller than the accuracy step-size [11]. That is, the step-size depends on stability requirement rather than the accuracy requirement [12], 2. the explicit method performs slowly or do not even work on it [13], 3. it has time scales that vary widely. In order words, some solution components decay more rapidly than others [12], 4. all its eigenvalues have negative real parts with large stiffness ratio [12], 5. Real(λ i ) < 0, i ∈ [1, m], where λ i is the eigenvalues of the matrix J = ∂ f /∂y called the Jacobian matrix of equation (1)  |Real(λ i )| is often called the stiffness ratio. The stiffness ratio measures the degree of stiffness of the system [14].
Conditions (5) and (6) as given by [14] shall be adopted as the definition of the stiff models that shall be considered in this research. The expression, is the general solution of the stiff model (1) where k i are the eigenvectors corresponding to the eigenvalues λ i , ∈ i are arbitrary constants and σ(x) is a particular integral. Taking x as time, the first term of y(x) = m i=0 ∈ i e λ i x k i is denoted as transient solution and the second term σ(x) as steady-state solution.
Assuming the stiff model (1) satisfies condition (5), then the term y(x) = m i=0 ∈ i e λ i x k i → 0 as x → ∞. Suppose |Real(λ u )| and |Real(λ v )| further satisfies the condition so that fastest and slowest transients are ∈ i e λ u x k i and ∈ i e λ v x k i respectively. If the stiff model (1) is solved numerically and aimed at achieving a steady-state phase, continuing integration is needed until the slowest transient is negligible [2]. Thus, the smaller |Real(λ v )| is, the longer the integration time. On the other hand, for the larger|Real(λ u )|, a sufficiently small step is required so that h will lie within the stability region of the method [12,15].
Many authors have made several attempts at solving stiff models. For instance, the authors in [16] proposed k-step hybrid methods for solving stiff models arising from chemical reactions. The authors also proved some relevant theorems regarding the regions of stability of the methods. The authors in Figure 1: Physical illustration of the implicit FPHBI [2] developed a block backward differentiation formula that is diagonally implicit for the approximation of stiff pharmacokinetics models. They further carried out convergence and stability analysis of the formula. Also, the authors in [17] derived an optimized hybrid block Adams method for the solution of first order differential equations. The A-stable method which is of order six was found to be convergence, zero-stable and consistent. The following authors also developed hybrid methods for solving different types of stiff differential equations [2,8,9,16,[18][19][20][21][22][23][24][25][26][27]. Furthermore, these authors also proposed different approaches to solving some special differential equations [28][29][30][31][32][33][34][35].
The predictor formulae that compute the previous block at the points x n−1 and x n is derived in similar fashion like the implicit FPHBI. It is given by the formulae,

Implementation
The newly formulated implicit FPHBI in equation (9) will be implemented in a predictor-corrector mode. Firstly, we set the data inputs namely, the problem to be solved, initial conditions, step size and tolerance level. The predictor in equation (10) is used to evaluate the previous block at the points x n−1 and x n . The implicit FPHBI which serves as the corrector and main method is then employed in solving the required stiff models. The codes for the implementation of the method were written in MATLAB 2021a while the method was derived using Scientific Workplace 5.5.

Analysis of the Implicit FPHBI
The analysis of the new implicit FPHBI shall be carried out in this section.

Order Definition 3.1 [36]
A method and its associated difference operator L given by Applying equation (12) on the implicit FPHBI (9), we obtain  Therefore, the implicit FPHBI is of uniform order 8 with the error constant Since the new implicit FPHBI is of order 8, it implies that it is consistent. 290 3.3. Zero-Stability Definition 3.3 [12] If no root of the characteristic polynomial has a modulus greater than one and every root with modulus one is simple, then such a method is called zero-stable. [1] proposed scalar test to ascertain the zero-stability status of a method. Thus substituting the scalar test problem into the implicit FPHBI equation (9), we obtain, Equation (14) is then written compactly in matrix form as, Equation (15) can be written as where Substituting H = 0 into (17) gives, Therefore, the zeros of equation (18) are t 1 = t 2 = t 3 = t 4 = t 5 = 0 and t 6 = 1. Since all the zeros lie within |t| ≤ 1, the implicit FPHBI is said to be zero-stable.

Convergence
Theorem 3.1 [12] Consistence and zero-stability are the necessary and sufficient conditions a method must satisfy in order to be convergent. Thus, the FPHBI is convergent since it satisfies the conditions for consistency and zero-stability.

Stability Regions 3.5.1. Stability Region of the Implicit FPHBI
The part of the complex plane where a method is absolutely stable when applied to the scalar test equation y = λy is called its stability region. Definition 3.4 [14] If the stability region of a method contains the whole left half-plane Re(hλ) < 0, such a method is referred to as A-stable.
The graphical plot of the implicit FPHBI is presented in Figure 2.

Comparison of Stability Regions and Intervals of Instability
The following figure shows the stability region of the Diagonally Implicit Block Backward Differentiation Formula (DIBBDF) derived by [2]. 292 Figure 3: Stability region of DIBBDF by [2]  The stability of the new implicit FPHBI (see Figure 2) shall be compared with that of the DIBBDF (in Figure 3). While both methods are A-stable, it is obvious from the two figures that the implicit FPHBI has a larger stability region than that of DIBBDF. Thus, the FPHBI is expected to be more accurate and efficient than the DIBBDF, since the larger the stability region of a method, the better the performance. Note that the stable region is the region that lies outside the closed curve. The interval of instability of the new implicit FPHBI is compared with that of DIBBDF as shown in Table 1. From Table 1, it is clear that the implicit FPHBI has a larger stability region than the DIBBDF.

Test Problems
The implicit FPHBI derived in equation (9) shall be used to simulate some rigid well-known stiff models in order to test its accuracy and efficiency. Numerical solution, maximum error, computation time, efficiency curves and solution curves of the said problems shall be presented.

Problem 1 (Pharmacokinetics Model)
Pharmacokinetics is a phenomenon that describes how an administered drug behaves in the body over a period of time. Consider the pharmacokinetics model (see Figure 4) composed of two components; that is, the drug concentration in the gastrointestinal (GI) tract y 1 (t)and drug concentration in the bloodstream y 2 (t) as a function of time t given by, y 1 = −k 1 y 1 + d, y 1 (0) = 1 y 2 = k 1 y 1 − k e y 2 , y 2 (0) = 0 (19) where k e (> 0) and k 1 are the clearance rate constant and rate constants from one compartment to another respectively. The parameter d is the regimen of drug intake. Taking k 1 = 2(ln 2)and k e = (ln 2)/5, the exact solution of equation (19) for t ∈ [0, 6] are given by Note that y 1 (x) = e −k 1 t is the drug's exponential decay in terms of absorption [37]. This implies that the behaviour in the long term of the concentration of the drug in the GI tract will reduce to level zero. Equation (18) is a two-compartmental pharmacokinetics model formulated by [38] to model the drug flow via the compartments of the body namely the GI tract and the circulatory system. The drug moves from the GI tract compartment into the bloodstream compartment at a rate relative to the drug's concentration in the GI tract [2]. The drug is finally metabolized and cleared from the bloodstream at a rate relative to its concentration there [39]. For more details on general two-compartment pharmacokinetics models, see [37,40,41].

Problem 2 (Robertson Model)
According to [4], the Robertson model defines an autocatalytic kinetics reaction. The Robertson problem, popularly called the ROBBER according to [13] is a set of three ordinary differential equations that can be put up under some ideal conditions. It is mathematically given by the stiff system with (y 1 (0), y 2 (0), y 3 (0)) T = (y 01 , y 02 , y 03 ) T , where y 1 , y 2 and y 3 are the concentrations of A, B and C respectively and y 01 , y 02 and y 03 are concentrations over a period of time t = 0. This problem was proposed by [42] and the file rober.f containing the software of the problem can be found in [43]. The Robertson problem is composed of the following reactions where A, B and Care chemical species and k 1 , k 2 and k 3 are rate constants. It is important to state that the cause of stiffness of this problem is the substantial variation in the reaction rate constants.

Problem 3 (Van der Pol Model)
The Van der Pol equation proposed originally in 1920s by the Dutch electrical engineer Balthazar Van der Pol (1889-1959 has been used to model systems that have self-excited oscillations of limit cycles. The equation depicts oscillatory processes in physics, electronics, biology, neurology, etc. [44]. Van der Pol himself used the equation to model the range of stability of the human heart dynamics. Furthermore, coupled neurons in gastric mill circuits of lobsters in the field of biology were successfully modelled using the Van der Pol equations [45,46]. The equation has also been used as a model in seismology for developing a model in a geological fault that shows interaction of two plates [47]. The generation of spike in giant squid axons can also be modelled using Van der Pol equation [48]. It is therefore important to explore ways of developing a deep understanding of the Van der Pol equation in view of its widespread applications. One of the ways to do this is by modeling and simulating it. The Van der Pol equation is a stiff nonlinear equation expressed as, Equation (23) is transformed to its equivalent system of first order system of the form, y 1 = y 2 , y 1 (0) = 2 y 2 = −y 1 + (1 − y 2 1 )y 2 ∈, y 2 (0) = 0 This transformation is achieved by substituting y 1 = ϕ(x), y 2 = µϕ (x) and x = x/µ, where ∈= 1 µ 2 is a parameter that controls stiffness. In this paper, the value of ∈= 500.

Results and Discussions
In this section, the implicit FPHBI is employed in simulating stiff models presented in Section 4. This is aimed at testing the accuracy, efficiency and computational reliability of the proposed integrator.
The following abbreviations shall be used in the Tables 2-4.
x: Point of evaluation h: Step size T /s: Computation/Execution time in seconds ABSE: Absolute error MAXE: Maximum error NUMSOL: Numerical (approximate) solution DIBBDF : Diagonally implicit block backward differentiation formula developed by [2] KSHM: k-step hybrid method derived by [16] Ode 15s: MATLAB inbuilt stiff solver  where NS is the total number of steps, y(x)is the exact solution and y n (x) is the computed (approximate) solution.
In Table 2, it is obvious that the FPHBI performed better than the DIBBDF developed by [2] and the ode 15s solver as the new method has lesser maximum error for Problem 2. It is also important to state that the FPHBI is more efficient in terms of computation time than those of DIBBDF and ode 15s. This is evident in Table 2 and the efficiency curves presented in Figure  5. Table 3 presents the approximate solutions of Problem 2 at points x = 0.4,x = 40 and x = 4000for the three solution components y 1 (x), y 2 (x) and y 3 (x). The results obtained clearly show that the implicit FPHBI effectively approximates Problem 2. To further buttress this point, accuracy (solution) curves were plotted for the problem, see Figure 6. Table 4 presents the numerical solution of the Van der Pol model in (24). Since the problem does not have exact solution, the approximate solution using the newly derived implicit FPHBI is compared with that of MATLAB in built stiff solver (ode 15s) at the end points x = 1, x = 5, x = 10 and x = 20. From the solution curve obtained, it is clear that the result of the FPHBI converges to that of the ode 15s solver, see Figure 7. Thus, the FPHBI is computationally reliable. 294

Conclusion
In this paper, an implicit FPHBI was derived for the simulations of first order stiff models. Special rigid stiff problems like the pharmacokinetics, Robertson and the Van der Pol models were considered. The results obtained clearly showed that the new integrator is computational reliable, as it performed creditably well. The paper further analysed the integrator on the basis on consistency, zero-stability and convergence. The stability region of the integrator was plotted and the plot generated shows that the integrator is A-stable, thus making it fit for simulating stiff models. Finally, the authors are optimistic that this study has added to the collective understanding of the nature of solutions of these stiff models.