INTEGRAL EQUATIONS FOR BIHARMONIC DATA COMPLETION

. A boundary integral based method for the stable reconstruction of missing boundary data is presented for the biharmonic equation. The solution (displacement) is known throughout the boundary of an annular domain whilst the normal derivative and bending moment are speciﬁed only on the outer boundary curve. A recent iterative method is applied for the data completion solving mixed problems throughout the iterations. The solution to each mixed problem is represented as a biharmonic single-layer potential. Matching against the given boundary data, a system of boundary integrals is obtained to be solved for densities over the boundary. This system is discretised using the Nystr¨om method. A direct approach is also given representing the solution of the ill-posed problem as a biharmonic single-layer potential and applying the similar techniques as for the mixed problems. Tikhonov regularization is employed for the solution of the corresponding discretised system. Numerical results are presented for several annular domains showing the eﬃciency of both data completion approaches.


1.
Introduction. Let u be a solution to the biharmonic equation (1) ∆ 2 u = 0 in D and suppose additionally that u satisfies the following conditions on the boundary, Here, D is an annular domain in IR 2 lying between the two simple closed nonintersecting curves Γ 2 and Γ 1 , with Γ 1 contained in the bounded interior of Γ 2 , and Γ = Γ 1 ∪ Γ 2 . The operators on the boundary are given by N u = ∂u ∂n , M u = ν∆u + (1 − ν) u x1x1 n 2 1 + 2u x1x2 n 1 n 2 + u x2x2 n 2 2 , where n = (n 1 , n 2 ) is the outward unit normal to Γ , = 1, 2, and 0 < ν < 1. In plate theory, M is the normal bending moment and ν the Poisson ratio. We assume that data is given such that there exists a classical or weak solution to (1)- (2). In [10], it is shown that (1)-(2) is ill-posed and that there exists at most one solution. Further in [10], an iterative method is proposed and analysed for the biharmonic data completion. At each iteration step, mixed boundary value problems are solved for (1) updating functions on the boundary curves. It is mentioned at the end of [10] that a single-layer approach can be used to solve the mixed problems numerically.
We shall undertake the laborious task of deriving and implementing a numerical solution strategy based on integral equations for solving (1) with mixed conditions. Moreover, it will be shown how to directly solve the biharmonic data completion problem using integral equations, that is without iterations. This in total constitute the novelty of the presented work. Related inverse problems for the biharmonic equation are mentioned in the introduction to [10]; for history and applications of the biharmonic equation, see [21].
For the Cauchy problem of the Laplace equation, boundary integral equations in combination with iteration schemes have been effective, see for example [3,19]. A direct integral approach with no iterations is given in [8] using ideas of [6]. We follow those works and extend the techniques to (1)- (2). The solution to the mixed problems in the iterative scheme is represented as a biharmonic single-layer potential with densities to be determined. Matching the given data, a system of boundary integral equations is obtained for the densities. This system is discretised using the Nyström method. We also investigate representing the solution to (1)-(2) directly as a single-layer potential following ideas in [8].
Layer potentials is a classical field with abundance of work (for direct problems), and we can therefore only give some guidance to where results and further references can be found for the biharmonic equation. The biharmonic single-layer potential that we use is studied for example in [11,16,22], see further [5,12,26]. Layer potential based methods for ill-posed Cauchy problems for other elliptic equations and the heat and wave equation are presented in [2,4,9,15,20,24,25].
For the outline of the work, in Section 2, the iterative method of [10] is briefly surveyed. In Section 3, the potential representation is given of the solution to the mixed biharmonic problems needed in the iterative procedure. Matching the boundary data renders a system of boundary integral equations for the densities, which is discretised using the Nyström method; details are in Section 4. In Section 5, a direct approach with no iterations is outlined. Properties of the obtained system of integral equations are shown in Theorem 5.1. Numerical examples are presented in Section 6, for both data completion approaches, showing that accurate reconstructions can be obtained. Some conclusions are stated in Section 7.
• Next, v 0 is constructed by solving (1) with the boundary conditions changed to v 0 = 0 on Γ, • Given that u k−1 and v k−1 are known, the approximation u k is determined from (1) with Here, The iterations continues repeating the last two steps until a suitable stopping rule has been satisfied. In the case of noisy data, a discrepancy principle can be applied to cease the iterations. Note that at each step, data is updated with data of the same kind (for example the normal bending moment is updated with a bending moment from a previous step). This is appealing from a physical point of view.
We recall from [10], has a solution u. Let the relaxation parameter γ > 0 be sufficiently small, and let u k be the k-th approximation in the given algorithm. Then lim k→∞ u − u k H 2 (D) = 0 for any initial function ζ 0 ∈ H −1/2 (Γ 1 ).
Precise statement on the size of γ is given in [10]. From the proof given in [10], it can be seen that the procedure is of Landweber type. Moreover, as a special case when the relaxation parameter γ = 1, the alternating method [17] is obtained.

