1 Introduction
Numerical integration(NI) is one of the most useful numerical tools that is routinely utilized in all scientific and engineering applications. There are general and specialized methods that efficiently compute definite integrals even for many pathological cases. Most of the integration algorithms are in software packages and are available in almost all scientific libraries. Ref. davis2007methods is an excellent book on the subject containing a comprehensive reference to scientific articles and the accompanying software. However, as efficient and advanced as those quadrature algorithms are, and perhaps for this very reason, the techniques implemented in those algorithms do not necessarily carry over to treat other related numerical problems.
Applications pertinent to NI are usually a lot more complex than simply a value at a single point at the upper integration limit. A broader and versatile advantage can be gained if integration is viewed from its general mathematical perspective as a specific case of an initial value problem. In this paper, we present an algorithm that solves ODEs that, with a slight modification, can be used on many relevant applications, one of the most important of which being NI. As far as numerical quadrature is concerned, the efficiency of the method given here is comparable to the common NI packages available in scientific libraries. Moreover, the present method is so simple to use that it can be readily modified and applied to a broader spectrum of numerical applications as long as they can be set up as initial value problems.
The algorithm in question has recently been applied to a second order ODE to solve and treat a challenging system–namely, the soft Coulomb problem, where its versatility and power is self–evident PhysRevE.89.053319 . In this work, we discuss in some detail how the same algorithm can be adapted to solve a firstorder ODE and thereby apply it, in a purely mathematical setting, to the calculation of the numerical integral of a given function. Illustrative examples that are best solved by the present method, and thereby highlight its important features, will also be included.
2 Description of Algorithm
The main intent of this paper is to present a finite element algorithm that numerically approximates a solution function to the following first order ODE
(1) 
with a given initial condition , where, both and are assumed to be finite. Generally, can be a simple function of such that the above equation may be a linear or non–linear eigenvalue problem. can also be a kernel function of homogeneous or inhomogeneous integral equation kythe2002computational . Particularly, if is a predetermined simple function and , then the above equation will be equivalent to a single integral given by
(2) 
Begin the numerical solution to eq. (1) by breaking the x–axis into finite elements and mapping into a local variable with domain
defined by a linear transformation given below.
(3) 
Here, labels the elements with , while is half the size of the element. In terms of the local variable eq. (1) can be re–written as
(4) 
The over–bar indicates the appropriate change in functional form while the element index is dropped for simplification of notation. At this point, we will expand in a polynomial basis set as follows.
(5) 
Notice that and similarly for . The main import of the above expansion is that it allows us to enforce continuity of both and across the boundary of two consecutive elements. This is because the basis functions and their derivatives (denoted by ) identically vanish at . These functions and are defined in terms of Legendre polynomials of the first kind citeulike:1816367 as
(6) 
They are discussed in PhysRevE.89.053319 in more detail. These polynomials satisfy the following recurrence identities presented here for the first time. It is surprising to see how these relations remain three–term, with no surface values, despite the fact that the polynomials and are sequentially primitives of Legendre polynomials.
(7)  
(8) 
The above two relations allow us to employ the Clenshaw recurrence formula Press:1992:NRF:141273 which is known to facilitate effective numerical evaluation of relevant summations as the one included in eq. (5).
Substituting the expansion given in eq. (5) into eq. (4), evaluating the resulting equation at Gauss–Legendre (GL) abscissas and after some rearrangement, we get the following set of simultaneous equations of size ,
(9) 
where, is a root of an order Legendre polynomial. This technique, known as collocation method, is an alternate way of constructing a linear system of equations compared to the more familiar projection integrals. are now elements of a square matrix which, for a given size , is constant and hence, merely needs to be constructed only once, LU decomposed and stored for all times. Thus solving for the unknown coefficients only involves back substitution. Construction of the right–hand side column, on the other hand, requires an evaluation of the function at points including at the beginning of the element. In the examples that follow we will denote the total number of function evaluations as .
After the coefficients are calculated this way, the value of the integral at the end point of the element can be obtained by evaluating eq. (5) at . The result is simply
(10) 
One of the advantages of representing the solution in the form given in eq. (5) is that it allows us to evaluate the value of the integral at any continuous point inside the element. In this sense, is simply a function whose domain extends into all of the solved elements as the propagation proceeds. Hence, from a numerical perspective, this algorithm is best implemented by way of object–oriented programming in order to incorporate and preserve all the necessary quantities of all the elements in a hierarchy of derived variables. The integrand function , for instance, can be saved (as an object) for later use by including, among other quantities, the grid containing the coordinates of the steps . Then locating the index to which a given point
belongs, is an interpolation exercise for which efficient codes already exist. See, for instance, subroutine
hunt and its accompanying notes in ref. Press:1992:NRF:141273 . For a sorted array, as is our case, this search generally takes about tries while further points within close proximity can be located rather easily.An even more important advantage of this form of the solution is its suitability for estimation of the size of the next element via the method described in
PhysRevE.89.053319 . This adaptive step size choice relies on the knowledge of the derivatives of up to order at the end of a previous element. Those derivatives can be directly computed from eq. (5). Instead of solving two quadratic equations as we did in the last paper, however, we will solve one cubic equation, which is more suited for the present purpose. The user is required to guess only the size of the very first element. Overestimation of the size of an element may cause an increase in the number of function evaluations, but there will be no compromise in the accuracy of the solution as it will be made clear in a moment.
One other attractive feature is that we can precisely measure the error of the calculated solution directly from eq. (4). Specifically, at the end of a solved element , the error between the exact and the calculated
(11) 
can be obtained. Notice that will be needed in the beginning point of the next element. If the resulting error is not satisfactory, a bisection step will be taken to reduce the size of the element by moving the upper limit closer to and re–solving. The main purpose of the adaptive step size choice is then to estimate a priory an optimum step size, thereby reducing the number of bisections and/or the number of function evaluations necessary. Since this error is that of the integrand , and not of the integral , it is not necessary to demand this error be as small as the machine precision. The reason is, higher derivatives of computed from eq. (5) will successively deteriorate in accuracy as the order increases. In a digit calculation a relative error of in often suffices, depending on how smooth the integrand is, for calculating the solution function correctly to within . Other quadrature methods do not directly take the integrand into account in their error estimation. They often use a formula, an outcome usually of a non–trivial analytic derivation, that estimates the upper bound of a residual term for a specific order Ehrich1999 .
In the present algorithm, the number of basis functions is kept constant. Other methods such as Clenshaw–Curtis quadrature, take advantage of a convenient property of the roots of Chebyshev polynomials that allows for conserving preceding function evaluations whenever the order of expansion is doubled. But this recursive doubling of basis set size is not necessarily effective as far as improving the accuracy of the integration is concerned. The reason is, for a given working precision, there is usually only a limited range of an optimum number of basis functions, say , whose half or double is either too small or excessive. Rather, a more effective way is to estimate an adaptive step size, fix an optimum order of expansion, and reduce the size of the element whenever necessary  which is what is done here.
Implementation of the algorithm starts by fixing the size of the first step and the size of the linear system . Also, define . The rest of the procedure is sketched below.

