The method of fundamental solutions for complex electrical impedance tomography

Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

The method of fundamental solutions for complex electrical impedance tomography Marjan Asadzadeh Heravi a, Liviu Marin b,c,n, Cristiana Sebu a a

Department of Mechanical Engineering and Mathematical Sciences, Oxford Brookes University, Wheatley Campus, Oxford OX33 1HX, UK Department of Mathematics, Faculty of Mathematics and Computer Science, University of Bucharest, 14 Academiei, 010014 Bucharest, Romania c Institute of Solid Mechanics, Romanian Academy, 15 Constantin Mille, 010141 Bucharest, Romania b

art ic l e i nf o

a b s t r a c t

Article history: Received 18 January 2014 Received in revised form 11 March 2014 Accepted 25 April 2014

The forward problem for complex electrical impedance tomography (EIT) is solved by means of a meshless method, namely the method of fundamental solutions (MFS). The MFS for the complex EIT direct problem is numerically implemented, and its efficiency and accuracy as well as the numerical convergence of the MFS solution are analysed when assuming the presence in the medium (i.e. background) of one or two inclusions with the physical properties different from those corresponding to the background. Four numerical examples with inclusion(s) of various convex and non-convex smooth shapes (e.g. circular, elliptic, peanut-shaped and acorn-shaped) and sizes are presented and thoroughly investigated. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Electrical impedance tomography (EIT) Multi-frequency Forward problem Method of fundamental solutions (MFS) Meshless method

1. Introduction Electrical impedance tomography (EIT) is a technique used for determining the admittivity distribution in the interior of an object, given simultaneous measurements of alternating electric currents (of frequencies varying from 10 Hz to 500 kHz) and of induced voltages on the boundary of the object [12]. For a given conductive object, the admittivity is a complex valued function whose real part is the electrical conductivity and whose imaginary part is the product of the frequency of the applied electric alternating current and the permittivity of the object. Since different materials display different electrical properties, a map of internal admittivity can be used to infer the internal structure of the object under consideration. Therefore, EIT can be used as a non-invasive and portable method of industrial, geophysical and medical imaging [1]. The reconstruction procedures proposed for static EIT include a wide range of iterative methods based on formulating the inverse problem in the framework of nonlinear optimization [2,24,40,41]. These approaches usually involve estimating the admittivity distribution of the object under consideration and then solve the forward problem (often using finite element methods [FEMs]) for the same

n Corresponding author at: Department of Mathematics, Faculty of Mathematics and Computer Science, University of Bucharest, 14 Academiei, 010014 Bucharest, Romania. E-mail addresses: [email protected] (M.A. Heravi), [email protected], [email protected] (L. Marin), [email protected] (C. Sebu).

input current patterns to compute the boundary voltages and then comparing the boundary data predicted by this estimate with the measured data. The discrepancy between these two data sets is then used to update the admittivity estimate and the procedure is repeated until a satisfactory agreement is achieved. However, the solution of the inverse problem in static EIT has suffered not only from its illposedness due to the inherent insensitivity of boundary measurements to any small changes of interior conductivity and permittivity values and its poor spatial resolution, but also from its reliance on accurate forward models which mimic every aspect of the imaging object (e.g. knowledge of boundary geometry, electrode positions and other sources of systematic artifacts in measured data). Hence, EIT has had limited applicability so far in clinical applications. This has encouraged the search of new reconstruction methods, such as the time-difference EIT (tdEIT) or frequency-difference EIT (fdEIT). Even though numerous tdEIT methods have been applied to image lung functions, stomach emptying or brain functions [12,31,32], there are medical applications where these time-reference data are not available (e.g. breast cancer or cerebral stroke detection). Since complex conductivity spectra of biological tissues show frequency-dependent changes [8,32], fdEIT methods have been proposed to image the changes in the admittivity distribution with respect to frequency [11,13,22,38]. It has been showed that, although there are no visible differences in the reconstructed frequency-difference images even though the true admittivity distributions have a strong frequency dependence, any anomaly can be clearly identified as long as its admittivity differs significantly from that of the background which is the case for tumour and stroke imaging.

http://dx.doi.org/10.1016/j.enganabound.2014.04.022 0955-7997/& 2014 Elsevier Ltd. All rights reserved.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

devices but also the development of a computer software fast enough to be used for real-time monitoring. However, optimized image reconstruction techniques for EIT rely on computationally efficient and numerically robust forward solvers. In this paper, we address this need by presenting an algorithm based on the method of fundamental solutions (MFS) which can be successfully used to find the numerical solution of the EIT forward problem for piecewise constant admittivity distributions to a high level of precision. The MFS is a meshless boundary collocation method applicable to boundary value problems in which a fundamental solution of the operator in the governing equation is known explicitly. The basic ideas of this method were first introduced, in the early 1960s, by Kupradze and Aleksidze [23], whilst its numerical formulation was first given by Mathon and Johnston [30] in the late 1970s. The main idea of the MFS consists of approximating the solution of the problem by a linear combination of fundamental solutions with respect to some singularities/source points which are located outside the domain. Consequently, the original problem is reduced to determining both the unknown coefficients of the fundamental solutions and the coordinates of the source points by requiring the approximation to satisfy the boundary conditions in some sense and hence solving a non-linear problem. If the source points are fixed a priori, then the coefficients of the MFS approximation are determined by solving a linear problem. The aforementioned MFS procedures are referred to in the literature as the dynamic and static approaches, respectively. Despite its constraint on the knowledge of a fundamental solution of the governing partial differential equation, the MFS

To achieve clinical acceptance, the theoretical developments of EIT reconstruction methods need to be closely connected with laboratory experiments and studies on real data. This implies not only the modelling of the real geometry and data collection

Fig. 1. Geometry of the problem. Possible placement of the sources for Ω0 ( ) and Ωj, j ¼ 1; 2, ( ), and the collocation points ().

20

30 4.14E-04

6.82E-03

9.59E-03

15

25

1.03E-03

1.74E-02 2.48E-03

20

2.91E-02

N01

N01

4.99E-03 5.10E-02

10

2.56E-02

2.03E-01

5

9.41E-03

15

1.04E-01

2.78E-01

10

7.32E-02

3.22E-01

2.04E-01 2.88E-01

5

3.22E-01

2.88E-01

2.78E-01

5

10

15

20

25

30

35

40

10

20

N00

30

40

50

60

N00

40 2.88E-05 9.30E-05

30

2.85E-04

N01

7.62E-04

20

2.54E-03 7.74E-03 2.56E-02

10

1.04E-01

2.80E-01 2.80E-01

10

20

30

40

50

60

70

80

N00 Fig. 2. The RMS error, errΓ 1 ðv0  v1 Þ, as a function of the number of sources N 00 A f1; 2; …; M 00 g and N 01 A f1; 2; …; M 01 g at t ¼0 s and various numbers of collocation points, namely (a) M ¼ 80, (b) M¼ 120 and (c) M ¼ 160, for Example 1.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

The MFS has been successfully applied to a large variety of physical problems, such as the Laplace equation [14,25,35,37], the biharmonic equation [16,33], the Helmholtz and modified Helmholtz [15,26,27], scattering and radiation problems [6], elastostatics [17,19,21,36], thermoelasticity [28,39], inhomogeneous elliptic equations [34], heat conduction problems in functionally graded materials [29], layered [3] and composite materials [4,18]. For a detailed and comprehensive description of the application of the MFS to direct and inverse problems we refer the reader to the survey papers by Fairweather and Karageorghis [5], Fairweather et al. [6], Golberg and Chen [9], and Karageorghis et al. [20], as well as the references therein. To the best of our knowledge, the application of the MFS to the forward problem for complex EIT has not been investigated yet. Motivated by the success of the application of the MFS to various physical problems, as well as the importance of the complex EIT forward problem, the purpose of the present study is to propose, implement and analyse the application of the MFS to the aforementioned problem. We also investigate the efficiency, accuracy and numerical convergence of the MFS solution of the complex EIT forward problem when assuming the presence in the medium (i.e. background) of one or two inclusions with physical properties (i.e. conductivity and permittivity) different from those of the background. The numerical analysis is studied for various convex and non-convex smooth shapes (e.g. circular, elliptic, peanut-shaped and acorn-shaped) and sizes of the inclusion(s). The paper is structured as follows. In Section 2 we formulate mathematically

