- Email: [email protected]
A numerical method for inverse source problems for Poisson and Helmholtz equations
Physics Letters A 380 (2016) 3707–3716
Contents lists available at ScienceDirect
Physics Letters A www.elsevier.com/locate/pla
A numerical method for inverse source problems for Poisson and Helmholtz equations A. Hamad, M. Tadi ∗ Mechanical Engineering Department, University of Colorado Denver, Denver, CO 80210, United States
a r t i c l e
i n f o
Article history: Received 20 June 2016 Received in revised form 14 August 2016 Accepted 16 August 2016 Available online 21 September 2016 Communicated by C.R. Doering Keywords: Inverse source problem Poisson’s equation Helmholtz equation Separable source function
a b s t r a c t This paper is concerned with an iterative algorithm for inverse evaluation of the source function for two elliptic systems. The algorithm starts with an initial guess for the unknown source function, obtains a background field and, obtains the working equations for the error field. The correction to the assumed value appears as a source term for the error field. It formulates two well-posed problems for the error field which makes it possible to obtain the correction term. The algorithm can also recover the source function with partial data at the boundary. We consider 2-D as well as 3-D domains. The method can be applied to both Poisson and Helmholtz operators. Numerical results indicate that the algorithm can recover close estimates of the unknown source functions based on measurements collected at the boundary. © 2016 Elsevier B.V. All rights reserved.
u (x) + 2 u = S (x) in ,
1. Introduction In this paper we develop a numerical method to recover the source function for an elliptic system in two and three dimensions. This problem falls within the general area of non-invasive imaging techniques that have seen increased amount of utility and applications in various fields of sciences. For example in diffuse optical tomography [1] one seeks to investigate neuronal activity in the living human brain. In bioluminesence imaging, biological entities (e.g. tumor cells) are tagged with luciferase enzymes and implanted in a small animal [2]. In geological applications, the method of self-potential is used where the interest is to recover the source function for a Poisson equation [3]. In electronic devices, the source function consists of piecewise constant functions on rectangular blocks, also called checkered functions [4]. In oil industry, the interest is to recover the pressure field which is often referred to as soil venting [5], to name a few. Consider a bounded domain ∈ R 2 (or R 3 ). In this paper, we consider the inverse problem of identifying a source function, S (x), for an elliptic system given by
u (x) = S (x) in .
(1)
The proposed method can also be applied to an inverse source problem for a Helmholtz operator given by
*
Corresponding author. E-mail address: [email protected] (M. Tadi).
http://dx.doi.org/10.1016/j.physleta.2016.08.057 0375-9601/© 2016 Elsevier B.V. All rights reserved.
where is the wave number. The boundary of the domain, i.e., ∂, is accessible and can be used to collect measurements. It is well known that such problems are highly ill-posed [6]. A number of investigators have developed various methods to deal with the inherent ill-posedness of such problems. Recent results for the Poisson operator includes a method based on boundary element method which uses Trefftz test function [7] and, methods based on the fundamental solution [8] which is a meshless boundary collocation method which belongs to the family of Trefftz methods. These methods use appropriate test functions that in some sense satisfy the differential equation under study. Additional methods include a method based on the Green’s function [9] which seeks to identify the strength and locations of point sources, a direct algorithm [10] which models the source function in terms of collections of point-sources and, a method based on reciprocity functional [11] in which the particular reciprocity functional maps harmonic functions to their integral defined over the support of the unknown source function. More general methods such as optimization based algorithms have also been applied to the inverse source problem [12,13]. A major difficulty with the inverse source problems for Poisson and Helmholtz equations is that the solutions are not unique [14]. It is possible to obtain unique solutions if one can introduce additional requirements on the unknown function. One of the early computational methods that was specifically devised for the inverse source problem for Poisson equation is the minimum support functional [15]. In this method, the algorithm seeks a source
3708
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
function with minimum volume (compactness) to restrict the solution of the inverse problem. In other words, the source is zero everywhere except within an unknown domain (support) and the minimum support functional seeks to minimize the volume enclosed by the function and the base. An a priori estimate of this volume was used to ensure the convergence of the algorithm to a nonzero function. The cost functional was later on modified to minimize the area over which the source function has a nonzero gradient, also known as minimum gradient support functional [16]. Another approach to deal with the non-uniqueness of this problem is to simply restrict the form of the source function. In [17], the authors developed a direct inversion algorithm to identify a one-point-mass-like source using boundary element approach. In [18], it was shown that a source function can be uniquely recovered if some a priori information is give. It was shown that a unique source function can be obtained, 1) if separation of variables are possible and one factor of the product is known; 2) if the domain and the source is in a cylindrical geometry. In another computational effort for the Poisson equation [19] the source function is assumed to be zero except within an unknown domain within which it’s value is equal to one. In other words the unknown source function is a characteristic function given by
χω (x) =
x∈ω
1
if
0
otherwise
,
on a domain ω ⊂⊂ . In fact uniqueness results for source functions with star shaped characteristic functions have been presented in [6]. In another formulation [20], the authors assume that the source function satisfies the homogeneous Laplace equation, i.e., the source is harmonic in . It is often the case that the unknown source function is composed of a uniform background and a small perturbation with compact support ([11], Eq. (1)). This is the case for the source functions that are under study in this paper. Such a source function appears in electromagnetic source imaging [21], and in particular in non-destructive evaluation of materials using lock-in vibrothermography ([22], Eq. (1)). The elliptic operator studied in [22] is the modified Helmholtz equation with a similar geometry that is under study in this paper (figure 1. in [22]). In this paper we consider source functions composed of a known uniform background and an unknown perturbation that is separable. This type of source function is similar to the unknown function in [11], except that the unknown perturbation part is a separable function. The value of the (known) uniform background can be zero, similar to the source functions in a number of references mentioned (for instance, [22,6]). However, here we can have a nonzero value, simply because our method treats the equation for the error field. Since the background field is assumed to be known then, the unknown source (the correction) for the error field has a zero background value. In this paper, we assume that the unknown part of the source function is separable, which is appropriate for some perturbations with small compact support. The present algorithm is iterative in nature and is composed of three steps: I: Assume a value (often the uniform value of the source) for the sought-after source function and use the Dirichlet conditions to obtain a background field. II: Use the collected data and the knowledge of the Dirichlet boundary condition and obtain working equation for the error field. This leads to an error field that is over-specified and therefore is ill-posed. In this paper, we introduce new methods to obtain solutions to the error field. In essence, the method considers two well-posed formulations of the error field. It is then possible to obtain corrections to the assumed values of the source function in step I.
Fig. 1. Error field with over-specified boundary conditions,
III: Update the assumed values and repeat steps I and II. In section 2 we develop the method for a 2-D domain. In section 3 we consider a 3-D domain, and develop a slightly different approach to achieve step II of the above algorithm. In section 4 we use a number of 2-D and 3-D examples to study the applicability of the method. 2. Inverse source problem in 2-D Consider a bounded domain in R 2 , with a continuous boundary ∂. Let S (x) = 1 + g (x) f ( y ) represent the source function in . For an inverse source problem, the goal is to recover the source term for an elliptic system given by
u xx + u y y = 1 + g (x) f ( y )
in ,
(2)
where is a unit square and the variable u (x) can be the electric potential or the material temperature. The boundary is accessible, and it is possible to measure both flux (or Neumann boundary) and Dirichlet conditions on the boundary ∂. Starting with an assumed value for the source term 1 + gˆ (x) ˆf ( y ) and using the Dirichlet condition, it is possible to obtain a background field, uˆ (x), given by
uˆ xx + uˆ y y = 1 + gˆ (x) ˆf ( y )
in .
(3)
Subtracting Eqn. (3) from Eqn. (2) leads to the error field, e (x), given by
e xx + e y y = g (x) f ( y ) − gˆ (x) ˆf ( y )
in
,
(4)
where, e (x) = u (x) − uˆ (x). The actual source functions are related to the assumed values according to
g (x) = gˆ (x) + σ (x),
and,
f ( y ) = ˆf ( y ) + ( y ),
(5)
where now, the correction terms σ (x) and ( y ) are unknowns. It is possible to recover these corrections one term at a time according to the following. This is the step II of the algorithm outlined in section 1. 2.1. Update in y direction First assume that g (x) ≈ gˆ (x) then, the above equation for the error field simplifies to
e xx + e y y = gˆ (x) ( y )
in .
(6)
The background field satisfies the Dirichlet condition, and therefore e = 0 on the boundary ∂. The data is in terms of flux at the boundaries as shown in Fig. 1. In order to recover the correction in the source term ( y ) it is possible to divide the domain in the y direction into ne equal intervals, i.e., y, and use the central differencing in the y direction according to
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
e xx +
1
y2
(e j −1 − 2e j + e j +1 ) = gˆ (x) ( y j ) = gˆ (x) j .
(7)
It is possible to write the above system in a vector form given by
⎡
−2
1
⎢ ⎢ 1 d e(x) 1 ⎢ ⎢ .. + dx2
y2 ⎢ ⎢ ⎣ .. .. 2
⎤
0
0
..
−2 1
0
0 ⎥
.. 0
..
n−2
(8)
T
the unknown source term, respectively. Note that dividing the domain y = [0 : 1] into ne equal elements leads to n = ne + 1 nodes. However, due to the zero boundary conditions for the error field e 1 = en = 0, there are only (n − 2) interior nodes at which the error field is unknown. Also, it is reasonable to assume that the unknown source functions are known at the boundary and as a result it is reasonable to assume that the corrections are also zero at the boundaries. It is due to the zero boundary condition for the error field that the coefficient matrix in Eq. (8) remains symmetric. The data is collected at the boundaries and as a result the above system is required to satisfy the boundary conditions given by
de(x)
e(x)|x=0 = 0,
dx de(x)
e(x)|x=1 = 0,
dx
|x=0 = w ( y j ), |x=1 = e ( y j ).
(9)
Since the matrix in Eq. (8) is symmetric, it is possible to decouple the above system using the matrix of eigenvectors according to e = v where ∈ R (n−2)×(n−2) is the orthogonal matrix of the eigenvectors of the matrix in Eq. (8) with T = I. The above system simplifies to
d2 v i (x) dx2
+ λi v i (x) = gˆ (x)c i
(10)
where v T = [ v 2 (x) v 3 (x), .., v ne (x)] T , λi are the eigenvalues of the matrix in Eq. (8), and c = T d. In addition, the boundary conditions in the modal space are given by
dx
|x=0 =
|x=1 =
T
dx de(x) dx
dx2
(11)
+ λi v i (x) = gˆ (x)c i ,
v i (x)|x=0 = 0,
dv i (x) dx
|x=1 = ˆi P. I,
v i (x)|x=1 = 0,
P. II,
(13)
⎢ 12 ⎢ x ⎣ ... ...
−2
x2
0
+ λi ... 1
2 x
1
x2
...
−4 2 x
A1
⎡
0
⎡
⎤
⎤ ⎡ v1 ⎤ ... i ⎢ 2⎥ ... ⎥ vi ⎥ ⎥⎢ ⎢ ⎥ ... ⎦⎣ ... ⎦
...
3 2 x
vn
i w
⎤
0 ⎢ gˆ (x2 ) ⎥
⎢0⎥ ⎢ ⎥ ⎥ : ⎥ c i . P. I, =⎢ ⎣ : ⎦+⎢ ⎣ gˆ (x ) ⎦ ne ˆi 0 b1
(14)
gˆ
where, for simplicity, we are also defining new variables according to w T = [ v 1i , v 2i , ..., v ni ] T , b1T = [0, 0, ..., ˆi ] T , and gˆ T = [0, gˆ (x2 ), .., gˆ (xne ), 0] T , and similarly, A1 . Note that the last row (equation) enforces the slope condition at x = 1. Similarly, a finite difference approximation to the second problem is given by
⎡
−3 2 x 1
x2
⎢ ⎢ ⎢ ⎢ ⎣ ... ...
x2
...
... 0
⎡
i
⎤
1
vn
i w
A2
⎡
0
1
...
⎤⎡
⎤ v1 ⎥⎢ i ⎥ ... ⎥ v 2i ⎥ ⎥⎢ ⎥ ⎥⎢ ⎣ ... ⎦ ... ⎦
−1
x
4 2 x −2 + λi
x2
0
⎤
gˆ (x2 ) ⎥ ⎢0⎥ ⎢ ⎥ ⎥+⎢ =⎢ c i P. II. ⎢ ⎣ : ⎦ ⎣ : ⎥ ⎦ gˆ (x )
(15)
ne
0
b2
where, we are also defining new variables and ˆ for simplicity. The above system is now fully decoupled. Every equation in the system is a second order differential equation that is required to satisfy two boundary conditions at each end. Each equation also includes an unknown constant c i . Therefore we are faced with a system of ill-posed linear boundary value problem. After solving the above equations and obtaining the values of unknown constants c i , the correction to the assumed value of the source term can be computed using the transformation d = c. It is possible to apply the method used in [23] to obtain solutions to the above boundary value problems. The approach is to consider two well-posed boundary value problems and note that they describe the same solution. In other words, instead of Eqs. (10) and (11), we consider two well-posed problems given by
d2 v i (x)
1
|x=0 = , ˆ |x=1 = ,
|x=0 = i
0
v(x)|x=0 = v(x)|x=1 = 0, dv(x) T de(x) dx dv(x)
+ λi v i (x) = gˆ (x)c i ,
for i = 1, 2, ..., (n − 2). Note that there are (n − 2) decouple equations. Both of the above problems are well-posed, and one can proceed to obtain solutions using finite difference approximations. If one divides the x domain into ne equal intervals and use central differencing to approximate the second derivative they lead to
⎡
where, e(x) = [e 2 (x) e 3 (x), ., ene (x) ] ∈ R and d = [ ( y 2 ) ( y 3 ), . ( yne )] T ∈ R n−2 are the unknown error field and T
dx2 dv i (x) dx
⎥ ⎥ .. .. .. ⎥ ⎥ e = gˆ (x)d, ⎥ 1 −2 1 ⎦ 0 1 −2
T
d2 v i (x)
3709
(12)
where, for simplicity, we are again defining new variables according to b2T = [i , 0, ..., 0] T , and, similarly, A2 . The above two problems are well-posed and the coefficient matrices can be stably inverted. Eliminating the unknown vector, w, from the above two equations leads to 1 −1 −1 −1 ˆ ˆ A− 1 b1 + A1 gc i = A2 b2 + A2 gc i ,
−1
−1
−1
−1
(A1 b1 − A2 b2 ) = (A2 − A1 )ˆg c i ,
r
or, (16)
s
where, we are also defining new vectors r and s for simplicity and space considerations. It is now possible to solve for the unknown scalar according to
ci = β
sT r sT s
,
(17)
where, the positive scalar β is used 0 < β and is set by the designer. We will explain this parameter in more details. Once the scalars c i , i = 1, 2, ..., (n − 2) are computed, then the correction to the assumed value of the source term can be solved for using the transformation d = c. We can apply a similar procedure to update the assumed value of the source term in the x direction.
3710
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
d2 v i (x)
+ λi v i (x) = gˆ (x)c i , dx2 v i (x)|x=0 = 0, v i (x)|x=1 = 0 d2 v i (x) dx2 dv i (x) dx
P. I,
(22)
+ λi v i (x) = gˆ (x)c i ,
|x=0 = ki
v i (x)|x=1 = 0,
P. II.
(23)
Both of the above problems are well-posed and one can follow the same formulation and obtain the unknown constant c i . The same procedure can be used to recover the correction in the x direction. We next proceed to consider inverse source problem in 3-D.
Fig. 2. Error field with over-specified boundary conditions.
3. Inverse source problem in 3-D Consider a bounded domain in R 3 , with a continuous boundary ∂. Let S (x) = 1 + g (x) f ( y )h( z) represent the source function in with the non-uniform portion being separable. The goal is to recover the source term for an elliptic system given by
u xx + u y y + u zz = 1 + g (x) f ( y )h( z)
2.2. Update in x directions After updating the source term in the y direction, it is possible to assume that f ( y ) = ˆf ( y ). Then equation (4) for the error field simplifies to
in .
(18)
The background field satisfies the Dirichlet condition, and therefore e = 0 on the boundary ∂. The data is in terms of flux at the boundaries as shown in Fig. 2. The same procedure can be used to recover the correction in the x direction. 2.3. Recovery with partial data It is possible to apply the present method when the data is collected at only two sides of the domain. Consider the same inverse problem presented in eq. (2) and assume that the data is collected at only two sides of the domain. Following similar steps and obtaining the error field leads to
e xx + e y y = g (x) f ( y ) − gˆ (x) ˆf ( y ),
in ,
e = 0 on
∂, (19)
uˆ xx + uˆ y y + uˆ zz = 1 + gˆ (x) ˆf ( y )hˆ ( z) in
(25)
.
Subtracting Eqn. (25) from (24) leads to the error field given by
e xx + e y y + e zz = g (x) f ( y )h( z) − gˆ (x) ˆf ( y )hˆ ( z) in
,
(26)
for which we have both Dirichlet and Neumann boundary conditions. The actual source functions are related to the assumed values according to
g (x) = gˆ (x) + σ (x),
and,
f ( y ) = ˆf ( y ) + ( y ),
h( z) = hˆ ( z) + η( z).
(27)
Similar to the problem in 2-D, it is possible to recover the corrections to the assumed values one term at a time. We can proceed as follows. 3.1. Update in z direction We can assume that g (x) ≈ gˆ (x) and f ( y ) ≈ ˆf ( y ). Then, the equation for the error simplifies to
e xx + e y y + e zz = gˆ (x) ˆf ( y )η( z)
in
.
(28)
Since e (x, y , z) = 0 at the boundaries x ∈ ∂, it is possible to assume an expansion of the form given by
in addition to
e x (0, y ) = κ ( y ),
(24)
,
where now is a unit cube. The boundary is accessible and it is possible to measure both Neumann and Dirichlet boundary conditions on the boundary ∂. Starting with an assumed value for the source term gˆ (x) ˆf ( y )hˆ ( z) and using the Dirichlet condition, it is possible to obtain a background field given by
Fig. 3. Error field with partial data at the boundary.
e xx + e y y = σ (x) ˆf ( y )
in
e y (x, 0) = ω(x).
(20)
For updating in y direction, one can assume that g (x) ≈ gˆ (x) and arrive at the system given in Fig. 3. After discretizing the second derivative and using the eigenvectors to decouple the working equations one can arrive at similar equations given by (10) with boundary conditions given by
e (x, y , z) =
∞
A (x, y )sin(π z),
=1
1 A (x, y ) = 2
e (x, y , z)sin(π z)dz,
(29)
0
v(x)|x=0 = v(x)|x=1 = 0,
dv(x) dx
|x=0 = T
de(x) dx
|x=0 = k.
(21)
Similar to eqs. (12) and (13), it is possible to formulate two wellposed problems for the error field with partial data according to
where the superscript simply denotes the component in the function expansion. By selecting sin(·) functions the Dirichlet boundary conditions in the z direction are satisfied. Substituting the above expansion in the error equation leads to
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
3711
Fig. 4. Two well-posed problems for the components of the error field.
∞
A xx + A y y − (π )2 A sin(π z) = gˆ (x) ˆf ( y )η( z)
in
.
=1
1 A y (x, 1) = ζ (x) = 2
(30) Using the orthogonality property of the sine functions, one can arrive at the following equation for each component
A xx + A y y − (π ) A = 2 gˆ (x) ˆf ( y )
2
1
η(z)sin(π z)dz in , 0
(31) for = 1, 2, ..∞. Retaining the first L components leads to
A xx + A y y − (π )2 A = gˆ (x) ˆf ( y )θ ,
1 θ = 2
η(z)sin(π z)dz, in ,
(32)
0
for = 1, 2, .., L. The collected data can provide Neumann boundary conditions for the above equations which are often called modified Helmholtz equation. The above equations are still ill-posed in the sense that every A (x, y ) is required to satisfy both Dirichlet (zero) and Neumann conditions. In addition, every equation includes an unknown scalar θ . Once these scalars are computed, then the correction to the assumed value of the source function in the z direction can be computed according to
η( z) =
∞
θ sin(π z) ≈
=1
L
θ sin(π z).
(33)
=1
In order to solve the above modified Helmholtz equations, one can follow the approach that was developed for the 2-D problem in section 2.1 and consider two well-posed problems given by where = ( ∂∂x2 + ∂ ∂y2 − (π )2 ) is the modified Helmholtz operator. The collected data appears as boundary conditions according to
e x (0, y , z)sin(π z)dz,
(34)
0
A x (1, y ) = ξ( y ) = 2
e x (1, y , z)sin(π z)dz,
(35)
e y (x, 0, z)sin(π z)dz,
(36)
0
1 A y (x, 0) = ν (x) = 2 0
The above problems are both well-posed, and one can obtain finite dimensional approximations using central differencing. Similar to Eqs. (14) and (15) finite difference approximations of the above problems lead to
A3 a = b3 + dˆ θ ,
(38)
A4 a = b4 + dˆ θ ,
(39)
where, A3 and A4 are the finite difference approximation of the Helmholtz operators in problem I and problem II shown in Fig. 4, respectively. The vector a is the nodal values of the unknown variable A (x, y ), b3 and b4 are the known right-hand sides of the finite difference approximations to problems I and II shown in Fig. 4, and the vector dˆ is the nodal values of the known product gˆ (x) ˆf ( y ). Note that the collected data, which appear as Neumann boundary conditions, show up in the vectors b3 and b4 . Also, the assumed values for the source function gˆ (x) ˆf ( y ) appear in the ˆ Both problems are well-posed, as a result A−1 and A−1 vector d. 3 4 exist. It is possible to eliminate e from the above equations and arrive at 1 1 (A−1 b3 − A− b4 ) = (A− − A−1 )dˆ θ , 3
4 4 3 w
(40)
v
where for simplicity and space considerations, we are redefining the vectors w and v. Solving for the scalar θ leads to
θ = β
wT v
(41)
vT v
where, the scalar β > 0 is chosen by the designer. After computing the scalars θ , for = 1, 2.., L, the correction to the assumed value can be obtained using Eq. (33). The updating in other directions can also be accomplished using a similar approach.
For updating in the y direction, we can assume that g (x) ≈ gˆ (x) and h( z) ≈ hˆ ( z). Then, the equation for the error simplifies to
e xx + e y y + e zz = gˆ (x) ( y )hˆ ( z)
1
(37)
3.2. Update in y and x directions
1 A x (0, y ) = μ( y ) = 2
e y (x, 1, z)sin(π z)dz. 0
in .
(42)
The above equation is similar to Eq. (28), and the same approach can be applied. Likewise, for updating in the x direction, we can assume that f ( y ) ≈ ˆf ( y ) and h( z) ≈ hˆ ( z). Then, the equation for the error simplifies to
e xx + e y y + e zz = σ (x) ˆf ( y )hˆ ( z) in
.
(43)
3712
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
Again, the above equation is similar to Eq. (28), and the same approach can be applied. 4. Helmholtz inverse source problem The same procedure described in section 3 can be applied to the inverse source problem for the Helmholtz equation. Starting with an initial guess and following the algorithm leads to the equation for the error field given by
e xx + e y y + e zz + 2 e = g (x) f ( y )h( z) − gˆ (x) ˆf ( y )hˆ ( z)
in . (44) Fig. 5. Reduction in the error for the Example 1.
Relating the actual values to the assumed values according to Eq. (27) and proceeding to recover the correction terms for one direction at a time lead to
e xx + e y y + e zz + 2 e = gˆ (x) ˆf ( y )η( z)
in .
(45)
Following the same procedure, the over-specified modified Helmholtz equation given in eq. (32) simplifies to
A xx + A y y + ( 2 − (π )2 ) A = gˆ (x) ˆf ( y )θ ,
1 θ = 2
η(z)sin(π z)dz,
(46)
0
for = 1, 2, .., L. The above equation is similar to Eq. (32) except for the nonzero wave number . For low to moderate wave numbers the above equations can be solved very accurately. The same procedure can be applied to solve the above ill-posed problem. We next apply the above algorithms to a number of numerical examples.
Fig. 6. Convergence of the recovered source function in the x direction for Example 1.
5. Numerical experiments In this section we use a number of numerical examples to study the applicability of the method. Example 1. First consider an inverse source problem in two dimensional domain. Assume that the actual source term in Eq. (2) is given by
g (x) f ( y ) = 1 + exp
(x − .3)4 0.000001
1 + 2exp
( y − .65)4 0.00001
. (47)
We can divide the domain in each direction into 100 equal intervals. The data is collected in the form of Neumann condition at the boundaries. The data is also contaminated with noise. The noise is a zero mean Gaussian white noise with signal-noise ratio of 6%. Before using the data a three point averaging is performed to smooth out the data. The actual source function has a nominal which is equal to 1. Since its value can be measured at the boundary, it is reasonable to start the iterations from gˆ (x) = ˆf ( y ) = 1. The algorithm seeks to recover the correction to the assumed values first in the y direction, and then in x direction. It is instructive to consider Eq. (6) which is used to update the source term in the y direction. It is presented here again
e xx + e y y =
gˆ (x)
( y ) in .
previous iteration
A similar equation, i.e., Eq. (18), is used to update the source term in the x direction. Considering the above equation it is clear that, in order to be able to recover the unknown function ( y ), one
Fig. 7. Convergence of the recovered source function in the y direction for Example 1.
needs to have the function gˆ (x) away from zero. For this example, the nominal value of the field is equal to one, and the inversion can be accomplished. As a result for this example, the parameter β in Eq. (17) is set equal to one, i.e., β = 1. If the nominal value of the field is close to zero, which is the case in most applications, then the above equation tends to over-estimate the unknown function. As a result a smaller value of β should be used. This will become clear in the next example. Fig. 5 presents the reduction in the error for this example. The error is the difference between the measured data and the computed value at the boundaries given by
Error =
∇n u − ∇n uˆ
2
dx.
∂
Figs. 6 and 7 show the recovered functions after a few iterations and compare them to the actual source functions. For this example, the algorithm can recover the source function very closely after only a few iterations. We next assume that only partial data is available for inverse evaluation of the source function.
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
3713
Fig. 11. Reduction in the error for the Example 2.
Fig. 8. Reduction in the error for the Example 1 with partial data.
Fig. 12. Actual source function for the Example 2. The background field is ≈ 0.0025.
Fig. 9. Convergence of the recovered source function in the x direction for Example 1 with partial data.
Fig. 13. Recovered source function for the Example 2 after 140 iterations.
Fig. 10. Convergence of the recovered source function in the y direction for Example 1 with partial data.
We consider the same problem and assume that data is available on only two sides given by (20). Fig. 8 shows the reduction in the total error. Only the error given in (20) is used in the recovery. With this limited data, it is still possible to formulate two well-posed problems given in eqs. (22) and (23). Figs. 9 and 10 show convergence of the source function when only partial data is available. Example 2. We next consider the recovering of a source function with a nominal value of the field close to zero. In this example, the actual source term is given by
(x − .3)4 0.000001 ( y − .65)4 . × 0.05 + 2exp 0.00001
g (x) f ( y ) = 0.05 + exp
(48)
For this example, since the nominal value of the field is very close to zero (0.05, and 0.0025 in 2-D) the algorithm over-estimates the
sought after function. It is then necessary to use a smaller value for the parameter β . For this case we use a value of β = 0.02. Using a smaller value for β will simply require more iterations. Fig. 11 shows the reduction in error for this case. The same amount of noise is used to model a realistic data. Figs. 12 and 13 show the actual source function and the recovered source function after 140 iterations. The present algorithm can recover a very close estimate of the source function. We next proceed to increase the noise in the data. Fig. 14 shows the reduction in error when the level of noise is increase by a factor of five. Fig. 15 shows the recovered source function after 120 iterations. The present algorithm can still recover a reasonably close estimate of the source function. Example 3. We next consider recovering of a source function in 3-D. Consider recovering a source function that is given by
g (x) f ( y )h( z)
(x − .3)4 ( y − .75)4 .1 + exp = .1 + exp 0.00002 0.00002 4 ( z − .5) . × .1 + exp 0.00002
(49)
For this case, the background field is 1.0 and the nominal (or apparent) value of the function is 0.001, and there is one target that is located at (x, y , z) ≈ (0.3, 0.75, 0.5). We can divide the domain into 50 equal intervals in each direction. This size mesh is required to calculate the field variables in 3-D with high accuracy.
3714
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
Fig. 14. Reduction in the error for the Example 2 with higher noise level. Fig. 18. Comparison of the actual function and the recovered source function in x for the Example 3.
Fig. 15. Recovered source function for the Example 2 after 120 iterations.
Fig. 19. Comparison of the actual function and the recovered source function in y for the Example 3.
In order to project the error down to the space of sine functions in Eq. (29), we use 40 terms in the series. However, when recovering the function back using Eq. (33) we use the first L = 18 terms. This is due to the fact that the noise in the data pollutes the higher frequency components and including them leads to instability. The parameter β = 30. Fig. 16 shows the location of the target within the 3-D domain. Fig. 17 presents the location of the recovered function. Figs. 18 through 20 present the individual components of the source function and compare them to the actual function. For this case, 60 iterations are used to recover the target presented in Fig. 16. Fig. 16. Actual target within a 3-D domain for the Example 3.
Example 4. We next consider recovering of multiple targets in a 3-D domain. Consider the recovering of a source function given by S (x) = 1 + g (x) f ( y )h( z) where,
(x − .3)4 (x − .65)4 + exp g (x) = .1 + exp 0.00002 0.000015 4 ( y − .75) ( y − .65)4 + exp f ( y ) = .1 + exp 0.00001 0.000015 4 ( z − .4) ( z − .7)4 + exp . h( z) = .1 + exp 0.00001 0.000015
(50) (51) (52)
The above source function represents four targets within the domain as is shown in Fig. 21. Fig. 22 shows the recovered targets after 60 iterations. All of the parameters are set to be the same as in Example 3. Also, the same amount of noise is used. The algorithm can recover multiple targets within the source function.
Fig. 17. Recovered target within a 3-D domain for the Example 3 after 60 iterations.
Example 5. We next consider recovering of multiple targets in a 3-D domain for the Helmholtz operator with the wave number
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
Fig. 20. Comparison of the actual function and the recovered source function in y for the Example 3.
Fig. 21. Actual source function with 4 targets for the Example 4.
3715
Fig. 23. Recovered source function after 80 iterations for the Helmholtz operator in Example 5.
A number of 2-D and 3-D examples were used to study the applicability of the proposed algorithm. We considered both Poisson and Helmholtz operators. The proposed method is particularly suitable for identifying a single anomaly with compact support. For the cases that the actual perturbation cannot be modelled by a separable function, the proposed method will converge to the closest approximation of that anomaly that can be described by a separable function. This is not a major weakness, simply because the uniqueness issues that are inherent for both of the above problems limits the actual source functions that can be recovered. The proposed method considers a unit square or cube. For other geometries the method can be modified in a number of ways. One possible approach is to use the method of immersed boundary to convert an irregular geometry into a cube. Another approach that can be used is the exact mapping of a general domain onto a square/cube. For this method the working equations in the transformed space will include a number of terms which will not be amendable to diagonalization. We envision that for these cases, it will be possible to set up an internal iteration to handle such terms. For geometries that are close to cube/square variable mesh finite difference method can be used. Again, an internal iteration can be used to handle the terms that are not amendable to diagonalization. All of these issues will be addresses in future works. References
Fig. 22. Recovered source function after 60 iterations for the Example 4.
2 = 4.5. The unknown source function is the same as in Example 4. Fig. 23 presents the recovered function after 80 iterations. All the parameters and the level of noise are set equal to their values in Example 4. 6. Conclusion In this paper we presented a computational algorithm for the inverse evaluation of the source function for Poisson and Helmholtz equations. The source term is assumed to be composed of a uniform background value and a separable function. The algorithm is iterative in nature. It updates the source function for each coordinate separately. The algorithm shows reasonable robustness to noise. It can also recover the source function with partial data.
[1] K. Ren, Recent developments in numerical techniques for transport-based medical imaging methods, Commun. Comput. Phys. 8 (1) (2010) 1–50. [2] S. Li, W. Driessen, S. Sullivan, H. Jiang, Bioluminescene tomography based on phantoms with different concentrations of bioluminescent cancer cells, J. Opt. A, Pure Appl. Opt. 8 (2006) 743–746. [3] B.J. Minsley, J. Sogade, F.D. Morgan, Three-dimensional source inversion of selfpotential data, J. Geophys. Res. 112 (2007) B02202. [4] A. Artemev, L. Parnovski, I. Polterovich, Inverse electrostatic and elasticity problems for checkered distributions, Inverse Probl. 29 (2013) 075010, 16 pp. [5] M. Slodicka, H.D. Schepper, On an inverse problem of pressure recovery arising from soil venting facilities, Appl. Math. Comput. 129 (2002) 469–480. [6] V. Isakov, Inverse Problems for Partial Differential Equations, Springer-Verlag, Berlin, 1998. [7] C.-S. Liu, A BIEM using the Trefftz test functions for solving the inverse Cauchy and source recovery problems, Eng. Anal. Bound. Elem. 62 (2016) 177–185. [8] A. Karageorghis, D. Lesnic, L. Marin, A survey of applications of the MFS to inverse problems, Inverse Probl. Sci. Eng. 19 (3) (2011) 309–336. [9] Y. Hon, M. Li, Y. Melnikov, Inverse source identification by Green’s function, Eng. Anal. Bound. Elem. 34 (2010) 352–358. [10] A.E. Badia, An inverse source problem in an anisotropic medium by boundary measurements, Inverse Probl. 21 (2005) 1487–1506. [11] N. Roberty, C. Alves, On the identification of star-shape sources from boundary measurements using a reciprocity functional, Inverse Probl. Sci. Eng. 17 (2) (2009) 187–202. [12] S.S. Adavani, G. Biros, Fast algorithms for source identification problems with elliptic PDE constraints, SIAM J. Imaging Sci. 3 (4) (2010) 791–808.
3716
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
[13] A. Canelas, A. Laurian, A. Novotny, A new reconstruction method for the inverse potential problem, J. Comput. Phys. 268 (2014) 417–431. [14] G. Wang, Y. Li, M. Jiang, Uniqueness theorems in bioluminescence tomography, Med. Phys. 31 (8) (2004) 2289–2299. [15] B.J. Last, K. Kubik, Compact gravity inversion, Geophysics 48 (6) (1983) 713–721. [16] O. Portniaguine, M.S. Zhdanov, Focusing geophysical inversion images, Geophysics 64 (3) (1999) 874–887. [17] T. Ohe, K. Ohmaka, Boundary element approach for an inverse source problem of the Poisson equation with a one-point-mass-like source, Appl. Math. Model. 18 (1994) 216–223. [18] A.E. Badia, T. Doung, An inverse source problem in potential analysis, Inverse Probl. 16 (2000) 651–663.
[19] N. Martins, An iterative shape reconstruction of source functions in a potential problem using the MFS, Inverse Probl. Sci. Eng. 20 (8) (2012) 1175–1193. [20] B. Jin, L. Marin, The method of fundamental solutions for inverse source problems associated with the steady-state heat conduction, Inverse Probl. Sci. Eng. 69 (2007) 1570–1589. [21] N.G. Gencer, I. Tanzer, Forward problem solution of electromagnetic source imaging using a new BEM formulation with high-order elements, Phys. Med. Biol. 44 (1999) 2275–2287. [22] R. Celorrio, A. Mendioroz, A. Salazar, Characterization of vertical buried defects using lock-in vibrothermography: II. Inverse problem, Meas. Sci. Technol. 24 (2013) 065602, 9 pp. [23] M. Tadi, A.K. Nandakumaran, S. Sritharan, An inverse problem for Helmholtz equation, Inverse Probl. Sci. Eng. 19 (6) (2011) 839–854.
Contents lists available at ScienceDirect
Physics Letters A www.elsevier.com/locate/pla
A numerical method for inverse source problems for Poisson and Helmholtz equations A. Hamad, M. Tadi ∗ Mechanical Engineering Department, University of Colorado Denver, Denver, CO 80210, United States
a r t i c l e
i n f o
Article history: Received 20 June 2016 Received in revised form 14 August 2016 Accepted 16 August 2016 Available online 21 September 2016 Communicated by C.R. Doering Keywords: Inverse source problem Poisson’s equation Helmholtz equation Separable source function
a b s t r a c t This paper is concerned with an iterative algorithm for inverse evaluation of the source function for two elliptic systems. The algorithm starts with an initial guess for the unknown source function, obtains a background field and, obtains the working equations for the error field. The correction to the assumed value appears as a source term for the error field. It formulates two well-posed problems for the error field which makes it possible to obtain the correction term. The algorithm can also recover the source function with partial data at the boundary. We consider 2-D as well as 3-D domains. The method can be applied to both Poisson and Helmholtz operators. Numerical results indicate that the algorithm can recover close estimates of the unknown source functions based on measurements collected at the boundary. © 2016 Elsevier B.V. All rights reserved.
u (x) + 2 u = S (x) in ,
1. Introduction In this paper we develop a numerical method to recover the source function for an elliptic system in two and three dimensions. This problem falls within the general area of non-invasive imaging techniques that have seen increased amount of utility and applications in various fields of sciences. For example in diffuse optical tomography [1] one seeks to investigate neuronal activity in the living human brain. In bioluminesence imaging, biological entities (e.g. tumor cells) are tagged with luciferase enzymes and implanted in a small animal [2]. In geological applications, the method of self-potential is used where the interest is to recover the source function for a Poisson equation [3]. In electronic devices, the source function consists of piecewise constant functions on rectangular blocks, also called checkered functions [4]. In oil industry, the interest is to recover the pressure field which is often referred to as soil venting [5], to name a few. Consider a bounded domain ∈ R 2 (or R 3 ). In this paper, we consider the inverse problem of identifying a source function, S (x), for an elliptic system given by
u (x) = S (x) in .
(1)
The proposed method can also be applied to an inverse source problem for a Helmholtz operator given by
*
Corresponding author. E-mail address: [email protected] (M. Tadi).
http://dx.doi.org/10.1016/j.physleta.2016.08.057 0375-9601/© 2016 Elsevier B.V. All rights reserved.
where is the wave number. The boundary of the domain, i.e., ∂, is accessible and can be used to collect measurements. It is well known that such problems are highly ill-posed [6]. A number of investigators have developed various methods to deal with the inherent ill-posedness of such problems. Recent results for the Poisson operator includes a method based on boundary element method which uses Trefftz test function [7] and, methods based on the fundamental solution [8] which is a meshless boundary collocation method which belongs to the family of Trefftz methods. These methods use appropriate test functions that in some sense satisfy the differential equation under study. Additional methods include a method based on the Green’s function [9] which seeks to identify the strength and locations of point sources, a direct algorithm [10] which models the source function in terms of collections of point-sources and, a method based on reciprocity functional [11] in which the particular reciprocity functional maps harmonic functions to their integral defined over the support of the unknown source function. More general methods such as optimization based algorithms have also been applied to the inverse source problem [12,13]. A major difficulty with the inverse source problems for Poisson and Helmholtz equations is that the solutions are not unique [14]. It is possible to obtain unique solutions if one can introduce additional requirements on the unknown function. One of the early computational methods that was specifically devised for the inverse source problem for Poisson equation is the minimum support functional [15]. In this method, the algorithm seeks a source
3708
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
function with minimum volume (compactness) to restrict the solution of the inverse problem. In other words, the source is zero everywhere except within an unknown domain (support) and the minimum support functional seeks to minimize the volume enclosed by the function and the base. An a priori estimate of this volume was used to ensure the convergence of the algorithm to a nonzero function. The cost functional was later on modified to minimize the area over which the source function has a nonzero gradient, also known as minimum gradient support functional [16]. Another approach to deal with the non-uniqueness of this problem is to simply restrict the form of the source function. In [17], the authors developed a direct inversion algorithm to identify a one-point-mass-like source using boundary element approach. In [18], it was shown that a source function can be uniquely recovered if some a priori information is give. It was shown that a unique source function can be obtained, 1) if separation of variables are possible and one factor of the product is known; 2) if the domain and the source is in a cylindrical geometry. In another computational effort for the Poisson equation [19] the source function is assumed to be zero except within an unknown domain within which it’s value is equal to one. In other words the unknown source function is a characteristic function given by
χω (x) =
x∈ω
1
if
0
otherwise
,
on a domain ω ⊂⊂ . In fact uniqueness results for source functions with star shaped characteristic functions have been presented in [6]. In another formulation [20], the authors assume that the source function satisfies the homogeneous Laplace equation, i.e., the source is harmonic in . It is often the case that the unknown source function is composed of a uniform background and a small perturbation with compact support ([11], Eq. (1)). This is the case for the source functions that are under study in this paper. Such a source function appears in electromagnetic source imaging [21], and in particular in non-destructive evaluation of materials using lock-in vibrothermography ([22], Eq. (1)). The elliptic operator studied in [22] is the modified Helmholtz equation with a similar geometry that is under study in this paper (figure 1. in [22]). In this paper we consider source functions composed of a known uniform background and an unknown perturbation that is separable. This type of source function is similar to the unknown function in [11], except that the unknown perturbation part is a separable function. The value of the (known) uniform background can be zero, similar to the source functions in a number of references mentioned (for instance, [22,6]). However, here we can have a nonzero value, simply because our method treats the equation for the error field. Since the background field is assumed to be known then, the unknown source (the correction) for the error field has a zero background value. In this paper, we assume that the unknown part of the source function is separable, which is appropriate for some perturbations with small compact support. The present algorithm is iterative in nature and is composed of three steps: I: Assume a value (often the uniform value of the source) for the sought-after source function and use the Dirichlet conditions to obtain a background field. II: Use the collected data and the knowledge of the Dirichlet boundary condition and obtain working equation for the error field. This leads to an error field that is over-specified and therefore is ill-posed. In this paper, we introduce new methods to obtain solutions to the error field. In essence, the method considers two well-posed formulations of the error field. It is then possible to obtain corrections to the assumed values of the source function in step I.
Fig. 1. Error field with over-specified boundary conditions,
III: Update the assumed values and repeat steps I and II. In section 2 we develop the method for a 2-D domain. In section 3 we consider a 3-D domain, and develop a slightly different approach to achieve step II of the above algorithm. In section 4 we use a number of 2-D and 3-D examples to study the applicability of the method. 2. Inverse source problem in 2-D Consider a bounded domain in R 2 , with a continuous boundary ∂. Let S (x) = 1 + g (x) f ( y ) represent the source function in . For an inverse source problem, the goal is to recover the source term for an elliptic system given by
u xx + u y y = 1 + g (x) f ( y )
in ,
(2)
where is a unit square and the variable u (x) can be the electric potential or the material temperature. The boundary is accessible, and it is possible to measure both flux (or Neumann boundary) and Dirichlet conditions on the boundary ∂. Starting with an assumed value for the source term 1 + gˆ (x) ˆf ( y ) and using the Dirichlet condition, it is possible to obtain a background field, uˆ (x), given by
uˆ xx + uˆ y y = 1 + gˆ (x) ˆf ( y )
in .
(3)
Subtracting Eqn. (3) from Eqn. (2) leads to the error field, e (x), given by
e xx + e y y = g (x) f ( y ) − gˆ (x) ˆf ( y )
in
,
(4)
where, e (x) = u (x) − uˆ (x). The actual source functions are related to the assumed values according to
g (x) = gˆ (x) + σ (x),
and,
f ( y ) = ˆf ( y ) + ( y ),
(5)
where now, the correction terms σ (x) and ( y ) are unknowns. It is possible to recover these corrections one term at a time according to the following. This is the step II of the algorithm outlined in section 1. 2.1. Update in y direction First assume that g (x) ≈ gˆ (x) then, the above equation for the error field simplifies to
e xx + e y y = gˆ (x) ( y )
in .
(6)
The background field satisfies the Dirichlet condition, and therefore e = 0 on the boundary ∂. The data is in terms of flux at the boundaries as shown in Fig. 1. In order to recover the correction in the source term ( y ) it is possible to divide the domain in the y direction into ne equal intervals, i.e., y, and use the central differencing in the y direction according to
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
e xx +
1
y2
(e j −1 − 2e j + e j +1 ) = gˆ (x) ( y j ) = gˆ (x) j .
(7)
It is possible to write the above system in a vector form given by
⎡
−2
1
⎢ ⎢ 1 d e(x) 1 ⎢ ⎢ .. + dx2
y2 ⎢ ⎢ ⎣ .. .. 2
⎤
0
0
..
−2 1
0
0 ⎥
.. 0
..
n−2
(8)
T
the unknown source term, respectively. Note that dividing the domain y = [0 : 1] into ne equal elements leads to n = ne + 1 nodes. However, due to the zero boundary conditions for the error field e 1 = en = 0, there are only (n − 2) interior nodes at which the error field is unknown. Also, it is reasonable to assume that the unknown source functions are known at the boundary and as a result it is reasonable to assume that the corrections are also zero at the boundaries. It is due to the zero boundary condition for the error field that the coefficient matrix in Eq. (8) remains symmetric. The data is collected at the boundaries and as a result the above system is required to satisfy the boundary conditions given by
de(x)
e(x)|x=0 = 0,
dx de(x)
e(x)|x=1 = 0,
dx
|x=0 = w ( y j ), |x=1 = e ( y j ).
(9)
Since the matrix in Eq. (8) is symmetric, it is possible to decouple the above system using the matrix of eigenvectors according to e = v where ∈ R (n−2)×(n−2) is the orthogonal matrix of the eigenvectors of the matrix in Eq. (8) with T = I. The above system simplifies to
d2 v i (x) dx2
+ λi v i (x) = gˆ (x)c i
(10)
where v T = [ v 2 (x) v 3 (x), .., v ne (x)] T , λi are the eigenvalues of the matrix in Eq. (8), and c = T d. In addition, the boundary conditions in the modal space are given by
dx
|x=0 =
|x=1 =
T
dx de(x) dx
dx2
(11)
+ λi v i (x) = gˆ (x)c i ,
v i (x)|x=0 = 0,
dv i (x) dx
|x=1 = ˆi P. I,
v i (x)|x=1 = 0,
P. II,
(13)
⎢ 12 ⎢ x ⎣ ... ...
−2
x2
0
+ λi ... 1
2 x
1
x2
...
−4 2 x
A1
⎡
0
⎡
⎤
⎤ ⎡ v1 ⎤ ... i ⎢ 2⎥ ... ⎥ vi ⎥ ⎥⎢ ⎢ ⎥ ... ⎦⎣ ... ⎦
...
3 2 x
vn
i w
⎤
0 ⎢ gˆ (x2 ) ⎥
⎢0⎥ ⎢ ⎥ ⎥ : ⎥ c i . P. I, =⎢ ⎣ : ⎦+⎢ ⎣ gˆ (x ) ⎦ ne ˆi 0 b1
(14)
gˆ
where, for simplicity, we are also defining new variables according to w T = [ v 1i , v 2i , ..., v ni ] T , b1T = [0, 0, ..., ˆi ] T , and gˆ T = [0, gˆ (x2 ), .., gˆ (xne ), 0] T , and similarly, A1 . Note that the last row (equation) enforces the slope condition at x = 1. Similarly, a finite difference approximation to the second problem is given by
⎡
−3 2 x 1
x2
⎢ ⎢ ⎢ ⎢ ⎣ ... ...
x2
...
... 0
⎡
i
⎤
1
vn
i w
A2
⎡
0
1
...
⎤⎡
⎤ v1 ⎥⎢ i ⎥ ... ⎥ v 2i ⎥ ⎥⎢ ⎥ ⎥⎢ ⎣ ... ⎦ ... ⎦
−1
x
4 2 x −2 + λi
x2
0
⎤
gˆ (x2 ) ⎥ ⎢0⎥ ⎢ ⎥ ⎥+⎢ =⎢ c i P. II. ⎢ ⎣ : ⎦ ⎣ : ⎥ ⎦ gˆ (x )
(15)
ne
0
b2
where, we are also defining new variables and ˆ for simplicity. The above system is now fully decoupled. Every equation in the system is a second order differential equation that is required to satisfy two boundary conditions at each end. Each equation also includes an unknown constant c i . Therefore we are faced with a system of ill-posed linear boundary value problem. After solving the above equations and obtaining the values of unknown constants c i , the correction to the assumed value of the source term can be computed using the transformation d = c. It is possible to apply the method used in [23] to obtain solutions to the above boundary value problems. The approach is to consider two well-posed boundary value problems and note that they describe the same solution. In other words, instead of Eqs. (10) and (11), we consider two well-posed problems given by
d2 v i (x)
1
|x=0 = , ˆ |x=1 = ,
|x=0 = i
0
v(x)|x=0 = v(x)|x=1 = 0, dv(x) T de(x) dx dv(x)
+ λi v i (x) = gˆ (x)c i ,
for i = 1, 2, ..., (n − 2). Note that there are (n − 2) decouple equations. Both of the above problems are well-posed, and one can proceed to obtain solutions using finite difference approximations. If one divides the x domain into ne equal intervals and use central differencing to approximate the second derivative they lead to
⎡
where, e(x) = [e 2 (x) e 3 (x), ., ene (x) ] ∈ R and d = [ ( y 2 ) ( y 3 ), . ( yne )] T ∈ R n−2 are the unknown error field and T
dx2 dv i (x) dx
⎥ ⎥ .. .. .. ⎥ ⎥ e = gˆ (x)d, ⎥ 1 −2 1 ⎦ 0 1 −2
T
d2 v i (x)
3709
(12)
where, for simplicity, we are again defining new variables according to b2T = [i , 0, ..., 0] T , and, similarly, A2 . The above two problems are well-posed and the coefficient matrices can be stably inverted. Eliminating the unknown vector, w, from the above two equations leads to 1 −1 −1 −1 ˆ ˆ A− 1 b1 + A1 gc i = A2 b2 + A2 gc i ,
−1
−1
−1
−1
(A1 b1 − A2 b2 ) = (A2 − A1 )ˆg c i ,
r
or, (16)
s
where, we are also defining new vectors r and s for simplicity and space considerations. It is now possible to solve for the unknown scalar according to
ci = β
sT r sT s
,
(17)
where, the positive scalar β is used 0 < β and is set by the designer. We will explain this parameter in more details. Once the scalars c i , i = 1, 2, ..., (n − 2) are computed, then the correction to the assumed value of the source term can be solved for using the transformation d = c. We can apply a similar procedure to update the assumed value of the source term in the x direction.
3710
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
d2 v i (x)
+ λi v i (x) = gˆ (x)c i , dx2 v i (x)|x=0 = 0, v i (x)|x=1 = 0 d2 v i (x) dx2 dv i (x) dx
P. I,
(22)
+ λi v i (x) = gˆ (x)c i ,
|x=0 = ki
v i (x)|x=1 = 0,
P. II.
(23)
Both of the above problems are well-posed and one can follow the same formulation and obtain the unknown constant c i . The same procedure can be used to recover the correction in the x direction. We next proceed to consider inverse source problem in 3-D.
Fig. 2. Error field with over-specified boundary conditions.
3. Inverse source problem in 3-D Consider a bounded domain in R 3 , with a continuous boundary ∂. Let S (x) = 1 + g (x) f ( y )h( z) represent the source function in with the non-uniform portion being separable. The goal is to recover the source term for an elliptic system given by
u xx + u y y + u zz = 1 + g (x) f ( y )h( z)
2.2. Update in x directions After updating the source term in the y direction, it is possible to assume that f ( y ) = ˆf ( y ). Then equation (4) for the error field simplifies to
in .
(18)
The background field satisfies the Dirichlet condition, and therefore e = 0 on the boundary ∂. The data is in terms of flux at the boundaries as shown in Fig. 2. The same procedure can be used to recover the correction in the x direction. 2.3. Recovery with partial data It is possible to apply the present method when the data is collected at only two sides of the domain. Consider the same inverse problem presented in eq. (2) and assume that the data is collected at only two sides of the domain. Following similar steps and obtaining the error field leads to
e xx + e y y = g (x) f ( y ) − gˆ (x) ˆf ( y ),
in ,
e = 0 on
∂, (19)
uˆ xx + uˆ y y + uˆ zz = 1 + gˆ (x) ˆf ( y )hˆ ( z) in
(25)
.
Subtracting Eqn. (25) from (24) leads to the error field given by
e xx + e y y + e zz = g (x) f ( y )h( z) − gˆ (x) ˆf ( y )hˆ ( z) in
,
(26)
for which we have both Dirichlet and Neumann boundary conditions. The actual source functions are related to the assumed values according to
g (x) = gˆ (x) + σ (x),
and,
f ( y ) = ˆf ( y ) + ( y ),
h( z) = hˆ ( z) + η( z).
(27)
Similar to the problem in 2-D, it is possible to recover the corrections to the assumed values one term at a time. We can proceed as follows. 3.1. Update in z direction We can assume that g (x) ≈ gˆ (x) and f ( y ) ≈ ˆf ( y ). Then, the equation for the error simplifies to
e xx + e y y + e zz = gˆ (x) ˆf ( y )η( z)
in
.
(28)
Since e (x, y , z) = 0 at the boundaries x ∈ ∂, it is possible to assume an expansion of the form given by
in addition to
e x (0, y ) = κ ( y ),
(24)
,
where now is a unit cube. The boundary is accessible and it is possible to measure both Neumann and Dirichlet boundary conditions on the boundary ∂. Starting with an assumed value for the source term gˆ (x) ˆf ( y )hˆ ( z) and using the Dirichlet condition, it is possible to obtain a background field given by
Fig. 3. Error field with partial data at the boundary.
e xx + e y y = σ (x) ˆf ( y )
in
e y (x, 0) = ω(x).
(20)
For updating in y direction, one can assume that g (x) ≈ gˆ (x) and arrive at the system given in Fig. 3. After discretizing the second derivative and using the eigenvectors to decouple the working equations one can arrive at similar equations given by (10) with boundary conditions given by
e (x, y , z) =
∞
A (x, y )sin(π z),
=1
1 A (x, y ) = 2
e (x, y , z)sin(π z)dz,
(29)
0
v(x)|x=0 = v(x)|x=1 = 0,
dv(x) dx
|x=0 = T
de(x) dx
|x=0 = k.
(21)
Similar to eqs. (12) and (13), it is possible to formulate two wellposed problems for the error field with partial data according to
where the superscript simply denotes the component in the function expansion. By selecting sin(·) functions the Dirichlet boundary conditions in the z direction are satisfied. Substituting the above expansion in the error equation leads to
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
3711
Fig. 4. Two well-posed problems for the components of the error field.
∞
A xx + A y y − (π )2 A sin(π z) = gˆ (x) ˆf ( y )η( z)
in
.
=1
1 A y (x, 1) = ζ (x) = 2
(30) Using the orthogonality property of the sine functions, one can arrive at the following equation for each component
A xx + A y y − (π ) A = 2 gˆ (x) ˆf ( y )
2
1
η(z)sin(π z)dz in , 0
(31) for = 1, 2, ..∞. Retaining the first L components leads to
A xx + A y y − (π )2 A = gˆ (x) ˆf ( y )θ ,
1 θ = 2
η(z)sin(π z)dz, in ,
(32)
0
for = 1, 2, .., L. The collected data can provide Neumann boundary conditions for the above equations which are often called modified Helmholtz equation. The above equations are still ill-posed in the sense that every A (x, y ) is required to satisfy both Dirichlet (zero) and Neumann conditions. In addition, every equation includes an unknown scalar θ . Once these scalars are computed, then the correction to the assumed value of the source function in the z direction can be computed according to
η( z) =
∞
θ sin(π z) ≈
=1
L
θ sin(π z).
(33)
=1
In order to solve the above modified Helmholtz equations, one can follow the approach that was developed for the 2-D problem in section 2.1 and consider two well-posed problems given by where = ( ∂∂x2 + ∂ ∂y2 − (π )2 ) is the modified Helmholtz operator. The collected data appears as boundary conditions according to
e x (0, y , z)sin(π z)dz,
(34)
0
A x (1, y ) = ξ( y ) = 2
e x (1, y , z)sin(π z)dz,
(35)
e y (x, 0, z)sin(π z)dz,
(36)
0
1 A y (x, 0) = ν (x) = 2 0
The above problems are both well-posed, and one can obtain finite dimensional approximations using central differencing. Similar to Eqs. (14) and (15) finite difference approximations of the above problems lead to
A3 a = b3 + dˆ θ ,
(38)
A4 a = b4 + dˆ θ ,
(39)
where, A3 and A4 are the finite difference approximation of the Helmholtz operators in problem I and problem II shown in Fig. 4, respectively. The vector a is the nodal values of the unknown variable A (x, y ), b3 and b4 are the known right-hand sides of the finite difference approximations to problems I and II shown in Fig. 4, and the vector dˆ is the nodal values of the known product gˆ (x) ˆf ( y ). Note that the collected data, which appear as Neumann boundary conditions, show up in the vectors b3 and b4 . Also, the assumed values for the source function gˆ (x) ˆf ( y ) appear in the ˆ Both problems are well-posed, as a result A−1 and A−1 vector d. 3 4 exist. It is possible to eliminate e from the above equations and arrive at 1 1 (A−1 b3 − A− b4 ) = (A− − A−1 )dˆ θ , 3
4 4 3 w
(40)
v
where for simplicity and space considerations, we are redefining the vectors w and v. Solving for the scalar θ leads to
θ = β
wT v
(41)
vT v
where, the scalar β > 0 is chosen by the designer. After computing the scalars θ , for = 1, 2.., L, the correction to the assumed value can be obtained using Eq. (33). The updating in other directions can also be accomplished using a similar approach.
For updating in the y direction, we can assume that g (x) ≈ gˆ (x) and h( z) ≈ hˆ ( z). Then, the equation for the error simplifies to
e xx + e y y + e zz = gˆ (x) ( y )hˆ ( z)
1
(37)
3.2. Update in y and x directions
1 A x (0, y ) = μ( y ) = 2
e y (x, 1, z)sin(π z)dz. 0
in .
(42)
The above equation is similar to Eq. (28), and the same approach can be applied. Likewise, for updating in the x direction, we can assume that f ( y ) ≈ ˆf ( y ) and h( z) ≈ hˆ ( z). Then, the equation for the error simplifies to
e xx + e y y + e zz = σ (x) ˆf ( y )hˆ ( z) in
.
(43)
3712
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
Again, the above equation is similar to Eq. (28), and the same approach can be applied. 4. Helmholtz inverse source problem The same procedure described in section 3 can be applied to the inverse source problem for the Helmholtz equation. Starting with an initial guess and following the algorithm leads to the equation for the error field given by
e xx + e y y + e zz + 2 e = g (x) f ( y )h( z) − gˆ (x) ˆf ( y )hˆ ( z)
in . (44) Fig. 5. Reduction in the error for the Example 1.
Relating the actual values to the assumed values according to Eq. (27) and proceeding to recover the correction terms for one direction at a time lead to
e xx + e y y + e zz + 2 e = gˆ (x) ˆf ( y )η( z)
in .
(45)
Following the same procedure, the over-specified modified Helmholtz equation given in eq. (32) simplifies to
A xx + A y y + ( 2 − (π )2 ) A = gˆ (x) ˆf ( y )θ ,
1 θ = 2
η(z)sin(π z)dz,
(46)
0
for = 1, 2, .., L. The above equation is similar to Eq. (32) except for the nonzero wave number . For low to moderate wave numbers the above equations can be solved very accurately. The same procedure can be applied to solve the above ill-posed problem. We next apply the above algorithms to a number of numerical examples.
Fig. 6. Convergence of the recovered source function in the x direction for Example 1.
5. Numerical experiments In this section we use a number of numerical examples to study the applicability of the method. Example 1. First consider an inverse source problem in two dimensional domain. Assume that the actual source term in Eq. (2) is given by
g (x) f ( y ) = 1 + exp
(x − .3)4 0.000001
1 + 2exp
( y − .65)4 0.00001
. (47)
We can divide the domain in each direction into 100 equal intervals. The data is collected in the form of Neumann condition at the boundaries. The data is also contaminated with noise. The noise is a zero mean Gaussian white noise with signal-noise ratio of 6%. Before using the data a three point averaging is performed to smooth out the data. The actual source function has a nominal which is equal to 1. Since its value can be measured at the boundary, it is reasonable to start the iterations from gˆ (x) = ˆf ( y ) = 1. The algorithm seeks to recover the correction to the assumed values first in the y direction, and then in x direction. It is instructive to consider Eq. (6) which is used to update the source term in the y direction. It is presented here again
e xx + e y y =
gˆ (x)
( y ) in .
previous iteration
A similar equation, i.e., Eq. (18), is used to update the source term in the x direction. Considering the above equation it is clear that, in order to be able to recover the unknown function ( y ), one
Fig. 7. Convergence of the recovered source function in the y direction for Example 1.
needs to have the function gˆ (x) away from zero. For this example, the nominal value of the field is equal to one, and the inversion can be accomplished. As a result for this example, the parameter β in Eq. (17) is set equal to one, i.e., β = 1. If the nominal value of the field is close to zero, which is the case in most applications, then the above equation tends to over-estimate the unknown function. As a result a smaller value of β should be used. This will become clear in the next example. Fig. 5 presents the reduction in the error for this example. The error is the difference between the measured data and the computed value at the boundaries given by
Error =
∇n u − ∇n uˆ
2
dx.
∂
Figs. 6 and 7 show the recovered functions after a few iterations and compare them to the actual source functions. For this example, the algorithm can recover the source function very closely after only a few iterations. We next assume that only partial data is available for inverse evaluation of the source function.
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
3713
Fig. 11. Reduction in the error for the Example 2.
Fig. 8. Reduction in the error for the Example 1 with partial data.
Fig. 12. Actual source function for the Example 2. The background field is ≈ 0.0025.
Fig. 9. Convergence of the recovered source function in the x direction for Example 1 with partial data.
Fig. 13. Recovered source function for the Example 2 after 140 iterations.
Fig. 10. Convergence of the recovered source function in the y direction for Example 1 with partial data.
We consider the same problem and assume that data is available on only two sides given by (20). Fig. 8 shows the reduction in the total error. Only the error given in (20) is used in the recovery. With this limited data, it is still possible to formulate two well-posed problems given in eqs. (22) and (23). Figs. 9 and 10 show convergence of the source function when only partial data is available. Example 2. We next consider the recovering of a source function with a nominal value of the field close to zero. In this example, the actual source term is given by
(x − .3)4 0.000001 ( y − .65)4 . × 0.05 + 2exp 0.00001
g (x) f ( y ) = 0.05 + exp
(48)
For this example, since the nominal value of the field is very close to zero (0.05, and 0.0025 in 2-D) the algorithm over-estimates the
sought after function. It is then necessary to use a smaller value for the parameter β . For this case we use a value of β = 0.02. Using a smaller value for β will simply require more iterations. Fig. 11 shows the reduction in error for this case. The same amount of noise is used to model a realistic data. Figs. 12 and 13 show the actual source function and the recovered source function after 140 iterations. The present algorithm can recover a very close estimate of the source function. We next proceed to increase the noise in the data. Fig. 14 shows the reduction in error when the level of noise is increase by a factor of five. Fig. 15 shows the recovered source function after 120 iterations. The present algorithm can still recover a reasonably close estimate of the source function. Example 3. We next consider recovering of a source function in 3-D. Consider recovering a source function that is given by
g (x) f ( y )h( z)
(x − .3)4 ( y − .75)4 .1 + exp = .1 + exp 0.00002 0.00002 4 ( z − .5) . × .1 + exp 0.00002
(49)
For this case, the background field is 1.0 and the nominal (or apparent) value of the function is 0.001, and there is one target that is located at (x, y , z) ≈ (0.3, 0.75, 0.5). We can divide the domain into 50 equal intervals in each direction. This size mesh is required to calculate the field variables in 3-D with high accuracy.
3714
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
Fig. 14. Reduction in the error for the Example 2 with higher noise level. Fig. 18. Comparison of the actual function and the recovered source function in x for the Example 3.
Fig. 15. Recovered source function for the Example 2 after 120 iterations.
Fig. 19. Comparison of the actual function and the recovered source function in y for the Example 3.
In order to project the error down to the space of sine functions in Eq. (29), we use 40 terms in the series. However, when recovering the function back using Eq. (33) we use the first L = 18 terms. This is due to the fact that the noise in the data pollutes the higher frequency components and including them leads to instability. The parameter β = 30. Fig. 16 shows the location of the target within the 3-D domain. Fig. 17 presents the location of the recovered function. Figs. 18 through 20 present the individual components of the source function and compare them to the actual function. For this case, 60 iterations are used to recover the target presented in Fig. 16. Fig. 16. Actual target within a 3-D domain for the Example 3.
Example 4. We next consider recovering of multiple targets in a 3-D domain. Consider the recovering of a source function given by S (x) = 1 + g (x) f ( y )h( z) where,
(x − .3)4 (x − .65)4 + exp g (x) = .1 + exp 0.00002 0.000015 4 ( y − .75) ( y − .65)4 + exp f ( y ) = .1 + exp 0.00001 0.000015 4 ( z − .4) ( z − .7)4 + exp . h( z) = .1 + exp 0.00001 0.000015
(50) (51) (52)
The above source function represents four targets within the domain as is shown in Fig. 21. Fig. 22 shows the recovered targets after 60 iterations. All of the parameters are set to be the same as in Example 3. Also, the same amount of noise is used. The algorithm can recover multiple targets within the source function.
Fig. 17. Recovered target within a 3-D domain for the Example 3 after 60 iterations.
Example 5. We next consider recovering of multiple targets in a 3-D domain for the Helmholtz operator with the wave number
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
Fig. 20. Comparison of the actual function and the recovered source function in y for the Example 3.
Fig. 21. Actual source function with 4 targets for the Example 4.
3715
Fig. 23. Recovered source function after 80 iterations for the Helmholtz operator in Example 5.
A number of 2-D and 3-D examples were used to study the applicability of the proposed algorithm. We considered both Poisson and Helmholtz operators. The proposed method is particularly suitable for identifying a single anomaly with compact support. For the cases that the actual perturbation cannot be modelled by a separable function, the proposed method will converge to the closest approximation of that anomaly that can be described by a separable function. This is not a major weakness, simply because the uniqueness issues that are inherent for both of the above problems limits the actual source functions that can be recovered. The proposed method considers a unit square or cube. For other geometries the method can be modified in a number of ways. One possible approach is to use the method of immersed boundary to convert an irregular geometry into a cube. Another approach that can be used is the exact mapping of a general domain onto a square/cube. For this method the working equations in the transformed space will include a number of terms which will not be amendable to diagonalization. We envision that for these cases, it will be possible to set up an internal iteration to handle such terms. For geometries that are close to cube/square variable mesh finite difference method can be used. Again, an internal iteration can be used to handle the terms that are not amendable to diagonalization. All of these issues will be addresses in future works. References
Fig. 22. Recovered source function after 60 iterations for the Example 4.
2 = 4.5. The unknown source function is the same as in Example 4. Fig. 23 presents the recovered function after 80 iterations. All the parameters and the level of noise are set equal to their values in Example 4. 6. Conclusion In this paper we presented a computational algorithm for the inverse evaluation of the source function for Poisson and Helmholtz equations. The source term is assumed to be composed of a uniform background value and a separable function. The algorithm is iterative in nature. It updates the source function for each coordinate separately. The algorithm shows reasonable robustness to noise. It can also recover the source function with partial data.
[1] K. Ren, Recent developments in numerical techniques for transport-based medical imaging methods, Commun. Comput. Phys. 8 (1) (2010) 1–50. [2] S. Li, W. Driessen, S. Sullivan, H. Jiang, Bioluminescene tomography based on phantoms with different concentrations of bioluminescent cancer cells, J. Opt. A, Pure Appl. Opt. 8 (2006) 743–746. [3] B.J. Minsley, J. Sogade, F.D. Morgan, Three-dimensional source inversion of selfpotential data, J. Geophys. Res. 112 (2007) B02202. [4] A. Artemev, L. Parnovski, I. Polterovich, Inverse electrostatic and elasticity problems for checkered distributions, Inverse Probl. 29 (2013) 075010, 16 pp. [5] M. Slodicka, H.D. Schepper, On an inverse problem of pressure recovery arising from soil venting facilities, Appl. Math. Comput. 129 (2002) 469–480. [6] V. Isakov, Inverse Problems for Partial Differential Equations, Springer-Verlag, Berlin, 1998. [7] C.-S. Liu, A BIEM using the Trefftz test functions for solving the inverse Cauchy and source recovery problems, Eng. Anal. Bound. Elem. 62 (2016) 177–185. [8] A. Karageorghis, D. Lesnic, L. Marin, A survey of applications of the MFS to inverse problems, Inverse Probl. Sci. Eng. 19 (3) (2011) 309–336. [9] Y. Hon, M. Li, Y. Melnikov, Inverse source identification by Green’s function, Eng. Anal. Bound. Elem. 34 (2010) 352–358. [10] A.E. Badia, An inverse source problem in an anisotropic medium by boundary measurements, Inverse Probl. 21 (2005) 1487–1506. [11] N. Roberty, C. Alves, On the identification of star-shape sources from boundary measurements using a reciprocity functional, Inverse Probl. Sci. Eng. 17 (2) (2009) 187–202. [12] S.S. Adavani, G. Biros, Fast algorithms for source identification problems with elliptic PDE constraints, SIAM J. Imaging Sci. 3 (4) (2010) 791–808.
3716
A. Hamad, M. Tadi / Physics Letters A 380 (2016) 3707–3716
[13] A. Canelas, A. Laurian, A. Novotny, A new reconstruction method for the inverse potential problem, J. Comput. Phys. 268 (2014) 417–431. [14] G. Wang, Y. Li, M. Jiang, Uniqueness theorems in bioluminescence tomography, Med. Phys. 31 (8) (2004) 2289–2299. [15] B.J. Last, K. Kubik, Compact gravity inversion, Geophysics 48 (6) (1983) 713–721. [16] O. Portniaguine, M.S. Zhdanov, Focusing geophysical inversion images, Geophysics 64 (3) (1999) 874–887. [17] T. Ohe, K. Ohmaka, Boundary element approach for an inverse source problem of the Poisson equation with a one-point-mass-like source, Appl. Math. Model. 18 (1994) 216–223. [18] A.E. Badia, T. Doung, An inverse source problem in potential analysis, Inverse Probl. 16 (2000) 651–663.
[19] N. Martins, An iterative shape reconstruction of source functions in a potential problem using the MFS, Inverse Probl. Sci. Eng. 20 (8) (2012) 1175–1193. [20] B. Jin, L. Marin, The method of fundamental solutions for inverse source problems associated with the steady-state heat conduction, Inverse Probl. Sci. Eng. 69 (2007) 1570–1589. [21] N.G. Gencer, I. Tanzer, Forward problem solution of electromagnetic source imaging using a new BEM formulation with high-order elements, Phys. Med. Biol. 44 (1999) 2275–2287. [22] R. Celorrio, A. Mendioroz, A. Salazar, Characterization of vertical buried defects using lock-in vibrothermography: II. Inverse problem, Meas. Sci. Technol. 24 (2013) 065602, 9 pp. [23] M. Tadi, A.K. Nandakumaran, S. Sritharan, An inverse problem for Helmholtz equation, Inverse Probl. Sci. Eng. 19 (6) (2011) 839–854.