A New Hybrid Block Method for Solving First-Order Differential System Models in Applied Sciences and Engineering

: This paper presents a new hybrid block method formulated in variable stepsize mode to solve some ﬁrst-order initial value problems of ODEs and time-dependent partial differential equations in applied sciences and engineering. The proposed method is implemented considering an adaptive stepsize strategy to maintain the estimated error in each step within a speciﬁed tolerance. In order to evaluate the performance and usefulness of the proposed technique in real-world applications, several differential problems from applied sciences and engineering, such as the SIR model, Jacobi elliptic function problem, and chemical reactions problems, are solved numerically. The results of numerical simulations in this work demonstrate that the proposed method is more efﬁcient than other existing numerical methods used for comparisons.


Introduction
Ordinary Differential Equations (ODEs) are essential for modeling real-life model problems in applied sciences and engineering.For instance, ODEs have found extensive application in modeling real-world problems in chemical, electrical, and mechanical engineering, physics, chemistry, and mathematical biology.Motivated by the diverse applications of ODEs in modeling physical phenomena, this paper introduces a new numerical method for integrating the following initial value problem (IVP): The above equation is defined over an integration interval [t 0 , t N ] and has initial conditions specified at t = t 0 .It is worth noting that the function f (t, v) in Equation ( 1) satisfies Lipchitz's condition, which guarantees the existence and uniqueness of the solution of the first-order IVP.
The first-order ODE given in Equation ( 1) has numerous applications in different engineering and applied sciences.One of the applications is chemical kinetics, where the rate of change of concentration of a chemical species in a reaction can be modeled using a first-order ODE.In this case, the function f (t, v) in Equation ( 1) represents the reaction rate, while the solution v(t) represents the concentration of the species at time t.Another application of first-order ODEs can be used to model the behavior of electrical circuits and systems in electrical engineering.In this case, the function f (t, v) in Equation (1) represents the current or voltage in the circuit, while the solution v(t) represents the state of the circuit or system at time t.In summary, a function f (t, v) in Equation (1) represents the rate of change, and the solution v(t) denotes the state of the system at time t.
A numerical approximation is important for understanding the behavior of differential systems when analytical solutions are known or unknown; numerical methods have been extensively applied to solving real-life model problems in many fields, including engineering and applied sciences.Numerical approximations play a critical role in developing and validating computational models by comparing approximate results with analytical solutions.Developing numerical techniques for solving differential systems is also a significant research activity since numerical methods provide valuable insights that enable scientists and engineers to make informed decisions and improve their understanding of real-world phenomena.For these reasons, much research is devoted to developing numerical techniques to solve Equation (1).
One of the numerical techniques for solving IVPs in ( 1) is the Runge-Kutta methods (RKMs).Explicit Runge-Kutta Methods (ERKMs) and Implicit Runge-Kutta Methods (IRKMs) are two main subclasses of RKMs.These methods form a wide range of numerical techniques that offer acceptable accuracy and efficiency in terms of both absolute and relative errors.ERKMs and IRKMs are one-step methods that use a Taylor series expansion to approximate the solution and are easy to implement and computationally efficient.However, explicit Runge-Kutta methods are frequently inappropriate for addressing stiff differential systems.Stiff problems are those where the solution varies rapidly over a short time.ERKMs can be unstable when applied to stiff problems, leading to inaccurate results.Therefore, researchers in numerical analysis often use IRKMs to integrate stiff differential problems.However, unlike ERKMs, IRKMs require solving nonlinear equations at each time step.This system of equations is usually solved using numerical methods such as Newton's method.IRKMs are more stable than ERKMs and can handle stiff problems more effectively.However, IRKMs are computationally more expensive than ERKMs and require more memory.In [1][2][3], the theory of RK-type methods, including the development, implementation, order condition, stability, convergence, and theoretical analysis of RK methods, are well discussed.
Multistep methods (MMs) are another numerical technique used to integrate (1).Unlike one-step methods, MMs utilize approximations of the solution at multiple mesh points to compute the subsequent approximate solutions, enhancing accuracy but demanding more function evaluations and computational time.Many experts in numerical solutions to differential equations have utilized various MMs to solve (1).The most common MMs for solving (1) are the linear multistep and the Adams-Bashforth methods.MMs can be used to solve stiff and nonstiff differential problems but require more memory and computational resources.Detailed discussions of MMs are reported in [4][5][6][7], including their development, implementation, stability, convergence, and theoretical analysis.
The main novelty of this paper is the introduction of a new hybrid block method (NHBM), which is capable of solving efficiently the IVP in (1).The NHBM is developed to overcome the limitations of the high number of function evaluations and stepsizes required by existing methods.The significance of the results presented in this paper is that the NHBM is an accurate and efficient approach for solving (1).The numerical experiments conducted in this study demonstrate that the NHBM outperforms other existing methods in terms of accuracy and efficiency.The NHBM also requires fewer total number of function evaluations and stepsizes, making it a more practical and computationally efficient method for solving (1).Overall, the improved accuracy and efficiency of the NHBM make it a valuable numerical method for researchers and practitioners in various fields, including engineering, physics, and applied mathematics.

