Numerical Solution Of Linear Systems Research Articles

Integral-algebraic equations with variable limits of integration. The case when the lower limit is at

The paper is devoted to the study and numerical solution of linear systems of integral equations with an identically singular matrix in the principal part (called integral-algebraic equations). The systems studied in the paper are fundamentally different from those previously studied in that the upper and lower limits of integration are variable. These systems of equations will be referred to as integral-algebraic equations with variable integration limits. The paper develops and justifies an algorithm for a numerical solution of the first order of accuracy and establishes conditions for the existence of a unique continuous solution to the given problem. There are numerical computations provided for model examples.

A nested divide-and-conquer method for tensor Sylvester equations with positive definite hierarchically semiseparable coefficients

Abstract Linear systems with a tensor product structure arise naturally when considering the discretization of Laplace-type differential equations or, more generally, multidimensional operators with separable coefficients. In this work, we focus on the numerical solution of linear systems of the form $$ \begin{align*} & \left(I\otimes \dots\otimes I \otimes A_1+\dots + A_d\otimes I \otimes\dots \otimes I\right)x=b, \end{align*}$$where the matrices $A_{t}\in \mathbb R^{n\times n}$ are symmetric positive definite and belong to the class of hierarchically semiseparable matrices. We propose and analyze a nested divide-and-conquer scheme, based on the technology of low-rank updates, which attains the quasi-optimal computational cost $\mathcal {O}(n^{d}\log (n))$. Our theoretical analysis highlights the role of inexactness in the nested calls of our algorithm and provides worst case estimates for the amplification of the residual norm. The performances are validated on 2D and 3D case studies.

The Reflection Method for the Numerical Solution of Linear Systems

The Reflection Method for the Numerical Solution of Linear Systems

A rational preconditioner for multi-dimensional Riesz fractional diffusion equations

We propose a rational preconditioner for an efficient numerical solution of linear systems arising from the discretization of multi-dimensional Riesz fractional diffusion equations. In particular, the discrete problem is obtained by employing finite difference or finite element methods to approximate the fractional derivatives of order α with α∈(1,2]. The proposed preconditioner is then defined as a rational approximation of the Riesz operator expressed as the integral of the standard heat diffusion semigroup. We show that, being the sum of k inverses of shifted Laplacian matrices, the resulting preconditioner belongs to the generalized locally Toeplitz class. As a consequence, we are able to provide the asymptotic description of the spectrum of the preconditioned matrices and we show that, despite the lack of clustering just as for the Laplacian, our preconditioner for α close to 1 and k≠1 reasonably small, provides better results than the Laplacian itself, while sharing the same computational complexity.

Algebraic Multigrid Preconditioner for Statically Condensed Systems Arising from Lowest-Order Hybrid Discretizations

We address the numerical solution of linear systems arising from the hybrid discretizations of second-order elliptic partial differential equations. Such discretizations hinge on a hybrid set of degrees of freedom (DoFs), respectively, defined in cells and faces, which naturally gives rise to a global hybrid system of linear equations. Assuming that the cell unknowns are only locally coupled, they can be efficiently eliminated from the system, leaving only face unknowns in the resulting Schur complement, which is also called the statically condensed matrix. We propose in this work an algebraic multigrid (AMG) preconditioner specifically targeting condensed systems corresponding to lowest-order discretizations (piecewise constant). Like traditional AMG methods, we retrieve geometric information on the coupling of the DoFs from algebraic data. However, as the condensed matrix only gives information on the faces, we use the uncondensed version to reconstruct the connectivity graph between elements and faces. An aggregation-based coarsening strategy mimicking a geometric coarsening or semicoarsening can then be set up to build coarse levels. Numerical experiments are performed on diffusion problems discretized by the hybrid high-order method at the lowest order. Our approach uses a K-cycle to precondition an outer flexible Krylov method. The results demonstrate similar performances, in most cases, compared to a standard AMG method and a notable improvement on anisotropic problems with Cartesian meshes.

Strong convergence analysis of iterative solvers for random operator equations

For the numerical solution of linear systems that arise from discretized linear partial differential equations, multigrid and domain decomposition methods are well established. Multigrid methods are known to have optimal complexity and domain decomposition methods are in particular useful for parallelization of the implemented algorithm. For linear random operator equations, the classical theory is not directly applicable, since condition numbers of system matrices may be close to degenerate due to non-uniform random input. It is shown that iterative methods converge in the strong, i.e. $$L^p$$, sense if the random input satisfies certain integrability conditions. As a main result, standard multigrid and domain decomposition methods are applicable in the case of linear elliptic partial differential equations with lognormal diffusion coefficients and converge strongly with deterministic bounds on the computational work which are essentially optimal. This enables the application of multilevel Monte Carlo methods with rigorous, deterministic bounds on the computational work.

