Boundary element methods for Dirichlet boundary control problems

In this article we discuss the application of boundary element methods for the solution of Dirichlet boundary control problems subject to the Poisson equation with box constraints on the control. The solutions of both the primal and adjoint boundary value problems are given by representation formulae, where the state enters the adjoint problem as volume density. To avoid the related volume potential we apply integration by parts to the representation formula of the adjoint problem. This results in a system of boundary integral equations which is related to the Bi‐Laplacian. For the related Dirichlet to Neumann map, we analyse two different boundary integral representations. The first one is based on the use of single and double layer potentials only, but requires some additional assumptions to ensure stability of the discrete scheme. As a second approach, we consider the symmetric formulation which is based on the use of the Calderon projector and which is stable for standard boundary element discretizations. For both methods, we prove stability and related error estimates which are confirmed by numerical examples. Copyright © 2010 John Wiley & Sons, Ltd.


Introduction
Optimal control problems of elliptic or parabolic partial differential equations with a Dirichlet boundary control play an important role [1]. In particular in fluid mechanics, the cost functional can be the integral over the strain tensor of the velocity field, where the Dirichlet control describes the inflow velocity [2]. A similar minimization problem is considered in [3], where the cost functional is the boundary integral to describe the work needed to overcome the drag exerted on a given body. A difficulty in handling the Dirichlet control problems is the choice of the control space, where the Sobolev trace space, H 1/2 ( ), appears as a natural choice [4]. To obtain smoother optimal solutions one may even consider H 2 ( ) as control space, where the Sobolev norm for sufficient regular boundaries can be realized by using the Laplace-Beltrami operator [5]. The most popular choice is to consider L 2 ( ) as control space. Although this choice allows the use of a piecewise constant control function, the associated partial differential equation has to be considered within an ultra-weak variational formulation, see, for example, [6] and [7] for an appropriate finite element approximation using standard piecewise linear basis functions. The use of the ultra-weak variational formulation of the primal Dirichlet boundary value problem in the context of an optimal control problem requires the adjoint variable p to be sufficiently regular, i.e. p ∈ H 2 ( )∩H 1 0 ( ). As the adjoint variable p itself is the unique solution of the adjoint partial differential equation with homogeneous Dirichlet boundary conditions, either a smooth boundary , or a polygonal or polyhedral but convex domain has to be assumed. For related finite element approximations, see, e.g., [2, 8--11] or [12] in the case of a finite dimensional Dirichlet control. To include a Dirichlet boundary condition, u = z ∈ L 2 ( ), in a standard variational formulation, alternatively one may consider a penalty approximation of the Dirichlet boundary condition by using a Robin boundary condition [4, 13--15]. Again, sufficient smoothness of the boundary has to be assumed.
In [16], a finite element approach was considered, where the energy norm was realized by using some hypersingular boundary integral operator which links the Dirichlet control with the normal derivative of the adjoint variable. The related optimality condition results in a higher regularity of the control and requires weaker assumptions on the smoothness of the adjoint variable, in fact, one may even consider general Lipschitz domains . Moreover, for polygonal-or polyhedral-bounded domains one also obtains higher order convergence results for the approximate finite element solution.

