Exponential integrators based on contour integral transforms lead to powerful numerical solvers for a variety of ODEs, PDEs, and other time-evolution equations. They are easy to parallelize and lead to global-in-time approximations that can be efficiently evaluated anywhere within a finite time horizon. However, there are theoretical challenges that restrict their use-cases to classes of smooth evolution equations called analytic semigroups. In this talk, we show how to use carefully regularized contour integral representations to construct high-order quadrature schemes for the much larger class of strongly continuous semigroups. Our algorithms are accompanied by explicit high-order error bounds and near-optimal parameter selection. We illustrate the method’s attractive features through several PDE examples associated with singular behavior, causality, and non-normality. Our approach is firmly rooted in contour integral techniques for matrix functions. Along the way, we highlight how simple ideas from semigroup theory can augment traditional techniques in numerical linear algebra to tackle common tensions that arise while computing functions of operators.

Contour integral methods. Contour integral methods approximate the action of the matrix exponential with a quadrature rule. Given \(\gamma \), a Jordan curve winding once clockwise around the spectrum of a square matrix \(A\), quadrature nodes \(z_1,\ldots ,z_N\), and weights \(w_1,\ldots ,w_N\), then

\begin{equation} \label {eqn:matrix_exp} \exp (At)x = \frac {1}{2\pi i}\int _\gamma e^{zt}(zI-A)^{-1}x\, dz \approx \frac {1}{2\pi i}\sum _{k=1}^N w_k e^{z_k t} (z_kI-A)^{-1}x. \end{equation}

The main cost of the scheme is solving a shifted linear equation at each quadrature node. These can be done in parallel and may be accelerated with iterative methods. After solving these linear systems, the approximation can be evaluated at any time \(t>0\) with only vector-vector operations.

When \(A\) is obtained by discretizing a differential operator, its spectrum may fill a large region of the complex plane. As the discretization is refined, this region typically grows and the quadrature approximation in (1) may rapidly deteriorate along with the performance of iterative methods for the linear systems. Famously, this effect can be mitigated through Talbot quadrature if the spectrum of \(A\) lies in a modest sector along the negative real axis [1]. These techniques have been used to develop highly parallelizable schemes for a broad class of parabolic equations [2]. Recently, Colbrook was awarded the Leslie Fox Prize for extending these techniques to analytic semigroups, which are associated with smooth evolution problems like parabolic and damped wave equations [3].

Regularization. Fundamentally, Talbot contour integral schemes exploit regularity in the operator exponential for analytic semigroups. Unfortunately, this is not possible for less regular evolution equations associated with challenging causal, singular, or non-normal behavior. Instead, we show how to regularize contour integral schemes to exploit regularity in the vector \(x\). This is common in time-stepping methods but global-in-time methods require new tools from semigroup theory.

To exploit regularity in \(x\) in our scheme and analysis, we develop the algorithm at the infinite-dimensional level and then carefully assess the impact of discretization. Consider a linear operator \(\mathcal {A}:D(\mathcal {A})\rightarrow \mathcal {X}\) on a separable Banach space \(\mathcal {X}\) and replace \(x\) by \(u\in \mathcal {X}\) and \(A\) by \(\mathcal {A}\). If \(\mathcal {A}\) generates a strongly continuous semigroup on \(\mathcal {X}\), then its spectrum is contained in a left half-plane \({\rm Im}(z)<\omega \) and its resolvent decays uniformly in the complimentary right-half plane. We use a regularized analogue of the contour integral in (1) with form, for some \(\delta >\omega \) and integer \(m\geq 2\),

\begin{equation} \label {eqn:reg_exp} \exp (\mathcal {A}t)u = (2\delta \mathcal {I-A})^m\frac {1}{2\pi i}\int _{\delta -i\infty }^{\delta +i\infty } \frac {e^{zt}}{(2\delta -z)^m}(z\mathcal {I-A})^{-1} x\,dz, \qquad u\in D(\mathcal {A}^2). \end{equation}

To approximate (2), we apply a truncated \(N\)-point trapezoidal rule with node spacing \(h>0\),

\begin{equation} \label {eqn:reg_quad} S_N^{(m)}(t)u = (2\delta \mathcal {I-A})^m \left [\frac {h}{2\pi }\sum _{k=-N}^N \frac {e^{(\delta +ihk)t}}{(\delta -ihk)^m} ((\delta +ihk)\mathcal {I-A})^{-1} \right ]u. \end{equation}

For a practical computational scheme, the shifted linear equations can be solved with a suitable discretization of smooth functions \(u\in D(\mathcal {A}^2)\subset \mathcal {X}\). The main requirement for convergence is that \(Ax\rightarrow \mathcal {A}u\) and \((z_kI-A)^{-1}x\rightarrow (z_k\mathcal {I-A})^{-1}u\), \(k=-N,\ldots ,N\), as the discretization is refined.

Explicit Error Bounds. At the operator level, the regularized quadrature approximation in (2)-(3) has several key advantages over (1). The integral converges absolutely and uniformly for \(u\in D(\mathcal {A}^2)\), while the amplification power of \((2\delta \mathcal {I-A})^m\) is controlled by the regularity of the vector \(u\). Moreover, the integrand decays at a controlled rate along the contour. These features allow us to derive simple explicit error bounds for smooth vectors \(u\in D(\mathcal {A}^m)\). In fact, we can derive closed-form expressions for quadrature parameters that achieve [4]

\begin{equation} \sup _{0\leq t\leq T}\|\exp (\mathcal {A}t)u - S_N^{(m)}(t)u\| \leq C \left (\frac {\delta }{hN}\right )^{m-1} ||(2\delta \mathcal {I-A})^m\|, \qquad \text {when}\qquad u\in D(A^m). \end{equation}

Here, \(C\) is an explicit constant depending on \(A\), the contour location \(\delta \), and the time horizon \(T\).

In practice, one must discretize the infinite-dimensional objects for a numerical algorithm. We show how semigroup theory provides simple, computable bounds for the action of the discretized resolvent that hold for a very broad class of discretization techniques including Galerkin-based schemes and modern spectral methods. These bounds are inspired by the residual bounds for linear systems commonly used in numerical linear algebra.

Applications and Outlook. High-order contour integral schemes for strongly continuous semigroups open the door to highly parallelizable methods for traditionally challenging PDEs and related time-evolution processes. They also provide a compact modal-like representation of semigroups that may be useful in operator learning, system identification, and other data-driven modeling techniques. We will illustrate some applications of our scheme to challenging simulation scenarios related to uncertainty quantification, where verified simulation of Koopman semigroups is a key ingredient. We will also discuss work in progress on applications of contour integral representations in data-driven settings. Finally, we will outline the potential of semigroup theory to resolve related tensions encountered while computing matrix functions of discretized operators.

• [1]  A. Talbot, The accurate numerical inversion of Laplace transforms, IMA J. Appl. Math., (1979).

• [2]  L.N. Trefethen, J.A.C. Weideman, and T. Schmelzer, Talbot quadratures and rational approximations, BIT Numer. Math., (2006).

• [3]  M.J. Colbrook, Computing semigroups with error control, SIAM J. Numer. Anal., (2022).

• [4]  A. Horning and A.R. Gerlach., A family of high-order accurate contour integral methods for strongly continuous semigroups, Arxiv preprint arXiv:2408.07691, (2024).