Runge--Kutta convolution coercivity and its use for time-dependent boundary integral equations
Lehel Banjai Christian Lubich



Runge–Kutta convolution coercivity



Lehel Banjai The Maxwell Institute for Mathematics in the Sciences; School of Mathematical & Computer Sciences, Heriot-Watt University, EH14 4AS Edinburgh, UK

L.Banjai@hw.ac.uk 
Christian Lubich Mathematisches Institut, Universität Tübingen, Auf der Morgenstelle,
  D-72076 Tübingen, Germany.

  Lubich@na.uni-tuebingen.de

L. Banjai and Ch. Lubich



Runge–Kutta convolution coercivity and its use for time-dependent boundary integral equations
    Lehel Banjai Christian Lubich
    
=============================================================================================

 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.



 65L05 65R20




§ 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 <cit.> 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 <cit.>.

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 <cit.>. 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 <cit.> 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 <ref>, which yields the discrete convolution coercivity.

Section 4 recapitulates error bounds of Runge–Kutta convolution quadrature shown in <cit.>.

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 <cit.>. 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.













§ COERCIVITY OF TEMPORAL CONVOLUTIONS


The following coercivity result is given in <cit.>, where it is used as a basic result in studying boundary integral operators for the acoustic wave equation; see also <cit.> for Maxwell's equation and <cit.> 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' its dual, and let ⟨·,·⟩ denote the anti-duality between V and V'.
Let L(s): V → V' and R(s) : V → V be analytic families of bounded linear operators for s > σ, continuous for s≥σ.  We assume the uniform bounds, with some real exponent μ,

    L(s)_V' ← V≤ M(1+ |s|)^μ,   R(s)_V ← V≤ M (1+|s|)^μ,      s > σ.











This polynomial bound guarantees that L is the Laplace transform of a distribution ℓ.
If we write L(s)=s^k L_k(s) with an integer k>μ+1, then the Laplace inversion formula

    ℓ_k(t) = 1/2π∫_σ'+i e^st  L_k(s)  ds,    t≥ 0    (σ'>σ)

defines a continuous and exponentially bounded function ℓ_k, which has ℓ as its kth distributional derivative. We write the convolution with ℓ as

    u(t) = L(∂_t)f (t) = (ℓ*f)(t)=(d/dt)^k ∫_0^t ℓ_k(t-τ)  f(τ)  dτ,   t>0,

for functions f on [0,T] whose extension to t<0 by 0 is k times continuously differentiable.
Similarly we consider the convolution R(∂_t) f.

 <cit.> Let α≥ 0.
  In the above situation, the following statements are equivalent:
  
  
  * ⟨ v, L(s) v⟩≥αR(s) v^2     for all  v ∈ V,   s > σ.
  


  * ∫_0^∞ e^-2σ t ⟨ f(t) , L(∂_t) f (t) ⟩  dt ≥α∫_0^∞ e^-2σ t R(∂_t) f(t)^2  dt
 for all  f ∈ C^k([0,∞), V) with finite support and f^(j)(0)=0 for 0≤ j<k, 

and for all t ≥ 0.
  


Property 1. is known to be satisfied for  the Laplace transforms of various boundary integral operators for wave equations <cit.>, 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(∂_t) and R(∂_t) that preserve this coercivity property. It was shown in <cit.> 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 <ref> 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 σ=0, it will be shown to be preserved for all algebraically stable Runge–Kutta methods.



§ PRESERVING COERCIVITY BY RUNGE–KUTTA CONVOLUTION QUADRATURE





 §.§ Runge–Kutta differentiation symbol and convolution quadrature


An m-stage Runge–Kutta  discretisation of the initial value problem y'
= f(t,y), y(0) = y_0, is given by

    Y_ni   = y_n + ∑_j = 1^m a_ij f(t_n+c_jh,Y_nj),      i =
    1,…,m,
    
    y_n+1    = y_n + ∑_j = 1^m b_j f(t_n+c_jh,Y_nj),

where >0 is the time step, t_n = n, and the internal stages Y_ni and grid values y_n are
approximations to y(t_n+c_i ) and y(t_n), respectively. 
In the  following we use the notation

    𝒜 = (a_ij)_i,j = 1^m,    b = (b_1,…,b_m)^T,
    
       = (1,1,…,1)^T.

We always assume that the Runge–Kutta matrix 𝒜 is invertible.

As has been shown in <cit.>, and in applications to wave propagation problems further  in <cit.>, Runge–Kutta methods can be used to construct convolution quadrature methods that enjoy favourable properties.
Here one uses the Runge–Kutta differentiation symbol

    Δ(ζ) = (𝒜+ζ/1-ζ b^T)^-1,     ζ∈ with  |ζ|<1.

This is well-defined for |ζ|<1 if 
R(∞)=1-b^T𝒜^-1 satisfies |R(∞)|≤ 1. In fact, the Sherman-Morrison formula then yields

    Δ(ζ) = 𝒜^-1 -ζ/1-R(∞)ζ𝒜^-1 b^T 𝒜^-1.

To formulate the Runge–Kutta convolution quadrature for L(∂_t )g, we formally replace in L(s) the differentiation symbol s by Δ(ζ)/ and expand the operator-valued matrix function

    L(Δ(ζ)/) = ∑_n=0^∞W_n(L) ζ^n,

