Runge--Kutta convolution coercivity and its use for time-dependent boundary integral equations

TL;DR

This paper investigates the inheritance of convolution coercivity by Runge--Kutta convolution quadrature methods, demonstrating that certain methods preserve this property and applying it to analyze stability and convergence in nonlinear boundary integral equations.

Contribution

It proves that all algebraically stable Runge--Kutta methods, with a shift in the Laplace domain, inherit convolution coercivity, extending previous results limited to multistep methods.

Findings

Radau IIA method of order three inherits coercivity without restrictions.

All algebraically stable Runge--Kutta methods inherit coercivity with a shifted Laplace variable.

Numerical experiments confirm the theoretical error and stability analysis.

Abstract

A coercivity property of temporal convolution operators is an essential tool in the analysis of time-dependent boundary integral equations and their space and time discretisations. It is known that this coercivity property is inherited by convolution quadrature time discretisation based on A-stable multistep methods, which are of order at most two. Here we study the question as to which Runge--Kutta-based convolution quadrature methods inherit the convolution coercivity property. It is shown that this holds without any restriction for the third-order Radau IIA method, and on permitting a shift in the Laplace domain variable, this holds for all algebraically stable Runge--Kutta methods and hence for methods of arbitrary order. As an illustration, the discrete convolution coercivity is used to analyse the stability and convergence properties of the time discretisation of a non-linear boundary integral equation that originates from a non-linear scattering problem for the linear wave equation. Numerical experiments illustrate the error behaviour of the Runge--Kutta convolution quadrature time discretisation.

Figures

Equations

$$\Vert L(s){\Vert }_{V^{\prime }\leftarrow V}\le M(1+|s|{)}^{\mu },\Vert R(s){\Vert }_{V\leftarrow V}\le M(1+|s|{)}^{\mu },\text{Re}s>\sigma .$$

$${\ell }_k(t)=\frac{1}{2\pi \mathrm{i}}{\int }_{{\sigma }^{\prime }+i\mathbb{R}}{e}^{st}{L}_k(s)ds,t\ge 0({\sigma }^{\prime }>\sigma )$$

$$u(t)=L(\partial_{t})f\,(t)=(\ell*f)(t)=\Bigl{(}\frac{d}{dt}\Bigr{)}^{k}\int_{0}^{t}\ell_{k}(t-\tau)\,f(\tau)\,d\tau,\quad\ t>0,$$

$$\begin{array}{rl}{Y}_{ni}& =y_n+\tau \sum _{j=1}^m{a}_{ij}f(t_n+c_{j}h,Y_{nj}),i=1,\dots ,m,\\ y_{n+1}& =y_n+\tau \sum _{j=1}^m{b}_{j}f(t_n+c_{j}h,Y_{nj}),\end{array}$$

$$\mathcal{A}=(a_{ij}{)}_{i,j=1}^m,b=(b_1,\dots ,b_m{)}^T,\text{mathbbm}1=(1,1,\dots ,1{)}^T.$$

$$\Delta(\zeta)=\Bigl{(}\mathscr{A}+\frac{\zeta}{1-\zeta}\mathbbm{1}b^{T}\Bigr{)}^{-1},\qquad\zeta\in\mathbb{C}\hbox{ with }|\zeta|<1.$$

$$\Delta (\zeta )={\mathcal{A}}^{-1}-\frac{\zeta }{1-R(\infty )\zeta }{\mathcal{A}}^{-1}\text{mathbbm}1b^T{\mathcal{A}}^{-1}.$$

$$L\Bigl{(}\frac{\Delta(\zeta)}{{\tau}}\Bigr{)}=\sum_{n=0}^{\infty}\boldsymbol{W}_{n}(L)\zeta^{n},$$

$${}{}\bigl{(}L(\boldsymbol{\partial}_{t}^{\tau})\boldsymbol{f}\bigr{)}_{n}=\sum_{j=0}^{n}\boldsymbol{W}_{n-j}(L)\boldsymbol{f}_{j}$$

$$\bigl{(}L(\partial_{t})f\bigr{)}(t_{n+1})\approx e_{m}^{T}\bigl{(}L(\boldsymbol{\partial}_{t}^{\tau})\boldsymbol{f}\bigr{)}_{n},$$

$$L_2({\partial }_t^{\tau })L_1({\partial }_t^{\tau })\mathbf{f}=(L_2{L}_1)({\partial }_t^{\tau })\mathbf{f}.$$

$$(u,v)=\sum _{i=1}^m{b}_i{\overline{u}}_i{v}_i,u,v\in {\mathbb{C}}^m.$$

$$\text{Re}(w,\Delta (\zeta )w)\ge 0for\; allw\in {\mathbb{C}}^m,|\zeta |<1.$$

$$\text{Re}(v,\Delta (\zeta {)}^{-1}v)\ge 0for\; allv\in {\mathbb{C}}^m,|\zeta |<1.$$

$$\Delta (\zeta {)}^{-1}=\mathcal{A}+\frac{\zeta }{1-\zeta }\text{mathbbm}1b^T=\mathcal{C}+\frac{1}{2}\frac{1+\zeta }{1-\zeta }\text{mathbbm}1b^T,with\mathcal{C}=\mathcal{A}-\frac{1}{2}\text{mathbbm}1b^T,$$

$$(\text{mathbbm}1b^{T}v,v)$$

$$2\text{Re}(\mathcal{C}v,v)$$

$$\mathcal{A}=\left(\begin{array}{cc}5/12& -1/12\\ 3/4& 1/4\end{array}\right),b^T=(3/4,1/4).$$

$$\text{Re}(w,\Delta (\zeta )w)\ge \frac{1}{2}(1-|\zeta {|}^2)|w{|}^{2}for\; allw\in {\mathbb{C}}^m,|\zeta |\le 1.$$

$$\text{Re}(w,\Delta (\zeta )w)\ge c\delta |w{|}^{2}for\; allw\in {\mathbb{C}}^m,|\zeta |\le e^{-\delta }.$$

$$\Delta (\zeta )=\frac{1}{2}\left(\begin{array}{cc}3& 1-4\zeta \\ -9& 5+4\zeta \end{array}\right).$$

$$\text{Re}(w,\Delta (\zeta )w)={\overline{w}}^T{\mathcal{B}}^{1/2}\cdot {\mathcal{B}}^{-1/2}\frac{1}{2}(\mathcal{B}\Delta (\zeta )+\Delta (\overline{\zeta }{)}^T\mathcal{B}){\mathcal{B}}^{-1/2}\cdot {\mathcal{B}}^{1/2}w.$$

