Adaptive rational Krylov methods for exponential Runge--Kutta integrators

We consider the solution of large stiff systems of ordinary differential equations with explicit exponential Runge--Kutta integrators. These problems arise from semi-discretized semi-linear parabolic partial differential equations on continuous domains or on inherently discrete graph domains. A series of results reduces the requirement of computing linear combinations of $\varphi$-functions in exponential integrators to the approximation of the action of a smaller number of matrix exponentials on certain vectors. State-of-the-art computational methods use polynomial Krylov subspaces of adaptive size for this task. They have the drawback that the required Krylov subspace iteration numbers to obtain a desired tolerance increase drastically with the spectral radius of the discrete linear differential operator, e.g., the problem size. We present an approach that leverages rational Krylov subspace methods promising superior approximation qualities. We prove a novel a-posteriori error estimate of rational Krylov approximations to the action of the matrix exponential on vectors for single time points, which allows for an adaptive approach similar to existing polynomial Krylov techniques. We discuss pole selection and the efficient solution of the arising sequences of shifted linear systems by direct and preconditioned iterative solvers. Numerical experiments show that our method outperforms the state of the art for sufficiently large spectral radii of the discrete linear differential operators. The key to this are approximately constant rational Krylov iteration numbers, which enable a near-linear scaling of the runtime with respect to the problem size

Improving Efficiency of Rational Krylov Subspace Methods

This thesis studies two classes of numerical linear algebra problems, approximating the product of a function of a matrix with a vector, and solving the linear eigenvalue problem $Av=\lambda Bv$ for a small number of eigenvalues. These problems are solved by rational Krylov subspace methods (RKSM). We present several improvements in two directions: pole selection and applying inexact methods. In Chapter 3, a flexible extended Krylov subspace method ($\mathcal{F}$-EKSM) is considered for numerical approximation of the action of a matrix function $f(A)$ to a vector $b$, where the function $f$ is of Markov type. $\mathcal{F}$-EKSM has the same framework as the extended Krylov subspace method (EKSM), replacing the zero pole in EKSM with a properly chosen fixed nonzero poles. For symmetric positive definite matrices, the optimal fixed pole is derived for $\mathcal{F}$-EKSM to achieve the lowest possible upper bound on the asymptotic convergence factor, which is lower than that of EKSM. The analysis is based on properties of Faber polynomials of $A$ and $(I-A/s)^{-1}$. For large and sparse matrices that can be handled efficiently by LU factorizations, numerical experiments show that $\mathcal{F}$-EKSM and a variant of RKSM based on a small number of fixed poles outperform EKSM in both storage and runtime, and they usually have advantage over adaptive RKSM in runtime. Chapter 4 concerns the theory and development of inexact RKSM for approximating the action of a function of matrix $f(A)$ to a column vector $b$. At each step of RKSM, a shifted linear system of equations needs to be solved to enlarge the subspace. For large-scale problems, arising from discretizations of PDEs in 3D domains, such a linear system is usually solved by an iterative method approximately. The main question is how to relax the accuracy of these linear solves without negatively affecting the convergence for approximating $f(A)b$. Our insight into this issue is obtained by exploring the residual bounds on the rational Krylov subspace approximations to $f(A)b$, based on the decaying behavior of the entries in the first column of the matrix function of the block Rayleigh quotient of $A$ with respect to the rational Krylov subspaces. The decay bounds on these entries for both analytic functions and Markov functions can be efficiently and accurately evaluated by appropriate quadrature rules. A heuristic based on these bounds is proposed to relax the tolerances of the linear solves arising from each step of RKSM. As the algorithm progresses toward convergence, the linear solves can be performed with increasingly lower accuracy and computational cost. Numerical experiments for large nonsymmetric matrices show the effectiveness of the tolerance relaxation strategy for the inexact linear solves of RKSM. In Chapter 5, inexact RKSM are studied to solve large-scale nonsymmetric eigenvalue problems. Similar to the problem setting in Chapter 4, each iteration (outer step) of RKSM requires solution to a shifted linear system to enlarge the subspace, but these linear solves by direct methods are prohibitive due to the problem scale. Errors are introduced at each outer step if these linear systems are solved approximately by iterative methods (inner step), and these errors accumulate in the rational Krylov subspace. In this thesis, we derive an upper bound on the errors that can be introduced at each outer step to maintain the same convergence as exact RKSM for approximating an invariant subspace. Since this bound is inversely proportional to the current eigenresidual norm of the desired invariant subspace, the tolerance of iterative linear solves at each outer step can be relaxed with the outer iteration progress. A restarted variant of the inexact RKSM is also proposed. Numerical experiments show the effectiveness of relaxing the inner tolerance to save computational cost

