Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Efficient computation of the sinc matrix function for the integration of second-order differential equations

Lidia Aceto, Fabio Durastante

TL;DR

The novelty contained here is that of using a suitable rational approximation formula for the sinc matrix function to apply a rational Krylov-like approximation method with suitable choices of poles to a finite element discretization of the wave equation.

Abstract

This work deals with the numerical solution of systems of oscillatory second-order differential equations which often arise from the semi-discretization in space of partial differential equations. Since these differential equations exhibit (pronounced or highly) oscillatory behavior, standard numerical methods are known to perform poorly. Our approach consists in directly discretizing the problem by means of Gautschi-type integrators based on $\operatorname{sinc}$ matrix functions. The novelty contained here is that of using a suitable rational approximation formula for the $\operatorname{sinc}$ matrix function to apply a rational Krylov-like approximation method with suitable choices of poles. In particular, we discuss the application of the whole strategy to a finite element discretization of the wave equation.

Efficient computation of the sinc matrix function for the integration of second-order differential equations

TL;DR

The novelty contained here is that of using a suitable rational approximation formula for the sinc matrix function to apply a rational Krylov-like approximation method with suitable choices of poles to a finite element discretization of the wave equation.

Abstract

This work deals with the numerical solution of systems of oscillatory second-order differential equations which often arise from the semi-discretization in space of partial differential equations. Since these differential equations exhibit (pronounced or highly) oscillatory behavior, standard numerical methods are known to perform poorly. Our approach consists in directly discretizing the problem by means of Gautschi-type integrators based on matrix functions. The novelty contained here is that of using a suitable rational approximation formula for the matrix function to apply a rational Krylov-like approximation method with suitable choices of poles. In particular, we discuss the application of the whole strategy to a finite element discretization of the wave equation.
Paper Structure (18 sections, 3 theorems, 79 equations, 8 figures, 4 tables)

This paper contains 18 sections, 3 theorems, 79 equations, 8 figures, 4 tables.

Table of Contents

  1. Introduction
  2. Notation
  3. Gautschi-type methods
  4. Rational Krylov methods
  5. Four approximations to the sinc function
  6. Padé-type approximants
  7. Padé approximants for the exponential function
  8. Padé approximants for
  9. Exponential sums from the inverse Fourier transform
  10. A method of lines for the linear wave equation
  11. FEM approximation spaces
  12. Numerical experiments
  13. Benchmark of the pole selection
  14. A synthetic example
  15. Solution of the linear wave equation
  16. ...and 3 more sections

Key Result

Theorem 1

Let $g$ be analytic in a neighborhood of a compact set $\Sigma \supseteq W(A)$ the field-of-value of $A$. Then the rational Krylov approximation eq:generic_krylov_approximation defined using the space eq:unseful_krylov satisfies with a constant $2 \leq C \leq 1 + \sqrt{2}$. If $A$ is Hermitian, the result holds even with $C = 1$ and $\Sigma \supseteq \Lambda(A) \cup \Lambda(V_k^T A V_k)$, for $\L

Figures (8)

  • Figure 1: Absolute error for the approximation \ref{['eq:sinc-exp-approximation']} and for the $[10/10]$-Padé approximation of the $\operatorname{sinc}$ function computed numerically in the $[-17,17]\times[-17i,17i]$ region. The color-map is $\log_{10}$-scale and the bold black dots represent the poles in the two cases; for the approximation \ref{['eq:sinc-exp-approximation']} the pole in $0$ is denoted in white to make it visible against its surrounding.
  • Figure 2: Absolute error for the approximation \ref{['eq:approx-pade-laguerre']} and for the $[10/10]$-Padé approximation of the $\operatorname{sinc}$ function computed numerically in the $[-17,17]\times[-17i,17i]$ region. The color-map is $\log_{10}$-scale and the bold black dots represent the poles.
  • Figure 3: Relative error for the computation of $\operatorname{sinc}(A)\mathbf{v}$ with $A$ the matrices described in Table \ref{['tab:test-matrices']} and $\mathbf{v}$ a randomly generated vector. We test all the pole choices described in Section \ref{['sec:rational-approx']}.
  • Figure 4: Computation of the exponential sums \ref{['eq:fourierapprox']} obtained using different approximations for the exponential, relative error versus time (s) graphs. To have a comparison with the results in Figure \ref{['fig:poles-benchmark']} we report in both pictures the results obtained with the $[n/n]$-Padé approximation, and the rational Krylov methods with poles $\mathcal{E}_n$, and $\overline{\mathcal{L}}_n$. The number of poles for the $[k/k]$-Padé approximation of the exponential is $k = 15$ for the 1D FD case and $k = 20$ for the 2D FD and Linear FEM cases. The number of Gauss-Legendre quadrature nodes goes from $\nu = 1,\ldots,15$ in all the cases.
  • Figure 5: Relative $2$-norm error at final time $T = 1$, for matrix size $N = 20$. We use the bounds from Proposition \ref{['pro:how-many-poles']} to determine the number of poles we need to achieve a given tolerance. The $\| \cdot \|_\Sigma$ norm is estimated by evaluating the $\| \cdot \|_\infty$ bound on an equally spaced grid of the interval $[0,h^2 \lambda_{\text{max}}(A)] = [0,h^2 ]$.
  • ...and 3 more figures

Theorems & Definitions (8)

  • Remark 1
  • Theorem 1: Near optimality, MR3095912MR3666309
  • Remark 2
  • Proposition 2
  • proof
  • Proposition 3
  • proof
  • Remark 3