Non-Iterative Solution Methods for Cauchy Problems for Laplace and Helmholtz Equation in Annulus Domain

This note is concerned with two new methods for the solution of a Cauchy problem. The first method is based on homotopy-perturbation approach which leads to solving a series of wellposed boundary value problems. No regularization is needed in this method. Laplace and Helmholtz equations are considered in an annular region. It is also proved that the homotopy solution for the Laplace operator converges to the actual exact solution. The second method is also non-iterative. It is based on the application of the Green’s second identity which leads to a moment problem for the unknown boundary condition. Tikhonov regularization is used to obtain a stable and close approximation of the missing boundary condition. A number of examples are used to study the applicability of the methods with the presence of noise.


Introduction
In this note, we introduce two direct (Non-iterative) methods for the solution of a particular ill-posed elliptic inverse problem. Specific applications of this problem in annulus domains appear in various fields including thermal systems [1], corrosion detection in pipes [2,3], interior boundary evaluation in Tokamak [4,5], reconstruction of interior voltage distribution [6], and continuation of magnetic field [7]. This problem falls within the general area of inverse problems (IHCP) that have numerous applications in various fields of engineering [8,9].
The specific problems studied in this note are the Cauchy problems of elliptic systems such as Laplace and Helmholtz [10,11] operators within annulus domains, where no information on the interior boundary is available.
Due to its importance the literature on this problem is vast. Recent results on this particular problem includes meshless methods [12,13], an optimal regularization method [14], a fitting algorithm [15], singular boundary method [16], discrete Fourier transform method [17], a method based on proper solution space [18], and an energy regularization method [19].
Almost all existing results are iterative [20] (and references therein). Iterative algorithms require an appropriate initial guess, and in case of minimization based algorithms [21], some level of convexity. The purpose of this note is to develop two direct (non-iterative) methods for inverse problems for annular domains. Analytical issues for this problem has been addressed in [22]. Iterative methods for this particular application in fusion research has been developed [4,23] (also references therein). However, fast direct (non-iterative) solution methods are more suitable because they can be used on-line for feedback control. The first method is based on Homotopy-perturbation technique [24].
Homotopy perturbation approach leads to solving a series of well-posed boundary value problems. It does not require any regularization and is direct. In case of partial data, it can also be used as an initial guess for optimization-based iterative methods [25]. It is also proved that this method can recover the exact boundary condition for the Laplace operator. This method is presented in Section 2. We also apply this method to the Helmholtz operator in an annular region which is presented in Section 3. The second method is based on the application of Green's identity. We have presented a different approach that uses Green's identity in [26]. Here, we are using point sources to excite the domain, and the formulation is simpler. It also does not require a fictitious domain. This method is also non-iterative and requires regularization. This method is presented in Section 4.
Notation: We use subscripts * x , * r , * z , * θ to denote differentiation with respect to the independent variable. We use integer subscript to denote an element in a finite dimensional space. For clarity, we specifically denote the normal derivative at the boundary by ∂ ∂n .

