Double Exponential Formula Research Articles

Computing the matrix exponential with the double exponential formula

Abstract This article considers the computation of the matrix exponential e A {{\rm{e}}}^{A} with numerical quadrature. Although several quadrature-based algorithms have been proposed, they focus on (near) Hermitian matrices. In order to deal with non-Hermitian matrices, we use another integral representation including an oscillatory term and consider applying the double exponential (DE) formula specialized to Fourier integrals. The DE formula transforms the given integral into another integral whose interval is infinite, and therefore, it is necessary to truncate the infinite interval. In this article, to utilize the DE formula, we analyze the truncation error and propose two algorithms. The first one approximates e A {{\rm{e}}}^{A} with the fixed mesh size, which is a parameter in the DE formula affecting the accuracy. The second one computes e A {{\rm{e}}}^{A} based on the first one with automatic selection of the mesh size depending on the given error tolerance.

Lightning Efficient Counterpoise Configurations for Transmission Line Grounding

Different configurations of counterpoises are often used to reduce the resistance and impulse impedance of the transmission line grounding. However, there are no general rules for an actual design. This paper analyses the performance of several configurations of counterpoises under typical lightning current impulses. The basis is a parametric analysis derived from time responses of 11,386 test cases for each configuration with parameters in broad ranges computed using a rigorous electromagnetic model based on the method of moments. The contributions of the paper are: Simple formulas for the resistance, impulse impedance, and effective length of counterpoise configurations are derived from and verified by the simulation results; a configuration with an 8-leg counterpoise bent parallel to the line is proposed, which extends the effective length and is efficient in high resistivity soil; and the method does not depend on the different lightning current impulse waveform representations based on Heidler, CIGRE, and double-exponential formulas, since the results can be easily converted. A simple procedure is provided to compare the capabilities and limitations of counterpoise configurations. New formulas and procedures make up a general and straightforward method for analyzing the optimal design, considering the cost-effective and best protection against lightning.

Double exponential quadrature for fractional diffusion

We introduce a novel discretization technique for both elliptic and parabolic fractional diffusion problems based on double exponential quadrature formulas and the Riesz–Dunford functional calculus. Compared to related schemes, the new method provides faster convergence with fewer parameters that need to be adjusted to the problem. The scheme takes advantage of any additional smoothness in the problem without requiring a-priori knowledge to tune parameters appropriately. We prove rigorous convergence results for both, the case of finite regularity data as well as for data in certain Gevrey-type classes. We confirm our findings with numerical tests.

Optimal selection formulas of mesh size and truncation numbers for the double-exponential formula

Optimal selection formulas of mesh size and truncation numbers for the double-exponential formula

A fast and accurate DE-formula algorithm to evaluate Ambartsumian-Chandrarasekhar $H$-function for isotropic scattering

In the present work, a compact, fast, and yet accurate algorithm is developed to calculate the numerical values of the Ambartsumian-Chandrasekhar’s \(H\)-function for isotropic scattering and its moments on the basis of the double exponential (DE) formula of Takahashi and Mori (RIMS, Kyoto Univ., 9:721, 1974). The main improvement made in the new method is an elimination of the iterative procedure for automatic adjustment of the step-size of integrations carried out with the DE-formula. Instead, a set of optimal values for the upper limit of integration \(T_{\text{max}}\) and the number of division points \(N_{\text{T}}\) to specify the step-size of the quadrature is predetermined for calculations of the \(H\)-function with a 15-digit accuracy, and also another for the evaluations of the moments with an accuracy of 14-digits or better. FORTRAN90 subroutines HFISCA for the \(H\)-function and HFMOMENT for the moments of arbitrary degrees are subsequently constructed (their source codes and a driver together with a sample set of output are shown in Appendix of this paper). Tables of sample calculations of the \(H\)-function and its moments of degree −1 through 6 carried out by these programs are also presented. The routines HFISCA and HFMOMENT should prove useful not only in astrophysical applications but also in other disciplines of science such as the electron transports in condensed matter and remote-sensing data analyses. A request for a copy of the Fortran 90 source code of the program can be made by writing to kawabata@rs.kagu.tus.ac.jp.

Double-exponential formula for infinite integrals of unilateral rapidly decreasing functions

Double-exponential formula for infinite integrals of unilateral rapidly decreasing functions

Numerical analysis of a self-similar turbulent flow in Bose–Einstein condensates

We study a self-similar solution of the kinetic equation describing weak wave turbulence in Bose–Einstein condensates. This solution presumably corresponds to an asymptotic behavior of a spectrum evolving from a broad class of initial data, and it features a non-equilibrium finite-time condensation of the wave spectrum n(ω) at the zero frequency ω. The self-similar solution is of the second kind, and it satisfies boundary conditions corresponding to a nonzero constant spectrum (with all its derivative being zero) at ω=0 and a power-law asymptotic n(ω)→ω−x at ω→∞x∈R+. Finding it amounts to solving a nonlinear eigenvalue problem, i.e. finding the value x* of the exponent x for which these two boundary conditions can be satisfied simultaneously. To solve this problem we develop a new high-precision algorithm based on Chebyshev approximations and double exponential formulas for evaluating the collision integral, as well as the iterative techniques for solving the integro-differential equation for the self-similar shape function. This procedures allow to achieve a solution with accuracy ≈4.7% which is realized for x*≈1.22.

