{'fraction_non_alphanumeric': 0.057160465295431784, 'fraction_numerical': 0.020530005392496727, 'mean_word_length': 3.9097957639939485, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Boundary integral solvers for an evolutionary exterior Stokes problem

This paper proposes and analyzes a full discretization of the exterior transient Stokes problem with Dirichlet boundary conditions. The method is based on a single layer boundary integral representation, using Galerkin semidiscretization in the space variables and multistep Convolution Quadrature in time. Convergence estimates are based on a Laplace domain analysis, which translates into a detailed study of the exterior Brinkman problem. Some numerical experiments are provided.


Introduction
In this paper we propose a fully discrete method based on an integral equation for the exterior Stokes problem with Dirichlet boundary conditions in two or three dimensions. The integral equation is based on a single layer potential representation of the velocity field. The numerical discretization uses a general Galerkin scheme for semidiscretization in space and Convolution Quadrature [17] for discretization in time. The analysis is carried out by combining ideas of Laplace domain analysis of integral operators [18] with the transformation of the Galerkin-BEM discretization in space into a set of exotic transmission conditions [16]. As part of the paper, we include a novel analysis of the single layer potential and operator for the Stokes resolvent equations (the Brinkman equations) on a general Lipschitz domain.
The literature on numerical methods for integral representations of parabolic problems has focused extensively on the heat equation. Most theoretical results are based on the single-layer representation, leading to a Volterrà-Fredholm integral equation that can be formally considered to be of the first kind. (We note that the mapping properties of the integral operators make the integral equations of the second kind for parabolic problems not to be a smooth perturbation of the identity, due to the mapping properties in the time variable. Additional complications arise when the boundary is not smooth.) This analysis was sparked by the work of Arnold and Noon [3] and Costabel [7], with some sequels as [12]. The work of Lubich and Schneider [18] offered a numerical treatment of the heat equation single layer operator equation. Other formulations, including fast multiplication techniques, appear in recent work of Tausch [22,23,19]. The mathematical literature for the unsteady exterior Stokes problem using integral equations seems to be quite limited: see, for instance, [10], [11]. A general overview of the state of the art of time domain integral equations one decade ago can be found in [6].
For our analysis we will rely on properties of the Brinkman single layer potential. We will however take a different approach than the one given in [13,14,15], since we need to study the behavior of all the bounds as functions of the parameter in the Brinkman model. We will adopt a Laplace domain approach similar to the one used in [4] for the wave equation. For some technical issues, we will rely on recent results on the Stokes potentials on general Lipschitz domains [20]. The passage to the time-domain will be done with a modification of a result in [18]. Following [16] we will analyze the semidiscretization in space in a systematic way, showing that a postprocessed solution (the velocity field) can have better properties than the preprocessed solution (the boundary density and, therefore, the pressure field, which is postprocessed with a steady-state operator). Finally, we will apply a general multistep-based Convolution Quadrature strategy and analyze it using the results in [17]. We note that this final step will be the only one where we will not be able to analyze how the constants that appear in the error estimates depend on time (as the latter grows to infinity).
The paper starts with two long sections (Sections 2 and 3) presenting the integral and variational forms of the single-layer potential for the Brinkman problem and providing bounds in terms of the parameter of the Brinkman equation. In Section 4 we transfer the Laplace domain estimates to estimates for the transient single layer potential for the Stokes equation, using a technical result that is proved in Appendix A. In Section 5 we introduce and analyze a general Galerkin semidiscretization in space of the integral equation. We provide bounds for the semidiscretization in space (Galerkin error operator) plus some stability bounds (Galerkin solver) that are needed for the analysis of the fully discrete method. In Section 6 we present and analyze the fully discrete scheme and show some numerical experiments. Finally, Appendix B shows an alternative integral formulation that can be used to eliminate some inconvenient Lagrange multipliers that are needed to impose conformity restrictions in the boundary element space. Foreword on background and notation. We will use basic properties of Sobolev spaces on bounded Lipschitz domains and on their boundaries [1]. All aspects related to integral operators will be proved using variational techniques [20]. The passage to the time-domain requires the momentary use of basic vector-valued distribution theory. It is important to remark that all brackets will be taken to be bilinear, even if they are employed in the context of complex-valued functions. In particular, for scalar fields complex-valued u, v, vector fields u, v and matrix-valued fields (tensors) U, V, and an open set O ⊂ R d , Given a Hilbert space X, we will write X := X d and immediately assume it to be endowed with the product topology.

The Brinkman single layer potential
In this section we present the variational theory for the Brinkman single layer potential as a holomorphic function of its parameter. This is equivalent to studying the single layer potential associated to the resolvent Stokes problem. In all the following arguments, the parameter s is a complex number not in the negative real axis The space of solenoidal vector fields will also play a key role. The geometric setting is as follows: we consider a bounded Lipschitz domain Ω − , with connected boundary Γ, and the associated unbounded exterior domain Ω + := R d \ Ω − . The superindices ± will be used to refer to limits/traces from Ω ± . We will use the angled bracket · , · Γ to denote the L 2 (Γ) and L 2 (Γ) inner products (with the above convention on not conjugating any component) as well as its extension to duality products between the spaces H ±1/2 (Γ), as well as between their vector-valued counterparts.
Jumps of traces and normal stresses The jump of the trace across Γ, for a locally H 1 function, is defined as be such that , ε(u) := 1 2 (Du + (Du) ). We can then define the functionals t ± (u, p) ∈ H −1/2 (Γ) given by the relations: We can thus define the jump of the normal stress [[t(u, p)]] := t − (u, p) − t + (u, p). In spite of the global definitions of t ± (u, p) using test functions in H 1 (Ω ± ), it is clear that these operators have a local behavior and can be extended to pairs (u, p) that only exhibit the required properties in a neighborhood of the boundary. This subtle distinction will be employed in the two dimensional case, where the pressure p is only locally in L 2 . In particular we will use the formula where a(u, v) := 2ν (ε(u), ε(v)) R d and D(R d ) is the set of C ∞ compactly supported functions.

Integral forms
The following definitions can be found in [14], [15], [13, p.81] The pressure potential. For a given density λ ∈ H −1/2 (Γ), we define where e p (r) := 1 2(d − 1)π 1 r d r is the negative gradient of the fundamental solution to the Laplace equation. The behavior at infinity of e p gives different properties for the operator S p in two and three dimensions. In the two dimensional case, the closed subspace plays an important role. As a simple fact that this pressure part of the single layer potential is the same for the Brinkman as for the Stokes problems we can show the following result [20, Propositions 5.2 and 7.2]: First order asymptotics in the two dimensional case. We note that the zero integral condition in the definition of the space H −1/2 0 (Γ) only affects the behavior at infinity of S p λ. We therefore explore the first order asymptotics at infinity of S p λ for general λ. Expanding the kernel function e p , we can write (S p λ)(z) = 1 2π as |z| → ∞.
Therefore, apart from the leading term, the pressure potential is in L 2 (R 2 ). Let then Here {e 1 , e 2 } is the canonical basis for R 2 . This leads to the proof of the following result.
is bounded. We note that while the result in Proposition 2.1 is a direct consequence of what is known for the Stokes operator (see [20] for a fully developed variational and integral theory), the decomposition of Proposition 2.2 (that subtracts an easily identifiable first order term from the potential) seems to be new.
The velocity potential. For λ ∈ H −1/2 (Γ), we define and K being the modified Bessel function of order . The square root in (2.6) is the one determination of the square root that is analytic in C . Note that A 3 and B 3 are entire functions with A 3 (0) = B 3 (0) = 1. The functions A 2 and B 2 are only analytic in C , and have logarithmic singularities in the cut (−∞, 0]. Comparing the integral expressions of the Brinkman potential with those of the Stokes potential, it is possible to prove that for any s ∈ C the operator S u (s) : We are however interested in the dependence on s of the bounds for this potential and some related integral operators.

Variational theory in three dimensions
Proposition 2.3 (Existence via potential theory). Let λ ∈ H −1/2 (Γ) and consider the functions u λ := S u (s)λ ∈ H 1 (R 3 ) and p λ : Proof. Regularity is guaranteed by the properties of the integral formulations of the layer potentials. The differential equations (2.7a)-(2.7b) are satisfied pointwise in a strong sense (this can be proved by differentiation directly in the fundamental solutions), and therefore, they are satisfied in a distributional sense. Condition (2.7c) is a direct consequence of the fact that u λ ∈ H 1 (R 3 ). Finally, condition (2.7d) follows from (2.7a) and (2.2). The equivalence of (2.7) and (2.8) is straightforward. (2.9) Problem (2.9) is well posed.
Proof. Taking v ∈ V(R 3 ) as test functions in (2.8) it is clear that u λ satisfies the equations (2.9). To prove well posedness, note that Applying a density argument, the bilinear form of (2.9) is shown to be coercive.
Proof. Uniqueness. If (u, p) is a solution of the corresponding homogeneous problem, by Proposition 2.4 it follows that u ≡ 0. Therefore p ∈ L 2 (R 3 ) satisfies ∇p = 0, which implies that p ≡ 0. Existence. Proposition 2.3 shows existence of solution using the integral form of the potentials.
Proof. It is clear that (0, −χ Ω − ) is a solution of (2.8) with λ = n. By uniqueness, this is the layer potential.

Variational theory in two dimensions
At this stage, the main difference between the two and three dimensional cases arises from the fact that [20,Proposition 3.2] and it is somehow due to the fact that constant functions are elements of the associated weighted Sobolev spaces. This is not the case for the Brinkman problem. The following approach uses the precise knowledge of the asymptotic behavior of the single layer potential for the pressure variable.
Proposition 2.7 (Existence via potential theory). If λ ∈ H −1/2 (Γ) and we consider the Proof. This is just a consequence of the properties of the associated integral operators.

Bounds in the Laplace domain
In this section we study properties of the operators S(s) := S u (s) (see (2.5)), V(s) := γS(s), V(s) −1 and S(s)V(s) −1 (s) as functions of s ∈ C . We start with two technical results.
Proof. In [20,Proposition 4.4] there is a right inverse whose range contains only compactly supported functions. The same right inverse is valid now. Then the decompositions are stable and there exists C > 0 such that Some technicalities Given s ∈ C = C \ (−∞, 0], we take its square root s 1/2 := |s| 1/2 exp( ı 2 Arg s) ∈ C + and denote We also consider the norms (depending on |s|) Note that The norm ||| · ||| (1) will be used as the standard norm in H 1 (R d ). Note finally that Proof. Let u λ := S(s)λ and u µ := S(s)µ and note that V(s)µ = γu µ . Taking v = u µ as test function in Propositions 2.4 and 2.8, it follows that which proves (a). Let now u λ := S(s)λ. By (a) and (3.5), we can write and therefore S(s)n = 0. If V(s)λ = 0, then, by (b) it follows that u λ = 0. Therefore, the associated pressure p λ = S p λ satisfies ∇p λ = 0 in R d \ Γ, and is decaying at infinity.
2) and the trace theorem.
Using Lemma 3.2 it follows that Therefore by (b), The remainder of the proof is straightforward, using (3.3) to get the final lower bound.
Remark. As part of the proof of Proposition 3.3(b) we have shown that Proposition 3.4 (Bounds for the single layer potential). There exists C > 0 such that Proof. Let u λ := S(s)λ. Then: by the trace theorem Proof. Note first that γu = φ. By Propositions 2.4 and 2.8, Therefore, using the lifting operator in (3.6) and defining u 0 := u − γ † φ, it follows that by (3.4). This implies that The remainder of the proof is straightforward.
Summary of bounds. In terms of ω (see (3.1)), α 1 (s) and α 2 (s) (see (3.3)), we can write The operator norms above are the natural ones using the spaces H 1 (R d ), H 1/2 (Γ), and H −1/2 (Γ), where it corresponds. Using (3.3), to bound we can obtain a new set of bounds, valid for all s ∈ C :

