Equispaced Grid Research Articles

Spatially Dispersionless, Unconditionally Stable FC–AD Solvers for Variable-Coefficient PDEs

We present fast, spatially dispersionless and unconditionally stable high-order solvers for partial differential equations (PDEs) with variable coefficients in general smooth domains. Our solvers, which are based on (i) A certain “Fourier continuation” (FC) method for the resolution of the Gibbs phenomenon on equi-spaced Cartesian grids, together with (ii) A new, preconditioned, FC-based solver for two-point boundary value problems for variable-coefficient Ordinary Differential Equations, and (iii) An Alternating Direction strategy, generalize significantly a class of FC-based solvers introduced recently for constant-coefficient PDEs. The present algorithms, which are applicable, with high-order accuracy, to variable-coefficient elliptic, parabolic and hyperbolic PDEs in general domains with smooth boundaries, are unconditionally stable, do not suffer from spatial numerical dispersion, and they run at Fast Fourier Transform speeds. The accuracy, efficiency and overall capabilities of our methods are demonstrated by means of applications to challenging problems of diffusion and wave propagation in heterogeneous media.

Discrete direct methods in the fractional calculus of variations

Finite differences, as a subclass of direct methods in the calculus of variations, consist in discretizing the objective functional using appropriate approximations for derivatives that appear in the problem. This article generalizes the same idea for fractional variational problems. We consider a minimization problem with a Lagrangian that depends on the left Riemann–Liouville fractional derivative. Using the Grünwald–Letnikov definition, we approximate the objective functional in an equispaced grid as a multi-variable function of the values of the unknown function on mesh points. The problem is then transformed to an ordinary static optimization problem. The solution to the latter problem gives an approximation to the original fractional problem on mesh points.

Nonuniform Fourier Transforms for Rigid-Body and Multidimensional Rotational Correlations

The task of evaluating correlations is central to computational structural biology. The rigid-body correlation problem seeks the rigid-body transformation (R, t), R ∈ SO(3), t ∈ ℝ3 that maximizes the correlation between a pair of input scalar-valued functions representing molecular structures. Exhaustive solutions to the rigid-body correlation problem take advantage of the fast Fourier transform to achieve a speedup either with respect to the sought translation or rotation. We present PFcorr, a new exhaustive solution, based on the non-equispaced SO(3) Fourier transform, to the rigid-body correlation problem; unlike previous solutions, ours achieves a combination of translational and rotational speedups without requiring equispaced grids. PFcorr can be straightforwardly applied to a variety of problems in protein structure prediction and refinement that involve correlations under rigid-body motions of the protein. Additionally, we show how it applies, along with an appropriate flexibility model, to analogs of the above problems in which the flexibility of the protein is relevant.

Spectral Dynamics and Regularization of Incompletely and Irregularly Measured Data

Constrained by practical and economical aspects, in many applications, one often deals with data sampled irregularly and incompletely. The use of irregularly sampled data may result in some artifacts and poor spatial resolution. Therefore, the preprocessing of the measurements onto a regular grid plays an important step. One of the methods achieving this objective is based on the Fourier reconstruction, which involves an underdetermined system of equations. The recent Uniform Uncertainty Principle (UUP) uses convex optimization through l1 minimization for solving underdetermined systems. The l1 minimization admits certain theoretical guarantees and simpler implementation. The present work applies UUP to the Fourier-based data regularization problem. For the signals having sparse Fourier spectra, our method replaces the incomplete and irregular coordinate grid with the grid that is a subset of equispaced complete grid. It then generates error resulting from the stated replacement. Finally, it applies UUP to realize its objective. To justify the applicability of our method, we present the empirical performance of it on different sets of measurement coordinates as a function of number of nonzero Fourier coefficients.

O(dlog N)-Quantics Approximation of N-d Tensors in High-Dimensional Numerical Modeling

