PDE-CONSTRAINED OPTIMAL CONTROL APPROACH FOR THE APPROXIMATION OF AN INVERSE CAUCHY PROBLEM

. This paper concerns the approximation of a Cauchy problem for the elliptic equation. The inverse problem is transformed into a PDE-constrain-ed optimal control problem and these two problems are equivalent under some assumptions. Diﬀerent from the existing literature which is also based on the optimal control theory, we consider the state equation in the sense of very weak solution deﬁned by the transposition technique. In this way, it does not need to impose any regularity requirement on the given data. Moreover, this method can yield theoretical analysis simply and numerical computation conveniently. To deal with the ill-posedness of the control problem, Tikhonov regularization term is introduced. The regularized problem is well-posed and its solution converges to the non-regularized counterpart as the regularization parameter approaches zero. We establish the ﬁnite element approximation to the regularized control problem and the convergence of the discrete problem is also investigated. Then we discuss the ﬁrst order optimality condition of the control problem further and obtain an eﬃcient numerical scheme for the Cauchy problem via the adjoint state equation. The paper is ended with numerical experiments. 2010 Mathematics Subject Classiﬁcation. Primary: 49J20, 65N21; Secondary: 65N30.

1. Introduction. The inverse Cauchy problem for elliptic partial differential equations, which consists in recovering data on the whole boundary of the domain from partial but over-determined measurements, is a typical ill-posed problem of mathematical physics. Specially, we consider a regular bounded domain with Lipschitz boundary and assume that the boundary of this domain is partitioned into two parts. One part is over-specified since both the Dirichlet and Neumann conditions of an elliptic PDE are given on it. The other part is under-specified because of no conditions are prescribed on this part of the boundary. One needs to solve the elliptic PDE and then reconstruct the corresponding conditions on the second part of the boundary. The Cauchy problem, a special boundary value problem for elliptic PDE, arises in electrostatic and has a large number of applications in many fields such as thermal inspection, electrical prospecting, geophysics and medical imaging and so on. Taking the detection of pipeline corrosion by electrical field as an example, this process can be modeled by a Cauchy problem for a Laplace equation in an annulus domain under the assumption that the thickness of the pipeline is relatively small compared with its radius. The purpose is to determine the information about the corrosion that occurs on the interior surface of the pipeline, but the only accessible data are the electrostatic measurements on the exterior surface. In the field of image processing, image inpainting techniques based on crack detection can reconstruct a damaged part of an image through solving an inverse Cauchy problem. More practical examples may be listed in a broader catalogue (see, e.g., [2,3,4,10,16,20,23,27,28,32,38]). Before introducing a model Cauchy problem for a linear second order elliptic PDE, we firstly give some notations. Let Ω ⊂ R d (d = 2, 3) be a bounded domain with a Lipschitz boundary Γ := ∂Ω and consider the elliptic differential operator where a 0 ∈ L ∞ (Ω), a 0 (x) ≥ 0 a.e. x ∈ Ω; a ij is Lipschitz on Ω (1 ≤ i, j ≤ d); a ij (x)ξ i ξ j ≥ c|ξ| 2 , c > 0, ∀ x ∈ Ω, ∀ ξ ∈ R d ; a ij n i ∂ xj , n ∈ R d is the unit outward normal to Γ; Γ = Γ 1 ∪ Γ 2 , Γ 1 , Γ 2 are nonempty closed subsets of Γ,Γ 1 ∩Γ 2 = ∅.
The model Cauchy problem we consider in the present paper takes the form Here f, g D and g N are given functions. The Cauchy problem for elliptic PDE is well known to be highly ill-posed. In many cases, this kind of problem does not have a solution. Even though the solution exists it is not stable with respect to the given data. That means the solution does not depend continuously on the given data and any small perturbation in the given data may cause large change to the solution. Theories on such kind of problem have been extensively studied for decades and are the subject of numerious references (see, e.g., [6,21,22,28,40]). On the other hand, investigating numerical issues of this problem also be a common concern due to the fact that many engineering processes can be modeled by the Cauchy problem. There are two large groups of numerical approaches to the solution of this kind of problem. The first kind comprises the methods based on regularization which transforms the ill-posed problem into a well-posed one using some regularization techniques. The solution of the regularized problem can approximate the corresponding ill-posed problem well by tuning the regularisation parameters. The lack of stability has promoted an important activity to develop relevant regularization schemes. In [41], the authors developed a scheme which is now widely used and known as Tikhonov regularization method (see [5,15,39] and the references therein). Another regularization method known as quasi-reversibility attributes to Latter and Lions [34]. Pertinent literatures dedicated to this topic are available in [8,31]. The other group consists of iterative methods which have received much attention in recent years. Among them, we should mention a convergent iterative method based on an alternating procedure proposed by Kozlov et al. [33]. One may refer to [13,18,29,42] for more details about the iterative methods.
In this paper, we study a numerical method for solving the Cauchy problem. The numerical method is based on a reformation of the Cauchy problem via PDEconstrained optimal control theory [26] coupling with Tikhonov regularization method. The objective functional, the state equation and the control variable in a boundary control problem provide an appropriate way to handle the over-specified and under-specified parts of the boundary in the Cauchy problem. The first result based on optimal control approach for problem (1) is given in [30]. The authors used the unknown solution value on Γ 2 as a control variable u and then the following Dirichlet boundary control problem is derived in Ω, where U ad is the admissible set for the control variable. Note that the compatibility of state equation (3) requires that u − g D ∈ H 1 2 00 (Γ 2 ) (the set of functions in H 1 2 (Γ) that vanish on Γ 1 ). This restriction causes difficulty for the purpose of developing efficient numerical methods. After that Chakib et al. [12] adopted the following Neumann boundary control approach to solve the Cauchy problem. In this reformation, the unknown flux ∂ n A y on Γ 2 is used as the control variable min u∈U ad They also proved that problem (4)-(6) is equivalent to the original problem (1) under some regularity assumptions on the data, and thus the ill-posedness of the original problem is inherited. Subsequently, Han et al. [24] improved the result of [12] by introducing a regularization term to deal with the ill-posedness, which minimizes instead of (4). It is evident that two state equations are involved in this scheme, which leads to more complexities in theoretical analysis and numerical computation. Moreover, the above two mentioned approaches all require that g D ∈ H 1 2 (Γ 1 ) to ensure the well-definition of the optimal control problems.
A more natural formulation which will be exploited in this paper is to use the unknown solution value on Γ 2 as the control variable and (4) as the objective functional. In this case the state equation should not be understood in the classical weak sense, which requires H 1 2 (Γ 2 )-Dirichlet boundary data to ensure the existence of a solution. Instead, we consider the state equation in the context of very weak solution defined by the transposition technique (see e.g., [7,37]). Compared with the above two schemes, our scheme has two obvious advantages. Firstly, it does not need to impose more regularity assumptions on the data f, g D and g N . For g D , g N ∈ L 2 (Γ 1 ), which are commonly encountered in applications, the new scheme is well-defined. Secondly, only one state equation is involved in the new scheme which makes the analysis and computation conveniently.
The rest of the paper is organized as follows. In section 2, we introduce the very weak solution to the elliptic Dirichlet boundary value problem and formulate our PDE-constrained optimal control problem. The equivalence between the control problem and the Cauchy problem is proved. To obtain a well-posed control problem, a Tikhonov regularization term is introduced. The convergence of the regularized control problem to the non-regularized one is also proved as the regularization parameter approaches zero. In section 3, we establish the finite element approximation to the state equation and the regularized control problem. The a priori error analysis for the finite element approximation posed on convex, polygonal or polyhedral domains is studied. In section 4, an efficient numerical scheme for the Cauchy problem is obtained through an adjoint state equation. Finally, some numerical experiments are carried out in section 5 to support the theoretical results.
Some notations. Let Ω ⊂ R d , 2 ≤ d ≤ 3 be either convex and polyhedral or ∂Ω belong to C 2 . For a bounded domain Ω we adopt the standard notation W m,s (Ω) for Sobolev spaces (see [1]) on Ω with norm · m,s,Ω and seminorm | · | m,s,Ω , and the standard modification for s = ∞. We denote by H m (Ω) with norm · m,Ω and seminorm | · | m,Ω for s = 2. The space H . Denote the L 2inner product in Ω and on Γ by (·, ·) and ·, · Γ , respectively. Furthermore, the space Here ∂ n A * is similar to ∂ n A and A * is the adjoint operator of A such that In addition, c and C denote generic positive constants.
2. Problem reformation. The formulation based on optimal control approach we utilize in this paper is to use the unknown solution value on Γ 2 as the control variable. The objective functional is the same as (4). Thus, we can obtain a Dirichlet boundary control problem which takes the form in Ω, Here the state equation (9) should not be understood in the classical weak sense, where u ∈ H 1 2 (Γ 2 ) is required to ensure the existence of a solution. Instead, we consider its very weak solution based on the transposition technique (see e.g., [7,37]). The very weak solution can be defined as follows.
Lemma 2.1. Given f ∈ L 2 (Ω), g N ∈ L 2 (Γ 1 ) and u ∈ L 2 (Γ 2 ), the boundary value problem (9) admits a unique solution y ∈ L 2 (Ω) which is denoted by y = y(u) in the sense that Moreover, y ∈ H 1 2 (Ω) and the following Lipschitz property holds: where L 1 is the Lipschitz constant. Finally, if u n u in L 2 (Γ 2 ), then y n := y(u n ) y in H Proof. For the proof we refer to [11] and [37]. Here a proof following the idea of [17] is included. Let us briefly sketch the existence of y in the case of f ≡ 0 and g N ≡ 0. For each ψ ∈ L 2 (Ω), let w ∈ H 2 Γ1 (Ω) ∩ H 1 Γ2 (Ω) be the solution of the following equation and T : L 2 (Ω) → L 2 (Γ 2 ) the linear operator which is defined by T ψ := −∂ n A * w| Γ2 . Denote T * the adjoint of T and y = T * u, then which verifies that y satisfies (10). Note that L 2 (Ω) embeds continuously into H − 1 2 (Ω), from the trace theorem and regularity of elliptic equation, we have (Ω) follows immediately from the inequality above for each ψ ∈ L 2 (Ω). Moreover, we have The uniqueness of the solution can be obtained from (12) and the proof of (11) is similar to (12).
Note that the solution operator y = y(u) from L 2 (Γ 2 ) to H 1 2 (Ω) is affine linear and continuous, and thus is weakly sequentially continuous, we conclude that y n := y(u n ) In the remainder of this section, we explore the properties of the Cauchy problem (1), the control problem (8)-(9) and its regularization as well as the relations among them. In general, problem (8)-(9) may not admit a solution, and if solution exists it may be not unique. However, we can show the following two theorems about the solution of the control problem which are similar to the results in [24].
Theorem 2.2. Let K be the solution set of the problem (8)- (9). Suppose U ad ⊂ L 2 (Γ 2 ) is bounded and weakly closed, then K is nonempty, closed and convex in U ad . Moreover, if U ad is bounded, closed and convex, then any solution u ∈ K ⊂ U ad with corresponding state y = y(u) satisfies the following first order optimality condition Proof. The problem (8)-(9) can be equivalently rewritten in the form By the definition of infimum, there is a sequence {u n } ⊂ U ad such that The sequence {u n } ⊂ U ad is bounded, thus a subsequence can be extracted which is still denoted by {u n } such that Since U ad is weakly closed in L 2 (Γ 2 ), we have u * ∈ U ad . Note that the norm is weakly lower semicontinuous (see, [43, Thm. 1 in Chap. V]) and thus the functional J is also weakly lower semi-continuous such that which implies that u * ∈ K and K is nonempty.
Under the additional assumption that U ad is closed and convex, the solution set K is nonempty. Let u ∈ K, we havê Note that the solution operator u → y(u) is affine, we can obtain (13) by a simple calculation.
On the other hand, if u ∈ U ad is a solution of problem (8)-(9) with associated y = y(u), we have 0 ≤ J(y, u) ≤ J(ȳ,ū) = 0, so y = g D on Γ 1 and y is a solution of (1).
From the above equivalence we know that the optimal control problem (8)- (9) inherits the ill-posedness of the Cauchy problem (1). To stabilize the control problem a Tikhonov regularization term is introduced, then the regularized control problem reads in Ω, In regard to the existence and uniqueness of solution to the regularized problem (15)-(16), we have the following result.
is closed and convex, then for each fixed > 0, problem (15)-(16) admits a unique solution u . Moreover, the solution u with corresponding state y := y(u ) can be characterized by the following first order optimality condition Proof. It is evident that (15)- (16) can be equivalently rewritten in the form min u∈U adĴ Let m = min u∈U adĴ (u). For a fixed > 0, we havê If U ad is unbounded, from the above condition, there exists a constant C > 0 such thatĴ We can thus define a bounded, closed and convex subset Similar to the proof of Theorem 2.2, the existence of a solution to the problem (19) can be obtained. By the strict convexity ofĴ , the solution is unique. In fact, let u 1 , u 2 be two different solutions of (19), then 1 2 (u 1 + u 2 ) ∈ U ad and This is in contrast to the strict convexity ofĴ . The first order optimality condition (17) follows immediately from the standard arguments (see e.g., [36]), we omit it here.
We can also show that the solution u with corresponding state y of the regularized problem (15)-(16) is stable with respect to the perturbations of data f, g D and g N .
Let u n and u be the unique solutions of problem (15) - (16) with the functions f n , g D,n , g N,n and f, g D , g N , respectively. The corresponding states are y n = y(u n ) and y = y(u ), respectively. Then the following result holds.
Let us show thatū is the solution of problem (15)- (16). It follows from (17) Taking n → ∞, (21) and (23) Thus,ū is a solution of problem (15)- (16) with the functions f, g D and g N . Therefore,ū = u with corresponding stateȳ = y because that (15)-(16) has a unique solution. Since the limit u is unique, the entire sequence converges, i.e., (20) holds. Now we are in the position to investigate the solution behavior of problem (15)-(16) as → 0. The following convergence result can be proved.
Theorem 2.6. Assume that U ad ⊂ L 2 (Γ 2 ) is bounded, closed and convex. Let u be the solution of problem (15)- (16). Then where u 0 ∈ K is defined by Here K is the solution set of problem (8)-(9).
3. Finite element approximation. In this section, we assume that Ω is a polygonal or polyhedral domain for simplicity. Let T h be a triangulation of Ω such that Ω = ∪ τ ∈T hτ , and any element side on ∂Ω belongs entirely to Γ 1 or Γ 2 . We suppose that the triangulation is shape regular and quasi-uniform, such that h τ ≤ cρ τ , where h τ is the diameter of element τ ∈ T h and ρ τ the diameters of inscribed circle of τ , Then one can show that (see [9]) The present section is divided into two subsections. Firstly, we analyze the finite element approximation to the state equation. The a priori error bound for the solution of the state equation and its approximation is derived. Based on these results, we investigate the discrete scheme of the regularized optimal control problem in second subsection.
3.1. Finite element approximation to the state equation. We introduce the finite element approximation to the state equation (10) as follows. For a given function u ∈ L 2 (Γ 2 ), let y h (u) ∈ Y h be the unique solution of where The goal of this subsection is to obtain the error estimates of the approximation y h (u) given by (29) to the solution y(u) of (10). Firstly, we present two lemmas which will be used later.
Then, the following result holds for 1 < r < ∞.
Proof. The statement (32) is a direct consequence of (31) (see [19]) and the interpolation error estimates (see e.g., [9, Thm. 14.3.12]). Here we only consider (33) and (34). In view of the linearity of A, it is sufficient to consider the problems where either f ≡ 0 and g N ≡ 0 or u ≡ 0. Let us first assume that f ≡ 0 and g N ≡ 0.
From the above theorem, it is obvious that there is no convergence order under norm · 1 2 ,Ω when the solution y = y(u) ∈ H 1 2 (Ω) of (10). However, we have following convergence property which will be indispensable for further analysis. (Ω) and y h (u) ∈ Y h be the solutions of (10) and (29), respectively. Then Proof. Let L = max{L 1 , L 2 } and L 1 , L 2 are Lipschitz constants in (11), (34).