has become very popular especially due to the ease with which it can be implemented, in particular for problems in complex geometries. It should be mentioned that the MFS can be used in conjunction with the so-called singularity subtraction technique (SST) for problems with boundary and/or solution singularities, in the sense that the standard MFS solution of the partial differential equation investigated is augmented by suitable singular functions/ eigenfunctions associated with the corresponding partial differential operator as given by the asymptotic expansion of the solution near the singular point. For details on the MFS-SST procedure applied to the numerical solution of direct and inverse problems associated with the Laplace, biharmonic, Helmholtz and modified Helmholtz, and Cauchy-Navier equations, we refer the reader to Karageorghis [14], Karageorghis et al. [21], and Marin [25–27]. In the MFS, there are no integrations. This clearly represents an advantage over the boundary element method (BEM) where integrations could be potentially troublesome and complicated. Furthermore, the MFS is a boundary collocation method, and hence only the boundary of the solution domain needs to be discretised, which is another advantage of the MFS over domain discretisation methods, such as FEMs or finite difference methods (FDMs). The major disadvantages of the MFS are related to the optimal location of the singularities, which is an important issue especially from the numerical point of view, and its inability to be directly applied to non-homogeneous or nonlinear partial differential equations, for which a fundamental solution is not available.

30

20

2.29E-03

2.39E-02

25

3.97E-02

15

6.04E-02

N01

1.51E-01

15

3.69E-02 6.98E-02

2.00E-01

1.64E-0

1

N01

1.05E-02 1.74E-02

1.04E-01

10

5.02E-03

20

7.48E-02

5

3

2.34E-01

10

1.29E-01

1.62E-01

2.21E-01

2.34E-01

5

2.00E-01

2.21E-01 1.62E-01

1.51E-01

5

10

15

20

25

30

35

40

10

20

N00

30

40

50

60

N00

40 2.18E-04 6.20E-04

30

1.70E-03

N01

3.58E-03

20

1.01E-02 2.85E-02

6.91E-02

10

1.50E-01

1.98E-01

1.50E-01

10

20

30

40

50

60

70

80

N00 Fig. 3. The RMS error, errΓ 1 ½nð0Þ  ðs0 ∇v0  ωε0 ∇h0 Þ þ nð1Þ  ðs1 ∇v1  ωε1 ∇h1 Þ, as a function of the number of sources N 00 A f1; 2; …; M 00 g and N 01 A f1; 2; …; M 01 g at t¼ 0 s and various numbers of collocation points, namely (a) M ¼ 80, (b) M ¼ 120 and (c) M ¼160, for Example 1.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

additive constant [7], which we fix by choosing the ground as Z ϕðxÞ dx ¼ 0: ð3Þ

the problem investigated herein. The MFS application to the complex EIT forward problem is presented in Section 3. The efficiency, accuracy and numerical convergence of the proposed method are validated by considering four examples in Section 4. Finally, some conclusions and future related problems are presented in Section 5.

∂Ω

In this paper, we are interested to find the numerical solution to the EIT forward problem for a piecewise constant admittivity distribution. To this end, we let Ωj, j ¼ 1; J , be disconnected subdomains compactly contained in Ω, such that Ωk \ Ωl ¼ ∅ for all k; l ¼ 1; J , k al, and Ω0 ¼ Ω\⋃Jj ¼ 1 Ωj is connected, see Fig. 1. The following notations are also introduced

2. Mathematical formulation Let Ω be a bounded, simply connected domain in R2 . The forward complex EIT problem can be stated as follows: for a given complex admittivity scalar-valued function γ ðxÞ ¼ sðxÞ þ iωεðxÞ, x A Ω , with strictly positive real part sðxÞ Zc 4 0, at a frequency ω of the applied current, find the complex potential distribution ϕðxÞ ¼ vðxÞ þ ihðxÞ, x A Ω , satisfying the elliptic partial differential equation ∇  ðγ ðxÞ∇ϕðxÞÞ ¼ 0;

x A Ω;

J

Z such that

∂Ω

Γ j ¼ ∂ Ωj ;

ð4bÞ

j ¼ 1; J ;

ϕj ðxÞ ¼ vj ðxÞ þ ihj ðxÞ; γ j ¼ sj þiωεj ;

ð1Þ

x A ∂Ω

ð4aÞ

j¼1

