Search

This paper derives a priori residual-type bounds for the Arnoldi approximation of a matrix function together with a strategy for setting the iteration accuracies in the inexact Arnoldi approximation of matrix functions. Such results are based on the decay behavior of the entries of functions of banded matrices. Specifically, a priori decay bounds for the entries of functions of banded non-Hermitian matrices will be exploited, using Faber polynomial approximation. Numerical experiments illustrate the quality of the results. [ABSTRACT FROM AUTHOR]

Golub and Meurant describe how pairs of Gauss and Gauss–Radau quadrature rules can be applied to determine inexpensively computable upper and lower bounds for certain real-valued matrix functionals defined by a symmetric matrix. However, there are many matrix functionals for which their technique is not guaranteed to furnish upper and lower bounds. In this situation, it may be possible to determine upper and lower bounds by evaluating pairs of Gauss and anti-Gauss rules. Unfortunately, it is difficult to ascertain whether the values determined by Gauss and anti-Gauss rules bracket the value of the given real-valued matrix functional. Therefore, generalizations of anti-Gauss rules have recently been described, such that pairs of Gauss and generalized anti-Gauss rules may determine upper and lower bounds for real-valued matrix functionals also when pairs of Gauss and (standard) anti-Gauss rules do not. The available generalization requires the matrix that defines the functional to be real and symmetric. The present paper reviews available anti-Gauss and generalized anti-Gauss rules and extends them in several ways that allow applications in new situations. In particular, the genarlized anti-Gauss rules for a real-valued non-negative measure described in Pranić and Reichel (J Comput Appl Math 284:235–243, 2015) are extended to allow the estimation of the error in matrix functionals defined by a non-symmetric matrix, as well as to matrix-valued matrix functions. Modifications that give simpler formulas and thereby make the application of the rules both easier and applicable to a larger class of problems also are described. [ABSTRACT FROM AUTHOR]

Some important applicative problems require the evaluation of functions Ψ of large and sparse and/or localized matrices A. Popular and interesting techniques for computing and , where is a vector, are based on partial fraction expansions. However, some of these techniques require solving several linear systems whose matrices differ from A by a complex multiple of the identity matrix I for computing or require inverting sequences of matrices with the same characteristics for computing. Here we study the use and the convergence of a recent technique for generating sequences of incomplete factorizations of matrices in order to face with both these issues. The solution of the sequences of linear systems and approximate matrix inversions above can be computed efficiently provided that shows certain decay properties. These strategies have good parallel potentialities. Our claims are confirmed by numerical tests. [ABSTRACT FROM AUTHOR]

We investigate the rational approximation of fractional powers of unbounded positive operators attainable with a specific integral representation of the operator function. We provide accurate error bounds by exploiting classical results in approximation theory involving Padé approximants. The analysis improves some existing results and the numerical experiments proves its accuracy. [ABSTRACT FROM AUTHOR]

The product of a matrix function and a vector is used to solve evolution equations numerically. Hashimoto and Nodera (ANZIAM J 58:C149–C161, 2016) proposed the Shift-invert Rational Krylov method for computing these products. However, since matrices produced by evolution equations behave like unbounded operators in infinite-dimensional spaces, an analysis with the unbounded operator is essential. In this paper, the Shift-invert Rational Krylov method is extended to be applied to unbounded operators. [ABSTRACT FROM AUTHOR]

Matrix functions have a major role in science and engineering. One of the fundamental matrix functions, which is particularly important due to its connections with certain matrix differential equations and other special matrix functions, is the matrix gamma function. This research article focus on the numerical computation of this function. Well-known techniques for the scalar gamma function, such as Lanczos, Spouge and Stirling approximations, are extended to the matrix case. This extension raises many challenging issues and several strategies used in the computation of matrix functions, like Schur decomposition and block Parlett recurrences, need to be incorporated to make the methods more effective. We also propose a fourth technique based on the reciprocal gamma function that is shown to be competitive with the other three methods in terms of accuracy, with the advantage of being rich in matrix multiplications. Strengths and weaknesses of the proposed methods are illustrated with a set of numerical examples. Bounds for truncation errors and other bounds related with the matrix gamma function will be discussed as well. [ABSTRACT FROM AUTHOR]