Dirichlet control problems
In boundary control problems, a Dirichlet boundary data z is to be determined such that the state u as the solution of the related Dirichlet boundary value problem matches a desired target u taking into account the costs of the control. As a model problem, we consider, for a bounded domain ⊂ R n , n = 2, 3, with Lipschitz boundary = * , the Dirichlet boundary control problem to minimize the cost functional subject to the constraint and where the control z satisfies the box constraints We assume f , u ∈ L 2 ( ), ∈ R + , and g a , g b ∈ H 1/2 ( ) satisfying g a g b on . Moreover, we use the hypersingular boundary integral operator D : H 1/2 ( ) → H −1/2 ( ) to describe the cost, or some regularization term, via a semi-norm in H 1/2 ( ). In particular, for z ∈ H 1/2 ( ) we have is the fundamental solution of the Laplace operator. Note that the hypersingular integral operator D is spectrally equivalent to the Steklov-Poincaré operator, which describes the energy of the harmonic extension of the Dirichlet control. Let u f ∈ H 1 0 ( ) be the weak solution of the homogeneous Dirichlet boundary value problem The solution of the Dirichlet boundary value problem (2) is then given by u = u z +u f , where u z ∈ H 1 ( ) is the unique solution of the Dirichlet boundary value problem Note that the solution of the Dirichlet boundary value problem (5) defines a linear map u z = Sz with S : H 1/2 ( ) → H 1 ( ) ⊂ L 2 ( ). Then, by using u = Sz +u f , we consider the problem to find the minimizer, z ∈ U ⊂ H 1/2 ( ), of the reduced cost functional To characterize the minimizer, z ∈ U, of the reduced cost functional (6) we introduce the self-adjoint, bounded and H 1/2 ( )-elliptic operator satisfying [16] T z, z c T Moreover, we define Hence we can rewrite the reduced cost functional (6) as As U ⊂ H 1/2 ( ) is convex and closed, and since T is self-adjoint and H 1/2 ( )-elliptic, the minimization of the reduced cost functional, J(z), is equivalent to solving the variational inequality to find z ∈ U such that As (9) is an elliptic variational inequality of the first kind, we can use standard arguments as given, for example in [6, 25--27], to establish unique solvability of the variational inequality (9). By using the primal variable u = Sz +u f , and by introducing the adjoint variable = S * (u−u) ∈ H −1/2 ( ), we can rewrite the variational inequality (9) as Note that for given z ∈ H 1/2 ( ) and f ∈ L 2 ( ) the application of u = Sz +u f corresponds to the solution of the Dirichlet boundary value problem (2). The application of the adjoint operator, = S * (u−u), is characterized by the Neumann datum where p ∈ H 1 0 ( ) is the unique solution of the adjoint Dirichlet boundary value problem As the unknown control z ∈ U ⊂ H 1/2 ( ) is considered on the boundary = * , the use of boundary integral equations to solve both the primal boundary value problem (2) and the adjoint boundary value problem (12) seems to be a natural choice. In what follows, we will describe and analyse two different boundary element methods to solve the variational inequality (9) numerically. This will be based on the use of appropriate boundary integral operator representations of T and g as introduced in (7) and (8), respectively.

Primal boundary value problem
The solution of the Dirichlet boundary value problem (2), is given by the representation formula for x ∈ , u( x) = U * ( x, y)t(y)ds y − * *n y U * ( x, y)z(y)ds y + U * ( x, y)f (y)dy, (Vt)(x) = U * (x, y)t(y)ds y for x ∈ is the Laplace single layer integral operator V : H −1/2 ( ) → H 1/2 ( ), and (Kz)(x) = * *n y U * (x, y)z(y)ds y for x ∈ is the Laplace double layer integral operator K : H 1/2 ( ) → H 1/2 ( ). Moreover, is the related Newton potential. The single layer integral operator V is H −1/2 ( )-elliptic; for n = 2 we assume the scaling condition diam < 1 to ensure this. For the solution of the boundary integral Equation (14) we therefore obtain

Adjoint boundary value problem
The solution of the adjoint Dirichlet boundary value problem (12), is given correspondingly by the representation formula for x ∈ , where q = (* / *n)p ∈ H −1/2 ( ) is the unique solution of the boundary integral equation Remark 3.1 While the boundary integral equation (14) can be used to determine the unknown Neumann datum, t ∈ H −1/2 ( ), of the primal Dirichlet boundary value problem (2), the unknown Neumann datum, q ∈ H −1/2 ( ), of the adjoint Dirichlet boundary value problem (12) is given as the solution of the boundary integral equation (17). Then, by using =−q, the control z ∈ H 1/2 ( ) is determined by the variational inequality (10). However, since the solution u of the primal Dirichlet boundary value problem (2) enters the volume potential N 0 u in the boundary integral equation (17), it seems to be necessary to include the representation formula (13).
In this case, we would have to solve a coupled system of boundary and domain integral equations, which still would require some domain mesh. Instead, we will now describe a system of only boundary integral equations to solve the adjoint boundary value problem (12).
To end up with a system of boundary integral equations only, instead of (16), we will introduce a modified representation formula for the adjoint state p as follows. First we note that is a solution of the Poisson equation i.e., V * (x, y) corresponds to the fundamental solution of the Bi-Laplacian. Hence we can rewrite the volume integral for u in (16), by using Green's second formula, as follows: Therefore, we now obtain from (16) the modified representation formula for x ∈ , where the volume potentials involve given data only.
The representation formula (20) results, when taking the limit x → x ∈ , in the boundary integral equation for x ∈ , which can be written as Note that In addition, we have introduced a second Newton potential, which is related to the fundamental solution of the Bi-Laplace operator, With (15), we conclude from (21) the boundary integral equation By replacing =−q in (10) we, therefore, obtain a boundary integral representation of the operator T as defined in (7), and a related representation for g as defined in (8), To investigate the unique solvability of the variational inequality (9) based on the boundary integral representations (23) and (24), we will first recall some properties of boundary integral operators which are related to the Bi-Laplace partial differential equation, see also [19,21,22].

