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

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.


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 [6] 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 [8].
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 [14,19]. 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 [6,9,12] 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 [5].
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 [8]. 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 [6], where it is used as a basic result in studying boundary integral operators for the acoustic wave equation; see also [12] for Maxwell's equation and [9] 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 Re s > σ, continuous for Re s ≥ σ. We assume the uniform bounds, with some real exponent µ, 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 defines a continuous and exponentially bounded function k , which has as its kth distributional derivative. We write the convolution with as 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 .
Property 1. is known to be satisfied for the Laplace transforms of various boundary integral operators for wave equations [1,6,12,13,19], 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 [6] 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 σ = 0, it will be shown to be preserved for all algebraically stable Runge-Kutta methods.

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 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 We always assume that the Runge-Kutta matrix A is invertible.
As has been shown in [16,20,4,5], and in applications to wave propagation problems further in [3,7,8,21], Runge-Kutta methods can be used to construct convolution quadrature methods that enjoy favourable properties. Here one uses the Runge-Kutta differentiation symbol In fact, the Sherman-Morrison formula then yields 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 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 , 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: where e m = (0, . . . , 0, 1) T is the mth unit vector. We recall the composition rule For λ ∈ C, the convolution quadrature (∂ τ t −λ) −1 f (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 [4,5,16]. The result of [5], which is relevant for operators L(s) arising in wave propagation, will be restated and extended to the internal stages in Section 4.

Numerical range of the Runge-Kutta differentiation symbol
We now consider methods that are algebraically stable: Gauss methods and Radau IIA methods are widely used classes of methods that satisfy this condition. We refer the reader to [11] for background literature on Runge-Kutta methods and their stability notions. We consider the weighted inner product on C m , We have the following characterisation.
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 [15, p. 232, (6.19)]. For the convenience of the reader we include the short proof.
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 Lemma 3.2 (a) For the two-stage Radau IIA method and for the b-weighted inner product (3.2) and corresponding norm | · | we have 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, Clearly, (3.3) implies (3.4) 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 (3.4).
Proof (a) For the two-stage Radau IIA method we find Denoting the diagonal matrix of the weights by B = diag(3/4, 1/4), we note or in other words, W T B 1/2 is an orthogonal matrix (with respect to the Euclidean inner product), and m is a skew-symmetric matrix with β m = 0 for the Gauss method and β m > 0 for the Radau IIA method. We write Now, by the definition of ∆(ζ) and the above-mentioned property of W −1 A W = X together with W T b = e 1 , the matrix W −1 ∆(ζ) −1 W is the sum of a skewhermitian 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: 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 B, which is of codimension 1 or 2 for Gauss or Radau, respectively, then we obtain from the above formulas that Re (w, ∆(ζ)w) = 0 in contradiction to (3.4).
As we will show in Theorem 3.1 below, Runge-Kutta convolution quadrature with (3.4) preserves the coercivity property of Theorem 2.1 for arbitrary abscissa σ ≥ 0, while general algebraically stable methods preserve it in the case σ = 0. Before we state and prove this theorem in Section 3.4, 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 antiduality denoted by ·, · . On C 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 ∈ C 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 .
Let the matrix S ∈ C m×m be such that Then, This result can be viewed as an extension of a theorem of von Neumann [18] (see also [11, p. 179]), 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 ).
Proof The proof adapts Michel Crouzeix's proof of von Neumann's theorem as given in [11, p. 179 f.]. 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 and we observe that S = S(1) and Re (w, S(z)w) = 1 2 (Re z) Re (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 Re 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 Re z ≥ 0.
If Re z = 0, then the matrix S(z) is normal, and hence the desired inequality is valid for S(z) with Re z = 0. The function 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), which is the desired inequality.
Let the matrix S ∈ C m×m be such that all its eigenvalues either have positive real part or are zero, and Then, 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
Theorem 3.1 Let the m-stage Runge-Kutta method satisfy (3.4) for some inner product (·, ·), 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 τ > 0 and with σ = σ/c, 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  Re ρ n f n , Here (3.5) used in Lemma 3.3 yields Moreover, again by Parseval's formula, which yields the result.

Error bounds of Runge-Kutta convolution quadrature
In this section we restate the result of [5] 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 Re 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: for all Re s > σ.
The constant C is independent of τ and f , but does depend onτ , T , and the constants in (4.1).
The proof in [5] 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.) Theorem 4.2 Let K satisfy (4.1) 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 ], The constant C is independent of τ and f , but does depend onτ , T , and the constants in (4.1).

A non-linear scattering problem
We consider the wave equation on an exterior smooth domain Ω + ⊂ R 3 . Following [8], we search for a function u(·, t) ∈ H 1 (Ω + ) satisfying the weak form of the wave equationü = ∆u in Ω + (5.1) with zero initial conditions and with the non-linear boundary condition where ∂ + ν is the outer normal derivative on the boundary Γ of Ω + , where g : R → 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 Γ 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. [8].

Boundary integral equation and representation formula
Using the standard notation of the boundary integral operators for the Helmholtz equation s 2 u − ∆u = 0 (Re s > 0) as used, for example, in [13,19] and [6,8], we denote by 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

