ReLaTIve. An Ansi C90 software package for the Real Laplace Transform Inversion

A software package for numerical inversion of Laplace transforms computable everywhere on the real axis is described. Besides the function to invert the user has only to provide the numerical value (even if it is an approximate value) of the abscissa of convergence and the accuracy required for the inverse function. The software provides a controlled accuracy, i.e. it dynamically computes the so-called maximum attainable accuracy such that numerical results are provided within the greatest value between the user’s required accuracy and the maximum attainable accuracy. This is done because the intrinsic ill posedness of the real inversion problem sometime may prevent to reach the desired accuracy. The method implemented is based on a Laguerre polynomial series expansion of the inverse function and belongs to the class of polynomial-type methods of inversion of the Laplace transform, formally characterized as Collocation methods (C-methods).     

Inverse two-sided Laplace transform for probability density functions

We present a method for the numerical inversion of two-sided Laplace transform of a probability density function. The method assumes the knowledge of the first M derivatives at the origin of the function to be antitransformed. The approximate analytical form is obtained by resorting to maximum entropy principle. Both entropy and L₁-norm convergence are proved. Some numerical examples are illustrated.     

红黑排序混合算法收敛速度分析

The algorithm of applying the block Gauss elimination to the Red-Black or-dering matrix to reduce the order of the system then solve the reduced system byiterative methods is called Hybrid Red-Black Ordering algorithm.In this paper,we discuss the convergence rate of the hybrid methods combined with JACOBI,CG,GMRES(m).Theoretical analysis shows that without preconditioner thesethree hybrid methods converge about 2 times as fast as the corresponding natural ordering methods.For the case that all the eigenvalues is near the real axis, the GMRES(m) algorithm converges about 3 times faster than the natural ordering GMRES(m).Various numerical experiments are presented.For large scale prob-lem with preconditioners, numerical experiments show that the GMRES(m) hybrid methods converge from about 3 times to even 5 times as fast as the natural order-ing methods and the computing time is reduced to about 1/3 even 1/6 of that of the natural ordering methods.     

Error analysis of a Collocation method for numerically inverting a Laplace transform in case of real samples

In [S. Cuomo, L. D’Amore, A. Murli, M.R. Rizzardi, Computation of the inverse Laplace transform based on a collocation method which uses only real values, J. Comput. Appl. Math., 198 (1) (2007) 98–115] the authors proposed a Collocation method (C-method) for real inversion of Laplace transforms (Lt), based on the truncated Laguerre expansion of the inverse function:

where σ, b are parameters and c_k, kN, are the MacLaurin coefficients of a function depending on the Lt. The computational kernel of a C-method is the solution of a Vandermonde linear system, where the right hand side is obtained evaluating the Lt on the real axis. The Bjorck Pereira algorithm has been used for solving the Vandermonde linear system, providing a computable componentwise error bound on the solution.

For an inversion problem on discrete data F is known on a pre-assigned set of points (we refer to these points as samples of F) only and the major challenge is to deal with a significative loss of information. A natural approach to overcome this intrinsic difficulty is to construct a suitable fitting model that approximates the given data. In this case, we show that such approach leads to a C-method with perturbed right hand side, and then we use again the Bjorck Pereira algorithm.

Starting from the error introduced by the fitting model, we study its propagation in order to determine the maximum attainable accuracy on f_N. Moreover we derive a computable error bound that allows to get the suitable value of the parameter N that gives the maximum attainable accuracy.     

Approximate inverse preconditioner by computing approximate solution of Sylvester equation

This paper presents a new method for obtaining a matrix M which is an approximate inverse preconditioner for a given matrix A, where the eigenvalues of A all either have negative real parts or all have positive real parts. This method is based on the approximate solution of the special Sylvester equation AX + XA = 2I. We use a Krylov subspace method for obtaining an approximate solution of this Sylvester matrix equation which is based on the Arnoldi algorithm and on an integral formula. The computation of the preconditioner can be carried out in parallel and its implementation requires only the solution of very simple and small Sylvester equations. The sparsity of the preconditioner is preserved by using a proper dropping strategy. Some numerical experiments on test matrices from Harwell–Boing collection for comparing the numerical performance of the new method with an available well-known algorithm are presented.     

