High-accuracy approximation of the complex probability function by Fourier expansion of exponential multiplier

The real K(x,y) and imaginary L(x,y) parts of the complex probability function are approximated as rapidly convergent series, based on the Fourier expansion of the exponential multiplier. This approach provides rapid and accurate calculations of the Voigt and complex error functions in the most challenging Humlí?ek regions 3 and 4.     

Computing Bessel functions of the second kind in extreme parameter regimes

A useful method of computing the integral order Bessel functions of the second kind Y_n(x+iy) when either, the absolute value of the real part, or the imaginary part of the argument z=x+iy is small, is described. This method is based on computing the Bessel functions for extreme parameter regimes when x∼0 (or y∼0) and is useful because a number existing algorithms and methods fail to give correct results for small x or small y. The approximating equations are derived by expanding the Bessel function in Taylor series, are tested and discussed. The present work is a continuation of the previous one conducted in regard to the Bessel function of the first kind. The results of our formalism are compared to the available existing numerical methods used in Mathematica, IMSL, MATLAB, and the Amos library. Our numerical method is easy to implement, efficient, and produces reliable results. In addition, this method reduces the computation of the Bessel functions of the second complex argument to that of real argument which simplify the computation considerably.     

Cylindrical decomposition for systems transcendental in the first variable

We present a generalization of the Cylindrical Algebraic Decomposition (CAD) algorithm to systems of equations and inequalities in functions of the form p(x,f₁(x),…,f_m(x),y₁,…,y_n), where p∈Q[x,t₁,…,t_m,y₁,…,y_n] and f₁(x),…,f_m(x) are real univariate functions such that there exists a real root isolation algorithm for functions from the algebra Q[x,f₁(x),…,f_m(x)]. In particular, the algorithm applies when f₁(x),…,f_m(x) are real exp-log functions or tame elementary functions.     

Predictive complexity and information

The notions of predictive complexity and of corresponding amount of information are considered. Predictive complexity is a generalization of Kolmogorov complexity which bounds the ability of any algorithm to predict elements of a sequence of outcomes. We consider predictive complexity for a wide class of bounded loss functions which are generalizations of square-loss function. Relations between unconditional KG(x) and conditional KG(x|y) predictive complexities are studied. We define an algorithm which has some “expanding property”. It transforms with positive probability sequences of given predictive complexity into sequences of essentially bigger predictive complexity. A concept of amount of predictive information IG(y:x) is studied. We show that this information is noncommutative in a very strong sense and present asymptotic relations between values IG(y:x), IG(x:y), KG(x) and KG(y).     

Curve generation of implicit functions by incremental computers

The conventional numerical solution of an implicit function f(x, y) = 0 is substantially complicated for calculating by any computer. We propose a new method representing the argument of the implicit function as a unary function of a parameter, t, if the continuous and unique solution of f(x, y) = 0 exists. The total differential $\text{d}\text{f}\text{d}\text{t}$ constitutes simultaneous differential equations of which the solution about x and y is unique. The Newton-Raphson method must be used to calculate the values near singular points of an implicit function and then the sign of dt has to be decided according to four special cases. Incremental computers are suitable for curve generation of implicit functions by the new method, because the incremental computer can perform more complex algorithms than the analog computer and can calculate faster than the digital computer. This method is easily applicable to curve generation in three-dimensional space.     

An algorithm for regular and irregular Coulomb and Bessel functions of real order to machine accuracy

We describe an algorithm to evaluate a wide class of functions and their derivatives, to extreme precision (25–30S) if required, which does not use any function calls other than square root. The functions are the Coulomb functions of positive argument (F_λ(x, η), G_λ(x, η), x > 0, η, λ real) and hence, as special cases with η = 0, the cylindrical Bessel functions (J_μ(x), Y_μ(x), x > 0, μ real), the spherical Bessel functions (i_λ(x), y_λ(x), x > 0, λ real), Airy functions of negative argument Ai(-x), Bi(-x) and others. The present method has a number of attractive features: both the regular and irregular solution are calculated, all others of the functions can be produced from a specified minimum (not necessarily zero) to a specified maximum, functions of a single order can be found without all of the orders from zero, the derivatives of the functions arise naturally in the solution and are readily available, the results are available to different precisions from the same subroutine (in contrast to rational approximation techniques) and the methods can be used for estimating final accuracies. In addition, the sole constant required in the algorithm is π, no precalculated arrays of coefficients are needed, and the final accuracy is not dependent on that of other subroutines. The method works most efficiently in the region x ≈ 0.5 to x ≈ 1000 but outside this region the results are still reliable, even though the number of iterations within the subroutine rises. Even in these more asymptotic regions the unchanged algorithm can be used with known accuracy to test other specific subroutines more appropriate to these regions. The algorithm uses the recursion relations satisfied by the Coulomb functions and contains a significant advance over Miller's method for evaluating the ratio of successive minimal solutions (F_λ+1/F_λ). It relies on the evaluation of two continued fractions and no infinite series is required for normalisation: instead the Wronskian is used.     