In the present paper, we discuss the novel concept of super-compressed tensor-structured data formats in high-dimensional applications. We describe the multifolding or quantics-based tensor approximation method of O(dlog N)-complexity (logarithmic scaling in the volume size), applied to the discrete functions over the product index set {1,…,N}⊗d, or briefly N-d tensors of size Nd, and to the respective discretized differential-integral operators in ℝd. As the basic approximation result, we prove that a complex exponential sampled on an equispaced grid has quantics rank 1. Moreover, a Chebyshev polynomial, sampled over a Chebyshev Gauss–Lobatto grid, has separation rank 2 in the quantics tensor format, while for the polynomial of degree m over a Chebyshev grid the respective quantics rank is at most 2m+1. For N-d tensors generated by certain analytic functions, we give a constructive proof of the O(dlog Nlog e−1)-complexity bound for their approximation by low-rank 2-(dlog N) quantics tensors up to the accuracy e>0. In the case e=O(N−α), α>0, our approach leads to the quantics tensor numerical method in dimension d, with the nearly optimal asymptotic complexity O(d/αlog 2e−1). From numerical examples presented here, we observe that the quantics tensor method has proved its value in application to various function related tensors/matrices arising in computational quantum chemistry and in the traditional finite element method/boundary element method (FEM/BEM). The tool apparently works.

Synthesis of Planar Arrays With Elements in Concentric Rings

The circular symmetry of concentric ring arrays becomes an advantage for situations where the radiation patterns are symmetric in the azimuthal direction. For planar arrays with elements in an equispaced grid, the fast Fourier transform can be applied in the calculations involved between the array factor and the source distribution. However, for circular arrays the elements are in a non-equispaced grid, making the application of fast computation algorithms difficult. Therefore, two techniques are proposed to synthesize circular arrays. The first one provides the control of the pattern by sampling the array factor. In this case, the fast Fourier transform can be applied. The windows technique is also utilized to control the sidelobe levels and the ripple structure of the shaped beam. The second technique considers the symmetric nature of concentric circular arrays and the control of the array factor is performed imposing some specific points of the pattern.

Spectral analysis and preconditioning techniques for radial basis function collocation matrices

SUMMARYMeshless collocation methods based on radial basis functions lead to structured linear systems, which, for equispaced grid points, have almost a multilevel Toeplitz structure. In particular, if we consider partial differential equations (PDEs) in two dimensions, then we find almost (up to a ‘low‐rank’ correction given by the boundary conditions) two‐level Toeplitz matrices, i.e. block Toeplitz with Toeplitz blocks structures, where both the number of blocks and the block‐size grow with the number of collocation points. In Bini et al. (Linear Algebra Appl. 2008; 428:508–519), upper bounds for the condition number of the Toeplitz matrices approximating a one‐dimensional model problem were proved. Here, we refine the one‐dimensional results, by explaining some numerics reported in the previous paper, and we show a preliminary analysis concerning conditioning, extremal spectral behavior, and global spectral results in the two‐dimensional case for the structured part. By exploiting the recent tools in the literature, a global distribution theorem in the sense of Weyl is proved also for the complete matrix‐sequence, where the low‐rank correction due to the boundary conditions is taken into consideration. The provided spectral analysis is then applied to design effective preconditioning techniques in order to overcome the ill‐conditioning of the matrices. A wide numerical experimentation, both in the one‐ and two‐dimensional cases, confirms our analysis and the robustness of the proposed preconditioners. Copyright © 2011 John Wiley & Sons, Ltd.

Effects of Discrete Charge Clustering in Simulations of Charged Interfaces.

A system of counterions between charged surfaces is investigated, with the surfaces represented by uniform charged planes and three different arrangements of discrete surface charges - an equispaced grid and two different clustered arrangements. The behaviors of a series of systems with identical net surface charge density are examined, with particular emphasis placed on the long ranged corrections via the method of "charged slabs" and the effects of the simulation cell size. Marked differences are observed in counterion distributions and the osmotic pressure dependent on the particular representation of the charged surfaces; the uniformly charged surfaces and equispaced grids of discrete charge behave in a broadly similar manner, but the clustered systems display a pronounced decrease in osmotic pressure as the simulation size is increased. The influence of the long ranged correction is shown to be minimal for all but the very smallest of system sizes.

Two-dimensional simulation of Rayleigh waves with staggered sine/cosine transforms and variable grid spacing