$${\mathcal{B}}^{-1/2}\frac{1}{2}(\mathcal{B}\Delta (\zeta )+\Delta (\overline{\zeta }{)}^T\mathcal{B}){\mathcal{B}}^{-1/2}=\frac{1}{2}\left(\begin{array}{cc}3& -\sqrt{3}(1+2\zeta )\\ -\sqrt{3}(1+2\overline{\zeta })& 5+4\text{Re}\zeta \end{array}\right),$$

$$W^T\mathcal{B}W=I_m$$

$$\mathcal{A}=WX{W}^{-1},$$

$$\text{Re}(w,\Delta (\zeta )w)=\text{Re}{\overline{w}}^T\mathcal{B}\Delta (\zeta )w$$

$$=\text{Re}{\overline{w}}^T{\mathcal{B}}^{1/2}\cdot {\mathcal{B}}^{1/2}W\cdot W^T{\mathcal{B}}^{1/2}\cdot {\mathcal{B}}^{1/2}\Delta (\zeta ){\mathcal{B}}^{-1/2}\cdot {\mathcal{B}}^{1/2}W\cdot W^T{\mathcal{B}}^{1/2}\cdot {\mathcal{B}}^{1/2}w$$

$$=\text{Re}{\overline{w}}^T{\mathcal{B}}^{1/2}\cdot {\mathcal{B}}^{1/2}W\cdot W^T\mathcal{B}\Delta (\zeta )W\cdot W^T{\mathcal{B}}^{1/2}\cdot {\mathcal{B}}^{1/2}w,$$

$$W^T\mathcal{B}\Delta (\zeta )W=W^{-1}\Delta (\zeta )W=(W^{-1}\Delta (\zeta {)}^{-1}W{)}^{-1}.$$

$$W^T\mathcal{B}\Delta (\zeta )W=Y+Z(\zeta ),$$

$$\text{Re}(w,\Delta (\zeta )w)=0$$

$$\text{Re}⟨v,L(s)v⟩\ge \alpha \Vert R(s)v{\Vert }^2,for\; allv\in V,\text{Re}s\ge \sigma .$$

Full text

††journal: BIT

1 Introduction

This paper is concerned with a discrete coercivity property that ensures the stability of time discretisations of boundary integral equations for wave equations, also in situations such as

• non-linear boundary integral equations;

• boundary integral equations coupled with a wave equation in an interior domain, with an explicit time discretisation in the domain.

For convolution quadrature based on A-stable multistep methods (which have approximation order at most two), it is known from BanLS that the coercivity property is preserved under time discretisation, uniformly in the temporal stepsize. Here we study the preservation of convolution coercivity under time discretisation by Runge–Kutta convolution quadrature. Up to a shift in the Laplace variable and a corresponding reformulation of the boundary integral equation for an exponentially scaled solution function, we show that the convolution coercivity property is preserved by all convolution quadratures based on algebraically stable Runge–Kutta methods, which include in particular Radau IIA methods of arbitrary order. Without any such shift and exponential scaling, the convolution coercivity is shown to be preserved by the two-stage Radau IIA method of order three.

We illustrate the use of the discrete convolution coercivity by the stability and convergence analysis of the Runge–Kutta convolution quadrature time discretisation of a non-linear boundary integral equation for a non-linear scattering problem for the acoustic wave equation. This problem has been studied with different numerical methods in BanR .

The discrete convolution coercivity is not needed for the corresponding linear scattering problem, because there the convolution quadrature time discretisation of the linear boundary integral equation can be interpreted as a convolution quadrature discretisation of the convolution operator that maps the data to the solution. Therefore known bounds of the Laplace transform of the solution operator and known error bounds of convolution quadrature yield stability and error bounds Lub94 ; Say16 . The same argument can also be used for the coupling of a linear wave equation in an interior domain with the boundary integral equation that describes transparent boundary conditions, provided that the convolution quadrature for the boundary integral equation is based on the same (implicit) time discretisation method as for the wave equation in the interior domain. This precludes explicit time-stepping in the interior. For the coupling of convolution quadrature on the boundary with an explicit time discretisation in the interior, the discrete convolution coercivity as considered in the present paper is needed; see BanLS ; Ebe ; KovL for the coupling of implicit BDF2 convolution quadrature on the boundary with explicit leapfrog time-stepping in the domain for acoustic, elastic and electro-magnetic wave equations, respectively.

The paper is organised as follows:

In Section 2 we recall the continuous-time convolution coercivity, which is related to a coercivity property of the Laplace transform of the (distributional) convolution kernel that holds uniformly for all values of the Laplace-domain frequency variable in a (possibly shifted) right half-plane.

In Section 3 we study the preservation of the convolution coercivity under time discretisation by Runge–Kutta convolution quadrature. This preservation depends on the numerical range of the Runge–Kutta differentiation symbol, which is shown to lie in the right half-plane for algebraically stable Runge–Kutta methods. With a matrix-function inequality that is obtained as an extension of a theorem of von Neumann, we then prove our main result, Theorem 3.1, which yields the discrete convolution coercivity.

Section 4 recapitulates error bounds of Runge–Kutta convolution quadrature shown in BanLM .

In Section 5 we apply our results to the time discretisation of the wave equation with a non-linear impedance boundary condition. We study only semi-discretisation in time, but note that this could be extended to full discretisation with the techniques of BanR . The error behaviour is illustrated by numerical experiments in Section 6. In the numerical experiments it is observed that the convolution quadrature based on the three-stage Radau IIA method performs well even without the shift and exponential scaling, which is more favourable than our theoretical results.

2 Coercivity of temporal convolutions

The following coercivity result is given in BanLS , where it is used as a basic result in studying boundary integral operators for the acoustic wave equation; see also KovL for Maxwell’s equation and Ebe for elastic wave equations. The result can be viewed as a time-continuous operator-valued extension of a theorem of Herglotz from 1911, which states that an analytic function has positive real part on the unit disk if and only if convolution with its coefficient sequence is a positive semi-definite operation.

Let $V$ be a complex Hilbert space and $V^{\prime}$ its dual, and let $\langle\cdot,\cdot\rangle$ denote the anti-duality between $V$ and $V^{\prime}$. Let $L(s):V\rightarrow V^{\prime}$ and $R(s):V\rightarrow V$ be analytic families of bounded linear operators for $\text{\rm Re}\,s>\sigma$, continuous for $\text{\rm Re}\,s\geq\sigma$. We assume the uniform bounds, with some real exponent $\mu$,

[TABLE]

