Machine Roundoff Research Articles

Parameter Identification of the Discrete-Time Stochastic Systems with Multiplicative and Additive Noises Using the UD-Based State Sensitivity Evaluation

The paper proposes a new method for solving the parameter identification problem for a class of discrete-time linear stochastic systems with multiplicative and additive noises using a numerical gradient-based optimization. The constructed method is based on the application of a covariance UD filter for the above systems and an original method for evaluating state sensitivities within the numerically stable, matrix-orthogonal MWGS transformation. In addition to the numerical stability of the proposed algorithm to machine roundoff errors due to the application of the MWGS-UD orthogonalization procedure at each step, the main advantage of the obtained results is the possibility of analytical calculation of derivatives at a given value of the identified parameter without the need to use finite-difference methods. Numerical experiments demonstrate how the obtained results can be applied to solve the parameter identification problem for the considered stochastic system model.

A general approach for designing the MWGS-based information-form Kalman filtering methods

The paper addresses a general approach to MWGS (Modified Weighted Gram-Schmidt) orthogonalization based Kalman filtering (KF) implementation methods. We propose two new numerically favored and convenient array information formulations of the MWGS-based KF that are the MWGS-based array Information Filter (algorithm MWGS-aIF) and the extended MWGS-based array Information Filter (algorithm eMWGS-aIF). To confirm the correctness of our results, we have proved that the newly constructed MWGS-based array computational schemes are algebraically equivalent to the “straight” (conventional) information filter. Although all these information-type algorithms are theoretically equivalent, their computational properties are different. The newly proposed algorithms are numerically robust to machine roundoff errors due to the numerically stable orthogonal transformations applied on each iteration. The obtained numerical results confirm this statement. Additionally, algorithm eMWGS-aIF has the extended array form, i. e., it allows for updating all required filter quantities with the use of the numerically stable MWGS orthogonalization procedure, only. Thus, our results extend the existing class of numerically efficient KF implementation methods and can be used in practical applications.

Asymmetric Density Fitting with Modified Cholesky Decomposition Applied to Second-Order Electron Propagator.

Computation of molecular orbital electron repulsion integrals (MO-ERIs) as a transformation from atomic orbital ERIs (AO-ERIs) is the bottleneck of second-order electron propagator calculations when a single orbital is studied. In this contribution, asymmetric density fitting is combined with modified Cholesky decomposition to generate efficiently the required MO-ERIs. The key point of the presented algorithms is to keep track of integrals through partial contractions performed on three-center AO-ERIs; these contractions are stored in RAM instead of the AO-ERIs. Two implementations are provided, an in-core, which reduces the arithmetic and memory scaling factors as compared to the four-center AO-ERIs contraction method, and a semidirect, which overcomes memory limitations by evaluating antisymmetrized MO-ERIs in batches. On the numerical side, the proposed approach is fast and stable. The bad effects due to ill conditioning, namely, several negative and close to zero eigenvalues due to machine round off errors of the matrix associated with the density fitting process, are effectively controlled by means of a modified Cholesky factorization that avoids the matrix inversion needed to perform the asymmetrical density fitting implementation. The numerical experience presented shows that the in-core implementation is highly competitive to perform calculations on medium and large basis sets, while the semidirect implementation has small variations in time by changes in the available memory. The general applicability is illustrated on a set of selected relatively large-size molecules.

Compatible, total energy conserving and symmetry preserving arbitrary Lagrangian–Eulerian hydrodynamics in 2D rz – Cylindrical coordinates