Derivation of the Proposed Method
To approximate the theoretical solution of the problem defined in Equation ( 1) on any subinterval [t n , t n+1 ], we can use a polynomial of degree five denoted by P(t) and can be expressed as follows: where coefficients k j are unknown real numbers that satisfy collocation conditions at specific nodes.The v(t) derivative is obtained as follows: To develop the NHBM, we select three hybrid points b 1 , b 2 , and b 3 within the interval We then substitute t n into Equation (2) and t n , t n+b 1 , t n+b 2 , t n+b 3 , and t n+1 into Equation ( 3) to obtain the following system of equations: where the coefficients in Equation ( 2) are represented by k n , where n ranges from 0 to 5. Additionally, v n+j and f n+j represent the approximations of v(t n+j ) and v (t n+j ) = f (t n+j , y(t n+j )), respectively.
Once we have determined the values of k n for n = 0, 1, 2, 3, 4 and 5, we substitute them into Equation ( 2) and replace the variable t with t n + xh.This substitution allows us to obtain a continuous approximation using the polynomial in Equation ( 2), which is given by the following equation: where To obtain the proposed NHBM, we substitute x = 1, x = 1 5 , 1 2 , and x = 4 5 into Equation (4).After inserting b and simplifying, we obtain the following formulas: The above formulas provide the approximations for v(t n+ 1 ), v(t n+ 1 ), v(t n+ 4 5 ), and v(t n+1 ), respectively.

Analysis of the Proposed Method
To evaluate the accuracy of the proposed method, we need to determine its local truncation errors and order of accuracy.This is achieved by expressing the NHBM method given in (5) in a different form by rewriting it as: where, Ā and B are constant matrices that are obtained from the formulas in Equation ( 5), and This new form of the NHBM method is a recurrence relation that relates the values of V n and V n to the matrices Ā and B and the stepsize h.The matrix Ā is derived from the coefficients of the polynomial approximation in Equation ( 5), while the matrix B is obtained by multiplying the coefficients of the polynomial by the corresponding values of b 1 , b 2 , b 3 , and 1.This matrix form of the NHBM method is useful for numerical implementation and analysis, as it allows us to compute the values of V n and V n recursively using matrix operations.
We then proceed to define the operator L associated with the NHBM method given in Equation ( 5) as follows: where the operator definition in (7), involves computing a summation for all i belonging to set i = {0, b 1 , b 2 , b 3 , 1}.We assume that v(t n ) is sufficiently differentiable.To obtain the local truncation errors of the NHBM method, we expand the term in Equation ( 5) as a Taylor series about the point t n in powers of h.The following local truncation errors (lte) are obtained by comparing the Taylor series expansion of the exact solution with the Taylor series expansion of the numerical approximation: Due to the fact that our method is a one-step method, it exhibits zero stability; furthermore, the NHBM has an accuracy greater than one, which means it is consistent with the differential systems under consideration.As a result, the NHBM is convergent since it is both zero-stable and consistent.

