Optimization of the richardson integration over fluctuations of its step sizes

In this paper, we examine the optimization of Richardson numerical integration of an arbitrary real valued function in the space of step sizes. Namely, as one of the most efficient numerical integrations of an integrable function, the Richardson method is optimized under the variations of its step sizes. Subsequently, we classify the stability domains of the Richardson integration of real valued functions. We discuss stability criteria of the Richardson integration via the sign of the fluctuation discriminant as a quintic or lower degree polynomials as a function of the step size parameter. As special cases, our proposal optimizes the trapezoidal, Romberg and other numerical integrations. Hereby, we consider the optimization of the Richardson schemes as a weighted estimation in the light of extrapolation techniques. Finally, optimal Richardson integrations are discussed towards prospective theoretical and experimental applications and their industrial counterparts. Subjects: Science; Mathematics & Statistics; Applied Mathematics; Technology; Engineering & Technology; Mathematics & Statistics forEngineers; Engineering Mathematics; Design Bhupendra Nath Tiwari ABOUT THE AUTHOR Dr. Bhupendra Nath Tiwari an alumnus of Jawaharlal Nehru University New Delhi, and Indian Institute of Technology Kanpur, India is an Assistant Professor at the University of Information Science and Technology, “St. Paul the Apostle”, Ohrid, Macedonia; act as Ricercatore, Visiting Scientist, PDF at INFNLaboratori Nazionali di Frascati, Rome, Italy. He has published research papers, reviews, books and monographs (in English & French). Having qualified UGC-FRP, CSIR-JRF/NET/SRF, GATE, GRE, TOEFL and CILS, he is a professional body member of the AMS, IOP, APS, MDPI, IEEE, Springer, Elsevier, World Scientific Publications and Hindawi Publishing Corporation, Worldwide. His research interests include theoretical and mathematical aspects of S&T. Er. Amarasingha Arachchige Mihiri Chathurika is a master candidate at the University of Information Science and Technology, “St. Paul the Apostle”, Ohrid, Macedonia. She is Machine Learning trainee at Lightblue Technology, Tokyo, Japan who aspires to work as a Machine Learning Engineer. PUBLIC INTEREST STATEMENT This paper examines the optimization of the Richardson integration of an arbitrary real valued function under variations of its step sizes. Richardson integration is one of the most efficient methods to integrate an integrable function numerically. Our method offers minimized value of the errors under variations of step sizes of the integration. We classify the stability domains of the Richardson integration of arbitrary real valued functions. Subsequently, we discuss stability criteria of the Richardson integration through the sign of the fluctuation discriminant that arises as a quintic polynomial or lower degree polynomials as a function of the step size parameters. As special cases, our proposal optimizes the trapezoidal, Romberg and other numerical integrations. Hereby, we consider the optimization of Richardson schemes as a weighted estimation in the light of extrapolation techniques. Finally, optimal Richardson integrations are discussed towards the prospective theoretical and experimental applications and their industrial counterparts. Tiwari & Chathurika, Cogent Mathematics & Statistics (2019), 6: 1643438 https://doi.org/10.1080/25742558.2019.1643438 © 2019 The Author(s). This open access article is distributed under a Creative Commons Attribution (CC-BY) 4.0 license. Received: 14 January 2019 Accepted: 24 June 2019 First Published :17 July 2019 Corresponding author: Bhupendra Nath Tiwari, INFN-Laboratori Nazionali di Frascati, Via E. Fermi 40, Frascati 00044, Rome, Italy E-mail: bhupendray2.tiwari.phd@iitkalumni.org Reviewing editor: Lutz Angermann, Mathematics, Clausthal University of Technology, Clausthal-Zellerfeld Germany Additional information is available at the end of the article


PUBLIC INTEREST STATEMENT
This paper examines the optimization of the Richardson integration of an arbitrary real valued function under variations of its step sizes. Richardson integration is one of the most efficient methods to integrate an integrable function numerically. Our method offers minimized value of the errors under variations of step sizes of the integration. We classify the stability domains of the Richardson integration of arbitrary real valued functions. Subsequently, we discuss stability criteria of the Richardson integration through the sign of the fluctuation discriminant that arises as a quintic polynomial or lower degree polynomials as a function of the step size parameters. As special cases, our proposal optimizes the trapezoidal, Romberg and other numerical integrations. Hereby, we consider the optimization of Richardson schemes as a weighted estimation in the light of extrapolation techniques. Finally, optimal Richardson integrations are discussed towards the prospective theoretical and experimental applications and their industrial counterparts. Keywords: Optimization theory; fluctuation theory; stability analysis; richardson integration; numerical techniques

