Numerical Integration as an Initial Value Problem

Numerical integration (NI) packages commonly used in scientific research are limited to returning the value of a definite integral at the upper integration limit, also commonly referred to as numerical quadrature. These quadrature algorithms are typically of a fixed accuracy and have only limited ability to adapt to the application. In this article, we will present a highly adaptive algorithm that not only can efficiently compute definite integrals encountered in physical problems but also can be applied to other problems such as indefinite integrals, integral equations and linear and non-linear eigenvalue problems. More specifically, a finite element based algorithm is presented that numerically solves first order ordinary differential equations (ODE) by propagating the solution function from a given initial value (lower integration value). The algorithm incorporates powerful techniques including, adaptive step size choice of elements, local error checking and enforces continuity of both the integral and the integrand across consecutive elements.



page 1

page 2

page 3

page 4


Least-Squares Finite Element Method for Ordinary Differential Equations

We consider the least-squares finite element method (lsfem) for systems ...

Comparison of Numerical Solvers for Differential Equations for Holonomic Gradient Method in Statistics

Definite integrals with parameters of holonomic functions satisfy holono...

Positivity-Preserving Adaptive Runge-Kutta Methods

Many important differential equations model quantities whose value must ...

Fenrir: Physics-Enhanced Regression for Initial Value Problems

We show how probabilistic numerics can be used to convert an initial val...

Variationally consistent mass scaling for explicit time-integration schemes of lower- and higher-order finite element methods

In this paper, we propose a variationally consistent technique for decre...

Parareal Algorithm Implementation and Simulation in Julia

We present a full implementation of the parareal algorithm---an integrat...

Linear multistep methods and global Richardson extrapolation

In this work, we study the application the classical Richardson extrapol...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 first-order 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


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


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.


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


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.


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


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.


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 ,


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


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


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.

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

  2. Calculate the solution coefficients by back substitution.

  3. Calculate the error between the two sides of eq. (2) at the end point using eq. (11) for the left–hand side and by evaluating directly. Retain the value of which will be needed at the beginning of the next element.

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

  5. 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 .

  6. 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.


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 111 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


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
Table 1: Results for closed integrals taken from doi:10.1080/10586458.2005.10128931 . The values of the integral from the present method are shown in the second column. and of both the present method and DQAG are also shown.

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.


The integral of the above function is


with . Both functions and are plotted in Fig. 1.

Figure 1: (Color online) Plots of the functions given in Eqs. (14) and (15). varies very rapidly towards the upper integration limit causing numerical difficulty.

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.

Figure 2: (Color online) Step size of the finite elements (blue .) and of the integral at the end of the elements (red *) are shown. The algorithm spends more time and picks up more errors towards the upper integration limit. Vertical axis values are of the actual.

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 1751-8121-47-23-235204 .


Since this is a nonlinear equation, it has to be solved in an iterative fashion as:


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 .

Figure 3: (Color online) Calculated solution functions to eq. (16) for with initial condition are shown. The functions are oscillatory near the origin before they evolve into an asymptotic discrete bundle.

Similar plots have been reported in ref. 1751-8121-47-23-235204 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
Table 2: Results for to eq. (16) for the respective initial conditions . The rest of the columns are: total number of steps, average step size, average number of functions evaluations per element and time elapsed.

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


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 .


Hence, eq. (18) can be written as


where, now represents the inner integral given below which needs to be propagated only once from the origin until the integral converges.


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.

Table 3 shows a sample calculation for the integral in eq. (20) for .

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]
Table 3: Exact and calculated values of the integral are shown for & . Columns 3 and 4 show the number of function evaluations taken in the integrals Eq. (20) and Eq. (19) respectively. Numbers in square bracket signify powers of ten.

The rest of the parameters are set as . The Bessel functions and have been computed using the subroutines in GNU Scientific Library (GSL) 222 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) DE-NA0002630. CAW was also supported in part by the Defense Threat Reduction Agency.

6 References