A Direct Method Based on Homotopy-Perturbation
Consider a Cauchy problem for a ring shown in Figure 1, and assume that the outer boundary, i.e., (r = 1), is accessible and can be used to collect measurements, and is given by Measurements in the form of flux can be collected and provided according to The Cauchy problem of interest is to solve for the temperature u(r, θ), and in particular the inner boundary condition u(w, θ). Assuming a homotopy where the scalar p : 0 → 1, we can consider an elliptic system given by [24] where for p = 1 we recover (1). We can then assume a perturbation for the temperature given by where p is the perturbation parameter, and u k is the temperature at different orders. For various orders it leads to p 0 : r 2 u 0 rr + ru 0 r = 0, p 1 : The boundary conditions can also be imposed according to Analytical solutions can be obtained according to where * ( ) denotes the -th order derivative d dθ . For p = 1, we have the solution given by In a more compact form, it is given by Applying the ratio test to the first and second summations leads to respectively. Therefore, if the functions f (θ) and g(θ) are analytic, and if the ratios g (2n+2) /g (2n) and f (2n+2) / f (2n) are bounded for all n, then the two power series converge for all r with 0 < w ≤ r ≤ 1. These two conditions are very restrictive, because the functions g (2n) , f (2n) can go through zero. Therefore the ratio test can only be used if the functions f (θ) and g(θ) and their derivatives satisfy required bounds. The ratio test is sufficient but not necessary. We next proceed to show that u(r, θ) in Equation (7) is indeed the solution to the elliptic problem posed in Equation (1). The Dirichlet problem for Equation (1) with specified boundary conditions at r = 1 and r = w is well-posed and has a unique solution [27]. We can proceed as follows. We can use Equation (7) to obtain the missing boundary condition for the problem. Using this condition and the given Dirichlet condition in Equation (1), we can write down the exact solution for the given Dirichlet problem ( [27], Section 9.5). Using this solution we can then evaluate the gradient of the temperature at r = 1, and show that, it is indeed the given gradient condition in Equation (2). (1) and (2) is u(r, θ) given in Equation (7).

Theorem 1. (main result) The solution to the Cauchy problem given in Equations
Proof. Consider the Dirichlet problem Equation (1) ( [27], Section 9.5) for the ring shown in Figure 1. The Dirichlet boundary conditions are given as u(1, θ) = f (θ) and u(w, θ) = φ(θ). We can use the solution obtained in Equation (7) to provide the Dirichlet condition at r = w according to The exact solution for this particular Dirichlet problem in given by [27] u(r, θ) [(a n r n + b n r −n ) cos nθ + (c n r n + d n r −n ) sin nθ].
For our specific geometry, the outer radius is r = 1 and the inner radius is r = w, and the coefficients are given by a n + b n = 1 π 2π 0 f (τ) cos nτdτ, w n a n + w −n b n = 1 We next proceed to use the value of u(w, θ) = φ(θ) in Equation (9) in Equation (11) and obtain the exact solution for the temperature u(r, θ) with the Dirichlet conditions u(1, θ) = f (θ) and u(w, θ) = φ(θ). We can then obtain normal derivative u r (1, θ) and show that it is indeed the specified u r (1, θ) = g(θ) in Equation (2). Evaluating the gradient u r (1, θ) leads to It is easy to note that Equation (12) is the Fourier representation of a function. We can proceed to show that it is indeed the Fourier representation of g(θ). Considering Equation (11) and substituting for φ(τ) leads to Due to the periodic boundary conditions, the terms involving 2π 0 g (2n) (τ)dτ and 2π 0 f (2n) (τ)dτ vanish, and the above relation leads to We next consider Equation (14) and substitute for φ(θ) from Equation (9) which leads to Various terms involving 2π 0 g (2n) (τ) cos nτdτ and 2π 0 f (2n) (τ) cos nτdτ can be integrated by parts leading to n ln(w) + (n ln(w)) 3 3! + (n ln(w)) 5 5! + ... g(τ) cos nτdτ.
The terms in the second and third integrals on the right hand side converge to 1 + (n ln(w)) 2 2! + (n ln(w)) 4 4! + (n ln(w)) 6 6! ... = cosh(n ln(w)) = e n ln(w) + e −n ln(w) 2 = w n + w −n 2 n ln(w) + (n ln(w)) 3 3! + (n ln(w)) 5 Replacing the above terms and simplifying leads to Following a similar procedure we can show that It follows that Equation (12) is indeed the Fourier representation of the gradient of the temperature u r (1, θ) = g(θ). This completes the proof.
We next consider a numerical example. Example 1. Consider the elliptic system given in Equation (1) and assume that the external boundary condition is given by u(1, θ) = f (θ) = 0.5 + 0.5e −(θ−π) 2 /1.5 , with the normal gradient at the boundary as the measurement g(θ). The actual boundary condition is given by We first use this give or actual (or the exact) boundary condition and generate the data, i.e., the gradient at the outer boundary. To simulate a realistic data, we add 3% noise, which is a zero mean, Gaussian, white noise generated with a simple random number generator. We can then provide this as the data. Figure 2 shows the recovered internal boundary at w = 0.63 with no noise in the data and, Figure 3 presents the results for the same problem with 3% noise.