where in the case of L(s):V→ V' we have the convolution quadrature matrices W_n(L):V^m → (V')^m. 
For the discrete convolution with these matrices we use the notation

    (L(∂_t^)f)_n = ∑_j=0^n W_n-j(L) f_j

for any sequence f=(f_n) in V^m. For vectors of function values of a function f:[0,T]→ V given as f_n=(f(t_n+c_i))_i=1^m, the ith component of the vector (L(∂_t^)f)_n is considered as an approximation to
(L(∂_t)f)(t_n+c_i).

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:

    (L(∂_t)f)(t_n+1) ≈ e_m^T (L(∂_t^)f)_n,

where e_m = (0, …, 0, 1)^T is the mth unit vector.
We recall the composition rule 

    L_2(∂_t^) L_1(∂_t^)f=(L_2L_1)(∂_t^)f.

For λ∈, the convolution quadrature (∂_t^-λ)^-1f (which is to be interpreted as L(∂_t^)f for the multiplication operator L(s)=(s-λ)^-1) contains the internal stages of the Runge–Kutta approximation to the linear differential equation y'=λ y +f with initial value y(0)=0. 

Results on the order of convergence of this approximation are given in <cit.>. The result of <cit.>, which is relevant for operators L(s) arising in wave propagation, will be restated and extended to the internal stages in Section <ref>.





 §.§ 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_ia_ij+b_ja_ji- b_ib_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 <cit.> for background literature on Runge–Kutta methods and their stability notions.

We consider the weighted inner product on ^m,

    (u,v) = ∑_i=1^m b_i u_i v_i,      u,v∈^m.

We have the following characterisation.

 
For an algebraically stable Runge–Kutta method and for the b-weighted inner product (<ref>),

    (w,Δ(ζ)w) ≥ 0    for all w∈^m, |ζ|<1.

Conversely, if the differentiation symbol of a Runge–Kutta method with positive weights b_i satisfies this inequality, then the method is algebraically stable.


 With a different notation, this is shown in <cit.>. For the convenience of the reader we include the short proof. 
Since for v=Δ(ζ)w we have (w,Δ(ζ)w)=(Δ(ζ)^-1v,v), it suffices to show that 

    (v,Δ(ζ)^-1v) ≥ 0    for all v∈^m, |ζ|<1.

We rewrite

    Δ(ζ)^-1 = 𝒜+ζ/1-ζ b^T = 𝒞+ 1/21+ζ/1-ζ b^T,     with  𝒞= 𝒜-12 b^T,

and observe that (cf.<cit.>)

    ( b^Tv,v)    = | ∑_i=1^m b_i v_i |^2
    
    
    2 (𝒞v,v)    = ∑_i,j=1^m (b_ia_ij+b_ja_ji-b_ib_j) v̅_i v_j.

Since (1+ζ)/(1-ζ)>0 for |ζ|<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

    𝒜 = [  5/12 -1/12;   3/4   1/4 ],     b^T = (3/4,1/4).


(a) For the two-stage Radau IIA method and for the b-weighted inner product (<ref>) and corresponding norm |·| we have

    (w,Δ(ζ)w) ≥12 (1-|ζ|^2) |w|^2    for all  w∈^m,  |ζ|≤ 1.

