# Numerical Solution of Continuation Problem for 3D Steady-State Diffusion in Cylindrically Layered Medium

Authors: Magira Kulbay; Saule Maussumbekova; Balgaisha Mukanova
Journal: Abstract and Applied Analysis
 (2015)

---

## Abstract

This work is based on the application of Fourier and quasi-solution methods for solving the continuation inverse problem for 3D steady-state diffusion model inside a cylindrical layered medium. The diffusion coefficient is supposed to be a piecewise constant function, Cauchy data are given on the outer boundary of the cylinder, and we seek to recover the temperature at the inner boundary of the cylinder. Numerical experiments are investigated and show the capacity of proposed method only for smooth boundary condition. Under the suitable choice of regularization parameters we recover the distribution of temperature on the inner boundary with satisfactory quality for noisy data.


---

## Body

## 1. Introduction

The Cauchy problem for Laplace’s equation arises from many branches of science and engineering such as medicine [1, 2], geophysics [3], and nondestructive testing [4]. It is so important to reconstruct data on the boundary because the object of interest of the boundary is not available for measurement in practice. As an additional information all possible measurements are used on the available part of the boundary. Numerical computations become very difficult if there is no a priory information on the solution. A small error in the data deduces large error in the solution. Usually the regularization technique is required to find a stable approximate solution. In order to obtain a stable numerical solution for these kinds of ill-posed problems many regularization methods have been proposed [5–21]. Application of spectral methods to the Cauchy problem for Laplace’s equation was used for the first time by Lavrent’ev in [5]. A lot of regularization methods are implemented for a Cauchy problem of Laplace’s equation by many authors in [6–9]. A new regularization method based on a finite dimensional subspace generated from fundamental solution for solving a Cauchy problem of Laplace’s equation in a simply connected bounded domain is proposed in [6]. A spectral regularization method for a Cauchy problem to Laplace’s equation in a rectangle is implemented by authors of [7]. A spectral method together with choice of regularization parameter is presented and error estimate is obtained in [8]. A semidiscrete difference scheme together with a choice of regularization parameter is presented and error estimate is obtained for the Cauchy problem of Laplace’s equation in [9].The application of the method of fundamental solution with invariant condition to determine an unknown portion of the boundary in the domain from the Cauchy data connected with the Laplace equation is studied in [10].A Cauchy problem of Laplace’s equation in a multiconnected domain, that is, determining the temperature or heat flux on the inner inaccessible boundary from Cauchy data on the outer accessible boundary, is considered by authors in [11]. The method of fundamental solutions is applied to solve considered problem. They study the application of the method of fundamental solutions to solve the Cauchy problem of Laplace’s equation based on Tikhonov regularization method with L-curve choice strategy for choosing the regularization parameter.To obtain stable numerical solution for a Cauchy problem to the Laplace equation the authors of paper [12] proposed a mollification method. The idea of this method is very simple: if the data are given inexactly, then they try to find a sequence of mollification operators which map the improper data into well-posed classes of the problem.A variational method to Cauchy problem for elliptic equations is widely considered in [13, 14]. A variational-type regularized method of fundamental solution for solving Cauchy problem of Laplace’s equation in annular domain with noisy Cauchy data is studied in [14]. Under the suitable choice of regularization parameter, a convergence result for the regularized solution is obtained.In [15] the authors applied a wavelet method for solving a Cauchy problem of elliptic equations with variable coefficients. They proved the error estimates for wavelet method.Conjugate gradient method is implemented for the Cauchy problem for Laplace equation by Háo and Lesnic in [16].Application of quasi-reversibility method to solve Cauchy problem for Laplace’s equation is investigated in [17–20]. In [17] the authors described a mixed formulation of the method and its relationship with a classical formulation. In [18] the quasi-reversibility method to solve the Cauchy problem for Laplace’s equation in a smooth bounded domain is considered. They obtained two results concerning the convergence rate of the method of quasi-reversibility in the presence of noisy data.In [19] the authors introduce a new version of the method of quasi-reversibility to solve the ill-posed Cauchy problems for the Laplace’s equation in the presence of noisy data. It enables one to regularize the noisy Cauchy data and to select a relevant value of the regularization parameter in order to use the standard method of quasi-reversibility. Their method is based on duality in optimization and is inspired by Morozov’s discrepancy principle.The work [20] concerns the use of the method of quasi-reversibility for solving Cauchy problems for Laplace’s equation. The paper begins with a derivation of an error estimate for this method using Carleman’s estimates. Next, a discretization of the method using finite differences is considered. Carleman-type estimates for the discrete scheme are derived and used to establish convergence of the numerical method.In [21] the Cauchy problem for the Laplace equation in a multiply connected region was solved. Iterative algorithm for solving inverse heat conduction problems was considered. This problem was solved in a multiply connected region and solving this problem was replaced by solving the Poisson equation in a simply connected region with an unknown source function different from zero in the adjoined region. A useful feature of this approach is that the optimal step length of the line minimum search along the vector representing the gradient of the error functional can be computed exactly in the case of linear equations.In [22] the alternating iterative algorithms are used for solving Cauchy problem for elliptic equations. The algorithm works by iteratively changing boundary conditions until a satisfactory result is obtained. These methods have been applied by Kozlov et al. [23] to solve the Cauchy problem for the Laplace equation and the Lame system. The authors also proved the convergence of the algorithms and established the regularizing properties.New method of reconstruction of unknown boundary data in the Cauchy problem for Laplace’s equation is proposed by Mukanova in [24]. The problem is reduced to a system consisting of mutually dependent direct and adjoint problems. The explicit formula for the quasi-solution has been obtained using the Fourier method.In this paper the continuation inverse problem for a 3D steady-state diffusion model inside a cylindrical layered medium is considered. Methods presented in [5, 24] for the first time are used for the numerical solution of continuation problem for 3D elliptic equation with piecewise constant coefficient. They are compared on numerical examples in three dimensional cases. This model is used to describe the process of heat conduction and diffusion in protective coatings of gas and oil pipelines, the shell of reactors in the chemical industry, and other industrial processes. In such cases, only the outer side of coverings is available for measurements, when it is necessary to recover the conditions inside the coverings in order to control the equipment state.The diffusion coefficient is supposed to be a piecewise constant function, depending on radius, Cauchy data are given on the outer boundary of the cylinder, and we seek to recover the temperature at the inner boundary of the cylinder. The problem is ill-posed in the sense that small errors in the Cauchy data may lead to large errors in the recovered solution. Thus the problem is solved on the base of quasi-solution and residual functional regularization method. The system consisting of necessary minimality conditions of the residual functional and analytical formula for the minimizer is derived. The formula is numerically implemented, and series of numerical experiments are performed.The paper is organized as follows: In Section2 the statement of the problem is described. The solution of the problem is displayed in Section 3. In Section 4 numerical algorithm and examples are provided to show the effect of method.


## 2. The Statement of the Problem