An evolutionary exterior Stokes problem
Given a Banach space X, we consider the set of causal C k X-valued functions The space of bounded linear operators from X to Y (for Hilbert spaces X and Y ) will be denoted B(X, Y ).) Abstract setting. The starting point is an operator valued holomorphic function F : In particular, there exists an B(X, Y )-valued casual distribution f , whose Laplace transform is F(s). The following result is based on [18, Lemma 2.2]. Its proof is given in an appendix. We note that in comparison with [18] we are more demanding in terms of regularity of g, but we pay attention to behavior of constants as t grows.
Proposition 4.1. Let f be such that its Laplace transform F : C → B(X, Y ) satisfies (4.1)-(4.2), and let g ∈ C 1 + (R; X). Then f * g ∈ C + (R; Y ) and where C F satisfies (4.2). Then, for all g ∈ C 2 + (R; X), we have that f * g ∈ C + (R; Y ) and Proof. Let ∂ −1 t f be the distribution whose transform is s −1 F(s). Then f * g = ∂ −1 t f * g and we can apply Proposition 4.1 to ∂ −1 t f and g .

Estimates for the single layer potential and operator
Proof. The distribution S satisfies the hypotheses of Proposition 4.1 (i.e. (4.1)-(4.2)) with µ = 1/2 and = 3. The result is then a direct consequence of Proposition 4.1. (Γ) such that V * λ = φ. Moreover λ ∈ C + (R; H −1/2 (Γ) and the associated potential u = S * λ is also continuous as a function of t. Finally, we have the bounds: Proof. There is a slightly delicate argument to show uniqueness. By causality, we can look at the equation assuming that φ is compactly supported. This means that φ has a Laplace transform Φ, and using the Payley-Wiener theorem and the bounds (3.9), there is a unique solution whose Laplace transform is V(s) −1 Φ(s).
Once existence and uniqueness is settled, the bounds follow from Corollary 4.2. For V −1 we have the hypotheses of Corollary 4.2 with µ = 1/2 and = 3. For S * V −1 , we have the hypotheses with µ = 0 and = 3.

The exterior Dirichlet problem
Our starting point is the velocity field on Γ at all times φ ∈ C 2 + (R; H We finally construct the pressure field, by applying the (time-independent) pressure part of the single layer operator for the steady-state Stokes equation: In the remainder of this section, it is necessary to clarify that all differential operators in the space variables will be used in the sense of distributions in R d \ Γ. If ∂ x i is the differentiation operator with respect to the i-th variable, it is well known that ∂ x i : For any t ≥ 0, the equation (4.7a) is to be understood in the sense of distributions in Proof. Note first that (4.7b) is satisfied because u is a continuous function with values in the space of solenoidal fields V(R d ). The initial condition (4.7d) is a consequence of the fact that u is continuous and causal. By causality, we can assume that φ (k) (t) is bounded in t for k ≤ 2 (this does not affect the generality of the result), and therefore, the Laplace transforms of φ, u and p exist for s ∈ C . Moreover, they satisfy and therefore sU(s) − ν∆U(s) + ∇P (s) = 0 ∀s ∈ C (4.9) and γU(s) = Φ(s) ∀s ∈ C . (4.10) The equality (4.10) proves the boundary condition (4.7c) in the time domain. Note now that ∆u ∈ C + (R; H −1 (R d \ Γ)). In the three dimensional case, it is clear from (4.6) that ∇p ∈ C + (R; H −1 (R 3 \ Γ)). In the two dimensional case, we have to use the decomposition of Proposition 2.2 and the fact that ∇(p ∞ · e ) ∈ L 2 (R 2 ) in order to prove that ∇p ∈ C + (R; H −1 (R 2 \ Γ)). In third place u ∈ C + (R; H −1 (R d \ Γ) and thereforė u is a causal distribution with values in H −1 (R d \ Γ). Taking Laplace transforms of ∆u, u and ∇p and using (4.9) we show thaṫ u = ν∆u − ∇p. (4.11) This equation is to be understood in the sense of causal distributions with values in H −1 (R 3 \Γ). However, as we have seen above, the right hand side of (4.11) is a continuous causal H −1 (R d \ Γ)-valued function. This proves (4.8) and the equality (4.11) is satisfied pointwise in time, that is, we have proved (4.7a) as equality of elements of H −1 (R d \ Γ) for all t. In its turn, this can be understood as a distributional equation in R d \ Γ for all t.

Galerkin semidiscretization in space
Γ) be a finite dimensional space. The semidiscretized BIE for the exterior Dirichlet problem starts with causal Dirichlet data φ : R → H 1/2 n (Γ), looks for a causal function λ h : R → X h (i.e. λ h (t) = 0 for all t ≤ 0) such that and finally constructs The semi-discretized integral equation (5.1) can be also written in the following abstract

The Galerkin solver
We first study properties of the Galerkin solver, i.e., the operator G h (s) : We will also be interested in the associated velocity field u h = S(s)λ h = S(s)G h (s)φ. Note that the space X h was chosen to work in the time domain and can be taken to be realvalued. In the context of Laplace transforms, it has to be closed by conjugation, which is equivalent to taking the same space with complex scalars to create linear combinations.
Proposition 5.1 (Bound for the Galerkin solver). There exists a constant independent of h such that Proof. The first estimate is a direct consequence of the coercivity estimate of Proposition 3.3(d). To show the second one, we need to replicate the proof of Proposition 3.5. Set λ h = G h (s)φ. By Propositions 2.4 and 2.8 (solenoidal variational form for the single layer potential) Consider now the closed space (5.5) Using the lifting (3.6), we can proceed as in the proof of Proposition 3.5 and decompose

and finally show that
from where the bound follows.

The Galerkin error operator
The Galerkin projector looks at the discrete problem from the point of view of the exact solution. We can define it as the operator G h (s)V(s) : H −1/2 m (Γ) → X h , or equivalently by setting λ h as the solution of the discrete equations Instead of studying this projection we will study the complementary projection, that corresponds to the error of the Galerkin semidiscretization. We thus consider the operators E h (s) := G h (s)V(s) − I and S(s)E h (s). Note that while G h (s)V(s) is a projection onto X h for all s ∈ C , the range of E h (s) varies with s.
Proposition 5.2 (Bounds for the Galerkin error operator). There exists a constant independent of h such that where λ h is the solution of (5.6). Proceeding as in the proof of Proposition 3.4 we can show that = −Re s 1/2 λ, V(s)(λ h − λ) Γ by Galerkin orthogonality (1) by the trace theorem.
Therefore, by (3.2) it follows that The second bound gives the estimate for S(s)E h (s) . Using now (3.8), we can bound which finishes the proof.
This decomposition will be used to derive error estimates.

Bounds in the time domain
Proposition 5.3 (Stability; bounds with respect to data). Let φ ∈ C 2 + (R; H 1/2 n (Γ)) and let λ h be the solution of (5.1) and u h be given by (5.2). Then and Proof. It is a direct consequence of Corollary 4.2 (abstract result in the time domain) and Proposition 5.1 (Laplace domain bounds for the Galerkin solver).
For the error estimates we use the orthogonal projection operator Π h : H −1/2 m (Γ) → X h . We first give an estimate of the velocity field, which requires much less regularity in time. We will next give an estimate for the density, which will in turn give an estimate for the pressure field.
Proposition 5.4 (Error estimate for the velocity field). Assume that λ ∈ C 1 Proof. This is the time-domain version of the second bound of Proposition 5.2, using Proposition 4.1 (with µ = 1/2 and = 2) and the identity (5.7) in order to introduce the orthogonal projector.

Full discretization and numerical experiments
We finally do a full discretization of equations (5.1) and (5.2) using Lubich's multistepbased Convolution Quadrature [17]. We next give a short introduction to this black-box technology applied to our particular problem. More implementation details can be found in [5] and [9] (although for wave propagation problems). Before we introduce the method, let us also mention that there is a faster version (which changes the implementation, but not the method itself) called the fast and oblivious CQ method [21] that we will not deal with in this paper. Let us choose a basis {µ j : j = 1, . . . , N } for X h , a time-step κ > 0, and let us consider the uniform grid in time t n := nκ, for n ≥ 0. The data are sampled in time and tested to define vectors The transfer operator corresponding to the convolution with V is defined as a matrixvalued function of s ∈ C : The CQ discretization of (5.1) starts with a Taylor expansion The function δ(ζ) is the characteristic function of the BDF method of order p. Lubich's theoretical results hold for p ≤ 6. Note that most of the matrices V n (κ) ∈ R N ×N do not have to be computed in the practical implementation of the method. The discretization of (5.1) looks for the sequence of vectors λ n ∈ R N given by the recurrence: If λ n = (λ n,1 , . . . , λ n,N ), we then reconstruct the discrete function λ h n := N j=1 λ j,n µ j ∈ X h . The discrete densities provide the discrete pressure field p h n := S p λ h n . (6.2) To compute the discrete velocity field we use another postprocessing of the discrete densities The convergence result follows from [17, Theorem 5.1] by using Proposition 5.1.
Proposition 6.1. Let p be the order of the BDF method used for the CQ discretization. Assume that φ ∈ C p+1

5)
where B = R 3 or B is any bounded open set in R 2 . The constants C 1 and C 2 depend on t, and C 2 depends on B in the two-dimensional case. For small t, C 1 ≤ Ct and C 2 ≤ Ct 1/2 .
A first numerical experiment. In order to be able to compare our method with an exact solution we will solve problem (4.7) in the domain Ω − = (−1, 1) 2 . We choose the data so that the exact solution is where H is the Heaviside function. The exact density λ(t) is not known. Note that even if the exact solution is smooth, there is no guarantee that λ(t) will be a smooth function in the space variable. We integrate from t = 0 to t = 1. For discretization in space we choose a uniform partition of Γ = ∂Ω − in N equally sized elements {e 1 , . . . , e N } where N is a multiple of four. We then consider the spaces where P 1 is the space of polynomials of degree less than or equal to one. Instead of building a basis for X h , we will enforce densities to be in X h using two Lagrange multipliers. This only affects the matrix V 0 (κ) in (6.1). Time discretization is carried out with CQ using BDF(3) as ODE solver in the background, using M time-steps to reach t = 1. We then compute errors for the pressure and the velocity If λ were smooth as a function of the space variable (which we do not know), the expected convergence order predicted by the theory would be h 2.5 + κ 3 , where h = 1/N and κ = 1/M . This does not take into account the possible regularization effects of the potentials. We note that, to the best of our knowledge, there is no theory of time-domain integral equations that is able to predict higher order convergence in weaker norms. The results are shown in A second experiment. We deal with the same exact solution but now use the unit circle as the domain. We measure the same errors, based now on three observation points placed at (0, 0), (1/2, 1/2) and (−.6, .1). We used BDF(3) as the time stepping method, taking M time steps to reach t = 1. For space discretization we use piecewise constant functions on a uniform grid (in parameter space), with N elements, and reduced integration. The fully discrete method that we obtain is equivalent to a Nyström method of the class given in [8]. Because of the smoothness of the domain, the density is a smooth function of the space variables and it is to be expected that order three convergence can be observed for potential postprocessings, although this has never been proved for problems in the time domain. The results are reported in Table 2 1636e-06 2.99 6.1078e-06 3.02 320 320 2.7053e-07 3.00 7.5887e-07 3.01 640 640 3.3819e-08 3.00 9.4578e-08 3.00 Table 2: Results at time t = 1 measured on three points interior to the unit circle. Timestepping is carried out with BDF(3)-based CQ. Piecewise constant functions with reduced integration are used in the space variable.

B An equivalent integral equation
We start by describing the formulation (at the continuous and semidiscrete level) for the Brinkman equation. The aim of this formulation is to incorportate the restrictions for test and trial functions to be in H