Bi-Laplace boundary integral equations and properties of T .
In this section, we consider a representation formula and related boundary integral equations for the Bi-Laplace equation which can be written as a system, As for the Laplace equation we can first write the boundary integral equations and where is the adjoint Laplace double layer integral operator K : are the associated Cauchy data on .
To obtain a representation formula for the solution u of the Bi-Laplace Equation (25), we first consider the related Green's first formula and in the sequel Green's second formula, When choosing v(y) = V * ( x, y) for x ∈ , i.e. the fundamental solution (18) of the Bi-Laplace operator, the solution of the Bi-Laplace partial differential equation (25) is given by the representation formula for x ∈ by By using (19), this can be written as Hence we obtain the boundary integral equation Moreover, when taking the normal derivative of the representation formula (30), this gives another boundary integral equation for x ∈ , where is the adjoint Bi-Laplace double layer integral operator K 1 : The boundary integral Equations (27), (28), (31), and (32) can now be written as a system, including the so-called Calderon

Lemma 4.1
The Calderon projection C as defined in (33) is a projection, i.e. C 2 = C.

Proof
The proof follows as in the case of the Laplace equation [20,24], for the Bi-Laplace equation, see also [22].
From the projection property as stated in Lemma 4.1 we obtain some well-known relations of all boundary integral operators, which were introduced for both the Laplace and the Bi-Laplace equations. and

Proof
The relations of (34) for the Laplace operator are well known, see, e.g., [24], for the Bi-Laplace operator, see also [22].
To prove the ellipticity of the boundary integral operator T as defined in (23), we use the following result.

Lemma 4.3
For any t ∈ H −1/2 ( ) there holds the equality where Proof For x ∈ and t ∈ H −1/2 ( ), we define the Bi-Laplace single layer potential which is a solution of the Bi-Laplace differential equation (25). Then, the related Cauchy data are given by On the other hand, for x ∈ , is a solution of the Laplace equation. Hence, the related Cauchy data are given by Now, for u = v = u t , Green's first formula (29) reads as and therefore we conclude The assertion now follows with w t = Vt. Now we are able to state the mapping properties of the boundary integral operator T as defined in (23), see also the properties of T as introduced in (7).

A non-symmetric boundary element method
be a boundary element space of piecewise linear and continuous basis functions i , which are defined with respect to a globally quasi-uniform and shape regular boundary mesh H of mesh size H. Define the discrete convex set Then the Galerkin discretization of the variational inequality (9) is to find z H ∈ U H such that

Proof
The error estimate (42) in the energy norm follows from the general abstract theory as presented, e.g. in [28,29], see also [26]. The error estimate (43) follows from the Aubin-Nitsche trick for variational inequalities, see [30] for the case U H ⊂ U and [31] for the more general case U H ⊂ U.
Although the error estimates (42) and (43) seem to be optimal, the operator T as considered in the variational inequality (41) does not allow a practical implementation, since this would require the discretization of the operator T as defined in (23), which is not possible in general. Hence, instead of (41) we need to consider a perturbed variational inequality to find z H ∈ U H such that where T and g are appropriate approximations of T and g, respectively. The following theorem, see, e.g., [16], presents an abstract consistency result, which will later be used to analyse the boundary element approximation of both the primal and adjoint boundary value problems.
where z H ∈ U H is the unique solution of the discrete variational inequality (41).
In the remaining part of this paper we will analyse approximations T and g which are based on the use of boundary element methods. This is related to the boundary element approximation of Steklov-Poincaré operators [32--35].

