Automatic Numerical Integration Research Articles

Highly accurate, efficient, and automatic computation of reversible cyclic voltammograms, using double exponential formulas for numerical integration

Highly accurate reference values of theoretical cyclic voltammograms are indispensable for those who develop and examine simulation methods or other modelling approaches for electroanalytical chemistry. Reference values having relative errors close to ± 10−16 (which is the level of the representation errors of standard double precision variables) are desirable, but they cannot be obtained efficiently using conventional simulations. A class of reversible cyclic voltammograms can be expressed in terms of improper convolution integrals that have to be evaluated numerically. By employing the efficient automatic numerical integration package INTDE of Ooura, based on double exponential formulas, and written in “C”, it is shown that such voltammograms can be calculated with relative errors close to ± 10−16, in a time of ca. 1ms per single current function value, on a contemporary laptop computer. Furthermore, this method appears resistant to failures that may accompany highly accurate numerical integration by alternative methods and software. This is demonstrated for three kinetic examples under conditions of pure diffusion transport: reversible charge transfer at a planar electrode; reversible charge transfer at a spherical electrode; and ErevCirr’ mechanism at a planar electrode.

Numerical integration and extrapolation for finite and UV-divergent 3-loop Feynman integrals

Automatic numerical integration attempts to alleviate the detailed application of solution methods to the problem. We use automatic multivariate integration with PARINT combined with numerical extrapolation for classes of 3-loop UV-divergent self-energy diagrams with massive internal lines, and finite or UV-divergent massless vertex diagrams. The extrapolation is performed with respect to the dimensional regularization parameter. It is shown that the PARINT software handles the integrations accurately and efficiently.

Feynman loop numerical integral expansions for 3-loop vertex diagrams

The inclusion of higher order corrections by multi-loop diagrams is required in perturbation methods for precise theoretical calculations of high energy interactions. We address 3-loop vertex Feynman diagrams with massless internal lines, and which may exhibit ultra-violet (UV) singularities. The computational methods target automatic numerical integration and extrapolation to approximate the leading coefficients of the integral expansion with respect to the dimensional regularization parameter. Convergence acceleration is achieved by applying linear extrapolation on the integral expansion. Multivariate integration is performed with the ParInt software package, layered over MPI (Message Passing Interface) to speed up the computations. Integrand transformations result in diminishing the effect of singular behavior in the integrand.

Stress Intensity Factor Equations for Internal Semi-Elliptical Cracks in Pressurized Cylinders

In order to calculate the fatigue life of cylinders subjected to internal pressure using the fracture mechanics approach, the stress intensity factors (SIFs) for internal semi-elliptical cracks are needed. Nowadays, the most accurate procedure for fatigue life calculation consists in starting from a postulated internal semi-elliptical crack and updating the flaw aspect ratio during the crack propagation. In this sense, assuming a semi-elliptical shape during crack propagation, SIFs both at the deepest crack point and at the surface point must be calculated in order to update the crack aspect ratio through its fatigue propagation. This continuous crack shape updating cannot be done using the conventional tabulated solutions, as those provided in the main design codes and SIF handbooks. This paper presents a number of closed-form equations, which very accurately fit the tabulated results in ASME Boiler and Pressure Vessel Code, Section VIII, Division 3, to calculate the SIF for internal semi-elliptical cracks in cylinders subjected to internal pressure. These equations can be used to avoid the use of the tabulated solutions, and, as a consequence, an automatic numerical integration of the propagation law can be done.

Extrapolation and Adaptivity in Software for Automatic Numerical Integration on a Cube

Asymptotic expansions related to the integration of well-behaved functions on simplices and cubes have been known for several decades. Extensions of these results to classes of vertex and line singularities are also known. The nature of these expansions justifies extrapolation using the e-algorithm of Wynn. In principle this requires a uniform subdivision of the region. This was implemented in QUADPACK for finite intervals and in TRIEX for triangles about 15 years ago. In this paper its incorporation in CUBPACK, a software package for automatic integration over a collection of cubes and simplices, is described and some results are reported. We also report on a special subdivision strategy that offers an alternative approach for higher-dimensional problems.

A Numerical Method for the Analysis of a Lifting Body in the Potential Flow

A numerical method is proposed for analyzing the potential flow around a lifting body such as a wing, a rudder, and a propeller. The method is different from the panel method in which boundary surface of the lifting body is divided into panels. In the present method the boundary surface is treated as continuous surface and the unknown (potential on the surface) is also expressed as the continuous function. So we call the method “continuous method”. In the method the boundary integral equation is solved as accurately as necessary for the purpose using automatic numerical integration routine.Numerical results by the present method are compared with the analytical solution in case of a sphere, an ellipsoid, and a Kármán-Trefftz aerofoil. The comparison shows the accuracy of the method. Finally numerical results are shown for a circular wing and the results are compared with other methods such as the panel method and the lifting surface method, and also compared with experimental results. Those comparisons show the validity of the present method.