Error Estimation of the Proposed Method
For optimal performance, variable stepsize methods are preferred.Here, we utilize an approach to estimate local errors at step endpoints.This involves employing a lowerorder formula to approximate v(t n+1 ) from the NHBM method, then comparing this to an approximation with a smaller stepsize.The difference helps estimate NHBM's error in each step, aiding stepsize adjustment for accuracy within tolerance.This balances larger steps for smooth solutions and smaller ones for rapid changes, enhancing accuracy while saving computational effort.In order to determine a reasonable estimate of the local error using the proposed NHBM, we use the following fourth-order method: with (lte) = − 17h 5 v (5) (t n )

4800
+ O h 6 .The following EREST formula estimates local error: In order to continue with our calculations, it is necessary to ensure that the absolute value of EREST is equal to or less than the specified absolute tolerances (AT).If this condition is met, we accept the results and proceed by increasing the stepsize.However, if EREST exceeds AT, we must reject the results and repeat the computations using a different strategy to improve accuracy.This computation involves adjusting the stepsize using the following formula: where AT denotes the absolute tolerance, and 0 < Λ < 1 is an adjustment factor that prevents stepsize failures during computation.This formula allows us to adjust the stepsize based on the magnitude of the error to achieve the desired accuracy level while minimizing computational effort.An additional p is also included in (9), which denotes the accuracy of the lower-order method presented in (8).

Implementation
Outlined steps for the proposed NHBM to integrate problem (1) using the strategy in Section 4 are as follows: Solution and Stopping: Solve resulting system of 4d nonlinear equations for d-dimensional IVPs.Utilize NM similarly for scalar IVPs.Implement a stopping criterion (| Ṽi+1 − Ṽi | ≤ 2 × AT) and set maximum iteration limit to 120 for NM execution.

Numerical Simulations
In this section, we demonstrate the numerical performance of the NHBM on several real-life models of the form (1). The formula used to measure the accuracy of the proposed method is given by: where v(t j ) represents the exact solution, while v j denotes the solution obtained by the NHBM at every time point represented by the discrete grid t j .This formula is used to compare the accuracy of the NHBM with existing numerical methods.Throughout this paper, different methods and performance measures are represented using the following abbreviations: • NHBM denotes the new hybrid block method defined in (5).• OOBM denotes one-step optimized block method in [23].• ode15s is the numerical solver available in MATLAB version R2021b.• RADAUIIA denotes an implicit Radau method of the fifth order in [1].• IGRKM represents the implicit Gauss-Legendre method of the sixth order in [4].• THBM is the two-step hybrid block method in [20].• FE denotes the total number of function evaluations.

Example 1
We first consider the SIR model problem from Kashkari and Syam [23], which arises in mathematical biology.This SIR model describes infectious disease dynamics over time.It involves the following three connected differential equations representing susceptible S(t), infected I(t), and recovered R(t) individuals: where ∆, θ, and α are positive parameters.Additionally, we introduce the function v(t), which is defined as the sum of the susceptible, infected, and recovered individuals at time t.This function is given by: By adding Equations ( 10)-( 12), we obtain the following expression for the derivative of v(t): For a particular closed population, assuming an initial condition of v(0) = 0.5 and ∆ = 0.5, we can rewrite the above equation as: The results presented in Table 1 provide evidence of the effectiveness and better performance of the proposed method.To further demonstrate the accuracy of the method, we include graphical representations in Figures 1 and 2 of the absolute errors of the NHBM, as well as a comparison of the exact and discrete solutions of the NHBM using AT = 10 −8 , h 0 = 10 −2 and t ∈ [0, 60].The plots clearly indicate that the NHBM produces highly accurate numerical solutions that closely match the exact solutions.These results highlight the effectiveness of the proposed method in accurately approximating solutions to the SIR model in mathematical biology.