Boundary element approximation of T
For an arbitrary but fixed z ∈ H 1/2 ( ), the application of T z reads as where q z ∈ H −1/2 ( ) is the unique solution of the boundary integral equation and t z ∈ H −1/2 ( ) solves For a Galerkin approximation of the above boundary integral equations, let be another boundary element space of, e.g., piecewise constant basis functions k , which are defined with respect to a second globally quasi-uniform and shape regular boundary element mesh of mesh size h. Now, q z,h ∈ S 0 h ( ) is the unique solution of the Galerkin formulation  (46) is bounded, i.e.

Proof
The assertion is a direct consequence of the mapping properties of all boundary integral operators involved, we skip the details.

Proof
For an arbitrary chosen but fixed z ∈ H 1/2 ( ) we have, by definition, By using (46), we also have and therefore Let us further define q z,h ∈ S 0 h ( ) as the unique solution of the variational problem Then the perturbed Galerkin orthogonality follows. From this we further conclude The assertion now follows from the triangle inequality, and by applying Cea's lemma.
By using the approximation property of the trial space S 0 h ( ) and the Aubin-Nitsche trick, we conclude an error estimate from (47) when assuming some regularity of q z and t z , respectively.

Corollary 5.5
Assume q z and t z to be piecewise smooth, i.e. q z ,t z ∈ H s pw ( ) for some s ∈ [0, 1]. Then there holds the error estimate

Boundary element approximation of g
As in (46), we may also define a boundary element approximation of the right-hand side g as defined in (24), In particular, g ∈ H −1/2 ( ) is the unique solution of the variational problem where t f = V −1 N 0 f ∈ H −1/2 ( ) solves the variational problem Hence we can define a boundary element approximation, g h ∈ S 0 h ( ), as the unique solution of the Galerkin variational problem where t f,h ∈ S 0 h ( ) is the unique solution of the Galerkin problem

Lemma 5.6
Let g be the right-hand side as defined in (24), and let g h be the boundary element approximation as defined in (50). Then there holds the error estimate

Proof
The assertion follows as in the proof of Lemma 5.4, we skip the details.
By using the approximation property of the trial space S 0 h ( ) and the Aubin-Nitsche trick, we conclude an error estimate from (52) when assuming some regularity of g and t f , respectively.

Corollary 5.7
Assume g, t f ∈ H s pw ( ) for some s ∈ [0, 1]. Then there holds the error estimate

Approximate variational inequality
We consider the variational inequality (10) with =−q to find z ∈ U such that where q ∈ H −1/2 ( ) is the unique solution of the boundary integral equation and t ∈ H −1/2 ( ) is the unique solution of The Galerkin boundary element approximation of the variational inequality (54), and therefore the boundary element discretization of the perturbed variational inequality (44), is to find z H ∈ U H such that where q h ∈ S 0 h ( ) is the unique solution of the Galerkin formulation and t h ∈ S 0 h ( ) solves The Galerkin formulation (58) is equivalent to the linear system and (59) is equivalent to where Moreover, we are also able to derive an error estimate in L 2 ( ), i.e.
when applying the Aubin-Nitsche trick.
In the particular case of a non-constrained minimization problem, instead of the discrete variational inequality (62) we have to solve the linear system which is equivalent to the system ⎛

Remark 5.10
The error estimates (64) and (65) provide optimal convergence rates when approximating the control z by using piecewise linear basis functions. However, we have to assume h c 0 H to ensure the unique solvability of the perturbed Galerkin variational inequality (44), where the constant c 0 is in general unknown. Moreover, the matrix T ,H as given in (63) defines a non-symmetric approximation of the self-adjoint operator T . Hence we are interested in deriving a symmetric boundary element method which is stable without any additional constraints in the choice of the boundary element trial spaces.