A pharmacodynamic model introduced earlier in the literature for in silico prediction of rifampicin-induced CYP3A4 enzyme production is described and some aspects of the involved curve-fitting based parameter estimation are discussed. Validation with our own laboratory data shows that the quality of the fit is particularly sensitive with respect to an unknown parameter representing the concentration of the nuclear receptor PXR (pregnane X receptor). A detailed analysis of the influence of that parameter on the solution of the model's system of ordinary differential equations is given and it is pointed out that some ingredients of the analysis might be useful for more general pharmacodynamic models. Numerical experiments are presented to illustrate the performance of related parameter estimation procedures based on least-squares minimization. [ABSTRACT FROM AUTHOR]

We consider Arnoldi-like processes to obtain symplectic subspaces for Hamiltonian systems. Large dimensional systems are locally approximated by ones living in low dimensional subspaces, and we especially consider Krylov subspaces and some of their extensions. These subspaces can be utilized in two ways: by solving numerically local small dimensional systems and then mapping back to the large dimension, or by using them for the approximation of necessary functions in exponential integrators applied to large dimensional systems. In the former case one can expect an excellent energy preservation and in the latter this is so for linear systems. We consider second order exponential integrators which solve linear systems exactly and for which these two approaches are in a certain sense equivalent. We also consider the time symmetry preservation properties of the integrators. In numerical experiments these methods combined with symplectic subspaces show promising behavior also when applied to nonlinear Hamiltonian problems. [ABSTRACT FROM AUTHOR]

Two numerical algorithms for computing interval matrices containing the matrix exponential are proposed. The first algorithm is based on a numerical spectral decomposition and requires only cubic complexity under some assumptions. The second algorithm is based on a numerical Jordan decomposition and applicable even for defective matrices. Numerical results show the effectiveness and robustness of the algorithms. [ABSTRACT FROM AUTHOR]

Let a(z)=∑i∈Zaizi be a complex valued function defined for |z|=1, such that ∑i∈Z|ai|<∞; define T(a)=(ti,j)i,j∈Z+,ti,j=aj-i for i,j∈Z+, the semi-infinite Toeplitz matrix associated with the symbol a(z); let E=(ei,j)i,j∈Z+ be a compact operator in ℓp, with 1≤p≤∞. A semi-infinite matrix of the kind A=T(a)+E is said quasi-Toeplitz (QT). The problem of the computation of exp(A) or exp(A)v, with A quasi-Toeplitz and v a vector, arises in many applications. We prove that the exponential of a QT-matrix A is QT, that is, exp(A)=T(exp(a))+F where F is a compact operator in ℓp. This property allows the design of an algorithm for computing exp(A) and exp(A)v up to any precision. The case of families of n×n matrices obtained by truncating infinite QT-matrices to finite size is also considered. Numerical experiments show the effectiveness of this approach. [ABSTRACT FROM AUTHOR]

Exponential integrators are special time discretization methods where the traditional linear system solves used by implicit schemes are replaced with computing the action of matrix exponential-like functions on a vector. A very general formulation of exponential integrators is offered by the Exponential Propagation Iterative methods of Runge-Kutta type (EPIRK) family of schemes. The use of Jacobian approximations is an important strategy to drastically reduce the overall computational costs of implicit schemes while maintaining the quality of their solutions. This paper extends the EPIRK class to allow the use of inexact Jacobians as arguments of the matrix exponential-like functions. Specifically, we develop two new families of methods: EPIRK-W integrators that can accommodate any approximation of the Jacobian, and EPIRK-K integrators that rely on a specific Krylov-subspace projection of the exact Jacobian. Classical order conditions theories are constructed for these families. Practical EPIRK-W methods of order three and EPIRK-K methods of order four are developed. Numerical experiments indicate that the methods proposed herein are computationally favorable when compared to a representative state-of-the-art exponential integrator, and a Rosenbrock-Krylov integrator. [ABSTRACT FROM AUTHOR]