Example 2
Furthermore, we solve Robertson's stiff problem, which is a well-known differential system that arises in chemical sciences.This system describes the kinetics of an autocatalytic reaction and is known for its stiffness, making it challenging to solve numerically.This problem involves a set of three nonlinear differential equations, which are given by the following equations: where, v 1 (t), v 2 (t), and v 3 (t) represent the concentrations of three chemical species at time t.The system is subject to initial conditions v 1 (0) = 1, v 2 (0) = 0, and v 3 (0) = 0, and is defined over the time interval 0 ≤ t ≤ 40.To evaluate the accuracy of our method, we compare the numerical results with a reference solution obtained from [24].Specifically, we use the maximum absolute errors to compare the obtained numerical results with the reference solution at the endpoint t N = 40.The reference solution is given by: Using the above reference solution, we determine the accuracy of our method and gain insights into the behavior of Robertson's stiff problem.
In Example 2, we compare the accuracy of the proposed NHBM method with two other commonly used methods, namely ode15s and RADAUIIA.The results are presented in Table 2, which clearly shows that the NHBM method outperforms the other two methods in terms of accuracy.By providing a more accurate and efficient approach for stiff differential systems, our study contributes to the development of new numerical methods for solving challenging problems in chemistry and chemical engineering.

Example 3
Let us consider a Jacobi elliptic function problem that arises in physics.The problem is defined by the following system of differential equations: The exact solution to this problem is given in terms of Jacobi elliptic functions SN, CN, DN, which can be expressed by series expansions.Specifically, the exact solution is given by: where SN, CN, DN are Jacobi elliptic functions.
In Example 3, we compared the proposed NHBM method with two existing numerical methods, namely ode15s and RADAUIIA.Our comparison results are presented in Table 3, which clearly demonstrates that the NHBM method outperforms the other two methods in terms of FE and MAE.Our study contributes to developing a new efficient numerical method for solving challenging problems in physics by providing a more accurate and efficient approach for Jacobi elliptic differential systems.For our last example, we consider a singularly perturbed parabolic time-dependent partial differential equation that arises in chemical engineering.This problem is a reaction-diffusion process in a catalytic pellet, where the concentration of the reactant is modeled by v(x, t), x ∈ [0, 1] and t ∈ [0, T].The reaction rate is assumed to be proportional to the concentration of the reactant, and the rate of diffusion is assumed to be small.The problem is given as follows: where k is a positive constant (reaction rate) and is a small positive constant (diffusion rate).Boundary conditions are at x = 0 and x = 1, with an initial condition at t = 0.However, the exact solution is unknown.
To solve Example 4 using the NHBM, the second-order spatial derivative is discretized using the method of lines.The following centered difference approximation is used to discretize v xx : where , and N represents the internal number of spatial nodes.
To evaluate the accuracy of our proposed method, we compare NHBM with the IGRKM method.We used NDSolve, a numerical method for solving differential equations implemented in Mathematica, as the reference solution to determine the MAE.Our comparison results are presented in Table 4, which clearly demonstrates that our proposed NHBM method outperforms the IGRKM method in terms of accuracy.To further demonstrate the effectiveness of our approach, we present surface plots of the proposed NHBM and NDSolve solutions for Example 4 in Figure 3.The results in Table 4 and Figure 3 distinctly showcase the proposed method's efficacy in accurately approximating solutions to the reaction-diffusion process that arises in the chemical engineering domain.

Concluding Remarks
This paper proposes a new numerical method that utilizes a variable stepsize formulation to solve differential systems in (1) efficiently.We also discuss a technique to control the stepsize of the method to enhance solution accuracy.The efficacy of the proposed approach is validated through its application to four distinct model differential problems within applied sciences and engineering, and the results are compared to existing solutions.The variable stepsize implementation used in this paper produces better performance, as evidenced by the computed results.We conclude that the NHBM can be a valuable numerical method for solving first-order differential model problems in various fields.

Table 2 .
Comparison of numerical solutions for Example 2.

Table 3 .
Comparison of numerical solutions for Example 3.