This polynomial bound guarantees that $L$ is the Laplace transform of a distribution $\ell$. If we write $L(s)=s^{k}L_{k}(s)$ with an integer $k>\mu+1$, then the Laplace inversion formula

[TABLE]

defines a continuous and exponentially bounded function $\ell_{k}$, which has $\ell$ as its $k$th distributional derivative. We write the convolution with $\ell$ as

[TABLE]

for functions $f$ on $[0,T]$ whose extension to $t<0$ by [math] is $k$ times continuously differentiable. Similarly we consider the convolution $R(\partial_{t})f$.

Theorem 2.1

(BanLS*, *, Lemma 2.2)** Let $\alpha\geq 0$. In the above situation, the following statements are equivalent:

$\text{\rm Re}\,\langle v,L(s)v\rangle\geq\alpha\|R(s)v\|^{2}\qquad\hbox{for all }v\in V,\;\text{\rm Re}\,s>\sigma$*. * 2. 2.

$\displaystyle\int_{0}^{\infty}e^{-2\sigma t}\,\text{\rm Re}\,\langle f(t),L(\partial_{t})f(t)\rangle\,dt\geq\alpha\int_{0}^{\infty}e^{-2\sigma t}\,\|R(\partial_{t})f(t)\|^{2}\,dt$

for all $f\in C^{k}([0,\infty),V)$ with finite support and $f^{(j)}(0)=0$ for $0\leq j<k$, and for all $t\geq 0$.

Property 1. is known to be satisfied for the Laplace transforms of various boundary integral operators for wave equations BamH ; BanLS ; KovL ; LalS ; Say16 , and it is a fundamental property in the study of boundary integral equations for wave equations.

We are interested in time discretisations of the convolution operators $L(\partial_{t})$ and $R(\partial_{t})$ that preserve this coercivity property. It was shown in BanLS that this is achieved by convolution quadrature based on A-stable multistep methods such as the first- and second-order backward differentiation formulae. In Theorem 3.1 below we will show that the coercivity property is also preserved by convolution quadrature based on certain Runge–Kutta methods such as the third-order, two-stage Radau IIA method. For the particular case $\sigma=0$, it will be shown to be preserved for all algebraically stable Runge–Kutta methods.

3 Preserving coercivity by Runge–Kutta convolution quadrature

3.1 Runge–Kutta differentiation symbol and convolution quadrature

An $m$-stage Runge–Kutta discretisation of the initial value problem $y^{\prime}=f(t,y)$, $y(0)=y_{0}$, is given by

[TABLE]

where ${\tau}>0$ is the time step, $t_{n}=n{\tau}$, and the internal stages $Y_{ni}$ and grid values $y_{n}$ are approximations to $y(t_{n}+c_{i}{\tau})$ and $y(t_{n})$, respectively. In the following we use the notation

[TABLE]

We always assume that the Runge–Kutta matrix $\mathscr{A}$ is invertible.

As has been shown in LubO ; SchLL ; BanL ; BanLM , and in applications to wave propagation problems further in BanK ; BanMS ; BanR ; WanW , Runge–Kutta methods can be used to construct convolution quadrature methods that enjoy favourable properties. Here one uses the Runge–Kutta differentiation symbol

[TABLE]

This is well-defined for $|\zeta|<1$ if $R(\infty)=1-b^{T}\mathscr{A}^{-1}\mathbbm{1}$ satisfies $|R(\infty)|\leq 1$. In fact, the Sherman-Morrison formula then yields

[TABLE]

To formulate the Runge–Kutta convolution quadrature for $L(\partial_{t})g$, we formally replace in $L(s)$ the differentiation symbol $s$ by $\Delta(\zeta)/{\tau}$ and expand the operator-valued matrix function

[TABLE]

where in the case of $L(s):V\to V^{\prime}$ we have the convolution quadrature matrices $\boldsymbol{W}_{n}(L):V^{m}\to(V^{\prime})^{m}$. For the discrete convolution with these matrices we use the notation

[TABLE]

for any sequence $\boldsymbol{f}=(\boldsymbol{f}_{n})$ in $V^{m}$. For vectors of function values of a function $f:[0,T]\to V$ given as \boldsymbol{f}_{n}=\bigl{(}f(t_{n}+c_{i}{\tau})\bigr{)}_{i=1}^{m}, the $i$th component of the vector \bigl{(}L(\boldsymbol{\partial}_{t}^{\tau})\boldsymbol{f}\bigr{)}_{n} is considered as an approximation to \bigl{(}L(\partial_{t})f\bigr{)}(t_{n}+c_{i}{\tau}).

In particular, if $c_{m}=1$, as is the case with Radau IIA methods, then the continuous convolution at $t_{n+1}$ is approximated by the last component of the discrete block convolution:

[TABLE]

where $e_{m}=(0,\dots,0,1)^{T}$ is the $m$th unit vector. We recall the composition rule

[TABLE]

For $\lambda\in\mathbb{C}$, the convolution quadrature $(\boldsymbol{\partial}_{t}^{\tau}-\lambda)^{-1}\boldsymbol{f}$ (which is to be interpreted as $L(\boldsymbol{\partial}_{t}^{\tau})\boldsymbol{f}$ for the multiplication operator $L(s)=(s-\lambda)^{-1}$) contains the internal stages of the Runge–Kutta approximation to the linear differential equation $y^{\prime}=\lambda y+f$ with initial value $y(0)=0$.

Results on the order of convergence of this approximation are given in BanL ; BanLM ; LubO . The result of BanLM , which is relevant for operators $L(s)$ arising in wave propagation, will be restated and extended to the internal stages in Section 4.

3.2 Numerical range of the Runge–Kutta differentiation symbol

We now consider methods that are algebraically stable:

• •

All weights $b_{i}$ are positive.

• •

The symmetric matrix with entries $b_{i}a_{ij}+b_{j}a_{ji}-b_{i}b_{j}$ is positive semi-definite.

Gauss methods and Radau IIA methods are widely used classes of methods that satisfy this condition. We refer the reader to HaiW for background literature on Runge–Kutta methods and their stability notions.

We consider the weighted inner product on $\mathbb{C}^{m}$,

[TABLE]

We have the following characterisation.

Lemma 1

For an algebraically stable Runge–Kutta method and for the $b$-weighted inner product (3),

[TABLE]

Conversely, if the differentiation symbol of a Runge–Kutta method with positive weights $b_{i}$ satisfies this inequality, then the method is algebraically stable.

Proof

With a different notation, this is shown in (LubO87, , p. 232, (6.19)). For the convenience of the reader we include the short proof. Since for $v=\Delta(\zeta)w$ we have $(w,\Delta(\zeta)w)=(\Delta(\zeta)^{-1}v,v)$, it suffices to show that

