Accurate Numerical Integration Research Articles

Fast and Accurate Numerical Integration of the Langevin Equation with Multiplicative Gaussian White Noise.

A univariate stochastic system driven by multiplicative Gaussian white noise is considered. The standard method for simulating its Langevin equation of motion involves incrementing the system's state variable by a biased Gaussian random number at each time step. It is shown that the efficiency of such simulations can be significantly enhanced by incorporating the skewness of the distribution of the updated state variable. A new algorithm based on this principle is introduced, and its superior performance is demonstrated using a model of free diffusion of a Brownian particle with a friction coefficient that decreases exponentially with the kinetic energy. The proposed simulation technique proves to be accurate over time steps that are an order of magnitude longer than those required by standard algorithms. The model used to test the new numerical technique is known to exhibit a transition from normal diffusion to superdiffusion as the environmental temperature rises above a certain critical value. A simple empirical formula for the time-dependent diffusion coefficient, which covers both diffusion regimes, is introduced, and its accuracy is confirmed through comparison with the simulation results.

On an accurate numerical integration for the triangular and tetrahedral spectral finite elements

On an accurate numerical integration for the triangular and tetrahedral spectral finite elements

Tiresia: A code for molecular electronic continuum states and photoionization

The Tiresia program [1] provides access to numerically accurate solutions of the one-particle Schrödinger equation for highly excited states of complex polyatomic molecules, both bound and continuum, that cannot be described by conventional Quantum Chemistry approaches. It is based on an expansion of the required solution in a local multicentric basis set, with primitive functions built as products of a radial B-spline times a real spherical harmonic. In conjunction with Density Functional Theory (DFT), it has been extensively employed in a large variety of photoionization studies, also for rather large systems. Highly excited bound states as well as wavepacket propagation can also be accurately described. In fact, the flexibility of the basis essentially allows accurate solutions of linear operator equations, like Poisson or inhomogeneous perturbative equations, which are employed in the code. The program is parallelized with standard MPI-I instructions and makes extensive use of the Scalapack linear algebra library. Ancillary programs are available for the evaluation of photoionization cross sections and angular distributions from randomly to fully oriented molecules. Program summaryProgram Title: TiresiaCPC Library link to program files:https://doi.org/10.17632/fcrjxwgjxh.1Licensing provisions: GPLv3Programming language: Fortran77, Fortran90, MPISupplementary material: Program manual documentNature of problem: Accurate solutions for highly excited and continuum electronic states in complex polyatomic molecules. Molecular photoionization cross sections and angular distributions under high energy, high-intensity radiation pulses from Synchrotron radiation and laser sources, photoelectron imaging in pump-probe experiments, basis for electronic wavepackets under ultrafast or nonperturbative excitation.Solution method: Solution of the Schrödinger and similar linear operator equations in a finite domain is obtained via basis set expansion. Flexible basis set, that may approach practical completeness within the domain, is obtained as a multicenter set of B-spline radial functions times spherical harmonics. Accurate numerical integration is employed for the evaluation of matrix elements, and conventional diagonalization for bound states, or Galerkin approach for the full multichannel solution in the continuum. Full hamiltonian and dipole matrices in the spectral basis are available for time propagation. DFT many-body description is available, and strong correlations in the bound states may be incorporated via Dyson orbitals.Additional comments including restrictions and unusual features: The code is noted for computational efficiency, which allows fast yet reasonably accurate photoionization calculations for medium-sized molecules, allowing, e.g., calculations at many molecular geometries as required to describe time-resolved photoelectron spectra in pump-probe experiments.

Robust numerical integration of embedded solids described in boundary representation