LetR1, R2 be inner and outer radii of the cylinder with height H, the cylinder material is nonuniform, and the thermal conductivity is piecewise constant function of the radius. Finding the unknown temperature on the inner surface of the cylinder r=R1 is required, if the measurements of the heat flux and the temperature at the outer boundary are given. Note that all of our subsequent arguments are also valid if there are data at the inner boundary and need to find a solution on the outer boundary.Stationary temperature distribution in cylindrical coordinates is described by an elliptic equation in the form(1)Lu≡1r∂∂rrar∂u∂r+arr2∂2u∂φ2+ar∂2u∂z2=0,r,φ,z∈Qin domain(2)Q=r,φ,z∣r∈R1,R2,φ∈0,2π,z∈0,H.Here the function a(r) describes the thermal diffusivity of the medium and is assumed to be known. We consider that the case of a(r) is the piecewise constant function with a finite number of discontinuities at r=r1,r2,…,rn-1. For consistency, we assume that r0=R1, rn=R2, then, taking the thermal conductivity a(r) presenting as a(r)=ai+1, ri≤r≤ri+1, i=0,…,n-1, (1) reduces to the Laplace equation within each layer, and at the junction of layers it is required to satisfy the condition of continuity of the temperature field and heat flux:(3)Lu≡1r∂∂rr∂u∂r+1r2∂2u∂φ2+∂2u∂z2=0,r,φ,z∈Q,r≠ri,i=1,n-1¯,(4)aiur=ri-0=ai+1ur=ri+0,ai∂u∂rr=ri-0=ai+1∂u∂rr=ri+0,i=1,n-1¯with initial and boundary conditions:(5)ur,0=0,ur,H=0,(6)urn,φ,z=0,(7)∂u∂rrn,φ,z=pφ,zand periodicity conditions on the corner:(8)ur,0,z=ur,2π,z,uφr,0,z=uφr,2π,z;wherein this function on the boundary r=r0 is to be determined:(9)ur0,φ,z=qφ,z-?.So, we now formulate the continuation problem statement:It is required to find the boundary valuesq(φ,z)=u(r0,φ,z) of u(r,φ,z) satisfying (3), the boundary conditions (4), (5), and (8), and the Cauchy data (6)-(7).


## 3. The Solution of the Problem

### 3.1. The Solution of the Direct Problem

Denote byur,φ,z;q the direct problem solution, if the function qφ,z on the inner boundary is defined, that is, the solution of (3) with the conditions (4)–(6) and (8)-(9).Let us seek a solution of the direct problem in the form of double series:(10)ur,φ,z=∑k=0∞∑l=1∞uklrcos⁡kφ+tklrsin⁡kφsin⁡πlzH.The coefficientsuklr and tklr must be solutions of the following ordinary differential equations in order for the functions (10) to satisfy (3),(11)1rrukl′′-k2r2+π2l2H2ukl=0,k=0,1,…;l=1,2,…,where λl=πl/H, l=1,2,….With conjugation conditions at the internal boundaries:(12)aiukl′r=ri-0=ai+1ukl′r=ri+0,i=0,n-1¯.The problem fortkl(r) is written similarly. Equation (11) leads to a modified Bessel equation by the variable transform x=λlr:(13)∂2Rk∂x2+1x∂Rk∂x-k2x2+1Rk=0,k=0,1,2….⁡Equation (13) coincides with the modified Bessel equation [25]. A general solution of that equation is represented as a linear combination of the modified Bessel functions of kth order of first and second kind. Then the functions Rk(x) can be written in the following form:(14)Rklx=CiIkx+DiKkx,x∈λlri,λlri+1,i=0,…,n,l=1,2,…,k=0,1,…here an additional index “l” is introduced into the notation of the function Rkl(x) because it depends on the corresponding interval for variable x. Below these indexes are omitted just for brevity.We now construct an auxiliary solution of (13) satisfying the identity boundary condition at the inner boundary and the homogeneous conditions on the outside:(15)Rλlr0=CnIkλlr0+DnKkλlr0=1,Rλlrn=CnIk′λrn+DnKk′λrn=0and with the conjugation conditions on the internal boundaries r=ri(i=1,2,…,n-1):(16)CiIkλlri+DiKkλlri=Ci-1Ikλlri+Di-1Kkλlri,CiIk′λri+DiKk′λlri=aiai+1Ci-1Ik′λlri+Di-1Kk′λlri.Conditions (15)-(16) define the 2n arbitrary constants C(i), D(i) for each pair of indices of harmonic numbers (k,l) in formulas of (14). Then coefficients of the series (14) for each pair (k,l) are presented as a product of undetermined coefficients and auxiliary solutions:(17)uklr=ARklλlr,tklr=BRklλlr.By construction, the partial sums of series (14) with coefficients in the form (17) satisfy (11) and coupling conditions (12). Now, obtain expressions for the undetermined coefficients A, B in (17). Boundary data (9) on the inner boundary are expanded in a Fourier series by calculating the coefficients qkl and tkl as follows:(18)qφ,z=ur0,φ,z=∑k=0∞∑l=1∞qklcos⁡kφ+tklsin⁡kφsin⁡πlzH.We put in formula (17) undetermined coefficients that are equal to the values:(19)A=qkl,B=tkl.Then the solution of the direct problem is represented as follows:(20)ur,φ,z=∑k=0∞∑l=1∞Rklλlrqklcos⁡kφ+tklsin⁡kφsin⁡λlz.It can be shown that the functions Rkl(λlr) are uniformly bounded on the interval (r0,rn), so series (20) converges, as soon as series (18) converge.


### 3.2. The Solution of the Inverse Problem