[TABLE]

We rewrite

[TABLE]

and observe that (cf.HaiL )

[TABLE]

Since $\text{\rm Re}\,(1+\zeta)/(1-\zeta)>0$ for $|\zeta|<1$, the result follows. ∎

In this paper we will need a stronger positivity property, for which we show the following order barrier and a positive result for the two-stage Radau IIA method, which is of order 3 and has the coefficients

[TABLE]

Lemma 2

(a) For the two-stage Radau IIA method and for the $b$-weighted inner product (3) and corresponding norm $|\cdot|$ we have

[TABLE]

(b) For none of the Gauss methods with two or more stages and none of the Radau IIA methods with three or more stages, there exists $c>0$ such that for all sufficiently small $\delta>0$,

[TABLE]

Clearly, (4) implies (5) with $c$ arbitrarily close to 1 for small $\delta$. We further note that the implicit Euler method and the implicit midpoint rule (which are the one-stage Radau IIA and Gauss methods, respectively) also satisfy (5).

Proof

(a) For the two-stage Radau IIA method we find

[TABLE]

Denoting the diagonal matrix of the weights by $\mathscr{B}={\rm diag}(3/4,1/4)$, we note

[TABLE]

We obtain the hermitian matrix

[TABLE]

which has the trace $4+2\text{\rm Re}\,\zeta$ and the determinant $3(1-|\zeta|^{2})$. It follows that both eigenvalues are positive and bounded by $6$, and hence the smaller eigenvalue is bounded from below by $3(1-|\zeta|^{2})/6=(1-|\zeta|^{2})/2$. This yields the inequality (4).

(b) The proof uses the W-transformation of Hairer & Wanner, see (HaiW, , p. 77). For each of the $m$-stage Gauss and Radau IIA methods, there exists an invertible real $m\times m$ matrix $W$ with first column $\mathbbm{1}$ such that, with the diagonal matrix $\mathscr{B}$ of the weights $b_{i}$,

[TABLE]

or in other words, $W^{T}\mathscr{B}^{1/2}$ is an orthogonal matrix (with respect to the Euclidean inner product), and

[TABLE]

where $X-\tfrac{1}{2}e_{1}e_{1}^{T}-\beta_{m}e_{m}e_{m}^{T}$ is a skew-symmetric matrix with $\beta_{m}=0$ for the Gauss method and $\beta_{m}>0$ for the Radau IIA method. We write

[TABLE]

where we note

[TABLE]

Now, by the definition of $\Delta(\zeta)$ and the above-mentioned property of $W^{-1}\mathscr{A}W=X$ together with $W^{T}b=e_{1}$, the matrix $W^{-1}\Delta(\zeta)^{-1}W$ is the sum of a skew-hermitian matrix plus a rank-1 or rank-2 matrix for Gauss or Radau IIA methods, respectively, and by the Sherman-Morrison-Woodbury formula so is its inverse:

[TABLE]

where $Y$ is skew-hermitian and $Z(\zeta)$ is of rank 1 or 2 for Gauss or Radau IIA methods, respectively. If $w\neq 0$ is in the null-space of $Z(\zeta)W^{T}\mathscr{B}$, which is of codimension 1 or 2 for Gauss or Radau, respectively, then we obtain from the above formulas that

[TABLE]

in contradiction to (5). ∎

As we will show in Theorem 3.1 below, Runge–Kutta convolution quadrature with (5) preserves the coercivity property of Theorem 2.1 for arbitrary abscissa $\sigma\geq 0$, while general algebraically stable methods preserve it in the case $\sigma=0$. Before we state and prove this theorem in Section 3.4, we need an auxiliary result of independent interest.

3.3 A matrix-function inequality related to a theorem by von Neumann

We consider again a complex Hilbert space $V$ and its dual $V^{\prime}$, with the anti-duality denoted by $\langle\cdot,\cdot\rangle$. On $\mathbb{C}^{m}$ we consider an inner product $(\cdot,\cdot)$ and associated norm $|\cdot|$. An inner product on $V^{m}$ and the anti-duality between $V^{m}$ and $(V^{\prime})^{m}$ are induced in the usual way: for Kronecker products $a\otimes u$ and $b\otimes v$ with $a,b\in\mathbb{C}^{m}$ and $u,v\in V$ one defines $(a\otimes u,b\otimes v)=(a,b)\,(u,v)$ and extends this to a sequilinear form on $V^{m}\times V^{m}$, and in the same way one proceeds for the anti-duality $\langle\cdot,\cdot\rangle$ on $V^{m}\times(V^{\prime})^{m}$.

Lemma 3

On the Hilbert space $V$, let $L(s):V\rightarrow V^{\prime}$ and $R(s):V\rightarrow V$ be analytic families of bounded linear operators for $\text{\rm Re}\,s>\sigma$, continuous for $\text{\rm Re}\,s\geq\sigma$, such that (1) is satisfied and for some $\alpha\geq 0$,

[TABLE]

Let the matrix $S\in\mathbb{C}^{m\times m}$ be such that

[TABLE]

Then,

[TABLE]

This result can be viewed as an extension of a theorem of von Neumann Neu (see also (HaiW, , p. 179)), which corresponds to the particular case where $L(s)$ is the identity operator on $V$ (when $V$ is identified with $V^{\prime}$ with the anti-duality given by the inner product on $V$).

Proof

The proof adapts Michel Crouzeix’s proof of von Neumann’s theorem as given in (HaiW, , p. 179 f.). Without loss of generality we assume here $\sigma=0$.

First we note that for a diagonal matrix $S$ the result holds trivially, and so it does for a normal matrix $S$, which is diagonalised by a similarity transformation with a unitary matrix.

For a non-normal matrix $S$ we consider the matrix-valued complex function

[TABLE]

and we observe that $S=S(1)$ and

[TABLE]

Together with the condition on $S$ this shows that the numerical range of $S(z)$ is in the right complex half-plane for $\text{\rm Re}\,z\geq 0$, and hence all eigenvalues of $S(z)$ have non-negative real part. Therefore, the operator functions $L(S(z))$ and $R(S(z))$ are well-defined for $\text{\rm Re}\,z\geq 0$.

If $\text{\rm Re}\,z=0$, then the matrix $S(z)$ is normal, and hence the desired inequality is valid for $S(z)$ with $\text{\rm Re}\,z=0$. The function

[TABLE]