Introduction
Optimization theory enables an understanding of fluctuation phenomena and their stability properties for a wide class of objective functions and searching algorithms (Parikh & Boyd, 2014). In this paper, we concentrate on optimization of numerical integration techniques over variations of their step sizes, see (Wen, Wang, Jin, Wong, & Ting, 2016) towards the Bayes-optimal joint channel-anddata estimation. Following the fact that the Richardson method incorporates most of the known numerical integration schemes (Chapra & Canale, 1998), we optimize the Richardson estimation of a definite integral of an arbitrary integrable function under variation of its step sizes. In the light of scientific computing (Kincaid, Kincaid, & Cheney, 2009), extrapolation techniques play an important role in trading-off the speed and accuracy of a given numerical integration or a proximal algorithm. In this regard, there are interesting developments towards the fractional calculus (Changpin & Zeng, 2015), critical exponents and series extrapolation techniques for a sequence of data (Cöster & Schmidt, 2016), and associated numerical extrapolation samples (Liu, 2006).
Historically, the Richardson method was found by Lewis Fry Richardson in 1927 who called it as the deferred approach to the limit (Schmidt & Henrici, 1966). This concerned a discrete variable analysis for a system of ordinary differential equations. Further, the Richardson method finds applications in a number of areas of numerical methods, statistical analysis and applied mathematics, see (Schmidt & Henrici, 1966) for analyzing a system of ordinary differential equations, (Little & Rubin, 2014) statistics of missing data, (Hundsdorfer & Verwer, 2013) numerical solution of diffusion equations and (Lambert, 1991) differential systems as initial value problems. As the result, it has become a powerful computational tool (Ouannas & Odibat, 2015), which is massively used in discrete dynamical systems and their evolutions as the initial/boundary value problems. For a system of ordinary differential equations, the Richardson method gives the best approximate solutions with an improved accuracy of the results, see (Ruyer-Quil & Manneville, 2002) for an instance. First of all, for a given real valued function, an extrapolation of the numerical integrations to the limiting zero interval spacing is realized on the basis of algebraic methods (Kappe & Warren, 1989), and associated numerical evaluations and their error analysis (Igoe, 1967). Hereby, it is demonstrated that some numerical evolutions confirm an essential validity towards their optimization procedure (Deb, 2005). In contrast to the above, some of the numerical schemes however show that their extrapolation technique must be used with a caution. As a matter of the fact, the implementation can vary into two types as the active and passive Richardson extrapolations, for an instance, see (Faragó, Havasi, & Zlatev, 2010) towards the stability analysis, optimization algorithms and reliability, (Mona, Lagzi, & Havasi, 2015) towards advection problems, and (Dimov, Zlatev, & Faragó, 2017) for the associated engineering applications and their practical aspects.
In order to reduce the step size, many methods can be used where the higher order terms are considered as the floating points while rounding-off of the errors in a given numerical integration scheme. In this concern, Refs. (Maharani, Salhi, & Suharto, 2018) provide an optimization analysis of floating locations and robust domain determination in the light of the linear system theory and artificial intelligence based algorithm. As an extrapolation model, such problems can be solved by using a weighted average (Liu, 2006), also see (Gao, Sun, & Zhang, 2014) for fractional numerical differentiations as an approximation of the Caputo fractional derivative (Podlubny, 1998) and its applications towards computational studies. Namely, a finite sequence fF n g can be realized as the values of the function FðhÞ. As h ! 0 þ , such an extrapolation model can be expressed (Liu, 2006) as the following asymptotic series where the first term F is understood as the limit F ¼ lim h!0 FðhÞ; α k are constants that are independent of the step size h, and the exponents fp k g satisfy p 1 < p 2 < Á Á Á < p sþ1 .
In this concern, the Richardson extrapolation arises as a type of numerical differentiations that includes highly accurate finite difference formulae and reduces various errors (Chapra & Canale, 1998;Keller, 2018;Liu, 2006). In order to get an accurate approximation, we offer the optimal Richardson method that combines various numerical integral estimates such as the trapezium, Simpsons and Romberg techniques (Chapra & Canale, 1998). As a special case of the Richardson integration scheme, the Romberg method arises as an initial value theory (Bauer, Rutishauser, & Stiefel, 1963), whereupon various numerical and computational aspects are established in the light of numerical integrations of an integrable real valued function (Samoradnitsky, 2017). In early time analysis, such evolutions correspond to stable non-Gaussian random processes (Lambers & Sumner, 2016).
Fluctuation theory finds its relations to the stochastic modelling, dynamical systems, recurring events, control systems and queuing theory (Asmussen & Albrecher, 2010). In this regard, there exists a concrete formula (Feller, 2015) for stable distributions and its applications to Markov chains and their ergodic properties. As per this introduction, we consider fluctuations of the Richardson numerical integration of an arbitrary real valued function in the space of step sizes. In practice, fluctuations introduce some fuzziness in the parameters, whereby they are measured probabilistically, see (Zimmermann, 1975) towards the optimization of fuzzy systems. Fluctuation theory has an intrinsic geometric origin, see (Ruppeiner, 1995) for a Riemannian geometric description of thermodynamic fluctuations, and (Tiwari, 2008) for its relation to black objects and their geometric invariants at different orders of perturbations. Indeed, fluctuations are closely intertwined with optimization theory, see for instance (Xie, Zhang, & Yang, 2002) towards the particle swarm optimization, evolutionary computing and dissipative models and (Harkey & Kenny, 2000) for conductance fluctuations (termed as 1=f noise) in microelectromechanical systems. System optimizations have their deep-rooted origin in the theory of estimation, for example, the second-order homogenization estimates (Castan˜eda, 2002) as well as in the theory of computing and information science, see (Saitou, Izui, Nishiwaki, & Papalambros, 2005) towards product developments and structural optimization techniques. Furthermore, it is worth mentioning that the fluctuation theory analysis plays an important role in examining the Brownian motion of a particle moving in a periodic potential (Bao-Quan et al., 2003) and associated numerical integral, also see (Hairer, Lubich, & Wanner, 2006) towards its geometric description, structure-preserving algorithms and numerical solutions of a system of ordinary differential equations.
In the light of embedded systems, our model does not stop here, but it finds its perspectives in the realm of computational methods, Monte Carlo simulations, random processors, viz. Gaussian and non-Gaussian configurations (Hairer et al., 2006). With this motivation, we focus on the stability structures of an arbitrary embedded system and its flow analysis, see (Tiwari, 2011a) in the light of convexity theory and its geometric perspectives. Perspectives of fluctuation theory include dynamic optical optimization , multicomponent constrained optimizations , power networks and their geometric stability analysis (Bellucci, Tiwari, & Gupta, 2012). Further engineering developments exist, see (Gupta, Tiwari, & Bellucci, 2013) towards an intrinsic geometric understating of the network reliability and voltage stability towards the optimal solution of power system planning and its operation. Such problems are highly nonlinear and large-scale in their nature. Notice that fluctuation theorems arise as a consequence of the time reversal symmetry, and symmetries of the hyperbolic spaces and transitive maps. Namely, such a theorem provides asymptotic observations in time. In the framework of approximation theory, a quantitative and parameter free relation between the stationary probabilities of observing a value emerges as the average entropy production rate and vice-versa. Further interests include its connections towards the chaotic hypothesis and chaotic dynamical systems (Ott, 2002).
Consequently, as per the above notion of the system reliability, optimizations of such systems as a state model yield different thermodynamic properties (Kozioł & Wiśniewski, 2008), see (Angelis, 1994) also towards a variation theory investigation and (Lei, Chengjun, & Chen, 2018) for associated flow properties in the space of parameters. In this paper, as far as the stability of the Richardson scheme is concerned, we find that the fluctuation discriminant arises as a quintic equation in the step size parameter, whereby we discuss the statistical nature of the undermining ensemble both in general and various limiting cases. This offers an intriguing relationship with fundamental objects in the algebraic theory of equations [47]. For example, finite degree rational curves show an interesting property in relation to compact Calabi-Yau three-folds (Candelas, Xenia, Green, & Parkes, 1991) and their moduli spaces that can be viewed as a discrete finite set, that is, a rational curve of its highest degree nine (Johnsen & Kleiman, 1996) or that of a curve of it degree eleven (Cotterill, 2012) on a generic quintic 3-fold. In this regard, the associated Donaldson-Thomas invariants render interesting platforms to study algebraic equations and rational curves in the light of the algebraic geometry (Clemens, 1984), Picard-Fuchs equations and mirror maps (Cox & Katz, 1999).
Dynamics of nonlinear configurations possess non-integrable and asymptotic characteristics. In general, such systems cannot be solved by using an analytic method. Therefore, numerical perturbation and reduction techniques are used to optimize the evolution of such systems as a function of the time. For example, the Hamilton-Jacobi-Bellman equation is one approach that optimizes such nonlinear systems (Clarkson, 1993). This offers an intrinsic interrelation of the optimization theory to nonlinear dynamical systems and associated geometric methods. Furthermore, as a nonlinear system, nonlinear adaptive controls enable us to identify the unknown parameters of the model (Marino, 1997). As a stabilization problem, the uncertainties can enter via either the state or input variables of a given system as a pair of dual variables, such as the coordinate and momentum or the energy and time, see (Tiwari, 2011b) towards generalized uncertainty principles and their relation to physical systems. In the light of the control theory, a complete set of solutions are known for certain limiting quadratic systems (Khargonekar, Petersen, & Zhou, 1990). Such robust stabilization problems exist in the realm of the control theory (Khargonekar et al., 1990) and stochastic models (Kushner, 1967).
Notice that considerations concerning the Richardson integration that focus on the optimization of coefficients in Romberg-like algorithms are of a fundamental importance, see for instance the classic book by Store and Bulirsch (Stoer & Bulirsch, 2013) for an overview of the literature. In this concern, there exist a number of computational analysis (Kiusalaas, 2013) in the light of the Romberg integration scheme and extrapolation methods for a smooth nonperiodic function. However, as far as the convergence of the solution f of an ordinary differential equation is concerned, it is worth mentioning that the existing Romberg quadratures do not provide an improvement (Ciarlet, 2001), when the derivative of f is discontinuous at certain point(s) in the solution space. The optimization of such solutions and their application arise in diverse engineering and system designing problems. Following the same, we have provided an intrinsic classification scheme and computed stability properties of the Richardson integration of an arbitrary real valued smooth function in terms of algebraic relations, namely, as the polynomials in the step-size parameter. As an innovation in the optimization theory, our proposal offers stable solutions of ordinary differential equations that are function independent in their nature. Hereby, the above proposed scheme finds its importance in engineering applications and interdisciplinary research. With this motivation, we explore the optimization of an arbitrary Richardson integral as a weighted extrapolation model, see (Cabrera & Parks, 1991) towards the understanding of spectral estimation techniques, data analysis and signal processing.
In this paper, we give the optimal Richardson integration of a real valued function under variations of its step sizes. Namely, we have provided the optimization of the Richardson numerical integration in the space of its step sizes. Herewith, about its critical points, we find that the flow components defined as the first order partial derivative of the integration have the same ratio as that of the step sizes. For a given real smooth function, we have shown that the critical points of the Richardson numerical integration are expressible as a function of the flow components. In this setting, the critical points of the Richardson integration arise as the roots of the corresponding flow equations. As a function of the step sizes, this results into a set of algebraic equations that need to be simultaneously solved in order to get the critical points. Hereby, we have obtained all possible critical points of the Richardson numerical integration of an arbitrary integrable real valued function. Further, in the limit of small step sizes, it is shown that the origin emerges as a repeated root of the flow equations. Apart from the origin, we find that there is a nontrivial repeated root along with two other simple roots. In the squeal, we have computed the fluctuation capacities as the pure second order partial derivatives of the Richardson numerical integration of an arbitrary real valued smooth function. In addition, we have calculated the undermining correlation between the step sizes as its mixed second order partial derivative. At a given critical point, the global stability of the Richardson numerical integration is governed by the sign of the discriminant. It is explicitly shown that the numerator of the discriminant arises as a fifth order polynomial in the step size parameter.
With the above motivation, the rest of the paper reads as follows. In section 2, we provide a brief review of the model and give highlights of the Richardson method by invoking the role of the trapezium rule. In the sequel, we discuss the Romberg method as a special case of the Richardson numerical integration. In section 3, considering the fact that the Richardson integration depends on a pair of step sizes, we give an overview of the stability analysis of a two variable function in the space of step sizes. In section 4, we optimize the Richardson integration under fluctuations of its step sizes and depict the system stability via a degree five polynomial in the square of the ratio of the step sizes around the unity. Hereby, we obtain all possible critical points of the Richardson numerical integration of arbitrary integrable function. In section 5, following the above results as in the section 4, we discuss stability criteria for the Richardson integration scheme in various limiting cases. In order to assimilate the stability structures of the Richardson integration, we offer the flowchart of our analysis, discuss various limiting cases to characterize the stability of the Richardson numerical integration of an arbitrary real valued function and provide summary of the concerned results in a tabular format. In section 6, we offer a qualitative discussion of the Richardson theory with the existing models that are commonly explored in the light of system optimizations, fluids engineering, business media and fast sea transportation. Hereby, perspective applications are anticipated in the domain of the sensitivity analysis, finite difference schemes, and non-uniform grids. In reference to the original Richardson's dream towards the atmospheric circulation and forecasting of the weather as a time integration, our study supports modeling of discrete dynamical systems, numerical simulations, system engineering and computational mechanics. Finally, in section 7, we conclude the paper with a brief summary and perspective directions for future research and developments.