Obviously, because of (15) and (19), the partial sums of (20) satisfy all the conditions of the direct problem. We suppose also the function q(φ,z) such that series (20) converges well with their second derivatives, so that its term-by-term differentiation of  (20) is permissible to calculate the normal derivative of the function at the outer boundary at r=rn:(21)∂u∂rr,φ,zr=rn=∂∂r∑k=0K∑l=1Luklrcos⁡kφ+tklrsin⁡kφsin⁡λlzr=rn=∑k=0∞∑l=1∞λlRk′λlrnqklcos⁡kφ+tklsin⁡kφsin⁡λlz.Expand the measured data (7) in a Fourier series:(22)pφ,z=∑k=0∞∑l=1∞pklcos⁡kφ+sklsin⁡kφsin⁡λlz.Equating (21) and (22), because of the orthogonality of the basic functions, we see that the Fourier coefficients of the measured data and boundary values at the inner boundary are related linearly via the following formulas:(23)pkl=γklqkl,skl=γkltkl,qkl=ηklpkl,tkl=ηklskl,where(24)γkl=λlRkl′λlrn,ηkl=γkl-1.Formulas (23) represent the discrete analogue of operator of inverse problem:(25)Λqφ,z=pφ,z.This operator is a diagonal and self-adjoint in discrete form; that is, the Fourier coefficient of the unknown function is multiplied by a diagonal matrix with elements γkl. Then the inverse operator is also recorded in a diagonal form for the discrete case:(26)qkl=γkl-1pkl,tkl=γkl-1skl.The coefficientsγkl=λlRk′(λlrn) do not depend on the measured data and can be calculated independently for each problem with different parameters. In addition, depending on the data of the problem, we can estimate what harmonic coefficients of the function q(φ,z) can be recovered.Let us estimate the values ofγkl. Due to the asymptotic properties of Bessel functions the coefficients γkl have the asymptotics by k,l→∞ as follows:(27)Kkλlrn~exp⁡-λlrnλlrn,Ikλlrn~exp⁡λlrnλlrnand also their first derivatives:(28)Kk′λlrn=-Kk-1λlrn-kKkλlrn~max⁡exp⁡-λlrnl,kexp⁡-λlrnll,Ik′λlrn=Ik-1λlrn-kIlλlrn~max⁡exp⁡λlrnl,kexp⁡λlrnll.By virtue of (14) and (20) the calculation error of Fourier coefficients of measured data in the recovered function q(φ,z) is multiplied by ηkl, which is inverse to γkl. Consider(29)ηkl~max⁡exp⁡λlrnl,kexp⁡λlrnll.The total recovery error inL2 norm is equal to the sum of the squared errors of the Fourier coefficients:(30)ε=∑k,lK,Lηkl2δpkl2+δtkl2.As the coefficientsηkl grow exponentially, errors of high-order harmonics give the greatest contribution to the recovery error. This implies a limit on the number of harmonics that can be recovered with the desired accuracy. Let the relative accuracy of restoring of the harmonics with the number (k,l) be less than the preset number δqkl/qkl≤ε. Then it follows from (26) that the inequality δqkl=ηklδpklqkl≤ε should be satisfied.Hence(31)δpkl≤εγkl=Cεmax⁡exp⁡-λlrnl,kexp⁡-λlrnll.Execution of (31) is possible, if we restrict ourselves to some finite number of harmonics, that is, as the set of admissible q(φ,z) is taken as functions which are linear combination of this number of harmonics. This is justified, if the measured data are very smooth, so that the contribution of higher harmonics is negligible.Restriction (31) gives us the conditions for choosing of the regularization parameter like the maximum possible number of harmonics K, L.It follows that one can neglect the harmonics, which will be submitted in the form of measured data with an amplitude which satisfies the following estimate:(32)exp⁡λlrnlpkl<ε.After calculating functionsur,φ,z by (20) and (26) the partial sum of series (14) for r=r0 can be calculated and the unknown function q(φ,z) is found approximately as(33)qφ,z=ur0,φ,z=∑k=0K∑l=1Lqklcos⁡kφ+tklsin⁡kφsin⁡πlzH.To get the Cauchy data for this problem, we need to generate synthetic data—the solution of the direct problem with a given functionqφ,z. First of all we have to solve the direct problem with a known qφ,z, then to obtain the values pφ,z. After that the inverse problem is solved, as if we had unknown function q(φ,z), and the function p(φ,z) had been measured. The direct problem is easily solved by Fourier method as it is described in Section 3.1. The only feature is that at the inner boundary r=ri the conditions of conjugation should be satisfied.In the case of a discontinuous coefficient of thermal conductivity, we use the representation of the solution in each layer similar to (22):(34)uklλr=CiIkλr+DiKkλr,r∈ri,ri+1,i=0,n¯,λ=πlH,l=1,…,L,k=0,…K.The conjugation conditions on the inner layers with r=ri are(35)CiIkλlri+DiKkλlri=Ci-1Ikλlri+Di-1Kkλlri,ai+1CiIk′λlri+DiKk′λlri=aiCi-1Ik′λlri+Di-1Kk′λlri,and the boundary conditions at r=r0 and r=rn are the following:(36)C0Ikλlr0+D0Kkλlr0=qkl,CnIkλlrn+1+DnKkλlrn+1=0.System (35)-(36) defines constants for each layer and allows one to build the coefficients of series (14) in form (34) for each layer. Obtaining the solution of the direct problem, we also obtain measured data and synthetic function psyn(φ,z) in the form(37)psynφ,z=∂∂r∑k=0K∑l=1Luklrcos⁡kφ+tklrsin⁡kφsin⁡πlzHr=r0=∑k=0K∑l=1Lpklcos⁡kφ+sklsin⁡kφsin⁡πlzH.The relative difference of synthetic data in normL20,2π×0,H for different harmonics number has been calculated for given q(φ,z). It turns out that for cases K=L=32 and K=L=50 it was found to be above 0.18%. Further increase of the cut-off parameters K and L almost does not affect that value. Then to generate synthetic data we eliminate these values by K=L=50.


### 3.3. Construction of the Quasi-Solution to the Inverse Problem