is subharmonic, since the last term is harmonic as the real part of an analytic function and the first term is the inner product of an analytic function with itself, which is subharmonic (as is readily seen by computing the Laplacian and noting that the real and imaginary parts of the analytic function are harmonic). Hence, by the maximum principle (or its Phragmén-Lindelöf-type extension to polynomially bounded subharmonic functions on the half-plane),

[TABLE]

which is the desired inequality. ∎

Remark 1

There is a slightly weaker variant of Lemma 3. We formulate it for $\sigma=0$. Let $L(s):V\rightarrow V^{\prime}$ and $R(s):V\rightarrow V$ be analytic families of bounded linear operators for $\text{\rm Re}\,s>0$, continuous for $s=0$, such that (1) is satisfied and for some $\alpha\geq 0$,

[TABLE]

Let the matrix $S\in\mathbb{C}^{m\times m}$ be such that all its eigenvalues either have positive real part or are zero, and

[TABLE]

Then,

[TABLE]

This is proved by continuity, using the previous result for $S+\epsilon I$ and letting $\epsilon\to 0$.

3.4 Preserving the convolution coercivity under discretisation

Theorem 3.1

Let the $m$-stage Runge–Kutta method satisfy (5) for some inner product $(\cdot,\cdot)$, as in particular is the case for the two-stage Radau IIA method. In the situation of Theorem 2.1, condition 1. of that theorem implies, for sufficiently small stepsize ${\tau}>0$ and with $\widetilde{\sigma}=\sigma/c$,

[TABLE]

for every sequence $\boldsymbol{f}=(\boldsymbol{f}_{n})_{n\geq 0}$ in $V^{m}$ with finitely many non-zero entries. Moreover, in the case $\sigma=0$ this inequality holds for every algebraically stable Runge–Kutta method, with $\widetilde{\sigma}=0$ and with respect to the $b$-weighted inner product (3) on $\mathbb{C}^{m}$.

Proof

The proof uses Parseval’s formula and combines Lemma 3 with Lemmas 1 and 2. By (5), with $\widetilde{\sigma}=\sigma/c$ and $\rho=e^{-\widetilde{\sigma}{\tau}}$ we have with respect to the inner product weighted by the $b_{i}$ that

[TABLE]

We abbreviate

[TABLE]

and similarly $\widehat{\boldsymbol{R}}(\theta)$. We denote the Fourier series

[TABLE]

By Parseval’s formula and the definition of the convolution quadrature weights $\boldsymbol{W}_{n}(L)$,

[TABLE]

Here (6) used in Lemma 3 yields

[TABLE]

Moreover, again by Parseval’s formula,

[TABLE]

which yields the result. ∎

4 Error bounds of Runge–Kutta convolution quadrature

In this section we restate the result of BanLM and extend it to cover the approximation properties of the internal stages, which will be needed in the next section. To avoid restating the list of properties required for the underlying Runge–Kutta method, we state the results just for the Radau IIA methods, which appear to be the practically most important class of Runge–Kutta methods to be used for convolution quadrature.

Let $K(s)$, for $\text{\rm Re}\,s>\sigma>0$, be an analytic family of operators between Hilbert spaces $V$ and $W$ (or just Banach spaces are sufficient here), such that for some real exponent $\mu$ and $\nu\geq 0$ the operator norm is bounded as follows:

[TABLE]

Theorem 4.1

(BanLM*, **, Theorem 3)**

Let $K$ satisfy (7) and consider the Runge-Kutta convolution quadrature based on the Radau IIA method with $m$ stages. Let $r>\max(2m-1+\mu,2m-1,m+1)$ and $f\in C^{r}([0,T],V)$ satisfy $f(0)=f^{\prime}(0)=\ldots=f^{(r-1)}(0)=0$. Then, there exists $\bar{{\tau}}>0$ such that for $0<{\tau}\leq\bar{{\tau}}$ and $t_{n}=n{\tau}\in[0,T]$,*

[TABLE]

The constant $C$ is independent of ${\tau}$ and $f$, but does depend on $\bar{{\tau}}$, $T$, and the constants in (7).

The proof in BanLM is readily extended to yield the following error bound for the internal stages. Note that here the full order $2m-1$ is replaced by the stage order plus one, $m+1$. We give the result for $m\geq 2$ stages, so that $m+1\leq 2m-1$. (For $m=1$, the implicit Euler method, one can use the previous result.)

Theorem 4.2

Let $K$ satisfy (7) and consider the Runge-Kutta convolution quadrature based on the Radau IIA method with $m\geq 2$ stages. Let $r>\max(m+1+\mu,m+1)$ and $f\in C^{r}([0,T],V)$ satisfy $f(0)=f^{\prime}(0)=\ldots=f^{(r-1)}(0)=0$. Then, there exists $\bar{{\tau}}>0$ such that for $0<{\tau}\leq\bar{{\tau}}$ and $t_{n}=n{\tau}\in[0,T]$,

[TABLE]

The constant $C$ is independent of ${\tau}$ and $f$, but does depend on $\bar{{\tau}}$, $T$, and the constants in (7).

5 Application to the time discretisation of the wave equation with a non-linear impedance boundary condition

5.1 A non-linear scattering problem

We consider the wave equation on an exterior smooth domain $\Omega^{+}\subset\mathbb{R}^{3}$. Following BanR , we search for a function $u(\cdot,t)\in H^{1}(\Omega^{+})$ satisfying the weak form of the wave equation

[TABLE]

with zero initial conditions and with the non-linear boundary condition

[TABLE]

where $\partial_{\nu}^{+}$ is the outer normal derivative on the boundary $\Gamma$ of $\Omega^{+}$, where $g:\mathbb{R}\to\mathbb{R}$ is a given monotonically increasing function, and where $u^{inc}(x,t)$ is a given solution of the wave equation. The interpretation is that the total wave $u^{tot}=u+u^{inc}$ is composed of the incident wave $u^{inc}$ and the unknown scattered wave $u$.

One approach to solve this exterior problem is to determine the Dirichlet and boundary data from boundary integral equations on $\Gamma$ and then to compute the solution at points of interest $x\in\Omega^{+}$ from the Kirchhoff representation formula. Here we are interested in the stability and convergence properties of the numerical approximation when the time discretisation in the boundary integral equation and in the representation formula is done by Runge–Kutta convolution quadrature. Since our interest in this paper is the aspect of time discretisation, we will not address the space discretisation by boundary elements, though with some effort this could also be included; cf. BanR .

5.2 Boundary integral equation and representation formula

Using the standard notation of the boundary integral operators for the Helmholtz equation $s^{2}u-\Delta u=0$ $(\text{\rm Re}\,s>0)$ as used, for example, in LalS ; Say16 and BanLS ; BanR , we denote by

[TABLE]