3.
A single-layer approach for mixed biharmonic problems. We consider the mixed problem: with the boundary operators N and M as in (3). This mixed biharmonic problem has a unique solution u ∈ H 2 (D) for f ∈ H 3/2 (Γ), h 1 ∈ H −1/2 (Γ 1 ) and g 2 ∈ H 1/2 (Γ 2 ) that depends continuously on the data (see for example [10]). The fundamental solution for the biharmonic equation in IR 2 is and satisfies ∆ 2 x G(x, y) = δ(x − y), in IR 2 , with δ the Dirac delta function. We consider the single-layer potential for the biharmonic equation, (7) v(x) = where ϕ ∈ H r (Γ) and ψ ∈ H r+1 (Γ) for f ∈ H r+3 (Γ), g 2 ∈ H r+2 (Γ 2 ) and h 1 ∈ H r+1 (Γ 2 ), and bounds for r given in [5,11,16]. This single-layer potential has the following continuity and jump properties (see [11,Chapt. 8 , For the kernels in these relations, we have the representations (see [16,Chapt. 10.4.4]) Matching the representation (7) against the data (5) and taking into account the above continuity and jump properties, a system of boundary integral equations is obtained to be solved for the densities ϕ and ψ. To have a unique solution to this system, (7) has to be slightly modified. According to [11,Chapt. 8.4] the following holds.
Theorem 3.1. The solution u ∈ H 2 (D) of the mixed problem (4)-(5) can be represented as with the elements w(x) = a 0 + a 1 x 1 + a 2 x 2 ((a 0 , a 1 , a 2 ) ∈ IR 3 ), ϕ k ∈ H r (Γ k ) and ψ k ∈ H r+1 (Γ k ), k = 1, 2, being the unique solution of the system of integral equations consisting of for a given triplet Sobolev spaces for the boundary data are given in [11,Chapt. 8.4]. Note that the constants A 0 , A 1 and A 2 can be chosen arbitrarily, but the solution of (9) depends on the choice of them. The need for those constants and the additional linear term w in the biharmonic single-layer potential in order to have existence and uniqueness of solution to the above system is explained in [11,Chap. 8.3]. As is shown in [11, Ex. 8.2] it is not in general possible to put the constants a 0 , a 1 , a 2 , A 0 , A 1 and A 2 simultaneously equal to zero.
For the sake of completeness, we give the similar result for the second mixed problem used in the iterative method, From [11,Chapt. 8.4] we have, of the mixed problem (11)-(12) is given by with the elements w(x) = a 0 + a 1 x 1 + a 2 x 2 ((a 0 , a 1 , a 2 ) ∈ IR 3 ), and ϕ k ∈ H r (Γ k ) and ψ k ∈ H r+1 (Γ k ), k = 1, 2, being the unique solution of the system of integral equations consisting of together with (10) for a given triplet We shall outline in the next section how to discretise (9) and (14) taking into account the singularities in the kernels, together with discretisation of (10). 4. Discretisation of the systems (9), (10) and (14). We assume that the boundary curves Γ , = 1, 2, are sufficiently smooth and given by a parametric representation The system (9) can then be written in parametric form, and is expressed as with s ∈ [0, 2π], together with the parametric form of (10), The system (15)-(16) is well-posed in the corresponding Hölder or Sobolev space setting (see [5,11,16]). We use the notation ϕ (s) := ϕ (x (s))|x (s)| and ψ (s) := ψ (x (s))|x (s)| for the densities, and for the kernels, These kernels are all continuous functions, but derivatives of H ,H , L and Q do contain a logarithmic singularity. It is advantageous to make these singularities explicit when applying quadratures for the numerical solution of the system (15). Note here that the kernelQ does not contain any logarithmic term and its derivatives are smooth. We therefore rewrite the following kernels, where and Moreover, we have the following representation The other kernels can be written analogously. The diagonal terms are in general zero apart from We can then, taking into account (18), apply the Nyström method to (15)-(16) based on the following trigonometric quadrature rules [18], with mesh points (20) s k = kh, k = 0, . . . , 2m − 1, h = π/m, and the weight functions As a result, we obtain the linear system consisting of , N w 2i = N w(x 2 (s i )) and R j = R j (0), to be solved for the values ϕ kj ≈ ϕ k (s j ), ψ kj ≈ ψ k (s j ), k = 1, 2, j = 0, . . . , 2m − 1.
Combining the representation (8) with the jump relation from Section 3 for the bending moment M u, we obtain Inverse Problems and Imaging Volume 13, No. 5 (2019), 1095-1111 Therefore, using (17)-(18) and the above quadratures together with the densities from (22)-(23), the sought numerical values of M u on Γ 2 can be calculated as A similar expression can be derived for the normal derivative N u on Γ 2 . The numerical solution of the system (14) can be realized analogously, and the representation (13) can be used to find the requested data on the boundary. The quadratures and discretisation strategies applied are all well-studied. Thus, in principle, a full error analysis can be carried out following, for example, [18,Chapt. 12]. However, this is better done in a separate work, and instead we present some numerical examples in Section 6 that highlights that the method has the expected properties, in particular exponential convergence.

5.
A "direct" integral equation approach for the Cauchy problem (1)- (2). As an alternative to the iterative regularizing method introduced above, it is possible to use a representation of the form (8) and match not against the data of a mixed problem but the original data (2). This renders an ill-posed system of integral equations.
We thus search for the solution of (1)-(2) as a (modified) single-layer potential, where w(x) = a 0 + a 1 x 1 + a 2 x 2 and (a 0 , a 1 , a 2 ) ∈ IR 3 , ϕ k , ψ k ∈ C(Γ k ), k = 1, 2, together solve the system of integral equations consisting of (26) for a given triplet (A 0 , A 1 , A 2 ) ∈ IR 3 . To analyse this system, we introduce an operator formulation. Let Define the operator Here, it is assumed without loss of generality that Γ ds(y) = 1, and we made the trivial extension of g and h to Γ by putting g = h = 0 on Γ 1 . Note that ϕ| Γ k = ϕ k and ψ| Γ k = ψ k , k = 1, 2. The operator A is considered as a mapping The system (26)-(27) can be written We show properties of the mapping A following [6,7].
Theorem 5.1. The operator A corresponding to the system (26)-(27) and defined in (32) is injective and has dense range.
Proof. We start with injectivity. Let therefore Aζ = 0. Using the components of this elementζ, denoted as in (30), define u by (25); the restriction of ϕ to Γ k generates ϕ k and similarly ψ k is obtained. Then, from the assumption that Aζ = 0 and since the system (26)-(27) can be written as (33), u has the data (2) being zero. Due to the uniqueness of a solution to (1)-(2) shown in [10], u is therefore identically zero inD. In particular, u and its normal derivative are zero on the boundary Γ. To show that u being identically zero implies thatζ is zero, we need an additional operator. Let the matrixK consist of the first five rows of K from (28), and define Then from [11,Thm. 8.6],Ã is an isomorphism for the biharmonic equation (1) supplied with Dirichlet and Neumann data on Γ. Hence, since this data is zero for u, the constants and densities in the elementζ are all identically zero. Thus,ζ = 0 and we have shown injectivity of A.
To prove that A has a dense range, we show that the adjoint operator, A * , is injective. Formally, A * is obtained by (35) [ with notation as in (31). The transpose of K from (28)-(29) needed in (35) is when x ∈ Γ 2 and (37) Assume then that A * f = 0, withf as in (31). Define the layer potential for y ∈ IR 2 \Γ, and w(y) = A 0 +A 1 y 1 +A 2 y 2 . The element v is a well-defined solution to the biharmonic equation in the domain D. Due to the assumption A * f = 0, it follows using the last two rows of K T in (36)-(37), that v has zero Dirichlet data and normal derivative on the boundary Γ of D. Hence, v is identically zero in D.
The element v is zero also in a region containing Γ 2 in its interior (for an explicit extension formula for the biharmonic equation with zero Dirichlet and Neumann condition, see [13,Theorem 1] and [23, Theorem 3.2]; see further [1] and [14]). We now apply the normal derivative N y in (38). This derivative has a jump across Γ 2 due to the integral involving h, according to the jump relations stated in [16,Chapt. 2.4.1]. This and since v is identically zero, imply h = 0.
Then v in (38) with the third integral removed is a modified biharmonic singlelayer potential over Γ 2 . This potential generates an isomorphism (generated as in the first part of the proof, see (34), but using the first five columns of K T from (36)) for the biharmonic equation with Dirichlet and Neumann data specified on Γ, see [11,Thm. 8.6]. Using that and since the data on Γ is zero, it follows that the remaining constants and densities in (38) are identically zero, that isf = 0. Hence, A * is injective implying that A has dense range.  (19) with mesh points (20) and weights (21), lead to the system of linear equations 22 (s i , s j )]ϕ 2j for i = 0, . . . , 2m − 1, to be solved for ϕ kj ≈ ϕ k (s j ), ψ kj ≈ ψ k (s j ), k = 1, 2, j = 0, . . . , 2m − 1. As above, f i = f (x (s i )) and R j = R j (0). The system (39)-(40) has a high-condition number since (1)- (2) is ill-posed, and therefore Tikhonov regularization is incorporated. Denoting the matrix corresponding to this system by A and the right-hand side by F , the regularization means solving where A * is the adjoint (transpose) of A, with α > 0 a regularization parameter to be chosen appropriately. Although there are optimal choices for the regularization parameter (the discrepancy principle), we have chosen α for our numerical experiments by trial and error. We consider the source point z 1 = (3, 0) and use as boundary data the restriction of the corresponding fundamental solution G(x, z) and the values N x G(x, z) and M x G(x, z). In the systems (9) and (14) supplemented with (10) corresponding to the two mixed problems, we take A 0 = A 1 = A 2 = 1. Table 1 contains the discrete L 2 relative errors for the normal bending moment M u on Γ 2 for the mixed problem (4)-(5) and correspondingly M v on Γ 1 for (11)- (12), that is The element u ex and v ex are the respectively exact solution andũ andṽ the corresponding numerical ones. Here, m is the parameter for the number of collocation points on each boundary curve, see further (20). The sub-index 2M indicates that it is the L 2 error of the bending moment M . To calculate M u on the boundary Γ 2 the expression (24) is used, and similarly M v on Γ 1 is generated. The expected exponential convergence order is clearly exhibited in Table 1. We have tested our approach for the mixed problems on other domains and data with the similar results. The quadratures applied are well-established with error estimates both in spaces of smooth functions and in Sobolev spaces. Thus, no real surprise is to be expected. We therefore move on to the biharmonic data completion problem. Ex. 2 (data completion, ill-posed): We consider again an arbitrary source point z 1 ∈ IR 2 \ D, and construct the needed boundary data functions The exact solution corresponding to (42) is u ex given by (6) and the computed one is denoted byũ. Thus, we can compare the exact and the numerically calculated data on the inner boundary Γ 1 . We choose the boundary curves as in the first example. The initial guess ζ 0 that starts the procedure is taken to be zero.
The discrete relative L 2 errors of the reconstructions for the data completion problem (1)-(2) on Γ 1 with the proposed iterative approach, that is Fig. 1. The sub-index 2N means the (discrete) L 2 norm of the normal derivative (N ). In generating the collocation (grid) points m = 32 was chosen (see further (20)), and in the iterations the relaxation parameter is γ = 0.5, and the source point is z 1 = (3, 0).
For noisy data, random pointwise errors are added to the corresponding boundary function, with the percentage given in terms of the L 2 -norm. The minimal values obtained for the errors are e 2N = 0.00328 and e 2M = 0.04886 for exact data after 300 iterations, and e 2N = 0.06675 and e 2M = 0.44542 for 3% noise in the data after 18 and 11 iterations, respectively. It is natural that the reconstruction of the normal bending moment M u is less accurate than the normal derivative N u,  Table 2. The errors in the second example calculated on Γ 1 for the "direct" single-layer approach with exact (δ = 0%) and noisy data (δ = 3%).
since M u contains higher derivatives. It is an advantage with the integral approach that also boundary data with derivatives can be completed accurately; it is due to the exact expressions for M u and N u that can be derived and discretised on the boundary for the mixed problems.   Table 2 contains the corresponding errors obtained for the data completion (1)-(2) using the "direct" single-layer approach from Section 5, with α being the Tikhonov regularization parameter in (41). The choice of the constants A 0 , A 1 and A 2 in (27) is as above. These constants can be chosen freely in both methods.
Comparing with the results of the iterative scheme, for both for exact and noisy data the layer approach is more accurate. However, for noisy data, the discrepancy between the two methods are smaller. This is mainly due to the high condition number of the discretised system (39) making it hard to regularize, whilst in the iterations well-posed mixed problems and well-conditioned systems are solved. We can report the similar conclusion between the two methods also for the other values of m.
The reconstructions of N u and M u corresponding to Fig. 1 and Table 2 follow the respectively exact solution well as indicated by the reported error levels. Plots of the reconstructions versus exact solutions are of the form similar to what have been reported in for example [3,7,8], and figures of the approximations are therefore left out. 7. Conclusion. A boundary data completion problem has been investigated for the biharmonic equation in planar annular domains. The missing normal derivative and bending moment on the inner boundary are reconstructed from known deflection on the boundary and the normal derivative and bending moment on the outer boundary curve. An iterative method was employed solving mixed problems at each iteration step. An efficient solver for these mixed problems based on a single-layer representation was developed. Numerical experiments included showed that accurate and stable numerical solutions can be obtained. As an alternative, a direct single-layer approach was developed for data completion without iterations. With exact data the direct layer approach produces reconstructions with higher accuracy than the iterative procedure, whilst for noisy data the situation is reversed.