Another way to obtain an explicit formula for the quasi-solution of problem (3)–(9) is to derive it from minimality conditions of the residual functional:(38)Jq=12∫02π∫0HurR2,z;q-pφ,z2dzdφ+12β∫02π∫0Hq2φ,zdzdφ⟶minq·,·.The idea of that method is described in detail in [24], so we describe here only key points of that method. Let us introduce formally a functional (Lagrangian):(39)Lq,v=Jq-∭Qvr,φ,zLur,φ,zdQ,where ur,φ,z is a solution of the direct problem for given qφ,z. Let the function vr,φ,z be a solution of the following adjoint problem:(40)1r∂∂rr∂v∂r+1r2∂2v∂φ2+∂2v∂z2=0,r,φ,z∈Q,r≠ri,i=1,n-1¯subject to boundary and periodicity conditions:(41)vr=ri-0=vr=ri+0,ai∂v∂rr=ri-0=ai+1∂v∂rr=ri+0,i=1,n-1¯,(42)vr,φ,0=0,vr,φ,H=0,vr0,φ,z=0,rnanvrn,φ,z=urrn,z-pz,vr,0,z=vr,2π,z,vφr,0,z=vφr,2π,z.Now calculate the first variation of the residual using the variation of the Lagrangian:(43)δJq=δLq,v=δJq-∭Qvr,φ,zLδur,φ,zdQ=∬ΩurR2,z;q-pφ,zδur+βqδqdΩ-∬Ω∫r0rnv∂∂rarrδurdrdΩ-∭QvLφ,zδurdQ,where Lφ,z is a part of operator (1) acting with respect to variables φ and z. Let us consider second term of expression (43). Integrating by part and taking into account conditions at discontinuity points r=ri of coefficient ar, we obtain(44)∫r0rnv∂∂rarrδurdr=rvarδurr=rn-rvarδurr=r0+∑i=1n-1rivri,φ,zaiδurr=ri-0r=ri+0-∫r0rnrvrarδurdr.Here arδurri-0ri+0=ai+1δurr=ri+0-aiδurr=ri-0 is the jump of the function at r=ri. By conditions (4) these values vanish. Then(45)∫r0rnv∂∂rarrδurdr=anrnvδurr=rn-a0r0vδurr=r0-anrnvrδur=rn+a0r0vrδur=r0-∑i=1n-1riuri,φ,zaiδvrr=ri-0r=ri+0+∫r0rnδu∂∂rarrvrdr.Using boundary conditions (6) and (42) we have(46)∫r0rnv∂∂rarrδurdr=anrnvδurr=rn+a0r0vrr0,φ,zδq+∫r0rnδu∂∂rarrvrdr.Finally, substituting (46) into expression (43), integrating by part over Ω, and applying conditions for the function vr,φ,z we get(47)δJq=∬Ωurrn,φ,z;q-pφ,zδur+βqδq-rnanvrn,φ,zδur-a0r0vrr0,φ,zδqdΩ-∭QδuLvdQ=∬Ωβq-a1r0vrr0,φ,zδqdΩ.Note that functional (38) is represented in the form of a sum of two squared norms in the space L2([0,2π]×[0,H]):(48)Jq=12Λq-p2+12βqz2.First item of (48) is a squared norm of a linear self-adjoint operator; second one is a strongly convex square functional. Then, by [26] functional (38) is strongly convex and has a unique minimizer in L2([0,2π]×[0,H]). This minimizer can be found via necessary minimality condition that requires that the first variation at the minimum point should vanish. Then, by equalizing to zero the first variation (47), we get(49)qφ,z=β-1a1r0vr0,φ,z≡β1-1vr0,φ,z,β1=βa1r0.Therefore, if we substitute expression (49) into the direct problem and take into account the adjoint problem, we get necessary minimality conditions of residual (38). These conditions are expressed as the following coupled system consisting of direct and adjoint problems subject to corresponding boundary conditions:(50)1r∂∂rr∂u∂r+1r2∂2u∂φ2+∂2u∂z2=0,1r∂∂rr∂v∂r+1r2∂2v∂φ2+∂2v∂z2=0,r,φ,z∈Q,r≠ri,i=1,n-1¯,ur=ri-0=ur=ri+0,ai∂u∂rr=ri-0=ai+1∂u∂rr=ri+0,i=1,n-1¯,vr=ri-0=vr=ri+0,ai∂v∂rr=ri-0=ai+1∂v∂rr=ri+0,i=1,n-1¯,ur,0,z=ur,2π,z,uφr,0,z=uφr,2π,z,ur0,φ,z=1β1vrr0,φ,z,urn,φ,z=0,vr,0,z=vr,2π,z,vφr,0,z=vφr,2π,z,anrnvrn,φ,z=urrn,z-pz.Let us apply the Fourier method to solve system (50). Represent the functions ur,φ,z, vr,φ,z in the form of series similar to (10):(51)ur,φ,z=∑k=0∞∑l=1∞uklrcos⁡kφ+tklrsin⁡kφsin⁡πlzH≡∑j=0∞ujrΦjφ,z,(52)vr,φ,z=∑j=0∞vjrΦjφ,z.Here for simplicity the set of basic functions is enumerated anew sequentially:(53)cos⁡kφsin⁡λlz,cos⁡kφsin⁡λlz,k=0,1,…,l=1,2,…≡Φjφ,z,j=0,1,…,so that each pair of k,l corresponds to the unique index j and vice versa.Let us firstly find functionsujr and vjr formally and then prove the convergence of the series (51) and (52). Substituting series (51) and (52) into boundary conditions in (50) we get the following equalities for coefficients:(54)ujr0=β1-1vj′r0,anvjrn=uj′rn-pj.As it is shown in Section3.1 functions uj(r) satisfy (11) and relations (23); that is, uj′rn,z=ηjujr0. Then we obtain (55)anrnvjrn=ηjujr0-pj.Consider now the functionsvjr. They satisfy the same equation as (11), which can be transformed in the form of Bessel equation (13). Now we consider its auxiliary solutions Rj¯x subject to the following boundary conditions:(56)Rj¯λlr0=0,Rj¯λlrn=1and the conjugation conditions (16) on the internal boundaries r=ri(i=1,2,…,n-1). Then coefficients vjr of series (52) for each index j are presented as a product of undetermined coefficients and auxiliary solutions:(57)vjr=VjRj¯λlr.Then first derivative at inner boundary is calculated as(58)vj′r0=λlRkl′¯λlr0Vj≡μjVj=μjvjrn.Let us compare values μj=λlRkl′¯λlr0 and ηj=λlRkl′λlrn. By direct computations one can see that the variable transforms (59)x1=x-0.5λlrn+r0,x2=0.5λlrn+r0-xapplied to the functions Rj and Rj¯, respectively, lead to the same equation and the same boundary conditions for the transformed functions y1,jx1=Rj(x1+0.5λlrn+r0) and y2,j(x2)=Rj¯(0.5λlrn+r0-x2). Then by unicity y1,jx1≡y2,jx2.Then(60)μkl=λlddrRkl¯λlr0=-λldy2,jdx2x2=0.5λlrn-r0,ηkl=λlddrRklλlrn=λldy1,jdx1x1=0.5λlrn-r0.Because y1,jx1≡y2,jx2, this yields to equality ηkl=-μkl. Therefore (61)ujr0=β1-1vj′r0=β1-1μjvjrn=β1-1μjηjujr0-pj=-β1-1ηjηjujr0-pj.This yields(62)ujr0=ηjpjβ+ηj2.Finally, the approximate regularized solution of the considered inverse problem is obtained in explicit form:(63)qβφ,z=∑j=0∞ηjpjβ+ηj2Φjφ,z=∑k=0∞∑l=1∞ηklpklβ+ηkl2cos⁡kφ+ηklsklβ+ηkl2sin⁡kφsin⁡πlzH,where ηkl are defined by expression (24). Via standard analysis of the function x/(β1+x2) we get an estimate ηj/(β1+ηj2)≤β1-1/2. Then the convergence of series (63) follows from the estimate(64)qβφ,z2=∑j=0∞ηj2pj2β+ηj22Φj2≤1β1∑j=0∞pj2Φj2=1β1pφ,z2<∞.It follows from (64) that the regularized problem is well-posed, because minimizer qβφ,z of functional (38) exists, is unique, and depends continuously on input data pφ,z.Since the direct problem with functionqβφ,z at inner boundary is well-posed that verifies presentation (51) of the solution in the form of Fourier series, as well as solution (52) of adjoint problem.As one can see, here three regularization parameters appear, namely, harmonic numbersK, L and a parameter β. The solutions for both methods, which are described in Sections 3.2 and 3.3, are calculated numerically and results are compared in Table 2.


## 3.1. The Solution of the Direct Problem