Spatio-temporal integral equation methods with applications

Electromagnetic interactions are vital in many applications including physics, chemistry, material sciences and so on. Thus, a central problem in physical modeling is the electromagnetic analysis of materials. Here, we consider the numerical solution of the Maxwell equation for the evolution of the electromagnetic field given the charges, and the Newton or Schr\\"odinger equation for the evolution of particles. By combining integral equation techniques with new spectral deferred correction algorithms in time and hierarchical methods in space, we develop fast solvers for the calculation of electromagnetism with relaxations of the model in different scenarios. The dissertation consists of two parts, aiming to resolve the challenges in the temporal and spatial direction, respectively. In the first part, we study a new class of time stepping methods for time-dependent differential equations. The core algorithm uses the pseudo-spectral collocation formulation to discretize the Picard type integral equation reformulation, producing a highly accurate and stable representation, which is then solved via the deferred correction technique. By exploiting the mathematical properties of the formulation and the convergence procedure, we develop some new preconditioning techniques from different perspectives that are accurate, robust, and can be much more efficient than existing methods. As is typical of spectral methods, the solution to the discretization is spectral accurate and the time step-size is optimal, though the cost of solving the system can be high. Thus, the solver is particularly suited to problems where very accurate solutions are sought or large time-step is required, e.g., chaotic systems or long-time simulation. In the second part, we study the hierarchical methods with emphasis on the spatial integral equations. In the first application, we implement a parallel version of the adaptive recursive solver for two-point boundary value problem by Cilk multithreaded runtime system based on the integral equation formulation. In the second application, we apply the hierarchical method to two-layered media Helmholtz equations in the acoustic and electromagnetic scattering problems. With the method of images and integral representations, the spatially heterogeneous translation operators are derived with rigorous error analysis, and the information is then compressed and spread in a fashion similar to fast multipole methods. The preliminary results suggest that our approach can be faster than existing algorithms with several orders of magnitude. We demonstrate our solver on a number of examples and discuss various useful extensions. Preliminary results are favorable and show the viability of our techniques for integral equations. Such integral equation methods could well have a broad impact on many areas of computational science and engineering. We describe further applications in biology, chemistry, and physics, and outline some directions for future work.Doctor of Philosoph

Structure-Preserving Model Reduction of Physical Network Systems

This paper considers physical network systems where the energy storage is naturally associated to the nodes of the graph, while the edges of the graph correspond to static couplings. The first sections deal with the linear case, covering examples such as mass-damper and hydraulic systems, which have a structure that is similar to symmetric consensus dynamics. The last section is concerned with a specific class of nonlinear physical network systems; namely detailed-balanced chemical reaction networks governed by mass action kinetics. In both cases, linear and nonlinear, the structure of the dynamics is similar, and is based on a weighted Laplacian matrix, together with an energy function capturing the energy storage at the nodes. We discuss two methods for structure-preserving model reduction. The first one is clustering; aggregating the nodes of the underlying graph to obtain a reduced graph. The second approach is based on neglecting the energy storage at some of the nodes, and subsequently eliminating those nodes (called Kron reduction).</p

Engineering Education and Research Using MATLAB

MATLAB is a software package used primarily in the field of engineering for signal processing, numerical data analysis, modeling, programming, simulation, and computer graphic visualization. In the last few years, it has become widely accepted as an efficient tool, and, therefore, its use has significantly increased in scientific communities and academic institutions. This book consists of 20 chapters presenting research works using MATLAB tools. Chapters include techniques for programming and developing Graphical User Interfaces (GUIs), dynamic systems, electric machines, signal and image processing, power electronics, mixed signal circuits, genetic programming, digital watermarking, control systems, time-series regression modeling, and artificial neural networks