Matrices, Moments and Quadrature: Applications to Time-Dependent Partial Differential Equations

The numerical solution of a time-dependent PDE generally involves the solution of a stiff system of ODEs arising from spatial discretization of the PDE. There are many methods in the literature for solving such systems, such as exponential propagation iterative (EPI) methods, that rely on Krylov projection to compute matrix function-vector products. Unfortunately, as spatial resolution increases, these products require an increasing number of Krylov projection steps, thus drastically increasing computational expense. This paper describes a modification of EPI methods that uses Krylov subspace spectral (KSS) methods, to compute these matrix function-vector products. KSS methods represent a balance between the efficiency of explicit methods and the stability of implicit methods. This balance is achieved by approximating the matrix exponential with different polynomials for each Fourier coefficient of the solution. These polynomials arise from techniques due to Golub and Meurant for computing bilinear forms involving matrix functions by treating them as Riemann-Stieltjes integrals, which are then approximated using Gaussian quadrature rules. This paper describes how the nodes for the quadrature rules required by KSS methods can be estimated very rapidly through asymptotic analysis of block Lanczos iteration, thus drastically reducing computational expense without sacrificing accuracy. Numeri‐ cal experiments demonstrate that this modification causes the number of Krylov projection steps to become bounded independently of the grid size, thus dramatically improving efficiency and scalability.