Showing papers in "arXiv: Numerical Analysis in 2001"

Adaptive Mollifiers for High Resolution Recovery of Piecewise Smooth Data from its Spectral Information

Abstract: We discuss the reconstruction of piecewise smooth data from its (pseudo-) spectral information. Spectral projections enjoy superior resolution provided the data is globally smooth, while the presence of jump discontinuities is responsible for spurious ${\cal O}(1)$ Gibbs oscillations in the neighborhood of edges and an overall deterioration to the unacceptable first-order convergence rate. The purpose is to regain the superior accuracy in the piecewise smooth case, and this is achieved by mollification. Here we utilize a modified version of the two-parameter family of spectral mollifiers introduced by Gottlieb & Tadmor [GoTa85]. The ubiquitous one-parameter, finite-order mollifiers are based on dilation. In contrast, our mollifiers achieve their high resolution by an intricate process of high-order cancelation. To this end, we first implement a localization step using edge detection procedure, [GeTa00a, GeTa00b]. The accurate recovery of piecewise smooth data is then carried out in the direction of smoothness away from the edges, and adaptivity is responsible for the high resolution. The resulting adaptive mollifier greatly accelerates the convergence rate, recovering piecewise analytic data within exponential accuracy while removing spurious oscillations that remained in [GoTa85]. Thus, these adaptive mollifiers offer a robust, general-purpose ``black box'' procedure for accurate post processing of piecewise smooth data.

High resolution conjugate filters for the simulation of flows

Abstract: This paper proposes a Hermite-kernel realization of the conjugate filter oscillation reduction (CFOR) scheme for the simulation of fluid flows. The Hermite kernel is constructed by using the discrete singular convolution (DSC) algorithm, which provides a systematic generation of low-pass filter and its conjugate high-pass filters. The high-pass filters are utilized for approximating spatial derivatives in solving flow equations, while the conjugate low-pass filter is activated to eliminate spurious oscillations accumulated during the time evolution of a flow. As both low-pass and high-pass filters are derived from the Hermite kernel, they have similar regularity, time-frequency localization, effective frequency band and compact support. Fourier analysis indicates that the CFOR-Hermite scheme yields a nearly optimal resolution and has a better approximation to the ideal low-pass filter than previously CFOR schemes. Thus, it has better potential for resolving natural high frequency oscillations from a shock. Extensive one- and two-dimensional numerical examples, including both incompressible and compressible flows, with or without shocks, are employed to explore the utility, test the resolution, and examine the stability of the present CFOR-Hermite scheme. Extremely small ratio of point-per-wavelength (PPW) is achieved in solving the Taylor problem, advancing a wavepacket and resolving a shock/entropy wave interaction. The present results for the advection of an isentropic vortex compare very favorably to those in the literature.

Partial Riemann problem, boundary conditions, and gas dynamics

Abstract: We introduce in this contribution the notion of partial Riemann problem. Recall that the Riemann problem describes a shock tube interaction between two given states ; the partial Riemann problem is a generalization of the previous concept and introduces the notion of boundary manifold. In what follows, we first recall very classical notions concerning gas dynamics and the associated Riemann problem. In a second part, we introduce the partial Riemann problem for general systems of conservation laws and proves that this problem admits a solution in some class of appropriate nonlinear waves. In section 3, we recall the linearized analysis with the method of characteristics, introduce the weak formulation of the Dirichlet boundary condition for nonlinear situations in terms of the partial Riemann problem and show that lot of physically relevant situations are described with this theoretical framework. In the last paragraph, we propose a practical implementation of the previous onsiderations with the finite volume method.

A unifying approach to software and hardware design for scientific calculations and idempotent mathematics

Abstract: A unifying approach to software and hardware design generated by ideas of Idempotent Mathematics is discussed. The so-called idempotent correspondence principle for algorithms, programs and hardware units is described. A software project based on this approach is presented.

Bayesian Estimation of Time Series Lags and Structure

Abstract: This paper derives practical algorithms, based on Bayesian inference methods, for several data analysis problems common in time series analysis of astronomical and other data. One problem is the determination of the lag between two time series, for which the cross-correlation function is a sufficient statistic. The second problem is the estimation of structure in a time series of measurements which are a weighted integral over a finite range of the independent variable.