Review of the model
In this section, we provide a brief review of the Richardson integration of an arbitrary real valued function. Considering a numerical method that performs the integration of a real valued function FðxÞ on a computer or a machine, it is known (Battaglia et al., 2015) that the accuracy depends on the step size h as in Eqn. (1). In practice, it follows that when the estimated value of a numerical integral tends to its true value, the step size of the integration tends to zero (Battaglia et al., 2015). In this case, let p ? be the order of the accuracy, then using step size h, it is known (Battaglia et al., 2015) that the function FðhÞ can be expressed as the following approximation Here, the true value of F is found as the 0 th order coefficient In this case, the exact value I of the integral of the function F in the trapezoidal rule (Battaglia et al., 2015), the estimated value IðhÞ and the associated true error EðhÞ (Kaw & Krivanek, 2018) can be represented as Notice that the step size h of the integral I can be obtained from its limits and the number of steps. Namely, if the integral is performed from a to b in n steps, then the step size h is given by ðb À aÞ=n. For example, the integration IðhÞ is realized as an estimated value of I and EðhÞ, that is, the true error for IðhÞ. To be precise, if we consider two different step sizes fh 1 ; h 2 g, then the true errors fEðh 1 Þ; Eðh 2 Þg and the respective estimated values fIðh 1 Þ; Iðh 2 Þg may vary, however, the exact value of I remains the same, see (Chapra & Canale, 1998) for an instance. Herewith, from Eqn. (4), it follows that we have the equality In the numerical analysis, it is known (Chapra & Canale, 1998) that the true error in the trapezoidal integration of an arbitrary real valued function F is given by Assuming that the second derivate f 00 evaluated at a given point is a constant, it follows from Eqn. (6) that the respective errors Eðh 1 Þ and Eðh 2 Þ corresponding to the step sizes h 1 and h 2 satisfy Substituting Eqn. (7) into Eqn. (5), we get the following estimation of the subsequent error From Eqns. (4,8), it follows that the integral I can be expressed as Hereby, we see that the integral I can be expressed solely as a function of the step sizes fh 1 ; h 2 g. In this paper, our motivation is to find the optimal value of I as in the above Eqn. (9). Further, as a special case, our method incorporates optimization of the trapezium integration and Romberg method, see (Chapra & Canale, 1998) for detail. Namely, from Eqn. (9), this is performed as a weighted estimation using In particular, Eqn. (10) offers an estimate of the trapezoidal scheme when the step sizes are repeatedly halved, see (Chapra & Canale, 1998) for an introduction. Moreover, in the light of the extrapolation theory (Battaglia et al., 2015), we get the Romberg Integration as the following recursion equation where j; k ¼ 1; 2; . . . ; N with N as the number of step sizes used to evaluate a given definite Richardson integral. Herewith, the Richardson method incorporates the Romberg extrapolation methods, wherefore our method inherits an intrinsic optimization of the Richardson extrapolation techniques.
With the above motivations, we consider the optimization of well-known Richardson integration with respect to the choice of the consecutive step sizes. Here, from the perspective of the intrinsic fluctuations, the objective function in our approach is to optimize the Richardson integration of an arbitrary real valued function F in the space of its step sizes. In particular, in this formulation, the task proposed is realized by combining different parts of the integral that is defined as a definite integral of the type I ¼ FðxÞdx that is sought to be computed numerically.
In this paper, one of the main motivations is to find the optimal value of the above I. Notice that I attains a fixed value when it is performed in the continuous way, that is, when there is no error. In the other words, when I is not evaluated by a numerical integration method, then there is nothing to be optimized. However, when the integration I is performed numerically, one deals with the quadrature formula. Therefore, the obtained numerical value of I needs to be optimized, that is, the exact error E governing the precision of I that depends on certain parameters should be minimized. In this setting, these parameters are considered as the consecutive step sizes that are used in the numerical evaluation of I satisfying From the above Eqns. (9, 12), it follows that in the limit of the quadrature formula (Battaglia et al., 2015), the optimization of I is the same as the minimization of the above exact error E. In this sense, for a given F as above, we optimize the Richardson integration I in the space of execution parameters fh 1 ; h 2 g, as depicted above in the Eqn. (9).
Following the same, we discuss the critical points of the quadrature formula, which is considered as a function of the two step sizes fh 1 ; h 2 g. In the light of numerical practices, the significance of the results of this discussion is that it gives the optimal value of the Richardson integration of a smooth real valued function when it is performed by a machine of finite memory storage capacity. To be precise, for a given f whose Richardson integration I is to be numerically computed on a given machine of finite memory capacity n, the inverse of n can be used to define the error of an evaluation of the Richardson integration I. Hereby, in the light of numerical practices our scheme offers an intrinsic computational significance that is equally supported by a machine. In this direction, numerical tests of the proposed approach as well as its comparison with other quadrature methods are among the promising directions that we leave open for prospective research and developments. In particular, the practical effects of the proposed optimization scheme are anticipated to be demonstrated in the future.
In nutshell, our proposal gives an intrinsic methodology to determine the optimal values of the Richardson integration of a smooth real valued function when it is numerically performed by a machine with its initial and final step size h 1 and h 2 . In addition, the condition that the quantities fðb À aÞ=h j j j ¼ 1; 2; . . . ; Ng should be natural numbers is considered in the light of derandomization theory, see (Agrawal & Biswas, 2003) for an introduction. Namely, in a derandomized system based formulation of the optimization task, these would emerge as the integer valued numerical constraints. Herewith, for a given sequence of step-sizes fh i j i ¼ 1; 2; . . . ; Ng, both step sizes are anticipated to give correct partitions of the integration interval ða; bÞ up to a computational error, that is, we have step-size spacing h iþ1 À h i in the sense of the forward difference and that of h i À h iÀ1 in the sense of the backward difference in the light of numerical studies.