A fast finite-difference algorithm for solving space-fractional filtration equation with a generalised Caputo derivative

The paper considers the problem of finite-difference simulation of filtration processes on the base of a space-fractional-differential equation with a generalised Caputo derivative with respect to a function. To make calculations faster, we propose an algorithm of fractional derivative approximation which reduces the computational complexity of matrix operations that arise while performing simulations on uniform and non-uniform grids. The algorithm is based on Taylor’s series expansions of fractional integral kernel and special representation of the matrices of linear systems that arise during finite-difference approximation. Theoretical estimates of the proposed algorithm’s complexity and accuracy are presented. We describe and experimentally compare computational schemes for matrix-vector multiplication used for numerical solution of linear systems by the BicgStab algorithm and analyse factors that influence its convergence. The results of computational experiments which demonstrate an efficiency of the proposed computational scheme for solving a filtration equation with a generalised Caputo derivative with respect to a space variable on a large-sized grids are given.

Application of the Jacobi–Davidson method for spectral low-rank preconditioning in computational electromagnetics problems

We consider the numerical solution of linear systems arising from computational electromagnetics applications. For large scale problems the solution is usually obtained iteratively with a Krylov subspace method. It is well known that for ill conditioned problems the convergence of these methods can be very slow or even it may be impossible to obtain a satisfactory solution. To improve the convergence a preconditioner can be used, but in some cases additional strategies are needed. In this work we study the application of spectral low-rank updates (SLRU) to a previously computed sparse approximate inverse preconditioner. The updates are based on the computation of a small subset of the eigenpairs closest to the origin. Thus, the performance of the SLRU technique depends on the method available to compute the eigenpairs of interest. The SLRU method was first used using the IRA’s method implemented in ARPACK. In this work we investigate the use of a Jacobi–Davidson method, in particular its JDQR variant. The results of the numerical experiments show that the application of the JDQR method to obtain the spectral low-rank updates can be quite competitive compared with the IRA’s method.

Numerical method for solving system of Fredhlom integral equations using Chebyshev cardinal functions

The focus of this paper is on the numerical solution of linear systems of Fredhlom integral equations of the second kind. For this purpose, the Chebyshev cardinal functions with Gauss-Lobatto points are used. By combination of properties of these functions and the effective Clenshaw-Curtis quadrature rule, an applicable numerical method for solving the mentioned systems is formulated. Some error bounds for the method are computed and its convergence rate is estimated. The method is numerically evaluated by solving some test problems caught from the literature by which the accuracy and computational efficiency of the method will be demonstrated.

Preconditioning of complex symmetric linear systems with applications in optical tomography

We consider the numerical solution of linear systems of the form (A+iκB)x=y, which arise in many applications, e.g., in time-harmonic acoustics, electromagnetics, or radiative transfer. We propose and analyze a class of preconditioners leading to complex symmetric iteration operators and investigate convergence of corresponding preconditioned iterative methods. Under mild assumptions on the operators A and B, we establish parameter and dimension independent convergence. The proposed methods are then applied to the solution of even-parity formulations of time-harmonic radiative transfer. For this application, we verify all assumptions required for our convergence analysis. The performance of the preconditioned iterations is then demonstrated by numerical tests supporting the theoretical results.

Adaptive Arnoldi-Tikhonov regularization for image restoration

In the framework of the numerical solution of linear systems arising from image restoration, in this paper we present an adaptive approach based on the reordering of the image approximations obtained with the Arnoldi-Tikhonov method. The reordering results in a modified regularization operator, so that the corresponding regularization can be interpreted as problem dependent. Numerical experiments are presented.

A projection method to solve linear systems in tensor format

SUMMARYIn this paper, we propose a method for the numerical solution of linear systems of equations in low rank tensor format. Such systems may arise from the discretisation of PDEs in high dimensions, but our method is not limited to this type of application. We present an iterative scheme, which is based on the projection of the residual to a low dimensional subspace. The subspace is spanned by vectors in low rank tensor format which—similarly to Krylov subspace methods—stem from the subsequent (approximate) application of the given matrix to the residual. All calculations are performed in hierarchical Tucker format, which allows for applications in high dimensions. The mode size dependency is treated by a multilevel method. We present numerical examples that include high‐dimensional convection–diffusion equations.Copyright © 2012 John Wiley & Sons, Ltd.