Adaptive numerical integration of exponential finite elements for a phase field fracture model

Phase field models for fracture are energy-based and employ a continuous field variable, the phase field, to indicate cracks. The width of the transition zone of this field variable between damaged and intact regions is controlled by a regularization parameter. Narrow transition zones are required for a good approximation of the fracture energy which involves steep gradients of the phase field. This demands a high mesh density in finite element simulations if 4-node elements with standard bilinear shape functions are used. In order to improve the quality of the results with coarser meshes, exponential shape functions derived from the analytic solution of the 1D model are introduced for the discretization of the phase field variable. Compared to the bilinear shape functions these special shape functions allow for a better approximation of the fracture field. Unfortunately, lower-order Gauss-Legendre quadrature schemes, which are sufficiently accurate for the integration of bilinear shape functions, are not sufficient for an accurate integration of the exponential shape functions. Therefore in this work, the numerical accuracy of higher-order Gauss-Legendre formulas and a double exponential formula for numerical integration is analyzed.

Computing the matrix fractional power with the double exponential formula

Two quadrature-based algorithms for computing the matrix fractional power $A^\alpha$ are presented in this paper. These algorithms are based on the double exponential (DE) formula, which is well-known for its effectiveness in computing improper integrals as well as in treating nearly arbitrary endpoint singularities. The DE formula transforms a given integral into another integral that is suited for the trapezoidal rule; in this process, the integral interval is transformed into an infinite interval. Therefore, it is necessary to truncate the infinite interval to an appropriate finite interval. In this paper, a truncation method, which is based on a truncation error analysis specialized to the computation of $A^\alpha$, is proposed. Then, two algorithms are presented---one where $A^\alpha$ is computed with a fixed number of abscissa points and one with $A^\alpha$ computed adaptively. Subsequently, the convergence rate of the DE formula for Hermitian positive definite matrices is analyzed. The convergence rate analysis shows that the DE formula converges faster than Gaussian quadrature when $A$ is ill-conditioned and $\alpha$ is a non-unit fraction. Numerical results show that our algorithms achieve the required accuracy and are faster than other algorithms in several situations.

Computing the matrix fractional power with the double exponential formula

Two quadrature-based algorithms for computing the matrix fractional power $A^\alpha$ are presented in this paper. These algorithms are based on the double exponential (DE) formula, which is well-known for its effectiveness in computing improper integrals as well as in treating nearly arbitrary endpoint singularities. The DE formula transforms a given integral into another integral that is suited for the trapezoidal rule; in this process, the integral interval is transformed to the infinite interval. Therefore, it is necessary to truncate the infinite interval into an appropriate finite interval. In this paper, a truncation method, which is based on a truncation error analysis specialized to the computation of $A^\alpha$, is proposed. Then, two algorithms are presented -- one computes $A^\alpha$ with a fixed number of abscissas, and the other computes $A^\alpha$ adaptively. Subsequently, the convergence rate of the DE formula for Hermitian positive definite matrices is analyzed. The convergence rate analysis shows that the DE formula converges faster than the Gaussian quadrature when $A$ is ill-conditioned and $\alpha$ is a non-unit fraction. Numerical results show that our algorithms achieved the required accuracy and were faster than other algorithms in several situations.

Numerical Computation of Matrix Sign Function by Double Exponential Formula

Numerical Computation of Matrix Sign Function by Double Exponential Formula

Linear Transport in Porous Media

The linear transport theory is developed to describe the time dependence of the number density of tracer particles in porous media. The advection is taken into account. The transport equation is numerically solved by the analytical discrete ordinates method. For the inverse Laplace transform, the double-exponential formula is employed. In this paper, we consider the travel distance of tracer particles whereas the half-space geometry was assumed in our previous paper [Amagai et al. (2020). Trans. Porous Media 132:311–331].

Kernel representation of long-wave dynamics on a uniform slope.

Long-wave propagation on a uniformly sloping beach is formulated as a transient-response problem, with initially stationary water subjected to an incident wave. Both the water surface elevation and the horizontal flow velocity on the slope can be represented as convolutions of the rate of displacement of the water surface at the toe of the slope with a singular kernel function of time and space. The kernel, which is typically expressed in the form of an infinite series, accommodates the dynamic processes of long waves, such as shoaling, reflection, and multiple reflections over the slope and yields exact solutions of the linear shallow water equations for any smooth incident wave. The kernel convolution can be implemented numerically by using double exponential formulas to avoid the kernel singularity. The kernel formulation can be extended readily to nonlinear dynamics via the hodograph transform, which in turn enables the instantaneous prediction of nonlinear wave properties and of the occurrence of wave breaking in the near-shore area. This general description of long-wave dynamics provides new insights into the long-studied problem.