We present a new discretization for 2D arbitrary Lagrangian–Eulerian hydrodynamics in rz geometry (cylindrical coordinates) that is compatible, total energy conserving and symmetry preserving.In the first part of the paper, we describe the discretization of the basic Lagrangian hydrodynamics equations in axisymmetric 2D rz geometry on general polygonal meshes. It exactly preserves planar, cylindrical and spherical symmetry of the flow on meshes aligned with the flow. In particular, spherical symmetry is preserved on polar equiangular meshes. The discretization conserves total energy exactly up to machine round-off on any mesh. It has a consistent definition of kinetic energy in the zone that is exact for a velocity field with constant magnitude.The method for discretization of the Lagrangian equations is based on ideas presented in [2,3,7], where the authors use a special procedure to distribute zonal mass to corners of the zone (subzonal masses). The momentum equation is discretized in its “Cartesian” form with a special definition of “planar” masses (area-weighted).The principal contributions of this part of the paper are as follows: a definition of “planar” subzonal mass for nodes on the z axis (r=0) that does not require a special procedure for movement of these nodes; proof of conservation of the total energy; formulated for general polygonal meshes.We present numerical examples that demonstrate the robustness of the new method for Lagrangian equations on a variety of grids and test problems including polygonal meshes. In particular, we demonstrate the importance of conservation of total energy for correctly modeling shock waves.In the second part of the paper we describe the remapping stage of the arbitrary Lagrangian–Eulerian algorithm. The general idea is based on the following papers [25–28], where it was described for Cartesian coordinates. We describe a distribution-based algorithm for the definition of remapped subzonal densities and a local constrained-optimization-based approach for each zone to find the subzonal mass fluxes.In this paper we give a systematic and complete description of the algorithm for the axisymmetric case and provide justification for our approach.1 The ALE algorithm conserves total energy on arbitrary meshes and preserves symmetry when remapping from one equiangular polar mesh to another.The principal contributions of this part of the paper are the extension of this algorithm to general polygonal meshes and 2D rz geometry with requirement of symmetry preservation on special meshes.We present numerical examples that demonstrate the robustness of the new ALE method on a variety of grids and test problems including polygonal meshes and some realistic experiments. We confirm the importance of conservation of total energy for correctly modeling shock waves.

Spherical Harmonics Power-spectrum of Global Geopotential Field of Gaussian-bell Type

Spherical harmonics power spectrum of the geopotential field of Gaussian-bell type on the sphere was investigated using integral formula that is associated with Legendre polynomials. The geopotential field of Gaussian-bell type is defined as a function of sine of angular distance from the bell's center in order to guarantee the continuity on the global domain. Since the integral-formula associated with the Legendre polynomials was represented with infinite series of polynomial, an estimation method was developed to make the procedure computationally efficient while preserving the accuracy. The spherical harmonics power spectrum was shown to vary significantly depending on the scale parameter of the Gaussian bell. Due to the accurate procedure of the new method, the power (degree variance) spanning over orders that were far higher than machine roundoff was well explored. When the scale parameter (or width) of the Gaussian bell is large, the spectrum drops sharply with the total wavenumber. On the other hand, in case of small scale parameter the spectrum tends to be flat, showing very slow decaying with the total wavenumber. The accuracy of the new method was compared with theoretical values for various scale parameters. The new method was found advantageous over discrete numerical methods, such as Gaussian quadrature and Fourier method, in that it can produce the power spectrum with accuracy and computational efficiency for all range of total wavenumber. The results of present study help to determine the allowable maximum scale parameter of the geopotential field when a Gaussian-bell type is adopted as a localized function.

A Cell-Integrated Semi-Lagrangian Semi-Implicit Shallow-Water Model (CSLAM-SW) with Conservative and Consistent Transport

Abstract A Cartesian semi-implicit solver using the Conservative Semi-Lagrangian Multitracer (CSLAM) transport scheme is constructed and tested for shallow-water (SW) flows. The SW equations solver (CSLAM-SW) uses a discrete semi-implicit continuity equation specifically designed to ensure a conservative and consistent transport of constituents by avoiding the use of a constant mean reference state. The algorithm is constructed to be similar to typical conservative semi-Lagrangian semi-implicit schemes, requiring at each time step a single linear Helmholtz equation solution and a single application of CSLAM. The accuracy and stability of the solver is tested using four test cases for a radially propagating gravity wave and two barotropically unstable jets. In a consistency test using the new solver, the specific concentration constancy is preserved up to machine roundoff, whereas a typical formulation can have errors many orders of magnitude larger. In addition to mass conservation and consistency, CSLAM-SW also ensures shape preservation by combining the new scheme with existing shape-preserving filters. With promising SW test results, CSLAM-SW shows potential for extension to a nonhydrostatic, fully compressible system solver for numerical weather prediction and climate models.

Fast and Accurate Con-Eigenvalue Algorithm for Optimal Rational Approximations