A preconditioned GMRES for complex dense linear systems from electromagnetic wave scattering problems

We consider the numerical solution of linear systems arising from the discretization of the electric field integral equation (EFIE). For some geometries the associated matrix can be poorly conditioned making the use of a preconditioner mandatory to obtain convergence. The electromagnetic scattering problem is here solved by means of a preconditioned GMRES in the context of the multilevel fast multipole method (MLFMM). The novelty of this work is the construction of an approximate hierarchically semiseparable (HSS) representation of the near-field matrix, the part of the matrix capturing interactions among nearby groups in the MLFMM, as preconditioner for the GMRES iterations. As experience shows, the efficiency of an ILU preconditioning for such systems essentially depends on a sufficient fill-in, which apparently sacrifices the sparsity of the near-field matrix. In the light of this experience we propose a multilevel near-field matrix and its corresponding HSS representation as a hierarchical preconditioner in order to substantially reduce the number of iterations in the solution of the resulting system of equations.     

GMRES, L-Curves, and Discrete Ill-Posed Problems

The GMRES method is a popular iterative method for the solution of large linear systems of equations with a nonsymmetric nonsingular matrix. This paper discusses application of the GMRES method to the solution of large linear systems of equations that arise from the discretization of linear ill-posed problems. These linear systems are severely ill-conditioned and are referred to as discrete ill-posed problems. We are concerned with the situation when the right-hand side vector is contaminated by measurement errors, and we discuss how a meaningful approximate solution of the discrete ill-posed problem can be determined by early termination of the iterations with the GMRES method. We propose a termination criterion based on the condition number of the projected matrices defined by the GMRES method. Under certain conditions on the linear system, the termination index corresponds to the vertex of an L-shaped curve.     

The Worst-Case GMRES for Normal Matrices

We study the convergence of GMRES for linear algebraic systems with normal matrices. In particular, we explore the standard bound based on a min-max approximation problem on the discrete set of the matrix eigenvalues. This bound is sharp, i.e. it is attainable by the GMRES residual norm. The question is how to evaluate or estimate the standard bound, and if it is possible to characterize the GMRES-related quantities for which this bound is attained (worst-case GMRES). In this paper we completely characterize the worst-case GMRES-related quantities in the next-to-last iteration step and evaluate the standard bound in terms of explicit polynomials involving the matrix eigenvalues. For a general iteration step, we develop a computable lower and upper bound on the standard bound. Our bounds allow us to study the worst-case GMRES residual norm as a function of the eigenvalue distribution. For hermitian matrices the lower bound is equal to the worst-case residual norm. In addition, numerical experiments show that the lower bound is generally very tight, and support our conjecture that it is to within a factor of 4/π of the actual worst-case residual norm. Since the worst-case residual norm in each step is to within a factor of the square root of the matrix size to what is considered an “average” residual norm, our results are of relevance beyond the worst case. This revised version was published online in July 2006 with corrections to the Cover Date.     

Weighted and deflated global GMRES algorithms for solving large Sylvester matrix equations

The solution of a large-scale Sylvester matrix equation plays an important role in control and large scientific computations. In this paper, we are interested in the large Sylvester matrix equation with large dimensionA and small dimension B, and a popular approach is to use the global Krylov subspace method. In this paper, we propose three new algorithms for this problem. We first consider the global GMRES algorithm with weighting strategy, which can be viewed as a precondition method. We present three new schemes to update the weighting matrix during iterations. Due to the growth of memory requirements and computational cost, it is necessary to restart the algorithm effectively. The deflation strategy is efficient for the solution of large linear systems and large eigenvalue problems; to the best of our knowledge, little work is done on applying deflation to the (weighted) global GMRES algorithm for large Sylvester matrix equations. We then consider how to combine the weighting strategy with deflated restarting, and propose a weighted global GMRES algorithm with deflation for solving large Sylvester matrix equations. In particular, we are interested in the global GMRES algorithm with deflation, which can be viewed as a special case when the weighted matrix is chosen as the identity. Theoretical analysis is given to show rationality of the new algorithms. Numerical experiments illustrate the numerical behavior of the proposed algorithms.

     

