Fast numerical solution of fredholm integral equations with stationary kernels

A fast recursive matrix method for the numerical solution of Fredholm integral equations with stationary kernels is derived. IfN denotes the number of nodal points, the complexity of the algorithm isO(N ²), which should be compared toO(N ³) for conventional algorithms for solving such problems. The method is related to fast algorithms for inverting Toeplitz matrices.Applications to equations of the first and second kind as well as miscellaneous problems are discussed and illustrated with numerical examples. These show that the theoretical improvement in efficiency is indeed obtained, and that no problems with numerical stability or accuracy are encountered.     

Diamond a sorting method for vector machines

A method for parallel sorting in vector machines is described. It is shown that it has complexityO(N(logN)²).     

Convergence analysis of the LDG method applied to singularly perturbed problems

Considering a two‐dimensional singularly perturbed convection–diffusion problem with exponential boundary layers, we analyze the local discontinuous Galerkin (DG) method that uses piecewise bilinear polynomials on Shishkin mesh. A convergence rate O(N^‐1 lnN) in a DG‐norm is established under the regularity assumptions, while the total number of mesh points is O(N²). The rate of convergence is uniformly valid with respect to the singular perturbation parameter ε. Numerical experiments indicate that the theoretical error estimate is sharp. © 2012 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 2013     

Improved approximations of the solution and derivatives to a singularly perturbed reaction-diffusion equation based on the solution decomposition method

In the case of the Dirichlet problem for a singularly perturbed ordinary differential reaction-diffusion equation, a new approach is used to the construction of finite difference schemes such that their solutions and their normalized first- and second-order derivatives converge in the maximum norm uniformly with respect to a perturbation parameter ɛ ∈(0, 1]; the normalized derivatives are ɛ-uniformly bounded. The key idea of this approach to the construction of ɛ-uniformly convergent finite difference schemes is the use of uniform grids for solving grid subproblems for the regular and singular components of the grid solution. Based on the asymptotic construction technique, a scheme of the solution decomposition method is constructed such that its solution and its normalized first- and second-order derivatives converge ɛ-uniformly at the rate of O(N ⁻²ln² N), where N + 1 is the number of points in the uniform grids. Using the Richardson technique, an improved scheme of the solution decomposition method is constructed such that its solution and its normalized first and second derivatives converge ɛ-uniformly in the maximum norm at the same rate of O(N ⁻⁴ln⁴ N).     

On Shellsort and the frobenius problem

A boundO(N ^1+1/k ) for the running time of Shellsort, withO(logN) passes, is proved very simply by application of a Frobenius basis withk elements.     

A fast algorithm for the evaluation of heat potentials

Numerical methods for solving the heat equation via potential theory have been hampered by the high cost of evaluating heat potentials. When M points are used in the discretization of the boundary and N time steps are computed, an amount of work of the order O(N²M²) has traditionally been required. In this paper, we present an algorithm which requires an amount of work of the order O(NM), and we observe speedups of five orders of magnitude for large-scale problems. Thus, the method makes it possible to solve the heat equation by potential theory in practical situations.     

Discrete fourier transform computation using prime Ramanujan numbers

Ramanujan numbers were introduced in [2] to implement discrete fourier transform (DFT) without using any multiplication operation. Ramanujan numbers are related to π and integers which are powers of 2. If the transform sizeN, is a Ramanujan number, then the computational complexity of the algorithms used for computing isO(N ²) addition and shift operations, and no multiplications. In these algorithms, the transform can be computed sequentially with a single adder inO(N ²) addition times. Parallel implementation of the algorithm can be executed inO(N) addition times, withO(N) number of adders. Some of these Ramanujan numbers of order-2 are related to the Biblical and Babylonian values of π [1]. In this paper, we analytically obtain upper bounds on the degree of approximation in the computation of DFT if JV is a prime Ramanujan number.     

Fast Computation of High‐Frequency Dirichlet Eigenmodes via Spectral Flow of the Interior Neumann‐to‐Dirichlet Map

We present a new algorithm for numerical computation of large eigenvalues and associated eigenfunctions of the Dirichlet Laplacian in a smooth, star‐shaped domain in ?^d, d ≥ 2. Conventional boundary‐based methods require a root search in eigenfrequency k, hence take O(N³) effort per eigenpair found, where N = O(k^d?1) is the number of unknowns required to discretize the boundary. Our method is O(N) faster, achieved by linearizing with respect to k the spectrum of a weighted interior Neumann‐to‐Dirichlet (NtD) operator for the Helmholtz equation. Approximations to the square roots k_j of all O(N) eigenvalues lying in [k ? ?, k], where ? = O(1), are found with O(N³) effort. We prove an error estimate with C independent of k. We present a higher‐order variant with eigenvalue error scaling empirically as O(?⁵) and eigenfunction error as O(?³), the former improving upon the “scaling method” of Vergini and Saraceno. For planar domains (d = 2), with an assumption of absence of spectral concentration, we also prove rigorous error bounds that are close to those numerically observed. For d = 2 we compute robustly the spectrum of the NtD operator via potential theory, Nyström discretization, and the Cayley transform. At high frequencies (400 wavelengths across), with eigenfrequency relative error 10^?10, we show that the method is 10³ times faster than standard ones based upon a root search. © 2014 Wiley Periodicals, Inc.     

