Runge–Kutta approximation for C0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_0$$\end{document}-semigroups in the graph norm with applications to time domain boundary integral equations

We consider the approximation of an abstract evolution problem with inhomogeneous side constraint using A-stable Runge–Kutta methods. We derive a priori estimates in norms other than the underlying Banach space. Most notably, we derive estimates in the graph norm of the generator. These results are used to study convolution quadrature based discretizations of a wave scattering and a heat conduction problem.

case of l 2 ½0; 1 was adressed, which is the case needed for PDEs with inhomogeneous boundary conditions. We point out that in the case of analytic semigroups, Lubich and Ostermann [19] had already established convergence also for inhomogeneous boundary conditions. All of these works focus on establishing convergence rates with respect to the norm of the underlying Banach space. In many applications one needs to establish convergence with respect to other norms, for example, in order to be able to bound boundary traces of the solution. Most notably, one might be interested in convergence of A H u, where A H is an extension of the generator that disregards boundary conditions. If u is assumed to be in domðAÞ, we get A H u ¼ Au and the convergence result can be easily established by using the fact that the time evolution commutes with the generator of the underlying semigroup (both in the continuous and discrete settings). If the boundary conditions are inhomogeneous, such a strategy cannot be pursued. It is the goal of this paper to establish convergence results for A H u also for the case uðtÞ 2 domðA l Þ for l 2 ½0; 1, again using the theory of interpolation spaces.
Similarly it is sometimes useful to compute discrete integrals of the time evolution by reusing the same Runge-Kutta method. Also in this case, we establish rigorous convergence rates.
Our interest in such estimates originally arose from the study of time domain boundary integral equations and their discretization using convolution quadrature (CQ). It has already been noticed in the early works (see e.g. [19]) that such discretizations have a strong relation to the Runge-Kutta approximation of the underlying semigroup. This approach of studying TDBIEs in a strictly time-domain way has recently garnered a lot of interest, see [3,13,15] and the monograph [31], as it potentially allows sharper bounds than the more standard Laplace domain based approach. Similar techniques have even been extended to the case of certain nonlinear problems in [4]. This paper can be seen as our latest addition to this effort. While the convergence rates provided by the Laplacedomain approach in [2] and the results in this current paper are essentially the same, the present new approach provides better insight into the dependence on the end-time of the computation (quadratic vs. general unknown polynomial behavior). This suggest that the present approach might be better suited for analyzing long term computations. It also fits more naturally with the time-domain analysis of the continuous problem and space discretization, as for example presented in [13].
The paper is structured as follows. Section 2 introduces the abstract setting and fixes notation, most notably for working with Runge-Kutta methods. Section 3 then contains the main estimates. Starting by summarizing known results from [1] in Sect. 3.1, we then formulate the main new results of this article in Sect. 3.2. After proving some preparatory lemmas related to Runge-Kutta methods in Sects. 4 and 5, we provide the proofs of the main estimates in Sect. 6. In Sect. 7, we show how our setting simplifies if we restrict our view to a subclass of admissible operators. In Sect. 8, to showcase how the theory developed in this paper is useful for this class of problems, we consider a simple exterior scattering problem in Sect. 8.3 and a heat transmission problem in Sect. 8.5. We note that Sect. 8.3 showcases the need for the bound on the discrete integral of the result, whereas Sect. 8.5 was chosen because, in order to bound the main quantity of interest on the boundary, we need to apply a trace theorem. This necessitates the use of the graph norm estimate.

Problem setting
We start by fixing the general setting used for the rest of the paper, first with respect to the equation to be solved and then with respect to its discretization. We assume that A :¼ A H j ker B generates a C 0 -semigroup and that B admits a bounded right inverse E such that range E & kerðI À A H Þ, where I : X ! X is the identity operator.
We are given u 0 2 domðAÞ and data functions F 2 C 1 ð½0; T; XÞ, N 2 C 1 ð½0; T; MÞ, and we consider the problem: find u 2 C 1 ð½0; T; XÞ such that For conditions on the well-posedness of this problem, see [13]. We start by recalling the following consequence of the Hille-Yosida theorem. When working with Runge-Kutta methods, it is useful to use a calculus that allows one to apply rational functions to (unbounded) operators, as long as the poles of the function are compatible with the spectrum of the operator.
Definition 2.2 (Rational functions of operators) Let q be a rational function that is bounded at infinity. Let K be the set of poles of q, which we can write in the form (note that we allow for some of the factors in the numerator to be constant) If A : domðAÞ & X ! X is a linear operator such that rðAÞ \ K ¼ ;, we define SN Partial Differential Equations and Applications qðAÞ :¼ c 0 ðc 1 I þ ðc 1 k 1 À 1ÞðA À k 1 IÞ À1 Þ Á Á Á ðc n I þ ðc n k n À 1ÞðA À k n IÞ À1 Þ: ð2:3Þ It is easy to see that different reorderings of the factors in the numerator and denominator of q produce the same result and that each factor in the definition of q(A) is a bounded linear operator in X since k i 6 2 rðAÞ: The bounded linear operator qðAÞ : X ! X satisfies kqðAÞk X!X C q 1 þ À max k2K kðA À kIÞ À1 k X!X Á n : ð2:4Þ The error estimates of this paper use the theory of interpolation spaces. For Banach spaces X 1 & X 0 with continuous embedding and l 2 ð0; 1Þ, we define the space ½X 0 ; X 1 l;1 using real interpolation with the following norm: We will not go into details of the definitions and instead refer to [35,36] or [22,Appendix B]. For simplicity of notation we often drop the second parameter 1 and just write ½X 0 ; X 1 l . The most important property is the following: a bounded linear operator T : X 0 ! Y 0 and X 1 ! Y 1 with X 1 X 0 and Y 1 Y 0 is also a bounded operator mapping ½X 0 ; X 1 l ! ½Y 0 ; Y 1 l with the following norm bound We also note that for l 1 l 2 , the spaces are nested, i.e., ½X 0 ; X 1 l 2 ½X 0 ; X 1 l 1 with continuous embedding. For notational convenience we write ½X 0 ; X 1 0 :¼ X 0 and ½X 0 ; X 1 1 :¼ X 1 . We will be interested in a collection of spaces defined by interpolating the domains of the powers of the operator A. The details of this construction can be found, for example in [11]. Definition 2.3 (Sobolev towers) Let A be a closed operator on a Banach space X . For l 2 N 0 , we define the following spaces X 0 :¼ dom A 0 ð Þ :¼ X and X l :¼ dom A l ð Þ, equipped with the following norm For l 2 ½0; 1Þ, we define X l :¼ X blc ; X blcþ1 Â Ã lÀblc by interpolation.
We sometimes consider domðAÞ as a Banach space. It is to be understood carrying the graph norm, same as X 1 .