In this work we derive new alternatives for efficient computation of the matrix exponential which is useful when solving Linear Initial Value Problems, vibratory systems or after semidiscretization of PDEs. We focus especially on the two classes of normal and nonnegative matrices and we present intervals of applications for rational L ∞ approximations of various degrees for these types of matrices in the lines of [7, 8]. Our method relies on Remez algorithm for rational approximation while the innovation here is the choice of the starting set of non-symmetrical Chebyshev points. Only one Remez iteration is then usually enough to quickly approach the actual L ∞ approximant. [ABSTRACT FROM AUTHOR]

We describe how to perform the backward error analysis for the approximation of exp(A)v by p(s-1A)sv, for any given polynomial p(x). The result of this analysis is an optimal choice of the scaling parameter s which assures a bound on the backward error, i.e. the equivalence of the approximation with the exponential of a slightly perturbed matrix. Thanks to the SageMath package expbea we have developed, one can optimize the performance of the given polynomial approximation. On the other hand, we employ the package for the analysis of polynomials interpolating the exponential function at so called Leja-Hermite points. The resulting method for the action of the matrix exponential can be considered an extension of both Taylor series approximation and Leja point interpolation. We illustrate the behavior of the new approximation with several numerical examples. [ABSTRACT FROM AUTHOR]

Multiple orthogonal polynomials generalize standard orthogonal polynomials by requiring orthogonality with respect to several inner products. This paper discusses an application to the approximation of matrix functions and presents quadrature rules that generalize the anti-Gauss rules proposed by Laurie. [ABSTRACT FROM AUTHOR]

In this article we investigate model order reduction of large-scale systems using time-limited balanced truncation, which restricts the well known balanced truncation framework to prescribed finite time intervals. The main emphasis is on the efficient numerical realization of this model reduction approach in case of large system dimensions. We discuss numerical methods to deal with the resulting matrix exponential functions and Lyapunov equations which are solved for low-rank approximations. Our main tool for this purpose are rational Krylov subspace methods. We also discuss the eigenvalue decay and numerical rank of the solutions of the Lyapunov equations. These results, and also numerical experiments, will show that depending on the final time horizon, the numerical rank of the Lyapunov solutions in time-limited balanced truncation can be smaller compared to standard balanced truncation. In numerical experiments we test the approaches for computing low-rank factors of the involved Lyapunov solutions and illustrate that time-limited balanced truncation can generate reduced order models having a higher accuracy in the considered time region. [ABSTRACT FROM AUTHOR]

If A is a matrix with all eigenvalues in the disk | z − 1 | < 1 , the principal p th root of A can be computed by Newton's method or Halley's method. The study of Newton's method and Halley's method for the matrix p th root can be done through a study of power series expansions of some sequences of scalar functions. In this paper, we prove monotonicity results for the coefficients in the power series expansions for both Newton's method and Halley's method, and for all integer p ≥ 2 . The sign patterns of these coefficients can be seen directly from the monotonicity results. We then use these monotonicity results and their consequences to obtain some new error estimates in the matrix case, to obtain a monotonic convergence result when A is a nonsingular M -matrix, and to obtain a structure preserving result when A is a nonsingular M -matrix or a real nonsingular H -matrix with positive diagonal entries. [ABSTRACT FROM AUTHOR]

In this research, an efficient variant of mid-point iterative method is given for computing the sign of a square complex matrix having no pure imaginary eigenvalues. It is proven that the method is new and has global convergence with high order of convergence seven. To justify the effectiveness of the new scheme, several comparisons for matrices of different sizes are worked out to show that the new method is efficient. [ABSTRACT FROM AUTHOR]

Light incident on a layer of scattering material such as a piece of sugar or white paper forms a characteristic speckle pattern in transmission and reflection. The information hidden in the correlations of the speckle pattern with varying frequency, polarization and angle of the incident light can be exploited for applications such as biomedical imaging and high-resolution microscopy. Conventional computational models for multi-frequency optical response involve multiple solution runs of Maxwell’s equations with monochromatic sources. Exponential Krylov subspace time solvers are promising candidates for improving efficiency of such models, as single monochromatic solution can be reused for the other frequencies without performing full time-domain computations at each frequency. However, we show that the straightforward implementation appears to have serious limitations. We further propose alternative ways for efficient solution through Krylov subspace methods. Our methods are based on two different splittings of the unknown solution into different parts, each of which can be computed efficiently. Experiments demonstrate a significant gain in computation time with respect to the standard solvers. [ABSTRACT FROM AUTHOR]