3.2.
Finite element approximation to the control problem. Now we introduce the finite element approximation to the regularized optimal control problem (15)- (16). For simplicity, the variational discretization approach developed in [25] is used to discretize the control problem. The discrete scheme takes the form (39) min Similar to Theorem 2.4, we have the following result for the discrete control problem (39)- (40).
Theorem 3.4. Suppose U ad ∈ L 2 (Γ 2 ) is closed and convex, then for each fixed > 0, problem (39)-(40) admits a unique solution u h . Moreover, the solution u h with corresponding state y h := y h (u h ) can be characterized by the following first order optimality condition The goal of this subsection is to obtain the convergence of the problem (39)- (40) to the problem (15)- (16). The main results are given in the following theorem.
Using embedding theorem, (44) can be deduced from (43). Thus, we only need to prove (43). For each fixed > 0, (42) where we have used (38) and (34). This gives the results. 4. Discrete scheme for the Cauchy problem. In this section, we discuss the numerical approximation to the Cauchy problem (1) through the discrete control problem (39)- (40). The superscript is omitted for convenience in the next two sections. To derive a numerical algorithm, further investigations of the optimality conditions (17) and (41) are required. Given a control u ∈ U ad , there is a corresponding state y(u) in control problem (8)- (9). Now, we introduce the so-called adjoint state variable p := p(u) ∈ H 1 Γ2 (Ω) as the solution of the following adjoint equation, With the help of the adjoint state equation and integration by part, it is not difficult to see that the optimality condition (17) has a simpler expression: Let P U ad denote the orthogonal projection in L 2 (Γ 2 ) onto the admissible set of controls U ad (see, [26, Lemma 1.10 in Chap. 1]). Then the variational inequality (47) can be equivalently written as a nonsmooth operator equation Similar to (46), we can introduce the discrete adjoint state variable. For discrete state variable y h = y h (u h ), let the corresponding discrete adjoint state Γ2 be the solution of the discrete adjoint equation Note that the discrete state y h and adjoint state p h are piecewise linear polynomials, so it might be not meaningful to prescribe boundary values for y h by piecewise constant functions. However, the discrete analogous of optimality condition (47) exactly proposes this, since ∂ n A * p h is a piecewise constant function on Γ 2 . In order to circumvent this difficulty and obtain the analogous expression to (47), we define a discrete normal derivative ∂ h n A * p h as follows It is obvious that ∂ h n A * p h | Γ2 ∈ V h . From the above equation and (40), we have where y h (v) is the solution of (29) with u substituted by v. It follows from (41) and (51) that the discrete optimality condition possesses the form Similar to (48), (52) is equivalent to the nonsmooth operator equation With above preparations, the solution u h ∈ U ad of problem (39)-(40) with corresponding state y h ∈ Y h and adjoint state p h ∈ Y h Γ2 can be characterized by the following system: A crucial issue in the above system is how to calculate u h , in other words, how to calculate the projection P U ad . From the viewpoint of Dirichlet boundary control problem with control constraints [11,17], the discrete problem (39)-(40) is solved on an infinite dimensional admissible control set U ad if one adopts the variational discretization concept to discretize the control variable u. It is worth emphasizing that the discrete control variable u h which depends on U ad is not usually a finite element function in most cases. Recall the discussions in section 2 and 3, to ensure the existence of a solution for the unregularized problem (8)- (9), the boundedness of U ad is indispensable (Theorem 2.2). However, this restriction is not necessary regarding the regularized scheme (15)-(16) as well as its discrete counterpart (39)- (40) (Theorem 2.4, 3.4). Therefore, there are two choices about U ad , i.e. it can be a bounded, closed, convex subset of L 2 (Γ 2 ) or U ad = L 2 (Γ 2 ). Generally speaking, it is difficult to give a concrete form of the projection P U ad if U ad is a general subset of L 2 (Γ 2 ). Fortunately, the projection has an explicit expression and is numerically implementable if U ad is of bilateral obstacle type. So in the computation of (57), we can set where a, b ∈ L ∞ (Γ 2 ). Then it is easy to verify that (see [25]) An alternative choice of the admissible control set is U ad = L 2 (Γ 2 ). In this case, P U ad is an identity operator and (57) takes the form It follows from (56), (59) and y h = Q h u h on Γ 2 that Thus, under this choice the system (54)-(57) can be rewritten as the following system: To sum up, the discrete systems (54)-(57) and (60)-(62) are the numerical schemes for the Cauchy problem. On a mesh with N nodes, we only need to solve an about 2N -order linear algebraic system to obtain the approximate solution y h of the Cauchy problem. Extensive numerical experiments show that one can get satisfactory results from the above two schemes. 5. Numerical experiments. In this section we carry out some numerical experiments to support our theoretical findings. For simplicity, the Cauchy problem for Laplace equation is considered as a model problem. Let Ω = (0, 1) 2 . The boundary Γ is split into two disjoint parts such that Γ 2 = {0} × [0, 1] and Γ 1 the rest. The model problem takes the form In the next two numerical examples, two different sets of data for f, g D , g N are described. The effect of noisy data is investigated in the second example. The admissible control set U ad is chosen as L 2 (Γ 2 ). For our computation the C++ software library AFEPack [35], supporting the creation of finite element codes, is used.
Example 5.1. This example can be found in [24]. We choose so that the model problem (63) has the unique solution The L 2 -norm error ||y −y h || 0,Ω and H 1 -norm error ||y −y h || 1,Ω between the exact solution and numerical solution for different values of mesh size h (number of mesh nodes N ) and regularization parameter are reported in TABLE 1-2. We explore the changes of errors from two aspects. Firstly, the regularization parameter is fixed. For a given , the error decreases slowly as the mesh size gets smaller if the regularization parameter has a larger value; while the trend downward accelerates if the regularization parameter is small. In particular when ≤ 10 −7 , the L 2norm error decreases with convergence order O(h 2 ) which is optimal for linear finite elements. Secondly, we investigate the error behavior when the mesh size h is fixed. For a given mesh size h, as the regularization parameter gets smaller, the error initially decreases but then increases once is smaller than a threshold value depending on h. How to select a proper regularization parameter is known to be an important and difficult problem for regularization method. In many cases, the regularization parameter is set as a matter of experience. It is worth further investigation on the choice of an optimal regularization parameter coupling with the mesh size to minimize the approximation error.
The influence of the regularization parameter and mesh size on the solution of the Cauchy problem on the under-specified part of boundary Γ 2 is illustrated in FIG. 1-4. It is shown that one can obtain a satisfactory approximate solution even on a coarse mesh (N = 2001) with = 10 −2 or 10 −3 which is sufficient for many practical applications. As the regularization parameter decreases and mesh refinement proceeds the approximate solution and exact solution are almost the same. Table 1. L 2 (Ω) error of the numerical solutions for Example 5.1.
The L 2 -norm and H 1 -norm errors of this example are given in TABLE 3-4. Similar phenomena as Example 5.1 can be observed. To further investigate the robustness of our method, the noise effect is considered in this example. We perturb g D and g N to get the noisy data g δ D and g δ N in the form g δ D (x) = g D (x) 1 + δ * rand(−1, 1) ∀ x ∈ Γ 1 , g δ N (x) = g N (x) 1 + δ * rand(−1, 1) ∀ x ∈ Γ 1 , where δ ≥ 0 denotes the noise level (the extent to which the measured data is disturbed) and rand(−1, 1) is the uniform distribution on the closed interval [−1, 1].
We mentioned the importance of choosing regularization parameter in the above example. This is particularly true for the data with noise. Extensive experiments and some references (see e.g., [8,24]) indicated that too small may lead to strong oscillation of the numerical solution and is not a good choice for the Cauchy problem in the presence of noisy data. Our numerical experiments reveal that > 10 −3 is appropriate which coincides with the results in [8,24]. We illustrate in FIG. 5-6 the influence of the regularization parameter and noise level on the numerical solution. From this figure, one can observe that the oscillations caused by noise is small. Numerical results show that the relative error is almost no more than 7% with regularization parameter = 10 −2 which is an appropriate choice. In addition, a striking scenario is worthy of note. Let y 1 h , y 2 h be numerical solutions in FIG. 5  and FIG. 6, respectively. It is apparent that y 1 h ≥ y 2 h initially, then y 1 h ≤ y 2 h , then y 1 h ≥ y 2 h and the exact solution y falls in between. In order to further reduce the error and to obtain higher accuracy, one can average these numerical solutions. The mean value of these two solutions is shown in FIG. 7. It is worth emphasizing that the method of average is often used in practical applications although the exact solution is unknown.