The need to compute small con-eigenvalues and the associated con-eigenvectors of positive-definite Cauchy matrices naturally arises when constructing rational approximations with a (near) optimally small $L^{\infty}$ error. Specifically, given a rational function with $n$ poles in the unit disk, a rational approximation with $m\ll n$ poles in the unit disk may be obtained from the $m$th con-eigenvector of an $n\times n$ Cauchy matrix, where the associated con-eigenvalue $\lambda_{m}>0$ gives the approximation error in the $L^{\infty}$ norm. Unfortunately, standard algorithms do not accurately compute small con-eigenvalues (and the associated con-eigenvectors) and, in particular, yield few or no correct digits for con-eigenvalues smaller than the machine roundoff. We develop a fast and accurate algorithm for computing con-eigenvalues and con-eigenvectors of positive-definite Cauchy matrices, yielding even the tiniest con-eigenvalues with high relative accuracy. The algorithm computes the $m$th con-eigenvalue in $\mathcal{O}\left(m^{2}n\right)$ operations and, since the con-eigenvalues of positive-definite Cauchy matrices decay exponentially fast, we obtain (near) optimal rational approximations in $\mathcal{O}(n\left(\log\delta^{-1}\right)^{2})$ operations, where $\delta$ is the approximation error in the $L^{\infty}$ norm. We provide error bounds demonstrating high relative accuracy of the computed con-eigenvalues and the high accuracy of the unit con-eigenvectors. We also provide examples of using the algorithm to compute (near) optimal rational approximations of functions with singularities and sharp transitions, where approximation errors close to machine roundoff are obtained. Finally, we present numerical tests on random (complex-valued) Cauchy matrices to show that the algorithm computes all the con-eigenvalues and con-eigenvectors with nearly full precision.

Communication: An exact short-time solver for the time-dependent Schrödinger equation.

The short-time integrator for propagating the time-dependent Schrödinger equation, which is exact to machine's round off accuracy when the Hamiltonian of the system is time-independent, was applied to solve dynamics processes. This integrator has the old Cayley's form [i.e., the Padé (1,1) approximation], but is implemented in a spectrally transformed Hamiltonian which was first introduced by Chen and Guo. Two examples are presented for illustration, including calculations of the collision energy-dependent probability passing over a barrier, and interaction process between pulse laser and the I(2) diatomic molecule.

A spurious evolution of turbulence originated from round-off error in pseudo-spectral simulation

We report a slowly-developing, spurious numerical solution in pseudo-spectral direct numerical simulation (DNS) of incompressible fluid turbulence. When the effect of machine round-off on the divergence-free condition is not carefully controlled, a problem can develop slowly (over about 50 large-eddy turnover times) and eventually leads to an unphysical flow field. The problem was found with a previously published, highly-compact algorithm for pseudo-spectral DNS and therefore it is important to document the contamination of this numerical artifact on simulated turbulence structure and statistics. This is a striking example since the problem is not easily noticeable due to its very long development time, and it does not lead to numerical instability but rather a different flow state. A theory is developed to explain the unphysical evolution and predicts the exponential growth of round-off error induced velocity divergence. The theory shows that any correlation of the large-scale forcing with the velocity field at the beginning of the time step could lead to amplification of the velocity divergence. For this reason, the problem is quite reproducible. Several simple remedies are tested and shown to correct the problem. It is shown that all revised algorithms are identical theoretically to the original algorithm, with the only difference in the level of control for the divergence-free condition of the simulated flow field. A general recommendation is that the pressure projection operation should be performed at the end of each time step to ensure that the divergence-free condition is not contaminated by machine round-off.

A study on the numerical stability of the two-fluid model near ill-posedness

The two-fluid model is widely used in studying gas–liquid flow inside pipelines because it can qualitatively predict the flow field at low computational cost. However, the two-fluid model becomes ill-posed when the slip velocity exceeds a critical value, and computations can be quite unstable before the flow reaches the ill-posed condition. In this work, computational stability of various convection schemes together with the Euler implicit method for the time derivatives in conjunction with the two-fluid model is analyzed. A pressure correction algorithm for the two-fluid model is carefully implemented to minimize its effect on numerical stability. von Neumann stability analysis shows that the central difference scheme is more accurate and more stable than the 1st-order upwind, 2nd-order upwind, and QUICK schemes. The 2nd-order upwind scheme is much more susceptible to instability than the 1st-order upwind scheme and is inaccurate for short waves. Excellent agreement is obtained between the predicted and computed growth rates of harmonic disturbances. The instability associated with the two-fluid model discretized system of equations is related to but quantitatively different from the instability associated with ill-posedness of the two-fluid model. When the computation becomes unstable due to the ill-posedness, the machine roundoff errors from a selected range of short wavelengths, which scale with the grid size, are amplified rapidly to render the computation of any targeted long wavelength variation useless. For the viscous two-fluid model with wall friction and interfacial drag, a small-amplitude long wavelength disturbance grows due to viscous Kelvin–Helmholtz instability without triggering the grid scale short waves when the system remains well posed. Under such a condition, central difference is found to be the most accurate discretization scheme among those investigated.

