Computational Krylov‐based methods for large‐scale differential Sylvester matrix problems

In the present paper, we propose Krylov‐based methods for solving large‐scale differential Sylvester matrix equations having a low‐rank constant term. We present two new approaches for solving such differential matrix equations. The first approach is based on the integral expression of the exact solution and a Krylov method for the computation of the exponential of a matrix times a block of vectors. In the second approach, we first project the initial problem onto a block (or extended block) Krylov subspace and get a low‐dimensional differential Sylvester matrix equation. The latter problem is then solved by some integration numerical methods such as the backward differentiation formula or Rosenbrock method, and the obtained solution is used to build the low‐rank approximate solution of the original problem. We give some new theoretical results such as a simple expression of the residual norm and upper bounds for the norm of the error. Some numerical experiments are given in order to compare the two approaches.     

A computational method for large‐scale differential symmetric Stein equation

We propose a numerical method for solving large‐scale differential symmetric Stein equations having low‐rank right constant term. Our approach is based on projection the given problem onto a Krylov subspace then solving the low dimensional matrix problem by using an integration method, and the original problem solution is built by using obtained low‐rank approximate solution. Using the extended block Arnoldi process and backward differentiation formula (BDF), we give statements of the approximate solution and corresponding residual. Some numerical results are given to show the efficiency of the proposed method.     

A structure‐preserving doubling algorithm for Lur'e equations

We introduce a numerical method for the numerical solution of the Lur'e equations, a system of matrix equations that arises, for instance, in linear‐quadratic infinite time horizon optimal control. We focus on small‐scale, dense problems. Via a Cayley transformation, the problem is transformed to the discrete‐time case, and the structural infinite eigenvalues of the associated matrix pencil are deflated. The deflated problem is associated with a symplectic pencil with several Jordan blocks of eigenvalue 1 and even size, which arise from the nontrivial Kronecker chains at infinity of the original problem. For the solution of this modified problem, we use the structure‐preserving doubling algorithm. Implementation issues such as the choice of the parameter γ in the Cayley transform are discussed. The most interesting feature of this method, with respect to the competing approaches, is the absence of arbitrary rank decisions, which may be ill‐posed and numerically troublesome. The numerical examples presented confirm the effectiveness of this method. Copyright © 2015 John Wiley & Sons, Ltd.     

On A‐stability for some two‐step W‐methods

In this paper we prove the A‐stability for some sequential and parallel two‐step W‐methods (TSW‐methods and PTSW‐methods) which have been constructed recently for the numerical solution of large stiff ODE systems. To show that the spectral radius of the amplification matrix of a fixed method is less than one we apply Schur's criterion recursively to the characteristic polynomial of the amplification matrix. The computations can be done automatically using computer algebra tools, e.g. Maple 7.     

A numerical approach for solving weakly singular partial integro‐differential equations via two‐dimensional‐orthonormal Bernstein polynomials with the convergence analysis

In this paper, we develop an efficient matrix method based on two‐dimensional orthonormal Bernstein polynomials (2D‐OBPs) to provide approximate solution of linear and nonlinear weakly singular partial integro‐differential equations (PIDEs). First, we approximate all functions involved in the considerable problem via 2D‐OBPs. Then, by using the operational matrices of integration, differentiation, and product, the solution of Volterra singular PIDEs is transformed to the solution of a linear or nonlinear system of algebraic equations which can be solved via some suitable numerical methods. With a small number of bases, we can find a reasonable approximate solution. Moreover, we establish some useful theorems for discussing convergence analysis and obtaining an error estimate associated with the proposed method. Finally, we solve some illustrative examples by employing the presented method to show the validity, efficiency, high accuracy, and applicability of the proposed technique.     

Hierarchical Schur complement preconditioner for the stochastic Galerkin finite element methods

