A note on the solution of not balanced banded Toeplitz systems

A direct algorithm is presented for the solution of linear systems having banded Toeplitz coefficient matrix with unbalanced bandwidths. It is derived from the cyclic reduction algorithm, it makes use of techniques based on the displacement rank and it relies on the Morrison–Sherman–Woodbury formula. The algorithm always equals and sometimes outperforms the already known direct ones in terms of asymptotic computational cost. The case where the coefficient matrix is a block banded block Toeplitz matrix in block Hessenberg form is analyzed as well. The algorithm is numerically stable if applied to M‐matrices that are point diagonally dominant by columns. Copyright © 2007 John Wiley & Sons, Ltd.     

Numerically stable deflation of hessenberg and symmetric tridiagonal matrices

While numerically stable techniques have been available for deflating a fulln byn matrix, no satisfactory finite technique has been known which preserves Hessenberg form. We describe a new algorithm which explicitly deflates a Hessenberg matrix in floating point arithmetic by means of a sequence of plane rotations. When applied to a symmetric tridiagonal matrix, the deflated matrix is again symmetric tridiagonal. Repeated deflation can be used to find an orthogonal set of eigenvectors associated with any selection of eigenvalues of a symmetric tridiagonal matrix.     

On the Linear Stability of Splitting Methods

A comprehensive linear stability analysis of splitting methods is carried out by means of a 2×2 matrix K(x) with polynomial entries (the stability matrix) and the stability polynomial p(x) (the trace of K(x) divided by two). An algorithm is provided for determining the coefficients of all possible time-reversible splitting schemes for a prescribed stability polynomial. It is shown that p(x) carries essentially all the information needed to construct processed splitting methods for numerically approximating the evolution of linear systems. By conveniently selecting the stability polynomial, new integrators with processing for linear equations are built which are orders of magnitude more efficient than other algorithms previously available. This paper is dedicated to Arieh Iserles on the occasion of his 60th anniversary.     

Parallel Implementation of Stability Analysis of Difference Schemes with MATHEMATICA

We consider a parallel algorithm for investigating the stability of the schemes of the finite-difference and finite-volume methods that approximate the two-dimensional Euler equations of compressible fluid on a curvilinear grid. The algorithm is implemented with the aid of the computer algebra system Mathematica 3.0. We apply a two-level parallelization process. At the first level, the symbolic computation of the amplification matrix is parallelized by a parallel computation of the matrix rows on different processors. At the second level, the values of the coordinates of points of the stability-region boundary are computed numerically. For the communication between the workstations, we apply a special program, LaunchSlave, which uses the MathLink communication protocol. Examples of application of the proposed parallel symbolic/numerical algorithm are presented. Bibliography: 15 titles.     

Spline smoothing under constraints on derivatives

An efficient algorithm for computing a smoothing polynomial splines under inequality constraints on derivatives is introduced where both order and breakpoints ofs can be prescribed arbitrarily. By using the B-spline representation ofs, the original semi-infinite constraints are replaced by stronger finite ones, leading to a least squares problem with linear inequality constraints. Then these constraints are transformed into simple box constraints by an appropriate substitution of variables so that efficient standard techniques for solving such problems can be applied. Moreover, the smoothing term commonly used is replaced by a cheaply computable approximation. All matrix transformations are realized by numerically stable Givens rotations, and the band structure of the problem is exploited as far as possible.     

A parallel radix‐4 block cyclic reduction algorithm

A conventional block cyclic reduction algorithm operates by halving the size of the linear system at each reduction step, that is, the algorithm is a radix‐2 method. An algorithm analogous to the block cyclic reduction known as the radix‐q partial solution variant of the cyclic reduction (PSCR) method allows the use of higher radix numbers and is thus more suitable for parallel architectures as it requires fever reduction steps. This paper presents an alternative and more intuitive way of deriving a radix‐4 block cyclic reduction method for systems with a coefficient matrix of the form tridiag{ ? I,D, ? I}. This is performed by modifying an existing radix‐2 block cyclic reduction method. The resulting algorithm is then parallelized by using the partial fraction technique. The parallel variant is demonstrated to be less computationally expensive when compared to the radix‐2 block cyclic reduction method in the sense that the total number of emerging subproblems is reduced. The method is also shown to be numerically stable and equivalent to the radix‐4 PSCR method. The numerical results archived correspond to the theoretical expectations. Copyright © 2013 John Wiley & Sons, Ltd.     

Numerical analysis of the method of fundamental solution for harmonic problems in annular domains

In this study, we investigate the application of the method of fundamental solutions (MFS) to the Dirichlet problem for Laplace's equation in an annular domain. We examine the properties of the resulting coefficient matrix and its eigenvalues. The convergence of the method is proved for analytic boundary data. An efficient matrix decomposition algorithm using fast Fourier transforms (FFTs) is developed for the computation of the MFS approximation. We also tested the algorithm numerically on several problems confirming the theoretical predictions. © 2005 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2005     

An SQP algorithm with cautious updating criteria for nonlinear degenerate problems

An efficient SQP algorithm for solving nonlinear degenerate problems is proposed in the paper. At each iteration of the algorithm, a quadratic programming subproblem, which is always feasible by introducing a slack variable, is solved to obtain a search direction. The steplength along this direction is computed by employing the 1∞ exact penalty function through Armijo-type line search scheme. The algorithm is proved to be convergent globally under mild conditions.     

A binary powering Schur algorithm for computing primary matrix roots