Application to a Helmholtz Equation
Consider a similar Cauchy problem for the Helmholtz equation given by where k is the wave number and, the flux at r = 1 is provided. Assuming a similar homotopy where the scalar p : 0 → 1, we can consider an elliptic system given by where (r, θ) ∈ Ω, and for p = 1 we recover (20). Similarly, the solution can be obtained in the form of a regular perturbation given by where various orders of the field variable must satify the following equations The boundary conditions can also be imposed according to The zero order is the homogenous Bessel equation of order zero. The solution is given by where J 0 (rk) and Y 0 (rk) are the Bessel functions of order zero. The above solution is unique because the determinant d 0 (k) = 2 πk is away from zero. For p 1 , we have a nonhomogeous Bessel equation given by where we are also denoting the known term on the right-hand-side by Γ 0 (r, θ). The above equation is subject to the boundary condition given in Equation (24). Using the method of variation of parameters [28] (p. 103) the solution is given by where D(sk) is the Wronskian and is given by for r ∈ Ω [29] (p. 302). Since D(rk) is away from zero for r ∈ Ω, the above solution is well defined. Higher order terms can also be obtained after solving similar nonhomogenous Bessel equations given by We next consider a numerical example for the Helmholtz operator. (20) and assume that the external boundary condition is given by u(1, θ) = f (θ) = 1 + 0.3e (−(θ−π)) 2 /0.7) . The actual boundary condition at r = w is given by u(w, θ) = 5 + e (−(θ−1.5) 4 /.02) + e (−(θ−4.51) 4 /0.02) , (29) and k = 5. Figure 4 shows the recovered internal boundary at w = 0.6. Figure 5 presents the unknown interior boundary for the same boundary condition at r = 1. The actual interior boundary condition is given by