Use of the stochastic Galerkin finite element methods leads to large systems of linear equations obtained by the discretization of tensor product solution spaces along their spatial and stochastic dimensions. These systems are typically solved iteratively by a Krylov subspace method. We propose a preconditioner, which takes an advantage of the recursive hierarchy in the structure of the global matrices. In particular, the matrices posses a recursive hierarchical two‐by‐two structure, with one of the submatrices block diagonal. Each of the diagonal blocks in this submatrix is closely related to the deterministic mean‐value problem, and the action of its inverse is in the implementation approximated by inner loops of Krylov iterations. Thus, our hierarchical Schur complement preconditioner combines, on each level in the approximation of the hierarchical structure of the global matrix, the idea of Schur complement with loops for a number of mutually independent inner Krylov iterations, and several matrix–vector multiplications for the off‐diagonal blocks. Neither the global matrix nor the matrix of the preconditioner need to be formed explicitly. The ingredients include only the number of stiffness matrices from the truncated Karhunen–Loève expansion and a good preconditioned for the mean‐value deterministic problem. We provide a condition number bound for a model elliptic problem, and the performance of the method is illustrated by numerical experiments. Copyright © 2013 John Wiley & Sons, Ltd.     

On greedy randomized coordinate descent methods for solving large linear least‐squares problems

For solving large scale linear least‐squares problem by iteration methods, we introduce an effective probability criterion for selecting the working columns from the coefficient matrix and construct a greedy randomized coordinate descent method. It is proved that this method converges to the unique solution of the linear least‐squares problem when its coefficient matrix is of full rank, with the number of rows being no less than the number of columns. Numerical results show that the greedy randomized coordinate descent method is more efficient than the randomized coordinate descent method.     

A matrix-splitting method for symmetric affine second-order cone complementarity problems

The affine second-order cone complementarity problem (SOCCP) is a wide class of problems that contains the linear complementarity problem (LCP) as a special case. The purpose of this paper is to propose an iterative method for the symmetric affine SOCCP that is based on the idea of matrix splitting. Matrix-splitting methods have originally been developed for the solution of the system of linear equations and have subsequently been extended to the LCP and the affine variational inequality problem. In this paper, we first give conditions under which the matrix-splitting method converges to a solution of the affine SOCCP. We then present, as a particular realization of the matrix-splitting method, the block successive overrelaxation (SOR) method for the affine SOCCP involving a positive definite matrix, and propose an efficient method for solving subproblems. Finally, we report some numerical results with the proposed algorithm, where promising results are obtained especially for problems with sparse matrices.     

Uzawa methods for a class of block three‐by‐three saddle‐point problems

In this work, we consider numerical methods for solving a class of block three‐by‐three saddle‐point problems, which arise from finite element methods for solving time‐dependent Maxwell equations and some other applications. The direct extension of the Uzawa method for solving this block three‐by‐three saddle‐point problem requires the exact solution of a symmetric indefinite system of linear equations at each step. To avoid heavy computations at each step, we propose an inexact Uzawa method, which solves the symmetric indefinite linear system in some inexact way. Under suitable assumptions, we show that the inexact Uzawa method converges to the unique solution of the saddle‐point problem within the approximation level. Two special algorithms are customized for the inexact Uzawa method combining the splitting iteration method and a preconditioning technique, respectively. Numerical experiments are presented, which demonstrated the usefulness of the inexact Uzawa method and the two customized algorithms.     

A new Jacobi operational matrix: An application for solving fractional differential equations

In this paper, we derived the shifted Jacobi operational matrix (JOM) of fractional derivatives which is applied together with spectral tau method for numerical solution of general linear multi-term fractional differential equations (FDEs). A new approach implementing shifted Jacobi operational matrix in combination with the shifted Jacobi collocation technique is introduced for the numerical solution of nonlinear multi-term FDEs. The main characteristic behind this approach is that it reduces such problems to those of solving a system of algebraic equations which greatly simplifying the problem. The proposed methods are applied for solving linear and nonlinear multi-term FDEs subject to initial or boundary conditions, and the exact solutions are obtained for some tested problems. Special attention is given to the comparison of the numerical results obtained by the new algorithm with those found by other known methods.     

