Krylov Approximation of Linear ODEs with Polynomial Parameterization

Abstract

Abstract. We propose a new numerical method to solve linear ordinary differential equations of the type ∂u ∂t (t, ε) = A(ε) u(t, ε), where A : C → C n×n is a matrix polynomial with large and sparse matrix coefficients. The algorithm computes an explicit parameterization of approximations of u(t, ε) such that approximations for many different values of ε and t can be obtained with a very small additional computational effort. The derivation of the algorithm is based on a reformulation of the parameterization as a linear parameter-free ordinary differential equation and on approximating the product of the matrix exponential and a vector with a Krylov method. The Krylov approximation is generated with Arnoldi's method and the structure of the coefficient matrix turns out to be independent of the truncation parameter so that it can also be interpreted as Arnoldi's method applied to an infinite dimensional matrix. We prove the superlinear convergence of the algorithm and provide a posteriori error estimates to be used as termination criteria. The behavior of the algorithm is illustrated with examples stemming from spatial discretizations of partial differential equations.

On a generalization of the Bessel function Neumann expansion

“…From the relations (6), we easily verify that the Bessel functions of the first kind are solutions to the infinite-dimensional ODE of the form (2), with H ∞ defined by (7). More precisely,…”

“…Lemma 2 in [6]). The conditions for the basis functions in (2) are satisfied with C = 2 for, (a) scaled monomials, i.e., ϕ i (t) = t!/i!, with H ∞ defined by (5); (b) Bessel functions, i.e., ϕ i (t) = J i (t), with H ∞ defined by (7); and (c) modified Bessel functions, i.e., ϕ i (t) = I i (t), with H ∞ defined by (9).…”

“…with α = 1/2. We use the Hessenberg matrices corresponding to scaled monomials (5), Bessel functions (7), modified Bessel functions (9), and an artifically constructed Hessenberg matrix given by H + αI where H is given by (7). The infinite matrix H + αI is bounded in Euclidean norm since H + αI ≤ H + |α|.…”

“…Although our work is of general character, it originated from the need of such expansions in iterative methods. More precisely, the expansion occurs in the infinite Arnoldi method both for linear ODEs [7,6] and similarly for nonlinear eigenvalue problems [5], where the fact that H ∞ is a Hessenberg matrix leads to the property that the m does not have to be determined a priori.…”