Step size fluctuations
In this section, we provide an overview of the fluctuation theory in the light of the Richardson numerical integration as defined above in Eqn. (9). Considering the step sizes h 1 and h 2 of the integration of F as the fluctuation parameters, we examine stability properties of the statistical system by considering the Richardson integration Iðh 1 ; h 2 Þ as the embedding function over the surface of step sizes fh 1 ; h 2 g. Namely, we have explained the nature of the embedding map from the space of step sizes h 1 and h 2 to the set of real numbers, its critical points whether they correspond to a minimum, maximum or a saddle point and discuss the convergence of the Richardson integral Iðh 1 ; h 2 Þ of F to its true value. To be precise, given the two embedding variables as the step sizes ðh 1 ; h 2 Þ 2 M 2 , the behavior of the integral Iðh 1 ; h 2 Þ about its critical points ðh 0 1 ; h 0 2 Þ can be determined by the Taylor series expansion where the symbols fI 1 ðh 1 ; h 2 Þ; I 2 ðh 1 ; h 2 Þg denote the first order partial derivatives of Iðh 1 ; h 2 Þ and fI 11 ðh 1 ; h 2 Þ; I 22 ðh 1 ; h 2 Þ; I 12 ðh 1 ; h 2 Þg are its second order partial derivatives evaluated at a critical point. About a given critical point ðh 0 1 ; h 0 2 Þ as roots of the flow equations I 1 ðh 1 ; h 2 Þ ¼ 0 and I 2 ðh 1 ; h 2 Þ ¼ 0, then up to the second order terms, the above Taylor series as in Eqn. (13) can be expressed as the below quadratic curve in h 1 and h 2 . Here, the expansion coefficients fa; b; c; dg are defined as a ¼ Iðh 0 In order to find the roots ðh 0 1 ; h 0 2 Þ, the partial derivatives fI 1 ðh 1 ; h 2 Þ; I 2 ðh 1 ; h 2 Þg with respect to h 1 and h 2 of the integral Iðh 1 ; h 2 Þ should be equated to zero, that is, we are required to simultaneously solve the equations I 1 ðh 1 ; h 2 Þ ¼ 0 and To be precise, let the solutions of the above equations be h 1 ¼ h 0 1 and h 2 ¼ h 0 2 , then the critical point of Iðh 1 ; h 2 Þ arises as a pair ðh 0 1 ; h 0 2 Þ. In order to check the behavior of a critical point, we need to find the signature of discriminant Physically, the behavior of a critical point ðh 0 1 ; h 0 2 Þ is summarized as per the following criteria (I) If Δ < 0 at ðh 0 1 ; h 0 2 Þ, then the critical point ðh 0 1 ; h 0 2 Þ is a saddle point of Iðh 1 ; h 2 Þ. (II) If d 1 4 < À n 0 =n 4 at Δ 2 < 0, then we have the following classification.
then the test fails and we say that the function Iðh 1 ; h 2 Þ is too flat. Therefore, we need a higher derivative test in order to examine the stability properties of Iðh 1 ; h 2 Þ under fluctuations of its step sizes fh 1 ; h 2 g for a chosen Richardson numerical integration.
In this regard, the aim is to optimize the Richardson integration I of an arbitrary integrable function under variations of its step sizes fh 1 ; h 2 g. In the next section, as fh 1 ; h 2 g vary, we study fluctuations of the Richardson integration Iðh 1 ; h 2 Þ about is critical points ðh 0 1 ; h 0 2 Þ. Hereby, the integration I reaches a fixed value when h 1 and h 2 are taken at their critical values ðh 0 1 ; h 0 2 Þ. This follows from the aforementioned fact that the numerical computation of a definite integral possesses an exact error. Further, the optimization of I is in correspondence with the optimization the exact error E. Moreover, up to a given nonzero maximum error e ¼ maxfEg, from Eqn. (12), it follows that the extrimization of I is equivalent to the extrimization of the exact error E. In the light of numerical practices, as a function of the step size parameter, the significance of the results is that our analysis gives a quintic type stability criterion as far as the optimization of the Richardson integration of an arbitrary integrable function is concerned under variations of its step sizes, see Eqn. (36) below.