A Direct Method in Terms of a Moment Problem
Consider the same bounded domain Ω shown in Figure 1. The domain is enclosed by smooth boundaries Γ 1 (or r = 1) and, Γ 2 (or r = w) with Γ 1 ∪ Γ 2 = ∂Ω and Γ 1 ∩ Γ 2 = 0. Consider a similar Cauchy problem given by Equation (1), i.e., This is similar to the problem given in (1) and (2), where Γ 2 is the interior boundary where no information is given. Consider two well-posed elliptic problems given by where δ(x − s j ) is a delta function centered at s j ∈ Ω, and σ j (θ) is a smooth function chosen by the designer. Applying the Green's 2-nd identity lead to where ∂ ∂n is the outward normal derivative. Using Equations (31)-(33) leads to The left hand sides are the field variable at the same point s j ∈ Ω. Equating the right hand sides leads to a moment problem given by where, the quantities on the right hand side are either known, or can be computed. By changing the location of the delta function s j ∈ Ω and the external sampling function σ j (θ), one can obtain a collection of moment problems for the unknown u(w, θ) given by The above moment problem can be used to solve for the unknown boundary condition u(w, θ), if the functions ∂γ ∂n j − ∂τ ∂n j are linearly independent. Before proceeding to prove the linearly independence of these functions, it is uselful to note that the solution to the elliptic system in Equation (33) (or (32)) is the combination of the solution to two problems, namely, where γ δ is the solution due to the point source, i.e., Equation (39). Also,γ is the solution due to the system being excited at the boundary only, i.e., Equation (40). Since the location of the point source is the same in both problems in Equations (32) and (33), this portion of the solution can be subtracted out. Also, since the boundary condition at Γ 1 is the same, i.e., f (θ), it is sufficient to state the following theorem. Proof. It is sufficient to consider two elliptic problems given by The solution to the above problems are given by [27] (pp. 340). Assume that σ j (θ) and σ k (θ) computing the normal gradient at r = w, multiplying the latter by µ and subtracting it from the former and simplifying leads to , then the righthand-side must be equal to zero for L =L, and L =L + 1, L =L + 2, L =L + 3,..., wherê L is a sufficiently large integer. This leads to σ j (θ) − µσ k (θ) = 0 ∀ θ ∈ [0, 2π], which is a contradiction. This completes the proof.
The moment problem given in Equation (38) is still ill-posed [30]. Assuming an expansion for the unknown boundary condition u(w, θ) according to u(w, θ) = ∑ L =1 ξ c (θ), Equation (38) leads to a non-square linear system given by where the unknown coefficients ξ are placed in the vector ξ and the known right-hand side entries λ j are placed in the vector λ. The above system can be solved after introducing Tikhonov regularization [31] which puts a bound on the slope of the unknown function. We next consider two specific examples.
Example 3. Consider the elliptic system given in Equation (1) and assume that the external boundary condition is given by f (θ) = 0.5(1 + exp(−(θ − π) 2 /1.5)). The problem is to recover the temperature at the interior boundary which is located at w = 0.7. One can choose 30 locations for the point sources, and 30 sampling functions that are denoted by σ j (θ). Dividing the θ direction into 600 equal intervals leads to dθ = (2π/600), and dividing the radial direction into 60 equal intervals leads to dr = (1 − w)/60. Appropriate locations for the point sources can be given by s j = (r j , θ j ) = ((w + (10.5dr) + (7j/6)), ((20dt) + j(2π/33))), for j = 1, 2, .., 30. This is to avoid the singularity associated with the solution of Equations in (32) and (33). To generate the sampling functions σ j (θ), j = 1, 2, ..., 30 we use a combination of cubic B-splines and sine functions. We essentially need smooth linearly independent functions that satisfy the periodic boundary condition of the problem. We are using Cubic-B splines, and the details are given in [26].
To approximate u(w, θ) at Γ 2 we use cubic B-splines with L = 28. Note that L, or the number of unknown coefficients, i.e., ξ , in the expansion u(w, θ) = ∑ L =1 ξ c (θ) needs to be lower than the number of point sources M = 30, or L < M. Once the location of the point sources and the sampling functions are selected, various terms can be computed and one can proceed to compute the non-square matrix given in (44). Figure 6 shows the normalized singular values of the coefficient matrix Γ. The actual boundary is given by There are about 12 significant (non-zero) eigenvalues out of 28, which is characteristic of a moment problem. In order to be able to invert the matrix, we need to introduce regularization by imposing a bound on the slope of the recovered function. We also need to require that the recovered function be periodic. We can impose that the recovered function be continuous by requiring that Therefore instead of Equation (44), we proceed to solve a least-square problem given by where Φ represents the first-derivative operator, and φ represents the condition in (45). For the numerical examples, we use β = 0.5 × 10 −4 , and α = 1 × 10 6 . To ensure a smooth inversion, we use singular-value-decomposition and ignore the singular values lower than 1 × 10 −4 . Figure 7 compares the recovered boundary to the exact value with 0.04% noise level.    A number factors can somewhat affect the results. In the work presented here, we did not require for the thickness of the annulus to be small. Methods presented here can recover the unknown boundary condition for all thicknesses. As the interior boundary gets closer to the outer boundary the results improve.

Conclusions
In this note we presented two non-iterative (direct) methods for a specific Cauchy problem. The first method is based on Homotopy-perturbation that leads to the solution in the form of a series. We applied this method to Laplace and Helmholtz operators. We also proved that the Homotopy perturbation solution is indeed the unique solusion of the Cauchy prblem for the Laplace operator. The second method is based on the application of the Green's second identity and leads to a moment problem. Both methods show good robustness to noise. The methods can be applied to a number of elliptic systems.
Author Contributions: These two authors contribute equally to this paper. All authors have read and agreed to the published version of the manuscript.