Efficient algorithms to compute Hankel transforms using wavelets

The aim of the paper is to propose two efficient algorithms for the numerical evaluation of Hankel transform of order ν, ν>−1 using Legendre and rationalized Haar (RH) wavelets. The philosophy behind the algorithms is to replace the part xf(x) of the integrand by its wavelet decomposition obtained by using Legendre wavelets for the first algorithm and RH wavelets for the second one, thus representing Fν(y) as a Fourier-Bessel series with coefficients depending strongly on the input function xf(x) in both the cases. Numerical evaluations of test functions with known analytical Hankel transforms illustrate the proposed algorithms.     

A recursive and a grammatical characterization of the exponential-time languages

We characterize the class of all languages which are acceptable in exponential time by means of recursive and grammatical methods. (i) The class of all languages which are acceptable in exponential time is uniquely characterized by the class of all (0-1)-functions which can be generated, starting with the initial functions of the Grzegorczyk-class $\text{E}$², by means of subtitution and limited recursion of the form f(x, y + 1) = h(x, y), f(x, y), f(x, l(x, y))), l(x, y) ? y. (ii) The class of all languages which are acceptable in exponential time is equal to the class of all languages generated by context-sensitive grammars with context-free control sets.     

Quadrature rules for the integration of the product of two oscillatory functions with different frequencies

We construct quadrature rules for the efficient computation of the integral of a product of two oscillatory functions y₁(x) and y₂(x), where , and the functions f_i,j(x) are smooth. The weights are evaluated by the exponential fitting technique of Ixaru [Comput. Phys. Comm. 105 (1997) 1-19], which is now extended to cover the case of two frequencies. We give a numerical illustration on how the new rules compare for accuracy with the one-frequency dependent rules and with the classical ones.     

Efficient computation of the cylindrical Bessel functions of complex argument

An algorithm that generates the cylindrical Bessel function very accurately for a wide range of complex arguments has been developed by Mason. The Mason algorithm consists of four different methods that apply to different portions of the complex plane. Experience with the Floating Point Systems FPS-364 minisupercomputer indicates several ways by which these methods can be made more efficient. Specific improvements relate to: 1) the method for determination of the point where backward recursion is initiated for the Bessel functions of the first kind; 2) the way that the Bessel functions of the first and second kind are normalized when |y| < 5 and |x| ⩽ 20; and 3) the extent that asymptotic expansions are used when |x| > 20 and |y| < 5. The first and third modifications will result in increased efficiency for all architectures. The second modification will be of value for many, but probably not all, architectures.     

f1: a code to compute Appell's F1 hypergeometric function

In this work we present the FORTRAN code to compute the hypergeometric function F₁(α,β₁,β₂,γ,x,y) of Appell. The program can compute the F₁ function for real values of the variables {x,y}, and complex values of the parameters {α,β₁,β₂,γ}. The code uses different strategies to calculate the function according to the ideas outlined in [F.D. Colavecchia et al., Comput. Phys. Comm. 138 (1) (2001) 29].

Program summary