An analogue of the Gelfand-Shilov estimate of the matrix exponential is proved for Green’s function of the problem of bounded solutions of the ordinary differential equation x′(t)-Ax(t)=f(t). [ABSTRACT FROM AUTHOR]

In this paper, we study the relative error propagation in the solution of linear autonomous ordinary differential equations with respect to perturbations in the initial value. We also consider equations with a constant forcing term and a nonzero equilibrium. The study is carried out for equations defined by normal matrices. [ABSTRACT FROM AUTHOR]

Computing the matrix exponential for large sparse matrices is essential for solving evolution equations numerically. Conventionally, the shift-invert Arnoldi method is used to compute a matrix exponential and a vector product. This method transforms the original matrix using a shift, and generates the Krylov subspace of the transformed matrix for approximation. The transformation makes the approximation converge faster. In this paper, a new method called the double-shift-invert Arnoldi method is explored. This method uses two shifts for transforming the matrix, and this transformation results in a faster convergence than the shift-invert Arnoldi method. [ABSTRACT FROM AUTHOR]

This paper deals with some numerical issues about the rational approximation to fractional differential operators provided by the Padé approximants. In particular, the attention is focused on the fractional Laplacian and on the Caputo’s derivative which, in this context, occur into the definition of anomalous diffusion problems and of time fractional differential equations (FDEs), respectively. The paper provides the algorithms for an efficient implementation of the IMEX schemes for semi-discrete anomalous diffusion problems and of the short-memory-FBDF methods for Caputo’s FDEs. [ABSTRACT FROM AUTHOR]

An efficient and robust restart strategy is important for any Krylov-based method for eigenvalue problems. The tensor infinite Arnoldi method (TIAR) is a Krylov-based method for solving nonlinear eigenvalue problems (NEPs). This method can be interpreted as an Arnoldi method applied to a linear and infinite dimensional eigenvalue problem where the Krylov basis consists of polynomials. We propose new restart techniques for TIAR and analyze efficiency and robustness. More precisely, we consider an extension of TIAR which corresponds to generating the Krylov space using not only polynomials, but also structured functions, which are sums of exponentials and polynomials, while maintaining a memory efficient tensor representation. We propose two restarting strategies, both derived from the specific structure of the infinite dimensional Arnoldi factorization. One restarting strategy, which we call semi-explicit TIAR restart, provides the possibility to carry out locking in a compact way. The other strategy, which we call implicit TIAR restart, is based on the Krylov-Schur restart method for the linear eigenvalue problem and preserves its robustness. Both restarting strategies involve approximations of the tensor structured factorization in order to reduce the complexity and the required memory resources. We bound the error introduced by some of the approximations in the infinite dimensional Arnoldi factorization showing that those approximations do not substantially influence the robustness of the restart approach. We illustrate the effectiveness of the approaches by applying them to solve large scale NEPs that arise from a delay differential equation and a wave propagation problem. The advantages in comparison to other restart methods are also illustrated. [ABSTRACT FROM AUTHOR]

A fast iterative algorithm for numerically computing an interval matrix containing the principal p th root of an n × n matrix A is proposed. This algorithm is based on a numerical spectral decomposition of A , and is applicable when a computed eigenvector matrix of A is not ill-conditioned. Particular emphasis is put on the computational efficiency of the algorithm which has only O ( n 3 + p n ) operations per iteration. The algorithm moreover verifies the uniqueness of the contained p th root. Numerical results show the efficiency of the algorithm. [ABSTRACT FROM AUTHOR]

