Automatic smoothness detection of the resolvent Krylov subspace method for the approximation of $C_0$-semigroups

TL;DR

This paper introduces an automatic smoothness detection feature in the resolvent Krylov subspace method, enabling faster convergence for smoother vectors without user intervention, demonstrated through numerical experiments.

Contribution

It presents a novel property of the resolvent Krylov method that accelerates convergence based on vector smoothness, eliminating the need for user adjustments.

Findings

Faster convergence for smoother vectors

No user modifications needed for smoothness adaptation

Numerical experiments validate the theoretical results

Abstract

The resolvent Krylov subspace method builds approximations to operator functions $f(A)$ times a vector $v$. For the semigroup and related operator functions, this method is proved to possess the favorable property that the convergence is automatically faster when the vector $v$ is smoother. The user of the method does not need to know the presented theory and alterations of the method are not necessary in order to adapt to the (possibly unknown) smoothness of $v$. The findings are illustrated by numerical experiments.

Equations

$$u(t)=e^{tA}v,v\in X,t\ge 0.$$

$$u^{\prime }(t)=Au(t),u(0)=v,t\in [0,\infty ),$$

$$u^{\prime }(t)=i\frac{{\partial }^2}{\partial x^2}u(t),u(0)=u_0\text{for}x\in (0,2\pi ),t\ge 0.$$

$$\sum _{k=-N/2}^{N/2-1}{\psi }_k^{\prime }(t){\varphi }_k=\sum _{k=-N/2}^{N/2-1}(-i{k}^2){\psi }_k(t){\varphi }_k.$$

$${\Psi }^{\prime }(t)=i{A}_N\Psi (t),\Psi (0)={\Psi }_0$$

$${\psi }_k(0)=\frac{1}{2\pi }{\int }_0^{2\pi }{u}_0(x)e^{-ikx}dx,k=-\frac{N}{2},\dots ,\frac{N}{2}-1.$$