the single-layer and double-layer potential operators, respectively, and by $V(s)$, $K(s)$, $K^{T}(s)$, $W(s)$ the corresponding boundary integral operators that form the Calderón operator

[TABLE]

and the related operator

[TABLE]

where the suffix imp stands for impedance. With these operators, the solution $u$ is determined by first solving, for

[TABLE]

(where $\gamma^{+}$ is the trace operator on $\Omega^{+}$), the time-dependent boundary integral equation (see BanR )

[TABLE]

The solution $u$ is then obtained from the representation formula

[TABLE]

We will address the question as to what are the approximation properties when the temporal convolutions in (12) and (13) are discretised by Runge–Kutta convolution quadrature. Since this will not turn out fully satisfactory, we will further consider time-differentiated versions of (12).

The following coercivity property was proved in BanLS .

Lemma 4

(BanLS*, *, Lemma 3.1)** Let $\langle\cdot,\cdot\rangle_{\Gamma}$ denote the anti-duality pairing between $H^{-1/2}(\Gamma)\times H^{1/2}(\Gamma)$ and $H^{1/2}(\Gamma)\times H^{-1/2}(\Gamma)$. There exists $\beta>0$ so that the Calderón operator (10) satisfies

[TABLE]

for $\text{\rm Re}\,s\geq\sigma>0$ and for all $\varphi\in H^{-1/2}(\Gamma)$ and $\psi\in H^{1/2}(\Gamma)$, with $c_{\sigma}=\min(1,\sigma^{2})\,\sigma$.

Since $B_{\mathrm{imp}}(s)$ differs from $B(s)$ by a skew-hermitian matrix, the same estimate then also holds for $B_{\mathrm{imp}}(s)$. Note that Lemma 4 yields property 1. of Theorem 2.1 for the Calderón operator $B(s)$ and for the multiplication operator $R(s)=s^{-1}$.

5.3 Time discretisation by Runge–Kutta convolution quadrature

Using the notation of Section 3.1, the boundary integral equation (12) is discretised in time with a stepsize ${\tau}>0$ over a time interval $(0,T)$ with $T=N{\tau}$ by

[TABLE]

where $(\boldsymbol{\varphi}^{\tau},\boldsymbol{\psi}^{\tau})=(\boldsymbol{\varphi}_{n},\boldsymbol{\psi}_{n})_{n=0}^{N-1}$ with $(\boldsymbol{\varphi}_{n},\boldsymbol{\psi}_{n})=(\varphi_{n,i},\psi_{n,i})_{i=1}^{m}$ and $(\varphi_{n,i},\psi_{n,i})\approx(\varphi(t_{n}+c_{i}{\tau}),\psi(t_{n}+c_{i}{\tau}))$ is the numerical approximation that is to be computed, and $\dot{\boldsymbol{u}}^{inc}=({\dot{\boldsymbol{u}}}^{inc}_{n})_{n=0}^{N-1}$ with ${\dot{\boldsymbol{u}}}^{inc}_{n}=(\dot{u}^{inc}(t_{n}+c_{i}{\tau}))_{i=1}^{m}$. The function $g$ acts componentwise. At the $n$th time step, a non-linear system of equations of the following form needs to be solved:

[TABLE]

where the dots represent known terms. This has a unique solution, because the eigenvalues of $\Delta(0)=\mathscr{A}^{-1}$ have positive real part, and Lemma 4 and the monotonicity of $g$ then yield the unique existence of the solution by the Browder–Minty theorem; cf. BanR for the analogous situation for multistep-based convolution quadrature.

As an alternative to (14), we further consider the time discretisation of the time-differentiated boundary integral equation:

[TABLE]

which is now solved for the approximations $(\dot{\boldsymbol{\varphi}}^{\tau},\dot{\boldsymbol{\psi}}^{\tau})$ (where the dot is just suggestive notation) to $(\dot{\varphi},\dot{\psi})$ (where the dot means again time derivative). Here we define $\boldsymbol{\psi}^{\tau}=(\boldsymbol{\partial}_{t}^{\tau})^{-1}\dot{\boldsymbol{\psi}}^{\tau}$ and the same for $\boldsymbol{\varphi}^{\tau}$. Furthermore, $\ddot{\boldsymbol{u}}^{inc}$ contains the values $\ddot{u}^{inc}(t_{n}+c_{i}\tau)$.

We can go even further and consider the time discretisation of the twice differentiated boundary integral equation:

[TABLE]

where again the dots on the approximation $(\ddot{\boldsymbol{\varphi}}^{\tau},\ddot{\boldsymbol{\psi}}^{\tau})$ are suggestive notation, and we set $\dot{\boldsymbol{\psi}}^{\tau}=(\boldsymbol{\partial}_{t}^{\tau})^{-1}\ddot{\boldsymbol{\psi}}^{\tau}$ and $\boldsymbol{\psi}^{\tau}=(\boldsymbol{\partial}_{t}^{\tau})^{-2}\ddot{\boldsymbol{\psi}}^{\tau}$, and the same for $\boldsymbol{\varphi}^{\tau}$.

Finally, at any point $x\in\Omega^{+}$ of interest we compute the approximation to the solution value $u(x,t_{n}+c_{i}\tau)$ by using the same Runge–Kutta convolution quadrature for discretizing the representation formula (13):

[TABLE]

5.4 Error bounds for the linear case

We consider first the case of a linear impedance function

[TABLE]

Let ${\boldsymbol{u}}^{\tau}=(\boldsymbol{u}_{n})_{n=0}^{N-1}$ with $\boldsymbol{u}_{n}=(u_{n,i})_{i=1}^{m}$ be the solution approximation obtained by the discretised representation formula (17) with either of the discretised boundary integral equations (14) or (15) or (16). The discretisation is done by Runge–Kutta convolution quadrature based on the Radau IIA method with $m$ stages.

Here we obtain the following optimal-order pointwise error bounds for $x$ bounded away from $\Gamma$.

Proposition 1

Suppose that in a neighbourhood of the boundary $\Gamma$, the incident wave $u_{inc}$ together with its extension by [math] to $t<0$ is sufficiently regular. For $x\in\Omega^{+}$ with $\operatorname{dist}(x,\Gamma)\geq\delta>0$, the following optimal-order error bound is satisfied in the linear situation described above: for $0\leq t_{n}=n\tau\leq T$,

[TABLE]

Proof

We denote

[TABLE]

By Lemma 4, $B_{\alpha}(s)$ is invertible for $\alpha\geq 0$ with the bound, for $\text{\rm Re}\,s\geq\sigma>0$,

[TABLE]