Shock-capturing with natural high frequency oscillations

Abstract: This paper explores the potential of a newly developed conjugate filter oscillation reduction (CFOR) scheme for shock-capturing under the influence of natural high-frequency oscillations. The conjugate low-pass and high-pass filters are constructed based on the principle of the discrete singular convolution. Two Euler systems, the advection of an isentropy vortex flow and the interaction of shock-entropy wave are considered to demonstrate the utility of the CFOR scheme. Computational accuracy and order of approximation are examined and compared with the literature. Some of the best numerical results are obtained for the shock-entropy wave interaction. Numerical experiments indicate that the proposed scheme is stable, conservative and reliable for the numerical simulation of hyperbolic conservation laws.

Bayesian Blocks in Two or More Dimensions: Image Segmentation and Cluster Analysis

Abstract: This paper describes an extension, to higher dimensions, of the Bayesian Blocks algorithm for estimating signals in noisy time series data (Scargle 1998, 2000). The mathematical problem is to find the partition of the data space with the maximum posterior probability for a model consisting of a homogeneous Poisson process for each partition element. For model M_{n}, attributing the data within region n of the data space to a Poisson process with a fixed event rate lambda_{n}, the global posterior is: P(M_{n}) = Phi(N,V) = Gamma(N+1)Gamma(V-N+1) / Gamma(V+2) = N!(V-N)!/(V+1)! . Note that lambda_{n} does not appear, since it has been marginalized, using a flat, improper prior. Other priors yield similar formulas. This expression is valid for a data space of any dimension. It depends on only N, the number of data points within the region, and V, the volume of the region. No information about the actual locations of the points enters this expression. Suppose two such regions, described by N_{1},V_{1} and N_{2},V_{2}, are candidates for being merged into one. From the above equation, construct a Bayes merge factor, giving the ratio of posteriors for the two regions merged and not merged, respectively: P(Merge) = Phi(N_{1}+N_{2},V_{1}+V_{2}) / Phi(N_{1},V_{1}) Phi(N_{2},V_{2}) . Then collect data points into blocks with a greedy cell coalescence algorithm.

New RBF collocation schemes and their applications

Abstract: The purpose of this study is to apply some new RBF collocation schemes and recently-developed kernel RBFs to various types of partial differential equation systems. By analogy with the Fasshauer's Hermite interpolation, we recently developed the symmetric BKM and boundary particle methods (BPM), where the latter is based on the multiple reciprocity principle. The resulting interpolation matrix of them is always symmetric irrespective of boundary geometry and conditions. Furthermore, the proposed direct BKM and BPM apply the practical physical variables rather than expansion coefficients and become very competitive alternative to the boundary element method. On the other hand, by using the Green integral we derive a new domain-type symmetrical RBF scheme called as the modified Kansa method (MKM), which differs from the Fasshaure's scheme in that the MKM discretizes both governing equation and boundary conditions on the same boundary nodes. Therefore, the MKM significantly reduces calculation errors at nodes adjacent to boundary with explicit mathematical basis. Experimenting these novel RBF schemes with 2D and 3D Laplace, Helmholtz, and convection-diffusion problems will be subject of this study. In addition, the nonsingular high-order fundamental or general solution will be employed as the kernel RBFs in the BKM and MKM.

A Priori Estimates for the Global Error Committed by Runge-Kutta Methods for a Nonlinear Oscillator