Title of the program: f1Catalogue identifier: ADSJProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADSJProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions: noneComputers: PC compatibles, SGI Origin2∗Operating system under which the program has been tested: Linux, IRIXProgramming language used: Fortran 90Memory required to execute with typical data: 4 kbytesNo. of bits in a word: 32No. of bytes in distributed program, including test data, etc.: 52 325Distribution format: tar gzip fileExternal subprograms used: Numerical Recipes hypgeo [W.H. Press et al., Numerical Recipes in Fortran 77, Cambridge Univ. Press, 1996] or chyp routine of R.C. Forrey [J. Comput. Phys. 137 (1997) 79], rkf45 [L.F. Shampine and H.H. Watts, Rep. SAND76-0585, 1976].Keywords: Numerical methods, special functions, hypergeometric functions, Appell functions, Gauss functionNature of the physical problem: Computing the Appell F₁ function is relevant in atomic collisions and elementary particle physics. It is usually the result of multidimensional integrals involving Coulomb continuum states.Method of solution: The F₁ function has a convergent-series definition for |x|<1 and |y|<1, and several analytic continuations for other regions of the variable space. The code tests the values of the variables and selects one of the precedent cases. In the convergence region the program uses the series definition near the origin of coordinates, and a numerical integration of the third-order differential parametric equation for the F₁ function. Also detects several special cases according to the values of the parameters.Restrictions on the complexity of the problem: The code is restricted to real values of the variables {x,y}. Also, there are some parameter domains that are not covered. These usually imply differences between integer parameters that lead to negative integer arguments of Gamma functions.Typical running time: Depends basically on the variables. The computation of Table 4 of [F.D. Colavecchia et al., Comput. Phys. Comm. 138 (1) (2001) 29] (64 functions) requires approximately 0.33 s in a Athlon 900 MHz processor.     

Bessel function of the first kind with complex argument

A new method of computing integral order Bessel functions of the first kind J_n(z) when either the absolute value of the real part or the imaginary part of the argument z = x + iy is small, is described. This method is based on computing the Bessel functions from asymptotic expressions when x∼ 0 (or y ∼ 0). These expansions are derived from the integral definition of Bessel functions. This method is necessary because some existing algorithms and methods fail to give correct results for small x small y. In addition, our overall method of computing Bessel functions of any order and argument is discussed and the logarithmic derivative is used in computing these functions. The starting point of the backward recurrence relations needed to evaluate the Bessel function and their logarithmic derivatives are investigated in order to obtain accurate numerical results. Our numerical method, together with established techniques of computing the Bessel functions, is easy to implement, efficient, and produces reliable results for all z.     

Gurland's ratio for the gamma function

We consider the ratio T(x, y) = г(x)г(y) / г²((x + y)/2) and its properties related to convexity, logarithmic convexity, Schur-convexity, and complete monotonicity. Several new bounds and asymptotic expansions for T are derived. Sharp bounds for the function x → x/(1 - e^−x) are presented, as well as bounds for the trigamma function. The results are applied to a problem related to the volume of the unit ball in Rⁿ and also to the problem of finding the inverse of the function x → T(1/x, 3/x), which is of importance in applied statistics.     

Some remarks on a conjecture in parameter adaptive control

A.S. Morse has raised the following question: Do there exist differentiable functions

$\text{f:R}{}^2\text{→ R}\text{and}\text{g:R}{}^2\text{→ R}$

with the property that for every nonzero real number λ and every (x₀, y₀) ∈ $\text{R}$² the solution (x(t),y(t)) of

$\text{x}\text{?}\text{(t) = x(t) + λf(x(t),y(t))}$

,

$\text{y}\text{?}\text{(t) = g(x(t),y(t))}$

,

$\text{x(0) = x}{}_0\text{, y(0) = y}{}_0$

, is defined for all t ? 0 and satisfies

$\text{lim}\text{t → + ∞}$

and y(t) is bounded on [0,∞)? We prove that the answer is yes, and we give explicit real analytic functions f and g which work. However, we prove that if f and g are restricted to be rational functions, the answer is no.     

Some procedures for function approximation based on the use of sample data and their application in heuristic methods for solving practical problems