Optimal Richardson schemes
In this section, we provide optimization properties of the Richardson integration Iðh 1 ; h 2 Þ of an arbitrary real valued function F with the integral Iðh 2 Þ when it is evaluated at the step size h 2 and the integral Iðh 1 Þ at that of the step size h 1 . In order to simplify the calculations, let us define d 1 and d 2 as With the above notations, the Richardson integration as in Eqn. (9) reads as the following function Hereby, the first order partial derivatives of the function Iðh 1 ; h 2 Þ with respect to the step sizes h 1 and h 2 give the flow components of Iðh 1 ; h 2 Þ. For a given pair of fd 1 ; d 2 g, the derivation of the Richardson integration Iðh 1 ; h 2 Þ give the flow components fI 1 ðh 1 ; h 2 Þ; I 2 ðh 1 ; h 2 Þg. Namely, we notice that the differentiation of Iðh 1 ; h 2 Þ with respect to h 1 results into the following expression Similarly, by differentiating Iðh 1 ; h 2 Þ with respect to h 2 , we find the associated flow component as Differentiations of the above Eqns. (19, 20) with respect to h 1 and h 2 give the respective fluctuation capacities that are define as the second order pure partial derivatives fI 11 ðh 1 ; h 2 Þ; I 22 ðh 1 ; h 2 Þg. In order to simplify the notations, we express the fluctuation capacities fI 11 ðh 1 ; h 2 Þ; I 22 ðh 1 ; h 2 Þg as the function of d 1 . Differentiating Eqn. (19) with respect to h 1 , we find the following h 1 fluctuation capacity where the coefficients fa 1 ; b 1 ; c 1 g are given by Similarly, by differentiating Eqn. (20) with respect to h 2 , we get the h 2 fluctuation capacity as where the coefficients fa 2 ; b 2 ; c 2 g are given by the expressions For the above Richardson integration as in Eqn. (18), the correlation between step sizes h 1 and h 2 is defined by the mixed partial derivative of I. In this case, by differentiating Eqn. (19) with respect to h 2 , as the function of d 1 , we find that the correlation I 12 arises as the following relation where the coefficients fa 3 ; b 3 ; c 3 g read as In this setting, the critical points of the Richardson numerical integration are given as roots of I 1 ðh 1 ; h 2 Þ and I 2 ðh 1 ; h 2 Þ. Therefore, we are required to simultaneously solve the flow equations: I 1 ðh 1 ; h 2 Þ ¼ 0 and I 2 ðh 1 ; h 2 Þ ¼ 0. From the above Eqn. (19), we see that I 1 ðh 1 ; h 2 Þ ¼ 0 yields Similarly, the other flow equation I 2 ðh 1 ; h 2 Þ ¼ 0 as in Eqn. (20) gives the following equalities For two distinct step sizes h 1 Þ h 2 and integral values Iðh 1 Þ Þ Iðh 2 Þ, by dividing Eqn. (27) by (28), we see that the step sizes fh 1 ; h 2 g satisfy the following ratio Therefore, we see that the ratio of step sizes fh 1 ; h 2 g arises as the ratio of the first derivatives fI 1 ðh 1 Þ; I 2 ðh 2 Þg at a critical point of I. This defines the equilibrium condition of the Richardson integration when fh 1 ; h 2 g vary. Squaring both side of Eqn. (29) and substituting the corresponding value of h 2 into Eqn. (27), we get the following values Substituting the above value of h 1 as in Eqn. (30) into Eqn. (28), we find that h 2 reads as follows In order to get all the solutions of the equations I 1 ðh 1 ; h 2 Þ ¼ 0 and I 2 ðh 1 ; h 2 Þ ¼ 0, for example, in the limit of h 1 ! h 2 , we need to directly solve the above Eqns. (27,28 In order to study the stability analysis of the Richardson numerical integration, let us consider a critical point ðh i 1 ; h i 2 Þ 2 <. Then, the type of fluctuations about the critical point ðh i 1 ; h i 2 Þ is governed by the sign of the discriminant as in Eqn. (15). Substituting the value of fluctuation capacities fI 11 ; I 22 g as in Eqns. (21,23), and the local correlation I 12 as in Eqn. (25), it is not difficult to show that the discriminant Δ reading as can be expressed as per the following ratio Here, the factors fn i j i ¼ 0; 1; 2; 3; 4; 5g as in Eqn. (34) are defined as From Eqn. (34), as a function of the step size parameter d 1 , it follows that the global stability of the Richardson integration of an arbitrary real valued function F is governed by the sign of the above quintic polynomial Namely, in the step size parameter d 1 , the roots of the above quintic polynomial nðd 1 Þ as in Eqn. (36) satisfy the quintic equation nðd 1 Þ ¼ 0. Here, as an algebraic polynomial, nðd 1 Þ ¼ 0 as in Eqn. (36) can be factorized by using Galois theory (Artin, 2007). One can solve such an algebraic equation nðd 1 Þ ¼ 0 when it is a degree five or reducible to lower degree equations, see (Kappe & Warren, 1989) for an introduction. In the sequel, we provide an explicit analysis of the stability of the Richardson integration of an arbitrary real valued function F under fluctuations of the step sizes h 1 and h 2 . In particular, in the next section, the stable configurations are categorized under various values of the model parameters fn i j i ¼ 0; 1; 2; 3; 4; 5g.

Discussion of the results
In this section, we provide a flow chart towards the classification of stability criteria of the Richardson numerical integration scheme and its limiting configurations under various cases of the modal parameters fn i j i ¼ 0; 1; 2; 3; 4; 5g. To offer a step by step picture of the same, we give the flowchart, discussion of special cases and summary of the results concerning the stability of the Richardson numerical integration of an arbitrary real valued function.

A flowchart of the stability analysis
To be precise, for the aforementioned Richardson numerical integration of an arbitrary real valued function, let us define the associated quantities as below: (I) Choose the step sizes h 1 ; h 2 2 R for a given I : R ! R.
Following the above steps, the flowchart of the stability analysis is given in the above Figure 1.
In the sequel, we study the stability of the Richardson integration, where we focus our attention on various special cases of the above discriminant Δ as in Eqn. (34). Namely, as a function of d 1 , we calculate its signature in various limiting cases as below.

Results in special cases
In this subsection, we give the stability properties of the Richardson numerical integration scheme when Δ has a specific set of non-vanishing fn i j i ¼ 0; 1; 2; 3; 4; 5g with a constant, linear, quadratic, cubic, quartic, quintic type of numerators as a function of the step size parameter d 1 .

Constant numerator
In the sequel, we discuss various possible cases of the discriminant Δ with a constant numerator.
In this case, from Eqn. (34), we see that the discriminant reads as follows Therefore, there is a positive discriminant Δ if n 0 > 0. In other words, the Richardson numerical integration remains stable if c 2 3 < À c 1 c 2 . In this case, we have either a maximum or a minimum according as I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0, or I 11 ðh 1 ; h 2 Þ > 0 or I 22 ðh 1 ; h 2 Þ > 0 respectively. We observe that there is a saddle point behavior for c 2 3 > À c 1 c 2 under fluctuations of fh 1 ; h 2 g.
Similar to the above case, we have the following fluctuating configurations with a different numerator fn i j i ¼ 0; 1; 2; 3; 4; 5g as below: Substituting the above values of n i 's, we find the following discriminant Since d 1 4 is always positive, the sign of the discriminant Δ remains the same as that of n 2 . In other words, the Richardson numerical integration is stable when b 1 b 2 >a 1 c 2 þ a 2 c 1 þ b 3 2 . In this case, there is either a maximum or a minimum according as I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0, or I 11 ðh 1 ; h 2 Þ>0 or I 22 ðh 1 ; h 2 Þ>0 respectively. Note that there is a saddle point behavior for n 2 < 0, that is, we have b 1 b 2 < a 1 c 2 þ a 2 c 1 þ b 3 2 , when the step sizes fh 1 ; h 2 g fluctuate.
Case 3. Configurations with n 0 ¼ 0 ¼ n 1 ¼ n 2 ¼ n 3 ¼ n 5 In this setup, we find that the discriminant reads as follows Hereby, we have a positive discriminant Δ for n 4 >0. Namely, the Richardson numerical integration remains stable if a 2 ðb 1 À a 1 Þ > 0. In this case, we have either a maximum or a minimum according as I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0, or I 11 ðh 1 ; h 2 Þ > 0 or I 22 ðh 1 ; h 2 Þ > 0 respectively. Equally, under fluctuations of fh 1 ; h 2 g, it follows that there is saddle point when b 1 < a 1 for a positive a 2 and b 1 > a 1 for a negative value of a 2 .
Case 4. Configurations with n 0 ¼ 0 ¼ n 2 ¼ n 3 ¼ n 4 ¼ n 5 With the above values of n i , we see that the discriminant is given by In this case, the discriminant Δ takes a positive value if n 1 and d 1 have the same sign. In particular, the Richardson numerical integration scheme remains stable if b 1 c 2 < b 2 c 1 þ 2b 3 c 3 holds for d 1 < 0 or b 1 c 2 > b 2 c 1 þ 2b 3 c 3 holds for d 1 > 0 with I 11 ðh 1 ; h 2 Þ > 0 or I 22 ðh 1 ; h 2 Þ > 0. Therefore, the integration I reaches a maximum value according as I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0 for b 1 c 2 < b 2 c 1 þ 2b 3 c 3 and d 1 < 0, or b 1 c 2 > b 2 c 1 þ 2b 3 c 3 holds for d 1 > 0, respectively. Further, the integration scheme attains a saddle point if Δ < 0, that is, the factors n 1 and d 1 have the opposite sign under fluctuations of fh 1 ; h 2 g.
Case 5. Configurations with n 0 ¼ 0 ¼ n 1 ¼ n 2 ¼ n 4 ¼ n 5 With the choice of n i 's as above, it follows that we have the following discriminant Therefore, there is a positive discriminant Δ if n 3 and d 1 have the same sign. In other words, the Richardson numerical integration scheme remains stable if b 1 a 2 < a 1 b 2 þ c 1 a 2 for d 1 < 0 or b 1 a 2 > a 1 b 2 þ c 1 a 2 for d 1 > 0 with I 11 ðh 1 ; h 2 Þ > 0 or I 22 ðh 1 ; h 2 Þ>0. In this case, we have a maximum according as I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0 for b 1 a 2 <a 1 b 2 þ c 1 a 2 and d 1 < 0, or b 1 a 2 >a 1 b 2 þ c 1 a 2 for d 1 > 0, respectively. Further, the configuration attains a saddle point if Δ < 0, that is, the factors n 3 and d 1 have the opposite sign under fluctuations of fh 1 ; h 2 g.
Case 6. Configurations with n 0 ¼ 0 In this case, it follows that we have the following discriminant Thus, there is a positive discriminant Δ if n 5 and d 1 have the same sign. In other words, for d 1 < 0 the Richardson numerical integration remains stable if a 1 a 2 > 0, viz. either a 1 and a 2 have the same sign for d 1 < 0, or a 1 and a 2 have the opposite sign for d 1 > 0 with I 11 ðh 1 ; h 2 Þ>0 or I 22 ðh 1 ; h 2 Þ > 0.
As a matter of the fact, the function I attains a maximum according as I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0 when a 1 and a 2 have the opposite sign with d 1 >0, or a 1 and a 2 have the same sign for d 1 <0, respectively. Indeed, the above configuration attains a saddle point behavior whenever Δ < 0, that is, when the factors n 5 and d 1 have the opposite sign under fluctuations of the step sizes fh 1 ; h 2 g of the integration, that is, either a 1 ; a 2 have the same sign for d 1 >0 or they have the opposite sign for d 1 < 0.

Linear numerator
In the sequel, we discuss stability properties of the Richardson numerical integration scheme when Δ has two of non-vanishing fn i j i ¼ 0; 1; 2; 3; 4; 5g with a linear type numerator in the factor d 1 . The undermining cases are as below.
Case 7. Configurations with n 2 ¼ 0 ¼ n 3 ¼ n 4 ¼ n 5 In this case, the numerator of the discriminant Δ becomes a linear function of the factor d 1 . Namely, it follows from Eqn. (34) that we have the following discriminant Thus, the sign of Δ will depend on positivity or negativity of n 1 d 1 þ n 0 . In other words, we have a maximum or a minimum according as d 1 > À n 0 =n 1 with I 11 ðh 1 ; h 2 Þ<0 or I 22 ðh 1 ; h 2 Þ<0, or I 11 ðh 1 ; h 2 Þ>0 or I 22 ðh 1 ; h 2 Þ>0 respectively.
On the other hand, we see that the Richardson integration of an arbitrary real valued function passes through a saddle point if the factor d 1 satisfies the inequality d 1 < À n 0 =n 1 where n 0 and n 1 are defined as in Eqn. (35).
Case 8. Configurations with n 0 ¼ 0 ¼ n 3 ¼ n 4 ¼ n 5 Following the above case, another linear equation in the numerator of Δ can be defined as With the above conditions, it follows that Δ has a positive sign when d 1 > 0 and n 2 d 1 þ n 1 > 0 or d 1 < 0 and n 2 d 1 þ n 1 < 0. In this case, we have either a maximum according as I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0, or vice-versa for a minimum, respectively. Similarly, the integration scheme passes through a saddle point whenever we have Δ < 0, that is, d 1 and n 2 d 1 þ n 1 have the opposite signature.
Case 9. Configurations with n 0 ¼ 0 ¼ n 1 ¼ n 4 ¼ n 5 In this case, we have the following discriminant With the above specifications, the sign of Δ solely depends on the sign of the numerator n 3 d 1 þ n 2 . In this concern, it follows that the Richardson integration scheme reaches a maximum or a minimum according as d 1 > À n 2 =n 3 with I 11 ðh 1 ; h 2 Þ<0 or I 22 ðh 1 ; h 2 Þ < 0, or I 11 ðh 1 ; h 2 Þ>0 or I 22 ðh 1 ; h 2 Þ>0 respectively. On the other hand, we notice that the Richardson integration goes through a saddle point if d 1 < À n 2 =n 3 in the space of step sizes fh 1 ; h 2 g.
Case 10. Configurations with n 0 ¼ 0 ¼ n 1 ¼ n 2 ¼ n 5 Substituting the values n 0 ¼ 0 ¼ n 1 ¼ n 2 ¼ n 5 in Eqn. (34), the discriminant is Therefore, it follows that there is a positive Δ when we have either d 1 > 0 and n 4 d 1 þ n 3 > 0 or d 1 < 0 and n 4 d 1 þ n 3 < 0. With the above conditions, we have a maximum according as I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0, or vice-versa for a minimum, respectively. Similarly, it follows that we have a saddle point behavior if the discriminant Δ < 0, that is, the factors d 1 and n 4 d 1 þ n 3 have the opposite sign under fluctuations of the step sizes fh 1 ; h 2 g.
Case 11. Configurations with n 0 ¼ 0 ¼ n 1 ¼ n 2 ¼ n 3 In the above limiting setting, the discriminant Δ simplifies as Therefore, the discriminant Δ remains positive if n 5 d 1 þ n 4 > 0. In this case, the Richardson integration scheme attains a maximum according as the inequality d 1 > À n 4 =n 5 with I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0, or vice-versa for a minimum, respectively. On the other hand, it follows that the Richardson integration scheme passes through a saddle point if we have d 1 < À n 4 =n 5 .

Quadratic numerator
In the next step, we examine stability properties of the Richardson numerical integration of an arbitrary real valued function when Δ has three of non-vanishing factors fn i j i ¼ 0; 1; 2; 3; 4; 5g with a quadratic type numerator as a function of d 1 . The undermining cases are as below.
Case 12. Configurations with n 3 ¼ 0 ¼ n 4 ¼ n 5 By substituting the above values of n 3 ; n 4 ; n 5 , we see that Δ simplifies as the quadratic equation In this case, in order to determine the sign of Δ, we may define d 1 AE ¼ Àn 1 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi n 2 1 À 4n 2 n 0 q 2n 2 : As in the case 13, it follows that we have a positive Δ, if On the other hand, for a negative Δ, that is, for Richardson integration Iðh 1 ; h 2 Þ approaches a saddle point in the space of fh 1 ; h 2 g.
Case 13. Configurations with n 1 ¼ 0 ¼ n 3 ¼ n 5 In this case, it follows from Eqn. (34) that the discriminant Δ reduces as In order to determine the sign of Δ, we may definẽ Here,d 1 AE are roots of the above Eqn. (50), when it is solved for d 1 2 . In an analogy with the previous case, we see that Δ is positive, if d 1 satisfies either d 1 2 <d 1 À or d 1 2 >d 1 þ . In particular, we have a stable Richardson integration scheme when either I 11 ðh 1 ; h 2 Þ > 0 or I 22 ðh 1 ; h 2 Þ > 0 holds. There are maxima, viz. we have is an unstable Richardson integration scheme when Δ > 0 and I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0. On the other hand, we have a saddle point when Δ takes a negative value, that is, we haved 1 À < d 1 2 <d 1 þ . In this case, as well as that of the case of a maximum, it is expected that the associated Richardson integration diverges.
Case 14. Configurations with n 0 ¼ 0 ¼ n 2 ¼ n 4 In the above special case, we see that the numerator of Δ simplifies as a quadratic equation in In order to determine the sign of Δ, let us define d AE 1 ¼ Àn 3 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 3 2 À 4n 5 n 1 p 2n 5 : (53) With the above notation, it follows that Δ can be factorized as Subsequently, we see that the Richardson integration Iðh 1 ; h 2 Þ attains a maximum value when Δ > 0 (that is, either for d 1 > 0 we have d 1 2 > d 1 þ or d 1 2 < d 1 À or for d 1 < 0 the condition d 1 À < d 1 2 <d 1 þ holds) such that I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0. Correspondingly, for a given Iðh 1 ; h 2 Þ as a function of ðh 1 ; h 2 Þ, it follows that there is a minimum when Δ > 0 (that is, either Case 15. Configurations with n 0 ¼ 0 ¼ n 4 ¼ n 5 In this specification, from Eqn. (34), we notice that the discriminant Δ reads as As the root of the numerator as above in Eqn. (55) in the step size parameter d 1 , this configuration follows the same analysis as in the case 14 with the following replacement d 1 AE ¼ Àn 2 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 2 2 À 4n 3 n 1 p Case 16 . Configurations with n 0 ¼ 0 ¼ n 1 ¼ n 5 In this case, the discriminant Δ reduces as per the following expression As discussed above, the corresponding stability analysis follows the case 12 with the replacement d 1 AE ¼ Àn 3 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 3 2 À 4n 4 n 2 p 2n 4 : Case 17. Configurations with n 0 ¼ 0 ¼ n 1 ¼ n 2 With the choice of vanishing coefficients fn 0 ; n 1 ; n 2 g, we observe the following quadratic discriminant Hereby, the corresponding stability analysis follows the case 14 by defining d 1 AE as follows d 1 AE ¼ Àn 4 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 4 2 À 4n 5 n 3 p 2n 4 : The stability analysis of the Richardson integration Iðh 1 ; h 2 Þ remains the same as in the case 14 with d 1 AE as the turning values of d 1 , that are, the roots of Δ as in Eqn. (59).
In the sequel, we discuss the stability properties of the Richardson numerical integration when Δ has a pair of non-vanishing fn i j i ¼ 0; 1; 2; 3; 4; 5g with a quadratic type numerator. The concerned cases are discussed as the below configurations.
Case 18. Configurations with n 1 ¼ 0 ¼ n 3 ¼ n 4 ¼ n 5 In this case, it is direct to see that we have the following discriminant This configuration follows the same analysis as in the case 7 with the replacement of d 1 2 as d 1 .
Namely, the sign of Δ depends on the positivity or the negativity of the resulting numerator n 2 d 1 2 þ n 0 . In other words, we have a maximum according as d 1 2 > À n 0 =n 2 with I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0, or a minimum if I 11 ðh 1 ; h 2 Þ > 0 or I 22 ðh 1 ; h 2 Þ > 0 for Δ > 0. On the other hand, we see that the Richardson integration passes through a saddle point if d 1 2 < À n 0 =n 2 . Thus, in order to have a stable Richardson integration scheme, we need to appropriately choose the factors fn 0 ; n 2 g that is achieved by optimally choosing the step sizes h 1 and h 2 .
Case 19. Configurations with n 1 ¼ 0 ¼ n 2 ¼ n 3 ¼ n 5 With the above specification, we have the following discriminant This configuration has the same analysis as in the case 7 with replacement of d 1 4 as d 1 . Namely, the sign of Δ will depend on positivity or negativity of n 4 d 1 4 þ n 0 . In other words, we have a maximum according as d 1 4 > À n 0 =n 4 with I 11 ðh 1 ; h 2 Þ < 0 or I 22 ðh 1 ; h 2 Þ < 0, or a minimum if I 11 ðh 1 ; h 2 Þ > 0 or I 22 ðh 1 ; h 2 Þ > 0 respectively. In contrast, we see that the Richardson integration scheme passes through a saddle point if d 1 4 < À n 0 =n 4 under fluctuations of the step sizes fh 1 ; h 2 g.
Case 20. Configurations with n 0 ¼ 0 ¼ n 2 ¼ n 4 ¼ n 5 With the above vanishing coefficients fn 0 ; n 2 ; n 4 ; n 5 g, we find the following discriminant The above configuration follows the same analysis as in the case 8 with the replacement of d 1 as.
Case 21. Configurations with n 0 ¼ 0 ¼ n 2 ¼ n 3 ¼ n 4 In this case, when the coefficients fn 0 ; n 2 ; n 3 ; n 4 g vanish, we have the following discriminant This follows the same stability analysis as in the case 8 with the replacement of d 1 as Given that the coefficients fn 0 ; n 1 ; n 3 ; n 5 g vanish, we have the following discriminant This system follows the same stability analysis as in the case 9 with the replacement of d 1 as d Case 23. Configurations with n 0 ¼ 0 ¼ n 1 ¼ n 2 ¼ n 4 In this case, substituting n 0 ¼ 0 ¼ n 1 ¼ n 2 ¼ n 4 in Eqn. (34), we find the following discriminant The concerning stability analysis can be carried out as in the case 10 with the replacement of d 1 as In the light of the Galois field theory, Galois groups, algebraic curves and commutative algebra (Igoe, 1967;Kappe & Warren, 1989), the other limiting cases of the discriminant Δ that are corresponding to the cubic (Nickalls, 1993), quartic (Lazard, 1988;Rees, 1922) and quintic (Alekseev, 2004;Křížek & Somer, 2015;Ruffini, 1799) numerators as a function of the factor d 1 are discussed as below.

Cubic numerator
For the given choice of fn 0 ; n 1 ; n 2 ; n 3 ; n 4 ; n 5 g as in Eqn. (35), the cubic numerator configurations can be expressed as the determinant withñ 3 Þ 0. For the above specific value of the discriminant Δ as in Eqn. (71), it follows (Nickalls, 1993) that the solutions fd ðkÞ 1 j k ¼ 0; 1; 2g of the corresponding cubic equation Here, the parameters Δ 0 ; C 1 and ω are referenced in the Appendix A. For k ¼ 0; 1; 2, the classification of roots d ðkÞ 1 as in Eqn. (72) of the above cubic equation is determined by the cubic discriminant Δ 2 as in Appendix A, viz. Eqn. (86). Namely, we have three distinct real roots of Δ when Δ 2 > 0. When Δ 2 < 0, we have two complex conjugate roots and a real root of Δ. In the case of Δ 2 ¼ 0, we have multiple real roots of Δ, see (Nickalls, 1993) for details. The concerned results concerning to a cubic numerator in d 1 are summarized in the Table 5 as in the next subsection.

Quartic numerator
For the choice of fn i ji ¼ 0; 1; . . . ; 5g as in Eqn. (35), the quartic numerator configurations can be described with the following fluctuation discriminant with n 4 Þ 0 as the limiting values of the discriminant Δ as in Eqn. (73). In this case, we see (Lazard, 1988;Rees, 1922) of the solutions fd ðiÞ 1 j i ¼ 1; 2; 3; 4g of the quartic equation ∑ 4 i¼0 n i d 1 i ¼ 0 are given by Here, the parameters fS; p; qg of the solutions as in the above Eqns. With the respective value of the parameters fP; R; D;Δ 0 ;Δ 2 g as given in Appendix B, the classification of its roots (Lazard, 1988;Rees, 1922) is summarized as below (I) ForΔ 2 < 0, the above quartic equation possesses two distinct real roots and a pair of nonreal complex conjugate roots as fd (i) For P < 0, there exist two real repeated roots.
(ii) For P > 0 and R ¼ 0, there exist two repeated complex conjugate roots.
The above behavior of the discriminant Δ is summarized in the next subsection.

Quintic numerator
In the sequel, we provide most general stability structure of the Richardson integration scheme, where it's fluctuation discriminant Δ is considered as a fifth degree polynomial in the step size parameter d 1 . In general, to discuss the stability structure, we need to factorized the quintic polynomial as in Eqn. (34) with n 5 Þ 0. As discussed before, because of the odd degree characteristic behavior of a quintic equation with n i as in Eqn. (35), the qualitative behavior of such a configuration appears to that of the cubic equation. Note that the most of the quintic equations cannot be solved in general because of the Abel-Ruffini theorem (Alekseev, 2004;Křížek & Somer, 2015;Ruffini, 1799). However, it is worth mentioning that some of the quintics can be solved explicitly. Such general solutions are highly non-trivial to be used in practice, wherefore one calculates them by an algorithmic prescription to obtain its roots.
Notice that solvable quintics are defined by a reducible polynomial that can be factorized as a product of irreducible radicals. For the solvable quintics, we may define d 1 ¼ y À n 4 =5n 5 . This results the above Eqn. (76) into a reduced quintic y 5 þ ty 3 þ uy 2 þ vy þ w ¼ 0: Hereby, the reduced parameters ft; u; v; wg are described in Appendix C. The above reduced quintic as in Eqn. (77) is solvable if it is factorizable into lower degree polynomials with rational coefficients. Namely, it is factorized whenever the Cayley-resolvent (Farb & Wolfson, 2018) defined as has a rational root in z. Here, the quintic discriminant Δ 0 and P are defined in Appendix C. In particular, the above quintic is said to be solvable if the parameter is a rational number. Hereby, a quintic is solved as a result of the Cayley-resolvent (Farb & Wolfson, 2018) that allows to test whether a given quintic is solvable or not.
Further, there exist implicit formulae such as Young/Lazard solutions for solving solvable quintics, see (Gray, 2018) for a historical introduction in the light of the abstract algebra. There are interesting solvable quintics such as Bring-Jerrard forms, Jacobi theta functions over a lattice index and number fields, elliptic modular forms, generalized hyper-geometric functions and their differential reductions, and computational aspects of a low genus curve (Bytev & Kniehl, 2016;Choie, Park, & Zagier, 2017;Skoruppa, 2015). In this regard, concerning solutions of quintic equations find their relations to algebraic geometry, viz. the rational curves, a pair of Calabi-Yau spaces, Abel-Jacobi maps, mirror symmetry, quintic 3-folds and others (Candelas et al., 1991;Clemens, 1984;Cotterill, 2012;Cox & Katz, 1999;Johnsen & Kleiman, 1996). As per this characterization, the behavior of the fluctuation discriminant Δ is summarized in the next subsection.

Summary of the results
In this section, we summarize the results concerning stability properties of the Richardson numerical integration of an arbitrary real valued function as below.

Configurations with a constant numerator
From the observation of the cases (1, 2, 3, 4, 5, 6) having a single nonzero n i , we find the following general expression of the discriminant The behavior of the above constant numerators with a nonzero n i is tabulated in the Table 1.

Configurations with a linear numerator
For the above cases (7,8,9,10,11) with two nonzero fn i ; n iþ1 g, it follows that we have The classification of the stability of the Richardson integration of an integrable function corresponding to the above linear cases of Δ as in Eqn. (81) arising with two nonzero n i is tabulated in Table 2 below

Configurations with a quadratic numerator
With the aforementioned observations as in the previous subsection, the quadratic cases of the numerator of Δ can be divided into two categories. The first arises as the case of the three nonzero fn i ; n iþ1 ; n iþ2 g or the three nonzero fn i ; n iþ2 ; n iþ4 g. The second arises as the case of the two nonzerofn i ; n iþ2 g or the two nonzero fn i ; n iþ4 g. The respective values of the discriminant Δ are given in Appendix D. In the sequel, we provide the tables corresponding to two different types of expressions of the discriminant. The stability properties of quadratic cases with three nonzero n i is tabulated in the Table 3 below.
The stability structures of the above quadratic cases with two nonzero n i is summarized in the below Table 4.

Configurations with a cubic numerator
As the cubic polynomial represented in Eqn. (71), there are four types of cubic discriminants that can be defined. They are the discriminants with four, three and two nonzero values of n i . Their explicit expressions are referred in Appendix E. The stability properties of the cubic cases for different nonzero n i are tabulated in the below Table 5.

Configurations with a quartic numerator
The quartic numerator configurations can be categorized via their four types of nonzero discriminant as described in Appendix F. For various nonzero n i , the stability structures of such quartic cases are summarized in the below Table 6.

Configurations with a quintic numerator
The quintic numerator configurations can be divided into sixteen subcases of a nonzero discriminant. Their explicit values are given in Appendix G. The stability properties of the quintic cases for various nonzero n i are summarized in the Table 7 below.

Future scope of the work
In this section, having provided the optimization analysis of the Richardson numerical integration for an arbitrary real valued function, we discuss its perspective applications under variations of its step sizes. Namely, we examine the behavior of flow components as the first derivative and fluctuation capacities as the pure second derivatives of the Richardson numerical integration for an arbitrary real valued function. The critical points of the Richardson integration scheme arise as the roots of its flow equations. Hereby, for a given real valued smooth function, we have shown that its critical points are expressible as a function of roots of the flow components. At a given pair of step sizes, this yields a pair of quadrature equations in the step size parameters d 1 . As a function of the step sizes, these equations are simultaneously solved in order to get critical points of the Richardson integration of a given real valued function.
Hitherto, from the perspective of reliability theories, we have obtained all possible critical points of the Richardson numerical integration of an arbitrary smooth real valued function. Furthermore, in the space of step sizes, we have shown that the origin is a repeated root of the flow equations. Namely, it is shown that the origin arises as a nontrivial repeated root along with four other simple roots as the critical points of the Richardson integration of a real valued function, see Eqn. (32). In order to discuss the stability of an arbitrary Richardson integration, we have computed the concerned fluctuation capacities as its pure second order partial derivatives.
In the sequel, we have calculated the undermining local correlation between the step sizes as its mixed second order partial derivative. At a given critical point, we have offered the global stability of the Richardson numerical integration via the sign of the fluctuation discriminant. It is explicitly demonstrated that the numerator of the fluctuation discriminant arises as a fifth order polynomial in the aforementioned step size parameter d 1 of the Richardson integration.

Sensitivity analysis
Optimization of the Richardson numerical integration of an arbitrary real valued smooth function finds applications in various domains of science and engineering. For example, in the light of the optimization theory and Richardson extrapolation models, a sensitivity analysis can be performed for analog circuits with a multi-objective function. Namely, by considering an evolutionary approach, a multi-parameter sensitivity analysis exists in a number of recent studies. Such investigations find their relevance in estimating the partial derivatives that govern the sensitivities of performances of a given amplifier. Thereby, one can compare the variations of sizes of the metal-oxide-semiconductor field-effect-transistors (MOSFET), see (Guerra-Gómez, Tlelo-Cuautle, & Luis, 2013) for an overview of the Richardson extrapolation analysis, multi-stage optimization and their applications to analog circuits.

Finite difference schemes
Applications of the Richardson extrapolation model are in practice as finite difference schemes towards the solution of ordinary and partial differential equations with an increased accuracy. The success of a finite linear Richardson model depends on the existence a rapidly converging asymptotic expansion over errors of the solution (Blum, Lin, & Rannacher, 1986) as a finite linear sum. Hereby, as an asymptotic error expansion, the Richardson extrapolation models offer a numerical understanding of a system of difference equations as a combination of finite linear elements. In this regard, the Ref. (Roos, Stynes, & Tobiska, 2008) offers robust numerical solutions for certain singularly perturbed differential equations and that of the Ref. (Křížek & Pekka, 1987) provides the associated super-convergence techniques.

Non-uniform grids solutions
As an application towards non-uniform grid solutions, the Richardson extrapolation model can be used to examine the grid independence of the solutions with a given grid refinement parameter. A solution with variable grid refinement factors in the coordinate directions is anticipated to offer the optimal numerical solution with the right degree of the convergence (Celik & Karatekin, 1997). The associated safety factors are discussed in (Xing & Stern, 2010) in the light of the fluid mechanics. There exists a number of engineering optimizations with non-uniform grid solutions. This includes the optimal shape design, dimensional reduction, geometric variability assessment, meta-models and deterministic particle swarm, see (Chen et al., 2015) towards global robust architecture techniques and high-fidelity optimizations.

Numerical weather forecasting
The fluid mechanical global optimization methods offer applications of our model towards the test of the system reliability (Chen et al., 2015;Xing & Stern, 2010) and robust design techniques for ships operability in an uncertain ocean environment in real time (Diez, Chen, Campana, & Stern, 2013). Hereby, our model supports the optimal numerical weather prediction (Bauer, Thorpe, & Brunet, 2015;Mason, 1986) in the light of the Richardson's original dream as the atmospheric circulation and forecasting of the weather as a time integration.

Discrete dynamical systems
Our analysis is well-suited towards the intrinsic stability examination of discrete dynamical systems such as a heuristic search planning over certain hybrid domains (Ramirez, Scala, Haslum, & Thiebaux, 2017), stochastic modeling in business industry (Zhi-Sheng & Xie, 2015) and their role in finance engineering (Lamberton & Lapeyre, 2011), charge carrier transport and ion vacancy dynamics, solar technology and their mathematical modeling (Courtier, Richardson, & Foster, 2018). Further studies exist towards the Runge-Kutta type numerical integration techniques (Kalogiratou, Monovasilis, Psihoyios, & Simos, 2014). This comprises of a system of ordinary or partial differential equations and their stability analysis as well as the error adjustment, case control analysis and applications in environmental sciences, e.g., Canadian INTEROCC study (Oraby et al., 2018).

Summary and conclusion
In this paper, we have focused on optimization of the Richardson numerical integration of an arbitrary real valued integrable function under the variations of its step sizes. Considering the fact that there always exists an error in a numerical integral of an arbitrary real valued function, we have investigated the nature of its optimization and convergence properties in the space of step sizes. In the sequel, we have shown that our method optimizes the value of a definite Richardson integral than its conventional evaluations. Specifically, in the space of step sizes, we have examined the nature of stabilities of the Richardson integral of an arbitrary real valued function about its critical points. Mathematically, our analysis shows that the optimal value of the Richardson integral arises directly via its extremization over step sizes. In general, we find that such an optimization holds for any smooth real valued function. To be precise, our investigation offers a stable Richardson numerical integration scheme with the optimal value of errors. This results in an intrinsic classification of the Richardson integration of arbitrary integrable functions in the space of step sizes.
Richardson method is among the most efficient numerical integrations that is analyzed under variations of its multiple step sizes fh 1 ; h 2 g. We have examined the optimization of the Richardson model that contains the trapezium rule and Romberg integration as its special cases. With the fact that the Richardson integration depends on a pair of step sizes, we have investigated its stability as a two variable function. Our analysis optimizes the Richardson integration of an arbitrary real valued function under fluctuations of its step sizes. Hereby, the critical points of the Richardson integration are defined as a pair of flow equations in the step sizes. For an arbitrary real valued smooth function, the critical points arise as roots of the flow equations. In particular, the flow equations that arise as a function of the step sizes need to be simultaneously solved in order to get the critical points of the Richardson integration. In general, this yields a pair of nonlinear equations in h 1 and h 2 .
To be specific, given two distinct step sizes fh 1 ; h 2 g, we have obtained all possible critical points of the Richardson integration of an arbitrary integrable function. Consequently, the fluctuation capacities as the auto-correlations between the step sizes are computed as the pure second order derivatives of the Richardson integration. At a given critical point, the global stability of the Richardson numerical integration is governed by the sign of the fluctuation discriminant. Hereby, we find that the discriminant arises as a fifth order polynomial in the step size parameter d 1 , viz. the square of the ratio of the step sizes around the unity. Following the above estimation, we have discussed the nature of critical points of the Richardson numerical integration and classified its stability domains in various limiting cases of the discriminant. In order to make our analysis comprehensible, we have provided a flowchart and tabular descriptions of the results toward an intrinsic classification of the stability structures of the Richardson integration of arbitrary real valued functions. Subsequently, under variations of the step sizes fh 1 ; h 2 g, a detailed discussion of the stability properties of the Richardson numerical integration is provided in the cases when the numerator of its fluctuation discriminant arises as an equation of the degree less than or equal to five in the factor d 1 . In the light of optimization theory, associated algebraic perspectives are considered under the agenda of the present research investigation, as well.
As a special case of the Richardson integration, we have discussed the optimization of the Romberg integration of an arbitrary real valued function as a weighted extrapolation model. As mentioned before, the notion of the Romberg integration arises as a special case of the present consideration. Herewith, by optimizing the Richardson integration, we have optimized the Simpson 1=3 or 3=8 values of a numerical integral, or the trapezium rule as a suitable convex combination of given intervals of the integration. It is worth emphasizing that the optimization of the Richardson integration equally optimizes various limiting numerical integrations, as well. Moreover, we have given qualitative analysis of the Richardson integration and its optimization with its relevance to existing theoretical and experimental models. As a future scope of this research, we have highlighted perspectives from the sensitivity analysis, finite difference schemes, non-uniform grids, forecasting of the weather as a time integration and the original Richardson's dream toward the atmospheric circulation as well as discrete dynamical systems.
In the light of numerical practices, the significance of the results is that our analysis gives polynomial type stability criteria for the Richardson integration of an arbitrary integrable function under variations of its step sizes. In this discussion, a comparative analysis of the proposed optimization technique with the other existing quadrature methods and practical practices are left open for a future study. On the other hand, we considered the number of steps ðb À aÞ=h j while performing the Richardson integration in the light of the randomized algorithms of Agrawal and Biswas. Hereby, it follows that both step sizes h 1 and h 2 are randomized and they give a correct partition of the integration interval ða; bÞ in the derandomization limit. Investigations of such processes to obtain a correctly derandomized optimization task with integer valued steps as the optimization constraints are among the future research scope.
Finally, it is worth mentioning that our analysis can be applied to both the classical numerical algorithms that are solutions of a given system of differential equations, and to other advanced numerical methods in the realm of different permutations, combinations and splitting procedures, e.g., label switching, multi-modality and population structures (Jakobsson & Rosenberg, 2007). From the perspective of the present analysis, classical numerical algorithms are devised and supplemented either as the corresponding supervised, semi-supervised or unsupervised solutions (Huang, Song, Gupta, & Cheng, 2014). Hereby, the following two issues are important: (a) a large time-step and two small time-steps need to be handled by using the same numerical method, and (b) the order of the selected numerical integration should be given prior to an initiation of the algorithm, see (Bousquet & Elisseeff, 2002) towards the machine learning and stable Richardson techniques. Industrial applications of the proposed study are the subject matter of future research and developments.