Construct the right–hand side of eq. (9) by evaluating the function at the GL nodes.

Calculate the solution coefficients by back substitution.

If the error is not satisfactory, reduce and go to step (1).

Otherwise calculate as per eq. (10) and its higher derivatives at using eq. (5) and estimate the size of the next step using the method presented in PhysRevE.89.053319 .

For finite domain systems with an upper limit
, make sure the resulting step size does not stride beyond
. i.e., take . For open integrals on the other hand, keep propagating until the value of the integral converges.
Finally, we may sometimes prefer not to evaluate the function at the lower limit of the integration. A good example is an integrand containing singular terms at the origin. In those instances, only for the very first element, the main expansion given in eq. (5) can be altered to be in terms of the polynomials as follows.
(12) 
The rest of the propagation can then resume normally and the other modifications that follow can be worked out straightforwardly.
3 Numerical Examples
We will now consider illustrative examples that demonstrate the typical behavior of the numerical algorithm discussed above. All of the examples chosen are difficult to calculate with an accuracy close to working precision, without breaking the range of integration into smaller intervals. Whenever applicable, comparisons will be made with DQAG Quadpackref , which is one of the integration subroutines compiled in QUADPACK ^{1}^{1}1http://www.netlib.org/quadpack. DQAG is an adaptive method that keeps bisecting the element with highest error estimate until the value of the overall integration is achieved to within the desired error.
In all of the examples that follow, the error at the end of the elements has been determined by
(13) 
where, and denote the absolute and relative errors respectively. The values of and will be used unless otherwise specified. The size of the first step is fixed to be 0.5–that is and the number of basis functions used is .
DQAG, on the other hand, has been run with absolute and relative errors of and the lowest order () has been used in order to keep the number of function evaluations at a minimum so as to provide a fair comparison. Our program is written in modern C++ and run on a late 2013, 2.8 GHz Intel Core i7 MacBook Pro laptop computer with Apple LLVM 8.1 compiler.
3.1 Closed Integral
We consider the closed integrals studied in doi:10.1080/10586458.2005.10128931 . These integrals exhibit different properties and were chosen by Bailey et al. to test three other methods of NI. Herein we attempt all of them and report the output. Table 1 shows the results for numerical values of those integrals calculated using the present method as well as DQAG.
label  (DQAG)  (DQAG)  

