Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 54, 3, pp. 935-944, Warsaw 2016 DOI: 10.15632/jtam-pl.54.3.935 TREFFTZ METHOD FOR POLYNOMIAL-BASED BOUNDARY IDENTIFICATION IN TWO-DIMENSIONAL LAPLACIAN PROBLEMS Leszek Hożejowski Kielce University of Technology, Faculty of Management and Computer Modelling, Kielce, Poland e-mail: hozej@tu.kielce.pl The paper addresses a two-dimensional boundary identification (reconstruction) problem in steady-state heat conduction. Having found the solution to the Laplace equation by super- positioning T-complete functions, the unknown boundary of a plane region is approximated by polynomials of an increasing degree. The provided examples indicate that sufficient ac- curacy can be obtained with a use of polynomials of a relatively low degree, which allows avoidance of large systems of nonlinear equations. Numerical simulations for assessing the performance of the proposed algorithm show better than 1% accuracy after a few iterations and very low sensitivity to small data errors. Keywords: inverse geometry problem, boundary identification, Trefftz method, Laplace equation 1. Introduction Inverse problems appear in many areas of engineering research. They contain, among others, problems of partly unknowngeometrywherewe seek for the shape and location of the boundary of a considered domain (or a part of the boundary) and try to determine it from the information available at the known portion of the boundary. A number of these problems, referred to as inverse geometry problems (IGPs), arise in engineering practice. The intensive study of IGPs in the past decades has resulted inmany propositions of efficient solution algorithms. Although purely computational aspects are crucial to many publications in this area, researchers usually indicate possible industrial applications of their theoretical results. Illustrative examples couldbe non-destructive detection of defects like voids, cracks or inclusions (Cheng andChang, 2003a,b; Karageorghis et al., 2014), locating the 1150◦C isotherm in a blast furnace hearth (Gonzalez andGoldschmit, 2006), or non-destructive evaluation of corrosion (Lesnic et al., 2002;Mera and Lesnic, 2005; Liu, 2008) – to name but a few. The problemof numerical identification (reconstruction) of geometrical characteristics of the studied object can be solved by a variety of computational methods. For example, Hsieh and Kassab (1986) presented a general solutionmethodwhich led to a nonlinear algebraic equation, or anonlinear ordinarydifferential equation solved numerically by theNewton-Raphsonmethod. Cheng and Chang (2003a,b) used a computational technique based on the conjugate-gradient method for the recovery of inner voids within a solid body, assuming different kinds of data prescribed on the outer surface. Huang and Liu (2010) estimated interfacial configurations in a 3D composite domain by the conjugate gradient method and commercial package CFD-ACE+. Lesnic et al. (2002) treated an inverse geometric problem for the Laplace equation by the bo- undary element method in conjunction with the Tikhonov first-order regularization. Mera et al. (2004) modelled the same problem as variational and minimized a defect functional by a real coded genetic algorithm combined with the function specification method. Kazemzadeh-Parsi andDaneshmand (2009) applied the smoothed fixed grid finite elementmethod for the problem 936 L. Hożejowski of cavity detection. Liu et al. (2008) proposed a new algorithm, named the regularized integral equation method, which employed Lavrentiev regularization and a second kind Fredholm inte- gral equation to ultimately attain the unknown boundary from a nonlinear algebraic equation by the bisection method. A good reason for the choice of a method is the claim to avoid mesh generation. One could point to several numerical schemeswhichmeet this requirement andperformwell in inverse pro- blems, among them the radial basis functionsmethod and themethod of fundamental solutions (MFS), both enjoying growing interest recently, as well as theTrefftzmethod (TM) dating back to the 1920s (Trefftz, 1926). MFShas been extensively used in scientific computation and simulation over the last two de- cades.A comprehensive surveyof its application to inverse problems (including inverse geometric problems) can be found in Karageorghis et al. (2011a,b). TM and MFS, which are equivalent for Laplace and biharmonic equations (Chen et al., 2007), share the same concept of expressing the solution to a differential equation by a linear combination of basis functions satisfying the governing equation. The coefficients of the bases (T-complete functions)have tobe selected tomake the solution satisfy, in somesense, all imposed conditions. A broad discussion of the method together with a tutorial on the applications can be found in Li et al. (2008). Clearly, TM suits for solving different kinds of inverse problems (Ciałkowski andGrysa, 2010). There have been quite a few Trefftz-type approaches to a boundary identification problem. Karageorghis et al. (2014) proposed collocation TM for reconstruction of the void(s) whose centre(s) is (are) unknown, solving the proceeding non-linear least squares problem with a use of a suitable predefined routine in MATLAB. The collocation Trefftz method, modified by Liu (2007) and further referred as to modified collocation Trefftz method (MCTM), was applied by Fan et al. (2012) to a Laplacian problem to infer the spatial position of the unknownpart of the boundary from theDirichlet condition. The authors adopted the exponentially convergent scalar homotopy algorithm (ECSHA) to solve the resulting system of nonlinear equations. A similar problem, butwith theNeumann condition on the unknown part of the boundary, was solved by Fan and Chan (2011) using the same technique. A combination of MCTM and ECSHA proved successful also in boundary detection problems for a modified Helmholtz equation (Chan and Fan, 2011); Fan et al., 2014) and biharmonic equation (Chan and Fan, 2013). In this study, we propose a fast converging computational algorithm based onTM to solve a two-dimensional boundary identification problem for the Laplace equation (excluding, however, highly complicated domains). In the first stage of computationmwe approximate the solution to the Laplace equation with a linear combination of harmonic polynomials so that it satisfies the prescribed conditions in the least-squares sense. The reconstruction of an unknown boundary is through the adjustment of parameters of an approximating polynomial curve. It is demon- strated by numerical examples that relatively low degree polynomials could provide, at a small computational cost, sufficient accuracy of boundary identification. Moreover, computation with simulated errors exhibited very low sensitivity to small disturbances of the data. The approach allows immediate extension to problems governed by other linear differential equations. 2. Boundary reconstruction problem for the Laplace equation Consider a two-dimensional inverse problem in a domainΩ whose boundary ∂Ω consists of two segments:Γ , which is known, and γwhose shape andposition is not knownand, therefore, being detected. The governing equation is ∆T =0 in Ω (2.1) Trefftz method for polynomial-based boundary identification... 937 where∆ is Laplace operator, and it satisfies corresponding conditions prescribed on the doma- in boundary. The problem could be regarded as a steady-state heat conduction problem. To compensate for the missing information about the geometry of the boundary, an overspecified condition will be imposed on the known part of the boundary. Hence, we have T =TΓ and ∂T ∂n = qΓ on Γ (2.2) where ∂/∂n denotes differentiation along the outward normal n. On the remaining part of the boundary, we assume αT +β ∂T ∂n = gγ on γ (2.3) where α and β denote arbitrary constants which are not zero simultaneously. In other words, a condition of Dirichlet, Neumann or Robin type can be imposed at the boundary γ. Rather than discuss the problem in full generality, let us look at a particular situation when theproblemdomainΩ is an areaunder the curveγ. It can always beassumed that in appropriate coordinates x varies from 0 to 1, so γ can be described by an explicit equation y= f(x) x∈ 〈0,1〉 (2.4) For this particular case, the problem specified by equations (2.1)-(2.3) is presented in Fig. 1. Fig. 1. A diagram of the boundary identification problem Both boundary curvesΓ (known) and γ (unknown),meet at points (0,y0) and (1,y1), hence, we can write f(0)= y0 f(1)= y1 (2.5) As stated before, the aim of the study is to determine the shape and location of the unknown part of the boundary γ. 3. Solution methodology In the first step, we approximate the solution to Laplace equation (2.1) with a combination of the T-complete functions Vk(x,y) T(x,y)≈Θ(x,y)= K ∑ 0 ckVk(x,y) (3.1) 938 L. Hożejowski with coefficients ck to be found. In the case under consideration, the functions Vk(x,y) which satisfy the Laplace equation will be taken from the set {1,Rezk,Imzk,k=1,2, . . .} (3.2) where z=x+iy is a complex number. These basis functions are termed harmonic polynomials. Coefficients ck in expansion (3.1) will be chosen to minimize, along the boundary Γ , squ- ared differences between true values of the function T and of its normal derivative and the corresponding values provided by the model. Hence, we minimize the following functional J(c0,c1, . . . ,cK)= ∫ Γ ( K ∑ 0 ckVk ∣ ∣ ∣ Γ −TΓ )2 dΓ + ∫ Γ ( K ∑ 0 ck ∂Vk ∂n ∣ ∣ ∣ ∣ ∣ Γ −qΓ )2 dΓ (3.3) The necessary condition for the minimum of J(c0,c1, . . . ,cK) is ∂J ∂ck =0 k=0,1, . . . ,K (3.4) which gives a system ofK+1 linear equations withK+1 variables. In order to avoid problems with numerical stability, the coordinates x and y of the points ofΩ should range from0 to 1. To meet this requirement, one can use a proper change of variables when building a mathematical model of the problem. Alternatively, basis functions (3.2) could be modified by taking z/L instead of z, where L denotes the characteristic length. The latter reminds modification of the corresponding T-complete functions in polar coordinates proposed by Liu (2007). Having found the coefficients ck, we can now focus our attention directly on identification of the unknown boundary γ assumed to be a graph of the function f(x). The function f(x) can be approximated with a polynomial of degree N, denoted by P(N)(x) or shortly P(x), if there is no danger of misunderstanding. Since P(x) has the form P(x)= a0+a1x+a2x 2+ . . .+aNx N (3.5) we only need to determine the coefficients an to know the approximate values of f(x) in the whole interval 〈0,1〉. Conditions (2.6) applied to P(x) give P(0)= y0 and P(1)= y1 (3.6) In consequence, any two an can be expressed as functions of the remaining coefficients. Further proceeding must include a prescribed boundary condition (2.3) on the sought cu- rve γ. Since Θ(x,y)≈ T(x,y) and P(x)≈ f(x), we can rewrite equation (2.3) in terms of the approximated functions as αΘ(x,P(x))+β ∂Θ ∂n (x,P(x)) = gγ(x) (3.7) where ∂Θ ∂n = 1 √ 1+[P ′(x)]2 ( −P ′(x) ∂Θ ∂x (x,P(x))+ ∂Θ ∂y (x,P(x)) ) (3.8) Equation (3.7) formally holds for all x ∈ (0,1) but depending on the amount of information available on γ could be given only for x ∈ {x1,x2, . . . ,xM}. The polynomial P(x) will be selected so as to minimize the functional J(a0,a1, . . . ,aN)= 1 ∫ 0 ( αΘ(x,P(x))+β ∂Θ ∂n (x,P(x))−gγ(x) )2 dx (3.9) Trefftz method for polynomial-based boundary identification... 939 In the case of discrete data, integration contained in functional (3.9) has to be changed into summation spanning over all xm. Minimization of J(a0,a1, . . . ,aN) results in solving simultaneous equations ∂J ∂an =0 n=0,1, . . . ,N (3.10) which are nonlinear.We propose an iterative algorithm for finding the coefficients of the polyno- mial P(N)(x). The computation always starts with N =1 (i.e. a polynomial of degree 1). Note that conditions (3.6) automatically determine the first approximation P(1)(x). We now check whether boundary condition (3.7) is satisfied.More precisely, we want to estimate the goodness of fit for (3.7). For this purpose, we introduce ameasure of inaccuracy defined by Er(1) = 1 ‖gγ‖ ∥ ∥ ∥ ∥ αΘ ∣ ∣ y=P(1)(x) +β ∂Θ ∂n ∣ ∣ ∣ y=P(1)(x) −gγ ∥ ∥ ∥ ∥ ‖ · ‖ − L2 norm (3.11) assuming either continuous or discrete range of variation of x. In the latter case, both the numerator and the denominator of (3.11) contain the norm of theM-dimensional vector. If the value of Er(1) is sufficiently close to zero, the computational procedure ends and P(1)(x) is the solution. Otherwise, we increase the polynomial degree by 1 and numerically solve (3.10) for a2 (since a0 and a1 can be expressed as functions of a2). Numerical root-finding methods require a first guess which we suggest by plotting J versus a2. Once the coefficients of P(2)(x) are determined, we compute the errorEr(2) defined similarly toEr(1). In general, when solving for the coefficients of a polynomialP(N+1)(x), our initial guess are an from the previous step and aN+1 = 0. Having determined P (N+1)(x), we compute Er(N+1). The procedure is continued to achieve the desired accuracy of fulfilment of (3.7), i.e. until Er(N+1) becomes less than the desired level. In geometrical terms, the algorithm allows one to approach the unknown boundary γ with successive polynomial curves. An alternative to the stopping criterion based on Er(N), may well be the rule saying that the iteration stops when the distance between the successive approximations P(N−1)(x) and P(N)(x) is sufficiently small, which gives that either ‖P(N)(x)−P(N−1)(x)‖ or relativemeasure ‖P(N)(x)−P(N−1)(x)‖·‖P(N−1)(x)‖−1 must be less than the user-specified tolerance. We must emphasise that the proposed method, when compared to the existing approaches using MCTM (Chan and Fan, 2011, 2013; Fan and Chan, 2011; Fan et al., 2012, 2014), differs not only in the representation of the sought boundary (a continuous curve rather than a set of discrete points) but it also uses a different computation scheme. We employ two systems of equations: a linear system to determine the coefficients of T-complete functions and a nonlinear system to adjust the coefficients of a polynomial approximation, the latter having relatively few unknowns. In the cited papers, both coefficients of T-complete functions and the unknown boundarypoints are obtained fromone nonlinear systemof equations. In practice, such a system mustbe large, especiallywhenweuseahigh-orderTrefftz solution for better accuracy anda large number of boundary points for more precise boundary reconstruction. Consequently, the more boundary points we wish to determine, the larger system we obtain, unlike using the present method which allows one to recover infinitely many boundary points at the cost of solving only a small system of equations. 4. Numerical examples In this Section, some examples will be shown to test the developed theoretical method. Altho- ugh the boundary detection problem originates from realistic applications, a number of studies proposing efficient solution algorithms, among them those based on MCTM which are cited in 940 L. Hożejowski the previous paragraph, validate the proposed computationalmethods on numerical simulations (synthetic data). Likewise, the presentmethodwill be tested on data coming from numerical si- mulations.The resultswill be shown for the functionT(x,y)= e2x cos(2y),which is ananalytical solution to the Laplace equation in the domainΩ defined by inequalities 0¬x¬ 1 and 0¬ y¬ f(x) (4.1) where the graph of y= f(x) represents the unknown boundary γ. In expansion (3.1), we tookK =14, and the values of Tγ, qγ and gγ are specified atM =49 distinct points. The resulting simultaneous nonlinear equations are solved using the Levenberg- Marquardt method. For quantitative evaluation of thefinal results,we introduce ε, a non-negative numberdefined by ε= ‖P(x)−f(x)‖ ‖f(x)‖ (4.2) which can be interpreted as amean error which ismadewhen replacing the original boundary γ with its polynomial approximation P(x). Example 1 The unknown boundary is inferred from the Dirichlet condition on γ. We take f(x) = (1+x−sin4x)e−x. In Fig. 2, the true boundary is comparedwith its polynomial approximation. Fig. 2. Boundary shape identification by polynomials of degreeN =2,3,4,5 Table 1 gives information about the accuracy of the presented approximations. Trefftz method for polynomial-based boundary identification... 941 Table 1.Mean error of boundary shape identification for f(x)= (1+x− sin4x)e−x Polynomial degree N =2 N =3 N =4 N =5 Error ε 32% 2.6% 1.9% 0.1% Example 2 The computations have beenperformed for the case of theNeumann condition on γ. Figure 3 shows the results of boundary reconstruction when f(x)= (1+x2)(2+x5)−1. Fig. 3. Boundary shape identification by polynomials of degreeN =3 andN =5 The approximation errors ε, concerning the case of the Neumann condition, are listed in Table 2. Table 2.Mean error of boundary shape identification for f(x)= (1+x2)(2+x5)−1 Polynomial degree N =2 N =3 N =4 N =5 Error ε 11.9% 1.2% 1.0% 0.4% Example 3 Finally, we present the results referring to the case of the Robin condition on γ with α=1 and β = −3 in (2.3). The true boundary γ is a curve f(x) = (2+ cos4x)(x3 +2)−1 and its approximation is presented in Fig. 4. The approximation errors ε, concerning the case of the Robin boundary condition, are listed in Table 3. All examples included in this Section prove the effectiveness of the proposed computational procedure for boundary identification. It has been enough to use a 5th degree polynomial to approach the true boundary with an error less than 1%. 5. Stability analysis The computation results shown in the previous Section have been performed on the exact func- tions Tγ, qγ and gγ. In order to examine the sensitivity of the method to changes in the inputs, one has to introduce random noise to the given functions and then evaluate the impact of such changes on the final boundary detection. 942 L. Hożejowski Fig. 4. Boundary shape identification by polynomials of degreeN =3 andN =5 Table 3.Mean error of boundary shape identification for f(x)= (2+cos4x)(x3+2)−1 Polynomial degree N =2 N =3 N =4 N =5 Error ε 18.0% 8.6% 1.8% 0.4% The values of Tγ, qγ and gγ are assumed to be given only at M discrete locations. In order to simulate measurement errors, we generate M random numbers having a normal distribution with a mean of 0 and a standard deviation of σ = 0.025. With such σ, approximately 95% of the simulated errors lie between −5% and 5%. One can evaluate the influence of the introduced errors on the final results of boundary detection in a variety of ways. For the purpose of stability analysis, we use ameasure∆, defined as ∆= ‖Pnoisy(x)−P(x)‖ ‖P(x)‖ (5.1) wherePnoisy(x) denotes a polynomial approximation of the unknown boundary γ, derived from the noisy data. Since we perform a sequence of computations, each giving a successive appro- ximation to γ, the input data errors might accumulate from step to step. Therefore, it seems reasonable to pay special attention to those solution inaccuracies which occur in the last itera- tion. Consequently, our discussion of numerical stability is based on the value of ∆ calculated forN =5. For a better insight into the problem, we have recorded ∆ in 10 runs of the computatio- nal procedure, each with different random noise. The values of ∆ in the respective numerical examples are presented in Table 4. Table 4.Relative changes (∆) between the solutions derived from exact and disturbed data Example 1 Example 2 Example 3 ∆ 0.13%-0.84% 0.24%-0.88% 0.06%-0.29% It turns out that changes in the final results are even smaller than changes in the inputs, as far as percent errors are concerned. A possible explanation is that a polynomial curve, which we fit to a large number of discrete points by least squares, is forced to follow the changes introduced by random errors. Therefore, it should be shifted up at some locations but, on the other hand, shifted down at others. A low-degree polynomial is not very ‘flexible’, so it remains almost unchanged. Trefftz method for polynomial-based boundary identification... 943 6. Conclusions and final remarks We proposed a method for the solution of a two-dimensional boundary identification problem governed by the Laplace equation. In contrast to approaches which derive coordinates of unk- nown boundary points from a large system of nonlinear equations, the proposed algorithm does not require solving large systems and yet it delivers a very accurate reconstruction of the unk- nown boundary. Its basic advantage is the reduction of the amount of computational work. The provided numerical results exhibit not only high accuracy, but fast convergence of the method. Testing the sensitivity of the algorithm to input data errors showed very low risk of large im- pact of possible errors in the input data on predicted model outputs. The presented solution procedure can be applied without changes to problems governed by other linear differential equ- ations, provided the appropriate Trefftz functions are known. The only limitation concerns the geometry of a domain, because global Trefftz solutions could provide poor results when used on very complicated domains. References 1. ChanH.-F., FanC.-M., 2013,Themodified collocationTrefftzmethodand exponentially conver- gent scalarhomotopyalgorithmfor the inverseboundarydeterminationproblem for thebiharmonic equation, Journal of Mechanics, 29, 2, 363-372 2. ChanH.-F., FanC.-M., andYeihW., 2011,Solution of inverseboundary optimizationproblem by Trefftz method and exponentially convergent scalar homotopy algorithm, CMC: Computers, Materials and Continua, 24, 2, 125-142 3. Chen J.-T., Wu Y.-T., Lee Y.-T., Chen K.-H., 2007, On the equivalence of the Trefftz me- thod andmethod of fundamental solutions for Laplace and biharmonic equations,Computers and Mathematics with Applications, 53, 6, 851-879 4. Cheng C.-H., Chang M.-H., 2003a, A simplified conjugate-gradient method for shape identifi- cation based on thermal data,Numerical Heat Transfer, Part B: Fundamentals, 43, 5, 489-507 5. ChengC.-H.,ChangM.-H., 2003b,Shape identificationby inverseheat transfermethod,Journal of Heat Transfer, 125, 2, 224-231 6. Ciałkowski M., Grysa K., 2010, Trefftz method in solving the inverse problems, Journal of Inverse and Ill-Posed Problems, 18, 6, 595-616 7. Fan C.-M., Chan H.-F., 2011, Modified collocation Trefftz method for the geometry boundary identification problem of heat conduction,Numerical Heat Transfer, Part B: Fundamentals, 59, 1, 58-75 8. FanC.-M., ChanH.-F., KuoC.-L., YeihW., 2012,Numerical solutions of boundary detection problems usingmodified collocationTrefftzmethod and exponentially convergent scalar homotopy algorithm,Engineering Analysis with Boundary Elements, 36, 1, 2-8 9. FanC.-M.,LiuY.-C.,ChanH.-F.,HsiaoS.-S., 2014,Solutionsofboundarydetectionproblems for modified Helmholz equation by Trefftz method, Inverse Problems in Science and Engineering, 22, 2, 267-281 10. Gonzalez M., Goldschmit M., 2006, Inverse geometry heat transfer problem based on radial basis functions geometry representation, International Journal for Numerical Methods in Engine- ering, 65, 8, 1243-1268 11. HsiehC.K.,KassabA.J., 1986,Ageneralmethod for the solutionof inverseheat conductionpro- blemswithpartiallyunknownsystemgeometries, International Journal of Heat andMassTransfer, 29, 1, 47-58 944 L. Hożejowski 12. Huang C.-H., Liu C.-Y., 2010, A three-dimensional inverse geometry problem in estimating simultaneously two interfacial configurations in a composite domain, International Journal of Heat and Mass Transfer, 53, 1/3, 48-57 13. Karageorghis A., Lesnic D.,Marin L., 2011a,A survey of applications of theMFS to inverse problems, Inverse Problems in Science and Engineering, 19, 3, 309-336 14. Karageorghis A., Lesnic D., Marin L., 2011b,TheMFS for inverse geometric problems, [In:] Inverse Problems and Computational Mechanics, L. Munteanu, L. Marin and V. Chiroiu (Eds.), Vol. 1, Editura Academiei, Bucharest, 191-216 15. Karageorghis A., Lesnic D., Marin L., 2014,Regularized collocationTrefftzmethod for void detection in two-dimensional steady-state heat conduction problems, Inverse Problems in Science and Engineering, 22, 3, 395-418 16. Kazemzadeh-ParsiM.J., DaneshmandF., 2009, Solution of geometric inverse conduction pro- blems by smoothed fixed grid finite elementmethod,Finite Elements in Analysis and Design, 45, 10, 599-611 17. Lesnic D., Berger J.R., Martin P.A., 2002, A boundary element regularization method for the boundary determination in potential corrosion damage, Inverse Problems in Engineering, 10, 2, 163-182 18. Li Z.-C., Lu T.-T., Hu H.-Y., Cheng A. H.-D., 2008, Trefftz and Collocation Methods, WIT Press, Southampton, Boston 19. Liu C.-S., 2007, Amodified Trefftzmethod for two-dimensional Laplace equation considering the domain’s characteristic length, CMES: Computer Modeling in Engineering and Sciences, 21, 1, 53-66 20. Liu C.-S., 2008, A highly accurate MCTM for inverse Cauchy problems of Laplace equation in arbitrary plane domains,CMES: Computer Modeling in Engineering and Sciences, 35, 2, 91-111 21. Liu C.-S., Chang C.-W., Chiang C.Y., 2008, A regularized integral equation method for the inverse heat conduction problem, International Journal of Heat and Mass Transfer, 51, 21/22, 5380-5388 22. Mera N.S., Lesnic D., 2005, A three-dimensional boundary determination problem in potential corrosion damage,Computational Mechanics, 36, 2, 129-138 23. MeraN.S., EliottL., InghamD.B., 2004,Numerical solution of a boundarydetectionproblem using genetic algorithms,Engineering Analysis with Boundary Elements, 28, 4, 405-411 24. Trefftz E., 1926, Ein Gegenstck zum Ritzschen Verfahren,Verhandlungen des II. Kongress für technische Mechanik, Zürich, 131-137 Manuscript received October 29, 2015; accepted for print December 22, 2015