Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Computing Linear Combinations of $\varphi$-Function Actions for Exponential Integrators

Awad H. Al-Mohy

TL;DR

This paper addresses the challenge of efficiently evaluating linear combinations of $φ$-function actions, $w_i=\sum_{j=0}^{p} α_i^{j} \, φ_j(t_i A) v_j$, within exponential integrators. It develops a matrix-free approach that combines scaling and recovering with a truncated Taylor series, and introduces a data-driven scheme to select a scaling parameter $s$ and a spectral shift $ξ^*$, reusing these across related calls. The authors present two algorithmic variants: a single-combination method and a block version that computes multiple combinations in parallel, decoupling the stage times from polynomial weights to improve performance on cache- and BLAS-efficient implementations. Thorough numerical experiments demonstrate near-machine accuracy for small steps and robust behavior for large steps across stiff and highly nonnormal matrices, often outperforming Krylov-based methods in reliability with competitive costs. The work offers a practical, scalable tool for exponential integrators with predictable cost and strong robustness, suitable for large-scale, matrix-free simulations.

Abstract

We propose a matrix-free algorithm for evaluating linear combinations of $\varphi$-function actions, $w_i := \sum_{j=0}^{p} α_i^{\,j}\,\varphi_j(t_i A)v_j$ for $i=1\colon r$, arising in exponential integrators. The method combines the scaling and recovering method with a truncated Taylor series, choosing a spectral shift and a scaling parameter by minimizing a power-based objective of the shifted operator. Accuracy is user-controlled and ultimately limited by the working precision. The algorithm decouples the stage abscissae $t_i$ from the polynomial weights $α_i^j$, and a block variant enables simultaneous evaluation of $\{w_i\}_{i=1}^r$. Across standard benchmarks, including stiff and highly nonnormal matrices, the algorithm attains near-machine accuracy (IEEE double precision in our tests) for small step sizes and maintains reliable accuracy for larger steps where several existing Krylov-based algorithms deteriorate, providing a favorable balance of reliability and computational cost.

Computing Linear Combinations of $\varphi$-Function Actions for Exponential Integrators

TL;DR

This paper addresses the challenge of efficiently evaluating linear combinations of -function actions, , within exponential integrators. It develops a matrix-free approach that combines scaling and recovering with a truncated Taylor series, and introduces a data-driven scheme to select a scaling parameter and a spectral shift , reusing these across related calls. The authors present two algorithmic variants: a single-combination method and a block version that computes multiple combinations in parallel, decoupling the stage times from polynomial weights to improve performance on cache- and BLAS-efficient implementations. Thorough numerical experiments demonstrate near-machine accuracy for small steps and robust behavior for large steps across stiff and highly nonnormal matrices, often outperforming Krylov-based methods in reliability with competitive costs. The work offers a practical, scalable tool for exponential integrators with predictable cost and strong robustness, suitable for large-scale, matrix-free simulations.

Abstract

We propose a matrix-free algorithm for evaluating linear combinations of -function actions, for , arising in exponential integrators. The method combines the scaling and recovering method with a truncated Taylor series, choosing a spectral shift and a scaling parameter by minimizing a power-based objective of the shifted operator. Accuracy is user-controlled and ultimately limited by the working precision. The algorithm decouples the stage abscissae from the polynomial weights , and a block variant enables simultaneous evaluation of . Across standard benchmarks, including stiff and highly nonnormal matrices, the algorithm attains near-machine accuracy (IEEE double precision in our tests) for small step sizes and maintains reliable accuracy for larger steps where several existing Krylov-based algorithms deteriorate, providing a favorable balance of reliability and computational cost.

Paper Structure

This paper contains 11 sections, 3 theorems, 31 equations, 3 figures, 4 tables, 3 algorithms.

Table of Contents

  1. Keywords.
  2. Introduction
  3. Method derivation
  4. Error analysis and parameter selection
  5. Algorithm
  6. Numerical experiment
  7. First experiment.
  8. Second experiment.
  9. Third experiment.
  10. Fourth experiment: Advection–Diffusion–Reaction (ADR).
  11. Conclusions.

Key Result

Theorem 1

Let $A\in\mathbb{C}^{n \times n}$, $V=\in\mathbb{C}^{n\times (p+1)}$, $J=[0]\oplus J_p(0)$, where $J_p(0)\in\mathbb{C}^{p\times p}$ is the Jordan block associated with zero, and Then the recurrence yields where the $j$th column of $F_s$ is $F_s(:,j) = \sum_{k=1}^{j-1}\varphi_k(A)v_{p-j+k+1}$, $2\le j\le p+1$. Consequently, we have $\mathrm{e}^Av_0 = AF_s(:,1)+v_0$.

Figures (3)

  • Figure 1: Relative forward errors for the computation $w=\sum_{j=0}^p\varphi_j(A)v_j$ using phimv and phi_funm. The solid line represents the condition number, $\kappa$, multiplied by $\mathrm{tol}=2^{-53}$.
  • Figure 2: The data in Figure \ref{['fig.test1']} presented as a performance profile.
  • Figure 3: ADR test with expRK4s6: execution time (s) vs. $\|\cdot\|_\infty$-relative error for bamphi, kiops, and phimv, all with $tol=10^{-7}$. Time steps are $h_i=2^{-8}(t_{\mathrm{end}}-t_0)/(i+1)$, $i=1,\ldots,7$. Reference solution: ode15s with RelTol$=$AbsTol$=10^{-12}$.

Theorems & Definitions (6)

  • Theorem 1
  • proof
  • Corollary 2
  • proof
  • Proposition 3
  • proof