Denote byur,φ,z;q the direct problem solution, if the function qφ,z on the inner boundary is defined, that is, the solution of (3) with the conditions (4)–(6) and (8)-(9).Let us seek a solution of the direct problem in the form of double series:(10)ur,φ,z=∑k=0∞∑l=1∞uklrcos⁡kφ+tklrsin⁡kφsin⁡πlzH.The coefficientsuklr and tklr must be solutions of the following ordinary differential equations in order for the functions (10) to satisfy (3),(11)1rrukl′′-k2r2+π2l2H2ukl=0,k=0,1,…;l=1,2,…,where λl=πl/H, l=1,2,….With conjugation conditions at the internal boundaries:(12)aiukl′r=ri-0=ai+1ukl′r=ri+0,i=0,n-1¯.The problem fortkl(r) is written similarly. Equation (11) leads to a modified Bessel equation by the variable transform x=λlr:(13)∂2Rk∂x2+1x∂Rk∂x-k2x2+1Rk=0,k=0,1,2….⁡Equation (13) coincides with the modified Bessel equation [25]. A general solution of that equation is represented as a linear combination of the modified Bessel functions of kth order of first and second kind. Then the functions Rk(x) can be written in the following form:(14)Rklx=CiIkx+DiKkx,x∈λlri,λlri+1,i=0,…,n,l=1,2,…,k=0,1,…here an additional index “l” is introduced into the notation of the function Rkl(x) because it depends on the corresponding interval for variable x. Below these indexes are omitted just for brevity.We now construct an auxiliary solution of (13) satisfying the identity boundary condition at the inner boundary and the homogeneous conditions on the outside:(15)Rλlr0=CnIkλlr0+DnKkλlr0=1,Rλlrn=CnIk′λrn+DnKk′λrn=0and with the conjugation conditions on the internal boundaries r=ri(i=1,2,…,n-1):(16)CiIkλlri+DiKkλlri=Ci-1Ikλlri+Di-1Kkλlri,CiIk′λri+DiKk′λlri=aiai+1Ci-1Ik′λlri+Di-1Kk′λlri.Conditions (15)-(16) define the 2n arbitrary constants C(i), D(i) for each pair of indices of harmonic numbers (k,l) in formulas of (14). Then coefficients of the series (14) for each pair (k,l) are presented as a product of undetermined coefficients and auxiliary solutions:(17)uklr=ARklλlr,tklr=BRklλlr.By construction, the partial sums of series (14) with coefficients in the form (17) satisfy (11) and coupling conditions (12). Now, obtain expressions for the undetermined coefficients A, B in (17). Boundary data (9) on the inner boundary are expanded in a Fourier series by calculating the coefficients qkl and tkl as follows:(18)qφ,z=ur0,φ,z=∑k=0∞∑l=1∞qklcos⁡kφ+tklsin⁡kφsin⁡πlzH.We put in formula (17) undetermined coefficients that are equal to the values:(19)A=qkl,B=tkl.Then the solution of the direct problem is represented as follows:(20)ur,φ,z=∑k=0∞∑l=1∞Rklλlrqklcos⁡kφ+tklsin⁡kφsin⁡λlz.It can be shown that the functions Rkl(λlr) are uniformly bounded on the interval (r0,rn), so series (20) converges, as soon as series (18) converge.


## 3.2. The Solution of the Inverse Problem

Obviously, because of (15) and (19), the partial sums of (20) satisfy all the conditions of the direct problem. We suppose also the function q(φ,z) such that series (20) converges well with their second derivatives, so that its term-by-term differentiation of  (20) is permissible to calculate the normal derivative of the function at the outer boundary at r=rn:(21)∂u∂rr,φ,zr=rn=∂∂r∑k=0K∑l=1Luklrcos⁡kφ+tklrsin⁡kφsin⁡λlzr=rn=∑k=0∞∑l=1∞λlRk′λlrnqklcos⁡kφ+tklsin⁡kφsin⁡λlz.Expand the measured data (7) in a Fourier series:(22)pφ,z=∑k=0∞∑l=1∞pklcos⁡kφ+sklsin⁡kφsin⁡λlz.Equating (21) and (22), because of the orthogonality of the basic functions, we see that the Fourier coefficients of the measured data and boundary values at the inner boundary are related linearly via the following formulas:(23)pkl=γklqkl,skl=γkltkl,qkl=ηklpkl,tkl=ηklskl,where(24)γkl=λlRkl′λlrn,ηkl=γkl-1.Formulas (23) represent the discrete analogue of operator of inverse problem:(25)Λqφ,z=pφ,z.This operator is a diagonal and self-adjoint in discrete form; that is, the Fourier coefficient of the unknown function is multiplied by a diagonal matrix with elements γkl. Then the inverse operator is also recorded in a diagonal form for the discrete case:(26)qkl=γkl-1pkl,tkl=γkl-1skl.The coefficientsγkl=λlRk′(λlrn) do not depend on the measured data and can be calculated independently for each problem with different parameters. In addition, depending on the data of the problem, we can estimate what harmonic coefficients of the function q(φ,z) can be recovered.Let us estimate the values ofγkl. Due to the asymptotic properties of Bessel functions the coefficients γkl have the asymptotics by k,l→∞ as follows:(27)Kkλlrn~exp⁡-λlrnλlrn,Ikλlrn~exp⁡λlrnλlrnand also their first derivatives:(28)Kk′λlrn=-Kk-1λlrn-kKkλlrn~max⁡exp⁡-λlrnl,kexp⁡-λlrnll,Ik′λlrn=Ik-1λlrn-kIlλlrn~max⁡exp⁡λlrnl,kexp⁡λlrnll.By virtue of (14) and (20) the calculation error of Fourier coefficients of measured data in the recovered function q(φ,z) is multiplied by ηkl, which is inverse to γkl. Consider(29)ηkl~max⁡exp⁡λlrnl,kexp⁡λlrnll.The total recovery error inL2 norm is equal to the sum of the squared errors of the Fourier coefficients:(30)ε=∑k,lK,Lηkl2δpkl2+δtkl2.As the coefficientsηkl grow exponentially, errors of high-order harmonics give the greatest contribution to the recovery error. This implies a limit on the number of harmonics that can be recovered with the desired accuracy. Let the relative accuracy of restoring of the harmonics with the number (k,l) be less than the preset number δqkl/qkl≤ε. Then it follows from (26) that the inequality δqkl=ηklδpklqkl≤ε should be satisfied.Hence(31)δpkl≤εγkl=Cεmax⁡exp⁡-λlrnl,kexp⁡-λlrnll.Execution of (31) is possible, if we restrict ourselves to some finite number of harmonics, that is, as the set of admissible q(φ,z) is taken as functions which are linear combination of this number of harmonics. This is justified, if the measured data are very smooth, so that the contribution of higher harmonics is negligible.Restriction (31) gives us the conditions for choosing of the regularization parameter like the maximum possible number of harmonics K, L.It follows that one can neglect the harmonics, which will be submitted in the form of measured data with an amplitude which satisfies the following estimate:(32)exp⁡λlrnlpkl<ε.After calculating functionsur,φ,z by (20) and (26) the partial sum of series (14) for r=r0 can be calculated and the unknown function q(φ,z) is found approximately as(33)qφ,z=ur0,φ,z=∑k=0K∑l=1Lqklcos⁡kφ+tklsin⁡kφsin⁡πlzH.To get the Cauchy data for this problem, we need to generate synthetic data—the solution of the direct problem with a given functionqφ,z. First of all we have to solve the direct problem with a known qφ,z, then to obtain the values pφ,z. After that the inverse problem is solved, as if we had unknown function q(φ,z), and the function p(φ,z) had been measured. The direct problem is easily solved by Fourier method as it is described in Section 3.1. The only feature is that at the inner boundary r=ri the conditions of conjugation should be satisfied.In the case of a discontinuous coefficient of thermal conductivity, we use the representation of the solution in each layer similar to (22):(34)uklλr=CiIkλr+DiKkλr,r∈ri,ri+1,i=0,n¯,λ=πlH,l=1,…,L,k=0,…K.The conjugation conditions on the inner layers with r=ri are(35)CiIkλlri+DiKkλlri=Ci-1Ikλlri+Di-1Kkλlri,ai+1CiIk′λlri+DiKk′λlri=aiCi-1Ik′λlri+Di-1Kk′λlri,and the boundary conditions at r=r0 and r=rn are the following:(36)C0Ikλlr0+D0Kkλlr0=qkl,CnIkλlrn+1+DnKkλlrn+1=0.System (35)-(36) defines constants for each layer and allows one to build the coefficients of series (14) in form (34) for each layer. Obtaining the solution of the direct problem, we also obtain measured data and synthetic function psyn(φ,z) in the form(37)psynφ,z=∂∂r∑k=0K∑l=1Luklrcos⁡kφ+tklrsin⁡kφsin⁡πlzHr=r0=∑k=0K∑l=1Lpklcos⁡kφ+sklsin⁡kφsin⁡πlzH.The relative difference of synthetic data in normL20,2π×0,H for different harmonics number has been calculated for given q(φ,z). It turns out that for cases K=L=32 and K=L=50 it was found to be above 0.18%. Further increase of the cut-off parameters K and L almost does not affect that value. Then to generate synthetic data we eliminate these values by K=L=50.