Review of inverse Laplace transform algorithms for Laplace-space numerical approaches

A boundary element method (BEM) simulation is used to compare the efficiency of numerical inverse Laplace transform strategies, considering general requirements of Laplace-space numerical approaches. The two-dimensional BEM solution is used to solve the Laplace-transformed diffusion equation, producing a time-domain solution after a numerical Laplace transform inversion. Motivated by the needs of numerical methods posed in Laplace-transformed space, we compare five inverse Laplace transform algorithms and discuss implementation techniques to minimize the number of Laplace-space function evaluations. We investigate the ability to calculate a sequence of time domain values using the fewest Laplace-space model evaluations. We find Fourier-series based inversion algorithms work for common time behaviors, are the most robust with respect to free parameters, and allow for straightforward image function evaluation re-use across at least a log cycle of time.     

Randomized Preprocessing of Homogeneous Linear Systems of Equations

Our randomized preprocessing enables pivoting-free and orthogonalization-free solution of homogeneous linear systems of equations. In the case of Toeplitz inputs, we decrease the estimated solution time from quadratic to nearly linear, and our tests show dramatic decrease of the CPU time as well. We prove numerical stability of our approach and extend it to solving nonsingular linear systems, inversion and generalized (Moore-Penrose) inversion of general and structured matrices by means of Newton’s iteration, approximation of a matrix by a nearby matrix that has a smaller rank or a smaller displacement rank, matrix eigen-solving, and root-finding for polynomial and secular equations and for polynomial systems of equations. Some by-products and extensions of our study can be of independent technical intersest, e.g., our extensions of the Sherman-Morrison-Woodbury formula for matrix inversion, our estimates for the condition number of randomized matrix products, and preprocessing via augmentation.     

An analysis of bilinear transform polynomial methods of inversion of Laplace transforms

Summary. Methods for the numerical inversion of a Laplace transform which use a special bilinear transformation of are particularly effective in many cases and are widely used. The main purpose of this paper is to analyze the convergence and conditioning properties of a special class of such methods, characterized by the use of Lagrange interpolation. The results derived apply both to complex and real inversion, and show that some known inversion methods are in fact in this class. Received June 21, 1993 / Revised version received March 10, 1994     

On the numerical solution of ill‐conditioned linear systems by regularization and iteration

We propose to reduce the (spectral) condition number of a given linear system by adding a suitable diagonal matrix to the system matrix, in particular by shifting its spectrum. Iterative procedures are then adopted to recover the solution of the original system. The case of real symmetric positive definite matrices is considered in particular, and several numerical examples are given. This approach has some close relations with Riley's method and with Tikhonov regularization. Moreover, we identify approximately the aforementioned procedure with a true action of preconditioning.     

Using successive approximations for improving the convergence of GMRES method

In this paper, our attention is concentrated on the GMRES method for the solution of the system (I–T)x=b of linear algebraic equations with a nonsymmetric matrix. We perform m pre-iterations y _l+1 =T _yl +b before starting GMRES and put y _m for the initial approximation in GMRES. We derive an upper estimate for the norm of the error vector in dependence on the mth powers of eigenvalues of the matrix T Further we study under what eigenvalues lay-out this upper estimate is the best one. The estimate shows and numerical experiments verify that it is advisable to perform pre-iterations before starting GMRES as they require fewer arithmetic operations than GMRES. Towards the end of the paper we present a numerical experiment for a system obtained by the finite difference approximation of convection-diffusion equations.     

Estimating the manipulation errors in finite element analysis, Part I

In Part I of this two-part work, the relative errors in representing and processing real numbers in digital computers are reviewed, and the computation of the digits lost from the relative errors is shown. The upper bound estimate of the relative error in the solution of Kξ = p for the displacements ξ is given as a function of the condition number of the stiffness matrix K, and the relative errors in K and the load vector p. The computation of relative errors for p and K by equilibrium checks is outlined, and various estimates for the condition number of the stiffness matrix are given.     

