Convergence analysis of the scaled boundary finite element method for the Laplace equation

The scaled boundary finite element method (SBFEM) is a relatively recent boundary element method that allows the approximation of solutions to PDEs without the need of a fundamental solution. A theoretical framework for the convergence analysis of SBFEM is proposed here. This is achieved by defining a space of semi-discrete functions and constructing an interpolation operator onto this space. We prove error estimates for this interpolation operator and show that optimal convergence to the solution can be obtained in SBFEM. These theoretical results are backed by a numerical example.


Introduction
The scaled boundary finite element method (SBFEM), first proposed by Song and Wolf, is a boundary element method that does not require a fundamental solution. It has proven to be particularly effective for problems with singularities or posed over unbounded media, see [6,7,9]. In SBFEM, a semi-analytical (or semidiscrete, as we shall call it here) solution to a PDE is constructed by transforming the weak formulation of the PDE into an ODE. Essentially, given a star-shaped domain Ω ⊂ R n , a coordinate transformation is performed (the scaled boundary transformation) in terms of a radial variable and n − 1 circumferential variables. Then, an approximate solution is sought in a space of functions discretized solely in the circumferential direction. The resulting weak formulation posed over this space is then transformed into an ODE which, under certain circumstances, can be solved exactly, yielding a semi-analytical approximation of the solution to the PDE.
The SBFEM has been applied to a wide range of problems that arise in science and engineering, such as crack propagation [8] and acoustic-structure interactions [3]. Moreover, the limitation to star-shaped domains has been overcome with the development of scaled-boundary polygon elements, in which the domain is broken into arbitrarily shaped polygons and shape functions are constructed over these polygons based on SBFEM [5,4,1].
The first author gratefully acknowledges support by the German Research Foundation (DFG) in the Priority Programme SPP 1748 Reliable simulation techniques in solid mechanics under grant number BE6511/1-1. The second author is a member of the INdAM Research group GNCS and his research is partially supported by IMATI/CNR and by PRIN/MIUR. The first and third author gratefully acknowledge the support by Mercator Research Center Ruhr (MERCUR) under grant Pr-2017-0017. We would also like to thank the project partners Prof. Carolin Birk (Universität Duisburg-Essen, Germany) and Prof. Christian Meyer (TU Dortmund, Germany) as well as Professor Gerhard Starke for the fruitful discussions. The objective of this paper is to introduce a rigorous framework in which the error of the approximate solution obtained by SBFEM can be estimated. In particular, the notion of a semi-discrete solution to a PDE is formalized by defining a space of semi-discrete functions and constructing an interpolation operator onto this space. Then, given a semi-analytical solution obtained in the framework of SBFEM, estimates of its error can be obtained by bounding the interpolation operator's error using Céa's lemma. We limit the analysis to Poisson's equation posed on a circular domain for simplicity; this setting is appropriate to highlight the main features of our theoretical setting.
The overview of this paper is as follows: in Section 2 we describe the continuous problem together with the polar coordinate change of variables. In Section 3 we introduce a semi-discretization of our problem, where the domain is discretized only in the angular coordinate. It is shown that the semi-discrete solution converges optimally to the continuous solution. Section 4, making use of the semidiscretization, transforms the original problem into an ODE. Finally, a numerical result reported in Section 5 shows that the method is performing optimally also in the presence of singularities.