On Computing the Subset Graph of a Collection of Sets

Let a given collection of sets have size N measured by the sum of the cardinalities. Yellin and Jutla presented an algorithm which constructed the partial order induced by the subset relation (a “subset graph”) in a claimed O(N²/log N) operations over a dictionary ADT, and exhibited a collection whose subset graph had Θ(N²/log² N) edges. This paper first establishes a matching upper bound on the number of edges in a subset graph. It also presents a finer analysis of the algorithm, which confirms the claimed upper bound and shows it to be tight. A simple implementation requiring O(1) bit-parallel operations per ADT operation is presented, along with a variant of the algorithm with an implementation requiring O(N²/log N) RAM operations.     

Fast fitting of radial basis functions: Methods based on preconditioned GMRES iteration

Solving large radial basis function (RBF) interpolation problems with non‐customised methods is computationally expensive and the matrices that occur are typically badly conditioned. For example, using the usual direct methods to fit an RBF with N centres requires O(N ²) storage and O(N ³) flops. Thus such an approach is not viable for large problems with N ≥10,000. In this paper we present preconditioning strategies which, in combination with fast matrix–vector multiplication and GMRES iteration, make the solution of large RBF interpolation problems orders of magnitude less expensive in storage and operations. In numerical experiments with thin‐plate spline and multiquadric RBFs the preconditioning typically results in dramatic clustering of eigenvalues and improves the condition numbers of the interpolation problem by several orders of magnitude. As a result of the eigenvalue clustering the number of GMRES iterations required to solve the preconditioned problem is of the order of 10-20. Taken together, the combination of a suitable approximate cardinal function preconditioner, the GMRES iterative method, and existing fast matrix–vector algorithms for RBFs [4,5] reduce the computational cost of solving an RBF interpolation problem to O(N) storage, and O(N \log N) operations. This revised version was published online in June 2006 with corrections to the Cover Date.     

A mixed‐type Galerkin variational formulation and fast algorithms for variable‐coefficient fractional diffusion equations

We consider the variable‐coefficient fractional diffusion equations with two‐sided fractional derivative. By introducing an intermediate variable, we propose a mixed‐type Galerkin variational formulation and prove the existence and uniqueness of the variational solution over . On the basis of the formulation, we develop a mixed‐type finite element procedure on commonly used finite element spaces and derive the solvability of the finite element solution and the error bounds for the unknown and the intermediate variable. For the Toeplitz‐like linear system generated by discretization, we design a fast conjugate gradient normal residual method to reduce the storage from O(N²) to O(N) and the computing cost from O(N³) to O(NlogN). Numerical experiments are included to verify our theoretical findings. Copyright © 2017 John Wiley & Sons, Ltd.     

Efficient Algorithms for Finding the Maximum Number of Disjoint Paths in Grids

In a rectangular grid, given two sets of nodes, (sources) and (sinks), of size each, the disjoint paths (DP) problem is to connect as many nodes in to the nodes in using a set of “disjoint” paths. (Both edge-disjoint and vertex-disjoint cases are considered in this paper.) Note that in this DP problem, a node in can be connected to any node in . Although in general the sizes of and do not have to be the same, algorithms presented in this paper can also find the maximum number of disjoint paths pairing nodes in and . We use the network flow approach to solve this DP problem. By exploiting all the properties of the network, such as planarity and regularity of a grid, integral flow, and unit capacity source/sink/flow, we can optimally compress the size of the working grid (to be defined) from O(N²) to O(N^1.5) and solve the problem in O(N^2.5) time for both the edge-disjoint and vertex-disjoint cases, an improvement over the straightforward approach which takes O(N³) time.     

Wave atoms and time upscaling of wave equations

We present a new geometric strategy for the numerical solution of hyperbolic wave equations in smoothly varying, two-dimensional time-independent periodic media. The method consists in representing the time-dependent Green’s function in wave atoms, a tight frame of multiscale, directional wave packets obeying a precise parabolic balance between oscillations and support size, namely wavelength ~(diameter).² Wave atoms offer a uniquely structured representation of the Green’s function in the sense that

A fast transform for spherical harmonics