## 3.3. Construction of the Quasi-Solution to the Inverse Problem

Another way to obtain an explicit formula for the quasi-solution of problem (3)–(9) is to derive it from minimality conditions of the residual functional:(38)Jq=12∫02π∫0HurR2,z;q-pφ,z2dzdφ+12β∫02π∫0Hq2φ,zdzdφ⟶minq·,·.The idea of that method is described in detail in [24], so we describe here only key points of that method. Let us introduce formally a functional (Lagrangian):(39)Lq,v=Jq-∭Qvr,φ,zLur,φ,zdQ,where ur,φ,z is a solution of the direct problem for given qφ,z. Let the function vr,φ,z be a solution of the following adjoint problem:(40)1r∂∂rr∂v∂r+1r2∂2v∂φ2+∂2v∂z2=0,r,φ,z∈Q,r≠ri,i=1,n-1¯subject to boundary and periodicity conditions:(41)vr=ri-0=vr=ri+0,ai∂v∂rr=ri-0=ai+1∂v∂rr=ri+0,i=1,n-1¯,(42)vr,φ,0=0,vr,φ,H=0,vr0,φ,z=0,rnanvrn,φ,z=urrn,z-pz,vr,0,z=vr,2π,z,vφr,0,z=vφr,2π,z.Now calculate the first variation of the residual using the variation of the Lagrangian:(43)δJq=δLq,v=δJq-∭Qvr,φ,zLδur,φ,zdQ=∬ΩurR2,z;q-pφ,zδur+βqδqdΩ-∬Ω∫r0rnv∂∂rarrδurdrdΩ-∭QvLφ,zδurdQ,where Lφ,z is a part of operator (1) acting with respect to variables φ and z. Let us consider second term of expression (43). Integrating by part and taking into account conditions at discontinuity points r=ri of coefficient ar, we obtain(44)∫r0rnv∂∂rarrδurdr=rvarδurr=rn-rvarδurr=r0+∑i=1n-1rivri,φ,zaiδurr=ri-0r=ri+0-∫r0rnrvrarδurdr.Here arδurri-0ri+0=ai+1δurr=ri+0-aiδurr=ri-0 is the jump of the function at r=ri. By conditions (4) these values vanish. Then(45)∫r0rnv∂∂rarrδurdr=anrnvδurr=rn-a0r0vδurr=r0-anrnvrδur=rn+a0r0vrδur=r0-∑i=1n-1riuri,φ,zaiδvrr=ri-0r=ri+0+∫r0rnδu∂∂rarrvrdr.Using boundary conditions (6) and (42) we have(46)∫r0rnv∂∂rarrδurdr=anrnvδurr=rn+a0r0vrr0,φ,zδq+∫r0rnδu∂∂rarrvrdr.Finally, substituting (46) into expression (43), integrating by part over Ω, and applying conditions for the function vr,φ,z we get(47)δJq=∬Ωurrn,φ,z;q-pφ,zδur+βqδq-rnanvrn,φ,zδur-a0r0vrr0,φ,zδqdΩ-∭QδuLvdQ=∬Ωβq-a1r0vrr0,φ,zδqdΩ.Note that functional (38) is represented in the form of a sum of two squared norms in the space L2([0,2π]×[0,H]):(48)Jq=12Λq-p2+12βqz2.First item of (48) is a squared norm of a linear self-adjoint operator; second one is a strongly convex square functional. Then, by [26] functional (38) is strongly convex and has a unique minimizer in L2([0,2π]×[0,H]). This minimizer can be found via necessary minimality condition that requires that the first variation at the minimum point should vanish. Then, by equalizing to zero the first variation (47), we get(49)qφ,z=β-1a1r0vr0,φ,z≡β1-1vr0,φ,z,β1=βa1r0.Therefore, if we substitute expression (49) into the direct problem and take into account the adjoint problem, we get necessary minimality conditions of residual (38). These conditions are expressed as the following coupled system consisting of direct and adjoint problems subject to corresponding boundary conditions:(50)1r∂∂rr∂u∂r+1r2∂2u∂φ2+∂2u∂z2=0,1r∂∂rr∂v∂r+1r2∂2v∂φ2+∂2v∂z2=0,r,φ,z∈Q,r≠ri,i=1,n-1¯,ur=ri-0=ur=ri+0,ai∂u∂rr=ri-0=ai+1∂u∂rr=ri+0,i=1,n-1¯,vr=ri-0=vr=ri+0,ai∂v∂rr=ri-0=ai+1∂v∂rr=ri+0,i=1,n-1¯,ur,0,z=ur,2π,z,uφr,0,z=uφr,2π,z,ur0,φ,z=1β1vrr0,φ,z,urn,φ,z=0,vr,0,z=vr,2π,z,vφr,0,z=vφr,2π,z,anrnvrn,φ,z=urrn,z-pz.Let us apply the Fourier method to solve system (50). Represent the functions ur,φ,z, vr,φ,z in the form of series similar to (10):(51)ur,φ,z=∑k=0∞∑l=1∞uklrcos⁡kφ+tklrsin⁡kφsin⁡πlzH≡∑j=0∞ujrΦjφ,z,(52)vr,φ,z=∑j=0∞vjrΦjφ,z.Here for simplicity the set of basic functions is enumerated anew sequentially:(53)cos⁡kφsin⁡λlz,cos⁡kφsin⁡λlz,k=0,1,…,l=1,2,…≡Φjφ,z,j=0,1,…,so that each pair of k,l corresponds to the unique index j and vice versa.Let us firstly find functionsujr and vjr formally and then prove the convergence of the series (51) and (52). Substituting series (51) and (52) into boundary conditions in (50) we get the following equalities for coefficients:(54)ujr0=β1-1vj′r0,anvjrn=uj′rn-pj.As it is shown in Section3.1 functions uj(r) satisfy (11) and relations (23); that is, uj′rn,z=ηjujr0. Then we obtain (55)anrnvjrn=ηjujr0-pj.Consider now the functionsvjr. They satisfy the same equation as (11), which can be transformed in the form of Bessel equation (13). Now we consider its auxiliary solutions Rj¯x subject to the following boundary conditions:(56)Rj¯λlr0=0,Rj¯λlrn=1and the conjugation conditions (16) on the internal boundaries r=ri(i=1,2,…,n-1). Then coefficients vjr of series (52) for each index j are presented as a product of undetermined coefficients and auxiliary solutions:(57)vjr=VjRj¯λlr.Then first derivative at inner boundary is calculated as(58)vj′r0=λlRkl′¯λlr0Vj≡μjVj=μjvjrn.Let us compare values μj=λlRkl′¯λlr0 and ηj=λlRkl′λlrn. By direct computations one can see that the variable transforms (59)x1=x-0.5λlrn+r0,x2=0.5λlrn+r0-xapplied to the functions Rj and Rj¯, respectively, lead to the same equation and the same boundary conditions for the transformed functions y1,jx1=Rj(x1+0.5λlrn+r0) and y2,j(x2)=Rj¯(0.5λlrn+r0-x2). Then by unicity y1,jx1≡y2,jx2.Then(60)μkl=λlddrRkl¯λlr0=-λldy2,jdx2x2=0.5λlrn-r0,ηkl=λlddrRklλlrn=λldy1,jdx1x1=0.5λlrn-r0.Because y1,jx1≡y2,jx2, this yields to equality ηkl=-μkl. Therefore (61)ujr0=β1-1vj′r0=β1-1μjvjrn=β1-1μjηjujr0-pj=-β1-1ηjηjujr0-pj.This yields(62)ujr0=ηjpjβ+ηj2.Finally, the approximate regularized solution of the considered inverse problem is obtained in explicit form:(63)qβφ,z=∑j=0∞ηjpjβ+ηj2Φjφ,z=∑k=0∞∑l=1∞ηklpklβ+ηkl2cos⁡kφ+ηklsklβ+ηkl2sin⁡kφsin⁡πlzH,where ηkl are defined by expression (24). Via standard analysis of the function x/(β1+x2) we get an estimate ηj/(β1+ηj2)≤β1-1/2. Then the convergence of series (63) follows from the estimate(64)qβφ,z2=∑j=0∞ηj2pj2β+ηj22Φj2≤1β1∑j=0∞pj2Φj2=1β1pφ,z2<∞.It follows from (64) that the regularized problem is well-posed, because minimizer qβφ,z of functional (38) exists, is unique, and depends continuously on input data pφ,z.Since the direct problem with functionqβφ,z at inner boundary is well-posed that verifies presentation (51) of the solution in the form of Fourier series, as well as solution (52) of adjoint problem.As one can see, here three regularization parameters appear, namely, harmonic numbersK, L and a parameter β. The solutions for both methods, which are described in Sections 3.2 and 3.3, are calculated numerically and results are compared in Table 2.