Runge-Kutta approximation and discrete stage derivative
An m-stage Runge-Kutta method is given by its Butcher tableau, characterized by Q 2 R mÂm and b, c 2 R m . The Runge-Kutta approximation of the problem (2.1a-2.1c) starts at u k 0 :¼ u 0 and then computes for n ! 0 the stage vector U k n 2 X m and the step approximation u k nþ1 2 X by solving We have used the following notation (the spaces Y and Z are generic): (a) For a function G : ½0; T ! Y we write Gðt n þ kcÞ :¼ Gðt n þ kc 1 Þ; . . . ; Gðt n þ kc m Þ ð Þ > 2 Y m : (b) For a matrix S 2 R mÂm and an operator C : Y ! Z we write (c) For the vector b and an operator C : (d) I is the m Â m identity matrix, and 1 ¼ ð1; Á Á Á ; 1Þ > . (e) We admit shortened expressions such as QFðt n þ kcÞ :¼ ðQ IÞFðt n þ kcÞ; b > Fðt n þ kcÞ :¼ ðb > IÞFðt n þ kcÞ: The following lemma involving inversion of matrices of operators associated to an operator can be proved by taking the Jordan canonical form of the matrix S.

Lemma 2.4 If A :
domðAÞ & X ! X is a linear operator on a Banach space X and S 2 C mÂm satisfies rðAÞ \ rðSÞ ¼ ;, then is invertible. Furthermore, there exists a constant C S , depending only on S, such that Under Assumption 2.I, the internal stage computation in the RK method can be decomposed in the following form: Y k n : ¼ ðI EÞNðt n þ k cÞ; ð2:8aÞ Z k n À kðQ AÞZ k n ¼ 1u k n À Y k n þ kQðY k n þ Fðt n þ k cÞÞ; ð2:8bÞ In (2.8b) we look for Z k n 2 ðdomðAÞÞ m . The stability function of the Runge-Kutta method is the rational function rðzÞ :¼ 1 þ zb > ðI À zQÞ À1 1. We will not consider the full class of Runge-Kutta methods, but will restrict our considerations to those satisfying the following Assumptions: The stability function r does not have poles in fz : Re z\0g, and rðitÞ j j 1 for all t 2 R (i.e., the method is A-stable). Equivalently, jrðzÞj\1 for all z with negative real part.
We note that Assumption 2.II (i) implies that the following limit exists lim z!1 and that r is a rational function with poles only in C þ and bounded at infinity.
The computation of the internal stages in the numerical approximation (2.7a-2.7c) requires the inversion of I I À kðQ AÞ ¼ ðQ IÞðQ À1 I À I ðk AÞÞ; as can be seen from the equivalent form (2.8a-2.8c).
If A is the infinitesimal generator of a C 0 -semigroup and x and M are given by Proposition 2.1 and if we choose (recall that rðQÞ & C þ ) then the RK method can be applied for any 0\k k 0 . By Proposition 2.1 and Lemma 2.4, it follows that This quantity is relevant for the study of the error propagation in the Runge-Kutta method.
Given an RK method, we consider the following matrix-valued rational function (The verification that these two formulas correspond to the same matrix is simple by using the Sherman-Morrison-Woodbury formula.) This matrix is related to the discrete differentiation process associated to an RK method satisfying Assumption 2.II: on the one hand k À1 dðzÞ is the discrete symbol associated to the discrete operational calculus built with the RK method [19]; on the other hand, a direct interpretation of this symbol is possible using the Z-transformation (see [14,Sect. 6]). Given a sequence U :¼ fU n g (tagged from n ! 0) on a space, its Z-transform Z is the formal series b For a detailed treatment on formal power series, see [12].
Definition 2.5 Let U :¼ fU n g and V :¼ fV n g be two sequences in X m and let b U and b V be their respective Z-transforms. If The above definition is consistent with the RK discrete operational calculus of Lubich and Ostermann, see Sect. 8.1 and [19]. We now show an explicit form of the computation of o k and its inverse.
Lemma 2.6 If U ¼ fU n g is a sequence in X m , then X :¼ ðo k Þ À1 U can be computed with the recurrence x 0 :¼ 0; X n :¼ 1x n þ kQU n ; x nþ1 :¼ x n þ kb > U n ¼ rð1Þx n þ b > Q À1 X n ; ð2:13Þ and V :¼ o k U can be computed with the inverse recurrence ð2:14Þ Proof The proof of (2.13) is a simple exercise in Z-transforms, while (2.14) follows from (2.13) by writing U n in terms of X n (and changing names to the sequences). h The first result of Lemma 2.6 expresses the fact that if we apply the RK method to the equation and rð1Þ ¼ 0: For stiffly accurate methods, taking the discrete derivative of a stage vector consisting of samples taken from a continuous function is particularly simple: Lemma 2.7 Let t7 !FðtÞ be a continuous function with Fð0Þ ¼ 0. For stiffly accurate RK methods the sequence G :¼ o k F with F n ¼ Fðt n þ kcÞ satisfies G n ¼ k À1 Q À1 ðFðt n þ kcÞ À 1Fðt n ÞÞ: Proof For stiffly accurate methods we have rð1Þ ¼ 0 and therefore However, since c m ¼ 1, we have e > m Fðt nÀ1 þ kcÞ ¼ Fðt nÀ1 þ kc m Þ ¼ Fðt n Þ; which proves the result. h We also make the following optional assumption, which allows us to increase the convergence order in some cases.

Error estimates
We are now in a position to formulate the main results of this article and put them into context with previous results, most notably from [1].
To simplify notation, we will write for v 2 Cð½0; T; X l Þ with l ! 0; For functions f : ½0; T ! Y, we will write ðo À1 f ÞðtÞ :¼ R t 0 f ðsÞds, where Y denotes a generic Banach space.