Weighted Convergence of Lagrange Interpolation Based on Gauss-Kronrod Nodes

The Gauss-Kronrod quadrature scheme, which is based on the zeros of Legendrepolynomials and Stieltjes polynomials, is a standard rule for automaticnumerical integration in mathematical software libraries. For a long time,very little was known about the underlying Lagrange interpolationprocesses. Recently, the authors proved new bounds and asymptoticproperties for the Stieltjes polynomials and, subsequently, appliedthese results to investigate the associated interpolation processes. Thepurpose of this paper is to survey the quality of these interpolationprocesses, with additional results that extend and complete the existingones. The principal new results in this paper are necessary and sufficientconditions for weighted convergence. In particular, we show that theLagrange interpolation polynomials associated with the above interpolationprocesses have the same speed of convergence as the polynomials of bestapproximation in certain weighted Besov spaces.

Automatic numerical integration techniques for polyatomic molecules

We describe a new algorithm for the generation of 3D grids for the numerical evaluation of multicenter molecular integrals in density functional theory. First, we use the nuclear weight functions method of Becke [A. D. Becke, J. Chem. Phys. 88, 2547 (1988)] to decompose a multicenter integral ∫F(r) dr into a sum of atomic-like single-center integrals. Then, we apply automatic numerical integration techniques to evaluate each of these atomic-like integrals, so that the total integral is approximated as ∫F(r) dr≊∑iωiF(ri). The set of abscissas ri and weights ωi constitutes the 3D grid. The 3D atomic-like integrals are arranged as three successive monodimensional integrals, each of which is computed according to a recently proposed monodimensional automatic numerical integration scheme which is able to determine how many points are needed to achieve a given accuracy. When this monodimensional algorithm is applied to 3D integration, the 3D grids obtained adapt themselves to the shape of the integrand F(r), and have more points in more difficult regions. The function F(r), which, upon numerical integration, yields the 3D grid, is called the generating function of the grid. We have used promolecule densities as generating functions, and have checked that grids generated from promolecule densities are also accurate for other integrands. Our scheme is very reliable in the sense that, given a relative tolerance ε, it generates 3D grids which are able to approximate multicenter integrals with relative errors smaller than ε for all the molecules tested in this work. Coarser or finer grids can be obtained using greater or smaller tolerances. For a series of 21 molecules, the average number of points per atom for ε=2.0⋅10−3, ε=2.0⋅10−4, ε=2.0⋅10−5, ε=2.0⋅10−6, and ε=2.0⋅10−7 is respectively 3141 (2.9⋅10−4), 10271 (2.4⋅10−5), 27184 (3.1⋅10−6), 72266 (1.9⋅10−7), and 164944 (5.2⋅10−9) (in parentheses are the maximum errors obtained when integrating the density). It is possible to reduce the number of points in the grid by taking advantage of molecular symmetry. It seems that our method achieves a given accuracy with fewer points than other recently proposed methods.

A simple, efficient and more reliable scheme for automatic numerical integration

We have recently proposed a simple, reliable and efficient scheme for automatic numerical integration that uses the change of variable x(z) = 1+(2/Π){[1+ 2 3 (1−z 2)]z√1−z 2−arccos z} to transform the integral to be computed, ∫ 1-1 f( x) dx, into (16/3 Π) ∫ 1-1 f( x( z)) (1− z 2)√1− z 2 dz,, which is approximated by successive n-point Gauss-Chebyshev quadrature formulas of the second kind ( I n ). The following sequence of formulas was generated: I 1, I 3, I 7,..., I ( n−1)/2 , I n , I 2 n+1 . In the present work we generate the same sequence, together with another one, I 2, I 5,..., I m , I 2 m +1. Both sequences are generated in an alternate way, I 2, I 1, I 5, 3,..., I ( n−1)/2 , I m , I n , I 2 m+1 , where m + 1 = 2 3 (n + 1) . The main advantage of the new scheme is that, unlike the previous one, the total number of points increases more moderately than doubling. This is possible because all the abscissas of I ( n−1)/2 are also abscissas of I n , all the abscissas of I m are also abscissas of I 2 m+1 , and, especially, all the abscissas of I n are also abscissas of I 2 m+1 . These allows us to implement a simple and efficient algorithm to generate both sequences alternatively (a FORTRAN 77 version is included). We propose an error estimate for the approximate integral that is more conservative than the one used in our previous work, which results in a more reliable automatic numerical integrator, very suitable for multiple integrals.

Automatic generation and numerical integration of differential-algebraic equations of multibody dynamics