A fast immersed boundary method using a nullspace approach and multi-domain far-field boundary conditions

We report on the continued development of a projection approach for implementing the immersed boundary method for incompressible flows in two and three dimensions. Boundary forces and pressure are regarded as Lagrange multipliers that enable the no-slip and divergence-free constraints to be implicitly determined to arbitrary precision with no associated time-step restrictions. In order to accelerate the method, we further implement a nullspace (discrete streamfunction) method that allows the divergence-free constraint to be automatically satisfied to machine roundoff. By employing a fast sine transform technique, the linear system to determine the forces can be solved efficiently with direct or iterative techniques. A multi-domain technique is developed in order to improve far-field boundary conditions that are compatible with the fast sine transform and account for the extensive potential flow induced by the body as well as vorticity that advects/diffuses to large distance from the body. The multi-domain and fast techniques are validated by comparing to the exact solutions for the potential flow induced by stationary and propagating Oseen vortices and by an impulsively-started circular cylinder. Speed-ups of more than an order-of-magnitude are achieved with the new method.

Constrained Transport Algorithms for Numerical Relativity. I. Development of a Finite‐Difference Scheme

A scheme is presented for accurately propagating the gravitational field constraints in finite difference implementations of numerical relativity. The method is based on similar techniques used in astrophysical magnetohydrodynamics and engineering electromagnetics, and has properties of a finite differential calculus on a four-dimensional manifold. It is motivated by the arguments that 1) an evolutionary scheme that naturally satisfies the Bianchi identities will propagate the constraints, and 2) methods in which temporal and spatial derivatives commute will satisfy the Bianchi identities implicitly. The proposed algorithm exactly propagates the constraints in a local Riemann normal coordinate system; {\it i.e.}, all terms in the Bianchi identities (which all vary as $\partial^3 g$) cancel to machine roundoff accuracy at each time step. In a general coordinate basis, these terms, and those that vary as $\partial g\partial^2 g$, also can be made to cancel, but differences of connection terms, proportional to $(\partial g)^3$, will remain, resulting in a net truncation error. Detailed and complex numerical experiments with four-dimensional staggered grids will be needed to completely examine the stability and convergence properties of this method. If such techniques are successful for finite difference implementations of numerical relativity, other implementations, such as finite element (and eventually pseudo-spectral) techniques, might benefit from schemes that use four-dimensional grids and that have temporal and spatial derivatives that commute.

Improvement of some bounds on the stability of fast Helmholtz solvers

Three fast Helmholtz solvers based on the Fast Fourier Transform are studied with respect to their stability. It is proved that they are strongly stable, and that the backward relative errors can grow at most like O(log 2n)ρ0, where ρ0 is the machine roundoff unit.

Limiting semantics of numerical programs

This paper treats semantics of numerical programs generally, but is principally concerned with limits as error tends to zero. We define three kinds of semantics for a simple programming language whose only data type is a real number type: a family of finite accuracy semantics, which explicitly represent numerical error, and two types of limiting program semantics. One limiting semantics models the limit as machine roundoff error becomes small; the other models the limit over programs as the mathematical error inherent in the algorithm used tends to zero. These three kinds of program semantics are based on domains that are similar to those commonly used in denotational semantics, but are continuous rather than algebraic, a slightly weaker property. Using these continuous bounded-complete domains, we are able to show that the two kinds of limiting semantics are in fact limits of the finite accuracy semantics.

Stability of the block cyclic reduction

The forward stability of the block cyclic reduction without back substitution for block tridiagonal systems is studied. The basic assumption is that the matrix of the system is block column diagonally dominant. Then it is shown that for nonstrictly diagonally dominant matrices the forward error is O( C n n 2log 2 nκϱ 0), and for strictly diagonally dominant matrices it is O( C n g( s) log 2 nκϱ 0), where n is the block size of the matrix, N is the size of each block, g( s) is a function which measures the diagonal dominance, κ is a condition number, and ϱ 0 is the machine roundoff unit. At the end some numerical evidence is presented.