Spherical harmonics arise on the sphere S² in the same way that the (Fourier) exponential functions {e^ik}k arise on the circle. Spherical harmonic series have many of the same wonderful properties as Fourier series, but have lacked one important thing: a numerically stable fast transform analogous to the Fast Fourier Transform (FFT). Without a fast transform, evaluating (or expanding in) spherical harmonic series on the computer is slow—for large computations probibitively slow. This paper provides a fast transform.For a grid ofO(N²) points on the sphere, a direct calculation has computational complexityO(N⁴), but a simple separation of variables and FFT reduce it toO(N³) time. Here we present algorithms with timesO(N^5/2 log N) andO(N²(log N)²). The problem quickly reduces to the fast application of matrices of associated Legendre functions of certain orders. The essential insight is that although these matrices are dense and oscillatory, locally they can be represented efficiently in trigonometric series.     

Optimized look-ahead recurrences for adjacent rows in the Padé table

A new look-ahead algorithm for recursively computing Padé approximants is introduced. It generates a subsequence of the Padé approximants on two adjacent rows (defined by fixed numerator degree) of the Padé table. Its two basic versions reduce to the classical Levinson and Schur algorithms if no look-ahead is required. The new algorithm can be viewed as a combination of the look-ahead sawtooth and the look-ahead Levinson and Schur algorithms that we proposed before, but now the look-ahead step size is minimal (as in the sawtooth version) and the computational costs are as low as in the least expensive competing algorithms (including our look-ahead Levinson and Schur algorithms). The underlying recurrences link well-conditioned basic pairs,i.e., pairs of sufficiently different neighboring Padé forms.The algorithm can be used to solve Toeplitz systems of equationsTx = b. In this application it comes in several versions: anO(N ²) Levinson-type form, anO(N ²) Schur-type form, and a superfastO(N log² N) Schur-type version. As an option of the first two versions, the corresponding block LDU decompositions ofT ^–1 orT, respectively, can be found.     

An effective direct solution method for certain boundary element equations in 3D

The boundary element method for the Dirichlet problem in a three-dimensional rotational domain leads to a system of linear equations with a full dense matrix having a special block structure. A direct solution method for such systems is presented, which requires O(N^3/2 ln N) arithmetical operations only, using a Fast Fourier Transformation (FFT), where N denotes the number of unknowns on the boundary surface.     

On a proposed divide-and-conquer minimal spanning tree algorithm

In a recent paper published in this journal, R. Chang and R. Lee purport to devise anO(N logN) time minimal spanning tree algorithm forN points in the plane that is based on a divide-and-conquer strategy using Voronoi diagrams. In this brief note, we present families of problem instances to show that the Chang-Lee worst-case timing analysis is incorrect, resulting in a time bound ofO(N ² logN). Since it is known that alternate, trulyO(N logN) time algorithms are available anyway, the general utility of the Chang-Lee algorithm is questionable.This author's research is supported in part by the Washington State Technology Center and by the National Science Foundation under grants ECS-8403859 and MIP-8603879.     

Character table of a controlling association scheme defined by the general orthogonal groupO(3,q).

LetN _m (q) be the set of nonisotropic lines in the vector space of dimensionm over a finite field of orderq. In a paper by Bannai, Hao, Song and Wei, it was shown that the association scheme character tableP(Sp(2n, q),N _2n (q)), withn 3 andq odd, is controlled byP(Sp(4,q),N ₄(q)) which is in turn controlled byP(O(3,q),O(3,q)/O ⁺ (2,q)). Our purpose in this paper is to compute the entries in the character tableP(O(3,q), O(3,q)/O ⁺ (2,q)) explicitly, which is left open in that paper.     

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 N ^d , 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 ε ⁻¹)-complexity bound for their approximation by low-rank 2-(dlog N) quantics tensors up to the accuracy ε>0. In the case ε=O(N ^−α ), α>0, our approach leads to the quantics tensor numerical method in dimension d, with the nearly optimal asymptotic complexity O(d/αlog ² ε ⁻¹). 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.     

Error analysis of the finite-strip method for parabolic equations

The finite-strip method (FSM) is a hybrid technique which combines spectral and finite-element methods. Finite-element approximations are made for each mode of a finite Fourier series expansion. The Galerkin formulated method is set apart from other weighted-residual techniques by the selection of two types of basis functions, a piecewise linear interpolating function and a trigonometric function. The efficiency of the FSM is due in part to the orthogonality of the complex exponential basis: The linear system which results from the weak formulation is decoupled into several smaller systems, each of which may be solved independently. An error analysis for the FSM applied to time-dependent, parabolic partial differential equations indicates the numerical solution error is O(h² + M^?r). M represents the Fourier truncation mode number and h represents the finite-element grid mesh. The exponent r ≥ 2 increases with the exact solution smoothness in the respective dimension. This error estimate is verified computationally. Extending the result to the finite-layer method, where a two-dimensional trigonometric basis is used, the numerical solution error is O(h² + M^?r + N^?q). The N and q represent the trucation mode number and degree of exact solution smoothness in the additional dimension. © 1993 John Wiley & Sons, Inc.