Abstract An algorithm for computing the matrix and vector functions in the numerical integration method of constrained equations of motion in multibody mechanical system dynamics is presented. Cartesian coordinates and Euler parameters have been used as the state variables to construct the mathematical model for a spatial multibody mechanical system. Four basic constraints, which describe the orthogonality and the parallelism of pairs of vectors, have also been written and will serve as building blocks of the kinematic joint library. Derivatives of these basic constraints are then derived by applying differential operators to the multivariable and vector-valued constraint functions. Verification is done by comparing the numerical solution with results computed by a symbolic computation software. In multibody dynamics, the acceleration of constrained equations of motion is a function of the state variables. Using the differentials of joint library functions and generalized force functions, the Jacobian of the acceleration with respect to the state variables is obtained. The efficiency and accuracy of computing this analytic Jacobian matrix are demonstrated by comparing execution time and numerical results with those obtained by applying a first order linearization method. The proposed computational technique is also suitable for various formulations of the equations of motion, for instance, the Euler angles instead of the Euler parameters can be used for the orientation variables. The final goal is to construct a standard library of joint constraints and force expressions for automatic generation of constrained equations of motion and to enhance the numerical integration methods for multibody dynamics.

A simple, reliable and efficient scheme for automatic numerical integration

A scheme for automatic numerical integration is presented. It uses the change of variable x=1+( 2 π ){[1+ 2 3 (1-z 2)]z√1-z 2 - arccosz} to transform the integral to be computed, ∫ 1 -1 f( x) d x, into ( 16 3π )∫ 1 -1f(x)(1-z 2)√1-z 2 dz , which is approximated by successive n-points Gauss-Chebyshev quadrature formulas of the second kind ( I n ). Due to the special nature of their abscissas and weights, a sequence of formulas I 1, I 3, I 7,…, I (n-1) 2 , I n, I 2n+1 may be generated, such that I 2 n+1 may be computed with only n + 1 new integrand evaluations, using the previous value of I n . An error estimation is proposed for I 2 n+1 , which only needs two previous values of the sequence ( I n and I (n+1) 2 ). The algorithm may be implemented by a very short program (a FORTRAN 77 version is included) that spends practically all its running time in integrand evaluations. It is compared with other methods for automatic numerical integration (trapezoidal rule, Simpson's rule, Romberg's method, an adaptive Gauss-Kronrod rule and Clenshaw-Curtis method) over a broad set of 20 functions. We conclude that the present method is very simple and reliable and is the most efficient among the methods tested here. Possible applications in density functional theory are explored.

A numerical study of three moving-grid methods for one-dimensional partial differential equations which are based on the method of lines

In recent years, several sophisticated packages based on the method of lines (MOL) have been developed for the automatic numerical integration of time-dependent problems in partial differential equations (PDEs), notably for problems in one space dimension. These packages greatly benefit from the very successful developments of automatic stiff ordinary differential equation solvers. However, from the PDE point of view, they integrate only in a semiautomatic way in the sense that they automatically adjust the time step sizes, but use just a fixed space grid, chosen a priori, for the entire calculation. For solutions possessing sharp spatial transitions that move, e.g., travelling wave fronts or emerging boundary and interior layers, a grid held fixed for the entire calculation is computationally inefficient, since for a good solution this grid often must contain a very large number of nodes. In such cases methods which attempt automatically to adjust the sizes of both the space and the time steps are likely to be more successful in efficiently resolving critical regions of high spatial and temporal activity. Methods and codes that operate this way belong to the realm of adaptive or moving-grid methods. Following the MOL approach, this paper is devoted to an evaluation and comparison, mainly based on extensive numerical tests, of three moving-grid methods for 1D problems, viz., the finite-element method of Miller and co-workers, the method published by Petzold, and a method based on ideas adopted from Dorfi and Drury. Our examination of these three methods is aimed at assessing which is the most suitable from the point of view of retaining the acknowledged features of reliability, robustness, and efficiency of the conventional MOL approach. Therefore, considerable attention is paid to the temporal performance of the methods.

Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region

A subroutine for automatic numerical integration is presented. The subroutine is developed by changing the basic rule used in an algorithm of P. van Dooren and L. de Ridder. Comparative test results are given.

An adaptive algorithm for numerical integration over an n-dimensional cube

A program is presented for automatic numerical integration over the n-dimensional cube. Interesting advantages of this program are the strong adaptivity of the algorithm combined with the use of a good basic integration rule. This is shown by some comparative tests with other programs.

Algorithm 468: algorithm for automatic numerical integration over a finite interval [D1

Algorithm 468: algorithm for automatic numerical integration over a finite interval [D1

Automatic integration of a function with a parameter

Two efficient methods for automatic numerical integration are Romberg integration and adaptive Simpson integration. For integrands of the form ƒ( x ) g ( x, α ) where α is a parameter, it is shown that Romberg's method is more efficient. A FORTRAN program shows how to achieve this greater efficiency.