The algorithm by Al-Mohy and Higham (2011) [2] computes an approximation to e A b for given A and b , where A is an n -by- n matrix and b is, for example, a vector of dimension n . It uses a scaling together with a truncated Taylor series approximation to the exponential of the scaled matrix. In this paper, a method is developed for estimating the roundoff error of the computed solution. An asymptotic expansion of this error for small values of the unit roundoff is the basis of the method. The roundoff error is further expressed in terms of sums of rounding errors, which occur during the computation. A second approximation to e A b , which is computed with a lower precision than the first one, is used to evaluate these rounding errors in practice. The result is an upper bound on the normwise relative roundoff error. Further, an algorithm is proposed for computing the error bound. The cost for performing this algorithm depends on the type of problem and the accuracy, which is required for the error estimate. In case that all computations are performed with standard precisions, this cost can be expressed in terms of the number of computed matrix-vector products and is bounded from above by two times the cost for computing e A b . [ABSTRACT FROM AUTHOR]

Perturbation bounds for Mostow's decomposition and the bipolar decomposition of matrices have been computed. To do so, expressions for the derivative of the geometric mean of two positive definite matrices have been derived. [ABSTRACT FROM AUTHOR]

We present an extended Krylov subspace analogue of the two-sided Lanczos method, i.e., a method which, given a nonsingular matrix A and vectors b, c with $\left \langle {{\mathbf {b}},{\mathbf {c}}}\right \rangle \neq 0$ , constructs bi-orthonormal bases of the extended Krylov subspaces ${\mathcal {E}}_{m}(A,{\mathbf {b}})$ and ${\mathcal {E}}_{m}(A^{T}\!,{\mathbf {c}})$ via short recurrences. We investigate the connection of the proposed method to rational moment matching for bilinear forms c f( A) b, similar to known results connecting the two-sided Lanczos method to moment matching. Numerical experiments demonstrate the quality of the resulting approximations and the numerical behavior of the new extended Krylov subspace method. [ABSTRACT FROM AUTHOR]

Although the principle of alternating maximization is well known in the optimization literature, it has not been used before in the context of calculating the hump of the matrix exponential. We propose a method that applies alternating maximization in this particular context, and we show that it has a number of advantages over traditional Newton-like methods. We establish convergence results that fit this context with mild assumptions than would otherwise be the case in general optimization problems. We conduct numerical tests to complement the theory and they show convergence in just a few iterations. [ABSTRACT FROM AUTHOR]

A numerical algorithm for enclosing solutions of nonlinear eigenvalue problems arising in delay differential equations is proposed. This algorithm verifies the uniqueness of the enclosed eigenvalue when it is simple, and is applicable even for multiple eigenvalues. Numerical results show properties of the algorithm. [ABSTRACT FROM AUTHOR]

We suggest a rational Krylov subspace approximation for products of matrix functions and a vector appearing in exponential integrators. We consider matrices with a field-of-values in a sector lying in the left complex half-plane. The choice of the poles for our method is suggested by a fixed rational approximation based on contour integration along a hyperbola around the sector. Compared to the fixed approximation, our rational Krylov subspace method exhibits an accelerated and more stable convergence of order O ( e − C n ) . [ABSTRACT FROM AUTHOR]

Given matrices A and B such that B = f ( A ) , where f ( z ) is a holomorphic function, we analyze the relation between the singular values of the off-diagonal submatrices of A and B . We provide a family of bounds which depend on the interplay between the spectrum of the argument A and the singularities of the function. In particular, these bounds guarantee the numerical preservation of quasiseparable structures under mild hypotheses. We extend the Dunford–Cauchy integral formula to the case in which some poles are contained inside the contour of integration. We use this tool together with the technology of hierarchical matrices ( H -matrices) for the effective computation of matrix functions with quasiseparable arguments. [ABSTRACT FROM AUTHOR]

We consider the numerical approximation of $$f(\mathcal{A})b$$ where $$b\in {\mathbb {R}}^{N}$$ and $$\mathcal A$$ is the sum of Kronecker products, that is $$\mathcal{A}=M_2 \otimes I + I \otimes M_1\in {\mathbb {R}}^{N\times N}$$ . Here f is a regular function such that $$f(\mathcal{A})$$ is well defined. We derive a computational strategy that significantly lowers the memory requirements and computational efforts of the standard approximations, with special emphasis on the exponential and on completely monotonic functions, for which the new procedure becomes particularly advantageous. Our findings are illustrated by numerical experiments with typical functions used in applications. [ABSTRACT FROM AUTHOR]