Simulation of Rayleigh waves requires high accuracy and an adequate spatial sampling at the surface. Discrete cosine and sine transforms are used to compute spatial derivatives along the direction perpendicular to the surface of the earth. Unlike the standard Fourier method, these transforms allow nonperiodic boundary conditions to be satisfied, in particular, the stress-free conditions at the surface. Because simulation of surface waves requires more points per minimum wavelength at the surface than simulation of body waves, the equispaced grid is not efficient. To overcome this problem, a grid compression is performed at the surface to obtain a denser spatial sampling. Grid size is minimal at the surface and increases with depth until reaching, at a relatively shallow depth, the grid points per wavelength required by the body waves. The stress-free boundary conditions are naturally handled by expanding the appropriate stress components in terms of the discrete sine transform. The wave equation is solved in the particle-velocity and stress formulation using a Runge-Kutta time integration and the convolutional PML (CPML) method to prevent reflections from the mesh boundaries. The simulations are very accurate for shallow sources and receivers and large offsets.

The MIPAS2D database of MIPAS/ENVISAT measurements retrieved with a multi-target 2-dimensional tomographic approach

Abstract. We present a multi-year database of atmospheric fields of the upper troposphere, stratosphere and lower mesosphere retrieved from satellite measurements adopting a 2-dimensional tomographic approach. The full mission of the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) instrument, on board the European Space Agency ENVISAT satellite, is analyzed with the Geofit Multi-Target Retrieval (GMTR) system to obtain the MIPAS2D database with atmospheric fields of pressure, temperature and volume mixing ratio of MIPAS main targets H2O, O3, HNO3, CH4, N2O, and NO2. The database covers both the MIPAS nominal observation mode measured at Full Resolution (FR) from July 2002 to March 2004 and the nominal observation mode of the new configuration, measured at Optimized Resolution (OR) and introduced in 2005. Further to the main targets, minor species N2O5, ClONO2, COF2, CFC-11, and CFC-12 for the FR mission only have been included in MIPAS2D to enhance its applicability in studies of stratospheric chemistry. The database is continuously updated with the analysis of the ongoing measurements that are planned to last until the end of 2013 and extended to other targets. The GMTR algorithm is operated on a fixed vertical grid coincident with the tangent altitudes of the FR nominal mode, spanning the altitude range from 6 to 68 km. In the horizontal domain, FR measurements are retrieved on both the observational grid and an equispaced 5 latitudinal-degrees grid which is made possible by the 2-dimensional retrieval algorithm. The analysis of MIPAS OR observations is operated on the same altitude-latitude fixed retrieval grid used for the FR measurements. This choice provides a database with a homogeneous altitude and latitude grid, over the whole globe, covering to date about seven years of measurements. The equispaced latitude grid provides a new and convenient layout for the much needed synergetic studies of data from various instrumental and modeling sources. MIPAS2D is available to the scientific community through the two web sites http://www.mbf.fci.unibo.it/mipas2d.html, and http://www.isac.cnr.it/~rss/mipas2d.htm

Stable Reconstruction of Piecewise Continuous Plane Stratified Biological Tissues via Electrical Impedance Tomography

Image reconstruction in electrical impedance tomography is, generally, an ill-posed nonlinear inverse problem. Regularization methods are widely used to ensure a stable solution. Herein, we present a case study, which uses a novel electrical impedance tomography method for reconstruction of layered biological tissues with piecewise continuous plane-stratified profiles. The algorithm implements the recently proposed reconstruction scheme for piecewise constant conductivity profiles, utilizing Legendre expansion in conjunction with improved Prony method. It is shown that the proposed algorithm is capable of successfully reconstructing piecewise continuous conductivity profiles with moderate slop. This reconstruction procedure, which calculates both the locations and the conductivities, repetitively provides inhomogeneous depth discretization, i.e., the depths grid is not equispaced. Incorporation of this specific inhomogeneous grid in the widely used mean least square reconstruction procedure results in a stable and accurate reconstruction, whereas, the commonly selected equispaced depth grid leads to unstable reconstruction. This observation establishes the main result of our investigation, highlighting the impact of physical phenomenon (the image series expansion) on electrical impedance tomography, leading to a physically motivated stabilization of the inverse problem, i.e., an inhomogeneous depth discretization renders an inherent regularization of the mean least square algorithm. The effectiveness and the significance of inhomogeneous discretization in electrical impedance tomography reconstruction procedure is further demonstrated and verified via numerical simulations.