Abstract: The Alekseev-Gr{\"o}bner lemma is combined with the theory of modified equations to obtain an \emph{a priori} estimate for the global error of numerical integrators. This estimate is correct up to a remainder term of order $h^{2p}$, where $h$ denotes the step size and $p$ the order of the method. It is applied to a class of nonautonomous linear oscillatory equations, which includes the Airy equation, thereby improving prior work which only gave the $h^p$ term. Next, nonlinear oscillators whose behaviour is described by the Emden-Fowler equation $y'' + t^ u y^n = 0$ are considered, and global errors committed by Runge-Kutta methods are calculated. Numerical experiments show that the resulting estimates are generally accurate. The main conclusion is that we need to do a full calculation to obtain good estimates: the behaviour is different from the linear case, it is not sufficient to look only at the leading term, and merely considering the local error does not provide an accurate picture either.

A conservative semi-Lagrangian method for oscillation-free computation of advection processes

Abstract: The semi-Lagrangian method using the hybrid-cubic-rational interpolation function [M. Ida, Comput. Fluid Dyn. J. 10 (2001) 159] is modified to a conservative method by incorporating the concept discussed in [R. Tanaka et al., Comput. Phys. Commun. 126 (2000) 232]. In the method due to Tanaka et al., not only a physical quantity but also its integrated quantity within a computational cell are used as dependent variables, and the mass conservation is completely achieved by giving a constraint to a forth-order polynomial used as an interpolation function. In the present method, a hybrid-cubic-rational function whose optimal mixing ratio was determined theoretically is employed for the interpolation, and its derivative is used for updating the physical quantity. The numerical oscillation appearing in results by the method due to Tanaka et al. is sufficiently eliminated by the use of the hybrid function.

Idempotent interval analysis and optimization problems

Abstract: Many problems in optimization theory are strongly nonlinear in the traditional sense but possess a hidden linear structure over suitable idempotent semirings. After an overview of `Idempotent Mathematics' with an emphasis on matrix theory, interval analysis over idempotent semirings is developed. The theory is applied to construction of exact interval solutions to the interval discrete stationary Bellman equation. Solution of an interval system is typically NP-hard in the traditional interval linear algebra; in the idempotent case it is polynomial. A generalization to the case of positive semirings is outlined.

A Generalization of Obreshkoff-Ehrlich Method for Multiple Roots of Algebraic, Trigonometric and Exponential Equations

Abstract: In this paper methods for simultaneous finding all roots of generalized polynomials are developed. These methods are related to the case when the roots are multiple. They possess cubic rate of convergence and they are as labour-consuming as the known methods related to the case of polynomials with simple roots only.

Detection of Edges in Spectral Data II. Nonlinear Enhancement

Abstract: We discuss a general framework for recovering edges in piecewise smooth functions with finitely many jump discontinuities, where $[f](x):=f(x+)-f(x-) eq 0$. Our approach is based on two main aspects--localization using appropriate concentration kernels and separation of scales by nonlinear enhancement. To detect such edges, one employs concentration kernels, $K_\epsilon(\cdot)$, depending on the small scale $\epsilon$. It is shown that odd kernels, properly scaled, and admissible (in the sense of having small $W^{-1,\infty}$-moments of order ${\cal O}(\epsilon)$) satisfy $K_\epsilon*f(x) = [f](x) +{\cal O}(\epsilon)$, thus recovering both the location and amplitudes of all this http URL an example we consider general concentration kernels of the form $K^\sigma_N(t)=\sum\sigma(k/N)\sin kt$ to detect edges from the first $1/\epsilon=N$ spectral modes of piecewise smooth f's. Here we improve in generality and simplicity over our previous study in [A. Gelb and E. Tadmor, Appl. Comput. Harmon. Anal., 7 (1999), pp. 101-135]. Both periodic and nonperiodic spectral projections are considered. We identify, in particular, a new family of exponential factors, $\sigma^{exp}(\cdot)$, with superior localization properties. The other aspect of our edge detection involves a nonlinear enhancement procedure which is based on separation of scales between the edges, where $K_\epsilon*f(x)\sim [f](x) eq 0$, and the smooth regions where $K_\epsilon*f = {\cal O}(\epsilon) \sim 0$. Numerical examples demonstrate that by coupling concentration kernels with nonlinear enhancement one arrives at effective edge detectors.

On solving fewnomials over intervals in fewnomial time

Abstract: Let f be a degree D univariate polynomial with real coefficients and exactlym monomial terms. We show that in the special case m=3 we can approximate within e all the roots of f in the interval (0, R) using just O log(D)log Dlog R e �� arithmetic operations. In particular, we can count the number of roots in any bounded interval using just O(log 2 D) arithmetic operations. Our speed-ups are significant and near-optimal: The asymptotically sharpest previous complexity upper bounds for both problems were super-linear in D, while our algorithm has complexity close to the respective complexity lower bounds. We also dis- cuss conditions under which our algorithms can be extended to general m, and a connection to a real analogue of Smale's 17 th Problem.

ABS Algorithms for Linear Systems and Optimization

Abstract: We present a review and bibliography of the main results obtained during a research on ABS (Abaffy, Broyden, Spedicato) methods

Interval straight line fitting

Abstract: I consider the task of experimental data fitting. Unlike the traditional approach I do not try to minimize any functional based on available experimental information, instead the minimization problem is replaced with constraint satisfaction procedure, which produces the interval hull of solutions of desired type. The method, called 'box slicing algorithm', is described in details. The results obtained this way need not to be labeled with confidence level of any kind, they are simply certain (guaranteed). The method easily handles the case with uncertainties in one or both variables. There is no need for, always more or less arbitrary, weighting the experimental data. The approach is directly applicable to other experimental data processing problems like outliers detection or finding the straight line, which is tangent to the experimental curve.

Holistically discretise the Swift-Hohenberg equation on a scale larger than its spatial pattern

Abstract: I introduce an innovative methodology for deriving numerical models of systems of partial differential equations which exhibit the evolution of spatial patterns. The new approach directly produces a discretisation for the evolution of the pattern amplitude, has the rigorous support of centre manifold theory at finite grid size $h$, and naturally incorporates physical boundaries. The results presented here for the Swift-Hohenberg equation suggest the approach will form a powerful method in computationally exploring pattern selection in general. With the aid of computer algebra, the techniques may be applied to a wide variety of equations to derive numerical models that accurately and stably capture the dynamics including the influence of possibly forced boundaries.

Direct linearization method for nonlinear PDE's and the related kernel RBFs

Abstract: The standard methodology handling nonlinear PDE's involves the two steps: numerical discretization to get a set of nonlinear algebraic equations, and then the application of the Newton iterative linearization or its variants to solve the nonlinear algebraic systems. Here we present an alternative strategy called direct linearization method (DLM). The DLM discretization algebraic equations of nonlinear PDE's is simply linear rather than nonlinear. The basic idea behind the DLM is that we see a nonlinear term as a new independent systematic variable and transfer a nonlinear PDE into a linear PDE with more than one independent variable. It is stressed that the DLM strategy can be applied combining any existing numerical discretization techniques. The resulting linear discretization equations can be either over-posed or well-posed. In particular, we also discuss how to create proper radial basis functions in conjunction with the DLM.

On Non-Oscillating Integrals for Computing Inhomogeneous Airy Functions

Abstract: Integral representations are considered of solutions of the inhomogeneous Airy differential equation $w''-z w=\pm1/\pi$. The solutions of these equations are also known as Scorer functions. Certain functional relations for these functions are used to confine the discussion to one function and to a certain sector in the complex plane. By using steepest descent methods from asymptotics, the standard integral representations of the Scorer functions are modified in order to obtain non-oscillating integrals for complex values of $z$. In this way stable representations for numerical evaluations of the functions are obtained. The methods are illustrated with numerical results.

Numerical Performance of Abs Codes for Nonlinear Systems of Equations

Abstract: Numerical experiments are performed in order to study the performance of ABS codes in solving non-linear systems of equations

A fast solver for systems of reaction-diffusion equations

Abstract: In this paper we present a fast algorithm for the numerical solution of systems of reaction-diffusion equations, $\partial_t u + a \cdot abla u = \Delta u + F (x, t, u)$, $x \in \Omega \subset \mathbf{R}^3$, $t > 0$. Here, $u$ is a vector-valued function, $u \equiv u(x, t) \in \mathbf{R}^m$, $m$ is large, and the corresponding system of ODEs, $\partial_t u = F(x, t, u)$, is stiff. Typical examples arise in air pollution studies, where $a$ is the given wind field and the nonlinear function $F$ models the atmospheric chemistry.

Derive boundary conditions for holistic discretisations of Burgers' equation

Abstract: I previously used Burgers' equation to introduce a new method of numerical discretisation of \pde{}s. The analysis is based upon centre manifold theory so we are assured that the discretisation accurately models all the processes and their subgrid scale interactions. Here I show how boundaries to the physical domain may be naturally incorporated into the numerical modelling of Burgers' equation. We investigate Neumann and Dirichlet boundary conditions. As well as modelling the nonlinear advection, the method naturally derives symmetric matrices with constant bandwidth to correspond to the self-adjoint diffusion operator. The techniques developed here may be used to accurately model the nonlinear evolution of quite general spatio-temporal dynamical systems on bounded domains.

Phase retrieval by iterated projections

Abstract: Several strategies in phase retrieval are unified by an iterative "difference map" constructed from a pair of elementary projections and a single real parameter $\beta$. For the standard application in optics, where the two projections implement Fourier modulus and object support constraints respectively, the difference map reproduces the "hybrid" form of Fienup's input-output map for $\beta = 1$. Other values of $\beta$ are equally effective in retrieving phases but have no input-output counterparts. The geometric construction of the difference map illuminates the distinction between its fixed points and the recovered object, as well as the mechanism whereby stagnation is avoided. When support constraints are replaced by object histogram or atomicity constraints, the difference map lends itself to crystallographic phase retrieval. Numerical experiments with synthetic data suggest that structures with hundreds of atoms can be solved.

A Generalization of Obreshkoff-Ehrlich Method for Multiple Roots of Polynomial Equations

Abstract: In this paper we develop a new method which is a generalization of the Obreshkoff -Ehrlich method for the cases of algebraic, trigonometric and exponential polynomials. This method has a cubic rate of convergence. It is efficient from the computational point of view and can be used for simultaneous finding all roots if the roots have known multiplicities. This new method in spite of the arbitrariness of multiplicities is of the same complexity as the methods for simultaneous finding all roots of simple roots. We do not use divided differences with multiple knots and this fact does not lead to calculation of derivatives of the given polynomial of higher order, but only of first ones.

Solving the difference initial-boundary value problems by the operator exponential method

Abstract: We suggest a modification of the operator exponential method for the numerical solving the difference linear initial boundary value problems. The scheme is based on the representation of the difference operator for given boundary conditions as the perturbation of the same operator for periodic ones. We analyze the error, stability and efficiency of the scheme for a model example of the one-dimensional operator of second difference.

A new algorithm for the volume of a convex polytope

Abstract: We provide two algorithms for computing the volume of a convex polytope with half-space representation {x>=0; Ax <=b} for some (m,n) matrix A and some m-vector b. Both algorithms have a O(n^m) computational complexity which makes them especially attractive for large n and relatively small m when the other methods with O(m^n) complexity fail. The methodology which differs from previous existing methods uses a Laplace transform technique that is well-suited to the half-space representation of the polytope.

Hybrid-Cubic-Rational Semi-Lagrangian Method with the Optimal Mixing

Abstract: A semi-Lagrangian method for advection equation with hybrid cubic-rational interpolation is introduced. In the present method, the spatial profile of physical quantities is interpolated with a combination of a cubic and a rational function. For achieving both high accuracy and convexity preserving of solution, the two functions are mixed in the optimal ratio which is given theoretically. Accuracy and validity of this method is demonstrated with some numerical experiments.

A semi-numerical computation for the added mass coefficients of an oscillating hemi-sphere at very low and very high frequencies

Abstract: A floating hemisphere under forced harmonic oscillation at very high and very low frequencies is considered. The problem is reduced to an elliptic one, that is, the Laplace operator in the exterior domain with standard Dirichlet and Neumann boundary conditions, so the flow problem is simplified to standard ones, with well known analytic solutions in some cases. The general procedure is based in the use of spherical harmonics and its derivation is based on a physics insight. The results can be used to test the accuracy achieved by numerical codes as, for example, by finite elements or boundary elements.

In search of an evolutionary coding style

Abstract: In the near future, all the human genes will be identified But understanding the functions coded in the genes is a much harder problem For example, by using block entropy, one has that the DNA code is closer to a random code then written text, which in turn is less ordered then an ordinary computer code; see \cite{schmitt} Instead of saying that the DNA is badly written, using our programming standards, we might say that it is written in a different style -- an evolutionary style We will suggest a way to search for such a style in a quantified manner by using an artificial life program, and by giving a definition of general codes and a definition of style for such codes

Approximate rational arithmetics and arbitrary precision computations

Abstract: We describe an approximate rational arithmetic with round-off errors (both absolute and relative) controlled by the user. The rounding procedure is based on the continued fraction expansion of real numbers. Results of computer experiments are given in order to compare efficiency and accuracy of different types of approximate arithmetics and rounding procedures.