0.250 000 000 000 000  1.110[16]  0.000  29  15  
0.210 657 251 225 807  0.000  1.318[16]  29  45  
1.905 238 690 482 68  0.000  0.000  191  15  
0.514 041 895 890 071  0.000  0.000  29  45  
0.444 444 444 444 445  8.743[16]  3.747[16]  871  885  
0.785 398 163 397 448  1.414[16]  2.827[16]  974  795  
1.198 140 227 142 81  6.337[9]  6.159[9]  2129  1725  
1.999 999 999 999 98  8.438[15]  4.441[16]  922  1545  
1.088 793 045 151 79  9.993[15]  2.855[15]  1243  1335  
2.221 441 454 672 65  6.485[9]  2.907[9]  2032  1725  
1.570 796 326 794 90  0.000  0.000  29  105  
1.772 453 840 168 93  6.057[9]  5.929[9]  2439  1455  
1.253 314 137 315 62  9.514[14]  0.000  96  255  
0.500 000 000 000 001  1.110[15]  0.000  231  375  
1.570 796 326 794 90  1.414[16]  2.218[11]  1523  210 
with respect to the exact values, and a number of function evaluations are also indicated. and of the output from DQAG are also shown for comparison. The two methods are essentially, qualitatively similar in performance. Notice that the labels of the integrals are taken from the article doi:10.1080/10586458.2005.10128931 .
Lets take a closer look at one of the integrals (index 6) in order to demonstrate the behavior of our algorithm. Example 6 has an integration limits and an integrand given below.
(14) 
The integral of the above function is
(15) 
with . Both functions and are plotted in Fig. 1.
The calculated value of the integral is which has a of compared to the exact value. It took function evaluations and steps to propagate the integral from to , where the size of the first step was set to as mentioned above. As can be seen from Fig. (2), the step sizes chosen by our algorithm kept decreasing to as low as , which is consistent with the singularity of the derivative of the integrand at the upper integration limit.
Fig. (2) also shows relative errors of the integral at the end of all the elements in comparison with the exact values from eq. (15).
Clearly, the maximum error occurs at the last step near , which is what motivated our choice of this example. Rapidly changing functions, such as those with integrable singularity, are generally troublesome to the algorithm because the step size may not be small enough and/or the order of polynomials high enough to accommodate portions of the integrand with (nearly) vertical shape.
In comparison, DQAG takes function evaluations and calculates the integral with a of .
3.2 Nonlinear Problem
Bender et al. have studied the following interesting nonlinear eigenvalue problem for which the numerical solution can be very challenging 175181214723235204 .
(16) 
Since this is a nonlinear equation, it has to be solved in an iterative fashion as:
(17) 
where and labels the levels of the iteration. Notice that at any stage of the iteration only the values of at the GL nodes are required. The first iterate of array
is seeded from the solution vector
of the final result in the preceding element. This is not only convenient but also an excellent approximation because the adaptive step choice implemented here is based on the assumption that the solution function between two consecutive elements remains constant up to a fourth order Taylor series expansion. The first element has been started by setting all the elements of the solution vector to unity. Apparently, only for this problem, the first two steps of the procedure given in Section 2, must be repeated until the iteration in eq. (17) converges.The results of our calculation for , are plotted in fig. 3 for initial values at the origin .
Similar plots have been reported in ref. 175181214723235204 which are a result of point–wise convergent calculations. Our solution function, on the other hand, has been propagated from the origin outward. The very large number of steps the solution required for such a modest distance of is quite remarkable. Table 1 summarizes the output of our program.
0.020 844 865 419 015 3  1 633 376  1.469  254.338  63  
0.104 224 327 270 128  1 701 378  1.411  253.336  66  
0.270 983 253 633 302  1 908 989  1.257  251.829  73  
0.437 742 187 280 245  2 068 413  1.160  250.609  79  
0.687 880 611 222 152  2 397 070  1.001  249.380  91  
0.938 019 076 811 230  2 636 633  0.9103  248.157  100  
1.271 537 122 002 93  3 041 709  0.7890  247.071  116  
1.688 434 875 810 57  3 586 933  0.6691  245.951  135  
2.105 332 915 403 23  4 021 621  0.5968  244.952  151  
2.605 611 041 676 66  4 626 563  0.5187  244.110  174 
It took millions of steps, with an average step size as low as and an average number of function evaluations up to per element. Notice that this includes all the iterations, which, along with the brevity of the elapsed times shown, indicates the efficiency of our implementation.
We have also included the final numerical values at for reference. For this example, the required error has been set lower as , for an obvious reason. In order to check the validity our results, we have re–run the program by still lowering the magnitude of the error to . With the lowered error, the indicated values of vary by an amount no larger than . The computational times in the last column of the table also increased to as high as seconds. Maintaining this much accuracy after millions of steps, and in an iterative calculation, shows how robust our method is. This non–linear eigenvalue problem is computationally challenging indeed.
3.3 Double–Range Integrals
In this example, we consider a double integral for which one of the limits of the inner integral is identical to the variable of the outer integral. These types of integrals are common in studies of many particle dynamical systems involving Green’s functions or double range addition theorems such as the Laplace expansion. Particularly, we will look at the following integral, which is the most prominent radial integral that is encountered in calculations that involve exponential type orbitals or Geminals in a spherical coordinate system. It stems from the addition theorem for given in QUA:QUA24319 and the resulting integrals are still topics of interest in recent research Jiao2015140 ; Rico2012 ; QUA:QUA21002 . Let us define the integral as
(18) 
where and are spherical modified Bessel functions of the first and second kind respectively citeulike:1816367 . The screening parameters are positive real numbers while all of the indices are integers. The correct composition of these parameters is such that both the inner and outer integrals remain finite as is the case in physical applications.
In order to propagate from the origin, both lower integration limits need to be set to zero. This can be attained by switching the order of the two integrals using the following identity which maintains identical – region of integration hassani2000mathematical .
(19) 
Hence, eq. (18) can be written as
(20) 
where, now represents the inner integral given below which needs to be propagated only once from the origin until the integral converges.
(21) 
Here signifies the end point of the last element where the result of the above integral converged to its value at infinity. Hence, beyond this point is considered to be constant function, i.e., for . Once the above integral is done and all the relevant parameters stored, it can be evaluated at any desired point , which is a significant computational gain since the double integral given in eq. (18) has essentially been reduced to two simple integrals. This demonstrates one of the main advantages contained in the present algorithm.
N of  N of  I – calculated  I – exact  