(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 δ>0,

    (w,Δ(ζ)w) ≥ cδ |w|^2    for all  w∈^m,  |ζ| ≤ e^-δ.



Clearly, (<ref>) implies (<ref>) with c arbitrarily close to 1 for small δ. 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  (<ref>).

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

    Δ(ζ) = 1/2[      3 1 - 4ζ;     -9 5 + 4ζ ].

Denoting the diagonal matrix of the weights by ℬ= diag(3/4,1/4), we 
note 

    (w,Δ(ζ)w) = w^T ℬ^1/2·ℬ^-1/212(ℬΔ(ζ)+Δ(ζ)^Tℬ) ℬ^-1/2·ℬ^1/2w.

We obtain the hermitian matrix

    ℬ^-1/212(ℬΔ(ζ)+Δ(ζ)^Tℬ) ℬ^-1/2 = 1/2[            3 -√(3) (1+2ζ); -√(3) (1+2ζ)      5 + 4 ζ ],

which has the trace 4+2ζ and the determinant
3(1-|ζ|^2). It follows that both eigenvalues are positive and bounded by 6, and hence the smaller eigenvalue is bounded from below by
3(1-|ζ|^2)/6=(1-|ζ|^2)/2. This yields the inequality (<ref>).

(b) The proof uses the W-transformation of Hairer & Wanner, see <cit.>. For each of the m-stage Gauss and Radau IIA methods, there exists an invertible real m× m matrix W with first column  such that, with the diagonal matrix ℬ of the weights b_i,

    W^T ℬ W=I_m

or in other words, W^T ℬ^1/2 is an orthogonal matrix (with respect to the Euclidean inner product), and

    𝒜= WXW^-1,

where X-12 e_1e_1^T - β_m e_me_m^T is a skew-symmetric matrix with β_m=0 for the Gauss method and β_m>0 for the Radau IIA method.
We write

    (w,Δ(ζ)w) = w^T ℬΔ(ζ) w 
    
          = w^T  ℬ^1/2·ℬ^1/2W · W^T ℬ^1/2·ℬ^1/2Δ(ζ) ℬ^-1/2·ℬ^1/2W · W^T ℬ^1/2·ℬ^1/2 w 
    
          = w^T  ℬ^1/2·ℬ^1/2W · 
    W^T ℬΔ(ζ) W · W^T ℬ^1/2·ℬ^1/2 w ,

where we note

    W^T ℬΔ(ζ) W = W^-1Δ(ζ) W = (W^-1Δ(ζ)^-1 W)^-1.

Now, by the definition of Δ(ζ) and the above-mentioned property of W^-1𝒜W=X together with W^Tb=e_1, the matrix W^-1Δ(ζ)^-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:

    W^T ℬΔ(ζ) W = Y+Z(ζ),

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

    (w,Δ(ζ)w) =0

in contradiction to (<ref>). 






As we will show in Theorem <ref> below, Runge–Kutta convolution quadrature with (<ref>) preserves the coercivity property of Theorem <ref> for arbitrary abscissa σ≥ 0, while general algebraically stable methods preserve it in the case σ=0.
Before we state and prove this theorem in Section <ref>, we need an  auxiliary result of independent interest.



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

We consider again a complex Hilbert space V and its dual V', with the anti-duality denoted by ⟨·, ·⟩. On ℂ^m we consider an inner product (·,·) and associated norm |·|. An inner product on V^m and the anti-duality between V^m and (V')^m are induced in the usual way: for Kronecker products a⊗ u and b⊗ v with a,b∈ℂ^m and u,v∈ V one defines
( a⊗ u, b⊗ v) =(a,b)  (u,v) and extends this to a sequilinear form on V^m× V^m, and in the same way one proceeds for the anti-duality ⟨·,·⟩ on V^m × (V')^m.

 
On the Hilbert space V,
let L(s): V → V' and R(s) : V → V be analytic families of bounded linear operators for s > σ, continuous for s ≥σ, such that (<ref>) is satisfied and for some α≥ 0,

    ⟨ v, L(s) v⟩≥αR(s) v^2,      for all  v ∈ V,   s ≥σ.

Let the matrix S∈ℂ^m× m be such that

    (w,Sw) ≥σ  |w|^2,      for all  w∈ℂ^m.

Then,

    ⟨v, L(S) v⟩≥αR(S) v^2,      for all v∈ V^m.



This result can be viewed as an extension of a theorem of von Neumann <cit.> (see also <cit.>), which corresponds to the particular case where L(s) is the identity operator on V (when V is identified with V' with the anti-duality given by the inner product on V).

 The proof adapts Michel Crouzeix's proof of von Neumann's theorem as given in <cit.>. Without loss of generality we assume here σ =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

    S(z) = z/2 (S+S^*) + 1/2 (S-S^*)

and we observe that S=S(1) and

    (w,S(z)w) = 12 ( z)   (w,Sw).

Together with the condition on S this shows that the numerical range of S(z) is in the right complex half-plane for z ≥ 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 z ≥ 0.

If z = 0, then the matrix S(z) is normal, and hence the desired inequality is valid for S(z) with z=0. The function

    φ(z) = α R(S(z)) v^2 - ⟨v, L(S(z))v⟩

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),

    φ(1) ≤sup_ z = 0φ(z) ≤ 0,

which is the desired inequality. 


 There is a slightly weaker variant of Lemma <ref>. We formulate it for σ=0. Let L(s): V → V' and R(s) : V → V be analytic families of bounded linear operators for s > 0, continuous for s=0, such that (<ref>) is satisfied and for some α≥ 0,

    ⟨ v, L(s) v⟩≥αR(s) v^2,      for all  v ∈ V,   s > 0.

Let the matrix S∈ℂ^m× m be such that all its eigenvalues either have positive real part or are zero, and

    (w,Sw) ≥ 0,      for all  w∈ℂ^m.

Then,

    ⟨v, L(S) v⟩≥αR(S) v^2,      for all v∈ V^m.

This is proved by continuity, using the previous result for S+ I and letting → 0.



































































 §.§ Preserving the convolution coercivity under discretisation
 


 Let the m-stage Runge–Kutta method satisfy (<ref>) for some inner product (·,·), as in particular is the case for the two-stage Radau IIA method.  In the situation of Theorem <ref>, condition 1. of that theorem implies, for sufficiently small stepsize >0 and with  σ= σ/c,

    ∑_n = 0^∞ e^-2σn⟨f_n, (L(∂_t^)f)_n ⟩≥α∑_n = 0^∞ e^-2σn(R(∂_t^)f)_n^2,

for every sequence f=(f_n)_n≥ 0 in V^m with finitely many non-zero entries. Moreover, in the case σ=0  this inequality holds for every algebraically stable Runge–Kutta method, with σ=0 and with respect to the b-weighted inner product (<ref>) on ^m.




The proof uses Parseval's formula and combines Lemma <ref> with Lemmas <ref> and <ref>.
By (<ref>), with σ=σ/c and ρ=e^-σ  we have with respect to the inner product weighted by the b_i that

    ( w,Δ(ρ e^iθ)/w ) ≥σ |w|^2      for all  w∈^m,  θ∈.

We abbreviate

    L(θ) = L( Δ(ρ e^iθ)/)

and similarly R(θ). We denote the Fourier series

    f(θ) = ∑_n=0^∞ρ^n e^inθf_n.

By Parseval's formula and the definition of the convolution quadrature weights W_n(L),

    ∑_n = 0^∞⟨ρ^n f_n , ∑_j = 0^n ρ^n-jW_n-j (L)ρ^j f_j⟩ = 1/2π∫_-π^π⟨f(θ), L(θ)  f(θ) ⟩ dθ.

Here (<ref>) used in Lemma <ref> yields

    ⟨f(θ), L(θ)  f(θ) ⟩≥αR(θ) f(θ)^2.

Moreover, again by Parseval's formula,

    1/2π∫_-π^πR(θ) f(θ)^2 dθ = ∑_n = 0^∞ρ^2n∑_j = 0^n W_n-j(R) f_j^2,

which yields the result. 






§ ERROR BOUNDS OF RUNGE–KUTTA CONVOLUTION QUADRATURE


In this section we restate the result of <cit.> 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  s > σ > 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 μ and ν≥ 0 the operator norm is bounded as follows:

    K(s)≤ M(σ) |s|^μ/( s)^ν   for all  s >  σ.


<cit.>
 
Let K satisfy (<ref>) and consider the Runge-Kutta 
convolution quadrature based on the Radau IIA method with m stages.
Let r > max(2m-1+μ, 2m-1, m+1) and f ∈ C^r([0,T],V) satisfy 
f(0) = f'(0) = … = f^(r-1)(0) = 0. Then, there exists > 0  such that for 0 < ≤ and t_n=n∈ [0,T],

    e_m^T (K(∂_t^)f)_n - (K(∂_t)f)(t_n+1)  
           ≤ C
     h^min(2m-1,m+1-μ+ν) (f^(r)(0)+∫_0^t f^(r+1)(τ) dτ).

The constant C is independent of  and f, but does depend on ,  T, and the constants in  (<ref>).


The proof in <cit.> 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≥ 2 stages, so that m+1≤ 2m-1. (For m=1, the implicit Euler method, one can use the previous result.)


 
Let K satisfy (<ref>) and consider the Runge-Kutta 
convolution quadrature based on the Radau IIA method with m≥ 2 stages.
Let r > max(m+1+μ,  m+1) and f ∈ C^r([0,T],V) satisfy 
f(0) = f'(0) = … = f^(r-1)(0) = 0. Then, there exists > 0  such that for 0 < ≤ and t_n=n∈ [0,T],

    (K(∂_t^)f)_n - (K(∂_t)f(t_n+c_i))_i=1^m  
           ≤ C
     h^min(m+1,m+1-μ+ν) (f^(r)(0)+∫_0^t f^(r+1)(τ) dτ).

The constant C is independent of  and f, but does depend on ,  T, and the constants in  (<ref>).






§ APPLICATION TO THE TIME DISCRETISATION OF THE WAVE EQUATION WITH A NON-LINEAR IMPEDANCE BOUNDARY CONDITION





 §.§ A non-linear scattering problem

We consider the wave equation on an exterior smooth domain Ω^+⊂^3. Following <cit.>, we search for a function u(·,t) ∈ H^1(Ω^+) satisfying the  weak form of the wave equation

    ü = Δ u      in Ω^+

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

    ∂_ν^+ u = g(u̇ + u̇^inc) - ∂_ν^+ u^inc     on Γ,

where ∂_ν^+ is the outer normal derivative on the boundary Γ of Ω^+, where g:→ 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 Γ and then to compute the solution at points of interest x∈Ω^+ 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. <cit.>.



 §.§ Boundary integral equation and representation formula


Using the standard notation of the boundary integral operators for the Helmholtz equation s^2 u-Δ u=0 ( s>0) as used, for example, in <cit.> and <cit.>, we denote by 

    S(s):H^-1/2(Γ) → H^1(^3∖Γ)    and   
    D(s):H^1/2(Γ) → H^1(^3∖Γ)
 
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

    B(s) = [     sV(s)      K(s);   -K^T(s) s^-1 W(s) ] :  H^-1/2(Γ) × H^1/2(Γ) → H^1/2(Γ) × H^-1/2(Γ)

and the related operator

    B_imp(s) = B(s) + [      0 -1/2 I;  1/2 I      0 ],

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

    φ=-∂_ν^+u  and  ψ=γ^+u̇
 
(where γ^+ is the trace operator on Ω^+), the time-dependent boundary integral equation (see <cit.>)

    B_imp(∂_t) [ φ; ψ ] +
    [             0; g(ψ + u̇^inc) ] =
    [           0; ∂_ν^+ u^inc ].

The solution u is then obtained from the representation formula

    u = S(∂_t)φ + D(∂_t)∂_t^-1ψ.

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

The following coercivity property was proved in <cit.>.

 <cit.>  Let
⟨·,·⟩_Γ denote the anti-duality pairing between H^-1/2(Γ) × H^1/2(Γ) and H^1/2(Γ) × H^-1/2(Γ).
There exists β > 0 so that the Calderón operator (<ref>) satisfies

    ⟨[ φ; ψ ], B(s) 
      [ φ; ψ ]⟩_Γ≥β c_σ (s^-1φ^2_H^-1/2(Γ) + s^-1ψ^2_H^1/2(Γ))

for  s ≥σ> 0 and for all φ∈ H^-1/2(Γ) and ψ∈ H^1/2(Γ),
with c_σ= min(1,σ^2)  σ.
 
Since B_imp(s) differs from B(s) by a skew-hermitian matrix, the same estimate then also holds for B_imp(s). Note that Lemma <ref> yields property 1. of Theorem <ref> for the Calderón operator B(s) and for the multiplication operator R(s)=s^-1.




 §.§ Time discretisation by Runge–Kutta convolution quadrature

Using the notation of Section <ref>, the boundary integral equation (<ref>) is discretised in time with a stepsize >0 over a time interval (0,T) with T=N by

    B_imp(∂_t^) [ φ^; ψ^ ] +
    [              0; g(ψ^ + u̇^inc) ] =
    [           0; ∂_ν^+ u^inc ],

where ( φ^ , ψ^)= ( φ_n , ψ_n)_n=0^N-1 with
( φ_n , ψ_n)=(φ_n,i,ψ_n,i)_i=1^m
 and (φ_n,i,ψ_n,i) ≈ (φ(t_n+c_i),ψ(t_n+c_i)) is the numerical approximation that is to be computed, and u̇^inc = (u̇^inc_n)_n=0^N-1 with u̇^inc_n=(u̇^inc(t_n+c_i))_i=1^m. The function g acts componentwise. At the nth time step, a non-linear system of equations of the following form needs to be solved:

    B_imp(Δ(0)/τ) [ φ_n; ψ_n ] +
    [                 0; g(ψ_n + u̇^inc_n) ] = …,

where the dots represent known terms. This has a unique solution, because the eigenvalues of Δ(0)=𝒜^-1 have positive real part, and Lemma <ref> and the monotonicity of g then yield the unique existence of the solution by the Browder–Minty theorem; cf. <cit.> for the analogous situation for multistep-based convolution quadrature.

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

    B_imp(∂_t^) [ φ̇^; ψ̇^ ] +
    [                            0; g'(ψ^ + u̇^inc)(ψ̇^ + ü^inc) ] =
    [            0; ∂_ν^+ u̇^inc ],

which is now solved for the approximations (φ̇^, ψ̇^ ) (where the dot is just suggestive notation) to (φ̇,ψ̇) (where the dot means again time derivative).
Here we define ψ^ = (∂_t^)^-1ψ̇^ and the same for φ^. Furthermore, ü^inc contains the values ü ^inc(t_n+c_iτ).

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

    B_imp(∂_t^) [ φ̈^; ψ̈^ ] +
    [                                                                0; g'(ψ^ + u̇^inc)(ψ̈^ + ⃛u^inc) + g”(ψ^ + u̇^inc)· (ψ̇^ + ü^inc)^2 ]
        =
    [           0; ∂_ν^+ ü^inc ],

where again the dots on the approximation (φ̈^, ψ̈^ ) are suggestive notation, and we set ψ̇^ = (∂_t^)^-1ψ̈^ and ψ^ = (∂_t^)^-2ψ̈^, and the same for φ^.

Finally, at any point x∈Ω^+ of interest we compute the approximation to the solution value u(x,t_n+c_iτ) by using the same Runge–Kutta convolution quadrature for discretizing the representation formula (<ref>):

    u^ = S(∂_t^)φ^ + D(∂_t^)(∂_t^)^-1ψ^.





























 §.§ Error bounds for the linear case

We consider first the case of a linear impedance function 

    g(ξ)=αξ   with   α≥ 0.

Let u^τ=(u_n)_n=0^N-1 with u_n=(u_n,i)_i=1^m be the solution approximation obtained by the discretised representation formula (<ref>) with either of the discretised boundary integral equations (<ref>) or  (<ref>) or  (<ref>). 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 Γ.
 Suppose that in a neighbourhood of the boundary Γ, the incident wave u_inc  together with its extension by 0 to t<0 is sufficiently regular. For x∈Ω^+ with (x,Γ)≥δ>0, the following optimal-order  error bound is satisfied in the linear situation described above: for 0≤ t_n=nτ≤ T,

    |u_n(x) - u(x,t_n) | ≤ C(δ,T) τ^2m-1.




We denote 

    B_α(s) = B_imp(s)+
    [   0   0;   0 α I ].

By Lemma <ref>, B_α(s) is invertible for α≥ 0 with the bound, for s≥σ>0,

    B_α(s)^-1≤ C(σ) |s|^2/ s.

The exact solution u(x,t) is given by the representation formula (<ref>) with

    [ φ; ψ ]
    = B_α^-1(∂_t) 
    [                    0; ∂_ν^+ u^inc -αu̇^inc ].

For x∈Ω^+ we define the operators S_x(s):H^-1/2(Γ)→ and 
D_x(s):H^1/2(Γ)→ by

    S_x(s)φ = (S(s)φ)(x)   and   D_x(s)ψ = (D(s)ψ)(x).

These operators are bounded for s≥σ >0 and dist(x,Γ)≥δ>0 by

    S_x(s)  _ H^-1/2(Γ)≤ C(σ,δ)  |s|  e^-δ  s
        D_x(s)  _ H^1/2(Γ) ≤ C(σ,δ)   |s|  e^-δ  s.

The first bound is proved in <cit.> and the second bound is proved similarly.

We thus have

    u(x,t) = (M_x(∂_t)f)(t)

with

    M_x(s)=(S_x(s),D_x(s)s^-1)B_α(s)^-1   and   
    f= [                    0; ∂_ν^+ u^inc -αu̇^inc ].

With the above operator bounds we obtain for s ≥σ >0 and dist(x,Γ)≥δ>0

    M_x(s) _ H^1/2(Γ) × H^-1/2(Γ)≤ C(σ,δ)  |s|^3/ s 
    e^-δ  s .

By the composition rule, the numerical solution obtained by (<ref>) and (<ref>) is given as

    u^(x) = M_x(∂_t^)f,

where f contains the values of f at the points t_n+c_i.
If we take instead (<ref>) or (<ref>), then we have

    u^(x) = M_x(∂_t^) (∂_t^)^-1ḟ   or   u^(x) = M_x(∂_t^) (∂_t^)^-2f̈,

respectively. In view of (<ref>), Theorem <ref> then yields the result.


The situation is different if we consider the H^1(Ω^+) norm of the error.

 Suppose that in a neighbourhood of the boundary Γ, the incident wave u_inc  together with its extension by 0 to t<0 is sufficiently regular. Then, the following   error bounds are satisfied in the linear situation described above: for 0≤ t_n=nτ≤ T,

    u_n - u(·,t_n) _H^1(Ω^+)≤ C(T) τ^k

with

    k= m+1/2    if (<ref>) is used,
       k=min(2m-1,m+3/2)    if (<ref>) is used,
       k=min(2m-1,m+5/2)    if (<ref>) is used.



 
Consider the Laplace transformed wave equation (<ref>)

    -Δû + s^2 û   = 0       in Ω^+,
    ∂_ν^+û -α s û   = f̂      on Γ^,

where  û is the Laplace transform of u and f̂ the Laplace transform of 

    f = ∂_ν^+ u^inc-αu̇^inc.

 We will require the estimate, see <cit.>, 

    ∂_ν^+ û_H^-1/2(Γ)≤ C(σ) |s|^1/2û_|s|,Ω^+,

with s ≥σ > 0 and the scaled H^1 norm

    û^2_|s|,Ω^+ = ∇û^2_L^2(Ω^+)+|s|^2 û^2_L^2(Ω^+).


Testing (<ref>) with sû, integrating by parts and taking the real part gives

    s û^2_|s|,Ω^+= -⟨∂_ν^+ û, s γ^+ û⟩_Γ≤  C(σ) |s|^1/2û_|s|,Ω^+ψ̂_H^1/2(Γ),

where ψ̂= s γ^+ û is the Laplace transform of ψ. Making use of û_H^1(Ω^+)≤ C(σ) û_|s|,Ω^+ and the bound (<ref>) gives

    û_H^1(Ω^+)≤ C(σ) |s|^5/2/( s)^2f̂_H^-1/2(Γ).
 The stated result then follows from Theorem <ref>.





 §.§ 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:→: there exists β>0 such that

    (ξ-η)( g(ξ)-g(η) ) ≥β  |ξ-η|^2    for all  ξ,η∈.

Furthermore, we assume that the pointwise application of g maps H^1/2(Γ)  to H^-1/2(Γ). As is shown in <cit.> by Sobolev embeddings, this is satisfied if g(ξ) grows at most cubically as |ξ|→∞.

In the following we write for a stepsize τ>0 and a sequence e=(e_n)_n=0^N-1 with e_n=(e_n,i)_i=1^m and e_n,i in a Hilbert space V

    e_ℓ_2^(0:N;V^m) = τ∑_n=0^N-1∑_i=1^m  e_n,i_V^2.

We denote the numerical solution by u^τ=(u_n,i) and the corresponding values of the exact solution by u=(u(t_n+c_iτ)), where in both cases n=0,…,N-1 and i=1,…,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 <ref>.

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

Consider the time discretisation (<ref>) and (<ref>) by the two-stage Radau IIA convolution quadrature method.
Then, there is τ̅>0 such that for stepsizes 0<τ≤τ̅, the error in the boundary values satisfies the bound

    γ^+u^τ - γ^+u_ℓ_2^(0:N;H^1/2(Γ)^2)≤ Cτ^3,

and the error in the exterior domain is bounded by

    u^τ - u_ℓ_2^(0:N;H^1(Ω^+)^2)≤ Cτ^3/2.

The constants C are independent of τ and N with 0<Nτ≤ T, but depend on T.


 
We eliminate φ in the system of boundary integral equations (<ref>) to arrive at a boundary integral equation for ψ,

    L(∂_t)ψ + g(ψ + u̇^inc)  =
     ∂_ν^+ u^inc,

where

    L(s) = s^-1(W(s) - (12 I - K^T(s)) V(s)^-1(12 I - K(s))) =
    -s^-1DtN^+(s)

with the exterior Dirichlet-to-Neumann operator DtN^+(s). It follows from Propositions 17 and 18 (and their proofs) in <cit.>  that, for s ≥σ >0, there exist C(σ) and α(σ)>0 such that

    L(s) _H^-1/2(Γ) H^1/2(Γ)≤ C(σ) |s|/ s,
    
       ⟨ψ, L(s) ψ⟩≥α(σ)   s/|s|^2 ψ_H^1/2(Γ)^2
       for all  ψ∈ H^1/2(Γ),

where ⟨·,·⟩ denotes the anti-duality pairing between H^-1/2(Γ) and H^1/2(Γ).

Thanks to the composition rule, we can do the same for the numerical discretisation (<ref>) and reduce the numerical system to an equation for ψ^, which is just the convolution quadrature time discretisation of (<ref>),

    L(∂_t^)ψ^ + g(ψ^ + u̇^inc)  =
     ∂_ν^+ u^inc.

For the error =ψ^ - ψ with ψ= ((ψ(t_n+c_i))_i=1^m)_n=0^N-1 we then have the error equation

    L(∂_t^) + g(ψ^ + u̇^inc)-g(ψ + u̇^inc)  =
    d

with the defect

    d = L(∂_t^)ψ - ((L(∂_t)ψ(t_n+c_i))_i=1^m)_n=0^N-1,

which is the convolution quadrature error for L(∂_t)ψ.
By Theorem <ref> and our assumption of a sufficiently regular ψ=γ^+u̇, this is bounded by

    d_n _H^-1/2(Γ)≤ C ^3    for  0≤ n≤ T.

Since we can apply the same argument also to spatial derivatives of ψ (in the assumed case of a smooth boundary Γ), we even have

    d_n _H^1/2(Γ)≤ C ^3.

We test (<ref>) with , multiply  with e^-2σt with σ=1/T and integrate from 0 to T. With (<ref>) and the Runge-Kutta convolution coercivity as given by Theorem <ref>, and with the strict monotonicity (<ref>) we conclude that

    α   ∑_n = 0^N e^-2σn((∂_t^)^-1)_n_H^1/2(Γ)^2 +β∑_n = 0^N e^-2σn_n_L_2(Γ)^2 
    
       ≤∑_n = 0^N e^-2σn⟨_n, d_n ⟩

and estimate further

    ⟨_n, d_n ⟩≤_n _L_2(Γ) d_n _L_2(Γ)≤β/2_n _L_2(Γ) +1/2βd_n _L_2(Γ)^2.

We thus find the stability estimate

    (∂_t^)^-1_ℓ_2^(0:N;H^1/2(Γ)^2) + 
    _ℓ_2^(0:N;L_2(Γ)^2)≤ C  d_ℓ_2^(0:N;L_2(Γ)^2).

Since (∂_t^)^-1 = γ^+u^τ - γ^+u, this proves (<ref>).

Let us denote by M(s) = S(s)V^-1(s):H^1/2(Γ) → H^1(Ω^+) the operator that maps Dirichlet data in H^1/2(Γ) to the corresponding solution 
û∈ H^1(Ω^+) of the Helmholtz equation s^2û-Δû=0. By <cit.>, this is bounded for s≥σ>0 by

    M(s) _H^1(Ω^+)  H^1/2(Γ)≤ C(σ) |s|^3/2/ s.

We then have

    u^τ - u = M(∂_t^)γ^+u^τ - 
    (( M(∂_t)γ^+u(t_n+c_i))_i=1^2 )_n=0^N-1
       = M(∂_t^) (γ^+u^τ - γ^+u) + 
     ( M(∂_t^)  γ^+u - (( M(∂_t)γ^+u(t_n+c_i))_i=1^2 )_n=0^N-1).

By Theorem <ref> and the bound for M, the last term is bounded by O(τ^5/2) in the H^1(Ω^+) norm. The first term is only O(τ^3/2), since we lose a factor τ^3/2 from the O(τ^3) error bound for γ^+ u because of the O(|s|^3/2) bound of M(s); this follows from Lemma 5.2 in <cit.> and Parseval's identity.



In a similar way we obtain the following results for the alternative discretisations (<ref>) and (<ref>):


(i) In addition to Proposition <ref>, assume that g has bounded second derivatives. With the discretisation (<ref>) instead of (<ref>), the error bound in the H^1(Ω^+) norm improves to O(τ^5/2), and the ℓ_2^ error in a point x bounded away from the boundary Γ is at most O(τ).

(ii) In addition to Proposition <ref>, assume that g has bounded second and third derivatives. With the discretisation (<ref>) instead of (<ref>), the error bound in the H^1(Ω^+) norm improves to O(τ^3), and the ℓ_2^ error in a point x bounded away from the boundary Γ is at most O(τ^2).


The proofs of these error bounds are very similar to that of Proposition <ref>, using in addition a discrete Gronwall inequality at the end of the estimation of , and an O(|s|^3) bound for the norm of the operator from H^1/2(Γ)→ that maps Dirichlet data to the solution of the Helmholtz equation s^2û-Δû=0 at a point x∈Ω^+ bounded away from Γ, 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.


 If we set L(s)=L(s+σ) and ψ(t)=e^-σ tψ(t) for some σ>0, then the boundary integral equation (<ref>) is equivalent to

    (L(∂_t)ψ)(t) + e^-σ t g(e^σ tψ(t)+ u̇^inc(t))  = e^-σ t∂_ν^+ u^inc(t).

By (<ref>), we then have the coercivity estimate for L(s) for all s ≥ 0 (and not just for s ≥σ):

    ⟨ψ, L(s) ψ⟩≥α(σ)   s+σ/|s+σ|^2 ψ_H^1/2(Γ)^2
       for all  ψ∈ H^1/2(Γ).

By Theorem <ref>, the coercivity estimate for the convolution quadrature approximation of L(∂_t)ψ 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 (<ref>) on an interval [0,T] with shift σ=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 <cit.>. 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.




§ NUMERICAL EXPERIMENTS




 §.§ Scattering by the unit sphere


In these experiments we let Ω^+ 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 <cit.>, the solution will also be constant in space. The eigenvalue for the combined operator L(s)  in (<ref>) is given by

    L(s)ψ̂= -s^-1DtN^+(s)ψ̂= (1+1/s)ψ̂,

for any ψ̂ 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

    L^-(s)ψ̂= s^-1DtN^-(s)ψ̂=  (-1/s+1+e^-2s/1-e^-2s)ψ̂,

again for ψ̂ 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 (<ref>) and (<ref>). Namely

    |L(s)| ≤ C(σ),      L(s) ≥ 1

and

    |L^-(s)| ≤ C(σ),      L^-(s) ≥α(σ).

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

    L^-(∂_t)ψ +g(ψ+u̇^inc) = 0.


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 ℓ^τ_2 norm. We have performed the numerical experiments with the following choices of g and u^inc

    g_1(ξ) = 14 ξ+ξ|ξ|,      g_2(ξ) = 14 ξ+ξ^3,      u^inc(t) = 2e^-10(t-5/2)^2

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 <ref> 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 <ref>. Note that the two solutions have a similar shape, but  a closer look at the derivative in Figure <ref> reveals that one is smooth and the other only once continuously differentiable.






For the interior problem, as L(s) ≥ 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 <ref>.





 §.§ 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 ψ and piecewise constant boundary element space to approximate the Neumann trace φ and the time-discretisation is performed using the two-stage Radau IIA method. The images of the solution are shown in Figure <ref>.








§ ACKNOWLEDGEMENT

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

1

BamH
A. Bamberger and T. Ha Duong.
Formulation variationelle espace-temps pour le calcul par potentiel
  retardé d'une onde acoustique.
Math. Meth. Appl. Sci., 8:405–435, 1986.

BamH2
A. Bamberger and T. Ha Duong.
Formulation variationnelle pour le calcul de la diffraction d'une
  onde acoustique par une surface rigide.
Math. Methods Appl. Sci., 8(4):598–608, 1986.

BanK
L. Banjai and M. Kachanovska. 
Fast convolution quadrature for the wave equation in three dimensions. 
J. Comp. Phys., 279, 103–126, 2014.

BanL
L. Banjai and C. Lubich.
An error analysis of Runge-Kutta convolution quadrature.
BIT 51:483–496, 2011.

BanLM
L. Banjai, C. Lubich, and J. M. Melenk.
Runge-Kutta convolution quadrature for operators arising in wave
  propagation.
Numer. Math., 119(1):1–20, 2011.

BanLS
L. Banjai, C. Lubich, and F.-J. Sayas.
Stable numerical coupling of exterior and interior problems for the
  wave equation.
Numer. Math., 129(4): 611–646, 2015.

BanMS
L. Banjai,  M. Messner, and M. Schanz. 
Runge-Kutta convolution quadrature for the boundary element method. 
Computer Methods in Applied Mechanics and Engineering, 245, 90–101, 2012.

BanR
L. Banjai and A. Rieder. 
Convolution quadrature for the wave equation with a non-linear impedance boundary condition.
arXiv preprint arXiv:1604.05212 (2016).

Ebe
S. Eberle.
The elastic wave equation and the stable numerical coupling of its interior and exterior problems.
Preprint, Univ. Tuebingen, na.uni-tuebingen.de/preprints.shtml, 2016.

HaiL
E. Hairer and C. Lubich.
On the stability of Volterra Runge–Kutta methods.
SIAM J. Numer. Anal. 21:123–135, 1984.

HaiW
E. Hairer and G. Wanner. Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Springer, 1996.

KovL
B. Kovács and C. Lubich.
Stable and convergent fully discrete interior-exterior coupling of Maxwell's equations.
Preprint,  arXiv:1605.04086 (2016). To appear in Numer. Math.

LalS
A. R. Laliena and F.-J. Sayas.
Theoretical aspects of the application of convolution quadrature to
  scattering of acoustic waves.
Numer. Math., 112(4):637–678, 2009.

Lub94
C. Lubich.
On the multistep time discretization of linear initial-boundary value
  problems and their boundary integral equations.
Numer. Math., 67:365–389, 1994.

LubO87 
C. Lubich and A. Ostermann.
Multi-grid dynamic iteration for parabolic equations.
BIT 27:216–234, 1987.

LubO
C. Lubich and A. Ostermann.
Runge-Kutta methods for parabolic equations and convolution
  quadrature.
Math. Comp., 60(201):105–131, 1993.

Ned
J.-C. Nédélec.
Acoustic and electromagnetic equations, volume 144 of Applied Mathematical Sciences.
Springer-Verlag, New York, 2001.


Neu
J. von Neumann.
Eine Spektraltheorie für allgemeine Operatoren eines unitären Raumes.
Math. Nachrichten, 4:258–281, 1951.

Say16 
F. Sayas. 
Retarded Potentials and Time Domain Boundary Integral Equations: A Road Map. 
Springer, 2016.

SchLL
A. Schädle, M. López-Fernández, and C. Lubich.
Fast and oblivious convolution quadrature.
SIAM J. Sci. Comput., 28(2):421–438, 2006.

WanW
X. Wang and D. Weile. 
Implicit Runge-Kutta methods for the discretisation of time domain integral equations. 
IEEE Transactions on Antennas and Propagation, 59(12): 4651–4663, 2011.