This article deduces geometric convergence rates for approximating matrix functions via inverse-free rational Krylov methods. In applications one frequently encounters matrix functions such as the matrix exponential or matrix logarithm; often the matrix under consideration is too large to compute the matrix function directly, and Krylov subspace methods are used to determine a reduced problem. If many evaluations of a matrix function of the form f ( A ) v with a large matrix A are required, then it may be advantageous to determine a reduced problem using rational Krylov subspaces. These methods may give more accurate approximations of f ( A ) v with subspaces of smaller dimension than standard Krylov subspace methods. Unfortunately, the system solves required to construct an orthogonal basis for a rational Krylov subspace may create numerical difficulties and/or require excessive computing time. This paper investigates a novel approach to determine an orthogonal basis of an approximation of a rational Krylov subspace of (small) dimension from a standard orthogonal Krylov subspace basis of larger dimension. The approximation error will depend on properties of the matrix A and on the dimension of the original standard Krylov subspace. We show that our inverse-free method for approximating the rational Krylov subspace converges geometrically (for increasing dimension of the standard Krylov subspace) to a rational Krylov subspace. The convergence rate may be used to predict the dimension of the standard Krylov subspace necessary to obtain a certain accuracy in the approximation. Computed examples illustrate the theory developed. [ABSTRACT FROM AUTHOR]

The extended Krylov subspace method is known to be very efficient in many cases in which one wants to approximate the action of a matrix function f ( A ) on a vector b , in particular when f belongs to the class of Laplace–Stieltjes functions. We prove that the Euclidean norm of the error decreases strictly monotonically in this situation when A is Hermitian positive definite. Similar results are known for the (polynomial) Lanczos method for f ( A ) b , and we demonstrate how the techniques of proof used in the polynomial Krylov case can be transferred to the extended Krylov case. [ABSTRACT FROM AUTHOR]

The need to compute expressions of the form $$u^*f(A)v$$ , where A is a large square matrix, u and v are vectors, and f is a function, arises in many applications, including network analysis, quantum chromodynamics, and the solution of linear discrete ill-posed problems. Commonly used approaches first reduce A to a small matrix by a few steps of the Hermitian or non-Hermitian Lanczos processes and then evaluate the reduced problem. This paper describes a new method to determine error estimates for computed quantities and shows how to achieve higher accuracy than available methods for essentially the same computational effort. Our methods are based on recently proposed generalized averaged Gauss quadrature formulas. [ABSTRACT FROM AUTHOR]

When using the Lanczos method to approximate $$f(A){\varvec{b}}$$ , the action of a matrix function on a vector, there is, in contrast to the solution of linear systems, no straightforward way to measure or estimate the error of the current iterate. Therefore, to be able to decide whether the desired accuracy has been reached, several different estimates and bounds for the error have been suggested, all of them specific to certain classes of functions. In this paper, we add to these results by developing a technique to compute error bounds for Stieltjes functions, using a recently suggested integral representation of the error, and we show how these bounds can be computed essentially for free, i.e., with cost independent of the iteration number and the dimension of the matrix A. [ABSTRACT FROM AUTHOR]

Golub and Meurant have shown how to use the symmetric block Lanczos algorithm to compute block Gauss quadrature rules for the approximation of certain matrix functions. We describe new block quadrature rules that can be computed by the symmetric or nonsymmetric block Lanczos algorithms and yield higher accuracy than standard block Gauss rules after the same number of steps of the symmetric or nonsymmetric block Lanczos algorithms. The new rules are block generalizations of the generalized averaged Gauss rules introduced by Spalević. Applications to network analysis are presented. [ABSTRACT FROM AUTHOR]

The Erlangian approximation of Markovian fluid queues leads to the problem of computing the matrix exponential of a subgenerator having a block-triangular, block-Toeplitz structure. To this end, we propose some algorithms which exploit the Toeplitz structure and the properties of generators. Such algorithms allow us to compute the exponential of very large matrices, which would otherwise be untreatable with standard methods. We also prove interesting decay properties of the exponential of a generator having a block-triangular, block-Toeplitz structure. [ABSTRACT FROM AUTHOR]