The exact solution $u(x,t)$ is given by the representation formula (13) with

[TABLE]

For $x\in\Omega^{+}$ we define the operators $S_{x}(s):H^{-1/2}(\Gamma)\to\mathbb{C}$ and $D_{x}(s):H^{1/2}(\Gamma)\to\mathbb{C}$ by

[TABLE]

These operators are bounded for $\text{\rm Re}\,s\geq\sigma>0$ and dist$(x,\Gamma)\geq\delta>0$ by

[TABLE]

The first bound is proved in (BanLM, , Lemma 6) and the second bound is proved similarly.

We thus have

[TABLE]

with

[TABLE]

With the above operator bounds we obtain for $\text{\rm Re}\,s\geq\sigma>0$ and dist$(x,\Gamma)\geq\delta>0$

[TABLE]

By the composition rule, the numerical solution obtained by (14) and (17) is given as

[TABLE]

where $\boldsymbol{f}$ contains the values of $f$ at the points $t_{n}+c_{i}{\tau}$. If we take instead (15) or (16), then we have

[TABLE]

respectively. In view of (21), Theorem 4.1 then yields the result. ∎

The situation is different if we consider the $H^{1}(\Omega^{+})$ norm of the error.

Proposition 2

Suppose that in a neighbourhood of the boundary $\Gamma$, the incident wave $u_{inc}$ together with its extension by [math] to $t<0$ is sufficiently regular. Then, the following error bounds are satisfied in the linear situation described above: for $0\leq t_{n}=n\tau\leq T$,

[TABLE]

with

[TABLE]

Proof

Consider the Laplace transformed wave equation (8)

[TABLE]

where $\hat{u}$ is the Laplace transform of $u$ and $\hat{f}$ the Laplace transform of

[TABLE]

We will require the estimate, see (BamH, , Equation (2.9)),

[TABLE]

with $\text{\rm Re}\,s\geq\sigma>0$ and the scaled $H^{1}$ norm

[TABLE]

Testing (22) with $s\hat{u}$, integrating by parts and taking the real part gives

[TABLE]

where $\hat{\psi}=s\gamma^{+}\hat{u}$ is the Laplace transform of $\psi$. Making use of $\|\hat{u}\|_{H^{1}(\Omega^{+})}\leq C(\sigma)\|\hat{u}\|_{|s|,\Omega^{+}}$ and the bound (18) gives

[TABLE]

The stated result then follows from Theorem 4.1. ∎

5.5 Convergence for the non-linear problem

There are several aspects which make the error analysis of the non-linear problem more intricate:

• •

The numerical solution can no longer be interpreted as a mere convolution quadrature for an appropriate operator $K(s)$ acting on the data (i.e., the incident wave).

• •

We need to impose regularity assumptions on the solution rather than the data.

• •

Convolution coercivity now plays an important role in ensuring the stability of the time discretisation.

We assume strict monotonicity of the non-linear function $g:\mathbb{R}\to\mathbb{R}$: there exists $\beta>0$ such that

[TABLE]

Furthermore, we assume that the pointwise application of $g$ maps $H^{1/2}(\Gamma)$ to $H^{-1/2}(\Gamma)$. As is shown in BanR by Sobolev embeddings, this is satisfied if $g(\xi)$ grows at most cubically as $|\xi|\to\infty$.

In the following we write for a stepsize $\tau>0$ and a sequence $\boldsymbol{e}=(\boldsymbol{e}_{n})_{n=0}^{N-1}$ with $\boldsymbol{e}_{n}=(e_{n,i})_{i=1}^{m}$ and $e_{n,i}$ in a Hilbert space $V$

[TABLE]

We denote the numerical solution by ${\boldsymbol{u}}^{\tau}=(u_{n,i})$ and the corresponding values of the exact solution by $\boldsymbol{u}=(u(t_{n}+c_{i}\tau))$, where in both cases $n=0,\dots,N-1$ and $i=1,\dots,m$.

We have the following error bound for the non-linear problem. Here the restriction to the two-stage Radau IIA method stems from Lemma 2.

Proposition 3

Let the non-linear function $g$ be continuous, strictly monotone and have at most cubic growth. Suppose that the solution $u$ to the problem (8)–(9) is sufficiently regular.

Consider the time discretisation (14) and (17) by the two-stage Radau IIA convolution quadrature method. Then, there is $\bar{\tau}>0$ such that for stepsizes $0<\tau\leq\bar{\tau}$, the error in the boundary values satisfies the bound

[TABLE]

and the error in the exterior domain is bounded by

[TABLE]

The constants $C$ are independent of $\tau$ and $N$ with $0<N\tau\leq T$, but depend on $T$.

Proof

We eliminate $\varphi$ in the system of boundary integral equations (12) to arrive at a boundary integral equation for $\psi$,

[TABLE]

where

[TABLE]

with the exterior Dirichlet-to-Neumann operator $\mathrm{DtN}^{+}(s)$. It follows from Propositions 17 and 18 (and their proofs) in LalS that, for $\text{\rm Re}\,s\geq\sigma>0$, there exist $C(\sigma)$ and $\alpha(\sigma)>0$ such that

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the anti-duality pairing between $H^{-1/2}(\Gamma)$ and $H^{1/2}(\Gamma)$.

Thanks to the composition rule, we can do the same for the numerical discretisation (14) and reduce the numerical system to an equation for $\boldsymbol{\psi}^{\tau}$, which is just the convolution quadrature time discretisation of (27),

[TABLE]

For the error $\boldsymbol{\epsilon}=\boldsymbol{\psi}^{\tau}-\boldsymbol{\psi}$ with \boldsymbol{\psi}=\bigl{(}(\psi(t_{n}+c_{i}{\tau}))_{i=1}^{m}\bigr{)}_{n=0}^{N-1} we then have the error equation

[TABLE]

with the defect

[TABLE]

which is the convolution quadrature error for $L(\partial_{t})\psi$. By Theorem 4.2 and our assumption of a sufficiently regular $\psi=\gamma^{+}\dot{u}$, this is bounded by

[TABLE]

Since we can apply the same argument also to spatial derivatives of $\psi$ (in the assumed case of a smooth boundary $\Gamma$), we even have

[TABLE]

We test (31) with $\boldsymbol{\epsilon}$, multiply with $e^{-2\widetilde{\sigma}t}$ with $\widetilde{\sigma}=1/T$ and integrate from [math] to $T$. With (29) and the Runge-Kutta convolution coercivity as given by Theorem 3.1, and with the strict monotonicity (31) we conclude that

[TABLE]

and estimate further

[TABLE]

We thus find the stability estimate