.3) and the related operator
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 [8]) The solution u is then obtained from the representation formula We will address the question as to what are the approximation properties when the temporal convolutions in (5.5) and (5.6) are discretised by Runge-Kutta convolution quadrature. Since this will not turn out fully satisfactory, we will further consider time-differentiated versions of (5.5).
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 5.1 yields property 1. of Theorem 2.1 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 3.1, the boundary integral equation (5.5) is discretised in time with a stepsize τ > 0 over a time interval (0, T ) with T = N τ by where (ϕ τ , ψ τ ) = (ϕ n , ψ n ) N −1 n=0 with (ϕ n , ψ n ) = (ϕ n,i , ψ n,i ) m i=1 and (ϕ n,i , ψ n,i ) ≈ (ϕ(t n + c i τ ), ψ(t n + c i τ )) is the numerical approximation that is to be computed, andu inc = (u inc n ) N −1 n=0 withu inc n = (u inc (t n + c i τ )) m i=1 . The function g acts componentwise. At the nth time step, a non-linear system of equations of the following form needs to be solved: where the dots represent known terms. This has a unique solution, because the eigenvalues of ∆(0) = A −1 have positive real part, and Lemma 5.1 and the monotonicity of g then yield the unique existence of the solution by the Browder-Minty theorem; cf. [8] for the analogous situation for multistep-based convolution quadrature.
As an alternative to (5.7), we further consider the time discretisation of the time-differentiated boundary integral equation: 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: 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 (5.6):

Error bounds for the linear case
We consider first the case of a linear impedance function be the solution approximation obtained by the discretised representation formula (5.10) with either of the discretised boundary integral equations (5.7) or (5.8) or (5.9). 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 Γ .
Proposition 5.1 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 dist(x, Γ ) ≥ δ > 0, the following optimal-order error bound is satisfied in the linear situation described above: for 0 ≤ t n = nτ ≤ T , The exact solution u(x, t) is given by the representation formula (5.6) with For x ∈ Ω + we define the operators S x (s) : H −1/2 (Γ ) → C and D x (s) : The first bound is proved in [5,Lemma 6] and the second bound is proved similarly. We thus have With the above operator bounds we obtain for Re s ≥ σ > 0 and dist(x, Re s e −δ Re s . (5.14) By the composition rule, the numerical solution obtained by (5.7) and (5.10) is given as contains the values of f at the points t n + c i τ . If we take instead (5.8) or (5.9), then we have , respectively. In view of (5.14), Theorem 4.1 then yields the result.
The situation is different if we consider the H 1 (Ω + ) norm of the error.

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 : R → R: there exists β > 0 such that Furthermore, we assume that the pointwise application of g maps H 1/2 (Γ ) to H −1/2 (Γ ). As is shown in [8] by Sobolev embeddings, this is satisfied if g(ξ) grows at most cubically as |ξ| → ∞. Consider the time discretisation (5.7) and (5.10) 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 18) and the error in the exterior domain is bounded by The constants C are independent of τ and N with 0 < N τ ≤ T , but depend on T .
Proof We eliminate ϕ in the system of boundary integral equations (5.5) to arrive at a boundary integral equation for ψ, with the exterior Dirichlet-to-Neumann operator DtN + (s). It follows from Propositions 17 and 18 (and their proofs) in [13] that, for Re s ≥ σ > 0, there exist C(σ) and α(σ) > 0 such that for all ψ ∈ H 1/2 (Γ ), (5.22) 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 (5.7) and reduce the numerical system to an equation for ψ τ , which is just the convolution quadrature time discretisation of (5.20), we then have the error equation , which is the convolution quadrature error for L(∂ t )ψ. By Theorem 4.2 and our assumption of a sufficiently regular ψ = γ +u , this is bounded by Since we can apply the same argument also to spatial derivatives of ψ (in the assumed case of a smooth boundary Γ ), we even have We test (5.24) with , multiply with e −2 σt with σ = 1/T and integrate from 0 to T . With (5.22) and the Runge-Kutta convolution coercivity as given by Theorem 3.1, and with the strict monotonicity (5.24) we conclude that Remark 5.1 If we set L(s) = L(s + σ) and ψ(t) = e −σt ψ(t) for some σ > 0, then the boundary integral equation (5.20) is equivalent to L(∂ t ) ψ (t) + e −σt g(e σt ψ(t) +u inc (t)) = e −σt ∂ + ν u inc (t). (5.26) By (5.22), we then have the coercivity estimate for L(s) for all Re s ≥ 0 (and not just for Re s ≥ σ): Re ψ, L(s)ψ ≥ α(σ) Re s + σ |s + σ| 2 ψ 2 for all ψ ∈ H 1/2 (Γ ).
By Theorem 3.1, 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 (5.26) 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 timedependent boundary integral equations [1]. 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.

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 [17], the solution will also be constant in space. The eigenvalue for the combined operator L(s) in (5.20) is given by 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 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 (ξ) = 1 4 ξ + ξ|ξ|, g 2 (ξ) = 1 4 ξ + ξ 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 6.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 6.2. Note that the two solutions have a similar shape, but a closer look at the derivative in Figure 6.3 reveals that one is smooth and the other only once continuously differentiable.
For the interior problem, as Re 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 6

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