x A Ω j ¼ Ωj [ ∂Ωj ; j ¼ 0; J ;

j ¼ 0; J ;

ð4cÞ ð4dÞ

where both the conductivity sj and the permittivity εj distribution in Ω j are assumed to be constant. Thus, Eq. (1) reduces to a set of Laplace equations

subject to the Neumann boundary condition

γ ðxÞðnðxÞ  ∇ϕðxÞÞ ¼ Iðx; ωÞ;

Γ 0 ¼ ∂Ω0 \ ⋃ ∂Ωj ;

Iðx; ωÞ dx ¼ 0;

ð2Þ

∇  ð∇ϕj ðxÞÞ ¼ 0;

x A Ωj ; j ¼ 0; J ;

ð5aÞ

which could be solved subject to the Neumann boundary condition

where nðxÞ is the outer normal at x A ∂Ω. If γ A L1 ðΩ Þ, the Neumann boundary value problem (1) and (2), for I A H  1=2 ð∂ΩÞ, has a unique solution ϕ A H 1 ðΩÞ up to an

x A Γ0

Z such that

Γ0

I 0 ðxÞ dx ¼ 0;

ð5bÞ

0.6

4.58E-04

0.4

8E

1.35E-06

1.29E-04

1.52E-04

d1 = d 2

d1 = d 2

4.09E-05

-05

0.8

1.66E-05

1.2

0.8

1.28E-05

0.6

1.28E-05 3.98E-05

0.4

1.81E-03

1.59E-04

8.20E-03

0.2

0.2

9.83E-03

6.31E-02

1.0

1.5

7.44E-05

7.60E-04

4.35E-02

0.5

1.59E-0

4

γ 0 ðnð0Þ ðxÞ  ∇ϕ0 ðxÞÞ ¼ I0 ðxÞeiωt ;

1.47E-02

2.0

2.5

3.0

0.5

1.0

2.0

2.5

3.0

d0

3.3

0E

05

5E-

-05

3.9

-06

8.4

7E-

1.30E

06

d1 = d 2

2.53

3.37E-05

0.8

0.6

E-0

4

d0

1.5

0.4

1.28E-05

3.37E-05 2.53E-04

0.2

8.56E-05

2.69E-03 4.03E-03

0.5

1.0

1.5

2.0

2.5

3.0

d0 Fig. 4. The RMS error, errΓ 1 ðh0  h1 Þ, as a function of d0 and d1 ¼ d2 at t ¼ 0:5  10  6 s and various numbers of collocation and source points, namely (a) M ¼ N¼ 80, (b) M ¼N ¼ 120 and (c) M ¼N ¼ 160, for Example 1.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Z

at a set frequency ω and time t. To guarantee uniqueness, we impose Z ϕ0 ðxÞ dx ¼ 0: ð5cÞ

Γ0

Due to the continuity of the potential and the current flux density across the interfaces Γj, the solutions ϕj, j ¼ 0; J , should also satisfy the following transmission conditions: x A Γ j ; j ¼ 1; J ;

h0 ðxÞ dx ¼ 0;

ð6fÞ

and the transmission conditions

Γ0

ϕ0 ðxÞ  ϕj ðxÞ ¼ 0;

ð5dÞ

v0 ðxÞ  vj ðxÞ ¼ 0;

x A Γ j ; j ¼ 1; J ;

ð6gÞ

h0 ðxÞ  hj ðxÞ ¼ 0;

x A Γ j ; j ¼ 1; J ;

ð6hÞ

nð0Þ ðxÞ  ðs0 ∇v0 ðxÞ  ωε0 ∇h0 ðxÞÞ þnðjÞ ðxÞ ðsj ∇vj ðxÞ  ωεj ∇hj ðxÞÞ ¼ 0;

γ 0 ðnð0Þ ðxÞ  ∇ϕ0 ðxÞÞ þ γ j ðnðjÞ ðxÞ  ∇ϕj ðxÞÞ ¼ 0; x A Γ j ; j ¼ 1; J ;

ðωεj ∇vj ðxÞ þ sj ∇hj ðxÞÞ ¼ 0;

where nðjÞ ðxÞ are the outer normals at x A ∂Ωj , j ¼ 0; J , see Fig. 1. Eqs. (5a)–(5e) can be rewritten in terms of vj ¼ Re ϕj and hj ¼ Im ϕj , j ¼ 0; J , as follows: x A Ωj ; j ¼ 0; J ;

ð6aÞ

Δhj ðxÞ ¼ 0;

x A Ωj ; j ¼ 0; J ;

ð6bÞ

x A Γ0;

ð6cÞ

nð0Þ ðxÞ  ðωε0 ∇v0 ðxÞ þ s0 ∇h0 ðxÞÞ ¼ I 0 ðxÞ sin ðωtÞ;

x A Γ0;

ð6dÞ

the uniqueness conditions Z v0 ðxÞ dx ¼ 0;

1 log J x  ξ J ; 2π

x A Ω j ; j ¼ 0; J ;

d1 = d 2

d1 = d 2

5.39E-04 1.89E-03

0.4

5 1.10E-0

0.6

4.37E-03

0.2

1.26E-01

4.73E-02

1.81E-01

0.5

1.0

1.5

7.03E-02

2.0

2.5

3.0

0.5

1.0

1.5

2.0

2.5

3.0

d0

2.31

2.31

0.8

5

E-04

3.12E-0

1.67E-05

E-04

d0

6.0

0.6

7E -06

3.85 E-0 5

d1 = d 2

5.60E-05

4.58E-04

3.00E-02

0.2

3.40E-05

3.40E-05 1.32E-04

0.4

7.23E-03

-05

0.6

0E

1.99E-04

1.1

5.55E-06

0.8

5.60E-05

4

5

ð7Þ

the Euclidean distance between the collocation point x and the singularity/source point ξ. In Fig. 1, we present the possible positions of the source and collocation points for a simple geometry.

1.49E-0

5.60E-0

0.8

ð6jÞ

where x ¼ ðx1 ; x2 Þ is a collocation point, ξ ¼ ðξ1 ; ξ2 Þ A R2 \Ω j is a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi singularity or source point and J x  ξ J ¼ ðx1  ξ1 Þ2 þ ðx2  ξ2 Þ2 is

ð6eÞ

Γ0

x A Γ j ; j ¼ 1; J :

The fundamental solution of the two-dimensional Laplace equation is given by, see e.g. [5] Fðx; ξÞ ¼ 

n ðxÞ  ðs0 ∇v0 ðxÞ  ωε0 ∇h0 ðxÞÞ ¼ I 0 ðxÞ cos ðωtÞ;

ð6iÞ

3. The method of fundamental solutions

subject to the Neumann boundary conditions ð0Þ

x A Γ j ; j ¼ 1; J ;

nð0Þ ðxÞ  ðωε0 ∇v0 ðxÞ þ s0 ∇h0 ðxÞÞ þnðjÞ ðxÞ

ð5eÞ

Δvj ðxÞ ¼ 0;

5

0.4

2.31E-04 6.74E-04

0.2

1.82E-02 2.73E-02

0.5

1.0

1.5

2.0

2.5

3.0

d0 Fig. 5. The RMS error, errΓ1 ½nð0Þ  ðωε0 ∇v0 þ s0 ∇h0 Þ þ nð1Þ  ðωε1 ∇v1 þ s1 ∇h1 Þ, as a function of the distances d0 and d1 ¼ d2 at t ¼ 0:5  10  6 s and various numbers of collocation and source points, namely (a) M ¼ N ¼80, (b) M ¼N ¼ 120 and (c) M ¼ N¼ 160, for Example 1.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

of the potential on ∂Ωj can be approximated by

In the MFS, both the real (vj) and the imaginary (hj) potential corresponding to the background (j¼0) and inclusions (j ¼ 1; J ) are approximated by a linear combination of fundamental solutions with respect to N0 ¼ N00 þ N 01 þ ⋯ þ N0J singularities, ðj;nÞ N j gn ¼ 1 ,

ð0;nÞ N 0

fξ gn ¼ 1 , and Nj singularities, fξ the form

h i Nj ðj;nÞ ðnðjÞ ðxÞ  ∇vj ðxÞÞ  ∑ cðj;nÞ nðjÞ ðxÞ  ∇x Fðx; ξ Þ ; n¼1 Nj

1; J , respectively, in

ðj;nÞ

ðnðjÞ ðxÞ  ∇hj ðxÞÞ  ∑ d

h

n¼1

ðj;nÞ

nðjÞ ðxÞ  ∇x Fðx; ξ

i Þ ;

x A ∂Ωj ; j ¼ 0; J ;

ð9aÞ

x A ∂Ωj ; j ¼ 0; J :

ð9bÞ

M

Nj

vj ðxÞ  ∑ cðj;nÞ Fðx; ξ

ðj;nÞ

Þ;

n¼1

Nj

hj ðxÞ  ∑ d n¼1

ðj;nÞ

Fðx; ξ

ðj;nÞ

Þ;

x A Ω j ; j ¼ 0; J ;

ð8aÞ

x A Ω j ; j ¼ 0; J ;

ð8bÞ

j Next, we select M j collocation points, fxðj;nÞ gn ¼ 1 , on each boundary Γj, j ¼ 0; J , such that the collocation points on each interface Γ j , j ¼ 1; J , are the same for both the boundary associated with the background, ∂Ω0 , and the boundary corresponding to the jth inclusion, ∂Ωj , i.e. the following relations hold: j1

xð0;∑ℓ ¼ 0 Mℓ þ kÞ ¼ xðj;kÞ ;

ðjÞ

d ¼ ½cðj;1Þ ; …; cðj;Nj Þ T A RNj , for j ¼0,..., J. It should be noted that, in the case of the background (j¼0), the total number of singularities (source points) was taken to be N 0 ¼ N 00 þ N 01 þ⋯ þ N 0J since ∂Ω0

is the outer boundary of the background and Γj, j ¼ 1; J , are the J inner boundaries of the background (interfaces between the

N0

M0

n¼1

k¼1

∑ cð0;nÞ

background and the inclusions), whilst N 00 and N 0j , j ¼ 1; J , are singularities corresponding to the outer and inner boundaries (interfaces), respectively. From Eqs. (6c), (6d), (6i), (6j), (8a) and (8b), it follows that the corresponding normal derivatives of the real and imaginary parts

ð0;nÞ

∑ Fðxð0;kÞ ; ξ

ð0;nÞ

M0

ð0;nÞ

ð11aÞ

#

∑ Fðxð0;kÞ ; ξ

n¼1

Þ  0;

Þ  0:

ð11bÞ

k¼1

1.0

+0 0

6E

0

x2

.9

+0

7E

+0

1

1

E+

00

83

E+

4. E+

00

91

00

E+

-0.5

-1

.5

4E

01

+0

0

E+

E+

98

0

1.

+0

36

00

15

01

2.

E+

0E

1.

E+

00

47

.0

.3

26

8.

E+

47

-2

-1

1.

0.0

00

01

E+

5.

11

9. 1.

01

01

-0.5

7E +0

0

1

1.

0.0

-1

5E

+0

+0

0.5

1

.0

0

6E

+0

.7

5E

+0

-5

3E

-5

.4

1E

.9

.7 -1

.1

-8

-2

-9

0.5

x2

"

N0

∑ d

1.0

-1.0 -1.0

ð10Þ

Consequently, the total number of boundary collocation points employed in the present approach is M 0 þ M 1 þ ⋯ þM J , i.e. the number of boundary collocation points associated with the backM 0 þ M1 þ ⋯ þ M J ground Ω0, namely fxð0;nÞ gn ¼ . 1 Then, the collocation of the uniqueness conditions (6e) and (6f) yields the following equations: " #

and we use the following notation cðjÞ ¼ ½cðj;1Þ ; …; cðj;Nj Þ T A RNj and

-0.5

k ¼ 1; M j ; j ¼ 1; J :

0.0

0.5

-1.0 -1.0

1.0

-0.5

x1

0.0

0.5

1.0

x1

1.0 -3

x2

00

E+

4E

+0

0

00

6E

00

-01

00

-0.5

E+

1.3 10

E+

1

6.8

E+

43

10

00

-1.0 -1.0

4.

0

-0.5

97

E+

2.

-0

01

0.0 2.

1E

+0

E-

.6

2E

86

.8

6.

-4.

.13

-1 -7

0.5

0.0

0.5

1.0

x1 Fig. 6. Distribution of the real part of the potential, v ¼ Re ϕ, for M ¼N ¼ 160 at various times, namely (a) t ¼0 s, (b) t ¼ 0:2  10  6 s and (c) t ¼ 0:5  10  6 s, for Example 1.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

with R¼ 1, and is characterised by the conductivity s0 ¼ 0:036 S m  1 and permittivity ε0 ¼ 1:9  10  9 F m  1 . In SI based 1 1 units, 1 S ¼ 1 A2 s3 m  2 kg and 1 F ¼ 1 A2 s4 m  2 kg . The inclusions considered in this study Ωj, j ¼ 1; 2, are assumed to have smooth boundaries ∂Ωj , j ¼ 1; 2, and be centered at xð1Þ ¼ ð0:45; 0Þ and xð2Þ ¼ ð  0:45; 0:35Þ, respectively, whilst their corresponding conductivities and permittivities are given by s1 ¼ 0:042 S m  1 and ε1 ¼ 9:7  10  9 F m  1 , and s2 ¼ 0:05 S m  1 and ε2 ¼ 1:2  10  8 F m  1 , respectively.

Finally, we collocate the Neumann boundary conditions (6c) and (6d), as well as the transmission conditions (6g)–(6j) to obtain the following ðjÞ system of linear equations for the unknown coefficients cðjÞ ; d A RNj , N j ¼ 0; J , grouped in the vector X A R , N ¼ 2ðN 0 þ N 1 þ ⋯ þ N J Þ: AX ¼ f;

ð12Þ ðM þ 2ÞN

is the corresponding MFS matrix whose elements where A A R are calculated from Eqs. (6c) to (6j), respectively, while f A RM þ 2 , M ¼ 2M 0 þ 4M 1 þ ⋯ þ 4M J , contains the corresponding discretised Neumann data from Eqs. (6c) and (6d), as well as zero entries associated with the transmission conditions (6g)–(6j) and uniqueness conditions (11a) and (11b). Clearly, in order for the discretised MFS system of linear equations (10) to have a unique solution X A RN , the numbers of collocation points M and singularities N should satisfy the inequality N r M þ 2. In this study, the system of linear equations (12) is solved in the least-squares sense by numerically inverting the normal system of equations associated with the system (12), namely X ¼ ðAT AÞ  1 ðAT fÞ:

Example 1. We consider the circular inclusion of radius r ¼0.25 and centered at xð1Þ , i.e. ð1Þ 2 2 2 Ω1 ¼ fx A R2 jðx1  xð1Þ 1 Þ þ ðx2  x2 Þ o r g:

ð13Þ

where θ A ½0; 2π Þ is the polar angle of point x with respect to the centre of the inclusion, a ¼1.0, b¼0.5, r ¼ 0.25, and X 1 and X 2 are given by

Without loss of generality, we assume that the background is described by the unit disk centered at the origin of the system of coordinates, i.e.

Ω ¼ fx A R

jx21 þ x22 oR2 g;

X 1 ¼ ðx1  x1ð1Þ Þ cos φ þ ðx2  x2ð1Þ Þ sin φ; ð1Þ X 2 ¼ ðx1  xð1Þ 1 Þ sin φ þ ðx2  x2 Þ cos φ:

ð14Þ

1.0

-8

3.1 0

00

E01

E+

00

x2

88

E+

00

.23

E+

00

00

00

0.0

-1

E+

01

0

E+

30

E-

-0.5

73

x2

+0

88

24 9.

86

+0

-0.5

1.

3. 6.

-6.

0

0

0E

0.0

00

+0

.1

E+

1.

00

+0

7E

0

00

0E

0 E+

E+

0E

.9

.9

+0 82

-2 -4

.34

.1

0E

+0

-6.

1.

1

-1

-2

0.5

3E

0

01

-2

4.4

+0

E-

00

-0

0.0

E+

6E

82

.8

61

.9

3E

1. 7.

-6

-1.0 -1.0

ð16bÞ

1.0

0.5

-0.5

ð15Þ

Example 2. We consider a peanut-shaped inclusion rotated by an angle φ ¼ 5π =6 about its centre xð1Þ  ( )  X 21 þ X 22 2 Ω1 ¼ x A R2  o r ; ð16aÞ a2 cos 2 θ þb2 sin 2 θ

4. Numerical results

2

7

0.5

-1.0 -1.0

1.0

-0.5

x1

0.0

0.5

1.0

x1

1.0 -2

0

+0

1

5E

+0

1

47

00

1.

E+

00

47

E+

9.

00

01

E+

5.

11

0.0

E+

15

01

2.

-0.5

1E

.4

1.

x2

00

E+

+0

6 .7 -5

0.5

3E

.7

-9

.1

-1

-2

.0

0E

+0

0

E+ 01

-1.0 -1.0

-0.5

0.0

0.5

1.0

x1 Fig. 7. Distribution of the imaginary part of the potential, h ¼ Im ϕ, for M ¼ N ¼160 at various times, namely (a) t ¼0 s, (b) t ¼ 0:2  10  6 s and (c) t ¼ 0:5  10  6 s, for Example 1.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

M

Example 3. We consider an acorn shaped inclusion centered at xð1Þ  ( ) ðx xð1Þ Þ2 þ ðx  xð1Þ Þ2 2 2 or 2 ; Ω1 ¼ x A R2  1 1 ð17Þ  a þ b sin ð3θÞ

j Generally speaking, the collocation points, fxðj;nÞ gn ¼ 1 , are uni-

Γj, j ¼ 0; J . The singularities/

formly distributed on each boundary sources fξ

ð0;nÞ N0 gn ¼ 1 ,

N 0 ¼ N00 þ N01 þ⋯ þ N 0J , associated with ∂Ω0 ¼

ðj;nÞ N j gn ¼ 1

⋃Jj ¼ 0 Γ j , and fξ

where θ A ½0; 2π Þ is the polar angle of point x with respect to the centre of the inclusion, a ¼4.25, b¼2.0 and r ¼0.1.

associated with

Γj, j ¼ 1; J , are also uni-

formly distributed. Moreover, we take M 0 ¼ N 00 and M j ¼ N j ¼ N0j , j ¼ 1; J . The aforementioned singularities/sources are preassigned and kept fixed throughout the solution process (i.e. the so-called static MFS approach is employed) on a pseudo-boundary Γ~ of a

Example 4. We consider two inclusions, i.e. J¼ 2, namely a circular inclusion of radius r 1 ¼ 0:25 and centered at xð1Þ

j

and an elliptical inclusion of semi-axes a ¼1.0 and b¼0.5, rotated by an angle φ ¼ π =4 about its centre xð2Þ     ( ) 2  X2 2 2 X1 2 Ω2 ¼ x A R  þ o r2 ; ð18bÞ  a b

similar shape to that of Γj, such that distðΓ~ j ; Γ j Þ, j A 0; J , are fixed constants, see e.g. [10]. Note that, for the examples considered in this paper, the number of inclusions is either J¼1 (Examples 1–3) or J¼2 (Example 4). To illustrate the positioning of the pseudo-boundaries used in ~ ¼ the present approach, we consider Example 1 and define ∂Ω 0 ~ ~ ~ Γ 0 [ Γ 01 and Γ 1 as follows:

where r 2 ¼ 0:25, whilst X 1 and X 2 are given by Eq. (16b).

Γ~ 0 ¼ fx A R2 jx21 þ x22 ¼ ½ð1 þ d0 ÞR2 g;

ð20aÞ

Γ~ 01 ¼ fx A R2 jðx1  x1ð1Þ Þ2 þ ðx2  x2ð1Þ Þ2 ¼ ½ð1  d1 Þr2 g;

ð20bÞ

Ω1 ¼ fx A R2 jðx1  x1ð1Þ Þ2 þ ðx2  x2ð1Þ Þ2 o r21 g;

ð18aÞ

In the present analysis, in the Neumann boundary conditions (6c) and (6d) the amplitude of the applied current was chosen as I 0 ðxÞ ¼  sin ðθðxÞ þ π =4Þ;

and

ð19Þ

ð1Þ 2 2 2 Γ~ 1 ¼ fx A R2 jðx1 xð1Þ 1 Þ þ ðx2  x2 Þ ¼ ½ð1 þ d2 Þr g;

while its frequency ω=2π was set to 500 kHz, and we only considered three different values for the time variable, namely t¼0 s, t ¼ 0:2  10  6 s and t ¼ 0:5  10  6 s, respectively.

where d0 4 0 and d1 ¼ d2 A ð0; 1Þ since Γ~ 01  Ω1 . A similar positioning of the pseudo-boundaries Γ~ j , j ¼ 0; J , is used in the case

1.0

1.0

8E 1

x2

+0 +0

1

E+

00

89

E+

x2

.9 7E

0 4.

00

E+ 00

37 E+

-0.5

01

0

99

98

0

+0

1.

01

+0

4E

1.

.0

96

00

E+ 1.

-1

-5 .5

1E

.0

8E

+0

0

+0

0

E+

E+

01

01

-0.5

.3

31

8.

-5 9E

-1

1.

0.0

00

E+

00

92

E+

36

.4

0

+0

E+

87

-1

+0

8E

1

32

4.

8. 1.

-1.0 -1.0

9E

.0

0.5

1.

0.0

-5

0

+0

-1

.9

+0

6E

1

.3

+0

-1

7E

0

4E

+0

.0

1E

-5

.9

.9

-8

-1

-8

0.5

-0.5

ð20cÞ

0.0

0.5

-1.0 -1.0

1.0

-0.5

x1

0.0

0.5

1.0

x1

1.0 -1

+0

7E

1

37

00

1.

E+

00

95

E+

8.

00

89

E+

4.

31

x2

0

+0

1.

0.0

E+ 01

1. 99

-0.5

1

0

8E

+0

.0

.3

+0

8E

-5

-1

8E

.9

.9

-8

0.5

-1

-5 .5

1E

.0

+0

8E

0

+0

0

E+ 01

-1.0 -1.0

-0.5

0.0

0.5

1.0

x1 Fig. 8. Distribution of the real part of the potential, v ¼ Re ϕ, at t ¼ 0:2  10  6 s and various numbers of collocation and source points, namely (a) M ¼N ¼ 80, (b) M ¼N ¼ 120 and (c) M ¼N ¼ 160, for Example 2.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

associated with the boundaries Γ0 and Γ1 (i.e. N 00 and N 01 ¼ N 1 ), for M¼80, M¼ 120 and M¼160, respectively. For all values of the number M of boundary collocation points considered, the following conclusions can be drawn from these figures:

where there are J inclusions present (i.e. d0 40 is constant associated with the outer boundary Γ0, and d2j  1 ¼ d2j associated with the interfaces Γj, j ¼ 1; J ). In order to assess the accuracy and convergence of the proposed MFS, for any real-valued function f : Γ j ⟶R, j ¼ 0; J , M and any set of points fxðj;mÞ gm j¼ 1  Γ j , we define the corresponding root mean square (RMS) error of f on Γj by errΓ j ðf Þ ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Mj ∑ jf ðxðj;mÞ Þj2 ; Mj m ¼ 1

j ¼ 0; J :

(i) For a fixed value of the number N 00 A f1; 2; …; M 00 g of source points corresponding to the outer boundary Γ0, the RMS error errΓ 1 ðv0  v1 Þ decreases as the number of singularities N 01 associated with the interface Γ1 increases. (ii) Although not illustrated in Figs. 2(a)–(c), it should be mentioned that, for fixed but large admissible values of d1 ¼ d2 (i.e. d1 ¼ d2 approaching 1.0) and a fixed value N00 A f1; 2; …; M 00 g, there exists a threshold value of the parameter N 01 , say N~ 01 , after which errΓ 1 ðv0  v1 Þ starts to increase as N 01 ⟶M 01 , the distribution of this RMS error becomes irregular and, consequently, the MFS approximations for the Dirichlet transmission condition on the interface Γ1 are more inaccurate. A possible reason for this increment in the inaccuracies of the MFS approximations for the real part of the potential transmission conditions on the interface Γ1 is the fact that increasing the number N00 of sources corresponding to Γ1 results in a highly ill-conditioned matrix A for the problem in Example 1. (iii) For a fixed value of the number N 01 A f1; 2; …; N~ 01 g of source points corresponding to the interface Γ1, the error in the Dirichlet transmission condition on Γ1, errΓ 1 ðv0  v1 Þ, has almost the same value for all N00 Z 5.

ð21Þ

We also investigate the influence of the number of source points N and of the distances from the boundaries of the solution domain, Γj, to the pseudo-boundaries on which the sources are located, Γ~ j , j ¼ 0; J , on the accuracy. Therefore, we analyse the dependence of the RMS error defined by Eq. (21) for the real and imaginary Dirichlet (potential) and Neumann (current) transmission conditions (6g)–(6j) on the interfaces Γj, j ¼ 1; J , with respect to the number of source points and the constants d0 and d2j  1 ¼ d2j , j ¼ 1; J , related to distðΓ~ j ; Γ j Þ, j A 0; J . We first consider the problem in Example 1 and set t¼0 s in Eq. (19), the number of collocation points M A f80; 120; 160g and distances d0 ¼ 2:0 and d1 ¼ d2 ¼ 0:3, at the same time varying the number of singularities N 0j A f1; 2; …; M 0j g, j ¼ 0; 1, i.e. such that N r M þ 2. Fig. 2(a), (b) and (c) presents the RMS error errΓ 1 ðv0  v1 Þ, which actually measures the accuracy of the potential transmission condition on the interface Γ1, as a function of the numbers of sources 1.0

1.0 -9

-6

0.5

.2

8E

.2

8E

+0

-9 +0

0

-4

0.5

.2

7E

.3

+0

6E

+0

0

0

81

01

E+ 00

E+

x2

-0.5

00

59

0

9.

00

-1

.0

5E

+0

0

00

00

E+

E+

53

+0

00

51

00

E+

1E

E+

14 6.

E+

00

.0

00

12 4.

01

E+

-1

E+

x2

0

.3

E-

2.

66

00

3.

E+

12 13 49 9.

-0.5

+0

-6

-2.

75

0.0

E-

2.

75

4. 6.

-1.0 -1.0

2E

0

-2. 3.

0.0

-0.5

9

0.0

0.5

-1.0 -1.0

1.0

-0.5

x1

0.0

0.5

1.0

x1

1.0 -9

-4

0.5

.2

0E

+0

-6 0

.3

5E

.3

+0

5E

+0

0

0

-2. E+ 00

E+

00

14

E+

01

12 4. 6. E+

00

51 58

00

9.

-0.5

E-

2.

71

x2

80 3.

0.0

-1

.0

4E

+0

0

E+ 00

-1.0 -1.0

-0.5

0.0

0.5

1.0

x1 Fig. 9. Distribution of the imaginary part of the potential, h ¼ Im ϕ, at t ¼ 0:2  10  6 s and various numbers of collocation and source points, namely (a) M ¼ N ¼80, (b) M ¼N ¼ 120 and (c) M ¼N ¼ 160, for Example 2.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

(iv) For all N 00 Z5, a fixed value of the number of source points corresponding to the interface Γ1 which is set to be proportional to the maximum admissible number of source points on Γ1, e.g. N 01 ¼ M 01 =2, the RMS error errΓ 1 ðv0  v1 Þ decreases with respect to the number of boundary collocation points, M, i.e. the MFS is convergent with respect to increasing M.

problems with exact boundary data, it is generally advised to place the singularities as far away from the boundary of the domain under investigation as possible, as much as the machine precision allows. Ramachandran [37] showed that the singular value decomposition (SVD) could mitigate this critical dependence. However, the direct problem investigated in this paper is somehow different from the situation presented above. More specifically, in the case of direct problems in EIT, the domain occupied by the background Ω0 is a multiply connected domain and, therefore, there exist source points associated with the interfaces Γj, j ¼ 1; J , and corresponding to the background which are located inside the simply connected domain Ω, more precisely, these sources are situated in the domains occupied by the inclusions Ωj  Ω, j ¼ 1; J . Consequently, the distances between the interfaces Γj and their corresponding pseudo-boundaries on which the MFS sources are located, and hence d2j  1 ¼ d2j , j ¼ 1; J , are constrained by the size of the inclusions Ωj  Ω0 , j ¼ 1; J . Nonetheless, it should be mentioned that there is no constraint on the distance between the outer boundary Γ0 and its associated pseudo-boundary Γ~ 0 and this can have any positive value. Fig. 4(a)–(c) present the RMS error errΓ 1 ðh0 h1 Þ as a function of the distances d0 A ½0:05; 3:0 and d1 ¼ d2 A ½0:01; 0:9 from the boundary Γj to the pseudo-boundary Γ~ j , j ¼ 0; 1, on which the sources are positioned, obtained using t ¼ 0:5  10  6 s and various numbers of collocation and source points, namely M¼ N ¼80, M¼N ¼ 120 and M¼ N ¼160, respectively, for Example 1. It can be

The same conclusions can be drawn from Fig. 3(a) to (c) for the RMS error corresponding to the transmission condition for the real part of the electric current on the interface Γ1, errΓ 1 ½nð0Þ  ðs0 ∇v0  ωε0 ∇h0 Þ þnð1Þ  ðs1 ∇v1  ωε1 ∇h1 Þ. The transmission conditions on Γ1 associated with the imaginary parts of the potential and current display similar behaviours, although these have not been presented here. As expected, errΓ 1 ðv0 v1 Þ o errΓ 1 ½nð0Þ  ðs0 ∇v0  ωε0 ∇h0 Þ þnð1Þ  ðs1 ∇v1  ωε1 ∇h1 Þ for the fixed values of the MFS parameters considered, namely t¼ 0 s in Eq. (19), the number of collocation points M A f80; 120; 160g and distances d0 ¼ 2:0 and d1 ¼ d2 ¼ 0:3, and all admissible values of the numbers of source points N 00 A f5; 6; …; M 00 g and N 01 A f1; 2; …; N~ 01 g, i.e. the MFS approximations for the Neumann (current) transmission conditions on the interface Γ1 are more inaccurate than the numerical results retrieved for the Dirichlet (potential) transmission conditions on Γ1. It is well-known that the accuracy of the MFS depends on the distance between the pseudo-boundary on which the singularities are located and the boundary of the solution domain. For direct

1.0

1.0 -2

x2

00

x2

E+

1

02 5.5

21

E+

E+

7E

-01

00

00

1.5

00

9E -0 2

-1.0 -1.0

1.0

18

0.5

E+

-4.

0.0

95

E+

0

-0.5

83

0

2. 3.

96

-0

+0

00

00

-0.5

1. 1.

00

00

4E

E-

0.0

1

03

E+

-0

E-

-1.0 -1.0

96

E+

2E

.01

-0.5

83

E+

E+

-3

2. 3.

95

5.3 20

.7

+0

1 1.

1.

59

5E

00

-0

0.0

.9

E+

0

0

7E

1.

8E

-7

0.5

.7

19

+0

.9

+0

-7

0E

7E

-4.

.9

.8

-1

-2

-1

0.5

00

-0.5

x1

0.0

0.5

1.0

x1

1.0 -2

E+ 00

x2

18

0

1

02

5.5 21

E+

E+

0E

-01

00

00

1.0

00

5E

E+

-4. -0 2

E+

83

+0

96

-0

0

2. 3.

96

5E

1.

1.

-1.0 -1.0

6E

+0

0.0

-0.5

.7

E-

8E

05

.7

1.

.9

-1 -7

0.5

00

-0.5

0.0

0.5

1.0

x1 Fig. 10. Distribution of the real part of the potential, v ¼ Re ϕ, at t ¼ 0:5  10  6 s and various numbers of collocation and source points, namely (a) M ¼ N¼ 80, (b) M ¼N ¼ 120 and (c) M ¼N ¼ 160, for Example 3.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

observed from these figures that, as the numbers of boundary collocation points M and singularities N increase, d0 increases and d1 ¼ d2 ⟶0:90, the distribution of the RMS error errΓ 1 ðh0  h1 Þ becomes irregular. This is a direct consequence of the fact that the corresponding MFS matrix A becomes highly ill-conditioned since not only the dimension of this matrix increases (i.e. M ¼N increase), but also the singularities on the pseudo-boundary Γ~

(ii) N0j ¼ M 0j for j ¼ 0; 1; 2, such that M ¼N for Examples 1–4; (iii) d0 ¼ 2:0 and d1 ¼ d2 ¼ 0:3 (Example 1); d0 ¼ 1:0 and d1 ¼ d2 ¼ 0:2 (Example 2); d0 ¼ 2:0 and d1 ¼ d2 ¼ 0:2 (Example 3); d0 ¼ 1:0, d1 ¼ d2 ¼ 0:3 and d3 ¼ d4 ¼ 0:3 (Example 4). For Examples 1–4, d1 and d2 set the positions of the pseudoboundaries corresponding to Γ1, whilst the distances between Γ2 and its corresponding pseudo-boundaries in Example 4 are determined by d3 and d4. Figs. 6 and 7 present the contour plots of the real and imaginary parts of the electric potential ϕ, v and h, respectively, obtained in the domain corresponding to Example 1 by applying the boundary current given by Eq. (19) with M¼N¼160 and various values of the time variable, namely t¼0 s, t ¼ 0:2  10  6 s and t ¼ 0:5  10  6 s. It can be seen from Figs. 6(a) and (b) that although the numerically reconstructed values of v ¼ Re ϕ at t¼0 s and t ¼ 0:2  10  6 s are different, the distributions of the real electric potential in Ω are very similar. However, the real electric potential numerically retrieved at t ¼ 0:5  10  6 s and illustrated in Fig. 6(c) clearly accounts for the presence of the circular inclusion Ω1 given by Eq. (15). As expected, the behaviour of the reconstructed imaginary part of the electric potential h ¼ Im ϕ is opposite to that of the real electric potential v ¼ Re ϕ, see Figs. 6 and 7. More precisely, at t¼0 s and t ¼ 0:2  10  6 s, although their numerical values are different, both imaginary electric potential distributions account for the presence of the circular inclusion Ω1, see Fig. 7(a) and (b). On the contrary, at t ¼ 0:5  10  6 s the presence of the inclusion Ω1

1

associated with the interface Γ1 cluster around the origin. Similar results have been obtained for the RMS error defined by Eq. (21) and corresponding to the imaginary part of the Neumann (current) transmission condition on Γ1 and these are displayed in Fig. 5(a)–(c). Although not illustrated, similar results were also obtained, in the case of Example 1, for the RMS errors errΓ 1 ðv0  v1 Þ and errΓ 1 ½nð0Þ  ðs0 ∇v0  ωε0 ∇h0 Þ þ nð1Þ  ðs1 ∇v1  ωε1 ∇h1 Þ, as well as the other examples analysed in this study. According to the RMS errors associated with the transmission conditions presented in Figs. 4 and 5, for stability and accuracy reasons it is recommended to choose the source points corresponding to the interfaces (inner boundaries) Γj, j ¼ 1; J , on the pseudo-boundaries Γ~ j , j ¼ 1; J , such that d2j  1 ¼ d2j A ½0:02; 0:35, j ¼ 1; J . The numerical results presented in the sequel have been obtained by setting the corresponding MFS parameters introduced in Section 3 as follows: (i) M 00 A f40; 60; 80g and M 0j ¼ M j A f20; 30; 40g for j ¼ 1; 2, such that M A f80; 120; 160g for Examples 1–3 and M A f120; 180; 240g for Example 4;

1.0

1.0 -2

0

x2

x2

+0

+0

1

1

00

E+

00

00 01

17

-1

.8

6E

+0

0

E+

01

01

E+

-0.5

E+

-0.5

E+

49

0

2.

+0

+0

5E

E+

17 1.

00

16

01

2.

E+

5E

8E

0

5. 64

00

E+

48

.8

.4

.1

24

9.

E+

60

-1

-1

1.

0.0

00

15

E+

5.

25

9. 1.

-1.0 -1.0

3E +0

1

1.

0.0

-2

5E

0

+0

0.5

1

.7

6E

+0

-5

.4

3E

.8

+0

0

7E

+0

.6

8E

-5

.6

-1

.1

-9

-9

0.5

-0.5

11

0.0

0.5

-1.0 -1.0

1.0

-0.5

x1

0.0

0.5

1.0

x1

1.0

-9

-2

+0

0

4E

+0

.7

9E

-5

.7

0.5

-1

.4

8E

+0

1

1

E+

9.

00

17

E+

5.

25

x2

0

+0

4E

1.

0.0

E+

00

64 49

00

1. E+

17

01

2.

-0.5

.1

-1

.8

5E

+0

0

E+ 01

-1.0 -1.0

-0.5

0.0

0.5

1.0

x1 Fig. 11. Distribution of the imaginary part of the potential, h ¼ Im ϕ, at t ¼ 0:5  10  6 s and various numbers of collocation and source points, namely (a) M ¼ N ¼80, (b) M ¼N ¼ 120 and (c) M ¼N ¼ 160, for Example 3.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

5. Conclusions

has little effect on the numerically reconstructed values of the imaginary part of the electric potential h ¼ Im ϕ presented in Fig. 7(c). The approach developed here was also shown to be efficient, accurate and convergent with respect to increasing the number of boundary collocation points when the inclusion Ω1 was assumed to be smooth, but non-convex. Two examples have been considered to illustrate this situation, namely the peanut-shaped inclusion considered in Example 2 for t ¼ 0:2  10  6 s and the acornshaped inclusion in Example 3 for t ¼ 0:5  10  6 s. The contour plots of the numerically retrieved values for v and h in Ω and corresponding to Examples 2 and 3 are illustrated in Figs. 8–11. The proposed procedure works equally well when two inclusions (J¼2) are present in the background medium as is the case in Example 4, in which we consider a circular inclusion and an elliptical one given by Eqs. (18a) and (18b), respectively. Figs. 12 and 13 present the contour plots of the retrieved real, v, and imaginary parts of the potential, h, respectively, using the boundary current given by Eq. (19) at t ¼ 0:5  10  6 s and various numbers of collocation and source points, i.e. M ¼ N A f80; 120; 160g. Overall, from Figs. 6–13 we can conclude that the MFS provides accurate and convergent numerical solutions for the real and imaginary parts of the electric potential in the background medium Ω0 as well as in the inclusions Ωj  Ω, j ¼ 1; J , at all time instants t Z 0.

In this paper, the MFS has been implemented and analysed for obtaining the numerical solution of the complex EIT forward problem in two dimensions. A generalized MFS formulation for this EIT direct problem has been developed for the case when a finite number of inclusions with constant physical properties (i.e. conductivity and permittivity) are assumed to be compactly contained in a given medium (background). The efficiency, accuracy and numerical convergence of the proposed method have been investigated and validated by considering four numerical examples corresponding to various smooth convex or non-convex shapes of the inclusion(s), namely circular, elliptic, peanut-shaped and acorn-shaped. The proposed technique is efficient and can be extended to other problems for complex EIT: three-dimensional forward problems, inverse geometric problems (i.e. reconstruction of the inclusion(s) in the background when the physical properties of both the inclusion (s) and the background are known a priori), parameter identification problems (i.e. identification of the physical properties of the inclusion (s) in the background) and combined geometric-parameter identification problems in two and three dimensions. However, these extensions will be considered in future work.

1.0

1.0 -3 -4.

x2

x2

1.

49

00

38

00

4.

E+

-01

00

26

-01

0E

E+

3.

9E

37

E+

7.0

2.

6.8 E+ 00

00

-0.5

00

0.0

-0.5

E+

-1.0 -1.0

E+

0

0

00

37

00

4.

E+

E+

24

00

0

+0

+0

36

E+

+0

3E

2E

2. 3.

-0.5

47

39E

.1

.1

-2

00

-2

0.0 1.

65

-3.

0

0.5

E+

+0

64

8E

-4.

.3

0.5

0.0

0.5

-1.0 -1.0

1.0

-0.5

x1

0.0

0.5

1.0

x1

1.0

-4.

-3

x2

00

0

E+

3.

-01

38

00

4.

E+

00

26

0E

38

00

+0

7.0

2.

E+

3E

0

49

+0

1.

.1

E+

9E -2

0.0

E+

-0.5

65

.3

0.5

00

-1.0 -1.0

-0.5

0.0

0.5

1.0

x1 Fig. 12. Distribution of the real part of the potential, v ¼ Re ϕ, at t ¼ 0:5  10  6 s and various numbers of collocation and source points, namely (a) M ¼N ¼ 120, (b) M ¼N ¼ 180 and (c) M ¼ N¼ 240, for Example 4.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1.0

1.0 -2

-1

0.5

+0

1

1

+0 0

x2

+0

7E

0 1. E+

00

00

73

E+

E+

4.

02

00

01 01

01

E+

E+

10

10

2.

01

E+

00

42

E+

00

1.

E+

00

E+

-0.5

15

9.

12

75

E+

4.

05 43

-0.5

0E

.0

+0

0

0.0

1.

9. 1. 2.

-1.0 -1.0

.4

3E

0

+0

x2

1

8E

1E

0.0

-1

0.5

1

.3

.4

+0

+0

-2 +0

-5

-5

2E

0E

7E

.3

.3

.4

.0

-9

-9

-0.5

13

0.0

0.5

-1.0 -1.0

1.0

-0.5

0.0

x1

0.5

1.0

x1

1.0 -2

-1

0.5

.4

-9

0E

.3 3E

-5

+0

1

1

0

8E +0

0

x2

+0

7E

+0

.3

0.0 1.

E+

E+

00

E+

00

42

00

72

E+

4.

01

15

9. 1. 10

01

2.

-0.5

.0

E+ 01

-1.0 -1.0

-0.5

0.0

0.5

1.0

x1 Fig. 13. Distribution of the imaginary part of the potential, h ¼ Im ϕ, at t ¼ 0:5  10  6 s and various numbers of collocation and source points, namely (a) M ¼N ¼ 120, (b) M¼ N ¼ 180 and (c) M ¼N ¼ 240, for Example 4.

Acknowledgments L. Marin acknowledges the financial support received from the Romanian National Authority for Scientific Research (CNCSUEFISCDI), Project number PN-II-ID-PCE-2011-3-0521. The authors also thank the reviewers for their constructive comments and suggestions, which significantly contributed to improving the quality of the publication. References [1] Borcea L. Electrical impedance tomography. Inverse Probl 2002;18:R99–136. [2] Cheney M, Isaacson D, Newell JC, Simske S, Goble J. NOSER: an algorithm for solving the inverse conductivity problem. Int J Imaging Syst Technol 1990;2:66–75. [3] Berger JA, Karageorghis A. The method of fundamental solutions for heat conduction in layered materials. Int J Numer Methods Eng 1999;45:1681–94. [4] Berger JA, Karageorghis A. The method of fundamental solutions for layered elastic materials. Eng Anal Bound Elem 2001;25:877–86. [5] Fairweather G, Karageorghis A. The method of fundamental solutions for elliptic boundary value problems. Adv Comput Math 1998;9:69–95. [6] Fairweather G, Karageorghis A, Martin PA. The method of fundamental solutions for scattering and radiation problems. Eng Anal Bound Elem 2003;27:759–69. [7] Folland GB. Introduction to partial differential equations. Princeton, NJ: Princeton University Press; 1995. [8] Gabriel S, Lau RW, Gabriel C. The dielectric properties of biological tissues: II. Measurements in the frequency range 10 Hz to 20 GHz. Phys Med Biol 1996;41:2251–69. [9] Golberg MA, Chen CS. The method of fundamental solutions for potential, Helmholtz and diffusion problems. In: Golberg MA, editor. Boundary integral

[10]

[11]

[12] [13]

[14]

[15] [16] [17] [18]

[19]

[20] [21]

methods: numerical and mathematical aspects. Computational engineering, vol. 1. Boston, MA: WIT Press/Computational Mechanics Publications; 1999. p. 103–76. Gorzelańczyk P, Kołodziej JA. Some remarks concerning the shape of the shape contour with application of the method of fundamental solutions to elastic torsion of prismatic rods. Eng Anal Bound Elem 2008;32:64–75. Harrach B, Seo JK, Woo EJ. Factorization method and its physical justification in frequency-difference electrical impedance tomography. IEEE Trans Med Imaging 2010;29:1918–26. Holder D. Electrical impedance tomography: methods, history and applications. Bristol, UK: Institute of Physics Publishing; 2005. Jun SC, Kuen J, Lee J, Woo EJ, Holder D, Seo JK. Frequency-difference EIT (fdEIT) using weighted difference and equivalent homogeneous admittivity: validation by simulation and tank experiment. Physiol Meas 2009;30:1087–99. Karageorghis A. Modified methods of fundamental solutions for harmonic and biharmonic problems with boundary singularities. Numer Methods Partial Differ Equ 1992;8:1–19. Karageorghis A. The method of fundamental solutions for the calculation of the eigenvalues of the Helmholtz equation. Appl Math Lett 2001;14:837–42. Karageorghis A, Fairweather G. The method of fundamental solutions for the numerical solution of the biharmonic equation. J Comput Phys 1987;69:434–59. Karageorghis A, Fairweather G. The method of fundamental solutions for axisymmetric elasticity problems. Comput Mech 2000;25:524–32. Karageorghis A, Lesnic D. Steady-state nonlinear heat conduction in composite materials using the method of fundamental solutions. Comput Methods Appl Mech Eng 2008;197:3122–37. Karageorghis A, Smyrlis Y-S. Matrix decomposition MFS algorithms for elasticity and thermo-elasticity problems in axisymmetric domains. J Comput Appl Math 2007;206:774–95. Karageorghis A, Lesnic D, Marin L. A survey of applications of the MFS to inverse problems. Inverse Probl Sci Eng 2011;19:309–36. Karageorghis A, Poullikkas A, Berger JR. Stress intensity factor computation using the method of fundamental solutions. Comput Mech 2006;37:445–54.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i

14

M.A. Heravi et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎

[22] Kim S, Lee EJ, Woo EJ, Seo JK. Asymptotic analysis of the membrane structure to sensitivity of frequency-difference electrical impedance tomography. Inverse Probl 2012;28:17. [23] Kupradze VD, Aleksidze MA. The method of functional equations for the approximate solution of certain boundary value problems. USSR Comput Math Math Phys 1964;4:82–126. [24] Lionheart W, Polydorides N, Borsic A. The reconstruction problem. In: Holder DS, editor. Electrical impedance tomography: methods, history and applications. Bristol, UK: Institute of Physics Publishing; 2005. [25] Marin L. Stable MFS solution to singular direct and inverse problems associated with the Laplace equation subjected to noisy data. Comput Model Eng Sci 2008;37:203–42. [26] Marin L. Treatment of singularities in the method of fundamental solutions for two-dimensional Helmholtz-type equations. Appl Math Model 2010;34: 1615–33. [27] Marin L. A meshless method for the stable solution of singular inverse problems for two-dimensional Helmholtz-type equations. Eng Anal Bound Elem 2010;34:274–88. [28] Marin L, Karageorghis A. The MFS-MPS for two-dimensional steady-state thermoelasticity problems. Eng Anal Bound Elem 2013;37:1004–20. [29] Marin L, Lesnic D. The method of fundamental solutions for nonlinear functionally graded materials. Int J Solids Struct 2007;44:6878–90. [30] Mathon R, Johnston RL. The approximate solution of elliptic boundary value problems by fundamental solutions. SIAM J Numer Anal 1977;14:638–50. [31] Metherall P, Barber DC, Smallwood RH, Brown BH. Three-dimensional electrical impedance tomography. Nature 1996;380:509–12.

[32] Oh TI, Koo W, Lee KH, Kim SM, Lee J, Kim SW, Seo JK, Woo EJ. Validation of multi-frequency electrical impedance tomography (mfEIT) KHU Mark 1: impedance spectroscopy and time-difference imaging. Physiol Meas 2008;28: 1175–88. [33] Poullikkas A, Karageorghis A, Georgiou G. Methods of fundamental solutions for harmonic and biharmonic boundary value problems. Comput Mech 1998;21:416–23. [34] Poullikkas A, Karageorghis A, Georgiou G. The method of fundamental solutions for inhomogeneous elliptic problems. Comput Mech 1998;22:100–7. [35] Poullikkas A, Karageorghis A, Georgiou G. The numerical solution of threedimensional Signorini problems with the method of fundamental solutions. Eng Anal Bound Elem 2001;25:221–7. [36] Poullikkas A, Karageorghis A, Georgiou G. The numerical solution for threedimensional elastostatics problems. Comput Struct 2002;80:365–70. [37] Ramachandran PA. Method of fundamental solutions: singular value decomposition analysis. Commun Numer Methods Eng 2002;18:789–801. [38] Seo JK, Lee J, Kim SW, Zribi H, Woo EJ. Frequency-difference electrical impedance tomography (fdEIT): algorithm development and feasibility study. Physiol Meas 2008;29:929–44. [39] Tsai CC. The method of fundamental solutions with dual reciprocity for threedimensional thermoelasticity under arbitrary forces. Int J Comput-Aided Eng Softw 2009;26:229–44. [40] Woo EJ, Hua P, Webster JG, Tompkins WJ. A robust image reconstruction algorithm and its parallel implementation in electrical impedance tomography. IEEE Trans Med Imaging 1993;12:137–46. [41] Yorkey TJ, Webster JG, Tompkins WJ. Comparing reconstruction algorithms for electrical impedance tomography. IEEE Trans Biomed Eng 1987;34:843–52.

Please cite this article as: Heravi MA, et al. The method of fundamental solutions for complex electrical impedance tomography. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.04.022i