## 4. Numerical Algorithm and Examples

### 4.1. Generation of Stochastic Noise

In practice, the data can never be measured precisely, due to their discreteness and measurement errors. Therefore, a simulation of noise in order to test our numerical method should be included into data (37). A noise is added as the sum of harmonics with random amplitudes. Then the noisy data for noise γ look like(65)pγφ,z=∑k=0K∑l=1Lpkl+dklcos⁡kφ+skl+eklsin⁡kφsin⁡πlzH,where dkl, ekl are random variables uniformly distributed on the interval -δ,δ. As a measure of the relative noise level we will take the value:(66)γ=pγ-psynpsyn≈2π∑l=1Ld0l2+∑k,l=1K,Ldkl2+ekl22π∑l=1Lp0l2+∑k,l=1K,Lpkl2+skl2.


### 4.2. Algorithm of Solving the Inverse Problem

The algorithm consists of the following steps:(i)
Step 1. It is necessary to set the number of harmonicsK, L over the angles φ and z and regularization parameter β. Generate synthetic and noisy data. Expand the boundary measurement with noise (38) to Fourier series by calculating the coefficients pkl,skl of series (37), k=0,…,K; l=1,…,L.(ii)
Step 2. Put in (36) value 1.0 instead of qkl. Solve the linear system (35)-(36) and define the coefficients C(i), D(i) for each pair of indexes (k,l) and i=0,…,n-1.(iii)
Step 3. Calculate first derivatives of functionsRkl(λlr)=C(0)Ik(λlr)+D(0)Kk(λlr) at r=r0 and find values ηkl defined by (24): ηkl=λlRkl′(λlrn)-1.(iv)
Step 4. Calculate the valuesqkl=ηklpkl, tkl=ηklskl.(v)
Step 5. Calculate the partial sum of series (33), which is the required approximate solution of the problem.Note that if parameterβ in formula (63) vanishes, then by (24) the solution is reduced to the solution obtained in Section 3.2. Therefore one can use the above algorithm for both cases with β=0 and β>0.


### 4.3. Examples of Numerical Calculation

We present here a numerical example to illustrate the performance of our algorithm. The thermal diffusivity coefficient is given by formula(67)ar=a1,1≤r<r1,a2,r1≤r<r2,a3,r2≤r≤r3.Example 1.
We reconstruct the smooth boundary condition as(68)fnzexp⁡-z2,q_exactφ,z=fnz-0.250.1+3fnz-0.750.1sin⁡πzHfnφ-π2+fnφ-3π2sin⁡φ2+fnz-0.252+φ-π20.15+3fnz-0.752+φ-π20.15,with the following parameters: K=16, L=16, H=1, a1=20, a2=1, a3=20, r1=1.25, r2=1.7, r3=2. The numerical reconstructions are obtained in this case for various percentages of noise γnoise∈0,2,5,8% are shown in Figures 1(b), 1(c), 1(d), and 1(e), respectively. L2 relative errors were computed to measure the quality of reconstructed boundary condition as follows:(69)εq=q_exactφ,z-q_recoveredφ,zL20,2π×0,Hq_exactφ,zL20,2π×0,H.
The comparison of relative errorsεq for different noise levels γnoise and admissible regularization parameter β for Example 1 is presented in Table 1.
The numbers that are highlighted in bold show preferable values of regularization parameterβ for a given noise level γnoise. In particular it follows from table that a large value of the noise level requires a larger value for the parameter of regularization.Table 1
Comparison of relative errors for different noise levelsγnoise and admissible regularization parameter β.
β
/
γ
noise
0%
2%
5%
8%
1
0
-
15
0.19
0.196
0.196
0.208
1
0
-
16
0.171
0.177
0.189
0.205
1
0
-
17
0.159
0.16
0.175
0.166
1
0
-
18
0.136
0.14
0.196
0.161
1
0
-
19
0.132
0.139
0.137
0.214
1
0
-
20
0.13
0.131
0.148
0.184
1
0
-
21
0.128
0.133
0.142
0.183Table 2
The values of relative errorsεq depending on different values of parameters K and L, with corresponding admissible values of β>0 and β=0 with noise 5%.
K
/
L
32
16
12
8
5
β
1
×
1
0
-
25
1
×
1
0
-
25
1
×
1
0
-
24
1
×
1
0
-
25
1
×
1
0
-
25
ε
q
0.05
0.08
0.07
0.12
0.22
β
0
0
0
0
0
ε
q
0.05
0.069
0.08
0.13
0.23Figure 1
Reconstruction of smooth boundary condition with various noise levelsγnoise∈0,2,5,8%.
(a)
Exact function
(b)
γ
noise
=
0
%, β=1⋅10-19, εq=0.132
(c)
γ
noise
=
2
%, β=1⋅10-19, εq=0.139
(d)
γ
noise
=
5
%, β=1⋅10-19, εq=0.137
(e)
γ
noise
=
8
%, β=1⋅10-19, εq=0.214Example 2.
In this example we presented the impact of harmonic numbers on reconstruction with noise 5% for parameters:H=1, a1=20, a2=1, a3=20, r1=1.25, r2=1.7, r3=2.
We attempt to reconstruct smooth boundary condition of the form:(70)fnz=exp⁡-z2,q_exactφ,z=fnz-0.250.11+3fnz-0.750.11sin⁡πzHfnφ-π2sin⁡φ2.The results for this example are illustrated in Figure 2.
From Figures1 and 2 we can see that when 5≤K, L≤32, we need to take β close to machine precision. According to Table 2 the results for β=0 and β=1×10-25 are similar with the same harmonics numbers K=32, L=32. However, the numerical experiments show that use of small value of the regularization parameter β≈1×10-25 is preferable to β=0.Figure 2
The influence of harmonic number on reconstruction.
(a)
K
=
L
=
32
(b)
K
=
L
=
16
(c)
K
=
L
=
12
(d)
K
=
L
=
8
(e)
K
=
L
=
5