Linear system solution by null‐space approximation and projection (SNAP)

Solutions of large sparse linear systems of equations are usually obtained iteratively by constructing a smaller dimensional subspace such as a Krylov subspace. The convergence of these methods is sometimes hampered by the presence of small eigenvalues, in which case, some form of deflation can help improve convergence. The method presented in this paper enables the solution to be approximated by focusing the attention directly on the ‘small’ eigenspace (‘singular vector’ space). It is based on embedding the solution of the linear system within the eigenvalue problem (singular value problem) in order to facilitate the direct use of methods such as implicitly restarted Arnoldi or Jacobi–Davidson for the linear system solution. The proposed method, called ‘solution by null‐space approximation and projection’ (SNAP), differs from other similar approaches in that it converts the non‐homogeneous system into a homogeneous one by constructing an annihilator of the right‐hand side. The solution then lies in the null space of the resulting matrix. We examine the construction of a sequence of approximate null spaces using a Jacobi–Davidson style singular value decomposition method, called restarted SNAP‐JD, from which an approximate solution can be obtained. Relevant theory is discussed and the method is illustrated by numerical examples where SNAP is compared with both GMRES and GMRES‐IR. Copyright © 2006 John Wiley & Sons, Ltd.     

A problem on heat conduction in an insulated wire solved by numerical inverse Laplace transform

A problem of transient heat conduction in an insulated wire is solved by use of Laplace transform and numerical inversion. The problem is solved for the radiation boundary condition and also for the boundary condition of no heat flux through the outer surface of the insulation. The results are presented both numerically with four significant figures and graphically. Asymptotic expansions are derived for small and large values of the time variable. The numerical inversion of the Laplace transform is checked by comparison with the asymptotic expansions and with the numerical results obtained by a numerical inversion formula utilizing one more abscissa than the previous one.     

Asymptotic estimation for the condition numbers in BEM

Summary. We apply the boundary element methods (BEM) to the interior Dirichlet problem of the two dimensional Laplace equation, and its discretization is carried out with the collocation method using piecewise linear elements. In this paper, some precise asymptotic estimations for the discretization matrix (where denotes the division number) are investigated. We show that the maximum norm of and the condition number of have the forms: and , respectively, as , where the constants and are explicitly given in the proof. Although these estimates indicate illconditionedness of this numerical computation, the -convergence of this scheme with maximum norm is proved as an application of the main results. Received January 25, 1993 / Revised version received March 13, 1995     

Some Observations on the TLM Numerical Solution of the Laplace Equation

This paper describes progress on the TLM modelling of the Laplace equation, in particular, how the rate of convergence is influenced by the choice of scattering parameter for a particular discretisation. The hypothesis that optimum convergence is achieved when the real and imaginary parts for the lowest harmonic in a Fourier solution cancel appears to be upheld. The Fourier solution for the problem has been advanced by a better understanding of the nature of the initial excitation. The relationship between the form of the initial condition used in this and many other numerical solutions of the Laplace equation and oscillatory behavior in the results is given a firmer theoretical basis. A correlation between TLM numerical results and those obtained from matrix spectral radius calculations has confirmed much previous work.     

Preconditioning for radial basis functions with domain decomposition methods

In our previous work, an effective preconditioning scheme that is based upon constructing least-squares approximation cardinal basis functions (ACBFs) from linear combinations of the RBF-PDE matrix elements has shown very attractive numerical results. The preconditioner costs O(N²) flops to set up and O(N) storage. The preconditioning technique is sufficiently general that it can be applied to different types of different operators. This was applied to the 2D multiquadric method, with c~1/√N on the Poisson test problem, the preconditioned GMRES converges in tens of iterations. In this paper, we combine the RBF methods and the ACBF preconditioning technique with the domain decomposition method (DDM). We studied different implementations of the ACBF-DDM scheme and provide numerical results for N > 10,000 nodes. We shall demonstrate that the efficiency of the ACBF-DDM scheme improves dramatically as successively finer partitions of the domain are considered.