In recent years, a great deal of attention has been focused on exponential integrators. The important ingredient to the implementation of exponential integrators is the efficient and accurate evaluation of the so called φ-functions on a given vector. The Krylov subspace method is an important technique for this problem. For this type of method, however, restarts become essential for the sake of storage requirements or due to computational complexities of evaluating matrix function on a reduced matrix of growing size. Another problem in computing φ-functions is the lack of a clear residual notion. The contribution of this work is threefold. First, we introduce a framework of the harmonic Arnoldi method for φ-functions, which is based on the residual and the oblique projection technique. Second, we establish the relationship between the harmonic Arnoldi approximation and the Arnoldi approximation, and compare the harmonic Arnoldi method and the Arnoldi method from a theoretical point of view. Third, we apply the thick-restarting strategy to the harmonic Arnoldi method, and propose a thick-restarted harmonic Arnoldi algorithm for evaluating φ-functions. An advantage of the new algorithm is that we can compute several φ-functions simultaneously in the same search subspace after restarting. The relationship between the error and the residual of the harmonic Arnoldi approximation is also investigated. Numerical experiments show the superiority of our new algorithm over many state-of-the-art algorithms for computing φ-functions. [ABSTRACT FROM AUTHOR]

An unsteady problem is considered for a space-fractional equation in a bounded domain. A first-order evolutionary equation involves the square root of an elliptic operator of second order. Finite element approximation in space is employed. To construct approximation in time, regularized two-level schemes are used. The numerical implementation is based on solving the equation with the square root of the elliptic operator using an auxiliary Cauchy problem for a pseudo-parabolic equation. The scheme of the second-order accuracy in time is based on a regularization of the three-level explicit Adams scheme. More general problems for the equation with convective terms are considered, too. The results of numerical experiments are presented for a model two-dimensional problem. [ABSTRACT FROM PUBLISHER]

An algorithm for computing an analytic function of a matrixis described. The algorithm is intended for the case wherehas some close eigenvalues, and clusters (subsets) of close eigenvalues are separated from each other. This algorithm is a modification of some well known and widely used algorithms. A novel feature is an approximate calculation of divided differences for the Newton interpolating polynomial in a special way. This modification does not require to reorder the Schur triangular form and to solve Sylvester equations. [ABSTRACT FROM AUTHOR]

The exponential time integration, i.e., time integration which involves the matrix exponential, is an attractive tool for time domain modeling involving Maxwell’s equations. However, its application in practice often requires a substantial knowledge of numerical linear algebra algorithms, such as Krylov subspace methods. In this note we discuss exponential Krylov subspace time integration methods and provide a simple guide on how to use these methods in practice. While specifically aiming at nanophotonics applications, we intentionally keep the presentation as general as possible and consider full vector Maxwell’s equations with damping (i.e., with nonzero conductivity terms). Efficient techniques such as the Krylov shift-and-invert method and residual-based stopping criteria are discussed in detail. Numerical experiments are presented to demonstrate the efficiency of the discussed methods and their mesh independent convergence. Some of the algorithms described here are available as Octave/Matlab codes from www.math.utwente.nl/~botchevma/expm/ . [ABSTRACT FROM AUTHOR]

The stream function–vorticity formulation of the (Navier–)Stokes equations yields a coupled system of a parabolic equation for the vorticity and an elliptic equation for the stream function. The essential coupling between them occurs through the boundary conditions which in case of a Dirichlet boundary involve only the stream function. Therefore, the boundary condition for the vorticity must be derived from them and thus the vorticity equation must be coupled to the stream function equation via its boundary condition. In this paper we propose an unconditionally stable splitting scheme for the unsteady Stokes equations in a stream function–vorticity formulation, that decouples the vorticity and stream function computations at each time step. The spatial discretization is based on a finite volume discretization on (generally) unstructured Delaunay grids and corresponding Voronoi finite volume cells. A generalization of the well-known Thom vorticity boundary condition is derived for such grids and the corresponding discrete problem is decoupled by a two-step splitting scheme which results in a decoupled discrete parabolic problem for the vorticity and an elliptic problem for the stream function. Furthermore, the scheme is extended to the unsteady Navier–Stokes equations. Finally, the stability and accuracy of the resulting schemes are demonstrated on numerical examples. [ABSTRACT FROM AUTHOR]