0.5  219  259  1.627 473 168 386 65 [27]  1.627 473 168 386 653 87 [27]  
1.0  219  218  2.559 085 779 949 79 [22]  2.559 085 779 949 794 01 [22]  
2.0  219  231  3.103 777 873 917 21 [17]  3.103 777 873 917 210 86 [17]  
0.5  232  259  2.946 389 365 576 82 [23]  2.946 389 365 576 741 23 [23]  
1.0  232  245  6.062 810 005 197 87 [18]  6.062 810 005 197 874 73 [18]  
2.0  232  204  9.533 337 428 978 80 [13]  9.533 337 428 978 808 27 [13]  
0.5  245  259  4.342 544 722 241 59 [19]  4.342 544 722 241 718 83 [19]  
1.0  245  245  1.097 615 571 907 40 [15]  1.097 615 571 907 438 80 [15]  
2.0  245  231  2.258 572 729 378 12 [10]  2.258 572 729 378 146 95 [10] 
The rest of the parameters are set as . The Bessel functions and have been computed using the subroutines in GNU Scientific Library (GSL) ^{2}^{2}2http://www.gnu.org/software/gsl/. For the first element, eq. (12) has been used in order to avoid evaluation at the origin. We also have calculated the exact values of the integrals using Mathematica Mathematica7 , the first 18 digits of which are displayed in the last column. Comparison with the the present calculated results reveals that the integral was done accurately. The table also further shows the total number of function evaluations for the inner and outer integrals.
4 Conclusion
There are many physical applications that can be modeled as initial value problems and the methods available to compute them are equally diverse. No single algorithm is known to address all of them at once, and hence, researchers usually digress from their main area of interest so as to familiarize themselves with many of the computational options. The present algorithm is by no means capable of solving all initial value problems, but it comes pragmatically close, especially for most physical applications. This is possible because it produces solutions on finite elements whose size is chosen to locally, checks the validity of the solution, and communicates the solution function and its first derivative to the next element, maintaining continiuity of both. More importantly, it can be applied on any ODE because it is easy to implement and modify, especially with the proper use of object–oriented programming techniques. The basis functions and discussed above are based on Legendre polynomials. In the future, we will use other classic orthogonal polynomials such as Chebyshev or Jacobi to determine if further advantages can be gained.
5 Acknowledgement
DHG and CAW were partially supported by the Department of Energy, National Nuclear Security Administration, under Award Number(s) DENA0002630. CAW was also supported in part by the Defense Threat Reduction Agency.
6 References
References