An algorithm for computing primary roots of a nonsingular matrix A is presented. In particular, it computes the principal root of a real matrix having no nonpositive real eigenvalues, using real arithmetic. The algorithm is based on the Schur decomposition of A and has an order of complexity lower than the customary Schur based algorithm, namely the Smith algorithm.     

Robust Parameter Estimation for Stochastic Differential Equations

We consider the estimation of parameters in stochastic differential equations (SDEs). The problem is treated in the setting of nonlinear filtering theory with a degenerate diffusion matrix. A robust stochastic Feynman–Kac representation for solutions of SDEs of Zakai-type is derived. It is verified that these solutions are conditional densities for the conditional measures defined by degenerate filtering problems. We show that the corresponding estimator for the parameters is robust in the following sense: It depends continuously on both the measurement path and on the intensity of the measurement noise. An algorithm based on a Monte-Carlo approach is given for the practical application of the estimator, and numerical results are reported. Mathematics Subject Classifications (2000) Primary: 62M05, 62M20; secondary: 62F15.     

A method for computing the Perron root for primitive matrices

Following the Perron theorem, the spectral radius of a primitive matrix is a simple eigenvalue. It is shown that for a primitive matrix A, there is a positive rank one matrix X such that B = A ° X , where ° denotes the Hadamard product of matrices, and such that the row (column) sums of matrix B are the same and equal to the Perron root. An iterative algorithm is presented to obtain matrix B without an explicit knowledge of X. The convergence rate of this algorithm is similar to that of the power method but it uses less computational load. A byproduct of the proposed algorithm is a new method for calculating the first eigenvector.     

The cauchy problem of a degenerate parabolic equation

This note is concerned with the existence of a weak solution for a degenerate Cauchy problem of parabolic type in then-dimensional spaceR ⁿ. The degenerate property is in the sense that the matrix (a _ij(t,x)) involved in the differential operator is not necessarily positive definite. The essential idea is the construction of a suitable function spaceH and to prove the existence of a weak solution inH.     

A saxpy formulation for plane rotations

An algorithm for pre- or post-multiplication of a matrix by a plane rotation, using only three vector saxpy operations instead of the four vector operations usually considered necessary, is described. No auxiliary storage for overwriting is required. The method is shown to be numerically stable.     

An efficient and numerically correct algorithm for the 2D convex hull problem

An efficient and numerically correct program called FastHull for computing the convex hulls of finite point sets in the plane is presented. It is based on the Akl-Toussaint algorithm and the MergeHull algorithm. Numerical correctness of the FastHull procedure is ensured by using special routines for interval arithmetic and multiple precision arithmetic. The FastHull algorithm guaranteesO(N logN) running time in the worst case and has linear time performance for many kinds of input patterns. It appears that the FastHull algorithm runs faster than any currently known 2D convex hull algorithm for many input point patterns.     

Numerical Evaluation of Integrals Involving the β‐PDF

The fast and accurate evaluation of expected values involving the probabilistic β‐distributions poses severe problems to standard numerical integration schemes. An efficient algorithm to evaluate such integrals is presented based on approximation and analytical evaluation rather than numerical integration. Starting from an extension of the 2‐parameter family of β‐distributions a criterion is derived to assess the correctness of any integration scheme in the numerically demanding limiting cases of this family.     

Construction of Matrices with Prescribed Singular Values and Eigenvalues

Two issues concerning the construction of square matrices with prescribe singular values an eigenvalues are addressed. First, a necessary and sufficient condition for the existence of an n × n complex matrix with n given nonnegative numbers as singular values an m ( n) given complex numbers to be m of the eigenvalues is determined. This extends the classical result of Weyl and Horn treating the case when m = n. Second, an algorithm is given to generate a triangular matrix with prescribe singular values an eigenvalues. Unlike earlier algorithms, the eigenvalues can be arranged in any prescribe order on the diagonal. A slight modification of this algorithm allows one to construct a real matrix with specified real an complex conjugate eigenvalues an specified singular values. The construction is done by multiplication by diagonal unitary matrices, permutation matrices and rotation matrices. It is numerically stable and may be useful in developing test software for numerical linear algebra packages.     

QR factorization of Toeplitz matrices

Summary This paper presents a new algorithm for computing theQR factorization of anm×n Toeplitz matrix inO(mn) operations. The algorithm exploits the procedure for the rank-1 modification and the fact that both principal (m–1)×(n–1) submatrices of the Toeplitz matrix are identical. An efficient parallel implementation of the algorithm is possible.     

The stability of linear potential gyroscopic systems

The stability of linear potential systems with a degenerate matrix of gyroscopic forces is investigated. Particular attention is devoted to the case of three degrees of freedom. In a development of existing results [Kozlov VV. Gyroscopic stabilization and parametric resonance. Prikl. Mat. Mekh. 2001; 65(5): 739–745], the sufficient conditions for gyroscopic stability are obtained. An algorithm for applying these conditions is proposed using the example of the problem of the motion of two mutually gravitating bodies, each of them being modelled by two equal point masses, connected by weightless inextensible rods.     

Finding all maximal efficient faces in multiobjective linear programming

An algorithm for finding the whole efficient set of a multiobjective linear program is proposed. From the set of efficient edges incident to a vertex, a characterization of maximal efficient faces containing the vertex is given. By means of the lexicographic selection rule of Dantzig, Orden and Wolfe, a connectedness property of the set of dual optimal bases associated to a degenerate vertex is proved. An application of this to the problem of enumerating all the efficient edges incident to a degenerate vertex is proposed. Our method is illustrated with numerical examples and comparisons with Armand—Malivert's algorithm show that this new algorithm uses less computer time.