Showing papers in "Numerical Algorithms in 1999"
TL;DR: A short Matlab implementation for P1-x1 finite elements on triangles and parallelograms is provided for the numerical solution of elliptic problems with mixed boundary conditions on unstructured grids to prove the flexibility of the Matlab tool.
Abstract: A short Matlab implementation forP1-Q1 finite elements on triangles and parallelograms is provided for the numerical solution of elliptic problems with mixed boundary conditions on unstructured grids. According to the shortness of the program and the given documenta- tion, any adaptation from simple model examples to more complex problems can easily be performed. Numerical examples prove the flexibility of the Matlab tool. Unlike complex black-box commercial computer codes, this paper provides a simple and short open-box Matlab implementation of combined Courant's P1 triangles and Q1 elements on parallelograms for the numerical solutions of elliptic problems with mixed Dirichlet and Neumann boundary conditions. Based on four or five data files, arbitrary regular triangulations are determined. Instead of covering all kinds of possible problems in one code, the proposed tool aims to be simple, easy to understand and modify. Therefore, only simple model examples are included to be adapted to whatever is needed. In further contributions we provide more complicated elements, a posteriori error estimators and flexible adaptive mesh-refining algorithms. Compared to the latest Matlab toolbox (3), our approach is shorter, allows more elements, is easily adapted to modified problems like convection terms, and is open to easy modifications for basically any type of finite element. The rest of the paper is organised as follows. As a model problem, the Laplace
223 citations
TL;DR: This work investigates the numerical solution of the stable generalized Lyapunov equation via the sign function method and considers some modifications and discusses how to solve generalized LyAPunov equations with semidefinite constant term for the Cholesky factor.
Abstract: We investigate the numerical solution of the stable generalized Lyapunov equation via the sign function method This approach has already been proposed to solve standard Lyapunov equations in several publications The extension to the generalized case is straightforward We consider some modifications and discuss how to solve generalized Lyapunov equations with semidefinite constant term for the Cholesky factor The basic computational tools of the method are basic linear algebra operations that can be implemented efficiently on modern computer architectures and in particular on parallel computers Hence, a considerable speed-up as compared to the Bartels–Stewart and Hammarling methods is to be expected We compare the algorithms by performing a variety of numerical tests
154 citations
TL;DR: It is shown that the three dimensional alignment problem of a histological data set of the human brain may be phrased in terms of a nonlinear PDE, and a fast direct solution technique is devised for the associated structured large systems of linear equations.
Abstract: In recent years, new nonlinear partial differential equation (PDE) based approaches have become popular for solving image processing problems. Although the outcome of these methods is often very promising, their actual realization is in general computationally intensive. Therefore, accurate and efficient schemes are needed. The aim of this paper is twofold. First, we will show that the three dimensional alignment problem of a histological data set of the human brain may be phrased in terms of a nonlinear PDE. Second, we will devise a fast direct solution technique for the associated structured large systems of linear equations. In addition, the potential of the derived method is demonstrated on real-life data.
96 citations
TL;DR: An iterative boundary element method based on the relaxed algorithm introduced in [8] is used to solve numerically a class of inverse boundary problems and improves the rate of convergence of Kozlov's scheme.
Abstract: In this paper, an iterative boundary element method based on our relaxed algorithm introduced in [8] is used to solve numerically a class of inverse boundary problems A dynamical choice of the relaxation parameter is presented and a stopping criterion based on our theoretical results is used The numerical results show that the algorithm produces a reasonably approximate solution and improves the rate of convergence of Kozlov's scheme [10]
58 citations
TL;DR: This communication describes Version 3.0 of Regularization Tools, a Matlab package for analysis and solution of discrete ill-posed problems.
Abstract: This communication describes Version 3.0 of Regularization Tools, a Matlab package for analysis and solution of discrete ill-posed problems.
54 citations
TL;DR: This paper constructs and analyzes two compact monotone finite difference methods to solve singularly perturbed problems of convection–diffusion type and proves that the methods are uniformly convergent of order two and three except for a logarithmic factor.
Abstract: In this paper we construct and analyze two compact monotone finite difference methods to solve singularly perturbed problems of convection–diffusion type. They are defined as HODIE methods of order two and three, i.e., the coefficients are determined by imposing that the local error be null on a polynomial space. For arbitrary meshes, these methods are not adequate for singularly perturbed problems, but using a Shishkin mesh we can prove that the methods are uniformly convergent of order two and three except for a logarithmic factor. Numerical examples support the theoretical results.
49 citations
TL;DR: A new method similar to QMR but based on the Hessenberg process instead of the Lanczos process is given, which is less expensive and requires slightly less storage than GMRES.
Abstract: The Generalized Minimal Residual (GMRES) method and the Quasi-Minimal Residual (QMR) method are two Krylov methods for solving linear systems. The main difference between these methods is the generation of the basis vectors for the Krylov subspace. The GMRES method uses the Arnoldi process while QMR uses the Lanczos algorithm for constructing a basis of the Krylov subspace. In this paper we give a new method similar to QMR but based on the Hessenberg process instead of the Lanczos process. We call the new method the CMRH method. The CMRH method is less expensive and requires slightly less storage than GMRES. Numerical experiments suggest that it has behaviour similar to GMRES.
47 citations
TL;DR: In this paper, the authors describe a Matlab 5.2 package for computing and modifying certain rank-revealing decompositions that have found widespread use in signal processing and other applications.
Abstract: We describe a Matlab 5.2 package for computing and modifying certain rank-revealing decompositions that have found widespread use in signal processing and other applications. The package focuses on algorithms for URV and ULV decompositions, collectively known as UTV decompositions. We include algorithms for the ULLV decomposition, which generalizes the ULV decomposition to a pair of matrices. For completeness a few algorithms for computation of the RRQR decomposition are also included. The software in this package can be used as is, or can be considered as templates for specialized implementations on signal processors and similar dedicated hardware platforms.
45 citations
TL;DR: This paper presents in a survey form the basic results of preconditioners of algebraic multilevel type and considers in particular additive methods, which has excellent parallelization properties.
Abstract: There exist two main versions of preconditioners of algebraic multilevel type, the additive and the multiplicative methods. They correspond to preconditioners in block diagonal and block matrix factorized form, respectively. Both can be defined and analysed as recursive two-by-two block methods. Although the analytical framework for such methods is simple, for many finite element approximations it still permits the derivation of the strongest results, such as optimal, or nearly optimal, rate of convergence and optimal, or nearly optimal order of computational complexity, when proper recursive global orderings of node points have been used or when they are applied for hierarchical basis function finite element methods for elliptic self-adjoint equations and stabilized in a certain way. This holds for general elliptic problems of second order, independent of the regularity of the problem, including independence of discontinuities of coefficients between elements and of anisotropy. Important ingredients in the methods are a proper balance of the size of the coarse mesh to the finest mesh and a proper solver on the coarse mesh. This paper presents in a survey form the basic results of such methods and considers in particular additive methods. This method has excellent parallelization properties.
38 citations
TL;DR: This paper considers algorithms to compute bounds of the A-norm of the error in the preconditioned conjugate gradient (PCG) algorithm and gives numerical experiments which show that good upper and lower bounds can be obtained provided estimates of the lowest and largest eigenvalues of the precONDitioned matrix are given or adaptively computed.
Abstract: In this paper we consider algorithms to compute bounds of the A-norm of the error in the preconditioned conjugate gradient (PCG) algorithm. We extend to PCG formulas that were given in an earlier paper [8]. We give numerical experiments which show that good upper and lower bounds can be obtained provided estimates of the lowest and largest eigenvalues of the preconditioned matrix are given or adaptively computed.
38 citations
TL;DR: In this article, backward and forward error analysis of corner cutting algorithms are performed and running error is also analyzed and as a consequence the general algorithm is modified to include the computation of an error bound.
Abstract: Corner cutting algorithms are used in different fields and, in particular, play a relevant role in Computer Aided Geometric Design Evaluation algorithms such as the de Casteljau algorithm for polynomials and the de Boor–Cox algorithm for B‐splines are examples of corner cutting algorithms Here backward and forward error analysis of corner cutting algorithms are performed The running error is also analyzed and as a consequence the general algorithm is modified to include the computation of an error bound
TL;DR: The concept of effective order allows for the possibility that the result computed in a Runge-Kutta step is an approximation to some quantity more general than the actual solution at a step point, and is applied here to singly-implicit methods.
Abstract: The concept of effective order allows for the possibility that the result computed in a Runge-Kutta step is an approximation to some quantity more general than the actual solution at a step point. This generalization is applied here to singly-implicit methods. The limitation that requires severe and inconvenient restrictions on the abscissae in the method is removed under this widening of the order requirement and all that is now needed is that the abscissae be distinct. Implementation questions, such as error estimation, stepsize change and dense output are also considered.
TL;DR: Numerical results demonstrate that the error estimation employed in the code is very reliable and that the step and order changing strategies are very robust, and this code outperforms the Matlab ode45 code for moderate and stringent tolerances.
Abstract: The issues related to the development of a new code for nonstiff ordinary differential equations are discussed. This code is based on the Nordsieck representation of type 1 DIMSIMs, implemented in a variable-step size variable-order mode. Numerical results demonstrate that the error estimation employed in the code is very reliable and that the step and order changing strategies are very robust. This code outperforms the Matlab ode45 code for moderate and stringent tolerances.
TL;DR: A deflation method is introduced that takes advantage of the IRA method, by extracting a GMRES solution from the Krylov basis computed within the Arnoldi process ofThe IRA method itself.
Abstract: We introduce a deflation method that takes advantage of the IRA method, by extracting a GMRES solution from the Krylov basis computed within the Arnoldi process of the IRA method itself. The deflation is well-suited because it is done with eigenvectors associated to the eigenvalues that are closest to zero, which are approximated by IRA very quickly. By a slight modification, we adapt it to the FOM algorithm, and then to GMRES enhanced by imposing constraints within the minimization condition. The use of IRA enables us to reduce the number of matrix-vector products, while keeping a low storage.
TL;DR: Algorithms in which alternating ways of re-using a given set of stored difference vectors are outlined in which the initial Hessian approximation is defined implicitly like the L-BFGS Hessian in terms of some stored vectors rather than the usual choice of a multiple of the unit matrix.
Abstract: This paper considers simple modifications of the limited memory BFGS (L-BFGS) method for large scale optimization It describes algorithms in which alternating ways of re-using a given set of stored difference vectors are outlined The proposed algorithms resemble the L-BFGS method, except that the initial Hessian approximation is defined implicitly like the L-BFGS Hessian in terms of some stored vectors rather than the usual choice of a multiple of the unit matrix Numerical experiments show that the new algorithms yield desirable improvement over the L-BFGS method
TL;DR: Sample numerical procedures for constructing one-coordinates, two-coordinate and arbitrary finite-coordinated generic differential quadrature models are presented.
Abstract: The differential quadrature (DQ) is generalized. Various methods for generating the weighting coefficients are developed. The design of a grid model is flexible. Weighting coefficients for general multi-coordinate grid models with arbitrary configurations can also be calculated. The calculation of weighting coefficients is easy. Sample numerical procedures for constructing one-coordinate, two-coordinate and arbitrary finite-coordinate generic differential quadrature models are presented.
TL;DR: The zeros of the hypergeometric polynomials lie asymptotically as n > \infty as well as Auxiliary results for the asymPTotic zero distribution of other functions related to hypergeometry are proved, including Jacobi polynmials with varying parameters and associated Legendre functions.
Abstract: We show that the zeros of the hypergeometric polynomials F ( n,kn + 1;kn + 2;z), k,n2N, cluster on the loop of the lemniscate {z: jz k (1 z)j = k k =(k + 1) k+1 ,R e(z) > k=(k + 1)} as n!1. We also state the equations of the curves on which the zeros of F ( n,kn+ 1; (k +‘)n + 2;z), k,‘,n2N, lie asymptotically as n!1. Auxiliary results for the asymptotic zero distribution of other functions related to hypergeometric polynomials are proved, including Jacobi polynomials with varying parameters and associated Legendre functions. Graphical evidence is provided using Mathematica.
TL;DR: It is shown that, with homogeneous Dirichlet boundary conditions, the condition number of finite element discretization matrices remains uniformly bounded independent of the size of the boundary elements provided that the sizes of the elements increases with their distance to the boundary.
Abstract: It is shown that, with homogeneous Dirichlet boundary conditions, the condition number of finite element discretization matrices remains uniformly bounded independent of the size of the boundary elements provided that the size of the elements increases with their distance to the boundary. This fact allows the construction of simple multigrid methods of optimal complexity for domains of nearly arbitrary shape.
TL;DR: The CGNR-method is considered, that is, the classical method of conjugate gradients by Hestenes and Stiefel applied to the associated normal equations of Ta = f* in Hilbert spaces, and two a posteriori stopping rules are introduced.
Abstract: We consider an ill-posed problem Ta = f* in Hilbert spaces and suppose that the linear bounded operator T is approximately available, with a known estimate for the operator perturbation at the solution. As a numerical scheme the CGNR-method is considered, that is, the classical method of conjugate gradients by Hestenes and Stiefel applied to the associated normal equations. Two a posteriori stopping rules are introduced, and convergence results are provided for the corresponding approximations, respectively. As a specific application, a parameter estimation problem is considered.
TL;DR: This paper determines the expressions of Bernstein type bases in spaces spanned by the constants and power functions with consecutive integer exponents in extended Chebyshev spaces.
Abstract: Extended Chebyshev spaces possess Bernstein type bases. In this paper, we determine the expressions of such bases in spaces spanned by the constants and power functions with consecutive integer exponents.
TL;DR: A biharmonic-type interpolation method is presented to solve 2D and 3D scattered data interpolation problems that makes it possible to significantly reduce the computational cost of the evaluation of the appearing domain integrals as well as the memory requirement of the procedure.
Abstract: A biharmonic-type interpolation method is presented to solve 2D and 3D scattered data interpolation problems. Unlike the methods based on radial basis functions, which produce a large linear system of equations with fully populated and often non-selfadjoint and ill-conditioned matrix, the presented method converts the interpolation problem to the solution of the biharmonic equation supplied with some non-usual boundary conditions at the interpolation points. To solve the biharmonic equation, fast multigrid techniques can be applied which are based on a non-uniform, non-equidistant but Cartesian grid generated by the quadtree/octtree algorithm. The biharmonic interpolation technique is applied to the multiple and dual reciprocity method of the BEM to convert domain integrals to the boundary. This makes it possible to significantly reduce the computational cost of the evaluation of the appearing domain integrals as well as the memory requirement of the procedure. The resulting method can be considered as a special grid-free technique, since it requires no domain discretisation.
TL;DR: A significant refinement of normal forms for differentiable maps near a fixed point is proposed, rational in the sense that if the coefficients of a map are in a field K, so is its normal form.
Abstract: We propose in this paper a significant refinement of normal forms for differentiable maps near a fixed point. We give a method to obtain further reduction of classical normal forms with concrete and interesting applications. Our method leads to unique normal forms either with respect to general diffeomorphisms in certain cases or with respect to near identity diffeomorphisms in other cases. Our approach is rational in the sense that if the coefficients of a map are in a field K, so is its normal form.
TL;DR: The behaviour of these iterative methods is investigated to define the effective strategy for this class of problems and to avoid breakdown in the numerical algorithm introduced by the use of Slotboom variables.
Abstract: The drift-diffusion model can be described by a nonlinear Poisson equation for the electrostatic potential coupled with a system of convection-reaction-diffusion equations for the transport of charge. We use a Gummel-like process [10] to decouple this system. Each of the obtained equations is discretised with the finite element method. We use a local scaling method to avoid breakdown in the numerical algorithm introduced by the use of Slotboom variables. Solution of the discrete nonlinear Poisson equation is accomplished with quasi-Newton methods. The nonsymmetric discrete transport equations are solved using an incomplete LU factorization preconditioner in conjunction with some robust linear solvers, such as (CGS), (BI-CGSTAB) and (GMRES). We investigate the behaviour of these iterative methods to define the effective strategy for this class of problems.
TL;DR: This paper describes a software package which implements a deflated and variable-block version of the Davidson method, a preconditioned eigenvalue technique aimed at computing a few of the extreme eigenpairs of large sparse symmetric matrices.
Abstract: The Davidson method is a preconditioned eigenvalue technique aimed at computing a few of the extreme (i.e., leftmost or rightmost) eigenpairs of large sparse symmetric matrices. This paper describes a software package which implements a deflated and variable-block version of the Davidson method. Information on how to use the software is provided. Guidelines for its upgrading or for its incorporation into existing packages are also included. Various experiments are performed on an SGI Power Challenge and comparisons with ARPACK are reported.
TL;DR: A package of computer programs for the unified treatment of initial-value problems for systems of ordinary differential equations which implement a numerical method which is efficient for a general class of differential equations.
Abstract: This paper describes a package of computer programs for the unified treatment of initial-value problems for systems of ordinary differential equations. The programs implement a numerical method which is efficient for a general class of differential equations. The user may determine the solutions over finite or infinite intervals. The solutions may have singularities at the end-points of the interval for which the solution is sought. Besides giving the initial values and the analytical expression for the differential equations to be solved the user needs to specify the nature of the singularities and give some other analytical information as described in the paper in order to take advantage of the speed and accuracy of the package described here.
TL;DR: This paper studies the Generalized Minimal Residual (GMRES) method for solving singular linear systems, particularly when the necessary and sufficient condition to obtain a Krylov solution is not satisfied and analyzes the convergence of GMRES and restarted GMRES.
Abstract: In this paper, we study the Generalized Minimal Residual (GMRES) method for solving singular linear systems, particularly when the necessary and sufficient condition to obtain a Krylov solution is not satisfied. Thanks to some new results which may be applied in exact arithmetic or in finite precision, we analyze the convergence of GMRES and restarted GMRES. These formulas can also be used in the case when the systems are nonsingular. In particular, it allows us to understand what is often referred to as stagnation of the residual norm of GMRES.
TL;DR: This paper deduces the recurrence relations for the Chebyshev and the Modified Chebyshv methods for every d∈ℕ and shows several formal and numerical tests realized with the software developed.
Abstract: This paper is concerned with the Shohat-Favard, Chebyshev and Modified Chebyshev methods for d-orthogonal polynomial sequences d∈ℕ. Shohat-Favard’s method is presented from the concept of dual sequence of a sequence of polynomials. We deduce the recurrence relations for the Chebyshev and the Modified Chebyshev methods for every d∈ℕ. The three methods are implemented in the Mathematica programming language. We show several formal and numerical tests realized with the software developed.
TL;DR: It is proved the existence of an asymptotic expansion of the error, which justifies the use of Richardson extrapolation, and how these expansions can be translated to a new expansion of potentials calculated with the numerical solution of a boundary integral equation such as those treated before.
Abstract: In this paper we propose a fully discretized version of the collocation method applied to integral equations of the first kind with logarithmic kernel. After a stability and convergence analysis is given, we prove the existence of an asymptotic expansion of the error, which justifies the use of Richardson extrapolation. We further show how these expansions can be translated to a new expansion of potentials calculated with the numerical solution of a boundary integral equation such as those treated before. Some numerical experiments, confirming our theoretical results, are given.
TL;DR: Finite difference approximations of Euler equations, the related solution algorithms, and applications to segmentation problems by using synthetic images are defined and discussed.
Abstract: This paper deals with finite-difference approximations of Euler equations arising in the variational formulation of image segmentation problems. We illustrate how they can be defined by the following steps: (a) definition of the minimization problem for the Mumford–Shah functional (MSf), (b) definition of a sequence of functionals Γ-convergent to the MSf, and (c) definition and numerical solution of the Euler equations associated to the kth functional of the sequence. We define finite difference approximations of the Euler equations, the related solution algorithms, and we present applications to segmentation problems by using synthetic images. We discuss application results, and we mainly analyze computed discontinuity contours and convergence histories of method executions.
TL;DR: A new iterative method, LAN/MGMRES, is introduced, which is designed to combine the reliability of GMRES with the reduced work of a Lanczos-type method and is comparable in terms of both the approximate number of iterations and the overall convergence behavior.
Abstract: For solving nonsymmetric linear systems, the well-known GMRES method is considered to be a stable method; however, the work per iteration increases as the number of iterations increases. We consider two new iterative methods GGMRES and MGMRES, which are a generalization and a modification of the GMRES method, respectively. Instead of using a minimization condition as in the derivation of GGMRES, we use a Galerkin condition to derive the MGMRES method. We also introduce another new iterative method, LAN/MGMRES, which is designed to combine the reliability of GMRES with the reduced work of a Lanczos-type method. A computer program has been written based on the use of the LAN/MGMRES algorithm for solving nonsymmetric linear systems arising from certain elliptic problems. Numerical tests are presented comparing this algorithm with some other commonly used iterative algorithms. These preliminary tests of the LAN/MGMRES algorithm show that it is comparable in terms of both the approximate number of iterations and the overall convergence behavior.