Setting of the problem
Given an angle Θ in (0, 2π), we are considering the Poisson problem in the following circular sector (see Figure 1): Since we are going to consider a change of variables when defining the scaled boundary method, we denote with• (with the hat symbol) quantities defined on Ω that correspond to quantities • defined on the reference domain. Hence our problem reads: findû : Ω → R such that wheref ∈ L 2 (Ω) and∆ = ∂ 2 x + ∂ 2 y is the Laplace operator in the Cartesian coordinates (x, y).
Let the curved part of the boundary of Ω be parametrized by the graph (x b (θ), y b (θ)) = (cos θ, sin θ) θ ∈ (0, Θ) and define the open rectangle Q := (0, 1) × (0, Θ). We consider the mapping F : Q → R 2 given by In this particular case, the scaled boundary transformation is given by the change of variables (r, θ) = F −1 (x, y) for (x, y) ∈ Ω, i.e., by the polar coordinate transformation. The Jacobian of F is given by Let u(r, θ) =û(F (r, θ)), then the relation between the gradient in Cartesian coor-dinates∇ = (∂ x , ∂ y ) ⊤ and the gradient in polar coordinates ∇ = (∂ r , ∂ θ ) ⊤ is given by∇û Moreover, the solutionû of (1) satisfies In order to consider the variational formulation of (1) in polar coordinates we have to consider appropriate weighted functional spaces. While this is pretty straightforward and well understood, we present the procedure in detail since the notation will be useful for the analysis of the numerical approximation.
Given a weight function w(r, θ) in Q, we define the weighted Lebesgue space We will use in particular w = r and w = 1/r; it is not difficult to see that we have ||u|| L 2 r (Q) ≤ ||u|| L 2 (Q) ≤ ||u|| L 2 1/r (Q) for all u ∈ L 2 1/r (Q). Furthermore, these spaces are complete [2].
The bound (3) motivates the definition of the following weighted Sobolev spacẽ The following lemma shows that H 1 (Ω) andH 1 (Q) are isometric.
Proof. Letû ∈ H 1 (Ω) and, for 0 < ρ < 1, let B ρ be the ball of radius ρ centred at the origin and B c ρ its complement. For Ω ρ = Ω ∩ B c ρ , the map F : Q ρ → Ω ρ with Q ρ := (ρ, 1) × (0, Θ) is a bi-Lipschitz map, i.e. there exist two constants C 1 , C 2 > 0 such that holds for all (r 1 , θ 1 ), (r 2 , θ 2 ) ∈ Q ρ . Indeed, by the mean value theorem we have and clearly ∇F ∞ ≤ 1. In the same way, F −1 : Ω ρ → Q ρ is a smooth, bijective map and F −1 ∞ ≤ 1/ρ; hence it is Lipschitz continuous and it follows that F and F −1 are bi-Lipschitz when restricted to Q ρ and Ω ρ respectively. As a result of [10, Theorem 2.2.2.], u = Φ(û) is weakly differentiable on Q ρ and the chain rule holds. For n ∈ N, define u n = Φ| Ω 1/n (û) on Q by extendingû by zero outside Ω ρ . For any 0 < ρ < 1 one has that so u n and its derivatives belong to the associated weighted Lebesgue spaces. As a result of the monotone convergence theorem we have that u ∈H 1 (Q) and Repeating the steps above, we can also show that for It is then natural to define the following space in order to take into account the boundary conditionsH 1 0 (Q) := Φ(H 1 0 (Ω)). We are now ready to state the variational formulation of (1) in both coordinate systems.