[TABLE]

Since $(\boldsymbol{\partial}_{t}^{{\tau}})^{-1}\boldsymbol{\epsilon}=\gamma^{+}{\boldsymbol{u}}^{\tau}-\gamma^{+}\boldsymbol{u}$, this proves (25).

Let us denote by $M(s)=S(s)V^{-1}(s):H^{1/2}(\Gamma)\to H^{1}(\Omega^{+})$ the operator that maps Dirichlet data in $H^{1/2}(\Gamma)$ to the corresponding solution $\hat{u}\in H^{1}(\Omega^{+})$ of the Helmholtz equation $s^{2}\hat{u}-\Delta\hat{u}=0$. By (Say16, , Equation (3.10)), this is bounded for $\text{\rm Re}\,s\geq\sigma>0$ by

[TABLE]

We then have

[TABLE]

By Theorem 4.2 and the bound for $M$, the last term is bounded by $O(\tau^{5/2})$ in the $H^{1}(\Omega^{+})$ norm. The first term is only $O(\tau^{3/2})$, since we lose a factor $\tau^{3/2}$ from the $O(\tau^{3})$ error bound for $\gamma^{+}u$ because of the $O(|s|^{3/2})$ bound of $M(s)$; this follows from Lemma 5.2 in BanL and Parseval’s identity. ∎

In a similar way we obtain the following results for the alternative discretisations (15) and (16):

(i) In addition to Proposition 3, assume that $g$ has bounded second derivatives. With the discretisation (15) instead of (14), the error bound in the $H^{1}(\Omega^{+})$ norm improves to $O(\tau^{5/2})$, and the $\ell_{2}^{\tau}$ error in a point $x$ bounded away from the boundary $\Gamma$ is at most $O(\tau)$.

(ii) In addition to Proposition 3, assume that $g$ has bounded second and third derivatives. With the discretisation (16) instead of (14), the error bound in the $H^{1}(\Omega^{+})$ norm improves to $O(\tau^{3})$, and the $\ell_{2}^{\tau}$ error in a point $x$ bounded away from the boundary $\Gamma$ is at most $O(\tau^{2})$.

The proofs of these error bounds are very similar to that of Proposition 3, using in addition a discrete Gronwall inequality at the end of the estimation of $\boldsymbol{\epsilon}$, and an $O(|s|^{3})$ bound for the norm of the operator from $H^{1/2}(\Gamma)\to\mathbb{C}$ that maps Dirichlet data to the solution of the Helmholtz equation $s^{2}\hat{u}-\Delta\hat{u}=0$ at a point $x\in\Omega^{+}$ bounded away from $\Gamma$, for $s$ in a right half-plane. Since our main concern here is to illustrate the use of the convolution coercivity, we omit the details of these extensions.

Remark 2

If we set $\widetilde{L}(s)=L(s+\sigma)$ and $\widetilde{\psi}(t)=e^{-\sigma t}\psi(t)$ for some $\sigma>0$, then the boundary integral equation (27) is equivalent to

[TABLE]

By (29), we then have the coercivity estimate for $\widetilde{L}(s)$ for all $\text{\rm Re}\,s\geq 0$ (and not just for $\text{\rm Re}\,s\geq\sigma$):

[TABLE]

By Theorem 3.1, the coercivity estimate for the convolution quadrature approximation of $\widetilde{L}(\partial_{t})\widetilde{\psi}$ is then obtained for every algebraically stable Runge-Kutta method (and not just the two-stage Radau IIA method). Hence, by discretising the shifted boundary integral equation (33) on an interval $[0,T]$ with shift $\sigma=1/T$, we obtain Runge–Kutta based convolution quadrature time discretisations of arbitrarily high order of convergence (assuming sufficient regularity of the exact solution). We remark that similar shifts are familiar in the convergence analysis of space-time Galerkin methods for time-dependent boundary integral equations BamH . As in that case, numerical experiments indicate that implementing the shift may not be necessary in practical computations, although this is not backed by theory.

6 Numerical experiments

6.1 Scattering by the unit sphere

In these experiments we let $\Omega^{+}$ be the exterior of the unit sphere and the trace of the incident wave $u^{inc}$ on the sphere be space independent. As constant functions are eigenfunctions of all the integral operators on the sphere Ned , the solution will also be constant in space. The eigenvalue for the combined operator $L(s)$ in (27) is given by

[TABLE]

for any $\hat{\psi}$ constant in space. This operator will reflect well the behaviour of scattering by a convex obstacle, but not that of a general scatterer. For this reason we concentrate on the corresponding interior problem with

[TABLE]

again for $\hat{\psi}$ constant. Treating both these operators as scalar, complex valued functions of $s$, we see that both have a better behaviour than the general operators, see (28) and (29). Namely

[TABLE]

and

[TABLE]

As the operator $L(s)$ is too simple, in the numerical experiments we only consider the scalar, non-linear equation

[TABLE]

Even though these operators are of such a simple form, due to the nonlinearity the exact solution is not available. Nevertheless, a highly accurate solution is not expensive to evaluate and can be used to compute the error in the $\ell^{\tau}_{2}$ norm. We have performed the numerical experiments with the following choices of $g$ and $u^{inc}$

[TABLE]

and with final time $T=6$. Note that $g_{1}$ is once continuously differentiable whereas $g_{2}$ is infinitely differentiable. The data $u^{inc}$ is not causal, but it is vanishingly small for $t<0$ and we have found that this discrepancy has no significant effect on the results.

In Figure 1 we show the convergence of the two-stage Radau IIA convolution quadrature. As expected, for the smooth non-linear condition we obtain full order of convergence. The solution and its first derivative are shown in Figure 2. Note that the two solutions have a similar shape, but a closer look at the derivative in Figure 3 reveals that one is smooth and the other only once continuously differentiable.

For the interior problem, as $\text{\rm Re}\,L(s)\geq 1$ the theory also applies to higher order Radau IIA methods. This is however not the case with $L^{-}(s)$. We nevertheless perform experiments with the three-stage Radau IIA method and obtain good results as shown in Figure 4.

6.2 A full non-scalar example

We end the paper with a 2D example that requires the full BEM discretisation in space. The domain is an L-shape and the incident wave is a plane wave. Piecewise linear boundary element space is used to approximate the Dirichlet trace $\psi$ and piecewise constant boundary element space to approximate the Neumann trace $\varphi$ and the time-discretisation is performed using the two-stage Radau IIA method. The images of the solution are shown in Figure 5.

Acknowledgement

We thank Ernst Hairer for helpful discussions. This work was partially supported by DFG, SFB 1173.