Numerical methods: principles, analyses, and algorithms

1. Background to Numerical Methods 2. Scope of Mathematical and Numerical Methods 3. Errors and their Propagation 4. Programming Tools and Techniques 5. Numerical Solution of One-Dimensional Non-Linear Equations 6. Numerical Solution of Linear Equations: Direct Methods 7. Numerical Solutions for Matrix Inversion 8. Numerical SOlution of Linear Systems of Equations: Iterative Methods 9. Eigenvalues and Eigenvalues 10. Numerical Solution of Systems of Nonlinear Equations 11. Linear Least Squares Problem 12. Interpolation and Extrapolation 13. Numerical Differentiation 14. Numeical Integration 15. Numerical Solution of Ordinary Differential Equations: Initial Value Problem 16. Numerical Solution of Ordinary Differential Equations: Boundary value Problem 17. Partial Differential Equations 18. Numerical Solution of Parabolic Partial Differential Equations 19. Numerical Solution of Hyperbolic Partial Differential Equations 20. Numerical Solution of Elliptic Partial Differential Equations 22. Advances in Numerical Methods Using Neurocomputing Paradigm 23. Numerical Solution of Differential Equations

On some Modifications of the Nekrassov Method for Numerical Solution of Linear Systems of Equations

A modification of the Nekrassov method for finding a solution of a linear system of algebraic equations is given and a numerical example is shown.

Krylov Subspace Methods for Linear Systems with Tensor Product Structure

The numerical solution of linear systems with certain tensor product structures is considered. Such structures arise, for example, from the finite element discretization of a linear PDE on a $d$-dimensional hypercube. Linear systems with tensor product structure can be regarded as linear matrix equations for $d=2$ and appear to be their most natural extension for $d>2$. A standard Krylov subspace method applied to such a linear system suffers from the curse of dimensionality and has a computational cost that grows exponentially with $d$. The key to breaking the curse is to note that the solution can often be very well approximated by a vector of low tensor rank. We propose and analyze a new class of methods, so-called tensor Krylov subspace methods, which exploit this fact and attain a computational cost that grows linearly with $d$.

Exponential Rosenbrock integrators for option pricing

In this paper, we are concerned with the time integration of differential equations modeling option pricing. In particular, we consider the Black–Scholes equation for American options. As an alternative to existing methods, we present exponential Rosenbrock integrators. These integrators require the evaluation of the exponential and related functions of the Jacobian matrix. The resulting methods have good stability properties. They are fully explicit and do not require the numerical solution of linear systems, in contrast to standard integrators. We have implemented some numerical experiments in Matlab showing the reliability of the new method.

Exponential Rosenbrock-Type Methods

We introduce a new class of exponential integrators for the numerical integration of large-scale systems of stiff differential equations. These so-called Rosenbrock-type methods linearize the flow in each time step and make use of the matrix exponential and related functions of the Jacobian. In contrast to standard integrators, the methods are fully explicit and do not require the numerical solution of linear systems. We analyze the convergence properties of these integrators in a semigroup framework of semilinear evolution equations in Banach spaces. In particular, we derive an abstract stability and convergence result for variable step sizes. This analysis further provides the required order conditions and thus allows us to construct pairs of embedded methods. We present a third-order method with two stages, and a fourth-order method with three stages, respectively. The application of the required matrix functions to vectors are computed by Krylov subspace approximations. We briefly discuss these implementation issues, and we give numerical examples that demonstrate the efficiency of the new integrators.

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.

Shadow block iteration for solving linear systems obtained from wavelet transforms

Matrices resulting from wavelet transforms have a special “shadow” block structure, that is, their small upper left blocks contain their lower frequency information. Numerical solutions of linear systems with such matrices require special care. We propose shadow block iterative methods for solving linear systems of this type. Convergence analysis for these algorithms are presented. We apply the algorithms to three applications: linear systems arising in the classical regularization with a single parameter for the signal de-blurring problem, multilevel regularization with multiple parameters for the same problem and the Galerkin method of solving differential equations. We also demonstrate the efficiency of these algorithms by numerical examples in these applications.

Semi-infinite multi-index perturbed block Toeplitz systems

In this article Banach algebra techniques are employed to study the numerical solution of linear systems with a semi-infinite multi-index suitably perturbed k× k block Toeplitz matrix. Decay properties of their solutions are studied by using suitably weighted Wiener algebras. A projection-type method for their numerical solution is introduced. Numerical results are presented illustrating both the accuracy of the method for certain perturbed Toeplitz systems and the need for generalizing the theoretical framework to other such systems.