Symplectic Householder transformations for a QR-like decomposition,a geometric and algebraic approaches

The aim of this paper is to show how geometric and algebraic approaches lead us to a new symplectic elementary transformations: the 2-D symplectic Householder transformations. Their features are studied in details. Their interesting properties allow us to construct a new algorithm for computing a SR factorization. This algorithm is based only on these 2-D symplectic Householder transformations. Its new features are highlighted. The study shows that, in the symplectic case, the new algorithm is the corresponding one to the classical QR factorization algorithm, via the Householder transformations. Some numerical experiments are given.     

Numerical enclosure for each eigenvalue in generalized eigenvalue problem

An algorithm for enclosing all eigenvalues in generalized eigenvalue problem Ax=λBx is proposed. This algorithm is applicable even if A∈C^n×n is not Hermitian and/or B∈C^n×n is not Hermitian positive definite, and supplies nerror bounds while the algorithm previously developed by the author supplies a single error bound. It is proved that the error bounds obtained by the proposed algorithm are equal or smaller than that by the previous algorithm. Computational cost for the proposed algorithm is similar to that for the previous algorithm. Numerical results show the property of the proposed algorithm.     

An approximate decomposition algorithm for convex minimization

For nonsmooth convex optimization, Robert Mifflin and Claudia Sagastizábal introduce a VU-space decomposition algorithm in Mifflin and Sagastizábal (2005) [11]. An attractive property of this algorithm is that if a primal-dual track exists, this algorithm uses a bundle subroutine. With the inclusion of a simple line search, it is proved to be globally and superlinearly convergent. However, a drawback is that it needs the exact subgradients of the objective function, which is expensive to compute. In this paper an approximate decomposition algorithm based on proximal bundle-type method is introduced that is capable to deal with approximate subgradients. It is shown that the sequence of iterates generated by the resulting algorithm converges to the optimal solutions of the problem. Numerical tests emphasize the theoretical findings.     

Jacobi matrices for sums of weight functions

Orthogonal polynomials are conveniently represented by the tridiagonal Jacobi matrix of coefficients of the recurrence relation which they satisfy. LetJ ₁ andJ ₂ be finite Jacobi matrices for the weight functionsw ₁ andw ₂, resp. Is it possible to determine a Jacobi matrix \(\tilde J\) , corresponding to the weight functions \(\tilde w\) =w ₁+w ₂ using onlyJ ₁ andJ ₂ and if so, what can be said about its dimension? Thus, it is important to clarify the connection between a finite Jacobi matrix and its corresponding weight function(s). This leads to the need for stable numerical processes that evaluate such matrices. Three newO(n ²) methods are derived that “merge” Jacobi matrices directly without using any information about the corresponding weight functions. The first can be implemented using any of the updating techniques developed earlier by the authors. The second new method, based on rotations, is the most stable. The third new method is closely related to the modified Chebyshev algorithm and, although it is the most economical of the three, suffers from instability for certain kinds of data. The concepts and the methods are illustrated by small numerical examples, the algorithms are outlined and the results of numerical tests are reported.     

Parallel defect control

How can small-scale parallelism best be exploited in the solution of nonstiff initial value problems? It is generally accepted that only modest gains inefficiency are possible, and it is often the case that “fast” parallel algorithms have quite crude error control and stepsize selection components. In this paper we consider the possibility of using parallelism to improvereliability andfunctionality rather than efficiency. We present an algorithm that can be used with any explicit Runge-Kutta formula. The basic idea is to take several smaller substeps in parallel with the main step. The substeps provide an interpolation facility that is essentially free, and the error control strategy can then be based on a defect (residual) sample. If the number of processors exceeds (p ? 1)/2, wherep is the order of the Runge-Kutta formula, then the interpolant and the error control scheme satisfy very strong reliability conditions. Further, for a given orderp, the asymptotically optimal values for the substep lengths are independent of the problem and formula and hence can be computed a priori. Theoretical comparisons between the parallel algorithm and optimal sequential algorithms at various orders are given. We also report on numerical tests of the reliability and efficiency of the new algorithm, and give some parallel timing statistics from a 4-processor machine.     

Factoring matrices with a tree-structured sparsity pattern

Let A be a matrix whose sparsity pattern is a tree with maximal degree d_max. We show that if the columns of A are ordered using minimum degree on |A|+|A|^∗, then factoring A using a sparse LU with partial pivoting algorithm generates only O(d_maxn) fill, requires only O(d_maxn) operations, and is much more stable than LU with partial pivoting on a general matrix. We also propose an even more efficient and just-as-stable algorithm called sibling-dominant pivoting. This algorithm is a strict partial pivoting algorithm that modifies the column preordering locally to minimize fill and work. It leads to only O(n) work and fill. More conventional column pre-ordering methods that are based (usually implicitly) on the sparsity pattern of |A|^∗|A| are not as efficient as the approaches that we propose in this paper.     