The estimates of Alonso-Mallo and Palencia
The following two propositions summarize the results of Alonso-Mallo and Palencia [1], rewritten with the notation of the present paper. The 'proofs' which we provide clarify how notation needs to be adapted and how the hypotheses of the main results of [1] are satisfied in our context. Proposition 3.1 ([1, Theorem 1]) Let Assumption 2.I hold and assume that the exact solution u satisfies u 2 C pþ1 0; T ½ ; X l À Á for some l ! 0. Let fu k n g denote the Runge-Kutta approximation from (2.7a-2.7c). Then there exist constants k 0 [ 0 and C [ 0 such that for 0\k k 0 and 0\nk T the following estimate holds: ku ð'Þ k T;l þ ku ðpþ1Þ k T;0 : ð3:1Þ The constant C depends on the Runge-Kutta method, l, and the constants M and x from (2.2). The constant k 0 depends only on x and the Runge-Kutta method.
Proof We only remark on the differences in notation. A different definition of interpolation spaces is given in [1], but the proof only relies on estimates of the form (2.6). The choice of k 0 follows from the fact that it is only needed to ensure that ðI À k Q AÞ is invertible, see (2.10). The assumption l p À q in [1, Theorem 1] can be replaced by using the rate minfp; q þ lg in (3.1) as the spaces X l X pÀq are nested for l ! p À q. We also lowered the regularity requirements on the highest derivative compared to their stated result. The fact that this holds true follows from inspection of the proof. See also Lemma 5.9 for the key ingredient. h For certain Runge-Kutta methods, these estimates can be improved: Theorem 2]) Let the assumptions of Proposition 3.1 hold and assume that, in addition, the RK method satisfies Assumption 2.III. Then there exist constants k 0 [ 0, C [ 0 such that for 0\k k 0 and 0\nk T the following improved estimate holds: The constant C depends on the Runge-Kutta method, l, and the constants M and x from (2.2); k 0 depends only on the constant x and the Runge-Kutta method.
Proof Again, this is just a reformulation of [1, Theorem 2]. We first note that, due to our assumption on rð1Þ, we are always in the case m ¼ 0 of [1]. Since we assumed that on the imaginary axis rðitÞ j j\1 for 0 6 ¼ t 2 R, we directly note that for sufficiently small k 0 , all the zeros of rðzÞ À 1 except z ¼ 0 satisfy Re z [ k 0 x. By the resolvent bound (2.2) we can therefore estimate for k k 0 ðzI À kAÞ À1 i.e., we have a uniform resolvent bound in the set Z a;d in [1]. We also note that we reformulated the convergence rate such that we do not have the restriction l p À q À 1, since the exceptional cases are already covered by Proposition 3.1. h

Remark 3.3
The assumption rðzÞ j j\1 for ReðzÞ 0 and rð1Þ 6 ¼ 1 is satisfied by the Radau IIA family of Runge-Kutta methods, but is violated by the Gauss methods, which satisfy rðzÞ j j ¼ 1 on the imaginary axis.

New results in this article
In this section we present some a priori estimates for the convergence of Runge-Kutta methods when applied to the abstract problem (2.1a-2.1c). These can be seen as a continuation of [1] to the case where the boundary conditions are not given exactly but stem from computing discrete integrals and differentials using the same Runge-Kutta method.
Theorem 3.4 (Integrated estimate) Let u solve (2.1a-2.1c) with u 0 ¼ 0 and assume that for some l ! 0 we have u 2 C p ð½0; T; X l Þ; EN; F 2 C pÀ1 ð½0; T; X l Þ \ C p ð½0; T; X 0 Þ: Let U k ¼ fU k n g and let u k ¼ fu k n g be the discrete approximation given by (2.7a-2.7c) for a method satisfying Assumption 2.II. If X k :¼ ðo k Þ À1 U k and we define x k ¼ fx k n g with the recurrence then there exists a constant k 0 [ 0 such that for all k\k 0 and n 2 N with nk T the following estimate holds: '¼q ku ð'Þ k T;l þ kEN ð'Þ k T;l þ kF ð'Þ k T;l þ ku ðpÞ k T;l þ kEN ðpÞ k T;0 þ kF ðpÞ k T;0 # þ C T 2 q k ðTÞk p kEN ðpÞ k T;0 þ kF ðpÞ k T;0 : If Assumption 2.III holds and if we assume the stronger regularities u 2 C pþ1 ð½0; T; X l Þ; F 2 C p ð½0; T; X l Þ; EN 2 C p ð½0; T; X l Þ; then kxðt n Þ À x k n k X Cð1 þ TÞq k ðTÞk minfqþlþ2;pg " X p '¼q ku ð'Þ k T;l þ kEN ð'Þ k T;l þ kF ð'Þ k T;l þ ku ðpþ1Þ k T;l # þ C T 2 q k ðTÞk p kEN ðpÞ k T;0 þ kF ðpÞ k T;0 : The constant k 0 depends only on x from (2.2) and the Runge-Kutta method. If x ¼ 0 then k 0 can be chosen arbitrarily large. C depends on x, M from (2.2), the Runge-Kutta method, and l.
, then there exists a constant k 0 [ 0 such that for all k\k 0 and n ! 1 such that nk T the following estimate holds: If, in addition, the method satisfies Assumption 2.III and u 2 C pþ2 ð½0; T; X l Þ; EN; F 2 C pþ1 ð½0; T; X l Þ \ C pþ2 ð½0; T; X 0 Þ; then kvðt n Þ À v k n k X þ kA H ðuðt n Þ À u k n Þk X Cð1 þ TÞq k ðTÞk minfqþl;pg X pþ1 '¼qþ1 À ku ð'Þ k T;l þ kEN ð'Þ k T;l þ kF ð'Þ k T;l Á þ ku ðpþ2Þ k T;l þ kEN ðpþ2Þ k T;0 þ kF ðpþ2Þ k T;0 : The constant k 0 depends only on x from (2.2) and the Runge-Kutta method. If x ¼ 0, then k 0 can be chosen arbitrarily large. C depends on x, M from (2.2), the Runge-Kutta method, and l.
Remark 3.6 Most of the effort in proving the above theorem is done in order to obtain a convergence rate higher than q, even though the constraint in the stages is approximated only with order q. This is possible by exploiting the additional structure of the discretization error of the side constraint.

Remark 3.7
We formulated all our results for homogeneous initial conditions, since it is sufficient for our purposes in time domain BEM and convolution quadrature. It should be possible to generalize these results to the case of u 0 2 domðA s Þ for sufficiently large s ! 1 by considering the evolution of the semigroup with inhomogeneous side constraint but homogeneous initial condition and the semigroup of homogeneous constraint but inhomogeneous u 0 separately.
Remark 3.8 The loss of order by 1 in Theorem 3.5 compared to Propositions 3.1 and 3.2 is to be expected. Indeed, if we look at the case u 2 domðA l Þ for l ! 1, this means A H u 2 domðA lÀ1 Þ. Applying Proposition 3.2 to this semigroup then also gives a reduced order of k minðqþl;pÞ .

Some computations related to the main theorems
We will collect the sampled data and the stage and step parts of the solutions in four formal series If the data functions are polynomially bounded in time, the series in (4.1a) are convergent (in X m and M m respectively) with at least unit radius of convergence. Because of the equivalent formulation of the numerical method in the form (2.8a-2.8c), and using (2.10), it follows that for k k 0 [with k 0 chosen using (2.9)], the numerical solution is at least bounded in the form U k n X .C n . Thus, the two series in (4.1b) also converge on a sufficiently small disk.
Proof Let us start by proving a simple result: the discrete equations (2.7a) and (2.7c) hold if and only if (2.7a) and hold. To see this, note that (2.7a) is equivalent to and therefore (2.7a) and (2.7c) imply or equivalently (4.3). The reciprocal statement is proved similarly. The recurrence (4.3) is equivalent to (4.2c). At the same time, the recurrence (2.7a) is equivalent to playing the role of a discrete Dirac delta at time t ¼ 0.
rðdðzÞÞ rðQ À1 Þ [ fw 2 C : rðwÞz ¼ 1g: In particular, if the Runge-Kutta method is A-stable (Assumption 2.II), then We need a corollary to the previous result: Proof In view of of Lemma 4.2, since rðQÞ is finite, independent of z, and contained in C þ , we are mainly concerned with the set fw 2 C : rðwÞz ¼ 1g. We first note that [ z j j r0 fw 2 C : rðwÞz ¼ 1g fw 2 C : rðwÞ j j! 1=r 0 g: Second, we observe that by taking d 0 small enough, we can ensure that w7 !rðwÞ is continuous for ReðwÞ d 0 and thus fw 2 C : rðwÞ j j! 1=r 0 g \ fw 2 C : ReðwÞ d 0 g ¼ rj fReðwÞ d0g À1 À ½1=r 0 ; 1Þ Á is a closed set. Third, by considering the limit along the imaginary axis, we get SN Partial Differential Equations and Applications rð1Þ j j¼ lim n!1 rðinÞ j j 1: Thus, for w j j sufficiently large, it holds that rðwÞ j j 1=r 0 . Overall, we get that is a compact set with empty intersection with the imaginary axis. Thus, it must have a positive distance from it. These observations and Lemma 4.2 conclude the proof. h has a unique solution for arbitrary b F 2 X m and b N 2 M m . If x ¼ 0 in Proposition 2.1, then there are no restrictions on k, and the results holds for all jzj\1.
Proof Assume first that S 2 C mÂm is such that rðSÞ & fz : Re z [ xg and consider the problem where E is the lifting operator of Assumption 2.I) and then seek b W 2 ðdom ðAÞÞ m satisfying This problem is uniquely solvable by Lemma 2.4, since rðAÞ & fz : Re z xg and therefore rðAÞ \ rðSÞ ¼ ;: W , which solves (4.6a, 4.6b). To see uniqueness, one observes that the difference of two solutions of (4.6a, 4.6b) solves the homogeneous problem ( b By Corollary 4.3, the union of the spectra of dðzÞ for jzj r 0 has a positive distance dðr 0 Þ [ 0 from the imaginary axis. If we take k 0 \dðr 0 Þ=x, then rðk À1 dðzÞÞ & fs : Re s [ xg for all jzj r 0 and k k 0 . When x ¼ 0, we can take any k 0 . By the previous considerations this implies unique solvability. h Proposition 4.5 Let U k ¼ fU k n g and u k ¼ fu k n g be sequences satisfying (2.7a-2.7c) with u k which proves that By Proposition 4.1, Eq. (4.9a-4.9c) are equivalent to (4.7a-4.7c). Finally (4.8a) follows from (4.2a), while (4.8b) follows from (4.4) and (4.8a). h Proposition 4.6 Let U k ¼ fU k n g and u k ¼ fu k n g be sequences satisfying (2.7a-2.7c) with for data x k Proof Follow the proof of Proposition 4.5. h

Some Lemmas regarding Runge-Kutta methods
In order to shorten the statements of the results of this section, in all of them we will understand that: (1) We have an RK method with coefficients Q; b; c satisfying Assumption 2.II (invertibility of Q and A-stability). The method has classical order p and stage order q. (2) We have an operator A in X that is the generator of a C 0 -semigroup, characterized by the quantities M and x of Proposition 2.1. The associated Sobolev tower fX l g, obtained by interpolation of domðA l Þ for positive integer values of l, will also be used.
The following lemma will be used at a key point in the arguments below.
Lemma 5.1 Let A be a linear operator in X and q be a rational function bounded at infinity whose poles are outside rðAÞ. The following properties hold: Using this result for each of the factors in the definition (2.3) the result follows. To prove (b) note first that p is rational, bounded at infinity, and that rðAÞ does not intersect the set of poles of p. Using Definition 2.2, we have pðAÞ ¼ A À' qðAÞ ¼ qðAÞA À' , and the result follows. h We start by recalling some simple facts about RK methods that we will need in the sequel. Using the notation c ' :¼ ðc ' 1 ; . . .; c ' m Þ > , the following equalities (order conditions) hold (see e.g. [1,25]): For a stiffly accurate method we have (2.15) and therefore The following result is well-known. We just summarize it for ease of reference later on.
Lemma 5.2 (Discrete antiderivative and RK quadrature) Let f : ½0; T ! X, g :¼ o À1 f , G k ¼ fG k n g ¼ ðo k Þ À1 ff ðt n þ kcÞg and fg k n g be given by the recursion For the errors d k n :¼ gðt n Þ À g k n , and for n such that nk T, we have the estimates kd k n k X CTk p kf ðpÞ k T;0 ; ð5:4aÞ Additionally, at the stage level we have Proof Follows from the fact that the Runge-Kutta method defines a quadrature formula of order p. h

Estimates on rational functions of the operator
The following results in this section are adaptations from [1]. While they focus on the case b ¼ 0, we present the necessary generalizations to b ¼ AE1. We will use the rational functions Note that these rational functions are bounded at infinity and that r ';b ð0Þ ¼ 0. We will also use the vector-valued rational function gðzÞ > :¼ zb > ðI À zQÞ À1 ; ð5:7Þ and note that gð0Þ ¼ 0 and rðzÞ ¼ 1 þ gðzÞ > 1.

Lemma 5.3
The rational functions (5.5) satisfy with q k ðTÞ defined in (2.11). If ' ¼ p and b ¼ 1, the left-hand side of (5.9) is bounded by Cq k ðTÞ. The constant C [ 0 in (5.9) depends only on the Runge-Kutta method, M and x, k 0 , ', and l, but is independent of n and k. If the Runge-Kutta method is stiffly accurate, then the estimate (5.9) also holds for b ¼ À1. If x ¼ 0, then k 0 can be chosen arbitrarily.
Proof We adapt the proof of [1, Lemma 6], which only covers the case b ¼ 0. Consider first the case p À ' À b ! 0 and take any integer l such that 0 l p À ' À b. Then By Lemma 5.3, the rational function r ';b has a zero of order p À ' À b þ 1 at z ¼ 0. The rational function ðrðzÞ À 1Þz l has a zero of order l þ 1 p À ' À b þ 1 at z ¼ 0, and all other zeros are in C þ by A-stability and Assumption 2.III. This implies that the rational function q ';b;l has its poles in Therefore, for k 0 [ 0 sufficiently small we get using (2.4): where C depends on M, x, k 0 , and the RK method. By Lemma 5.1 we have This, (5.11), and applying (2.11) to control rðkAÞ nþ1 by q k ðTÞ, proves (5.9) for integer l p À ' À b. For larger integer values of l, the result does not need to be proved as the maximum rate is already attained. We just have to estimate the X pÀ'Àb norm by the stronger X l norm. For real values of l, we use interpolation. We still need to prove the result when p À ' À b ¼ À1; which can only happen when ' ¼ p and b ¼ 1: We note that r p;1 ð0Þ ¼ 0 and we can therefore argue as in the previous case for l ¼ 0. h Lemma 5.5 If the RK method satisfies Assumption 2.III and k 0 is the value given in Lemma 5.4, then ks n ðk AÞgðk AÞ > k X m !X C q k ðTÞ; ð5:12Þ for all k k 0 and n such that nk T.
Proof Since gð0Þ ¼ 0, we can adapt the proof of Lemma 5.4 to each of the components of the vector-valued function g. The key step is to show that hðzÞ > :¼ ðrðzÞ À 1Þ À1 gðzÞ is bounded at infinity and has all its poles in the set defined in (5.10) and therefore Since the operator s n ðk AÞgðk AÞ > on the left-hand side of (5.12) can be rewritten as ðrðkAÞ nþ1 À IÞhðk AÞ > , the bound (5.12) follows readily. h When dealing with Runge-Kutta methods that do not satisfy the additional Assumption 2.III, we still have the following result: Lemma 5.6 For k 0 [ 0 taken as in Lemma 5.4, we can bound for all k k 0 for ' p, b 2 f0; 1g, and l ! 0. The constant C depends on M, x, k 0 , l, and the RK method. The estimate (5.13) also holds for b ¼ À1 if the method is stiffly accurate. Additionally Proof The argument to prove (5.13) is very similar to that of Lemma 5.4. By interpolation it is clear that we just need to prove the result for any integer l satisfying 0 l p þ 1 À ' À b. Consider then the rational function q ';b;l ðzÞ :¼ z Àl r ';b ðzÞ, which is bounded at infinity and has all its poles in rðQ À1 Þ [see (5.10)]. We can then use the same argument to prove (5.11) for this redefined new function q ';b;l . (Note that we do not use Assumption 2.III in this argument.) Using that r ';b ðk AÞ ¼ k l q ';b;l ðkAÞA l in dom A l ð Þ, the result follows. Stiff accuracy of the method is used in the case b ¼ À1 when we apply Lemma 5.3, dealing with the zeros of r ';À1 .
The proof of (5.14) is a similar adaptation of the proof of Lemma 5.5. h

Estimates on discrete convolutions
The RK error will naturally induce several types of discrete convolutions that we will need to estimate separately. In all of them we will have the structure We first deal with the simplest cases.
Lemma 5.7 For nk T, the sequence defined by (5.15) can be bounded by kx n k X nkq k ðTÞ max j n kg j k X Tq k ðTÞ max j n kg j k X : If g n :¼ gðkAÞ > n n for n n 2 X m , then kx n k X CTq k ðTÞ max j n kn n k X m :

SN Partial Differential Equations and Applications
Proof Follows by writing the recurrence (5.15) as a discrete convolution.
h The next estimate is related to the consistency error of the RK method in the sense of how the RK method approximates derivatives at the stage level. We introduce the operator D k ðy; tÞ :¼ yðt þ kcÞ À yðtÞ1 À kQ _ yðt þ kcÞ: ð5:16Þ The following well-known result about D k ðy; tÞ underlies the proofs of [  Proof We will use the function e k b defined in (5.19) and Abel's summation by parts: : Since e k b ðtÞ À e k b ðt À kÞ ¼ gðkAÞ > ðkQÞ b D k ðy À yðÁ À kÞ; tÞ; and using that y ðjÞ ðtÞ À y ðjÞ ðt À kÞ X l k max tÀk s tþk y ðjþ1Þ X l , a computation analogous to the above bound, but using y À yðÁ À kÞ as data implies Note that if b 2 f0; 1g we can make a simpler estimate for the term originating from R k , (i.e., the one containing the highest derivative) using less regularity for y by not taking advantage of the difference between y ðpþ1Þ ðt j Þ and y ðpþ1Þ ðt jÀ1 Þ and thus end up requiring less regularity. Using the estimate (5.20) for the last term in (5.21), we have thereby already derived estimates for all three terms in (5.21). h

Proofs
The two different cases (with or without Assumption 2.III) will be collected by using the parameter On the other hand, fX k n g ¼ ðo k Þ À1 fU k n g solves by Proposition 4.6: ð6:3bÞ Before we can estimate the difference between the functions x and x k n , we need one final lemma.
Assume that for some l ! 0 we have x 2 C pþ1 ð½0; T; X l Þ; H 2 C p ð½0; T; X l Þ; EC 2 C p ð½0; T; X l Þ: Proof We set y :¼ x À EC. By assumption we have y 2 C p ð½0; T; X l Þ and B x À EC ð Þ¼0. Since x 2 domðA H Þ and range E & domðA H Þ this implies yðtÞ 2 domðAÞ for all t 2 ½0; T. We further calculate using (6.4) and range E kerðI À AÞ H : Each of the terms is assumed in C p ð½0; T; X l Þ, thus y 2 C p ð½0; T; X lþ1 Þ. h We will need the sequences fc k n g and fh k n g with the scalar parts of the computations of fC k n g and fH k n g respectively, namely (see Lemma 2.6), We then consider D k n :¼ ðI EÞðCðt n þ kcÞ À C k n Þ; d k n :¼ EðCðt n Þ À c k n Þ: Using (6.5a), the definition C k ¼ ðo k Þ À1 N, and (2.13), we can write D k n À 1d k n ¼ ðI EÞCðt n þ kcÞ À 1ECðt n Þ À kQ E _ Cðt n þ kcÞ ¼ D k ðEC; t n Þ: ð6:6Þ Lemma 5.2 (take f ¼ EN for the first three inequalities and f ¼ F for the last one) proves that kd k n k X CTk p kEN ðpÞ k T;0 ; ð6:7aÞ kd k n À d k nÀ1 k X Ck pþ1 kEN ðpÞ k T;0 ; ð6:7bÞ kkb > D k n k X Ck pþ1 ðkEN ðpÀ1Þ k T;0 þ TkEN ðpÞ k T;0 Þ; ð6:7cÞ kHðt n Þ À h k n k X CTk p kF ðpÞ k T;0 : ð6:7dÞ The error analysis is derived by tracking the evolution of the following differences E k n :¼ xðt n þ kcÞ À X k n À D k n 2 ðdomðAÞÞ m ; e k n :¼ xðt n Þ À x k n À d k n ; (compare (6.2) and (6.3b) to see the vanishing boundary condition for E k n ) and note that by (6.7a) kxðt n Þ À x k n k X ke k n k X þ CTk p kEN ðpÞ k T;0 ; which shows that we only need to estimate e k n to prove Theorem 3.4. We start with the observation that x solves the following equation, as can be easily derived from Eq. (6.2): xðt n þ kcÞ ¼ 1xðt n Þ þ kQ A H xðt n þ kcÞ þ xðt n þ kcÞ À kQ _ xðt n þ kcÞ þ kQHðt n þ kcÞ À 1xðt n Þ ¼ 1xðt n Þ þ kQ A H xðt n þ kcÞ þ kQHðt n þ kcÞ þ D k ðx; t n Þ: ð6:8Þ Recalling that Assumption 2.I included the hypothesis range E & kerðI À A H Þ, we have ðQ A H ÞD k n ¼ QD k n . Combining (6.8) and (6.3a), we get Naive estimation of the terms D k ðx; t n Þ and D k n À 1d k n would yield convergence rates similar to Propositions 3.1 and 3.2. In order to get an increased rate, as stated in Theorem 3.4, we combine these two terms using the function YðtÞ :¼ xðtÞ À ECðtÞ. Lemma 6.1 and the assumptions of Theorem 3.4 ensure Y 2 C pþa ð½0; T; X lþ1 Þ \ C pþ1 ð½0; T; XÞ.

Proof of Theorem 3.5
This proof is very similar to the one of Theorem 3.4 but slightly simpler. We will point out the main steps of the proof. We first focus on showing the estimate for v À v k . Note that we use the simple form of o k for stiffly accurate RK methods given in Lemma 2.7. We define and fV k n g ¼ o k fU k n g satisfies (see Proposition 4.5 and Lemma 2.7, where we use stiff accuracy of the RK scheme, and recall that fG k n g ¼ o k fFðt n þ kcÞg and fH k n g ¼ o k fNðt n þ kcÞg) Let then D k n :¼ ðI EÞðH k n À Hðt n þ kcÞÞ ¼ k À1 Q À1 D k ðEN; t n Þ and [note (6.11b) and (6.11c)] E k n :¼ V k n À e V k n À D k n 2 ðdomðAÞÞ m ; e k n :¼ v k n À e v k n : By (6.11a) and (6.12a), using that ðQ A H ÞD k n ¼ QD k n (assumption on the lifting) and Lemma 2.7 to represent G k n , we have kðb > AÞE k n ¼ gðkAÞ > ð1e k n À D k n þ kQD k n þ D k ðF; t n ÞÞ and therefore, from (6.11c) and (6.12c) e k nþ1 ¼ rðkAÞe k n À gðkAÞ > ðkQÞ À1 D k ðEN; t n Þ þ gðkAÞ > D k ðEN þ F; t n Þ þ kb > Q À1 D k ðEN þ F; t n Þ: ð6:13Þ The final term can be shown to be of order Oðk pþ1 Þ by combining ( The estimate involving A H u can be proved as an easy corollary of the estimate on v. Since the last stage of a stiffly accurate method is the step, we have that (4.8a) implies that A H u k n ¼ v k n À Fðt n Þ and therefore A H uðt n Þ À A H u k n ¼ vðt n Þ À v k n :

Maximal dissipative operators in Hilbert space
In this short section we summarize some results that show that the hypotheses on the abstract equation and its discretization are simpler for maximal dissipative operators on Hilbert spaces. These results are well-known and will be needed when applying the theory developed in the previous sections to some model problems in Sect. 8

Applications
In this section X is a bounded Lipschitz open set in R d (d ¼ 2 or 3) with boundary C. We use the usual (fractional) Sobolev spaces H s ðXÞ for s ! 0 and introduce the space H 1 D ðXÞ :¼ fu 2 H 1 ðXÞ : Du 2 L 2 ðXÞg. On the boundary C, we also consider Sobolev spaces H s ðCÞ and their duals H Às ðCÞ. Details can, for example be found in [22].
For the trace operators, we make the convention that the index þ relates to exterior and -means the trace is taken from the interior of X. For example, the two bounded surjective trace operators c AE : H 1 ðR d n CÞ ! H 1=2 ðCÞ denote the trace from R d n X and X, respectively. and we will denote H À1=2 ðCÞ for the dual of the trace space. The angled bracket h Á ; Á i C will be used for the H À1=2 ðCÞ Â H 1=2 ðCÞ duality pairing and ðÁ; ÁÞ R d will be used for the inner product in L 2 ðR d Þ and Â L 2 ðR d Þ Ã d . We will also use the normal traces c AE m : Hðdiv; R d n CÞ ! H À1=2 ðCÞ and the normal derivative operators o AE m . Here we make the convention that the normal derivative points out of X for both interior and exterior trace.
We note that the applications in this section are chosen for their simplicity. More complicated applications, also involving full discretizations by convolution quadrature and boundary elements of systems of time domain boundary integral equations can be found in [29] and [27].

Boundary integral equations and convolution quadrature
In this section, we give a very brief introduction to boundary integral equations and their discretization using convolution quadrature. In that way, we can later easily state our methods for both the heat and wave equations in a concise and unified language. We present the result mostly formally, but note that they can be made rigorous under mild assumptions on the appearing functions. This theory can be found in most monographs on boundary element methods, see e.g. [22,32,33] or [31].
For s 2 C þ , we consider solutions u 2 H 1 ðR d n CÞ to the Helmholtz equation For this problem, the fundamental solution is given by where H ð1Þ 0 denotes the first kind Hankel function of order 0. Using the representation formula, u can be rewritten using only its boundary data: where the single layer and double layer potentials are given by We note that both SðsÞk and DðsÞw solve the Helmholtz equation for any given densities k 2 H À1=2 ðCÞ and w 2 H 1=2 ðCÞ.
We will need the following four boundary integral operators:  [21]). Given a Runge-Kutta method, it is then easy to define the convolution quadrature approximation to such operators, as was introduced in [19]. We just replace the Laplace transform by the Z-transform Z and s with the function d=k, i.e., we define: where g denotes a sequence in the shared domain of F(s) and k [ 0 denotes the step size. The matrix-valued function z7 !Fð dðzÞ k Þ is defined using the Riesz-Dunford calculus, but can be computed in practice by diagonalizing the argument.
Remark 8. 1 We note that our use of the notation o k and ðo k Þ À1 is consistent with this definition by using the functions FðsÞ :¼ s and FðsÞ :¼ s À1 .

An exotic transmission problem
In this section we show how to apply Theorems 3.4 and 3.5 to a transmission problem in free space associated to the infinitesimal generator of a group of isometries (both AEA are maximal dissipative) with some exotic transmission conditions which impose partial observation of a trace. In Sect. 8.3 we will explain how this problem is related to a boundary integral representation of a scattering problem and how the current results yield the analysis of a fully discrete method for that integral representation. We keep the presentation brief. For more details and exemplary applications we refer to [13].
Let Y h be a closed subspace of H 1=2 ðCÞ (in practice it will be finite-dimensional) and consider the spaces Hðdiv; R d n CÞ :¼fw 2 L 2 ðR d n CÞ d : r Á w 2 L 2 ðR d n CÞg; ð8:4aÞ The condition scvt 2 Y h is equivalent to We then set In X we use the natural inner product, in V we use the norm of H 1 ðR d n CÞ Â Hðdiv; R d n CÞ, and in M we use the usual norm. We will define A H : Let k 2 L 2 ðCÞ be given. By applying Theorem A.4 to the exterior and setting e w ¼ 0 inside, we can construct a function e w 2 Hðdiv; R d n CÞ satisfying sc m e wt ¼ k and e w k k Hðdiv;R d nCÞ þ e w k k ½L 2 ðXÞ;H0ðdiv;R d nXÞ 1=2 . k k k L 2 ðCÞ : ð8:9Þ Upon identifying the product of function spaces on X and R d n X with a function space on R d n C, we have The product of interpolation spaces equals the interpolation of product spaces (cf. Lemma A.5); we can therefore also estimate: If we consider ðv; wÞ :¼ Ek, then ðv; w À e wÞ 2 domðAÞ by construction of the lifting. Thus we have The continuity of E from Proposition 8.3 and (8.9) conclude the proof. h Proposition 8.5 If g 2 C 2 ð½0; 1Þ; H À1=2 ðCÞÞ satisfies gð0Þ ¼ _ gð0Þ ¼ 0, then (8.6a-8.6c) has a unique strong solution.
We also need a regularity result that allows us to bound time derivatives of the solution in terms of the data. The continuity condition for the ðs þ 2Þ-nd derivative of g in Proposition 8.6 can be relaxed to local integrability, but then the norms on the right-hand side of (8.10) have to be modified. We now consider the RK approximation of (8.6a-8.6c) in a finite time interval [0, T], which provides pairs of stage values ðV k h;n ; W k h;n Þ 2 X m and step approximations ðv k h;n ; w k h;n Þ 2 X. We then define fU k h;n g ¼ ðo k Þ À1 fV k h;n g; u k h;n ¼ rð1Þu k h;n þ b > Q À1 U k h;n ; n ! 0 ð8:11Þ with u k h;0 ¼ 0 (see Lemma 2.6) and w k h;n :¼ scu k h;n t .
Proposition 8.7 For sufficiently smooth g, with RK approximations using a method satisfying Assumption 2.II, and with a given by (6.1), for nk T we have the estimates follows from Propositions 3.1 and 3.2, using (8.10) for the estimate in terms of the data. The H 1 ðR d n CÞ estimate (8.13) is then a direct consequence of (8.12) and (8.15). The estimate for w h À w k h follows from the standard trace theorem. h

Scattering
As it does not incur much difficulty, we cover both the exterior scattering problem, which is an exterior Neumann problem, as well as the interior Neumann problem. In order to do so, we define the domain X þ :¼ R d n X and X À :¼ X and distinguish the problems by adding the superscripts þ or -to the functions involved. We stay in the geometric setting of the previous section. Assume that d 2 R d is a unit vector (direction of propagation) and that c 2 R is such that Let / : R ! R be a function such that /ðrÞ ¼ 0 for all r ! c: The incident wave u inc ðx; tÞ :¼ /ðx Á d À tÞ propagates in the direction d at unit speed and has not reached the scatterer given by X at time t ¼ 0. The data for our problem will be the function g : ½0; T ! L 2 ðCÞ given by gðtÞ :¼ Ào AE m u inc ðÁ; tÞ. The problems under consideration are: Find u AE : so that o AE m ðu AE þ u inc Þ ¼ 0. (Note that we can take the trace of the normal derivative of the incident wave, since it is locally smooth.) The exterior problem (posed on X þ ) is the classical sound soft scattering problem of the incident wave u inc .
A direct formulation for solving this problem is obtained by extended the function by zero to the complement of the domain of interest. That is, we solve: with u AE ð0Þ ¼ _ u AE ð0Þ ¼ 0: By imposing some additional hypotheses on the growth of g (which is needed to have a well-defined distributional Laplace transform), we can represent the solution to (8.16) as u AE ¼ ÇSðoÞg À DðoÞw AE , where w AE :¼ scu AE t. Note that, to be precise with the use of weak distributional definitions, all functions have to be extended by zero to t\0 (we say that they are causal) and the time interval is extended to infinity.
Taking the trace in this representation formula, the solution of (8.16) can be found by solving an equation for w AE and then postprocessing with the potential operators: ð8:17Þ and we still have that w AE ¼ scu AE t.
For simplicity of notation, we will skip the indices ± for the different functions from now on. We can equivalently write (8.16)  For the discretization, we consider a finite dimensional space Y h and the Galerkin approximation to (8.17), so that we look for w h : R ! X h causal such that ð8:18Þ The functions v h :¼ _ u h and w h :¼ ru h satisfy (8.6a-8.6c). The difference between the solutions of (8.16) and (8.18) can be studied by comparing the solutions to (8.6a-8.6c) when Y h ¼ H 1=2 ðCÞ and when Y h is a finite dimensional space, see [13] for details. For our purposes, it is sufficient to note that we get quasi-optimal estimates for the discretization in space.
Discretization in time is performed by applying convolution quadrature to (8.18). The fully discrete solution reads The approximations w k h and u k h are then computed by the usual post-processing, i.e., Lemma 8.8 The sequences u k h and w k h computed via (8.19) coincide with the Runge-Kutta approximations to (8.6a-8.6c) and their traces respectively.
Proof The details of the computation can be found in the appendix of [23]. The basic idea is to take the Z-transform and show that both approaches solve the matrix-valued Helmholtz problem (4.2a-4.2c). h This gives the following immediate corollary, representing an a priori bound for the fully discrete method: see [20,Appendix 2]. Applying the abstract theory of [2] then implies convergence rate minðq þ 1; pÞ for the boundary data w h . Modifying their proof, one can also get for b g 2 L 2 ðCÞ that which would yield the same convergence rate as Corollary 8.20, but without insight into the dependence on the end-time T. Ä

Numerical example
We solve (8.19) on a ''hollow square'', as depicted in Fig. 1, and focus on the interior Neumann problem, i.e. computing w À ¼: w. The geometry was chosen to be non-convex and not simply connected, in order to test if the rate observed is a general result, or if our estimates might prove sharp in some situation.
We prescribe the exact solution as a traveling wave, given by i.e., we compare to the L 2 -projection of the exact solution. Since we are interested in the convergence rate with respect to the timestep size k, we consider a fixed, but sufficiently fine mesh. We used 3 and 5 stage Radau IIA methods, with orders (q, p) of (3,5) and (5,9), respectively (see [16] for their definition). While their strong damping properties are not advantageous for wave propagation problems, they nevertheless are the standard method used with convolution quadrature. This is in part due to the fact that the standard theory (see, e.g., [2]) makes some assumptions not satisfied by the Gauss methods. A more detailed analysis of the dissipation and dispersion of the Radau methods was performed in [7,Sect. 4.3], showing that higher order Runge-Kutta methods posess favorable properties compared to their low order brethren.
Our theory predicts convergence rates of 4.5 and 6.5. In Fig. 2, we observe a rate that is closer to 5 and 8. This means that (just like the standard Laplace-domain estimates) our estimates do not appear to be sharp in this case. Further investigations into the cause of this phenomenon are required. Results trying to explain this phenomenon, initially prompted by the work on this article, can be found in [24] but with a different model problem.

The heat equation
In this section, as an example where our estimates turn out to be sharp, we consider a heat conduction problem and will apply Theorem 3.5 to get convergence of the boundary trace. The physical situation is a body X & R d that is held at a given temperature distribution and radiates heat into a medium X þ :¼ R d n X. We make the simplifying assumption that at t ¼ 0 the temperature is 0. Since the problem is posed on an unbounded domain, it is a Fig. 2 Performance of Radau IIA methods for the wave equation, cf. Sect. 8.4 good candidate for boundary integral equations, while being simple enough to showcase our more general results. We only briefly give the mathematical setting. More details and a more involved physical example can be found in [27]. The setting is as follows: find u : In order to derive the boundary integral formulation, we take the Laplace transform of (8.21a), giving for j :¼ ffiffi s p : ÀDb uðsÞ þ j 2 b uðsÞ ¼ 0; which is Helmholtz's equation for a complex wave number j. We make an ansatz of the form b u ¼ SðjÞ b k for some unknown density b k, which can be determined by applying the trace operator, giving the equation VðjÞ b k ¼ LðgÞ.
Transforming back, and using the definition V j ðsÞ :¼ Vð ffiffi s p Þ, we get the formulation: The solution u can then be recovered by computing u ¼ S j ðoÞ, where S j ðsÞ :¼ Sð ffiffi s p Þ.
The discrete version of this is then given by solving It can be shown that plugging the discrete solution into the representation formula U k :¼ S j ðo k ÞK k gives back the Runge-Kutta approximation of (8.21a-8.21c). The approximations at the endpoints t n ¼ n k, denoted by k k and u k respectively can be computed by the usual post-processing. We refer to the appendix of [23] for an analogous computation in the context of the Schrödinger equation, which easily transfers to our situation. For simplicity, we do not consider any discretization in space. A Galerkin approach could easily be included into the analysis, analogously to Sect. 8 Since for stiffly accurate RK-methods, u k also satisfies the boundary conditions (it is the just last entry of the stage vector) we get from the Dirichlet-boundary conditions that c þ uðt n Þ ¼ gðt n Þ ¼ c þ u k ðt n Þ. Therefore, integration by parts and the Cauchy-Schwarz inequality give: kru k ðt n Þ À ruðt n Þk 2 L 2 ðR d nXÞ ¼ À À Du k ðt n Þ À Duðt n Þ; u k ðt n Þ À uðt n Þ Á L 2 ðR d nXÞ kDu k ðt n Þ À Duðt n Þk L 2 ðR d nXÞ ku k ðt n Þ À uðt n Þk L 2 ðR d nXÞ : Estimate ( We are still free to choose the precise lifting v. Doing so as in [ The result then follows from the previous estimates. h Remark 8.13 Note that in the cases q ¼ p À 1 with a ¼ 1 and q ¼ p with a ¼ 0, the rates r 1 and r k in Theorem 8.12 are sharp from what can be extracted from Theorem 3.5 and Propositions 3.1 and 3.2. Nevertheless, we expect it to be possible to extract better rates from a more explicit investigation of these limiting cases.

Numerical example
In order to demonstrate that the estimate (8.25) is sharp, we consider a simple model problem. Following [34], we take X to be the unit sphere and consider a right-hand side g(x, t) of the form gðx; tÞ :¼ wðtÞY m n ðxÞ; where Y m n is the spherical harmonic of degree n and order m. It is well-known that the spherical harmonics are eigenfunctions of the pertinent boundary integral operators. Most notably for us, we have where j n denotes the spherical Bessel functions and h ð1Þ n is the spherical Hankel function of the first kind. Due to this relation, solving (8.22) becomes a purely one dimensional problem, i.e., we can write kðx; tÞ ¼ e kðtÞY m n ðxÞ and the solution can be easily computed to very high accuracy. For our experiments we chose n ¼ 2.
We compare the 3-stage and 5-stage Radau IIA methods (see [16] for their definitions). These methods have stage orders 3 and 5 respectively and both are stiffly accurate and satisfy Assumption 2.III. We therefore expect convergence rates for the density k of order 3.5 and 5.5. Since the exact solution is not available, we compute the difference to an approximation with step size k/4 and use this as an approximation to the discretization error. The results can be seen in Fig. 3. We observe that the results are in good agreement with our predictions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Fig. 3 Convergence for the density e k for the heat conduction problem (cf. Sect. 8.5), comparing Radau IIA methods

SN Partial Differential Equations and Applications
Funding Open access funding provided by Austrian Science Fund (FWF).

A Interpolation of Sobolev spaces
In this appendix we prove that in Lipschitz domains and for certain parameters l, the spaces ½L 2 ðXÞ; H 1 0 ðXÞ l contain functions with non-vanishing boundary conditions. Such estimates are the main ingredient when determining the convergence rate of Runge-Kutta methods using the theory developed in the previous sections. For l\1=2, it is well-known that the fractional Sobolev spaces H l ðXÞ ¼ ½L 2 ðXÞ; H 1 ðXÞ l;2 and e H l ðXÞ ¼ ½L 2 ðXÞ; H 1 0 ðXÞ l;2 coincide (see e.g. [22,Theorem 3.40] together with the results in [22, Appendix B] to identify the Sobolev spaces with the interpolation space). We prove that when interpolating using the index 1, the critical value l ¼ 1=2 is also admissible, provided that some further regularity is provided.
In order to state our result, we need additional notation, notably we define interpolation spaces for q 2 ½1; 1Þ as When considering the Neumann problem in Sect. 8.3, we need to devise a lifting to a vector field with a given normal jump in L 2 . In general, such liftings do not have B Proof We focus on a single chart in the parametrization of (a vicinity of) C. Let O X and D R dÀ1 be open, r 2 R n , u : D ! R, y 0 : D ! R such that we can write X t \ O ¼ È ðx; uðxÞ þ yrÞ : x 2 D; and y 2 ð0; y 0 ðxÞÞ É : By the Lipschitz assumption, we note that y 0 ðxÞ.Ct. Following the considerations in [ with an implied constant depending only on X.
Proof For simplicity, assume that X is bounded. By performing an appropriate cutoff away from oX, all arguments can be localized.
Step 1: Consider the case R C g ¼ 0. Let u be the solution of the Neumann problem Du ¼ 0 in X; o m u ¼ g on oX; Z X u ¼ 0: In addition to u 2 H 1 ðXÞ, by [17] (see also [10,Theorem A.6]), this harmonic function u also satisfies NðruÞ k k L 2 ðCÞ g k k L 2 ðCÞ : For fixed t [ 0 we again select a smooth cutoff function v t satisfying (A.4). We set w :¼ ru and calculate using Lemma A.3: w k k ½L 2 ðXÞ;H 0 ðdiv;XÞ 1=2;1 . esssup t ! 0 t À1=2 Â ð1 À v t Þw k k L 2 ðXÞ þt v t w k k Hðdiv;XÞ Ã . esssup t ! 0 t À1=2 kruk L 2 ðX 2t Þ . NðruÞ k k L 2 ðCÞ . g k k L 2 ðCÞ : Step 2: In the case R C g 6 ¼ 0, the harmonic Neumann problem does not have a solution. Instead, we define u as the solution to ÀDu þ u ¼ 0 in X; o m u ¼ g on oX; and again set w :¼ ru. By construction we have r Á w ¼ Àu 2 H 1 ðXÞ. We decompose u ¼ u 0 þ u 1 , where u 0 solves the full-space problem where u was extended by 0 outside of X. As u 2 L 2 ðR d Þ, standard regularity theory gives u 0 2 H 2 ðBÞ on any ball B and in particular u 0 2 H 2 ðXÞ. In turn this yields o m u 0 2 L 2 ðoXÞ. By construction u 1 then solves As g, o m u 0 2 L 2 ðoXÞ and R C g À o m u 0 ¼ 0, we can apply Step 1 to get ru 1 2 Â L 2 ðXÞ; H 0 ðdiv; XÞ Ã 1=2;1 : Since ru 0 2 ðH 1 ðXÞÞ d ðB to conclude the proof of (A.7). h The following lemma appears to be known in the community, see e.g. [5, Section 3.13, Exercise 4], but in order to be able to rigorously cite, we provide a short proof.
Lemma A.5 Let X :¼ ðX 1 ; . . .; X N Þ and Y :¼ ðY 1 ; . . .; Y N Þ, where X j ; Y j are Banach spaces with continuous embedding Y j & X j , and the product space carries any l p -norm. Fix q 2 ½1; 1 and h 2 ð0; 1Þ: Then, the product of the interpolation spaces coincides with the interpolation of the product spaces. Namely, the following estimate holds for all x :¼ ðx 1 ; . . .; x N Þ 2 ½X; Y h;q : Proof For j 2 f1; . . .; Ng, consider the operators S j and T j defined as SN Partial Differential Equations and Applications T j : X ! X j ðx 1 ; . . .; x N Þ7 !x j and S j : X j ! X x j 7 !ð0; . . .; x j ; . . .0Þ: It is easy to see using the interpolation estimate (2.6) that We therefore calculate: For the opposite direction, we observe that