(1)
P. Davis, P. Rabinowitz,
Methods of Numerical
Integration, Dover Books on Mathematics Series, Dover Publications, 2007.
URL https://books.google.com/books?id=gGCKdqka0HAC 
(2)
D. H. Gebremedhin, C. A. Weatherford,
Calculations for
the onedimensional soft coulomb problem and the hard coulomb limit, Phys.
Rev. E 89 (2014) 053319.
doi:10.1103/PhysRevE.89.053319.
URL http://link.aps.org/doi/10.1103/PhysRevE.89.053319 
(3)
P. Kythe, P. Puri,
Computational Methods
for Linear Integral Equations, Birkhäuser Boston, 2002.
URL https://books.google.com/books?id=CmX1aHur7jcC  (4) G. B. Arfken, H. J. Weber, F. E. Harris, Mathematical Methods for Physicists, Sixth Edition: A Comprehensive Guide, 6th Edition, Academic Press, 2005.
 (5) W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in FORTRAN (2Nd Ed.): The Art of Scientific Computing, Cambridge University Press, New York, NY, USA, 1992.

(6)
S. Ehrich, Applications
and Computation of Orthogonal Polynomials: Conference at the Mathematical
Research Institute Oberwolfach, Germany March 22–28, 1998, Birkhäuser
Basel, Basel, 1999, Ch. Stieltjes Polynomials and the Error of GaussKronrod
Quadrature Formulas, pp. 57–77.
doi:10.1007/9783034886857_4.
URL http://dx.doi.org/10.1007/9783034886857_4 
(7)
R. Piessens, E. d. DonckerKapenga, C. Überhuber, D. Kahaner,
Quadpack, Springer
Series in Computational Mathematics, SpringerVerlag Berlin Heidelberg, 1983.
URL http://www.springer.com/us/book/9783540125532 
(8)
D. H. Bailey, K. Jeyabalan, X. S. Li,
A comparison of three
highprecision quadrature schemes, Experimental Mathematics 14 (3) (2005)
317–329.
arXiv:http://dx.doi.org/10.1080/10586458.2005.10128931, doi:10.1080/10586458.2005.10128931.
URL http://dx.doi.org/10.1080/10586458.2005.10128931 
(9)
C. M. Bender, A. Fring, J. Komijani,
Nonlinear eigenvalue
problems, Journal of Physics A: Mathematical and Theoretical 47 (23) (2014)
235204.
URL http://stacks.iop.org/17518121/47/i=23/a=235204 
(10)
D. Gebremedhin, C. Weatherford,
Canonical tworange addition
theorem for slatertype orbitals, International Journal of Quantum Chemistry
113 (1) (2013) 71–75.
doi:10.1002/qua.24319.
URL http://dx.doi.org/10.1002/qua.24319 
(11)
L. G. Jiao, Y. K. Ho,
Computation
of twoelectron screened coulomb potential integrals in hylleraas basis
sets, Computer Physics Communications 188 (2015) 140 – 147.
doi:http://dx.doi.org/10.1016/j.cpc.2014.11.019.
URL http://www.sciencedirect.com/science/article/pii/S0010465514004020 
(12)
J. F. Rico, R. López, G. Ramírez, I. Ema,
Repulsion integrals
involving slatertype functions and yukawa potential, Theoretical Chemistry
Accounts 132 (1) (2012) 1–9.
doi:10.1007/s002140121304x.
URL http://dx.doi.org/10.1007/s002140121304x 
(13)
J. J. Fern ndez, R. L pez, I. Ema, G. Ram rez, J. Fern ndez Rico,
Auxiliary functions for molecular
integrals with slatertype orbitals. i. translation methods, International
Journal of Quantum Chemistry 106 (9) (2006) 1986–1997.
doi:10.1002/qua.21002.
URL http://dx.doi.org/10.1002/qua.21002  (14) S. Hassani, Mathematical Methods: For Students of Physics and Related Fields, Undergraduate Texts in Contemporary Physics, Springer New York, 2000.
 (15) Mathematica, Version 7.0 1.0, Wolfram Research, Champaign, Illinois, 2009.
Comments
There are no comments yet.