Parallel domain decomposition procedures of improved D-D type for parabolic problems

Two parallel domain decomposition procedures for solving initial-boundary value problems of parabolic partial differential equations are proposed. One is the extended D-D type algorithm, which extends the explicit/implicit conservative Galerkin domain decomposition procedures, given in [5], from a rectangle domain and its decomposition that consisted of a stripe of sub-rectangles into a general domain and its general decomposition with a net-like structure. An almost optimal error estimate, without the factor H^−1/2 given in Dawson-Dupont’s error estimate, is proved. Another is the parallel domain decomposition algorithm of improved D-D type, in which an additional term is introduced to produce an approximation of an optimal error accuracy in L²-norm.     

Generalized Newton’s method based on graphical derivatives

This paper concerns developing a numerical method of the Newton type to solve systems of nonlinear equations described by nonsmooth continuous functions. We propose and justify a new generalized Newton algorithm based on graphical derivatives, which have never been used to derive a Newton-type method for solving nonsmooth equations. Based on advanced techniques of variational analysis and generalized differentiation, we establish the well-posedness of the algorithm, its local superlinear convergence, and its global convergence of the Kantorovich type. Our convergence results hold with no semismoothness and Lipschitzian assumptions, which is illustrated by examples. The algorithm and main results obtained in the paper are compared with well-recognized semismooth and B-differentiable versions of Newton’s method for nonsmooth Lipschitzian equations.     

B-spline collocation for solution of two-point boundary value problems

A numerical method based on B-spline is developed to solve the general nonlinear two-point boundary value problems up to order 6. The standard formulation of sextic spline for the solution of boundary value problems leads to non-optimal approximations. In order to derive higher orders of accuracy, high order perturbations of the problem are generated and applied to construct the numerical algorithm. The error analysis and convergence properties of the method are studied via Green’s function approach. O(h⁶) global error estimates are obtained for numerical solution of these classes of problems. Numerical results are given to illustrate the efficiency of the proposed method. Results of numerical experiments verify the theoretical behavior of the orders of convergence.     

On the fast reduction of a quasiseparable matrix to Hessenberg and tridiagonal forms

In this paper we design a fast new algorithm for reducing an N × N quasiseparable matrix to upper Hessenberg form via a sequence of N − 2 unitary transformations. The new reduction is especially useful when it is followed by the QR algorithm to obtain a complete set of eigenvalues of the original matrix. In particular, it is shown that in a number of cases some recently devised fast adaptations of the QR method for quasiseparable matrices can benefit from using the proposed reduction as a preprocessing step, yielding lower cost and a simplification of implementation.     

A global harmonic Arnoldi method for large non-Hermitian eigenproblems with an application to multiple eigenvalue problems

The global Arnoldi method can be used to compute exterior eigenpairs of a large non-Hermitian matrix A, but it does not work well for interior eigenvalue problems. Based on the global Arnoldi process that generates an F-orthonormal basis of a matrix Krylov subspace, we propose a global harmonic Arnoldi method for computing certain harmonic F-Ritz pairs that are used to approximate some interior eigenpairs. We propose computing the F-Rayleigh quotients of the large non-Hermitian matrix with respect to harmonic F-Ritz vectors and taking them as new approximate eigenvalues. They are better and more reliable than the harmonic F-Ritz values. The global harmonic Arnoldi method inherits convergence properties of the harmonic Arnoldi method applied to a larger matrix whose distinct eigenvalues are the same as those of the original given matrix. Some properties of the harmonic F-Ritz vectors are presented. As an application, assuming that A is diagonalizable, we show that the global harmonic Arnoldi method is able to solve multiple eigenvalue problems both in theory and in practice. To be practical, we develop an implicitly restarted global harmonic Arnoldi algorithm with certain harmonic F-shifts suggested. In particular, this algorithm can be adaptively used to solve multiple eigenvalue problems. Numerical experiments show that the algorithm is efficient for the eigenproblem and is reliable for quite ill-conditioned multiple eigenproblems.     

Optimal adaptive solution of initial-value problems with unknown singularities

The optimal solution of initial-value problems in ODEs is well studied for smooth right-hand side functions. Much less is known about the optimality of algorithms for singular problems. In this paper, we study the (worst case) solution of scalar problems with a right-hand side function having r   continuous bounded derivatives in R$\mathbf{R}$, except for an unknown singular point. We establish the minimal worst case error for such problems (which depends on r similarly as in the smooth case), and define optimal adaptive algorithms. The crucial point is locating an unknown singularity of the solution by properly adapting the grid. We also study lower bounds on the error of an algorithm for classes of singular problems. In the case of a single singularity with nonadaptive information, or in the case of two or more singularities, the error of any algorithm is shown to be independent of r.     