Definition 2 (Weak form of the Poisson problem in Cartesian coordinates
.

Definition 3 (Weak form of the Poisson problem in polar coordinates
It is well-known that (5) is well posed and so is (6) thanks to the properties of the map Φ and of the isometry shown above.

The semi-discrete Poisson equation
The discretization of (1) with the scaled boundary finite element method is based on a spatial semi-discretization that is described in this section.
We introduce a partition of the parametrized boundary θ → (cos θ, sin θ) given by and consider a finite dimensional approximation of H 1 (0, Θ) generated by a basis at this point is arbitrary. It could be based on finite elements, splines, global Lagrange polynomials, etc.
Due to our choice of boundary conditions, we could also have defined the basis , but we prefer to avoid this in order to allow our analysis to be extended more easily to more general boundary conditions or to a situation where Θ = 2π.
The main idea behind the semi-discretization is to consider families of functions where the variables r and θ are separated formally as follows: Ideally, we would like to have u i (r) = u(r, θ i ) and this choice will be used later on in Subsection 3.1 for the error analysis; it will lead to the analogous of the interpolation operator for standard finite elements. In order to do so, we need to give sense to the radial trace u(r, θ i ). For the sake of readability, we now introduce an abstract setting and we postpone the actual definition of the involved functional spaces to Subsection 3.1. Ultimately, we want to define a semi-discrete space whereŨ is a suitable functional space on the interval (0, 1). We will then consider its subspace U s 0 = U s ∩H 1 0 (Q), so that the semi-discretization of problem (6) will read: find u s ∈ U s 0 such that (7) a(u s , v s ) = b(v s ) for all v s ∈ U s 0 . We will prove (see Theorem 6) that U s 0 is a closed subspace ofH 1 0 (Q) so that problem (7) is uniquely solvable and the error between u and u s is bounded as usual by the best approximation (Céa's lemma): The solution of problem (7) is actually computed by solving a system of ordinary differential equations where the unknowns are the coefficients u i (r) of u s (r, θ). This procedure is detailed in Section 4.
In order to show the convergence of this procedure, we need to estimate the right-hand side of (8).

3.1.
Error estimates for the interpolation operator. We plan to construct an interpolation operator Π with values in U s that, if applied to smooth functions, would act as follows Since we will work with Sobolev functions, it is useful to define an adequate trace-like operator that we are going to call the "radial trace operator". To this end, the following bound is required.
Reordering and integrating over r, we have , such that (9) holds for smooth functions.
The following space on the interval (0, 1) will be used for the definition ofŨ Given an angle 0 ≤ ϑ ≤ Θ, inequality (9) shows that the natural norm for a space U where the radial trace operator can be defined, is It is apparent that not all functions in C ∞ (Q) have a bounded U -norm because in general C ∞ (Q) is not included inH 1 (Q). This is due to the fact that ∂ θ u L 2 1/r (Q) might not be bounded for some u ∈ C ∞ (Q). Hence, we definẽ and extend it to the closure of C ∞ (Q) ∩H 1 (Q) with respect to the U -norm. We denote by U ⊂H 1 (Q) this space and by γ ϑ the extension of the trace operator, so that we have a bounded radial trace operator that extends the restriction operatorγ ϑ defined on smooth enough functions. It is then natural to chooseŨ = H 1 r (0, 1), so that the definition of U s reads as follows Theorem 6. The space of semi-discrete functions U s is a closed subspace ofH 1 (Q).
Proof. Let (u n ) be a Cauchy sequence in U s . By completeness, there is a functioñ u such that u n →ũ inH 1 (Q). By Lemma 5 we have as n, m → ∞, so (u n (·, θ i )) is a Cauchy sequence in L 2 r (0, 1) and by completeness there is a limit u n (·, as n → ∞, so u =ũ and therefore u n → u inH 1 (Q).
Given u ∈ U we can then define the interpolant as where u i (r) is defined as γ θi (u) in H 1 r (0, 1). In order to see that the interpolant is well defined, we need to show that Πu belongs toH 1 (Q). To limit the technicalities, from now on in this section we are assuming that {e i } is the basis of continuous piecewise linear finite elements on (0, Θ). The general case can be handled with similar arguments. We apply Lemma 5 and obtain (11) For the third term in (10), we fix r ∈ (0, 1) and observe that for i = 1, 2, . . . , N − 1. Taking into account that ∂ θ Πu is well defined in (θ i , θ i+1 ), we apply the mean value theorem and holds for someθ ∈ (θ i , θ i+1 ). Since Πu is linear in this interval, the following equality holds for all θ ∈ (θ i , θ i+1 ) and r ∈ (0, 1). After multiplying by 1/r, integrating, and applying Hölder's inequality, we have By integrating over each interval and summing up the terms, we have (12) ∂ θ Πu L 2 1/r (Q) ≤ ∂ θ u L 2 1/r (Q) . Finally, putting (11) and (12) together, we get In the next theorem we prove the approximation properties of U s . As usual, we need to assume suitable regularity that will be characterized by the following space The space U ′ requires extra regularity only in the angular variable θ. We will see in Section 5 that singular solutions (with respect to the Cartesian coordinates) can be in U ′ and be approximated optimally by SBFEM.

Constructing semi-discrete solutions with SBFEM
In order to solve our problem, the scaled boundary finite element method rewrites the formulation (7) as a system of ordinary differential equations.
Let us consider An immediate consequence of the definition of the space of semi-discrete functions U s 0 is that the bilinear form a : U s 0 × U s 0 → R and the linear form b : U s 0 → R may be rewritten as follows: We now proceed formally with the derivation of the differential equation. To this aim, we will use the following integration by parts formula which is clearly valid for smooth enough u i and v i . In order to simplify our notation we introduce a name for the space of the coefficients in U s as follows holds for all i = 1, . . . , N . Here, we have used that the matrix A is positive definite. Using a matrix-vector notation, (19) can be rewritten as a first order ODE by introducing the variable q(r) = rAu ′ (r) and defining the vector Then, it is straightforward to see that (19) is equivalent to (20) rx ′ (r) = Ex(r) + G(r) for r ∈ (0, 1) SBFEM provides a methodology for the construction of solutions for (20), see [9]. Whenever such solutions exist and satisfy the integration by parts formula for all v ∈ U s 0 , then they correspond with the unique solution in U s 0 .

Numerical examples
One of the main features of the method presented in the previous sections is that it can achieve high order of convergence also in presence of singular solutions. We are going to show this behavior with a simple example.
We take Θ = 3π/2 and consider the function u e (r, θ) = r  This function satisfies ∆u e = 0 on Ω, so it is a solution to the homogeneous problem: find u ∈ H 1 0 (Ω) such that − ∆u = 0 in Ω u = u e on ∂Ω.
Our theory can be easily extended to accommodate non-homogeneous boundary conditions. Moreover, it is well-known that u e belongs to H 1 (Ω) but not to H 2 (Ω). Indeed, this follows from the fact that the following inequality holds for any 0 < R < 1 and 0 < ε < R u e which tends to infinity as ε goes to 0. On the other hand, it can be easily seen that u e belongs to U ′ and this makes it possible to use the result of Theorem 9 which implies, in particular, that first order elements achieve second order of convergence in L 2 even in presence of a corner singularity.
We report in Figure 2 the results of our numerical test that confirm our theoretical findings. We include also higher order approximations (up to order six) and a convergence plot in H 1 .