A rational approximation of the inverse normal probability function

We present coefficients for a rational approximation to the inverse normal probability function. The approximants are well suited for the double precision computation of the function. The error of the rational approximation is much smaller than the machine roundoff.

A simple approach to mode analysis for parabolic waveguides

A simple method of mode analysis for parabolic waveguides, based on one-dimensional analytic continuation, is presented. The method gives essentially exact values for parabolic cylinder functions; i.e., the relative error in the computed result is on the order of the machine round-off. When supplemented with a Newton-Poisson shooting method and simple homotopy techniques, this continuation method can be used to find the transverse electric and transverse magnetic mode eigenvalues, and associated separation constants, for arbitrary parabolic domains. These methods are also used to compute a power-handling efficiency factor for a range of parabolic regions.< <ETX xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">&gt;</ETX>

Flux-corrected transport techniques for multidimensional compressible magnetohydrodynamics

A prescription is given for conservatively integrating generalized hydromagnetic equations using flux-corrected transport (FCT) techniques. By placing the magnetic-field components at the interface locations of the finite-difference grid, the field is kept divergence-free to within machine roundoff error. The use of FCT techniques allows an integration scheme of high accuracy to be employed, while the numerical ripples associated with large dispersion errors are avoided. The method is particularly well suited for problems involving magnetohydrodynamic shocks and other discontinuities.

Modeling one‐dimensional infiltration into very dry soils: 1. Model development and evaluation

With the increasing economic growth in the arid regions of the West, and with the growing need for waste disposal and storage, the ability to efficiently model water flow and pollutant transport through unsaturated soils is becoming more important. One of the more difficult water flow problems to model, from a numerical point of view, is infiltration into very dry soils. The presence of very steep pressure gradients combined with the large field scales leads to algorithms that are very CPU intensive. Here we develop a water content based algorithm that is suitable for modeling one‐dimensional unsaturated water flow into layered soils. We show that this algorithm is a numerical approximation of a general form of Richards equation. We compare the computational efficiency of this algorithm with that of two pressure head based finite difference water flow algorithms for several test problems. We find that the mass balance errors for the noniterative water content formulation are of the order of machine round off error for most applications. We also find that for soils with fairly wet initial conditions (h ≈ −100 cm H2O) the water content formulation requires approximately the same CPU time as the faster of the two pressure head formulations. For very dry soils (h=−1000 to −50,000 cm H2O) the CPU time required for the water content formulation is not a function of the initial water content of the soil, whereas the CPU time required for the pressure head formulation strongly increases with decreasing initial water content. Because of this lack of sensitivity to initial conditions, the water content based algorithm is from 1 to 3 orders of magnitude faster than the pressure head based algorithms when applied to infiltration into very dry soils. The water content algorithm is not suitable for combined saturated‐unsaturated or near‐saturated flow that may be present because of local heterogeneities in the soil. In addition, the water content algorithm cannot handle positive pressure upper boundary conditions such as those associated with ponded surface water.

Perpendicularly Propagating Plasma Cyclotron Instabilities Simulated with a One-Dimensional Computer Model

Perpendicular propagating (k‖ = 0) cyclotron instabilities of velocity distributions f(υ⊥) = δ(υ⊥ − υ0)/2πυ⊥ and similar loss cone distributions with velocity spreads are simulated with a one-dimensional computer model. Essentially noise-free starting conditions are created by loading phase space in a uniform way. Growth starting from machine roundoff levels can be observed at precisely the linear Vlasov theory predictions, up to the saturated level, many orders of magnitude above the starting point. When the growth is dominated by a single wavenumber, the post-saturation time behavior of the electric field is dominated by recurring peaks separated by the cyclotron period. This causes the major frequency components of the nonlinear signal to be at or near harmonics of the cyclotron frequency rather than the linear real frequency, which is nonharmonic. The amplitude of succeeding bursts decays slowly, suggesting the existence of quasistable nonlinear states. A nonlinear instability has also been observed. Furthermore, we have examined another k‖ = 0 instability, a double hump mode arising from the combination of a stable loss cone plus a small amount of cold ions.