Two class of synchronous matrix multisplitting schemes for solving linear complementarity problems

Many problems in the areas of scientific computing and engineering applications can lead to the solution of the linear complementarity problem LCP (M,q). It is well known that the matrix multisplitting methods have been found very useful for solving LCP (M,q). In this article, by applying the generalized accelerated overrelaxation (GAOR) and the symmetric successive overrelaxation (SSOR) techniques, we introduce two class of synchronous matrix multisplitting methods to solve LCP (M,q). Convergence results for these two methods are presented when M is an H-matrix (and also an M-matrix). Also the monotone convergence of the new methods is established. Finally, the numerical results show that the introduced methods are effective for solving the large and sparse linear complementary problems.     

Convergence analysis of a conforming adaptive finite element method for an obstacle problem

The adaptive algorithm for the obstacle problem presented in this paper relies on the jump residual contributions of a standard explicit residual-based a posteriori error estimator. Each cycle of the adaptive loop consists of the steps ‘SOLVE’, ‘ESTIMATE’, ‘MARK’, and ‘REFINE’. The techniques from the unrestricted variational problem are modified for the convergence analysis to overcome the lack of Galerkin orthogonality. We establish R-linear convergence of the part of the energy above its minimal value, if there is appropriate control of the data oscillations. Surprisingly, the adaptive mesh-refinement algorithm is the same as in the unconstrained case of a linear PDE—in fact, there is no modification near the discrete free boundary necessary for R-linear convergence. The arguments are presented for a model obstacle problem with an affine obstacle χ and homogeneous Dirichlet boundary conditions. The proof of the discrete local efficiency is more involved than in the unconstrained case. Numerical results are given to illustrate the performance of the error estimator.     

Computation of the inverse Laplace transform based on a collocation method which uses only real values

We develop a numerical algorithm for inverting a Laplace transform (LT), based on Laguerre polynomial series expansion of the inverse function under the assumption that the LT is known on the real axis only. The method belongs to the class of Collocation methods (C-methods), and is applicable when the LT function is regular at infinity. Difficulties associated with these problems are due to their intrinsic ill-posedness. The main contribution of this paper is to provide computable estimates of truncation, discretization, conditioning and roundoff errors introduced by numerical computations. Moreover, we introduce the pseudoaccuracy which will be used by the numerical algorithm in order to provide uniform scaled accuracy of the computed approximation for any x   with respect to e^σx${\mathrm{e}}^{\sigma x}$. These estimates are then employed to dynamically truncate the series expansion. In other words, the number of the terms of the series acts like the regularization parameter which provides the trade-off between errors.     

Determination of a matrix function using the divided difference method of Newton and the interpolation technique of Hermite

Computing a function f(A) of an n-by-n matrix A is a frequently occurring problem in control theory and other applications. In this paper we introduce an effective approach for the determination of matrix function f(A). We propose a new technique which is based on the extension of Newton divided difference and the interpolation technique of Hermite and using the eigenvalues of the given matrix A. The new algorithm is tested on several problems to show the efficiency of the presented method. Finally, the application of this method in control theory is highlighted.     

Evaluating the boundary and covering degree of planar Minkowski sums and other geometrical convolutions

Algorithms are developed, based on topological principles, to evaluate the boundary and “internal structure” of the Minkowski sum of two planar curves. A graph isotopic to the envelope curve is constructed by computing its characteristic points. The edges of this graph are in one-to-one correspondence with a set of monotone envelope segments. A simple formula allows a degree   to be assigned to each face defined by the graph, indicating the number of times its points are covered by the Minkowski sum. The boundary can then be identified with the set of edges that separate faces of zero and non-zero degree, and the boundary segments corresponding to these edges can be approximated to any desired geometrical accuracy. For applications that require only the Minkowski sum boundary, the algorithm minimizes geometrical computations on the “internal” envelope edges, that do not contribute to the final boundary. In other applications, this internal structure is of interest, and the algorithm provides comprehensive information on the covering degree for different regions within the Minkowski sum. Extensions of the algorithm to the computation of Minkowski sums in R³${\mathbb{R}}^3$, and other forms of geometrical convolution, are briefly discussed.     

A generalized unitary Hessenberg matrix

This paper describes the use of a generalized isometric Arnoldi algorithm to reduce a unitary matrix, via unitary similarity, to a product of elementary reflectors and permutations. The computation is analogous to the reduction of a unitary matrix to a unitary Hessenberg matrix using the isometric Arnoldi algorithm. In the case in which A is a shift matrix, the reduction provides a novel recurrence for the factor R in the QR factorization of a Toeplitz-like matrix.