Let f(x) be a member of a set of functions over a probability space. Samples of f(x) are 2-tuples (xⁱ,f(xⁱ) where xⁱ is a sample of the random variable X and f(xⁱ) is a sample of f(x) at x = xⁱ. Some procedures and analysis are presented for the approximation of such functions by systems of orthonormal functions. The approximations are based on the data samples. The analysis includes the case of error in the measurement of f(xⁱ). The properties of the expected square error in the approximation are examined for a number of different estimators for the coefficients in the expansion and these well-behaved and easily analyzed estimators are compared to those obtained using the method of least squares. The effectiveness of different sets of basis functions, those involved in the Karhunen-Loeve expansion and others, can be compared and an approach is suggested to adaptive basis selection in order to select that basis which is most efficient in approximating the particular function under examination. The connection between results and applications are discussed in the introduction and conclusion.     

Robust algorithms that locate local extrema of a function of one variable from interval measurement results: A remark

The problem of locating local maxima and minima of a function from approximate measurement results is vital for many physical applications: inspectral analysis, chemical species are identified by locating local maxima of the spectra; inradioastronomy, sources of celestial radio emission, and their subcomponents, are identified by locating local maxima of the measured brightness of the radio sky;elementary particles are identified by locating local maxima of the experimental curves. Since measurements are never absolutely precise, as a result of the measurements, we have aclass of possible functions. If we measuref(x _i ) with interval uncertainty, this class consists of all functionsf for whichf(x _i ) ε [y _i ??, y _i +?], wherey _i are the results of measuringf(x _i ), andε is the measurement accuracy. For this class, in [2], a linear-time algorithm was described. In real life, a measuring instrument can sometimes malfunction, leading to the so-calledoutliers, i.e., measurementsy _i that can be way offf(x _i ) (and thus do not restrict the actual valuesf(x _i ) at all). In this paper, we describerobust algorithms, i.e., algorithms that find the number of local extrema in the presence of possible outliers. These algorithms solve an important practical problem, but they are not based on any new mathematical results: they simply use algorithms from [2] and [3].     

Two-stage continuation algorithms for Bloch waves of Bose-Einstein condensates in optical lattices

Bloch waves of Bose-Einstein condensates (BEC) in optical lattices are extremum nonlinear eigenstates which satisfy the time-independent Gross-Pitaevskii equation (GPE). We describe an efficient Taylor predictor-Newton corrector continuation algorithm for tracing solution curves of parameter-dependent problems. Based on this algorithm, a novel two-stage continuation algorithm is developed for computing Bloch waves of 1D and 2D Bose-Einstein condensates (BEC) in optical lattices. We split the complex wave function into the sum of its real and imaginary parts. The original GPE becomes a couple of two nonlinear eigenvalue problems defined in the real domain with periodic boundary conditions. At the first stage we use the chemical potential μ as the continuation parameter. The Bloch wavenumber k(kx,ky), and the coefficient of the cubic term are treated as the second and third continuation parameters, respectively. Then we compute the Bloch bands/surfaces for the 1D/2D problem with linear counterparts. At the second stage we use μ and k/kx or ky as the continuation parameters simultaneously with two constraint conditions. The states without linear counterparts in the GPE can be obtained via states with linear counterparts. Numerical results are reported for both 1D and 2D problems.     

A sixth-order extrapolation method for special nonlinear fourth-order boundary value problems

A sixth-order convergent finite difference method is developed for the numerical solution of the special nonlinear fourth-order boundary value problem y^(iv)(x) = f(x, y), a < x < b, y(a) = A₀, y″(a) = B₀, y(b) = A₁ y′(b) = B₁, the simple-simple beam problem.The method is based on a second-order convergent method which is used on three grids, sixth-order convergence being obtained by taking a linear combination of the (second-order) numerical results calculated using the three individual grids.Special formulas are proposed for application to points of the discretization adjacent to the boundaries x = a and x= b, the first two terms of the local truncation errors of these formulas being the same as those of the second-order method used at the other points of each grid.Modifications to these two formulas are obtained for problems with boundary conditions of the form y(a) = A₀, y′(a) = C₀, y(b) = A₁, y′(b) = C₁, the clamped-clamped beam problem.The general boundary value problem, for which the differential equation is y^(iv)(x) = f(x, y, y′, y″, y‴), is also considered.     

Estimation of densities of probability and regression surfaces in one or two dimensions

The FORTRAN programs presented make it possible to build curves and surfaces of densities for the Lebesgue-measure, when one has a sample of n independent observations of a random variable in one or two dimensions and when this number n can be high (many thousands). The method uses kernel-estimators with varying window-parameters estimated via a modified maximum likelihood procedure. In the case of a three-dimensional variable, it is possible to estimate the function: (x, y) → E(Z/X> = x, Y = y) of conditional expectation. Studies based on simulations as well as on real data are pres ented.     

Asymptotic behaviour of the solutions of second order functional differential equations

This paper presents integral criteria to determine the asymptotic behaviour of the solutions of second order nonlinear differential equations of the type y^″(x)+q(x)f(y(x))=0, with q(x)>0 and f(y) odd and positive for y>0, as x tends to +∞. It also compares them with the results obtained by Chanturia (1975) in [11] for the same problem.