- Email: [email protected]
The boundary element method for the numerical recovery of a circular inhomogeneity in an elliptic equation
Engineering Analysis with Boundary Elements 28 (2004) 413–419 www.elsevier.com/locate/enganabound
The boundary element method for the numerical recovery of a circular inhomogeneity in an elliptic equation L. Marina,*, L. Elliottb, D.B. Inghamb, D. Lesnicb b
a School of the Environment, University of Leeds, Leeds LS2 9JT, UK Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK
Received 17 November 2002; revised 22 February 2003; accepted 26 April 2003
Abstract In this paper, the numerical recovery of a circular inhomogeneity in the Schro¨dinger (modified Helmholtz) equation is investigated numerically using the boundary element method, in conjunction with a least-squares minimisation. We obtain unique and stable numerical solutions. q 2003 Elsevier Ltd. All rights reserved. Keywords: Inverse problem; Schro¨dinger equation; Laplace equation; Boundary element method; Inhomogeneity
1. Introduction A direct problem in the mathematical modelling of physical systems is to determine the response of a system, provided that the governing partial differential equations, the geometry of the domain of interest, the complete boundary and initial conditions, the material properties and the external sources acting in the solution domain are given [1]. When one or more of the conditions for solving the direct problem are partially or entirely unknown then an inverse problem may be formulated to determine the unknowns from specified or measured system responses. It should be noted that most of the inverse problems are ill-posed and hence they are more difficult to solve than direct problems. It is well known that inverse problems are in general unstable [2], in the sense that small measurement errors in the input data may amplify significantly the errors in the solution. In recent years, inverse problems have been extensively treated in several branches of science, such as solid mechanics [1], heat transfer [3], acoustic and electromagnetic scattering [4], electrical impedance tomography [5], etc. The most common approach is to determine the optimal estimates of the model parameters by minimising a selected measure-to-fit between the responses of the system and the model. * Corresponding author. 0955-7997/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0955-7997(03)00088-2
An important class of inverse problems in engineering consists of the detection of defects or inhomogeneities in a given solution domain. Ikehata [5] has proposed an inversion algorithm for the electrical impedance tomography in order to recover the shape of an inclusion and the conductivity using the Dirichlet-to-Neumann map. A reconstruction formula for a three-dimensional sound-soft/sound-hard obstacle by employing the surface data of the scattering solution generated by a point source has been given by Ikehata [6]. The reconstruction of electromagnetic imperfections of small diameter from boundary measurements has been analysed in Ammari et al. [7]. Alessandrini et al. [8] have studied the identification of an inclusion embedded in an elastic material from boundary measurements of the traction and the displacement vectors and have given an estimate for the volume of the inclusion. The sensitivity analysis of structures immersed in an inviscide fluid and illuminated by harmonic incident waves has been investigated by Mallardo and Aliabadi [9]. Mallardo and Alessandri [10] have proposed a boundary element method (BEM) algorithm for identifying both deformable and rigid inclusions embedded in an elastic material. In this paper, we consider the inverse problem to recover numerically a circular inclusion in the Schro¨dinger (modified Helmholtz) equation using a single boundary measurement as input data. Referring to heat conduction applications, the problem investigated in this study consists
414
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
of determining the support of a temperature dependent heat source from a single boundary measurement of the temperature and heat flux. There have been several investigations related to the aforementioned inverse problem, which can also be motivated by the determination of the transistor contact resistivity and the contact window location corresponding to planar electronic devices [11 – 13]. A uniqueness result within a one-parameter monotone family from a one-point measurement has been obtained in Ref. [11]. Moreover, a global uniqueness result and a reconstruction scheme based on the Finite-Difference Method, in conjunction with constrained minimisation methods, within the class of two- or three-dimensional balls from a single boundary measurement have been provided by Hettlich and Rundell [12]. Further concrete practical applications arise in the determination of heat transfer coefficients, see Kraus et al. [14], and gravimetry, see Isakov [15] and the references therein. In this study, we propose a straightforward BEM to discretise the inverse problem under consideration. The BEM is further combined with a constrained least-squares minimisation routine.
2. Mathematical model We consider the mathematical model described by the following equations: ðD 2 k2 xD ÞTðxÞ ¼ 0; ~ TðxÞ ¼ TðxÞ;
FðxÞ ;
TðxÞl›D2 ¼ TðxÞl›Dþ ;
x [ V;
›T ðxÞ ¼ F~ ðxÞ; x [ ›V; ›n ›T ¼ 2 ›T ðxÞ ; ðxÞ ›Dþ 2 þ ›n ›n ›D2
Fig. 1. A schematic diagram of the domain V and the inclusion D:
In this study, we propose a straightforward BEM, as explained in Section 3, to discretise the transmission problem given by Eqs. (1) – (3). The BEM is further combined with a constrained least-squares minimisation routine. No regularisation is required since we only look for the location and the size of the circular inclusion D for which the uniqueness holds [13].
ð1Þ ð2Þ ð3Þ
in which given the data k [ R; T [ H 1=2 ð›VÞ and F~ [ H 21=2 ð›VÞ; one would like to find the inclusion D , V; where V is a bounded solution domain for the function T: In Eq. (1), xD is the characteristic function of D which is 1 inside D and 0 outside D: Also n; nþ and n2 are the outward unit normal vectors to the boundaries ›V; ›Dþ and ›D2 ; respectively, see Fig. 1. For example, referring to heat conduction applications, through the model (1) –(3) we would like to determine the support of a temperature dependent heat source, k2 xD T; ~ from a single boundary measurement of the temperature, T; and heat flux, F~ ; as given by Eq. (2). Otherwise, T could be a potential, T~ a voltage, F~ a current flux and k2 xD T a dielectric source [12]. Prior to this study, numerical methods to determine inhomogeneities in elliptic equations have been developed by Hettlich and Rundell [12] based on the Finite-Difference Method, in conjunction with constrained minimisation methods. Much closer to our study would be the recent method proposed by Kang et al. [13] which, however, requires volume integral evaluations.
3. Boundary integral formulation The problem under investigation given by Eqs. (1) and (2) and the corresponding transmission conditions (3) can be recast in a more convenient form by defining T1 ðxÞ ; TðxÞ and T2 ðxÞ ; TðxÞ for x [ D; where T1 and T2 for x [ Vw D satisfy x [ Vw D;
DT1 ðxÞ ¼ 0;
›T ~ F1 ðxÞ ; 1 ðxÞ ¼ F~ ðxÞ; T1 ðxÞ ¼ TðxÞ; ›n ðD 2 k2 ÞT2 ðxÞ ¼ 0; T1 ðxÞ ¼ T2 ðxÞ;
ð4Þ x [ ›V;
ð6Þ
x [ D;
›T1 ›T ðxÞ ¼ 2 þ2 ðxÞ; ›n2 ›n
ð5Þ
x [ ›D:
ð7Þ
Let G and E be the fundamental solutions of the Laplace equation (4) and the modified Helmholtz equation (6), respectively, which in two-dimensions are given by Ref. [16] Gðx;yxÞ ¼ 2
1 lnðlx 2 ylÞ; 2p
i Eðx;yÞ ¼ H0ð1Þ ðklx 2 ylÞ; 4
ð8Þ ð9Þ
where H0ð1Þ is the order zero Hankel function of the first kind. On applying the boundary condition (5_1) on ›V and
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
the interface transmission conditions (7) on ›D; we obtain the following boundary integral representations: ð ›Gðx;yÞ ›T ~ cðxÞT1 ðxÞ¼ Gðx;yÞ 1 ðyÞ2 dGðyÞ TðyÞ ›n ›n ›V ð ›Gðx;yÞ ›T 2 Gðx;yÞ þ2 ðyÞþ T ðyÞ dGðyÞ; ›n ›n2 2 ›D x[ Vw D; ð10Þ ð
›Eðx;yÞ ›T 2 cðxÞT2 ðxÞ ¼ Eðx;yÞ þ ðyÞ2 T ðyÞ dGðyÞ; ›n ›nþ 2 ›D x [ D; ð11Þ where cðxÞ ¼ 1 for x [ Vw ›D and cðxÞ ¼ 1=2 for x [ ›V < ›D (smooth). If we discretise uniformly the boundary ›V into N1 constant boundary elements in a counterclockwise sense and the boundary ›D into N2 constant boundary elements in both a counterclockwise and clockwise sense, and apply Eqs. (10) and (11) at the nodes on ›V < ›D and ›D; respectively, then we obtain a system of ðN1 þ 2N2 Þ nonlinear equations with ðN1 þ 4N2 Þ unknowns, namely Aðxð1Þ ; …; xðN2 Þ Þx ¼ F:
ð12Þ
Here Aðxð1Þ ; …; xðN2 Þ Þ [ RðN1 þ2N2 Þ£ðN1 þ4N2 Þ depends on the two-dimensional Cartesian coordinates of the ðjÞ boundary element nodes xðjÞ ¼ ðxðjÞ 1 ; x2 Þ; j ¼ 1; …; N2 ; on ›D; X [ RN1 þ4N2 contains the unspecified values of T1 ¼ T2 on ›D; F1 ; ›T1 =›n on ›V and F2 ; ›T2 =›nþ on ›D; and F [ RN1 þ2N2 contains the known boundary values. The resulting coefficient matrices from the BEM discretisation depend solely on the geometries of ›V and ›D and can be evaluated analytically, in the case of the boundary integral representation (10), and numerically, in the case of the boundary integral representation (11), for more details see Ref. [16]. For a given initial guess ðjÞ for the boundary element nodes xðjÞ ¼ ðxðjÞ 1 ; x2 Þ; j ¼ 1; …; N2 ; on ›D; the system of Eq. (12) becomes linear and we can solve it for the measured data F1 ; ›T1 =›n on ›V: Hence we can minimise the least-squares functional of the errors between the calculated and the measured fluxes on the outer boundary ›V; namely F ðxð1Þ ; …; xðN2 Þ Þ ¼
N1 X
ðFðnumÞ ðxðjÞ Þ 2 F~ ðxðjÞ ÞÞ2 : 1
ð13Þ
j¼1
However, even by minimising the functional (13), the problem is still underdetermined having ðN1 þ 2N2 Þ equations for ðN1 þ 3N2 Þ unknowns. Therefore, additional information is necessary in order to account for the ill-conditioned nature of the discretised inverse problem. Such constraints (additional information) may include: (i)
The inclusion D is always contained in V; i.e. xðjÞ [ V; j ¼ 1; …; N2 :
415
(ii) The introduction into the functional P 2 (13)ðjÞ of2 the penalty P 2 kx k ; l2 Nj¼1 regularising P terms, such as l1 Nj¼1 2 kxðjÞ0 k2 or l3 Nj¼1 kxðjÞ00 k2 ; where lj . 0; j ¼ 1; 2; 3; are regularisation parameters which may allow for C 0 ; C 1 or C 2 -boundaries. (iii) The boundary ›D is the union of two disjoint graphs of functions, for example, x2 ¼ f1 ðx1 Þ and x2 ¼ f2 ðx1 Þ; x1 [ ½a; b; such that the number of unknowns is reduced by N2 ; with only the components xðjÞ 1 ; j ¼ 1; …; N2 ; needed to be recovered. It should be mentioned that in the case of the problem studied in this paper, only the constraint ðiÞ has been used and no penalty regularising terms have been introduced into the objective functional, since we only look for the location and the size of a circular inclusion D:
4. Numerical results and discussion In order to test the convergence and the stability of the numerical method proposed, we assume, for simplicity, that the solution domain V is a smooth two-dimensional one, namely the disk V ¼ {x ¼ ðx1 ; x2 Þlx21 þ x22 , R2 }; R ¼ 3:0; and we restrict ourselves to the detection of a circular inclusion D ¼ {x ¼ ðx1 ; x2 Þlðx1 2 h1 Þ2 þ ðx2 2 h2 Þ2 , r 2 } embedded in V; see Fig. 1. However, it should be noted that the numerical results presented and discussed in this section remain valid even for other smooth two-dimensional domains V; such as ellipses. In this case, the inclusion D is completely determined by its location, i.e. its centre, h ¼ ðh1 ; h2 Þ; and its size, i.e. its radius, r: Consequently, the boundary element nodes ðjÞ xðjÞ ¼ ðxðjÞ 1 ; x2 Þ; j ¼ 1; …; N2 ; on ›D are parametrised by the location and the size of the inclusion, i.e. xðjÞ 1 ¼ h1 þ r cos uj and xðjÞ ¼ h þ r sin u ; j ¼ 1; …; N ; where uj is 2 j 2 2 the polar angle corresponding to the point xðjÞ : Hence the least-squares functional (13) depends only on the coordinates of the centre and the radius of the inclusion, i.e. F ¼ F ðh1 ; h2 ; rÞ: Then the following physical bounds and nonlinear constraint, respectively, are imposed on the geometrical parameters which determine the circular inclusion D : 2R , h1 , R; 2 R , h2 , R and 0 , r , R;
ð14Þ
0 , h21 þ h22 2 r 2 þ 2rR , R2 :
ð15Þ
Numerically, the least-squares functional (13) is minimised using the NAG subroutine E04UPF, which is designed to minimise an arbitrary smooth sum of squares subject to constraints. This may include simple bounds on the variables, linear constraints and smooth nonlinear constraints. Each iteration of the subroutine includes the following: (a) the solution of a quadratic programming subproblem; (b) a line search with an augmented
416
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
Lagrangian function; and (c) a quasi-Newton update of the approximate Hessian of the Lagrangian function, for more details see Ref. [17]. It should be noted that it is important to impose the constraints (14) and (15) in order to generate feasible solutions since unconstrained minimisation was found to produce physically meaningless solutions. 4.1. Example Since we assume that the solution exists, we consider the following analytical solution for the temperature distribution ( krJ 00 ðkrÞlnðrðxÞ=rÞ þ J0 ðkrÞ; x [ Vw D ðanÞ T ðxÞ ¼ ; ð16Þ J0 ðkrðxÞÞ; x[D where J0 is the Bessel function of order zero and rðxÞ ¼ lx 2 hl: The temperature (16) satisfies the Eqs. (4) –(7) and its corresponding flux on the outer boundary ›V and the interface ›D; respectively, is given by 8 ›T1ðanÞ R 0 ðx1 2 h1 Þx1 þ ðx2 2 h2 Þx2 > > ; > < ›n ðxÞ ¼ k r J 0 ðkrÞ r2 ðxÞ ðanÞ F1 ðxÞ ; > > ›T1ðanÞ > : ðxÞ ¼ kJ 00 ðkrÞ; ›n2
F2ðanÞ ðxÞ ;
›T2ðanÞ ðxÞ ¼ 2kJ 00 ðkrÞ; x [ ›D:: ›nþ
ð18Þ
F~ l›V ¼ FðexactÞ l›V are given by relation (16) and by 1 solving numerically the mixed boundary value problem (19), respectively, and in what follows they will be referred to as ‘exact data’. The numerical results presented in this section have been obtained using a discretisation of the outer boundary ›V and the inner boundary ›D corresponding to the inclusion with N1 ¼ N2 ¼ 40 boundary elements, such that N ¼ N1 þ 2N2 ¼ 120 boundary elements. These values were found to be sufficiently large such that any further refinement of the mesh size did not significantly improve the accuracy of the results. 4.2. Numerical results obtained with exact data Starting with an arbitrary initial guess for the inclusion, ðgÞ ðgÞ hðgÞ ¼ 1:5; we estimate their exact 1 ¼ 0:0; h2 ¼ 0:0 and r values, h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; using the numerical method proposed in Section 3, for exact input data and k ¼ 2:5: Fig. 2 shows the iterative convergence process as the initial guess for the inclusion D moves towards the
x [ ›V ;
ð17Þ
x [ ›D
target, for the various numbers of iterations performed. From this figure it can be seen that the numerically retrieved inclusion is a very good approximation of the target and the
It is worth noting that by considering the analytical data (16) and (17) as input data in the inverse problem under consideration, we inherently introduce numerical noise. An alternative way to fabricate the input data for the inverse problem is to solve a direct problem with complete boundary conditions and consider the numerically retrieved boundary temperature or flux on ›V as the input data for the inverse problem. More specifically, since analytical expressions for both the temperature and the flux are available in the domain V; we may solve the following mixed boundary value problem: 8 DT1 ðxÞ ¼ 0; > > > > ðanÞ > > < T1 ðxÞ ¼ T ðxÞ;
x [ Vw D x [ ›V
ðD2k2 ÞT2 ðxÞ ¼ 0; x[D ; > > > > > > : ›T1 ðxÞ ¼ 2 ›T2 ðxÞ ¼ FðanÞ ðxÞ; x [ ›D ›n2 ›nþ
ð19Þ
l›V : Hence in order to obtain the numerical flux FðexactÞ 1 the inverse problem under investigation is given by ~ ›V ¼ T ðanÞ l›V and Eqs. (4) – (7) in which the input data Tl
Fig. 2. The iterative convergence process for the circular inclusion as the ðgÞ ðgÞ initial guess ð· · ·Þ; given by hðgÞ ¼ 1:5; moves 1 ¼ 0:0; h2 ¼ 0:0 and r towards the target ð – – – Þ; i.e. the exact solution h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; for various numbers of iterations performed, namely n [ {1; 3; 6}:
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
Fig. 3. The iterative convergence process for the numerical flux FðnumÞ l›V (—) towards the analytical flux FðanÞ l›V ð – – – Þ for various numbers of iterations performed, namely n [ {1; 3; 6}:
convergence is achieved after a small number of iterations, namely n ¼ 6: Consequently, the numerical results for the location and the size of the inclusion, hðnumÞ ¼ 1:0434; 1 hðnumÞ ¼ 0:9932 and r ðnumÞ ¼ 1:0155 are in very good 2 agreement with their exact values. Fig. 3 shows the iterative convergence process for the numerical flux FðnumÞ l›V towards its exact value FðexactÞ l›V for various numbers of iterations performed, namely n [ {1; 3; 6}: It can be seen that the convergence is very rapid and at the end of the iterative process the numerical and the exact values of the flux on the outer boundary ›V are almost identical. Table 1 The values of the objective function F and the numerically retrieved values for the parameters characterising the position and the shape of the inclusion D given by h1 ¼ 1:0; h2 [ {0:25; 0:50; 0:75; 1:00; 1:25; 1:50} and r ¼ 1:0; ðgÞ ðgÞ ¼ 1:5 for k ¼ 2:5 and the initial guess hðgÞ 1 ¼ 0:0; h2 ¼ 0:0 and r Exact Numerical Error (%) F
h1 h2 R h1 h2 R h1 h2 R h1 h2 R h1 h2 R h1 h2 R
1.00 0.25 1.00 1.00 0.50 1.00 1.00 0.75 1.00 1.00 1.00 1.00 1.00 1.25 1.00 1.00 1.50 1.00
1.0430 0.2545 1.0160 1.0436 0.5018 1.0160 1.0438 0.7485 1.0163 1.0434 0.9943 1.0171 1.0423 1.2388 1.0185 1.0412 1.4837 1.0205
4.2960 1.7879 1.5974 4.3555 0.3528 1.6016 4.3755 0.2011 1.6330 4.3430 0.5674 1.7060 4.2285 0.8987 1.8456 4.1204 1.0847 2.0495
0.1048 £ 1024
No. of iterations 6
0.1273 £ 10minus;4 6
0.1726 £ 10minus;4 6
0.2639 £ 10minus;4 6
0.4657 £ 10minus;4 7
0.8803 £ 10minus;4 7
417
Table 1 presents the numerical results obtained for the parameters characterising the position and the size of the inclusion D in comparison with their exact values when the position of the centre of the inclusion is varied with respect to the x2 coordinate, while keeping the other parameters constant, i.e. h1 ¼ 1:0; r ¼ 1:0 and k ¼ 2:5: It can be seen that the numerical results retrieved using the numerical method proposed are very good approximations of their exact values, they are obtained after maximum n ¼ 7 iterations and the value of the objective function F is Oð1024 Þ in the worst case. Similar results have been obtained when varying h1 and therefore they are not presented here. In Table 2 we show the exact and the numerical results for the location and the size of the inclusion obtained using the BEM minimisation presented in Section 3 and exact input data, when the position of the inclusion and the constant k are kept constant, i.e. h1 ¼ 1:0; h2 ¼ 1:0 and k ¼ 2:5; whilst the radius of the inclusion is varied, namely r [ {0:1; 0:3; 0:5; 0:7; 0:9; 1:1; 1:3}: From this table it can be seen that better estimation of the parameters describing the position and the size of the inclusion is obtained, provided that the radius, r; of the inclusion is larger. However, it should be noted that the numerical results obtained using the proposed numerical method deteriorate as the radius, r; of the inclusion becomes very small, e.g. see the numerical results obtained in the case when h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 0:1; but they still represent reasonable estimations of their exact values. Table 2 The values of the objective function F and the numerically retrieved values for the parameters characterising the position and the shape of the inclusion D given by h1 ¼ 1:0; h2 ¼ 1:0 and r [ {0:1; 0:3; 0:5; 0:7; 0:9; 1:1; 1:3}; for ðgÞ ðgÞ ¼ 1:5 k ¼ 2:5 and the initial guess hðgÞ 1 ¼ 0:0; h2 ¼ 0:0 and r
h1 h2 r h1 h2 r h1 h2 r h1 h2 R h1 h2 r h1 h2 r h1 h2 r
Exact
Numerical
Error (%)
F
‘No. of iterations
1.00 1.00 0.10 1.00 1.00 0.30 1.00 1.00 0.50 1.00 1.00 0.70 1.00 1.00 0.90 1.00 1.00 1.10 1.00 1.00 1.30
0.9267 0.9933 0.0888 0.9314 0.9754 0.2708 0.9188 0.9751 0.4652 0.9420 0.9932 0.6778 1.0254 0.9972 0.9106 1.0286 0.9832 1.0970 0.9901 0.9767 1.2521
7.3305 0.6687 11.1575 6.8647 2.4604 9.7298 8.1194 2.4867 6.9539 5.7958 0.6787 3.1716 2.5368 0.2803 1.1812 2.8576 1.6811 0.2693 0.9943 2.3320 3.6821
0.2019 £ 1029
8
0.9622 £ 1027
7
0.2621 £ 1025
7
0.9751 £ 1025
7
0.5785 £ 1025
6
0.1251 £ 1023
8
0.2743 £ 1023
7
418
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
Table 3 The values of the objective function F and the numerically retrieved values for the parameters characterising the position and the shape of the inclusion D given by h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; for various values of k and the ðgÞ ðgÞ initial guess hðgÞ ¼ 1:5 1 ¼ 0:0; h2 ¼ 0:0 and r k
h1
Error (%)
h2
Error (%)
r
Error (%)
1.50 2.50 3.50 4.50 5.50 6.50 7.50
0.8486 1.0434 1.0048 0.9792 1.0570 0.9602 0.9618
15.1383 4.3430 0.4796 2.0850 5.6996 3.9760 3.8234
0.9754 0.9943 1.0026 1.0165 0.9887 1.0002 0.9993
2.4636 0.5674 0.2632 1.6503 1.1331 0.0210 0.0655
0.9478 1.0171 0.9636 0.9740 1.0241 0.9840 0.9831
5.2198 1.7060 2.6014 2.6014 2.4074 1.6030 1.6899
Table 3 shows the comparison of the exact and the numerical results for the location and the size of the inclusion D retrieved using the BEM algorithm described in Section 3 for exact input data, when the parameters determining the position and the size of the inclusion are kept constant, i.e. h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; and different materials are considered in D by varying the value of k in the interval [1.5,7.5]. Also in this case the numerical results obtained for the coordinates of the centre and the radius of the inclusion are in very good agreement with their exact values. It is worth noting that the larger the value of k; the more accurate are the numerical results. From Figs. 2 and 3 and Tables 1– 3 it can be concluded that the numerical method described in Section 3 is very efficient in retrieving circular inclusions of different sizes, locations and materials embedded in the domain V; provided that exact input data is available on the outer boundary of the domain under investigation. 4.3. Numerical results obtained with noisy data Exact data is seldom available in practice since measurement errors always include noise in the prescribed boundary conditions and this is studied next. In order to investigate the stability of the numerical method proposed, ~ ›V has been perturbed as the boundary data Tl ~ ›V þ dT; ~ ~T1 l›V ¼ Tl dT~ ¼ G05DDFð0; sÞ; ð20Þ ~ p ; s ¼ max lTl ›V 100 where dT is a Gaussian random variable with mean zero and standard deviation s; generated by the NAG subroutine G05DDF, and p is the percentage of additive noise included ~ ›V in order to simulate the inherent into the input data Tl measurement errors. Fig. 4 shows the exact and the numerical solution for the inclusion D given by h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; ðgÞ obtained starting with the initial guess hðgÞ 1 ¼ 0:0; h2 ¼ 0:0 ðgÞ and r ¼ 1:5; for k ¼ 2:5 and various levels of noise added
Fig. 4. The exact (—) and the numerical solution for the inclusion D; given by h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; for various levels of noise added into the temperature data Tl›V ; namely p ¼ 1% (· · ·), p ¼ 3% (– – –), p ¼ 5% ( – · · –) and p ¼ 10% ( – · · · –).
~ ›V ; namely p [ {1; 3; 5; 10}: into the input temperature Tl It can be seen that as p decreases, so the level of noise added ~ ›V also decreases, then the numerical into the input data Tl solution for the inclusion D more closely approximates the exact target, while remaining at the same time stable. The same conclusion can be drawn from Table 4 which presents the numerical results for the position and the size of the inclusion obtained with the noisy input data T~ 1 given by Eq. (20).
Table 4 The numerically retrieved values for the parameters characterising the position and the shape of the inclusion D given by h1 ¼ 1:0; h2 ¼ 1:0 and ðgÞ ðgÞ r ¼ 1:0; for k ¼ 2:5; the initial guess hðgÞ ¼ 1:5; 1 ¼ 0:0; h2 ¼ 0:0 and r and various levels p of noise added into the input temperature data T ðanÞ l›V p
h1
Error (%)
h2
Error (%)
r
Error (%)
0 1 2 3 4 5 6 7 8 9 10
1.0434 1.0482 1.0487 1.0336 1.0177 1.0005 0.9819 0.9621 0.9412 0.9200 0.9089
4.3430 4.8230 4.8712 3.3640 1.7667 0.0495 1.8053 3.7931 5.8806 7.9965 9.1066
0.9943 0.9780 0.9655 0.9595 0.9539 0.9488 0.9445 0.9411 0.9389 0.9381 0.9287
0.5674 2.2023 3.4476 4.0530 4.6112 5.1163 5.5506 5.8889 6.1061 6.1943 7.1308
1.0171 1.0250 1.0299 1.0282 1.0264 1.0244 1.0221 1.0195 1.0168 1.0140 1.0123
1.7060 2.4996 2.9873 2.8193 2.6380 2.4353 2.2060 1.9509 1.6779 1.4022 1.2346
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
5. Conclusions In this paper, the BEM has been developed for the numerical recovery of a circular inclusion in an inverse problem for the Schro¨dinger (modified Helmholtz) equation given by Eqs. (1) – (3). The BEM can also be easily implemented not only for the Schro¨dinger (modified Helmholtz) Eq. (1) with k [ R; but also for the Helmholtz equation (1) with k [ Cw R: It should be noted that the analysis performed in this study can be applied to other partial differential equations, such as convection diffusion equations, which can be reduced to Eq. (1) with k [ R or k [ Cw R [18]. From the numerical results presented and discussed in Section 4, it can be concluded that the proposed numerical method produces a stable numerical solution with respect to the location and the size of the circular inclusion D embedded in the two-dimensional domain V as the amount of noise added into the input data decreases. However, future work will involve the detection of other shapes, such as polygons for which the uniqueness of the inclusion has recently been reported by Sungwhan and Yamamoto [19], or ellipses for which no uniqueness result is as yet available in the literature.
Acknowledgements L. Marin would like to acknowledge the financial support received from EPSRC.
References [1] Kubo S. Inverse problems related to the mechanics and fracture of solids and structures. JSME Int J 1988;31:157 –66.
419
[2] Hadamard J. Lectures on Cauchy problem in linear partial differential equations. London: Oxford University Press; 1923. [3] Beck JV, Blackwell B, St Clair CR. Inverse heat conduction: ill-posed problems. New York: Wiley– Interscience; 1985. [4] Colton D, Kress R. Inverse acoustic and electromagnetic scattering. Berlin: Springer; 1992. [5] Ikehata M. How to draw a picture of an unknown inclusion from boundary measurements. Two mathematical inversion algorithms. J Inverse Ill-Posed Prob 1999;7:255 –71. [6] Ikehata M. Reconstruction of obstacle from boundary measurements. Wave Motion 1999;30:205–23. [7] Ammari H, Moskow S, Vogelius M. Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter. ESAIM Control Optim Calc Var 2003;9:49–66. [8] Alessandrini G, Morassi A, Rosset E. Detecting an inclusion in an elastic body. SIAM J Math Anal 2002;33:1247– 68. [9] Mallardo V, Aliabadi MH. A BEM sensitivity and shape identification analysis for acoustic scattering in fluid-solid problems. Int J Numer Meth Engng 1998;41:1527–41. [10] Mallardo V, Alessandri C. Inverse problems in the presence of inclusions and unilateral constraints: a boundary element approach. Comput Mech 2000;41:571 –81. [11] Fang W, Cumberbatch E. Inverse problems for metal oxide semiconductor field-effect transistor contact resistivity. SIAM J Appl Math 1992;52:699–702. [12] Hettlich F, Rundell W. Recovery of the support of a source term in an elliptic differential equation. Inverse Prob 1997;13:959–76. [13] Kang H, Kwon K, Yun K. Recovery of an inhomogeneity in an elliptic equation. Inverse Prob 2001;17:25 –44. [14] Kraus AD, Aziz A, Welty J. Extended surface heat transfer. New York: Wiley; 2001. [15] Isakov V. Inverse source problems. Providence, Rhode Island: AMS; 1990. [16] Chen G, Zhou J. Boundary element methods. London: Academic Press; 1992. [17] Gill PE, Murray W, Wright MH. Practical optimization. London: Academic Press; 1981. [18] Chen K. Improving the accuracy of DRBEM for convective partial differential equations. Engng Anal Bound Elem 1999;23: 639–44. [19] Sungwhan K, Yamamoto M. Uniqueness in identification of the support of a source term in an elliptic equation. Preprint; 2002.
The boundary element method for the numerical recovery of a circular inhomogeneity in an elliptic equation L. Marina,*, L. Elliottb, D.B. Inghamb, D. Lesnicb b
a School of the Environment, University of Leeds, Leeds LS2 9JT, UK Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK
Received 17 November 2002; revised 22 February 2003; accepted 26 April 2003
Abstract In this paper, the numerical recovery of a circular inhomogeneity in the Schro¨dinger (modified Helmholtz) equation is investigated numerically using the boundary element method, in conjunction with a least-squares minimisation. We obtain unique and stable numerical solutions. q 2003 Elsevier Ltd. All rights reserved. Keywords: Inverse problem; Schro¨dinger equation; Laplace equation; Boundary element method; Inhomogeneity
1. Introduction A direct problem in the mathematical modelling of physical systems is to determine the response of a system, provided that the governing partial differential equations, the geometry of the domain of interest, the complete boundary and initial conditions, the material properties and the external sources acting in the solution domain are given [1]. When one or more of the conditions for solving the direct problem are partially or entirely unknown then an inverse problem may be formulated to determine the unknowns from specified or measured system responses. It should be noted that most of the inverse problems are ill-posed and hence they are more difficult to solve than direct problems. It is well known that inverse problems are in general unstable [2], in the sense that small measurement errors in the input data may amplify significantly the errors in the solution. In recent years, inverse problems have been extensively treated in several branches of science, such as solid mechanics [1], heat transfer [3], acoustic and electromagnetic scattering [4], electrical impedance tomography [5], etc. The most common approach is to determine the optimal estimates of the model parameters by minimising a selected measure-to-fit between the responses of the system and the model. * Corresponding author. 0955-7997/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0955-7997(03)00088-2
An important class of inverse problems in engineering consists of the detection of defects or inhomogeneities in a given solution domain. Ikehata [5] has proposed an inversion algorithm for the electrical impedance tomography in order to recover the shape of an inclusion and the conductivity using the Dirichlet-to-Neumann map. A reconstruction formula for a three-dimensional sound-soft/sound-hard obstacle by employing the surface data of the scattering solution generated by a point source has been given by Ikehata [6]. The reconstruction of electromagnetic imperfections of small diameter from boundary measurements has been analysed in Ammari et al. [7]. Alessandrini et al. [8] have studied the identification of an inclusion embedded in an elastic material from boundary measurements of the traction and the displacement vectors and have given an estimate for the volume of the inclusion. The sensitivity analysis of structures immersed in an inviscide fluid and illuminated by harmonic incident waves has been investigated by Mallardo and Aliabadi [9]. Mallardo and Alessandri [10] have proposed a boundary element method (BEM) algorithm for identifying both deformable and rigid inclusions embedded in an elastic material. In this paper, we consider the inverse problem to recover numerically a circular inclusion in the Schro¨dinger (modified Helmholtz) equation using a single boundary measurement as input data. Referring to heat conduction applications, the problem investigated in this study consists
414
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
of determining the support of a temperature dependent heat source from a single boundary measurement of the temperature and heat flux. There have been several investigations related to the aforementioned inverse problem, which can also be motivated by the determination of the transistor contact resistivity and the contact window location corresponding to planar electronic devices [11 – 13]. A uniqueness result within a one-parameter monotone family from a one-point measurement has been obtained in Ref. [11]. Moreover, a global uniqueness result and a reconstruction scheme based on the Finite-Difference Method, in conjunction with constrained minimisation methods, within the class of two- or three-dimensional balls from a single boundary measurement have been provided by Hettlich and Rundell [12]. Further concrete practical applications arise in the determination of heat transfer coefficients, see Kraus et al. [14], and gravimetry, see Isakov [15] and the references therein. In this study, we propose a straightforward BEM to discretise the inverse problem under consideration. The BEM is further combined with a constrained least-squares minimisation routine.
2. Mathematical model We consider the mathematical model described by the following equations: ðD 2 k2 xD ÞTðxÞ ¼ 0; ~ TðxÞ ¼ TðxÞ;
FðxÞ ;
TðxÞl›D2 ¼ TðxÞl›Dþ ;
x [ V;
›T ðxÞ ¼ F~ ðxÞ; x [ ›V; ›n ›T ¼ 2 ›T ðxÞ ; ðxÞ ›Dþ 2 þ ›n ›n ›D2
Fig. 1. A schematic diagram of the domain V and the inclusion D:
In this study, we propose a straightforward BEM, as explained in Section 3, to discretise the transmission problem given by Eqs. (1) – (3). The BEM is further combined with a constrained least-squares minimisation routine. No regularisation is required since we only look for the location and the size of the circular inclusion D for which the uniqueness holds [13].
ð1Þ ð2Þ ð3Þ
in which given the data k [ R; T [ H 1=2 ð›VÞ and F~ [ H 21=2 ð›VÞ; one would like to find the inclusion D , V; where V is a bounded solution domain for the function T: In Eq. (1), xD is the characteristic function of D which is 1 inside D and 0 outside D: Also n; nþ and n2 are the outward unit normal vectors to the boundaries ›V; ›Dþ and ›D2 ; respectively, see Fig. 1. For example, referring to heat conduction applications, through the model (1) –(3) we would like to determine the support of a temperature dependent heat source, k2 xD T; ~ from a single boundary measurement of the temperature, T; and heat flux, F~ ; as given by Eq. (2). Otherwise, T could be a potential, T~ a voltage, F~ a current flux and k2 xD T a dielectric source [12]. Prior to this study, numerical methods to determine inhomogeneities in elliptic equations have been developed by Hettlich and Rundell [12] based on the Finite-Difference Method, in conjunction with constrained minimisation methods. Much closer to our study would be the recent method proposed by Kang et al. [13] which, however, requires volume integral evaluations.
3. Boundary integral formulation The problem under investigation given by Eqs. (1) and (2) and the corresponding transmission conditions (3) can be recast in a more convenient form by defining T1 ðxÞ ; TðxÞ and T2 ðxÞ ; TðxÞ for x [ D; where T1 and T2 for x [ Vw D satisfy x [ Vw D;
DT1 ðxÞ ¼ 0;
›T ~ F1 ðxÞ ; 1 ðxÞ ¼ F~ ðxÞ; T1 ðxÞ ¼ TðxÞ; ›n ðD 2 k2 ÞT2 ðxÞ ¼ 0; T1 ðxÞ ¼ T2 ðxÞ;
ð4Þ x [ ›V;
ð6Þ
x [ D;
›T1 ›T ðxÞ ¼ 2 þ2 ðxÞ; ›n2 ›n
ð5Þ
x [ ›D:
ð7Þ
Let G and E be the fundamental solutions of the Laplace equation (4) and the modified Helmholtz equation (6), respectively, which in two-dimensions are given by Ref. [16] Gðx;yxÞ ¼ 2
1 lnðlx 2 ylÞ; 2p
i Eðx;yÞ ¼ H0ð1Þ ðklx 2 ylÞ; 4
ð8Þ ð9Þ
where H0ð1Þ is the order zero Hankel function of the first kind. On applying the boundary condition (5_1) on ›V and
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
the interface transmission conditions (7) on ›D; we obtain the following boundary integral representations: ð ›Gðx;yÞ ›T ~ cðxÞT1 ðxÞ¼ Gðx;yÞ 1 ðyÞ2 dGðyÞ TðyÞ ›n ›n ›V ð ›Gðx;yÞ ›T 2 Gðx;yÞ þ2 ðyÞþ T ðyÞ dGðyÞ; ›n ›n2 2 ›D x[ Vw D; ð10Þ ð
›Eðx;yÞ ›T 2 cðxÞT2 ðxÞ ¼ Eðx;yÞ þ ðyÞ2 T ðyÞ dGðyÞ; ›n ›nþ 2 ›D x [ D; ð11Þ where cðxÞ ¼ 1 for x [ Vw ›D and cðxÞ ¼ 1=2 for x [ ›V < ›D (smooth). If we discretise uniformly the boundary ›V into N1 constant boundary elements in a counterclockwise sense and the boundary ›D into N2 constant boundary elements in both a counterclockwise and clockwise sense, and apply Eqs. (10) and (11) at the nodes on ›V < ›D and ›D; respectively, then we obtain a system of ðN1 þ 2N2 Þ nonlinear equations with ðN1 þ 4N2 Þ unknowns, namely Aðxð1Þ ; …; xðN2 Þ Þx ¼ F:
ð12Þ
Here Aðxð1Þ ; …; xðN2 Þ Þ [ RðN1 þ2N2 Þ£ðN1 þ4N2 Þ depends on the two-dimensional Cartesian coordinates of the ðjÞ boundary element nodes xðjÞ ¼ ðxðjÞ 1 ; x2 Þ; j ¼ 1; …; N2 ; on ›D; X [ RN1 þ4N2 contains the unspecified values of T1 ¼ T2 on ›D; F1 ; ›T1 =›n on ›V and F2 ; ›T2 =›nþ on ›D; and F [ RN1 þ2N2 contains the known boundary values. The resulting coefficient matrices from the BEM discretisation depend solely on the geometries of ›V and ›D and can be evaluated analytically, in the case of the boundary integral representation (10), and numerically, in the case of the boundary integral representation (11), for more details see Ref. [16]. For a given initial guess ðjÞ for the boundary element nodes xðjÞ ¼ ðxðjÞ 1 ; x2 Þ; j ¼ 1; …; N2 ; on ›D; the system of Eq. (12) becomes linear and we can solve it for the measured data F1 ; ›T1 =›n on ›V: Hence we can minimise the least-squares functional of the errors between the calculated and the measured fluxes on the outer boundary ›V; namely F ðxð1Þ ; …; xðN2 Þ Þ ¼
N1 X
ðFðnumÞ ðxðjÞ Þ 2 F~ ðxðjÞ ÞÞ2 : 1
ð13Þ
j¼1
However, even by minimising the functional (13), the problem is still underdetermined having ðN1 þ 2N2 Þ equations for ðN1 þ 3N2 Þ unknowns. Therefore, additional information is necessary in order to account for the ill-conditioned nature of the discretised inverse problem. Such constraints (additional information) may include: (i)
The inclusion D is always contained in V; i.e. xðjÞ [ V; j ¼ 1; …; N2 :
415
(ii) The introduction into the functional P 2 (13)ðjÞ of2 the penalty P 2 kx k ; l2 Nj¼1 regularising P terms, such as l1 Nj¼1 2 kxðjÞ0 k2 or l3 Nj¼1 kxðjÞ00 k2 ; where lj . 0; j ¼ 1; 2; 3; are regularisation parameters which may allow for C 0 ; C 1 or C 2 -boundaries. (iii) The boundary ›D is the union of two disjoint graphs of functions, for example, x2 ¼ f1 ðx1 Þ and x2 ¼ f2 ðx1 Þ; x1 [ ½a; b; such that the number of unknowns is reduced by N2 ; with only the components xðjÞ 1 ; j ¼ 1; …; N2 ; needed to be recovered. It should be mentioned that in the case of the problem studied in this paper, only the constraint ðiÞ has been used and no penalty regularising terms have been introduced into the objective functional, since we only look for the location and the size of a circular inclusion D:
4. Numerical results and discussion In order to test the convergence and the stability of the numerical method proposed, we assume, for simplicity, that the solution domain V is a smooth two-dimensional one, namely the disk V ¼ {x ¼ ðx1 ; x2 Þlx21 þ x22 , R2 }; R ¼ 3:0; and we restrict ourselves to the detection of a circular inclusion D ¼ {x ¼ ðx1 ; x2 Þlðx1 2 h1 Þ2 þ ðx2 2 h2 Þ2 , r 2 } embedded in V; see Fig. 1. However, it should be noted that the numerical results presented and discussed in this section remain valid even for other smooth two-dimensional domains V; such as ellipses. In this case, the inclusion D is completely determined by its location, i.e. its centre, h ¼ ðh1 ; h2 Þ; and its size, i.e. its radius, r: Consequently, the boundary element nodes ðjÞ xðjÞ ¼ ðxðjÞ 1 ; x2 Þ; j ¼ 1; …; N2 ; on ›D are parametrised by the location and the size of the inclusion, i.e. xðjÞ 1 ¼ h1 þ r cos uj and xðjÞ ¼ h þ r sin u ; j ¼ 1; …; N ; where uj is 2 j 2 2 the polar angle corresponding to the point xðjÞ : Hence the least-squares functional (13) depends only on the coordinates of the centre and the radius of the inclusion, i.e. F ¼ F ðh1 ; h2 ; rÞ: Then the following physical bounds and nonlinear constraint, respectively, are imposed on the geometrical parameters which determine the circular inclusion D : 2R , h1 , R; 2 R , h2 , R and 0 , r , R;
ð14Þ
0 , h21 þ h22 2 r 2 þ 2rR , R2 :
ð15Þ
Numerically, the least-squares functional (13) is minimised using the NAG subroutine E04UPF, which is designed to minimise an arbitrary smooth sum of squares subject to constraints. This may include simple bounds on the variables, linear constraints and smooth nonlinear constraints. Each iteration of the subroutine includes the following: (a) the solution of a quadratic programming subproblem; (b) a line search with an augmented
416
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
Lagrangian function; and (c) a quasi-Newton update of the approximate Hessian of the Lagrangian function, for more details see Ref. [17]. It should be noted that it is important to impose the constraints (14) and (15) in order to generate feasible solutions since unconstrained minimisation was found to produce physically meaningless solutions. 4.1. Example Since we assume that the solution exists, we consider the following analytical solution for the temperature distribution ( krJ 00 ðkrÞlnðrðxÞ=rÞ þ J0 ðkrÞ; x [ Vw D ðanÞ T ðxÞ ¼ ; ð16Þ J0 ðkrðxÞÞ; x[D where J0 is the Bessel function of order zero and rðxÞ ¼ lx 2 hl: The temperature (16) satisfies the Eqs. (4) –(7) and its corresponding flux on the outer boundary ›V and the interface ›D; respectively, is given by 8 ›T1ðanÞ R 0 ðx1 2 h1 Þx1 þ ðx2 2 h2 Þx2 > > ; > < ›n ðxÞ ¼ k r J 0 ðkrÞ r2 ðxÞ ðanÞ F1 ðxÞ ; > > ›T1ðanÞ > : ðxÞ ¼ kJ 00 ðkrÞ; ›n2
F2ðanÞ ðxÞ ;
›T2ðanÞ ðxÞ ¼ 2kJ 00 ðkrÞ; x [ ›D:: ›nþ
ð18Þ
F~ l›V ¼ FðexactÞ l›V are given by relation (16) and by 1 solving numerically the mixed boundary value problem (19), respectively, and in what follows they will be referred to as ‘exact data’. The numerical results presented in this section have been obtained using a discretisation of the outer boundary ›V and the inner boundary ›D corresponding to the inclusion with N1 ¼ N2 ¼ 40 boundary elements, such that N ¼ N1 þ 2N2 ¼ 120 boundary elements. These values were found to be sufficiently large such that any further refinement of the mesh size did not significantly improve the accuracy of the results. 4.2. Numerical results obtained with exact data Starting with an arbitrary initial guess for the inclusion, ðgÞ ðgÞ hðgÞ ¼ 1:5; we estimate their exact 1 ¼ 0:0; h2 ¼ 0:0 and r values, h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; using the numerical method proposed in Section 3, for exact input data and k ¼ 2:5: Fig. 2 shows the iterative convergence process as the initial guess for the inclusion D moves towards the
x [ ›V ;
ð17Þ
x [ ›D
target, for the various numbers of iterations performed. From this figure it can be seen that the numerically retrieved inclusion is a very good approximation of the target and the
It is worth noting that by considering the analytical data (16) and (17) as input data in the inverse problem under consideration, we inherently introduce numerical noise. An alternative way to fabricate the input data for the inverse problem is to solve a direct problem with complete boundary conditions and consider the numerically retrieved boundary temperature or flux on ›V as the input data for the inverse problem. More specifically, since analytical expressions for both the temperature and the flux are available in the domain V; we may solve the following mixed boundary value problem: 8 DT1 ðxÞ ¼ 0; > > > > ðanÞ > > < T1 ðxÞ ¼ T ðxÞ;
x [ Vw D x [ ›V
ðD2k2 ÞT2 ðxÞ ¼ 0; x[D ; > > > > > > : ›T1 ðxÞ ¼ 2 ›T2 ðxÞ ¼ FðanÞ ðxÞ; x [ ›D ›n2 ›nþ
ð19Þ
l›V : Hence in order to obtain the numerical flux FðexactÞ 1 the inverse problem under investigation is given by ~ ›V ¼ T ðanÞ l›V and Eqs. (4) – (7) in which the input data Tl
Fig. 2. The iterative convergence process for the circular inclusion as the ðgÞ ðgÞ initial guess ð· · ·Þ; given by hðgÞ ¼ 1:5; moves 1 ¼ 0:0; h2 ¼ 0:0 and r towards the target ð – – – Þ; i.e. the exact solution h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; for various numbers of iterations performed, namely n [ {1; 3; 6}:
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
Fig. 3. The iterative convergence process for the numerical flux FðnumÞ l›V (—) towards the analytical flux FðanÞ l›V ð – – – Þ for various numbers of iterations performed, namely n [ {1; 3; 6}:
convergence is achieved after a small number of iterations, namely n ¼ 6: Consequently, the numerical results for the location and the size of the inclusion, hðnumÞ ¼ 1:0434; 1 hðnumÞ ¼ 0:9932 and r ðnumÞ ¼ 1:0155 are in very good 2 agreement with their exact values. Fig. 3 shows the iterative convergence process for the numerical flux FðnumÞ l›V towards its exact value FðexactÞ l›V for various numbers of iterations performed, namely n [ {1; 3; 6}: It can be seen that the convergence is very rapid and at the end of the iterative process the numerical and the exact values of the flux on the outer boundary ›V are almost identical. Table 1 The values of the objective function F and the numerically retrieved values for the parameters characterising the position and the shape of the inclusion D given by h1 ¼ 1:0; h2 [ {0:25; 0:50; 0:75; 1:00; 1:25; 1:50} and r ¼ 1:0; ðgÞ ðgÞ ¼ 1:5 for k ¼ 2:5 and the initial guess hðgÞ 1 ¼ 0:0; h2 ¼ 0:0 and r Exact Numerical Error (%) F
h1 h2 R h1 h2 R h1 h2 R h1 h2 R h1 h2 R h1 h2 R
1.00 0.25 1.00 1.00 0.50 1.00 1.00 0.75 1.00 1.00 1.00 1.00 1.00 1.25 1.00 1.00 1.50 1.00
1.0430 0.2545 1.0160 1.0436 0.5018 1.0160 1.0438 0.7485 1.0163 1.0434 0.9943 1.0171 1.0423 1.2388 1.0185 1.0412 1.4837 1.0205
4.2960 1.7879 1.5974 4.3555 0.3528 1.6016 4.3755 0.2011 1.6330 4.3430 0.5674 1.7060 4.2285 0.8987 1.8456 4.1204 1.0847 2.0495
0.1048 £ 1024
No. of iterations 6
0.1273 £ 10minus;4 6
0.1726 £ 10minus;4 6
0.2639 £ 10minus;4 6
0.4657 £ 10minus;4 7
0.8803 £ 10minus;4 7
417
Table 1 presents the numerical results obtained for the parameters characterising the position and the size of the inclusion D in comparison with their exact values when the position of the centre of the inclusion is varied with respect to the x2 coordinate, while keeping the other parameters constant, i.e. h1 ¼ 1:0; r ¼ 1:0 and k ¼ 2:5: It can be seen that the numerical results retrieved using the numerical method proposed are very good approximations of their exact values, they are obtained after maximum n ¼ 7 iterations and the value of the objective function F is Oð1024 Þ in the worst case. Similar results have been obtained when varying h1 and therefore they are not presented here. In Table 2 we show the exact and the numerical results for the location and the size of the inclusion obtained using the BEM minimisation presented in Section 3 and exact input data, when the position of the inclusion and the constant k are kept constant, i.e. h1 ¼ 1:0; h2 ¼ 1:0 and k ¼ 2:5; whilst the radius of the inclusion is varied, namely r [ {0:1; 0:3; 0:5; 0:7; 0:9; 1:1; 1:3}: From this table it can be seen that better estimation of the parameters describing the position and the size of the inclusion is obtained, provided that the radius, r; of the inclusion is larger. However, it should be noted that the numerical results obtained using the proposed numerical method deteriorate as the radius, r; of the inclusion becomes very small, e.g. see the numerical results obtained in the case when h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 0:1; but they still represent reasonable estimations of their exact values. Table 2 The values of the objective function F and the numerically retrieved values for the parameters characterising the position and the shape of the inclusion D given by h1 ¼ 1:0; h2 ¼ 1:0 and r [ {0:1; 0:3; 0:5; 0:7; 0:9; 1:1; 1:3}; for ðgÞ ðgÞ ¼ 1:5 k ¼ 2:5 and the initial guess hðgÞ 1 ¼ 0:0; h2 ¼ 0:0 and r
h1 h2 r h1 h2 r h1 h2 r h1 h2 R h1 h2 r h1 h2 r h1 h2 r
Exact
Numerical
Error (%)
F
‘No. of iterations
1.00 1.00 0.10 1.00 1.00 0.30 1.00 1.00 0.50 1.00 1.00 0.70 1.00 1.00 0.90 1.00 1.00 1.10 1.00 1.00 1.30
0.9267 0.9933 0.0888 0.9314 0.9754 0.2708 0.9188 0.9751 0.4652 0.9420 0.9932 0.6778 1.0254 0.9972 0.9106 1.0286 0.9832 1.0970 0.9901 0.9767 1.2521
7.3305 0.6687 11.1575 6.8647 2.4604 9.7298 8.1194 2.4867 6.9539 5.7958 0.6787 3.1716 2.5368 0.2803 1.1812 2.8576 1.6811 0.2693 0.9943 2.3320 3.6821
0.2019 £ 1029
8
0.9622 £ 1027
7
0.2621 £ 1025
7
0.9751 £ 1025
7
0.5785 £ 1025
6
0.1251 £ 1023
8
0.2743 £ 1023
7
418
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
Table 3 The values of the objective function F and the numerically retrieved values for the parameters characterising the position and the shape of the inclusion D given by h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; for various values of k and the ðgÞ ðgÞ initial guess hðgÞ ¼ 1:5 1 ¼ 0:0; h2 ¼ 0:0 and r k
h1
Error (%)
h2
Error (%)
r
Error (%)
1.50 2.50 3.50 4.50 5.50 6.50 7.50
0.8486 1.0434 1.0048 0.9792 1.0570 0.9602 0.9618
15.1383 4.3430 0.4796 2.0850 5.6996 3.9760 3.8234
0.9754 0.9943 1.0026 1.0165 0.9887 1.0002 0.9993
2.4636 0.5674 0.2632 1.6503 1.1331 0.0210 0.0655
0.9478 1.0171 0.9636 0.9740 1.0241 0.9840 0.9831
5.2198 1.7060 2.6014 2.6014 2.4074 1.6030 1.6899
Table 3 shows the comparison of the exact and the numerical results for the location and the size of the inclusion D retrieved using the BEM algorithm described in Section 3 for exact input data, when the parameters determining the position and the size of the inclusion are kept constant, i.e. h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; and different materials are considered in D by varying the value of k in the interval [1.5,7.5]. Also in this case the numerical results obtained for the coordinates of the centre and the radius of the inclusion are in very good agreement with their exact values. It is worth noting that the larger the value of k; the more accurate are the numerical results. From Figs. 2 and 3 and Tables 1– 3 it can be concluded that the numerical method described in Section 3 is very efficient in retrieving circular inclusions of different sizes, locations and materials embedded in the domain V; provided that exact input data is available on the outer boundary of the domain under investigation. 4.3. Numerical results obtained with noisy data Exact data is seldom available in practice since measurement errors always include noise in the prescribed boundary conditions and this is studied next. In order to investigate the stability of the numerical method proposed, ~ ›V has been perturbed as the boundary data Tl ~ ›V þ dT; ~ ~T1 l›V ¼ Tl dT~ ¼ G05DDFð0; sÞ; ð20Þ ~ p ; s ¼ max lTl ›V 100 where dT is a Gaussian random variable with mean zero and standard deviation s; generated by the NAG subroutine G05DDF, and p is the percentage of additive noise included ~ ›V in order to simulate the inherent into the input data Tl measurement errors. Fig. 4 shows the exact and the numerical solution for the inclusion D given by h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; ðgÞ obtained starting with the initial guess hðgÞ 1 ¼ 0:0; h2 ¼ 0:0 ðgÞ and r ¼ 1:5; for k ¼ 2:5 and various levels of noise added
Fig. 4. The exact (—) and the numerical solution for the inclusion D; given by h1 ¼ 1:0; h2 ¼ 1:0 and r ¼ 1:0; for various levels of noise added into the temperature data Tl›V ; namely p ¼ 1% (· · ·), p ¼ 3% (– – –), p ¼ 5% ( – · · –) and p ¼ 10% ( – · · · –).
~ ›V ; namely p [ {1; 3; 5; 10}: into the input temperature Tl It can be seen that as p decreases, so the level of noise added ~ ›V also decreases, then the numerical into the input data Tl solution for the inclusion D more closely approximates the exact target, while remaining at the same time stable. The same conclusion can be drawn from Table 4 which presents the numerical results for the position and the size of the inclusion obtained with the noisy input data T~ 1 given by Eq. (20).
Table 4 The numerically retrieved values for the parameters characterising the position and the shape of the inclusion D given by h1 ¼ 1:0; h2 ¼ 1:0 and ðgÞ ðgÞ r ¼ 1:0; for k ¼ 2:5; the initial guess hðgÞ ¼ 1:5; 1 ¼ 0:0; h2 ¼ 0:0 and r and various levels p of noise added into the input temperature data T ðanÞ l›V p
h1
Error (%)
h2
Error (%)
r
Error (%)
0 1 2 3 4 5 6 7 8 9 10
1.0434 1.0482 1.0487 1.0336 1.0177 1.0005 0.9819 0.9621 0.9412 0.9200 0.9089
4.3430 4.8230 4.8712 3.3640 1.7667 0.0495 1.8053 3.7931 5.8806 7.9965 9.1066
0.9943 0.9780 0.9655 0.9595 0.9539 0.9488 0.9445 0.9411 0.9389 0.9381 0.9287
0.5674 2.2023 3.4476 4.0530 4.6112 5.1163 5.5506 5.8889 6.1061 6.1943 7.1308
1.0171 1.0250 1.0299 1.0282 1.0264 1.0244 1.0221 1.0195 1.0168 1.0140 1.0123
1.7060 2.4996 2.9873 2.8193 2.6380 2.4353 2.2060 1.9509 1.6779 1.4022 1.2346
L. Marin et al. / Engineering Analysis with Boundary Elements 28 (2004) 413–419
5. Conclusions In this paper, the BEM has been developed for the numerical recovery of a circular inclusion in an inverse problem for the Schro¨dinger (modified Helmholtz) equation given by Eqs. (1) – (3). The BEM can also be easily implemented not only for the Schro¨dinger (modified Helmholtz) Eq. (1) with k [ R; but also for the Helmholtz equation (1) with k [ Cw R: It should be noted that the analysis performed in this study can be applied to other partial differential equations, such as convection diffusion equations, which can be reduced to Eq. (1) with k [ R or k [ Cw R [18]. From the numerical results presented and discussed in Section 4, it can be concluded that the proposed numerical method produces a stable numerical solution with respect to the location and the size of the circular inclusion D embedded in the two-dimensional domain V as the amount of noise added into the input data decreases. However, future work will involve the detection of other shapes, such as polygons for which the uniqueness of the inclusion has recently been reported by Sungwhan and Yamamoto [19], or ellipses for which no uniqueness result is as yet available in the literature.
Acknowledgements L. Marin would like to acknowledge the financial support received from EPSRC.
References [1] Kubo S. Inverse problems related to the mechanics and fracture of solids and structures. JSME Int J 1988;31:157 –66.
419
[2] Hadamard J. Lectures on Cauchy problem in linear partial differential equations. London: Oxford University Press; 1923. [3] Beck JV, Blackwell B, St Clair CR. Inverse heat conduction: ill-posed problems. New York: Wiley– Interscience; 1985. [4] Colton D, Kress R. Inverse acoustic and electromagnetic scattering. Berlin: Springer; 1992. [5] Ikehata M. How to draw a picture of an unknown inclusion from boundary measurements. Two mathematical inversion algorithms. J Inverse Ill-Posed Prob 1999;7:255 –71. [6] Ikehata M. Reconstruction of obstacle from boundary measurements. Wave Motion 1999;30:205–23. [7] Ammari H, Moskow S, Vogelius M. Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter. ESAIM Control Optim Calc Var 2003;9:49–66. [8] Alessandrini G, Morassi A, Rosset E. Detecting an inclusion in an elastic body. SIAM J Math Anal 2002;33:1247– 68. [9] Mallardo V, Aliabadi MH. A BEM sensitivity and shape identification analysis for acoustic scattering in fluid-solid problems. Int J Numer Meth Engng 1998;41:1527–41. [10] Mallardo V, Alessandri C. Inverse problems in the presence of inclusions and unilateral constraints: a boundary element approach. Comput Mech 2000;41:571 –81. [11] Fang W, Cumberbatch E. Inverse problems for metal oxide semiconductor field-effect transistor contact resistivity. SIAM J Appl Math 1992;52:699–702. [12] Hettlich F, Rundell W. Recovery of the support of a source term in an elliptic differential equation. Inverse Prob 1997;13:959–76. [13] Kang H, Kwon K, Yun K. Recovery of an inhomogeneity in an elliptic equation. Inverse Prob 2001;17:25 –44. [14] Kraus AD, Aziz A, Welty J. Extended surface heat transfer. New York: Wiley; 2001. [15] Isakov V. Inverse source problems. Providence, Rhode Island: AMS; 1990. [16] Chen G, Zhou J. Boundary element methods. London: Academic Press; 1992. [17] Gill PE, Murray W, Wright MH. Practical optimization. London: Academic Press; 1981. [18] Chen K. Improving the accuracy of DRBEM for convective partial differential equations. Engng Anal Bound Elem 1999;23: 639–44. [19] Sungwhan K, Yamamoto M. Uniqueness in identification of the support of a source term in an elliptic equation. Preprint; 2002.