## 4.1. Generation of Stochastic Noise

In practice, the data can never be measured precisely, due to their discreteness and measurement errors. Therefore, a simulation of noise in order to test our numerical method should be included into data (37). A noise is added as the sum of harmonics with random amplitudes. Then the noisy data for noise γ look like(65)pγφ,z=∑k=0K∑l=1Lpkl+dklcos⁡kφ+skl+eklsin⁡kφsin⁡πlzH,where dkl, ekl are random variables uniformly distributed on the interval -δ,δ. As a measure of the relative noise level we will take the value:(66)γ=pγ-psynpsyn≈2π∑l=1Ld0l2+∑k,l=1K,Ldkl2+ekl22π∑l=1Lp0l2+∑k,l=1K,Lpkl2+skl2.


## 4.2. Algorithm of Solving the Inverse Problem

The algorithm consists of the following steps:(i)
Step 1. It is necessary to set the number of harmonicsK, L over the angles φ and z and regularization parameter β. Generate synthetic and noisy data. Expand the boundary measurement with noise (38) to Fourier series by calculating the coefficients pkl,skl of series (37), k=0,…,K; l=1,…,L.(ii)
Step 2. Put in (36) value 1.0 instead of qkl. Solve the linear system (35)-(36) and define the coefficients C(i), D(i) for each pair of indexes (k,l) and i=0,…,n-1.(iii)
Step 3. Calculate first derivatives of functionsRkl(λlr)=C(0)Ik(λlr)+D(0)Kk(λlr) at r=r0 and find values ηkl defined by (24): ηkl=λlRkl′(λlrn)-1.(iv)
Step 4. Calculate the valuesqkl=ηklpkl, tkl=ηklskl.(v)
Step 5. Calculate the partial sum of series (33), which is the required approximate solution of the problem.Note that if parameterβ in formula (63) vanishes, then by (24) the solution is reduced to the solution obtained in Section 3.2. Therefore one can use the above algorithm for both cases with β=0 and β>0.


## 4.3. Examples of Numerical Calculation

We present here a numerical example to illustrate the performance of our algorithm. The thermal diffusivity coefficient is given by formula(67)ar=a1,1≤r<r1,a2,r1≤r<r2,a3,r2≤r≤r3.Example 1.
We reconstruct the smooth boundary condition as(68)fnzexp⁡-z2,q_exactφ,z=fnz-0.250.1+3fnz-0.750.1sin⁡πzHfnφ-π2+fnφ-3π2sin⁡φ2+fnz-0.252+φ-π20.15+3fnz-0.752+φ-π20.15,with the following parameters: K=16, L=16, H=1, a1=20, a2=1, a3=20, r1=1.25, r2=1.7, r3=2. The numerical reconstructions are obtained in this case for various percentages of noise γnoise∈0,2,5,8% are shown in Figures 1(b), 1(c), 1(d), and 1(e), respectively. L2 relative errors were computed to measure the quality of reconstructed boundary condition as follows:(69)εq=q_exactφ,z-q_recoveredφ,zL20,2π×0,Hq_exactφ,zL20,2π×0,H.
The comparison of relative errorsεq for different noise levels γnoise and admissible regularization parameter β for Example 1 is presented in Table 1.
The numbers that are highlighted in bold show preferable values of regularization parameterβ for a given noise level γnoise. In particular it follows from table that a large value of the noise level requires a larger value for the parameter of regularization.Table 1
Comparison of relative errors for different noise levelsγnoise and admissible regularization parameter β.
β
/
γ
noise
0%
2%
5%
8%
1
0
-
15
0.19
0.196
0.196
0.208
1
0
-
16
0.171
0.177
0.189
0.205
1
0
-
17
0.159
0.16
0.175
0.166
1
0
-
18
0.136
0.14
0.196
0.161
1
0
-
19
0.132
0.139
0.137
0.214
1
0
-
20
0.13
0.131
0.148
0.184
1
0
-
21
0.128
0.133
0.142
0.183Table 2
The values of relative errorsεq depending on different values of parameters K and L, with corresponding admissible values of β>0 and β=0 with noise 5%.
K
/
L
32
16
12
8
5
β
1
×
1
0
-
25
1
×
1
0
-
25
1
×
1
0
-
24
1
×
1
0
-
25
1
×
1
0
-
25
ε
q
0.05
0.08
0.07
0.12
0.22
β
0
0
0
0
0
ε
q
0.05
0.069
0.08
0.13
0.23Figure 1
Reconstruction of smooth boundary condition with various noise levelsγnoise∈0,2,5,8%.
(a)
Exact function
(b)
γ
noise
=
0
%, β=1⋅10-19, εq=0.132
(c)
γ
noise
=
2
%, β=1⋅10-19, εq=0.139
(d)
γ
noise
=
5
%, β=1⋅10-19, εq=0.137
(e)
γ
noise
=
8
%, β=1⋅10-19, εq=0.214Example 2.
In this example we presented the impact of harmonic numbers on reconstruction with noise 5% for parameters:H=1, a1=20, a2=1, a3=20, r1=1.25, r2=1.7, r3=2.
We attempt to reconstruct smooth boundary condition of the form:(70)fnz=exp⁡-z2,q_exactφ,z=fnz-0.250.11+3fnz-0.750.11sin⁡πzHfnφ-π2sin⁡φ2.The results for this example are illustrated in Figure 2.
From Figures1 and 2 we can see that when 5≤K, L≤32, we need to take β close to machine precision. According to Table 2 the results for β=0 and β=1×10-25 are similar with the same harmonics numbers K=32, L=32. However, the numerical experiments show that use of small value of the regularization parameter β≈1×10-25 is preferable to β=0.Figure 2
The influence of harmonic number on reconstruction.
(a)
K
=
L
=
32
(b)
K
=
L
=
16
(c)
K
=
L
=
12
(d)
K
=
L
=
8
(e)
K
=
L
=
5


## 5. Conclusion

In this study we come to conclusion that the application of the proposed approach gives an efficient way to solve the Cauchy problem for the 3D Laplace equation in a cylindrical layered medium. Numerical experiments are investigated and show the capacity of proposed method for smooth boundary condition only. Under the suitable choice of regularization parameter we recover the distribution of temperature on the inner boundary with satisfactory quality for noisy data. We have determined more admissible values of the main parameters such as a regularization parameter and harmonic numbers for different noise level.



---

*Source: 329052-2015-09-10.xml*