Improvement of selection formulas of mesh size and truncation numbers for the double-exponential formula

Improvement of selection formulas of mesh size and truncation numbers for the double-exponential formula

Construction of conformal maps based on the locations of singularities for improving the double exponential formula

Abstract The double exponential formula, or DE formula, is a high-precision integration formula using a change of variables called a DE transformation; it has the disadvantage that it is sensitive to singularities of an integrand near the real axis. To overcome this disadvantage, Slevinsky & Olver (2015, On the use of conformal maps for the acceleration of convergence of the trapezoidal rule and Sinc numerical methods. SIAM J. Sci. Comput., 37, A676–A700) attempted to improve the formula by constructing conformal maps based on the locations of singularities. Based on their ideas, we construct a new transformation formula. Our method employs special types of the Schwarz–Christoffel transformation for which we can derive the explicit form. The new transformation formula can be regarded as a generalization of DE transformations. We confirm its effectiveness by numerical experiments.

Algorithms for the computation of the matrix logarithm based on the double exponential formula

We consider the computation of the matrix logarithm by using numerical quadrature. The efficiency of numerical quadrature depends on the integrand and the choice of quadrature formula. The Gauss–Legendre quadrature has been conventionally employed; however, the convergence could be slow for ill-conditioned matrices. This effect may stem from the rapid change of the integrand values. To avoid such situations, we focus on the double exponential formula, which has been developed to address integrands with endpoint singularity. In order to utilize the double exponential formula, we must determine a suitable finite integration interval, which provides the required accuracy and efficiency. In this paper, we present a method for selecting a suitable finite interval based on an error analysis as well as two algorithms, and one of these algorithms is an adaptive quadrature algorithm.

Improvement of the double exponential formula with conformal maps based on the locations of singularities

Improvement of the double exponential formula with conformal maps based on the locations of singularities

Accurate and efficient computation of the 3D radiative transfer equation in highly forward-peaked scattering media using a renormalization approach

This paper considers an accurate and efficient numerical scheme for solving the radiative transfer equation (RTE) based on the discrete ordinates method (DOM) in highly forward-peaked scattering media. The DOM approximates the scattering integral of the RTE including the phase function into a quadrature sum with the total number of discrete angular directions (TND) by a quadrature set. Due to large numerical errors of the scattering integral based on the DOM in highly forward-peaked scattering, the phase function is renormalized to satisfy its normalization conditions. Although the renormalization approaches of the phase function improve the accuracy of the numerical results of the RTE, the computational efficiency of the RTE is still required. This paper develops the first order renormalization approach using the double exponential formula for three quadrature sets: level symmetric even, even and odd, and Lebedev sets in a wide range of the TND from 48 to 1454. Numerical errors of the three-dimensional time-dependent RTE are investigated by the analytical solutions of the RTE. The investigation shows that the level symmetric even set with the TND of 48 using the developed approach provides the most accurate results of the RTE among the quadrature sets in the range of the TND, while to obtain the same accuracy by the conventional zeroth order renormalization approach, the TND needs to be larger than 360. The results suggest the large reduction of computational loads by the developed approach to less than 10 percent from those in the conventional approach.

A new theory and automatic computation of reversible cyclic voltammograms at an inlaid disk electrode

By generalising the recent theory of chronoamperometry [L. K. Bieniasz, Electrochim. Acta 199 (2016) 1], a new, semi-analytical description of reversible cyclic voltammetry at an inlaid disk electrode is obtained, assuming equal diffusion coefficients of the electroactive species. In contrast to some former modelling studies, the new theory does not involve simplified or low-accurate approximations to the chronoamperometric current or integral transformation kernel, but it is mathematically rigorous. An explicit expression for the voltammetric current function is derived, in the form of an inverse Laplace transform of an infinite series involving spheroidal wave functions. An equivalent formula involving convolution integrals, more suitable for calculations, is also obtained, as well as an integral equation for the current function. The voltammograms are calculated automatically with a prescribed accuracy, by using either the adaptive integrator INTDE based on double exponential formulas, or the adaptive Huber method for integral equations. INTDE is more efficient, allowing one to compute a current function value with a relative error as close to ±10−16 as is possible, in a computing time of ca. 1 ms on a contemporary laptop computer.

An upper bound estimate and stability for the global error of numerical integration using double exponential transformation

The double exponential formula was introduced for calculating definite integrals with singular point oscillation functions and Fourier integral. The double exponential transformation is not only useful for numerical computations but it is also used in different methods of Sinc theory. In this paper we give an upper bound estimate for the error of double exponential transformation. By improving integral estimates having singular final points, in Theorem 1 we prove that the method is convergent and the rate of convergence is O(h2) where h is a step size. Our main tool in the proof is DE formula in Sinc theory. The advantage of our method is that the time and space complexity is drastically reduced. Furthermore, we discovered an upper bound error in DE formula independent of N-truncated number, as a matter of fact, we proved stability. Numerical tests are presented to verify the theoretical predictions and confirm the convergence of the numerical solution.