Abstract and Applied Analysis

Volume 2015, Article ID 329052, 13 pages

http://dx.doi.org/10.1155/2015/329052

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

^{1}Department of Mathematical and Computing Modelling, Al-Farabi Kazakh National University, Almaty 050040, Kazakhstan^{2}Faculty of Mathematics and Mechanics, L.N.Gumilyov Eurasian National University, Astana 010008, Kazakhstan

Received 23 May 2015; Revised 17 August 2015; Accepted 18 August 2015

Academic Editor: Sergei V. Pereverzyev

Copyright © 2015 Magira Kulbay et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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 -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 Section 2 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

Let , be inner and outer radii of the cylinder with height , 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 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 formin domainHere the function describes the thermal diffusivity of the medium and is assumed to be known. We consider that the case of is the piecewise constant function with a finite number of discontinuities at . For consistency, we assume that , , then, taking the thermal conductivity presenting as , , , (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:with initial and boundary conditions:and periodicity conditions on the corner:wherein this function on the boundary is to be determined:So, we now formulate the continuation problem statement:

It is required to find the boundary values of 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 by the direct problem solution, if the function 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:

The coefficients and must be solutions of the following ordinary differential equations in order for the functions (10) to satisfy (3),where , .

With conjugation conditions at the internal boundaries:

The problem for is written similarly. Equation (11) leads to a modified Bessel equation by the variable transform :

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 th order of first and second kind. Then the functions can be written in the following form:here an additional index “” is introduced into the notation of the function because it depends on the corresponding interval for variable . 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:and with the conjugation conditions on the internal boundaries :Conditions (15)-(16) define the arbitrary constants , for each pair of indices of harmonic numbers in formulas of (14). Then coefficients of the series (14) for each pair are presented as a product of undetermined coefficients and auxiliary solutions:

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 , in (17). Boundary data (9) on the inner boundary are expanded in a Fourier series by calculating the coefficients and as follows:We put in formula (17) undetermined coefficients that are equal to the values:Then the solution of the direct problem is represented as follows:It can be shown that the functions are uniformly bounded on the interval , 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 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 :Expand the measured data (7) in a Fourier series: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:whereFormulas (23) represent the discrete analogue of operator of inverse problem: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 . Then the inverse operator is also recorded in a diagonal form for the discrete case:

The coefficients 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 can be recovered.

Let us estimate the values of . Due to the asymptotic properties of Bessel functions the coefficients have the asymptotics by as follows:and also their first derivatives:By virtue of (14) and (20) the calculation error of Fourier coefficients of measured data in the recovered function is multiplied by , which is inverse to . Consider

The total recovery error in norm is equal to the sum of the squared errors of the Fourier coefficients:

As the coefficients 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 be less than the preset number . Then it follows from (26) that the inequality should be satisfied.

Hence

Execution of (31) is possible, if we restrict ourselves to some finite number of harmonics, that is, as the set of admissible 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 , .

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:

After calculating functions by (20) and (26) the partial sum of series (14) for can be calculated and the unknown function is found approximately as

To get the Cauchy data for this problem, we need to generate synthetic data—the solution of the direct problem with a given function . First of all we have to solve the direct problem with a known , then to obtain the values . After that the inverse problem is solved, as if we had unknown function , and the function 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 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):The conjugation conditions on the inner layers with areand the boundary conditions at and are the following:

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 in the form

The relative difference of synthetic data in norm for different harmonics number has been calculated for given . It turns out that for cases and it was found to be above . Further increase of the cut-off parameters and almost does not affect that value. Then to generate synthetic data we eliminate these values by .

##### 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:

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):where is a solution of the direct problem for given . Let the function be a solution of the following adjoint problem:subject to boundary and periodicity conditions:Now calculate the first variation of the residual using the variation of the Lagrangian:where is a part of operator (1) acting with respect to variables and . Let us consider second term of expression (43). Integrating by part and taking into account conditions at discontinuity points of coefficient , we obtainHere is the jump of the function at . By conditions (4) these values vanish. ThenUsing boundary conditions (6) and (42) we haveFinally, substituting (46) into expression (43), integrating by part over , and applying conditions for the function we get

Note that functional (38) is represented in the form of a sum of two squared norms in the space :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 . 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

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:

Let us apply the Fourier method to solve system (50). Represent the functions , in the form of series similar to (10):

Here for simplicity the set of basic functions is enumerated anew sequentially:so that each pair of corresponds to the unique index and vice versa.

Let us firstly find functions and 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:

As it is shown in Section 3.1 functions satisfy (11) and relations (23); that is, . Then we obtain

Consider now the functions . They satisfy the same equation as (11), which can be transformed in the form of Bessel equation (13). Now we consider its auxiliary solutions subject to the following boundary conditions:and the conjugation conditions (16) on the internal boundaries . Then coefficients of series (52) for each index are presented as a product of undetermined coefficients and auxiliary solutions:Then first derivative at inner boundary is calculated asLet us compare values and . By direct computations one can see that the variable transforms applied to the functions and , respectively, lead to the same equation and the same boundary conditions for the transformed functions and . Then by unicity .

ThenBecause , this yields to equality . Therefore This yieldsFinally, the approximate regularized solution of the considered inverse problem is obtained in explicit form:where are defined by expression (24). Via standard analysis of the function we get an estimate . Then the convergence of series (63) follows from the estimateIt follows from (64) that the regularized problem is well-posed, because minimizer of functional (38) exists, is unique, and depends continuously on input data .

Since the direct problem with function 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 numbers , 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 likewhere , are random variables uniformly distributed on the interval . As a measure of the relative noise level we will take the value:

##### 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 harmonics , over the angles and and regularization parameter . Generate synthetic and noisy data. Expand the boundary measurement with noise (38) to Fourier series by calculating the coefficients of series (37), ; .(ii)Step 2. Put in (36) value 1.0 instead of . Solve the linear system (35)-(36) and define the coefficients , for each pair of indexes and .(iii)Step 3. Calculate first derivatives of functions at and find values defined by (24): .(iv)Step 4. Calculate the values , .(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 and .

##### 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

*Example 1. *We reconstruct the smooth boundary condition aswith the following parameters: , , , , , , , , . The numerical reconstructions are obtained in this case for various percentages of noise are shown in Figures 1(b), 1(c), 1(d), and 1(e), respectively. relative errors were computed to measure the quality of reconstructed boundary condition as follows:The comparison of relative errors for different noise levels 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 . In particular it follows from table that a large value of the noise level requires a larger value for the parameter of regularization.