Divergence (Runge Phenomenon) for least-squares polynomial approximation on an equispaced grid and Mock–Chebyshev subset interpolation

Runge showed more than a century ago that polynomial interpolation of a function f( x), using points evenly spaced on x ∈ [ - 1 , 1 ] , could diverge on parts of this interval even if f( x) was analytic everywhere on the interval. Least-squares fitting of a polynomial of degree N to an evenly spaced grid with P points should improve accuracy if P ≫ N . We show through numerical experiments that such an overdetermined fit reduces but does not eliminate the Runge Phenomenon. More precisely, define β ≡ ( N + 1 ) / P . The least-squares fit will fail to converge everywhere on [−1, 1] as N → ∞ for fixed β if f( x) has a singularity inside a curve D ( β ) in the complex-plane. The width of the region enclosed by the convergence boundary D shrinks as β diminishes (i. e., more points for a fixed polynomial degree), but shrinks to the interval [−1, 1] only when β → 0 . We also show that the Runge Phenomenon can be completely defeated by interpolation on a “mock–Chebyshev” grid: a subset of ( N + 1 ) points from an equispaced grid with O ( N 2 ) points chosen to mimic the non-uniform N + 1 -point Chebyshev–Lobatto grid.

On the efficient evaluation of Fourier patterns for nanoparticles and clusters

Samples made of an isotropically oriented ensemble of atomic clusters or structures that are not large crystals (i.e. extended less than 10 periods in each direction) are at the frontier of today's material science and chemistry. Examples are nanoparticles, nanotubes, amorphous matter, polymers, and macromolecules in suspension. For such systems the computation of powder diffraction patterns (which may provide an efficient characterization) is to be performed the hard way, by summing contributions from each atom pair. This work deals with performing such computation in the most practical and efficient way. Three main points are developed: how to encode the enormous array of interatomic distances (which increase as the square or higher powers of the cluster diameter) to a much smaller array of equispaced values on a coarse grid (whose size increases linearly with the diameter); how to perform a fast computation of the diffraction pattern from this equispaced grid; how to optimize the grid step to obtain an arbitrarily small error on the computed diffraction pattern. Theory and examples are jointly developed and presented.

Efficient Computation of Linearized Cross-Covariance and Auto-Covariance Matrices of Interdependent Quantities

In many geostatistical applications, spatially discretized unknowns are conditioned on observations that depend on the unknowns in a form that can be linearized. Conditioning takes several matrix–matrix multiplications to compute the cross-covariance matrix of the unknowns and the observations and the auto-covariance matrix of the observations. For large numbers n of discrete values of the unknown, the storage and computational costs for evaluating these matrices, proportional to n2, become strictly inhibiting. In this paper, we summarize and extend a collection of highly efficient spectral methods to compute these matrices, based on circulant embedding and the fast Fourier transform (FFT). These methods are applicable whenever the unknowns are a stationary random variable discretized on a regular equispaced grid, imposing an exploitable structure onto the auto-covariance matrix of the unknowns. Computational costs are reduced from O(n2) to O(nlog2n) and storage requirements are reduced from O(n2) to O(n).

Direct numerical simulation of turbulent spots

A group of high-order finite-difference schemes for incompressible flow was implemented to simulate the evolution of turbulent spots in channel flows. The long-time accuracy of these schemes was tested by comparing the evolution of small disturbances to a plane channel flow against the growth rate predicted by linear theory. When the perturbation is the unstable eigenfunction at a Reynolds number of 7500, the solution grows only if there are a comparatively large number of (equispaced) grid points across the channel. Fifth-order upwind biasing of convection terms is found to be worse than second-order central differencing. But, for a decaying mode at a Reynolds number of 1000, about a fourth of the points suffice to obtain the correct decay rate. We show that this is due to the comparatively high gradients in the unstable eigenfunction near the walls. So, high-wave-number dissipation of the high-order upwind biasing degrades the solution especially. But for a well-resolved calculation, the weak dissipation does not degrade solutions even over the very long times (O(100)) computed in these tests. Some new solutions of spot evolution in Couette flows with pressure gradients are presented. The approach to self-similarity at long times can be seen readily in contour plots.