Generating equations approach for quadratic matrix equations

We show how Van Loan's method for annulling the (2,1) block of skew‐Hamiltonian matrices by symplectic‐orthogonal similarity transformation generalizes to general matrices and provides a numerical algorithm for solving the general quadratic matrix equation: For skew‐Hamiltonian matrices we find their canonical form under a similarity transformation and find the class of all symplectic‐orthogonal similarity transformations for annulling the (2,1) block and simultaneously bringing the (1,1) block to Hessenberg form. We present a structure‐preserving algorithm for the solution of continuous‐time algebraic Riccati equation. Unlike other methods in the literature, the final transformed Hamiltonian matrix is not in Hamiltonian–Schur form. Three applications are presented: (a) for a special system of partial differential equations of second order for a single unknown function, we obtain the matrix of partial derivatives of second order of the unknown function by only algebraic operations and differentiation of functions; (b) for a similar transformation of a complex matrix into a symmetric (and three‐diagonal) one by applying only finite algebraic transformations; and (c) for finite‐step reduction of the eigenvalues–eigenvectors problem of a Hermitian matrix to the eigenvalues– eigenvectors problem of a real symmetric matrix of the same dimension. Copyright © 1999 John Wiley & Sons, Ltd.     

Spectral smoothed boundary methods: The role of external boundary conditions

The spectral smoothed boundary method (SSBM) is a recently proposed numerical method to approximate the solution of partial differential equations in irregular domains with no‐flux boundary conditions by means of Fourier spectral methods. In this article we explore the robustness and accuracy of the scheme under variations of the artificial boundary conditions that must be imposed on the boundary of the enlarged domain in which the problem is solved. As a test model, we present quantitative numerical results based on a problem of propagation of waves of electrical activity in cardiac tissue for which the method is relevant. © 2005 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2006     

Error estimations for some meshless boundary interpolation methods