$$u^{q}_{0}(x)=\left\{\begin{array}[]{l@{\qquad}l}\left(\tfrac{2}{\pi}\right)^{4q}(x-\pi)^{2q}x^{2q}\,,&x\in(0,\pi]\,,\\[8.61108pt] \left(\tfrac{2}{\pi}\right)^{4q}(x-\pi)^{2q}(x-2\pi)^{2q}\,,&x\in(\pi,2\pi]\,.\end{array}\right.$$

$$\begin{array}[]{l@{\qquad}l@{\hspace*{0.1cm}}l}\text{explicit Euler:}&\displaystyle\lim_{n\to\infty}\left(I+\tfrac{\tau}{n}A\right)^{n}v\hfil\hskip 2.84544pt&=e^{\tau A}v\,,\\[4.30554pt] \text{implicit Euler:}&\displaystyle\lim_{n\to\infty}\left(I-\tfrac{\tau}{n}A\right)^{-n}v\hfil\hskip 2.84544pt&=e^{\tau A}v\,.\end{array}$$

$$\left\|e^{\tau A}v-\left(1-\tfrac{\tau}{n}A\right)^{-n}v\right\|\leq\left\{\begin{array}[]{l@{\qquad}l}C\,\tfrac{\tau}{\sqrt{n}}\,\|Av\|\,,&v\in\mathcal{D}(A)\,,\\[8.61108pt] C\,\tfrac{\tau^{2}}{n}\,\|A^{2}v\|\,,&v\in\mathcal{D}(A^{2})\,.\end{array}\right.$$

$$e^{A}v\approx\bigl{(}1+\tfrac{\tau}{n-1}A\bigr{)}^{n-1}v=a_{0}v+a_{1}Av+\cdots+a_{n-1}A^{n-1}v=p(A)v\,,\quad p\in\mathcal{P}_{n-1}\,,$$

$${\mathcal{K}}_n(A,v)=\text{mbox}span\{v,Av,A^{2}v,\dots ,A^{n-1}v\}.$$

$$e^{A}v\approx a_{0}v+a_1(\gamma -A{)}^{-1}v+\cdots +a_{n-1}(\gamma -A{)}^{-(n-1)}v=r(A)v,r\in \frac{{\mathcal{P}}_{n-1}}{(\gamma -\cdot {)}^{n-1}},$$

$${\mathcal{K}}_n((\gamma -A{)}^{-1},v)=\mathrm{span}\{v,(\gamma -A{)}^{-1}v,(\gamma -A{)}^{-2}v,\dots ,(\gamma -A{)}^{-n+1}v\},\gamma >0.$$

$${\mathcal{M}}_{+}=\left\{f\in C(\mathbb{R})|\mathcal{F}f\in L^1(\mathbb{R})\text{mbox}and\text{operator@font}supp(\mathcal{F}f)\subseteq [0,\infty )\right\},$$

$$\mathcal{F}f(s)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{-ixs}f(x)dx\text{mbox}forf\in L^1(\mathbb{R}).$$

$$\tilde{\mathcal{M}}:=\left\{f\text{mbox}holomorphicandboundedfor\text{mbox}Rez\le 0|f_{(0)}\in {\mathcal{M}}_{+}\right\}.$$

$$f(A)={\int }_0^{\infty }{e}^{sA}\mathcal{F}{f}_{(0)}(s)ds.$$

$$\left(\frac{1}{(z_0-z{)}^k}\right)(A)={\int }_0^{\infty }{e}^{sA}{e}^{-s{z}_0}\cdot \frac{s^{k-1}}{(k-1)!}ds=(z_0-A{)}^{-k},$$

$${\mathcal{M}}_0:=\{f\text{mbox}holomorphicfor\text{mbox}Rez\le 0|\exists n\in {\mathbb{N}}_0:\frac{f(z)}{(1-z{)}^n}\in \tilde{M}\}.$$

$$f(A):=(1-A{)}^n\left(\frac{f(z)}{(1-z{)}^n}\right)(A),$$

$$\tilde{\mathcal{M}}\subseteq \mathcal{M}:=\{f\in {\mathcal{M}}_0|f(A):X\to X\text{mbox}isbounded\}$$

$$\left(e^{\tau z}\right)(A)=e^{\tau A}.$$

$$\left(\frac{e^{z\tau }}{1-z}\right)(A)$$

$$={\int }_0^{\infty }{e}^{(s+\tau )A}{e}^{-s}ds=e^{\tau A}{\int }_0^{\infty }{e}^{sA}{e}^{-s}ds=e^{\tau A}(1-A{)}^{-1}.$$

$$(1-A)\left(\frac{e^{\tau z}}{1-z}\right)(A)=(1-A)e^{\tau A}(1-A{)}^{-1}=e^{\tau A}$$

$$f(A)v=g(A)A^{l}v\text{for}v\in \mathcal{D}(A^l).$$

$${\mathcal{K}}_n((\gamma -A{)}^{-1},v):=\mathrm{span}\{v,(\gamma -A{)}^{-1}v,(\gamma -A{)}^{-2}v,\dots ,(\gamma -A{)}^{-n+1}v\},\gamma >0.$$

$${\mathcal{K}}_{\infty }((\gamma -A{)}^{-1},v):=\mathrm{span}\{v,(\gamma -A{)}^{-1}v,(\gamma -A{)}^{-2}v,\dots \}.$$

$$f(A)v\in \overline{{\mathcal{K}}_{\infty }((\gamma -A{)}^{-1},v)}.$$

$$Y:=\overline{{\mathcal{K}}_{\infty }((\gamma -A{)}^{-1},v)},$$

$$(\mu -A{)}^{-1}y\in Y\mathrm{for}\mathrm{all}y\in Y,\mu >0.$$

Full text

1 Introduction

Let $X$ be some Banach space with norm $\|\cdot\|$. For $t\geq 0$, we consider a $C_{0}$-semigroup $e^{tA}$, which is generated by $A$, applied to some initial data $v\in X$, or more exactly,

[TABLE]

Due to a standard rescaling argument (cf. Section 2.2 on page 60 in [8]), it suffices to study bounded semigroups, that is, semigroups satisfying $\|e^{tA}\|\leq N$ for all $t\geq 0$. The object of interest (1) is just the (mild) solution of the abstract linear evolution equation

[TABLE]

whose effective approximation is important in many applications, especially for the numerical solution of semilinear evolution equations by either splitting methods (e.g. [35, 25]) or exponential integrators (e.g. [20]). In order to approximate the solution (1) of the abstract evolution equation in an efficient and reliable way, one has to use a method which leads to an error reduction that is independent of the norm of the matrix representing the discretized operator $A$ (see [18]). Such error bounds can therefore be designated as grid-independent, since the refinement of the grid in space does not deteriorate the convergence in time (cf. [11, 10]).

In the case of a matrix or a bounded operator $A$, the basic importance of an efficient approximation to $e^{tA}$ and possible methods for this problem are well reflected in “Nineteen dubious ways to compute the exponential of a matrix” [27] by Moler and van Loan. The subsequent finding that the standard Krylov subspace approximation can be used for the approximation of the matrix exponential times a vector, $e^{tA}v$, led to an updated version “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later” with the Krylov subspace method as twentieth method (see [28]). Recently, it becomes more and more apparent, that rational Krylov subspace methods constitute a promising twenty-first possibility that is even suitable for matrices with a large norm or unbounded operators. The use of rational Krylov subspaces for the approximation of matrix/operator functions $f(A)$ times $v$ has been studied and promoted, e.g., in [2, 3, 5, 4, 7, 12, 9, 11, 15, 16, 21, 23, 24, 29, 31, 30, 32, 33, 36].

In this paper, we will study the approximation of $e^{tA}v$ and products of related operator functions, the so-called $\varphi$-functions, times $v$ in the resolvent Krylov subspace spanned by $(\gamma-A)^{-1}$ and $v$. An efficient approximation of these operator functions is of major importance particularly in the context of exponential integrators. Our error analysis provides sublinear error bounds for unbounded operators $A$ that translate to error bounds independent of the norm of the discretized operator. That is, the error bounds prove a grid-independent convergence for the discretized problem. Moreover, it turns out that the error reduction correlates with the smoothness of the initial value $v$. A favorable property is that the resolvent Krylov subspace method detects the smoothness of the initial vector by itself and converges the faster the smoother $v$ is. All of this happens automatically, the user of the method does not even need to know the precise smoothness of the initial value.

After this introduction and a motivation in Section 2, we briefly review a functional calculus in Section 3. In Section 4 we prove that any function of the presented functional calculus times a vector can be approximated in the resolvent Krylov subspace spanned by the resolvent and this vector. For the proof of our main results, some smoothing operators are introduced in Section 5. In Section 6, the approximation of the semigroup in Banach spaces is considered. Our main theorems can be found in Section 7, where the effect of the smoothness of the vector on the convergence rate for the approximation of the semigroup and related functions in Hilbert spaces is studied. Some numerical illustrations of our results are given in Section 8, followed by a conclusion.

2 Motivation

For a first illustration of this nice feature of the resolvent Krylov subspace method just mentioned above, we consider the one-dimensional Schrödinger equation on $L^{2}(0,2\pi)$

[TABLE]

With $A=\tfrac{\partial^{2}}{\partial x^{2}}$, we obtain the abstract equation $u^{\prime}(t)=iAu(t)$, where the domain of $A$ is the Sobolev space $H^{2}_{\pi}(0,2\pi)$ containing all $2\pi$-periodic functions that admit a second order weak derivative. We now discretize (3) by a pseudospectral method. Therefore, we approximate the unknown solution $u$ by a finite linear combination of the basis functions $\phi_{k}(x)=e^{ikx}$, that is $u(t)\approx\sum_{k=-N/2}^{N/2-1}\psi_{k}(t)\phi_{k}$ with $N$ even, and search for coefficients $\psi_{k}(t)$ such that

[TABLE]

This ansatz is equivalent to

[TABLE]

with solution $\Psi(\tau)=e^{i\tau A_{N}}\Psi_{0}$, where the vector $\Psi(t)\in\mathbb{C}^{N}$ contains the Fourier coefficients $\psi_{k}(t)$ for $k=-\frac{N}{2},\ldots,\frac{N}{2}-1$ and the matrix $A_{N}\in\mathbb{R}^{N\times N}$ is a diagonal matrix with entries $\left(-\frac{N}{2}\right)^{2},\left(-\tfrac{N}{2}+1\right)^{2},\ldots,\left(\tfrac{N}{2}-1\right)^{2}$. The discretized initial vector \Psi_{0}=\big{(}\psi_{k}(0)\big{)} is given by

[TABLE]

These coefficients $\psi_{k}(0)$ can be approximated by a discrete Fourier transform of the discretized function $u_{0}$. Here, we use the initial data

[TABLE]

Differentiating this function $2q+1$ times, $\frac{d^{2q+1}}{dx^{2q+1}}u_{0}^{q}$ becomes discontinuous at $x=\pi$ and at $x=2\pi$, if $u_{0}^{q}$ is considered as a $2\pi$-periodic function. So, we have $u_{0}^{q}\in\mathcal{D}(A^{q})$ but $u_{0}^{q}\not\in\mathcal{D}(A^{q+1})$. By $\Psi_{0}^{q}\in\mathbb{C}^{N}$, we denote the corresponding spectral discretizations of the initial value. The solution $e^{i\tau A_{N}}\Psi_{0}^{q}$ of the discretized initial value problem at time $\tau>0$ is now approximated in the rational Krylov subspace $\mathcal{K}_{n}((\gamma-i\tau A_{N})^{-1},\Psi_{0}^{q})$, where $\gamma=1$.

In Figure 1, the error of the rational Krylov subspace approximation is plotted against the dimension of the Krylov subspace (blue solid lines) for $N=131072$, $\tau=0.02$, and smoothness indices $q=2,4,6,8$. We can observe that $e^{i\tau A_{N}}\Psi_{0}^{q}$ is approximated the better the smoother the continuous initial value $u_{0}^{q}$ is, or more exactly, the higher the number $q$ with $u_{0}^{q}\in\mathcal{D}(A^{q})$ is. Furthermore, we applied for comparison the implicit Euler method to the discretized problem and added the obtained error curves to Figure 1 (red dashed lines).

In order to introduce the resolvent Krylov subspace and to get a first idea why this subspace might be a good choice for the approximation of the operator/matrix exponential, we consider the implicit Euler scheme, which is, besides the explicit Euler scheme, a standard method to approximate the matrix exponential times a vector. The two methods are based on the relations

[TABLE]

For matrices with a large norm, the explicit Euler method does not work efficiently. The discretized Schrödinger equation is a stiff system of ordinary differential equations. If we increase the number $N$ of basis functions, the norm of the discretization matrix $A_{N}$ grows. For the discretized problem considered here, the explicit Euler method is therefore not suitable. The situation is even worse for the continuous equation, since the explicit Euler scheme cannot be used unless the initial data is very smooth and lies in $\mathcal{D}(A^{\infty})=\cap_{n=1}^{\infty}\mathcal{D}(A^{n})$. The implicit Euler method, however, provides an approximation to the semigroup for all initial vectors $v$ in the associated Banach space $X$. The resolvent $(I-\frac{\tau}{n}A)^{-1}$ maps $X$ to $\mathcal{D}(A)$ and can thus be seen as a smoothing operator. While the explicit Euler method cannot be applied for initial values $v\not\in\mathcal{D}(A^{\infty})$, the implicit Euler method can be proven to possess the convergence rates (see [6])

[TABLE]

For even smoother data, the implicit Euler method does not converge faster. In order to improve both methods, we briefly review the basic idea of Krylov subspace methods. Consider $n-1$ steps of the explicit Euler method that can be seen as a product of a polynomial in $A$ and $v$, that is

[TABLE]

where $\mathcal{P}_{n-1}$ is the space of polynomials of maximum degree $n-1$. Instead of using a fixed polynomial approximation given by the explicit Euler scheme, it might be better to search a best approximation in the polynomial Krylov subspace

[TABLE]

Using this approximation space for the numerical solution of the discretized Schrödinger equation, it turns out that the approximation improves with respect to stability, but a substantial error reduction would just begin after nearly $\|\tau A_{N}\|$ iteration steps (see [19]), where $\|\tau A_{N}\|$ becomes large for fine space discretizations. Analogous to the explicit Euler method, the standard polynomial Krylov subspace method is not suitable for a grid-independent approximation of the Schrödinger equation.

Instead of applying the implicit Euler method, one can try to find a better approximation of the type

[TABLE]

that means, to search a best approximation in the rational Krylov subspace

[TABLE]

This so-called resolvent Krylov subspace has been proposed by Ruhe in [34] for eigenvalue computations and is by now a standard technique for this purpose (cf. [1]). We will use the resolvent Krylov subspace as approximation space for the approximation of $e^{A}v$ and related operator functions in the following. Analogous to the implicit Euler method, the approximation based on the resolvent Krylov subspace will be grid-independent but improve on the implicit Euler method with respect to the convergence rate dependent on the smoothness of the vector $v$, as illustrated above (cf. Figure 1).

3 Preliminaries

We briefly review a functional calculus that has been formerly used in [11] and [14]. The Lebesgue space of complex-valued integrable functions defined on $\mathbb{R}$ is denoted by $L^{1}(\mathbb{R})$ with norm $\|\cdot\|_{1}$. By $C(\mathbb{R})$, we designate the space of continuous functions $f\,:\,\mathbb{R}\to\mathbb{C}$. Moreover, let

[TABLE]

where $\mathcal{F}f$ is the Fourier transform of $f$ given as

[TABLE]

For $f\not\in L^{1}(\mathbb{R})$, the Fourier transform is understood in the sense of distributions. For each function holomorphic in the left half-plane, we denote by $f_{(0)}:\mathbb{R}\rightarrow\mathbb{C}$ the restriction of $f$ to $\mbox{Re}\,z=0$ so that $f_{(0)}(\xi)=f(i\xi)$, $\xi\in\mathbb{R}$, and we define the algebra

[TABLE]

Let $A$ generate a bounded strongly continuous semigroup with $\|e^{\tau A}\|\leq N$ on some Banach space $X$. For functions $f\in\widetilde{\mathcal{M}}$, we introduce a functional calculus via

[TABLE]

This defines a bounded linear operator $f(A)$ satisfying $\|f(A)\|\leq N\|\mathcal{F}f_{(0)}\|_{1}$. Until we know that the functional calculus is consistent with standard operator functions such as the resolvent and the semigroup, we write $\left(f(z)\right)(A)$, when the definition of the operator functions is according to the new calculus (6), instead of simply $f(A)$. For $f(z)=(z_{0}-z)^{-k}$ with $\mbox{Re}\,z_{0}>0$ and $k\geq 1$, we have by elementary semigroup theory (cf. Corollary 1.11, pp. 56–57 in [8]),

[TABLE]

that is, the definition via (6) coincides with the definition in terms of the resolvent. Analogously, all rational functions with a smaller degree of the numerator than the denominator and poles in the right complex half-plane are included by our functional calculus so far.

We will need another extension in order to include the semigroup, i.e., we want that the generator $A$ inserted in the exponential function $e^{tz}$, $t\geq 0$, coincides with the semigroup. Let

[TABLE]

For $f\in\mathcal{M}_{0}$, we set

[TABLE]

where $n$ is such that $\frac{f(z)}{(1-z)^{n}}\in\widetilde{M}$. Note that the definition does not depend on the choice of $n$ and that the definition results in a closed operator on $X$. Finally, we define the set

[TABLE]

which is sufficient for our purposes. The following lemma can be found as Proposition 1.12 in [17].

Lemma 1**.**

The mapping $f\rightarrow f(A)$ via (6) is a homomorphism of $\mathcal{M}$ into the algebra of bounded linear operators on $X$.

We can check, that the semigroup is now included in the extended functional calculus.

Lemma 2**.**

For $\tau\geq 0$, we have

[TABLE]

Proof.

For $n=1$, one can verify that $\frac{e^{\tau z}}{1-z}\in\widetilde{\mathcal{M}}$. Hence, we have by (6) that

[TABLE]

Finally, we conclude

[TABLE]

which proves the assertion. ∎

For all functions relevant to our discussion, the functional calculus (6) coincides with the definitions in semigroup theory. From now on, we therefore do not use different notations and simply write $f(A)$ for a function $f$ of an operator $A$ with respect to (6) . We will also need the following lemma of Brenner and Thomée (cf. Lemma 4 in [6]), whose proof extends to our case.

Lemma 3**.**

For $f,g\in\mathcal{M}$ with $f(z)=z^{l}g(z)$ for some $l>0$ and $\mbox{Re}\,z\leq 0$, we have

[TABLE]

4 Approximation in the resolvent Krylov subspace

Here and in the following, we always consider bounded semigroups with generator $A$ on some Banach space $X$ which satisfy $\|e^{tA}\|\leq N$. For bounded semigroups, it is well-known that the right complex half-plane belongs to the resolvent set of the generator $A$ (e.g. Theorem 1.10 on page 55 in [8]) which guarantees that the resolvent $(\gamma-A)^{-1}$ exists for all $\gamma>0$.

We are interested in the approximation of operator functions, especially the semigroup, times a vector $v\in X$ in the resolvent Krylov space

[TABLE]

For $n=1,2,3,\ldots$, these spaces form a nested sequence of subspaces. If there exists an index $n_{0}$ for which $\mathcal{K}_{n_{0}}((\gamma-A)^{-1},v)$ is invariant under $(\gamma-A)^{-1}$, we have $\mathcal{K}_{n_{0}}((\gamma-A)^{-1},v)=\mathcal{K}_{k}((\gamma-A)^{-1},v)$ for all $k\geq n_{0}$. For a Banach space $X$ of finite dimension, this always happens. At the latest, when $n$ reaches the dimension of $X$. For a Banach space of infinite dimension, this might happen or it might not. In most cases, the spaces build an infinite series of nested spaces that are different. We therefore first discuss the natural question, whether all functions of our functional calculus can be approximated to an arbitrary precision in the space (7) when $n$ tends to infinity. For this purpose, we define the maximal resolvent Krylov subspace.

Definition 4**.**

The maximal resolvent Krylov subspace for a given vector $v\in X$ and a fixed $\gamma>0$ is given as the space

[TABLE]

We also need the closure of this space that we designate by $\overline{\mathcal{K}_{\infty}((\gamma-A)^{-1},v)}\subseteq X$.

The following theorem states that all functions that are defined for $A$ via (6) times $v$ are in the closure of the maximal resolvent Krylov subspace (8), that is, $f(A)v$ can be approximated in the Krylov subspace (8) to any desired precision. Since the span designates all finite linear combinations, this also means that all functions in our functional calculus can be approximated in the space (7) to any arbitrary precision, if we let $n$ go to infinity.

Theorem 5**.**

For all $v\in X$ and all functions $f\in\mathcal{M}$, we have

[TABLE]

Proof.

If we define

[TABLE]

then $Y$ is an invariant subspace of $(\mu-A)^{-1}$ for all $\mu>0$ (cf. proof of Theorem 4.6.1 in [26]). Hence, we have

[TABLE]

Theorem 4.6.1 in [26] now states that $Y$ is an invariant subspace of our semigroup $e^{tA}$, $t\geq 0$, and $A$. Furthermore, the restriction $\left.e^{tA}\right|_{Y}$ of the semigroup $e^{tA}$ to $Y$ is again a semigroup with generator $\left.A\right|_{Y}$ and $\left.A\right|_{Y}y=Ay$ for all $y\in\mathcal{D}(A)\cap Y=:\mathcal{D}(\left.A\right|_{Y})$. For $f\in\widetilde{\mathcal{M}}\subseteq\mathcal{M}$, we thus find

[TABLE]

Since $v\in Y$, we obtain $f(A)v\in Y$. Now we proceed with the case of functions belonging to $\mathcal{M}$. By the definition of $\mathcal{M}$, we have for $f\in\mathcal{M}$ and all $y\in Y$ that

[TABLE]

where $n$ has been chosen appropriately. Because of $Y$ being $A$-invariant, we obtain

[TABLE]

By the first part of the proof, since $g\in\widetilde{\mathcal{M}}$, we can conclude for $y\in Y$ that

[TABLE]

Again, due to $v\in Y$, the statement $f(A)v\in Y$ follows. ∎

In the case that an index $n_{0}$ exists for which the resolvent Krylov subspace is invariant, we obtain from Theorem 5, that

[TABLE]

Thus, $f(A)v$ can be represented exactly in the finite-dimensional space (7) with index $n=n_{0}$.

We also study subspaces of $X$ of the type

[TABLE]

with a smoothed initial vector. These spaces are usually different. For example, if $v\in X\backslash\mathcal{D}(A)$ holds true, then $v$ is in the space (8), but $v$ is not in any of the spaces (9), which are all subsets of $\mathcal{D}(A)$. An intriguing fact is that the closures of the spaces (9) are identical and coincide with the closure of (8). Hence, from a numerical analyst’s point of view, if $w\in X$ can be approximated to an arbitrary precision in any of the spaces (8) or (9), then $w$ can be approximated in all spaces to an arbitrary precision.

Lemma 6**.**

For every $q=1,2,3,\ldots$ and $v\in X$, we have

[TABLE]

Proof.

If we set

[TABLE]

then we have, analogously to the previous proof, that

[TABLE]

Obviously, $(\gamma-A)^{-1}v\in Y$ and hence $(\gamma-\mu)(\mu-A)^{-1}(\gamma-A)^{-1}v\in Y$. By the resolvent equation (cf. (1.2) on page 239 in [8]), it follows

[TABLE]

Since $\mu>0$ has been arbitrarily chosen, we have $\mu(\mu-A)^{-1}v\in Y$ for all $\mu>0$. Due to the well-known fact that

[TABLE]

(e.g. Lemma 3.4, p. 73 in [8]), we find $v\in Y$, since $\mu(\mu-A)^{-1}v\in Y$ and $Y$ is closed. This immediately shows our assertion for $q=1$. The statement for $q>1$ now follows by induction. ∎

5 Smoothing operators with range in the resolvent Krylov subspace

We study in this section, how well an initial vector $v\in\mathcal{D}(A^{q})$ can be approximated in different resolvent Krylov subspaces. These bounds will be necessary to prove our main theorems. The constants occurring in the following bounds and estimates will always be generic constants denoted by $C(param_{1},param_{2},\ldots)$, where the terms in brackets indicate the parameters on which $C$ depends.

The first lemma provides an error estimate for the approximation of $v\in\mathcal{D}(A^{q})$ in the special resolvent Krylov subspace $\mathcal{K}_{q}((\sqrt{n}-A)^{-1},(\sqrt{n}-A)^{-q}v)$.

Lemma 7**.**

There are bounded operators $H_{n,q}$ of the form

[TABLE]

such that, for all $v\in\mathcal{D}(A^{q})$, we have

[TABLE]

where the constants depend only on $q$ and the bound $N$ of the bounded semigroup, but not on $n$.

Proof.

The coefficients $h_{k}^{q}$, $k=q,\ldots,2q-1$, are chosen such that

[TABLE]

can be continued to a holomorphic function on $\mathbb{C}\backslash\{1\}$. One can check that $g(z)$ is a linear combination of $(1-z)^{-q},(1-z)^{-(q+1)},\ldots,(1-z)^{-(2q-1)}$. For any generator $B$ of a bounded semigroup $e^{tB}$ with $\|e^{tB}\|\leq N$, we have $\|(1-B)^{-1}\|\leq N$ (cf. Theorem 1.10 on page 55 in [8]) and thus

[TABLE]

We will now use this estimate for $B=\frac{1}{\sqrt{n}}A$, which generates a semigroup satisfying $\|e^{\frac{t}{\sqrt{n}}A}\|\leq N$ to bound the difference $\|H_{n,q}v-v\|$. Since $g\in\mathcal{M}$, $1\in\mathcal{M}$ and $\sum_{k=q}^{2q-1}h_{k}^{q}(1-z)^{-k}\in\mathcal{M}$ are functions that belong to our extended functional calculus, we obtain according to Lemma 3

[TABLE]

The bound on $H_{n,q}$ follows, since $\|(\sqrt{n}-A)^{-1}\|\leq N/\sqrt{n}$ (e.g. Theorem 1.10 on page 55 in [8]) and therefore $\|\sqrt{n}(\sqrt{n}-A)^{-1}\|\leq N$. ∎

In the next lemma, we study the best approximation $W_{n}$ to $\sqrt{n}(\sqrt{n}-A)^{-1}$ in the resolvent subspace

[TABLE]

Lemma 8**.**

For any $p\in\mathbb{N}$ and any $n\in\mathbb{N}$, there are operators $W_{n}$ of the form

[TABLE]

such that

[TABLE]

with $n$-independent generic constants $C(\gamma,p,N)$.

Proof.

We estimate the best approximation in the space (11). This best approximation exists, since the space is finite. By standard calculations with the resolvent, we obtain

[TABLE]

We now proceed analogously to the proof of Theorem 4.1 in [13]. By using

[TABLE]

with the $m$-th Laguerre polynomial

[TABLE]

it follows

[TABLE]

where $\mathcal{F}f_{(0)}$ is the Fourier transform of $f(z)=\frac{\sqrt{n}}{\sqrt{n}-z}$ restricted to $\text{Re}~{}z=0$. It remains to estimate

[TABLE]

Let $L_{\omega_{0}}^{1}(\mathbb{R})$ be the space of Lebesgue integrable functions with respect to the weight function $\omega_{0}$. Moreover, we equip the space $W_{r}^{1}(\omega_{0}):=\left\{g\in L_{\omega_{0}}^{1}(\mathbb{R})\,:\,\left\|g^{(r)}\varphi^{r}\omega_{0}\right\|_{1}<\infty\right\}$ with the norm

[TABLE]

According to [22], there exists for any $r\geq 1$ a constant $C(r)$ such that we have for any $g\in W_{r}^{1}(\omega_{0})$

[TABLE]

where $C(r)$ neither depends on $m$ nor on $g$. Our function

[TABLE]

is even smoother and we obtain

[TABLE]

Hence, we have

[TABLE]

and therefore with the choice $r=2p$ and $m=n$

[TABLE]

Choosing the coefficients $a_{k}$ in accordance with the coefficients $w_{k}^{n}$ that belong to a best approximation, the first statement is proved. Since

[TABLE]

the second statement is an immediate consequence. ∎

The choice $\gamma=\sqrt{n}$ in Lemma 8 obviously gives error zero, since the resolvent $\sqrt{n}(\sqrt{n}-A)^{-1}$ can be represented exactly in the resolvent subspace $\mathcal{R}_{n}(\sqrt{n},A)$.

We now continue with two further lemmas that state how well a vector $v\in\mathcal{D}(A^{q})$ can be approximated in resolvent Krylov subspaces of type (9).

Lemma 9**.**

There exist operators $\tilde{S}_{n,q}$ with

[TABLE]

such that for all $v\in\mathcal{D}(A^{q})$ and $n\in\mathbb{N}$ the bounds

[TABLE]

hold true with constants $C(\gamma,q,N)$ not depending on $n$. Moreover, we have

[TABLE]

Proof.

We choose the coefficients $h_{k}^{q}$ as in Lemma 7 and set

[TABLE]

where the $W_{n}^{k}$ are the $k$-th powers of the operators $W_{n}$ in Lemma 8. Since these operators $W_{n}$ are uniformly bounded according to Lemma 8, it is clear that the operators $\tilde{S}_{n,q}$ are uniformly bounded with respect to $n$ as well. From the definition of $W_{n}$, we conclude

[TABLE]

and hence the same holds true for $\tilde{S}_{n,q}$. For fixed $q$, let now $v\in\mathcal{D}(A^{q})$ be arbitrary. It follows

[TABLE]

with Lemma 7. Since we can write

[TABLE]

and since $W_{n}$ and $\sqrt{n}(\sqrt{n}-A)^{-1}$ are bounded independently of $n$, we obtain with Lemma 8

[TABLE]

and therefore

[TABLE]

Altogether, the first bound in our lemma is proved. Due to the special form of the operators $W_{n}$, the statement $\tilde{S}_{n,q}v\in\mathcal{K}_{(2q-1)n-(q-1)}((\gamma-A)^{-1},(\gamma-A)^{-q}v)$ is clear. ∎

Lemma 10**.**

There exist operators $S_{n,q}$ with

[TABLE]

such that for all $v\in\mathcal{D}(A^{q})$ and $n\geq 2q$, we have

[TABLE]

with $n$-independent constants and

[TABLE]

Proof.

This is an immediate result of Lemma 9 by setting $S_{n,q}:=\tilde{S}_{\lfloor\frac{n-2q}{2(2q-1)}\rfloor+1,q}$. ∎

Lemma 10 basically says that one can approximate a vector $v\in\mathcal{D}(A^{q})$ in about half of the Krylov subspace $\mathcal{K}_{n}((\gamma-A)^{-1},(\gamma-A)^{-q}v)$ and retaining a convergence according to the smoothness of $v$. The operator $S_{n,q}$ might therefore be seen as a smoothing projecting operator because of $S_{n,q}v\in\mathcal{D}(A^{2q})$ for $v\in\mathcal{D}(A^{q})$.

6 Approximation properties of the semigroup in Banach spaces

We study in this section the best approximation of $e^{A}v$ in the resolvent Krylov subspace $\mathcal{K}_{n}((\gamma-A)^{-1},v)$ dependent on the smoothness of the vector $v$.

Theorem 11**.**

Let $A$ be the generator of a bounded semigroup with bound $N$. Then we have

[TABLE]

where $C(\gamma,q,N)$ does not depend on $n$.

Proof.

With the smoothing operator $\tilde{v}_{n}:=S_{n,q}v$, $n\geq 2q$, from Lemma 10, it follows

[TABLE]

and thus

[TABLE]

It remains to bound the second term on the right. Since $\tilde{v}_{n}\in\mathcal{K}_{\lfloor\frac{n}{2}\rfloor}((\gamma-A)^{-1},(\gamma-A)^{-q}v)$, we have for arbitrary coefficients $b_{0},\ldots,b_{q-1}$, $a_{1},\ldots,a_{\lfloor\frac{n}{2}\rfloor-q}\in\mathbb{C}$ that

[TABLE]

In the following, we use an approximation result given in [13] for the so-called $\varphi$-functions

[TABLE]

that belong to $\widetilde{\mathcal{M}}$ and fit in the functional calculus introduced above. We leave $a_{1},\ldots,a_{\lfloor\frac{n}{2}\rfloor-q}$ arbitrary and choose $b_{j}=1/j!$, $j=0,\ldots,q-1$, according to the sum in (12) to obtain

[TABLE]

Now we have for $n>2q+1$ the estimate

[TABLE]

where we used Theorem 4.2 in [13] for the second inequality. Due to

[TABLE]

we find for $n>2q+1$ that

[TABLE]

Since this bound holds true up to finitely many numbers, adjusting the constant renders the bound true for all $n\in\mathbb{N}$ and our theorem is proved. ∎

Remark 6.1**.**

A discussion on a suitable choice of the shift $\gamma$ in the resolvents of the rational approximations to the $\varphi$-functions can be found in [11].

7 Smoothness-detection by the resolvent Krylov subspace method

Let $A$ be a linear operator on a Hilbert space $H$ with $\mathrm{Range}(\lambda-A)=H$ for some $\lambda$ with $\mathrm{Re}\,\lambda>0$ and

[TABLE]

By the Lumer-Phillips theorem for Hilbert spaces (e.g. Corollary 4.3.11 on page 186 in [26]), $A$ is the generator of a contraction $C_{0}$-semigroup with

[TABLE]

i.e. a bounded semigroup with $N=1$. We designate by $P_{n}$ the orthogonal projection onto $\mathcal{K}_{n}((\gamma-A)^{-1},v)$ and by $L_{n}=P_{n}AP_{n}$ the restriction of $A$ to this subspace. For simplicity, we assume in this section that $v\in\mathcal{D}(A)$, and therefore $\mathcal{K}_{n}((\gamma-A)^{-1},v)\subseteq\mathcal{D}(A)$. Then $L_{n}$ also satisfies

[TABLE]

as well as $\mathrm{Range}(\lambda-L_{n})=H$ for some $\lambda$ with $\mathrm{Re}\,\lambda>0$, and therefore generates a contraction semigroup. Hence, our functional calculus can be applied to functions of $A$ as well as $L_{n}$.

The next theorem gives an upper bound for the approximation of the operator $\varphi$-function of $A$ times $v\in\mathcal{D}(A^{q})$ in the resolvent Krylov subspace $\mathcal{K}_{n}((\gamma-A)^{-1},v)$. More precisely, we approximate $\varphi_{j}(A)v$ by $\varphi_{j}(L_{n})v$, where $\varphi_{j}$ is for $j=1,2,\ldots$ the $\varphi$-function defined above in (12) and $\varphi_{0}$ denotes the exponential function.

Theorem 12**.**

Let $L_{n}=P_{n}AP_{n}$, where $P_{n}$ is the orthogonal projection onto the resolvent Krylov subspace $\mathcal{K}_{n}((\gamma-A)^{-1},v)$. For $v\in\mathcal{D}(A^{q})$ and $j\geq 0$, we have

[TABLE]

where $C(\gamma,q,j)$ does not depend on $n$.

Proof.

Let $v\in\mathcal{D}(A^{q})$ and $k\geq 0$ be arbitrary. To bound $\|\varphi_{j}(A)v-\varphi_{j}(L_{n})v\|$, we use as a first step our smoothing operator $\tilde{v}_{n}:=S_{n,q}v$, $n>2q+1$ w.l.o.g., and turn this difference into

[TABLE]

In the following, both terms will be bounded separately. We start with the second term. Since

[TABLE]

and since the resolvent Krylov subspace approximation is exact for every rational function $r(z)=p_{n-1}(z)/(\gamma-z)^{n-1}$, $p_{n-1}(z)\in\mathcal{P}_{n-1}$ (see [3]), we have for every function $r_{1}(z)$ of the form

[TABLE]

that

[TABLE]

Choosing $r_{1}$ as the best approximation to $\varphi_{j}$ in the sense of our functional calculus, we can estimate

[TABLE]

where we used Theorem 4.2 in [13] for the second inequality. Note that for the case $j=0$ the terms $\|e^{A}-r_{1}(A)\|$ and $\|e^{L_{n}}-r_{1}(L_{n})\|$ are just bounded by generic constants $C(\gamma,q)$. With Lemma 10, we finally obtain

[TABLE]

It remains to bound the first term on the right-hand side of (13). For this purpose, we first remark that we have

[TABLE]

due to the special choice of our smoothing operator $\tilde{v}_{n}$. We define a second function $r_{2}(z)$ by

[TABLE]

Then the exactness property of the resolvent Krylov subspace approximation yields for an arbitrary choice of the coefficients $b_{k}$ and $a_{k}$ that

[TABLE]

and therefore

[TABLE]

We now choose $b_{0},\ldots,b_{q-1}$ according to the Taylor series of $\varphi_{j}$ and obtain with

[TABLE]

and Lemma 3 that

[TABLE]

Hence, we have with our functional calculus

[TABLE]

Exactly the same calculation holds true for $L_{n}$ and we conclude altogether

[TABLE]

for all choices of the coefficients $a_{k}$. By the exactness property, we have $L_{n}^{q}\tilde{v}_{n}=A^{q}\tilde{v}_{n}$ and hence

[TABLE]

By using Theorem 4.2. in [13], we thus obtain

[TABLE]

Since

[TABLE]

we finally have

[TABLE]

With (14) and (15), our statement is proved. ∎

How $\varphi_{k}(L_{n})v$ can be represented with the help of quasi-matrices when an orthonormal basis of the approximation space $\mathcal{K}_{n}((\gamma-A)^{-1},v)$ is known, is described in [11].

In view of abstract evolution equations, we are usually interested in the approximation of $e^{\tau A}v$ or $\varphi_{j}(\tau A)v$ for $v\in\mathcal{D}(A^{q})$ and $\tau>0$. In this case, all the above results remain valid, we only have to replace $A$ by $\tau A$ anywhere.

Moreover, all presented results for the semigroup $e^{\tau A}$ applied to some initial data $v~{}\in~{}\mathcal{D}(A^{q})$ transfer to the discrete case, that is, to the approximation of $e^{\tau A_{N}}\Psi_{0}$, where $A_{N}$ is the discretization matrix of the differential operator $A$ and $\Psi_{0}$ is the discretized initial value $v$. Since the error bounds do not depend on $\|A\|$, we obtain a convergence rate that is independent of the spatial grid.

8 Numerical experiments

We illustrate our theoretical findings by a finite-difference discretization in Subsection 8.1 as well as by a finite-element discretization of a wave equation in Subsection 8.2. Besides these illustrations, our theory provides an explanation for the behavior observed in several applications. For example, the grid-independent convergence of the rational Krylov method for the solution of Maxwell’s Equations in photonic crystal modeling in [4] is explained by our analysis.

8.1 Finite-difference discretization of the wave equation

We consider the standard wave equation on the unit square $\Omega=(0,1)^{2}$ with homogeneous Dirichlet boundary conditions written as a system of first-order

[TABLE]

The operator $A$ on the Hilbert space $\mathcal{D}(\sqrt{-\Delta})\times L^{2}(\Omega)$ equipped with the inner product

[TABLE]

where $v=[v_{1},v_{2}]^{T}$ and $w=[w_{1},w_{2}]^{T}$, satisfies the properties of Section 7 and therefore generates a contraction semigroup (or more exactly, a $C_{0}$-group). The domain of the operator $A$ reads $\mathcal{D}(A)=\mathcal{D}(\Delta)\times\mathcal{D}(\sqrt{-\Delta})$ and, with standard Sobolev spaces, $\mathcal{D}(A)=(H^{2}(\Omega)\cap H_{0}^{1}(\Omega))\times H_{0}^{1}(\Omega)$. By $\|\cdot\|$, we denote the norm induced by the inner product defined in (17). The $q$-dependent initial values $y_{0}^{q}$ with $u_{0}=u^{\prime}_{0}=g_{0}^{q}/\|[g_{0}^{q},g_{0}^{q}]^{T}\|$, where

[TABLE]

satisfy $y_{0}^{q}\in\mathcal{D}(A^{q})$ and $y_{0}^{q}\not\in\mathcal{D}(A^{q+1})$.

In order to illustrate and verify our theory numerically, we discretize the operator $A$ via finite differences on the grid $(ih,jh)$, $i,j=1,\ldots,d$ with $h=\frac{1}{d+1}$, which leads to a system of the type above with the matrix

[TABLE]

where $\otimes$ is the Kronecker product. The matrix $\Delta_{N}\in\mathbb{R}^{N\times N}$ is the standard discretization with the five-point stencil for the Laplacian. We deal with the space $\mathbb{R}^{2N}$ equipped with the inner product

[TABLE]

where $\Psi=[\Psi_{1},\Psi_{2}]^{T}$ and $\Theta=[\Theta_{1},\Theta_{2}]^{T}$, $\Psi_{i},\Theta_{i}\in\mathbb{R}^{N}$, $i=1,2$, and $(\cdot,\cdot)_{2}$ designates the standard Euclidean inner product in $\mathbb{R}^{N}$. The matrix $A_{N}$ also satisfies the assumptions in Section 7. We define discretizations of the initial values $y_{0}^{q}$ by

[TABLE]

For these initial values, we have

[TABLE]

where $\|\cdot\|_{h}$ is the norm induced by the inner product (18) and $C$ is a generic constant that does not depend on $N$. Therefore, the error of the resolvent Krylov subspace approximation to the matrix exponential measured in the discrete norm is bounded independently of $N$ according to our Theorem 12 as

[TABLE]

with $L_{n}=P_{n}A_{N}P_{n}$, where $P_{n}$ is the orthogonal projection onto $\mathcal{K}_{n}((\gamma-\tau A_{N})^{-1},\Psi^{q})$. Note that the right-hand side with the continuous values, does not depend on $N$. This worst case sublinear convergence can be clearly observed in the numerical results in Figure 2 with $\tau=0.5$, where the error in the discrete norm is shown over the dimension of the Krylov subspace for discretizations of $A$ with $d=15,31,127,255,511,1023,2047$, leading to matrices $A_{N}$ from size $450\times 450$ to size $8,380,418\times 8,380,418$. For the smaller matrices with $d=15$ and $d=31$, the convergence is faster than predicted. For the remaining matrices up to size $8,380,418\times 8,380,418$, the predicted sublinear convergence can be seen. Furthermore, the convergence is faster for the smoother initial value $\Psi^{4}$ on the right-hand side of Figure 2 compared to $\Psi^{2}$ on the left-hand side, which also fits perfectly to our theorem. For a suitable space discretization, this behavior is always to be expected. In Figure 3, the error of the backward Euler method and the resolvent Krylov subspace method, respectively, is shown versus the computing time in minutes for the discretization with dimension $8,380,418$. The resolvents have been computed by a multigrid method and the exact solution for the computation of the error has been calculated by a discrete fast Fourier transform.

8.2 Finite-element discretization of a wave equation on a non-standard domain

For a trapezoidal domain $\Omega$ with a slit in the upper half, we consider the two-dimensional wave equation

[TABLE]

with homogeneous Neumann boundary conditions which can be represented by the first-oder system

[TABLE]

where $\Delta$ is the Laplacian including the boundary conditions. It can be shown that $A$ is the generator of a contraction semigroup with respect to the inner product analogous to (17), where the operator $\sqrt{-\Delta}$ is replaced by $\sqrt{I-\Delta}$. The domain of $A$ is given by $\mathcal{D}(A)=\mathcal{D}(\Delta)\times\mathcal{D}(\sqrt{I-\Delta})$. We use here the equivalent norms that are commonly used for finite elements. For $v=[v_{1},v_{2}]^{T}$ and $w=[w_{1},w_{2}]^{T}$ with $v_{1},w_{1}\in H^{1}(\Omega)$, $\nabla_{n}v_{1}=\nabla_{n}w_{1}=0$ on $\partial\Omega$ and $v_{2},w_{2}\in L^{2}(\Omega)$, the inner product reads

[TABLE]

We solve equation (19) numerically by using finite elements with $N$ linear nodal basis functions $\phi_{k}\in H^{1}(\Omega)$, $k=1,\ldots,N$. This leads to the semi-discrete formulation

[TABLE]

where the vectors $\Psi_{1}(t),\Psi_{2}(t)\in\mathbb{R}^{N}$ are the coordinate vectors for $u(t)$ and $u^{\prime}(t)$. The mass matrix $M\in\mathbb{R}^{N\times N}$ and the stiffness matrix $S\in\mathbb{R}^{N\times N}$ are defined by $(M)_{jk}=(\phi_{j},\phi_{k})_{L^{2}(\Omega)}$ and $(S)_{jk}=-(\nabla\phi_{j},\nabla\phi_{k})_{L^{2}(\Omega)}$ for $j,k=1,\ldots,N$. Multiplying (20) from the left with the inverse of the block diagonal matrix $\text{diag}(M,M)$, we end up with the initial value problem

[TABLE]

with $\Psi(t)=[\Psi_{1}(t),\Psi_{2}(t)]^{T}$ and $\Psi_{0}=[\Psi_{1,0},\Psi_{2,0}]^{T}$. For the initial data $y_{0}^{q}$ in (19), we choose $u_{0}=u^{\prime}_{0}=g_{0}^{q}/\|[g_{0}^{q},g_{0}^{q}]^{T}\|$ with

[TABLE]

where $y_{0}^{q}\in\mathcal{D}(A^{q})$ but $y_{0}^{q}\not\in\mathcal{D}(A^{q+1})$. For the initial data used in our numerical experiment, we evolved $y_{0}^{q}$ with the given wave equation for time $\tau=0.5$. Its discrete counterpart, depicted in Figure 4 on the left-hand side, we denote by $\Psi_{0}^{q}$. Moreover, we show on the right-hand side the numerical solution for time $\tau=3$. In Figure 5, the obtained error curves are plotted for a coarse grid with 31,232 triangles and 160,323 nodes as well as a fine grid with 499,712 triangles and 251,520 nodes. As parameters, we have chosen the time step size $\tau=0.05$ and the smoothness indices $q=1$ and $q=3$.

In Figure 5, the obtained error curves for the approximation of $e^{\tau A_{N}}\Psi_{0}^{q}$ in the Krylov subspace $\mathcal{K}_{n}((15-\tau A_{N})^{-1},\Psi_{0}^{q})$ are plotted for the coarse grid as well as for the a fine grid. As parameters, we have chosen the time step size $\tau=0.05$ and the smoothness indices $q=1$ and $q=3$. The linear systems with the matrix $15-\tau A_{N}$ were solved again by a multigrid method.

9 Conclusion

In this work, we could show that the resolvent Krylov subspace method is suitable for the approximation of a large set of operator functions. For the semigroup and related operator functions, convergence rates dependent on the smoothness of the initial data have been presented. In contrast to standard methods, the faster convergence for smoother initial data is automatic, that is, the method does not need to be altered in any way to achieve the faster convergence for smoother initial data. The theoretical findings have been illustrated by numerical experiments.

Acknowledgements

This work has been supported by the Deutsche Forschungsgemeinschaft (DFG) via GR 3787/1-1.