Embedded and immersed methods have become essential tools in computational mechanics, as they allow discretizing arbitrarily complex geometries without the need for boundary-fitted meshes. One of their main challenges is the accurate numerical integration of cut elements. Among the various integration schemes developed for this purpose, moment fitting has proven to be a powerful technique that provides highly efficient and accurate integration rules.This publication presents a framework for the robust and efficient numerical integration of embedded solids described by oriented boundary meshes using moment fitting. The developments include an intersection algorithm that aims to drastically accelerate the computation of the necessary moments while achieving high accuracy. A closed surface parameterization of each cut domain is computed to facilitate the direct application of the divergence theorem. The algorithm is subject to a single quality criterion that guarantees the accurate evaluation of boundary integrals. At the same time, it allows to disregard classical mesh criteria, such as high aspect ratios, strongly varying angles, etc., resulting in extremely fast runtimes. In addition, an existing robust flood fill-based element classification scheme is further developed to initiate filling from arbitrary seed elements and to enable parallel execution, increasing its flexibility and efficiency.The successful application of all proposed algorithms to 4948 valid and flawed STLs from the Thingi10K database (Zhou and Jacobson, 2016) demonstrates their extraordinary robustness. In all cases, the wall-clock time scales at most linearly with the number of elements in the background mesh. We show that higher-order quadrature rules on the boundary elements enable efficient computation of the moments via the divergence theorem with near-machine precision. Finally, the presented methodologies are used to perform direct FE analyses on clean and flawed B-Rep models. All proposed algorithms are publicly available in the open-source C++ framework QuESo – Quadrature for Embedded Solids (https://github.com/manuelmessmer/QuESo), where the moment fitting equations are assembled and solved.

Higher-continuity s-version of finite element method with B-spline functions

This paper proposes a strategy to solve the problems of the conventional s-version of finite element method (SFEM) fundamentally. Because SFEM can reasonably model an analytical domain by superimposing meshes with different spatial resolutions, it has intrinsic advantages of local high accuracy, low computation time, and simple meshing procedure. However, it has disadvantages such as accuracy of numerical integration and matrix singularity. Although several additional techniques have been proposed to mitigate these limitations, they are computationally expensive or ad-hoc, and detract from the method's strengths. To solve these issues, we propose a novel strategy called B-spline based SFEM. To improve the accuracy of numerical integration, we employed cubic B-spline basis functions with C2-continuity across element boundaries as the global basis functions. To avoid matrix singularity, we applied different basis functions to different meshes. Specifically, we employed the Lagrange basis functions as local basis functions. The numerical results indicate that using the proposed method, numerical integration can be calculated with sufficient accuracy without any additional techniques used in conventional SFEM. Furthermore, the proposed method avoids matrix singularity and is superior to conventional methods in terms of convergence for solving linear equations. Therefore, the proposed method has the potential to reduce computation time while maintaining a comparable accuracy to conventional SFEM.

Revisiting Hansen’s Ideal Frame Propagation with Special Perturbations—1: Basic Algorithms for Osculating Elements

A review of the basic Hansen’s ideal frame algorithms for accurate numerical integration of perturbed elliptic motion is carried out. The fundamental approaches rely on the use of nonsingular variables and differ in the ways in which the ellipse in the orbital plane is determined. It is well known that the accuracy of the propagation of the orbit geometry is notably increased when using time-regularization techniques to transform the independent variable. However, this is at the expense of adding a differential equation to compute the time, which gathers the Lyapunov-type instabilities that are removed from the coordinates. The asynchronism resulting from errors in the numerical integration of the time may be palliated with the use of time elements, to which end a constant and a linear nonsingular time element are presented, which are new to our knowledge.

Indirect Optimization of Optimal Continuous Trajectory for Electric Sails With Refined Thrust Model in the Multi-Target Exploration Scenario

Electric sail-based propulsion is an innovative propellant-less propulsion technology that generates continuous thrust through the interaction between an artificial electric field and the solar wind. In an electric sail, the propulsive acceleration is adjusted by controlling the attitude of its normal plane and the coefficient determining the maximum thrust. However, the attitude adjustment speed of the electric sail is relatively small and the direction of the normal vector is constrained. Consequently, the electric sail is required to maintain a continuous propulsive acceleration vector when flying by the intermediate targets in multi-target interplanetary exploration. Therefore, an indirect optimization of three-dimensional optimal continuous interplanetary trajectory for electric sails with refined thrust model is investigated in this study. First, the optimal propulsive angles and thrust adjustment coefficient of electric sails with a refined thrust model are derived using Pontryagin's minimum principle. Second, a homotopy function is introduced in the process of trajectory optimization with an indirect method to approximate the step of the thrust adjustment coefficient to improve the accuracy of numerical integration. Additionally, the initial costates of the electric sail are transformed from the numerical simulation results of the Bezier shaping approach (BSA) and integrative Bezier shaping approach (IBSA) using the Karush-Kuhn-Tucker (KKT) condition. The numerical simulation results for the Earth-Mars rendezvous mission and the multi-target flyby mission reveal that the initial values of the costates are effective to implement an indirect optimization process of the optimal continuous trajectory for the rapid convergence of the electric sail. According to the numerical simulation results for multi-target mission, the indirect method is on average 2.5% better than the Gauss pseudospectral method (GPM) in overall index, and the calculation time is only 1.06% of GPM.

Quaternion and Biquaternion Methods and Regular Models of Analytical Mechanics (Review)

The work is of a survey analytical nature. The first part of the work presents quaternion and biquaternion methods for describing motion, models of the theory of finite displacements and regular kinematics of a rigid body based on the use of four-dimensional real and dual Euler (Rodrigues–Hamilton) parameters. These models, in contrast to the classical models of kinematics in Euler–Krylov angles and their dual counterparts, do not have division-by-zero features and do not contain trigonometric functions, which increases the efficiency of analytical research and numerical solution of problems in mechanics, inertial navigation, and motion control. The problem of regularization of differential equations of the perturbed spatial two-body problem, which underlies celestial mechanics and space flight mechanics (astrodynamics), is discussed using the Euler parameters, four-dimensional Kustaanheimo–Stiefel variables, and Hamilton quaternions: the problem of eliminating singularities (division by zero), which are generated by the Newtonian gravitational forces acting on a celestial or cosmic body and which complicate the analytical and numerical study of the motion of a body near gravitating bodies or its motion along highly elongated orbits. The history of the regularization problem and the regular Kustaanheim–Stiefel equations, which have found wide application in celestial mechanics and astrodynamics, are presented. We present the quaternion methods of regularization, which have a number of advantages over Kustaanheimo–Stiefel matrix regularization, and various regular quaternion equations of the perturbed spatial two-body problem (for both absolute and relative motion). The results of a comparative study of the accuracy of numerical integration of various forms of regularized equations of celestial mechanics and astrodynamics in Kustaanheimo–Stiefel variables and Newtonian equations in Cartesian coordinates are presented, showing that the accuracy of numerical integration of regularized equations in Kustaanheimo–Stiefel variables is much higher (by several orders of magnitude) than the accuracy of numerical integration Newtonian equations.

Adaptive Mesh Generation Technique for Efficient Electromagnetic Computation in RFIC Designs

A novel adaptive mesh generation technique for efficient electromagnetic simulation of radio-frequency integrated circuits (RFICs) is herein presented. By exploring the geometrical and physical characteristics of RFICs, some adaptive mesh treatments, such as mesh projection, edge refinement, via polymerization, etc., are utilized to improve the accuracy and efficiency of electromagnetic computations. For strong coupling structures, such as two conductors in close proximity for a relatively large area, a projection-based mesh scheme is introduced to improve the accuracy of numerical integration. Moreover, the current most likely concentrates near the edges of conductors due to the edge effect. To better model the edge effect, an edge refinement scheme is applied. For via arrays that appear common in RFICs, an automatic via aggregation approach is adopted to improve computational efficiency yet still keep good computational accuracy. Finally, some numerical examples are given to validate the computational accuracy and efficiency of the novel adaptive mesh generation technique.

Numerical integration rules with improved accuracy close to discontinuities

This work is devoted to the construction and analysis of a new nonlinear technique that allows obtaining accurate numerical integrations of any order using data that contains discontinuities, and when the integrand is only known at grid points. The novelty of the technique consists in the inclusion of correction terms with a closed expression that depend on the size of the jumps of the function and its derivatives at the discontinuities, that are supposed to be known. The addition of these terms allows recovering the accuracy of classical numerical integration formulas close to the discontinuities, as these correction terms account for the error that the classical integration formulas commit up to their accuracy at smooth zones. Thus, the correction terms can be added during the integration or as post-processing, which is useful if the main calculation of the integral has been already done using classical formulas. We include several numerical experiments that confirm the theoretical conclusions reached in this article.

Forecast and Correction of the Orbital Motion of the Space Vehicle Using Regular Quaternion Equations and Their Solutions in the Kustaanheimo–Stiefels Variables and Isochronic Derivatives

The regular quaternion equations of the orbital motion of a spacecraft (SC) proposed by us earlier in four-dimensional Kustaanheimo–Stiefel variables (KS-variables) are considered. These equations use as a new independent variable a variable related to real time by a differential relation (Sundman time transformation) containing the distance to the center of gravity. Various new regular quaternion equations in these variables and equations in regular quaternion osculating elements (slowly varying variables) are also constructed, in which the half generalized eccentric anomaly, widely used in celestial mechanics and space flight mechanics, is used as a new independent variable. Keplerian energy and time are used as additional variables in these equations. These equations are used to construct quaternion equations and relations in variations of KS-variables and their first derivatives and in variations of Keplerian energy and real time; the isochronous derivatives of the KS-variables and of their first derivatives and the matrix of isochronous derivatives for the elliptical Keplerian motion of the spacecraft are found, which are necessary for solving the problems of predicting and correcting its orbital motion. The results of a comparative study of the accuracy of the numerical integration of the Newtonian equations of the spatial restricted three-body problem (Earth, Moon, and spacecraft) in Cartesian coordinates and the regular quaternion equations of this problem in KS-variables are presented, which show that the accuracy of the numerical integration of regular quaternion equations is much higher (by several orders) of the accuracy of numerical integration of equations in Cartesian coordinates. This substantiates the expediency of using regular quaternion equations of the spacecraft orbital motion and the quaternion equations and relations in variations constructed in the article on their basis for the prediction and correction of the orbital motion of a spacecraft.

Leaky-Wave Analysis of TM-, TE-, and Hybrid-Polarized Aperture-Fed Bessel-Beam Launchers for Wireless-Power-Transfer Links

Bessel beams (BBs) can be generated at microwave frequencies by means of compact, planar devices based on Fabry–Perot cavity leaky-wave resonators. Most designs are based on the excitation of a single transverse magnetic/transverse electric (TM/TE) leaky mode by means of a vertical electric/magnetic dipole source. In this work, we analyze the important case of a horizontal magnetic dipole, which excites both TM and TE modes, developing an original theoretical framework suitable for resonant radiators. Analytical expressions are provided for the dominant leaky-wave aperture field. The near-field distribution is then evaluated through an accurate numerical integration of the radiating currents and validated through full-wave simulations for the relevant case of a BB launcher when excited by a waveguide-fed slot. Different designs are presented and analyzed in order to obtain TM, TE, and hybrid polarizations, finding in all cases an excellent agreement between simulations and theoretical results. Finally, we evaluate and compare the wireless-power-transfer efficiency of two coupled TM-, TE-, or hybrid-polarized BB launchers when a wireless link is established in the radiative near-field region. Interestingly, full-wave results demonstrate the superior performance of resonant BB launchers in TM or TE polarization with respect to the nonresonant hybrid case.

Accurate numerical integration of highly stiff pharmacokinetics models using continuous block implicit hybrid one-step collocation methods

The differential equations of pharmacokinetic models, obtained from the formulation based on the Fick's perfusion principle and law of conservation of mass action, deals with absorption, distribution, and elimination of drugs by the body systems. In the formulation, the ordinary differential equations obtained may result into highly stiff systems whereby suitable implicit methods have to be employed to solve accurately drug concentration in the body systems, owing to the complex nature of the exact solutions, if they exist. Further, because of the large and increasing interest in the problems of drug kinetics, it is necessary to apply considerably accurate implicit numerical methods in evaluating and in applying them to specific cases. These solutions are to help determine the distribution of drug concentration in different compartments (parts) of the human body network systems as a fundamental step to help in understanding and improving treatments. This is necessary since drug concentration in each compartment is different from one another because of the differences in drug affinity to tissues. From the study, the solution curves obtained show that drug level in the gastrointestinal tract decreases with the passage of time, while drug concentration in the blood increases from zero and reaches its maximum level and then decreases steadily again. We further, compared the model curves with some experimental data plots published in scientific papers. The results obtained from the study can be used extensively for various drug diffusion problems arising in pharmaceutical studies. The presented results widen the applicability of the continuous block implicit hybrid one-step collocation methods to diffusion process which have good number of applications in the drug control, drug dosage and other related problems in pharmaceutical and biomedical industries.

A meshless fragile points method for the solution of the monodomain model for cardiac electrophysiology simulation

Meshless methods for in silico modeling and simulation of cardiac electrophysiology are gaining more and more popularity. These methods do not require a mesh and are more suitable than the Finite Element Method (FEM) to simulate the activity of complex geometrical structures like the human heart. However, challenges such as numerical integration accuracy and time efficiency remain, which limit their applicability. Recently, the Fragile Points Method (FPM) has been introduced in the meshless methods family. It uses local, simple, polynomial, discontinuous functions to construct trial and test functions in the Galerkin weak form. This allows for accurate integration and improved efficiency while enabling the imposition of essential and natural boundary conditions as in FEM. In this work, we consider the application of FPM for cardiac electrophysiology simulation. We derive the cardiac monodomain model using the FPM formulation and we solve several benchmark problems in 2D and 3D. We show that FPM leads to solutions of accuracy and efficiency similar to FEM while alleviating the need for a mesh. Additionally, FPM demonstrates comparable convergence to FEM in the considered benchmarks.

ARPIST: Provably accurate and stable numerical integration over spherical triangles

Numerical integration on spherical triangles, including the computation of their areas, is a core computation in geomathematics. The commonly used techniques sometimes suffer from instabilities and significant loss of accuracy. We describe a new algorithm, called ARPIST, for accurate and stable integration of functions on spherical triangles. ARPIST is based on an easy-to-implement transformation to the spherical triangle from its corresponding linear triangle via radial projection to achieve high accuracy and efficiency. More importantly, ARPIST overcomes potential instabilities in computing the Jacobian of the transformation, even for poorly shaped triangles that may occur at poles in regular longitude-latitude meshes, by avoiding potential catastrophic rounding errors. We compare our proposed technique with L’Huilier’s Theorem for computing the area of spherical triangles, and also compare it with the recently developed LSQST method (Beckmann et al., 2014) and a radial-basis-function-based technique (Reeger and Fornberg, 2016) for integration of smooth functions on spherical triangulations. Our results show that ARPIST enables better or comparable accuracy over previous methods while being easier to implement, significantly faster, and more tolerant of poor element shapes.

Hierarchical deep learning of multiscale differential equation time-steppers.

Nonlinear differential equations rarely admit closed-form solutions, thus requiring numerical time-stepping algorithms to approximate solutions. Further, many systems characterized by multiscale physics exhibit dynamics over a vast range of timescales, making numerical integration expensive. In this work, we develop a hierarchy of deep neural network time-steppers to approximate the dynamical system flow map over a range of time-scales. The model is purely data-driven, enabling accurate and efficient numerical integration and forecasting. Similar ideas can be used to couple neural network-based models with classical numerical time-steppers. Our hierarchical time-stepping scheme provides advantages over current time-stepping algorithms, including (i) capturing a range of timescales, (ii) improved accuracy in comparison with leading neural network architectures, (iii) efficiency in long-time forecasting due to explicit training of slow time-scale dynamics, and (iv) a flexible framework that is parallelizable and may be integrated with standard numerical time-stepping algorithms. The method is demonstrated on numerous nonlinear dynamical systems, including the Van der Pol oscillator, the Lorenz system, the Kuramoto–Sivashinsky equation, and fluid flow pass a cylinder; audio and video signals are also explored. On the sequence generation examples, we benchmark our algorithm against state-of-the-art methods, such as LSTM, reservoir computing and clockwork RNN.This article is part of the theme issue ‘Data-driven prediction in dynamical systems’.

An XFEM-based numerical strategy to model three-dimensional fracture propagation regarding crack front segmentation

Under mixed I-III loading, a crack frequently segments into several echelon-shaped daughter cracks rotated towards the direction of principal stress. This phenomenon is called crack front segmentation and has been widely observed at different length scales in both laboratory and nature. In this paper, based on the extended finite element method (XFEM), a three-dimensional numerical strategy is proposed, validated, and adopted to model quasi-static fracture propagation considering crack front segmentation. A new crack propagation criterion based on the weighted average principal stress in conjunction with a cubic spline curve fitting of front vertexes is introduced to update the crack geometry. Further, to get a higher resolution of the stress field around the crack front, a robust and efficient local mesh refinement scheme fully integrated with the XFEM is proposed. Besides, a partitioning scheme of enriched elements containing segmentation points is presented to achieve an accurate numerical integration. Several examples are presented to illustrate the effectiveness and robustness of the proposed strategy. The hydraulic fracturing example indicates that ignoring fracture front segmentation for a mixed-mode I-III fluid-driven fracture will severely underestimate the injection pressure due to the stress shadow effects between adjacent overlapped segments.

Artificial neural network assisted numerical quadrature in finite element analysis in mechanics

Finite element analysis (FEA) is a numerical tool that is broadly implemented in solving engineering problems. The Gaussian quadrature rule is commonly used to integrate elemental quantities like stiffness matrices over arbitrarily shaped elements in FEA. The Gaussian quadrature rule is a numerical technique in which the accuracy of numerical integration depends on the number of Gauss quadrature points. A higher number of quadrature points in the Gaussian quadrature rule implies higher accuracy but at the cost of higher computations. Thus, the number of Gauss points directly affects the computational cost and accuracy of the FE results. The main objective of the paper is to develop a deep learning model to predict the optimum number of Gauss points required in three-dimensional FEA utilizing hexahedral elements. In recent times some works have been proposed using deep learning in this direction. The present work further explores new possibilities in the area. In this paper, instead of the coordinates of the hexahedral element, its length and angle between two sides are taken as the input to the model, and the optimum number of Gauss quadrature points are predicted. The performance of the model is evaluated by the accuracy and confusion matrix. The present results are compared and are found to be better than the results available in the literature.

A numerical integration strategy of meshless numerical manifold method based on physical cover and applications to linear elastic fractures

The meshless numerical manifold method (MNMM) inherits two covers of the numerical manifold method. A mathematical cover is composed of nodes' influence domains and a physical cover consists of physical patches, which are produced through cutting mathematical cover by physical boundaries. Because two covers are adopted, MNMM can naturally and uniquely solve both the continuous and discontinuous problems under Galerkin's variational framework. However, Galerkin's meshless method needs background integration grids to realize solving, which often does not match the nodes' influence domains, so the accuracy of numerical integration is reduced. Consider that MNMM allows even distribution of nodes and the physical cover contains the characteristics of the boundary and nodes' influence domains, the study presents a new numerical integration strategy to ensure that the background integration grids match the nodes' influence domains. The method can be applied to continuous and discontinuous problems, and is proved to be equivalent to the influence domain integration. At the same time, a reasonable arrangement of mathematical nodes is made to assure the background integration grid that it is accordant and straightforward. In this way, the number of physical patches of each integral point is the same, which improves the accuracy of interpolation calculation. The effectiveness of the proposed method is verified by numerical examples of both continuous and discontinuous problems.

Forced Variational Integrators for the Formation Control of Multiagent Systems

Formation control of autonomous agents can be seen as a physical system of individuals interacting with local potentials, and whose evolution can be described by a Lagrangian function. In this article, we construct and implement forced variational integrators for the formation control of autonomous agents modeled by double integrators. In particular, we provide an accurate numerical integrator with a lower computational cost than traditional solutions. We find error estimations for the rate of the energy dissipated along with the agents' motion to achieve desired formations. Consequently, this permits to providing sufficient conditions on the time step for the convergence of discrete formation control systems such as the consensus problem in discrete systems. We present practical applications such as the rapid estimation of regions of attraction to desired shapes in distance-based formation control.