Meshless methods have become quite popular in numerical treatment of partial differential equations because of their simplicity and the fact that they require neither domain nor boundary mesh. In general, however, they convert the original problem to a highly ill‐conditioned linear system of algebraic equations with a dense matrix. Recently, a special technique has been proposed which circumvents this computational difficulty. This method, called Direct Multi‐Elliptic Interpolation Method, is based on a scattered data interpolation which defines the interpolation function as a solution of a higher order multi‐elliptic equation. Here the boundary version of this meshless method which is based on a multi‐elliptic boundary interpolation is considered. Error estimations are derived justifying the interpolation function to be a good approximation of the solution of the original boundary value problem as well. At the same time, the problem of large, dense and ill‐conditioned matrices as well as the mesh generation are completely avoided. (© 2004 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

On order‐reducible sinc discretizations and block‐diagonal preconditioning methods for linear third‐order ordinary differential equations

By introducing a variable substitution, we transform the two‐point boundary value problem of a third‐order ordinary differential equation into a system of two second‐order ordinary differential equations (ODEs). We discretize this order‐reduced system of ODEs by both sinc‐collocation and sinc‐Galerkin methods, and average these two discretized linear systems to obtain the target system of linear equations. We prove that the discrete solution resulting from the linear system converges exponentially to the true solution of the order‐reduced system of ODEs. The coefficient matrix of the linear system is of block two‐by‐two structure, and each of its blocks is a combination of Toeplitz and diagonal matrices. Because of its algebraic properties and matrix structures, the linear system can be effectively solved by Krylov subspace iteration methods such as GMRES preconditioned by block‐diagonal matrices. We demonstrate that the eigenvalues of certain approximation to the preconditioned matrix are uniformly bounded within a rectangle on the complex plane independent of the size of the discretized linear system, and we use numerical examples to illustrate the feasibility and effectiveness of this new approach. Copyright © 2013 John Wiley & Sons, Ltd.     

An inverse boundary value problem for the Oseen equation

Based on the two‐dimensional stationary Oseen equation we consider the problem to determine the shape of a cylindrical obstacle immersed in a fluid flow from a knowledge of the fluid velocity on some arc outside the obstacle. First, we obtain a uniqueness result for this ill‐posed and non‐linear inverse problem. Then, for the approximate solution we propose a regularized Newton iteration scheme based on a boundary integral equation of the first kind. For a foundation of Newton‐type methods we establish the Fréchet differentiability of the solution to the Dirichlet problem for the Oseen equation with respect to the boundary and investigate the injectivity of the linearized mapping. Some numerical examples for the feasibility of the method are presented. Copyright © 2000 John Wiley & Sons, Ltd.     

Analytical solution for the evolution of exit point along the sloping seepage face

This paper developed an analytical solution for the problem of exit point evolution on the seepage face in the unconfined aquifer with sloping interface. A theoretical model for the groundwater drawdown problem in a half‐infinite aquifer with a sloping boundary is built in accordance with the linearized one‐dimensional Boussinesq equation and the Neumann boundary condition at the seepage point. The homotopy analysis method is then adopted for solving this dynamic boundary problem. By constructing two continuous deformations, the original problem could be converted into a group of subproblems with the same physical essence and similar mathematical solutions. To compare this analytical solution, a numerical model based on the finite volume method is developed, which employs adaptive grids to settle the dynamic boundary condition. The comparisons show that the analytical solution agrees with the numerical model well. The results are useful for the quantification of various hydrological problems. The methodology applied in this study is referential for other dynamic boundary problems as well.     

An iterative method based on equation decomposition for the fourth‐order singular perturbation problem

In this article, we propose an iterative method based on the equation decomposition technique ( 1 ) for the numerical solution of a singular perturbation problem of fourth‐order elliptic equation. At each step of the given method, we only need to solve a boundary value problem of second‐order elliptic equation and a second‐order singular perturbation problem. We prove that our approximate solution converges to the exact solution when the domain is a disc. Our numerical examples show the efficiency and accuracy of our method. Our iterative method works very well for singular perturbation problems, that is, the case of 0 < ε ? 1, and the convergence rate is very fast. © 2012 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2013     

A wavelets method for plane elasticity problem

A Neumann boundary value problem of plane elasticity problem in the exterior circular domain is reduced into an equivalent natural boundary integral equation and a Poisson integral formula with the DtN method. Using the trigonometric wavelets and Galerkin method, we obtain a fast numerical method for the natural boundary integral equation which has an unique solution in the quotient space. We decompose the stiffness matrix in our numerical method into four circulant and symmetrical or antisymmetrical submatrices, and hence the solution of the associated linear algebraic system can be solved with the fast Fourier transform (FFT) and the inverse fast Fourier transform (IFFT) instead of the inverse matrix. Examples are given for demonstrating our method has good accuracy of our method even though the exact solution is almost singular.     

Properties and methods for finding the best rank-one approximation to higher-order tensors

The problem of finding the best rank-one approximation to higher-order tensors has extensive engineering and statistical applications. It is well-known that this problem is equivalent to a homogeneous polynomial optimization problem. In this paper, we study theoretical results and numerical methods of this problem, particularly focusing on the 4-th order symmetric tensor case. First, we reformulate the polynomial optimization problem to a matrix programming, and show the equivalence between these two problems. Then, we prove that there is no duality gap between the reformulation and its Lagrangian dual problem. Concerning the approaches to deal with the problem, we propose two relaxed models. The first one is a convex quadratic matrix optimization problem regularized by the nuclear norm, while the second one is a quadratic matrix programming regularized by a truncated nuclear norm, which is a D.C. function and therefore is nonconvex. To overcome the difficulty of solving this nonconvex problem, we approximate the nonconvex penalty by a convex term. We propose to use the proximal augmented Lagrangian method to solve these two relaxed models. In order to obtain a global solution, we propose an alternating least eigenvalue method after solving the relaxed models and prove its convergence. Numerical results presented in the last demonstrate, especially for nonpositive tensors, the effectiveness and efficiency of our proposed methods.