Turbulence characteristics along the path of a heavy particle

The turbulence along the path of a particle in an isotropic turbulent flow is investigated using direct numerical simulation. The turbulence is simulated using a collocation spectral method on an equi-spaced three-dimensional grid of 128 3 points with a Taylor micro-scale Reynolds number of 65. The particulate phase is treated as a simple system where the flow is not affected by the presence of the particles. A total of 8192 particles are tracked through the flow using a Lagrangian advancement scheme with Stokes drag as the only force applied to the particles. The implementation of the forcing method used to maintain a constant level of kinetic energy within the flow is studied. The total fluid turbulence kinetic energy was found to be a function of particle Stokes number with smaller particles (St<1) seeing slightly higher average velocity fluctuations than a Lagrangian average of the turbulence and intermediate particles (1<St<10) seeing slightly lower average fluctuations. Likewise, the turbulence integral time scales of the fluid along the path of a particle appear slightly larger than both the Eulerian and Lagrangian integral time scales for these intermediate particles. All preferentially concentrated particles (0.1<St<10) saw a lower average vorticity and higher deformation of the fluid than a volume average of the turbulence. The preferential concentration of the particles peaked at a Stokes number of 1 with particles being on average 20% closer to the nearest neighboring particle than in a randomly distributed case. These data suggest that it is the vorticity measured along the path of a particle which is the dominant factor for variations of statistical properties in the particle Lagrangian reference frame.

Spatial Finite Difference Approximations for Wave-Type Equations

The simplest finite difference approximations for spatial derivatives are centered, explicit, and applied to "regular" equispaced grids. Well-established generalizations include the use of implicit (compact) approximations and staggered grids. We find here that the combination of these two concepts, together with high formal order of accuracy, is very effective for approximating the first derivatives in space that occur in many wave-type PDEs.

Seismic Hazard Estimation in Northern Algeria

In the present study, the seismic hazard in northern Algeria is estimated using both physical strain energy release and Gumbel's extreme values approaches. For six of the most industrial and populated cities in Algeria, seismic hazard is assessed and examined in greater detail. Gumbel's extreme values approach has been used to estimate seismic hazard in terms of magnitude and P.G.A at each point of an equispaced grid all over the north of Algeria. An average attenuation relationship for PGA has been provided using known relations which have been established in regions with similar attenuation characteristics. The results are presented mainly in the form of graphs and contour maps of magnitudes (respectively PGA) with a 60% probability of not being exceeded in the next 100 and 200 years. Globally, they give main features of northern Algeria in terms of zoning (as well as in terms of magnitude and in terms of PGA). They corroborate the ones obtained through other works, especially in the basin areas (Mitidja, Cheliff, Soumam and Constantine Basin).

Defeating the Runge phenomenon for equispaced polynomial interpolation via Tikhonov regularization

Runge showed that polynomial interpolation using an equispaced grid often diverges as the degree N of the interpolating polynomial f X increases, even when f(x) is analytic over the whole interval. We suppress Runge divergence by defining the approximating polynomial as the minimizer of the sum of the interpolation residual plus a constant - times a smoothness norm. The Tikhonov parameter - can be determined easily through the method of the L-shaped curve. The resulting “Tikhonov approximant” is at least as accurate as the truncated, N th degree Chebyshev series for the same function.

FIR digital filter design using non-equispaced frequency sampling

Frequency-sampling design of FIR digital filters is a well established technique. The method would be of even greater utility if the constraint frequencies could be specified to lie at arbitrary points, rather than just on the standard equispaced grid. Here we show how non-equispaced sampling can easily be achieved, and offer an example of a typical filter design using this flexible new approach.