A symmetric boundary element method
The boundary integral formulation of the primal boundary value problem (2) is given by (14), while the adjoint boundary value problem (12) corresponds to the modified boundary integral equation (21). In what follows, we will use a second boundary integral equation of the adjoint boundary value problem to obtain an alternative representation for q and therefore of the adjoint operator S * . In particular, when computing the normal derivative of the representation formula, (20), this gives where Hence, from (11) we obtain = S * (u−u) =−q =−( 1 2 I+K )q+D 1 z +K 1 t +N 1 u+M 1 f, and, by using (15) and (22), we conclude the alternative representations and  In the particular case of a non-constrained minimization problem, instead of the discrete variational inequality (79) we have to solve the linear system which is equivalent to a system of linear equations, ⎛ Remark 6.7 The symmetric boundary element approximation T ,H is positive definite for any choice of conformal boundary element spaces S 1 H ( ) ⊂ H 1/2 ( ) and S 0 h ( ) ⊂ H −1/2 ( ). In particular we may use the same boundary element mesh with mesh size h = H to define the basis functions i and k , respectively. From a theoretical point of view, this is not possible in general when using the non-symmetric approximation T ,H .

Numerical results
As numerical example we first consider as in [8,11], see also [16], the unconstrained Dirichlet boundary control problem (1)-(2) for the domain = (0, 1 2 In this case we have to solve the coupled linear system (66) in the case of the non-symmetric boundary element approach, and (82) for the symmetric approach. For the boundary element discretization, we introduce a uniform triangulation of the boundary = * on several levels where the mesh size is h L = 2 −(L+1) . As the minimizer of (1) is not known in this case, we use the boundary element solution z h 9 of the ninth refinement level as reference solution. The boundary element discretization is done by using the trial space S 0 h ( ) of piecewise constant basis functions, and S 1 h ( ) of piecewise linear and continuous functions. In particular, we use the same boundary element mesh to approximate the control z by a piecewise linear approximation, and piecewise constant approximations for the fluxes t and q. Note that we have h = H in this case, and therefore we can not ensure the S 1 h ( )-ellipticity of the non-symmetric boundary element approximation, see Theorem 5.8. However, the numerical example shows stability in this case.
In Table I, we present the errors for the control z in the L 2 ( ) norm and the estimated order of convergence (eoc). These results correspond to the error estimate (65) of the non-symmetric boundary element approximation, and to the error estimate (81) of the symmetric boundary element approximation. Note that in this example we have u ∈ H s ( ) for s < 2 / 3 implying p ∈ H 2+s ( ), and (* / *n)p ∈ H 1/2+s ( ). As D : H 1/2+ ( ) → H −1/2+ ( ) for 2 / 3, see [36], we further conclude z ∈ H 7/6 ( ), and therefore we can expect 7 / 6 as order of convergence. For comparison, we also give the error of the related finite element solution, see [16]. From the numerical results we conclude, that all three different approaches behave almost similar, as predicted by the theory.
As a second example, we consider in addition the box constraints (3) with g a =−1, which is inactive, and g b = 2.23. In Figure 1 we give a comparison of the unconstrained and constrained solutions, and in Figure 2 we plot the related controls for x 1 ∈ (0, 0.5),x 2 = 0.

Concluding remarks
In this paper, we have shown that we can use boundary element methods to solve Dirichlet boundary control problems. The first approach is based on the first Bi-Laplace boundary integral equation only, but a stable boundary element discretization requires   the use of appropriate boundary element spaces. However, the symmetric formulation is stable for all standard boundary element spaces without any further condition. Moreover, the Galerkin discretization of the symmetric formulation results in a symmetric approximation, while the first approach does not. Hence, the symmetric formulation seems to be the method of choice. The numerical results coincide with those of a comparable finite element approach. The advantage of using boundary element methods lies in the fact that only a discretization of the boundary is required. In the case of smooth data we can prove, with respect to the used lowest order trial spaces, the best possible order of convergence for the boundary element approximation of the control z, whereas for a finite element approximation we are only able to prove some reduced order, see [16]. Moreover, optimal control problems subject to partial differential equations in unbounded exterior domains can be handled analogously. While this paper is on the stability and error analysis of boundary element methods for optimal control problems only, further research will be done for an efficient solution of the resulting discrete systems. Hereby, special focus will be on appropriate solution methods to solve the discrete variational inequalities. This also involves the construction of efficient preconditioners, as well as the use of fast boundary element methods.