Convergence of finite element approximation for electrical impedance tomography with the complete electrode model

The complete electrode model (Cheng et al 1989) is a standard model for electrical impedance tomography (EIT), when the conductivity distribution over a region is of interest and direct current is applied (Cheney et al 1999, Holder 2005). The forward problem of EIT with this model involves a boundary value problem (BVP) that has been proved to have an unique solution in an appropriate space (Somersalo et al 1992).

In practice, the forward problem is often solved by the finite element method (FEM) (Vauhkonen et al 1999, Polydorides and Lionheart 2002). The finite element approximation is desired to be as accurate as possible for some iterative reconstruction algorithms (Yorkey et al 1987) of EIT to perform effectively (Paulson et al 1992). High accuracy in finite element approximation is usually attempted to achieve by refining the underlying meshes.

It is worth studying whether the finite element approximation will converge to the true solution of BVP in EIT with the complete electrode model. First, the BVP for the complete electrode model has nontrivial boundary conditions, and is of different kind from the ones about which the FEM has been proved to converge in standard texts e.g. (Braess 2007, Brenner and Scott 2008). Second, the current density in EIT with the complete electrode model exhibited some kind of singularity around the edge of the electrodes (Pidcock et al 1995, Ciulli and Ispas 1996). This may complicate the convergence of finite element approximation for this BVP.

Some works (Gehre et al 2014, Dardé and Staboulis 2016, Jin et al 2017) involve the convergence analysis of finite element approximation in EIT. The (Gehre et al 2014, Jin et al 2017) are mainly focused on inverse problems in EIT. However, the convergence analysis cover also the case under consideration in this paper. The (Dardé and Staboulis 2016) also gives the convergence rate of finite element approximation.

We discuss the convergence of finite element approximation for the forward problem of EIT with the complete electrode model. Here the region in EIT is assumed to be a convex polygonal domain. Conforming finite element method with linear triangle or tetrahedron elements is employed. The resulting finite element solution is to be proved to converge to the true one in some norm as the sizes of elements in the underlying mesh all tend to zero.

The outline of this paper is as follow. Section 2 starts by reviewing the BVP for EIT with the complete electrode model, the weak solution of this BVP, the existence and the uniqueness of the weak solution. Then we discuss the finite element approximation for the solution of this BVP and its convergence which is the main result in this paper. In section 3, convergence of the finite element approximation were demonstrated with simulations on a toy 2-D EIT model. The conclusion of this paper is made in the last section.

2.1. Complete electrode model

Assume the conductivity distribution σ over a region Ω is of interest. On the boundary ∂Ω of this region, a number of L electrodes are mounted each of which is assumed to be a perfect conductor. The part of boundary that is covered by the i-th electrode is denoted by e_i, i = 1, 2, ..., L. Through each electrode, direct current is either injected into or extracted out of the region, with an amount denoted by I_i for the i-th electrode, i = 1, 2, ..., L. It is specified that the amount of current takes a positive value when the current is injected into the region through its corresponding electrode, and a negative value when being extracted out of the region. The resulting voltage on each electrode is measured, and is denoted by U_i for the i-th electrode, i = 1, 2, ..., L. Contact impedance is assumed to exist between each electrode and the region, and its value is z_i for the i-th electrode. Let u be the potential distribution over the region. All these quantities are related as follow:

$$\begin{eqnarray}&&-{\rm{\nabla }}\cdot (\sigma {\rm{\nabla }}u)=0\qquad \qquad u\in {\rm{\Omega }},\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&-\sigma \displaystyle \frac{\partial u}{\partial n}=0\ \qquad \mathrm{on}\ \partial {\rm{\Omega }}\setminus \,{\cup }_{i=1}^{L}{e}_{i},\end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&{\int }_{{e}_{i}}\sigma \displaystyle \frac{\partial u}{\partial n}{\rm{d}}s={I}_{i}\ \qquad \mathrm{for}\ i=1,2,\,\ldots ,\,L,\end{eqnarray} \tag{ 1c }$$

$$\begin{eqnarray}&&u+{z}_{i}\sigma \displaystyle \frac{\partial u}{\partial n}={U}_{i}\quad \qquad \mathrm{on}\ {e}_{i}\ \mathrm{for}\ i=1,2,\,\ldots ,\,L,\end{eqnarray} \tag{ 1d }$$

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{L}{I}_{i}=0.\end{eqnarray} \tag{ 1e }$$

The (1a) is because no source of current exists within the region. The (1b) is because currents go through the region only through the electrodes. The (1d) is due to that a thin layer of impedance exists beneath each electrode. The (1e) represents current conservation. The model (1a) is called the complete electrode model of EIT (Cheng et al 1989). It can predict experimentally measured voltages to within 0.1% (Somersalo et al 1992).

The forward problem of EIT is to solve BVP (1a) for u and U₁, U₂, ..., U_L, given that all the other quantities in this BVP are known. In practice, the BVP is often solved numerically by the finite element method (FEM) (Vauhkonen 1997, Vauhkonen et al 1999, Polydorides and Lionheart 2002). Because it can handle easily boundary conditions (1b), (1c), (1d), even for a region with complex geometry.

2.2. Weak Solution

The finite element approximation for the numerical solution of BVP (1a) could be based on its variational formulation (Somersalo et al 1992). Assume Ω is a open bounded domain in R² or R³ with Lipsize-continuous boundary which is defined in (Ciarlet 2002). The electrodes e_l, l = 1, 2, ..., L are open connected subsets of ∂Ω whose closures are disjoint. When n = 3, the boundary of the electrodes are in addition assumed to be smooth curves on ∂Ω. The conductivity σ is assume to be real-valued here. It is continuously differentiable in Ω up to the boundary as defined in (Somersalo et al 1992).

Let

$$\begin{eqnarray}&&H\mathop{=}\limits^{\mathrm{def}}{H}^{1}({\rm{\Omega }})\otimes {R}^{L}\end{eqnarray} \tag{ 2 }$$

i.e. the Cartesian product of L²-based Sobolev space H¹(Ω) and L-dimensional real Euclidean space R^L. For $(u,U),(v,V)\in H$ with $u,v\in {H}^{1}({\rm{\Omega }})$ and $U={({U}_{1},{U}_{2},\ldots ,{U}_{L})}^{T},V={({V}_{1},{V}_{2},\ldots ,{V}_{L})}^{T}\in {R}^{L}$,

$$\begin{eqnarray}&&a((u,U),(v,V))\mathop{=}\limits^{\mathrm{def}}{\int }_{{\rm{\Omega }}}\sigma {\rm{\nabla }}u\cdot {\rm{\nabla }}v{\rm{d}}{\rm{\Omega }}+\displaystyle \sum _{i=1}^{L}\displaystyle \frac{1}{{z}_{i}}{\int }_{{e}_{i}}(u-{U}_{i})(v-{V}_{i}){\rm{d}}s,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\langle {\ell },(v,V)\rangle \mathop{=}\limits^{\mathrm{def}}\displaystyle \sum _{i=1}^{L}{V}_{i}{I}_{i}.\end{eqnarray} \tag{ 4 }$$

Then $(u,U)\in H$ is called a weak solution of BVP (1a) if

$$\begin{eqnarray}&&a((u,U),(v,V))=\langle {\ell },(v,V)\rangle \quad \mathrm{for}\ \mathrm{all}\ (v,V)\in H.\end{eqnarray} \tag{ 5 }$$

Every weak solution that lies in ${C}^{2}({\rm{\Omega }})\cap C(\bar{{\rm{\Omega }}})$ is a classical solution of this BVP, and every classical solution is also a weak solution (Somersalo et al 1992). In H, the weak solution of BVP (1a) if it exists is not unique. This is because if $(u,U)\in H$ is a weak solution, so is (u, U) + c, for any $c\in R$, according to (3), (4) and (1e).

2.3. Existence and uniqueness

If σ and z₁, z₂, ..., z_L are positive-valued and σ is bounded from above, then the weak solution of BVP (1a) exists and is unique in the quotient space of H over R (Somersalo et al 1992). This quotient space is defined as

$$\begin{eqnarray}&&H/R\mathop{=}\limits^{\mathrm{def}}\{[(u,U)]:(u,U)\in H\}\end{eqnarray} \tag{ 6 }$$

with

$$\begin{eqnarray}&&[(u,U)]=\{(u,U)+c:c\in R\},\end{eqnarray} \tag{ 7 }$$

where $u\in {H}^{1}({\rm{\Omega }})$ and $U\in {R}^{L}$. By saying BVP (1a) has an unique weak solution in H/R, we mean there exists only one element $[(u,U)]\in H/R$ such that

$$\begin{eqnarray}&&a([(u,U)],[(v,V)])=\langle {\ell },[(v,V)]\rangle \quad \mathrm{for}\ \mathrm{all}\ [(v,V)]\in H/R,\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&a([(u,U)],[(v,V)])\mathop{=}\limits^{\mathrm{def}}{\int }_{{\rm{\Omega }}}\sigma {\rm{\nabla }}u\cdot {\rm{\nabla }}v{\rm{d}}{\rm{\Omega }}+\displaystyle \sum _{i=1}^{L}\displaystyle \frac{1}{{z}_{i}}{\int }_{{e}_{i}}(u-U)(v-V){\rm{d}}s,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\langle {\ell },[(v,V)]\rangle \mathop{=}\limits^{\mathrm{def}}\displaystyle \sum _{i=1}^{L}{V}_{i}{I}_{i}.\end{eqnarray} \tag{ 10 }$$

In other words, the BVP (1a) has weak solutions in H and all of its weak solutions are in the form $\{(u,U)+c,c\in R\}$ where (u,U) is one of its weak solutions.

The existence and uniqueness of the week solution in H/R could be proved as follow. First, the bilinear form (9) is elliptic with respect to H/R (Somersalo et al 1992). That is to say, there exists α > 0 such that

$$\begin{eqnarray}&&a([u,U],[(u,U)])\geqslant \alpha \parallel [(u,U)]{\parallel }_{H/R}\qquad \mathrm{for}\ \mathrm{all}\ [(u,U)]\in H/R,\end{eqnarray} \tag{ 11 }$$

where $\parallel {\parallel }_{H/R}$ is the natural quotient norm for the space H/R. Since bilinear form (9) is elliptic (11) and the linear form (10) is bounded, the term

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}a([(u,U)],[(u,U)])-\langle {\ell },[(u,U)]\rangle \end{eqnarray} \tag{ 12 }$$

has a unique minimizer in H/R according to the Lax-Milgram theorem in Braess (2007, Page 38). Because the bilinear form (9) is also symmetric and positive, the characterization theorem in Braess (2007, Page 35) ensures that the unique minimizer of (12) is also the unique one satisfying (8). So the weak solution of BVP (1a) exists and is unique in H/R.

It is expected that the weak solution is unique in the quotient space H/R. This is because the BVP (1a) determines the potential distribution and voltages on the electrodes only up to an additive constant, and reference has to be specified when voltages are measured.

2.4. Finite element approximation and its convergence

The approximations for the numerical solutions of BVP (1a) by the finite element method are the minimizer of (12) over finite element spaces. The finite element spaces are associated with a triangulation ${{ \mathcal T }}_{h}$ of the domain Ω, where the subscript h represents the size of largest element in the mesh.

As regard to the convergence of finite element approximations, we prove the following theorem.

Theorem. Assume the domain Ω in BVP (1a) is polygonal, and is discretized by a family of shape-regular triangulation ${{ \mathcal T }}_{h}$. Linear triangle or tetrahedron elements are employed. All the finite elements associated with this family of triangulations are affine-equivalent to a single reference finite element as defined in Ciarlet (2002). The resulting conforming finite element method produces estimated voltages ${U}_{h}\in {R}^{L}$ on all the L electrodes. $U\in {R}^{L}$ are true voltages on the electrodes. Then

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{h\to 0}\mathop{\inf }\limits_{c\in R}\parallel U-{U}_{h}+c{\parallel }_{{R}^{L}}=0.\end{eqnarray} \tag{ 13 }$$

Proof. Let $u\in {H}^{1}({\rm{\Omega }})$ be a true potential distribution over the region. The space where a finite element approximation u_h for u lies is denoted by ${H}_{h}^{1}({\rm{\Omega }})$. Because conforming finite element method is employ,

$$\begin{eqnarray}&&{H}_{h}^{1}({\rm{\Omega }})\subset {H}^{1}({\rm{\Omega }}).\end{eqnarray} \tag{ 14 }$$

A finite element solution (u_h, U_h) lies in

$$\begin{eqnarray}&&{H}_{h}\mathop{=}\limits^{\mathrm{def}}{H}_{h}^{1}({\rm{\Omega }})\otimes {R}^{L}.\end{eqnarray} \tag{ 15 }$$

i.e. the Cartesian product of ${H}_{h}^{1}({\rm{\Omega }})$ and R^L. The quotient space of H_h over R is

$$\begin{eqnarray}&&{H}_{h}/R\mathop{=}\limits^{\mathrm{def}}\{[({u}_{h},{U}_{h})]:({u}_{h},{U}_{h})\in {H}_{h}\}\end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray}&&[({u}_{h},{U}_{h})]\mathop{=}\limits^{\mathrm{def}}\{({u}_{h},{U}_{h})+c:c\in R\}.\end{eqnarray} \tag{ 17 }$$

Note that

$$\begin{eqnarray}&&{H}_{h}/R\subset H/R.\end{eqnarray} \tag{ 18 }$$

according to the definitions of these spaces.

First, we bound the difference between the finite element solutions [(u_h, U_h)] and the true solution [(u, U)] that is measured in the quotient norm $\parallel \cdot {\parallel }_{H/R}$:

$$\begin{eqnarray}&&{\parallel [(u,U)]-[({u}_{h},{U}_{h})]\parallel }_{H/R}\leqslant C\mathop{\inf }\limits_{[({v}_{h},{V}_{h})]\in {H}_{h}/R}{\parallel [(u,U)]-[({v}_{h},{V}_{h})]\parallel }_{H/R}\end{eqnarray} \tag{ 19 }$$

with some positive constant C. The quotient norm is defined as follow

$$\begin{eqnarray}&&{\parallel [(u,U)]\parallel }_{H/R}\mathop{=}\limits^{\mathrm{def}}\mathop{\inf }\limits_{c\in R}(\parallel u-c{\parallel }_{{H}^{1}{\rm{\Omega }}}+\parallel U-c{\parallel }_{{R}^{L}})\end{eqnarray} \tag{ 20 }$$

The (19) is proved as follow. A finite element solution (u_h, U_h) is a minimizer of (12) in H_h, and is characterized by

$$\begin{eqnarray}&&a(({u}_{h},{U}_{h}),({v}_{h},{V}_{h}))=\langle {\ell },({v}_{h},{V}_{h})\rangle \quad \mathrm{for}\ \mathrm{all}\ ({v}_{h},{V}_{h})\in {H}_{h}.\end{eqnarray} \tag{ 21 }$$

So [(u_h, U_h)] satisfies

$$\begin{eqnarray}&&a([({u}_{h},{U}_{h})],[({v}_{h},{V}_{h})])=\langle {\ell },[({v}_{h},{V}_{h})]\rangle \quad \mathrm{for}\ \mathrm{all}\ [({v}_{h},{V}_{h})]\in {H}_{h}/R,\end{eqnarray} \tag{ 22 }$$

because of (18), (9) and (10). The (22), together with (8) and (18), gives

$$\begin{eqnarray}&&a([(u-{u}_{h},U-{U}_{h})],[(v,V)])=0\quad \mathrm{for}\ \mathrm{all}\ [(v,V)]\in {H}_{h}/R.\end{eqnarray} \tag{ 23 }$$

From (23) and that the bilinear form (9) is H/R-elliptic, we have

$$\begin{eqnarray}&&{C}_{1}\parallel [(u,U)]-[({u}_{h},{U}_{h})]{\parallel }_{H/R}^{2}\leqslant a([(u,U)]-[({u}_{h},{U}_{h})],[(u,U)]-[({u}_{h},{U}_{h})])\\ &&\quad =a([(u,U)]-[({u}_{h},{U}_{h})],[(u,U)]-[({v}_{h},{V}_{h})])\\ &&\quad +a([(u,U)]-[({u}_{h},{U}_{h})],[({v}_{h},{V}_{h})]-[({u}_{h},{U}_{h})])\\ &&\quad \leqslant {C}_{2}{\parallel [(u,U)]-[({u}_{h},{U}_{h})]{\parallel }_{H/R}\parallel [(u,U)]-[({v}_{h},{V}_{h})]\parallel }_{H/R}\\ &&\quad \mathrm{for}\ \mathrm{any}[({v}_{h},{V}_{h})]\in {H}_{h}/R\end{eqnarray} \tag{ 24 }$$

with some positive constant C₁ and C₂. From (24),

$$\begin{eqnarray}&&{\parallel [(u,U)]-[({u}_{h},{U}_{h})]\parallel }_{H/R}\leqslant \displaystyle \frac{{C}_{2}}{{C}_{1}}{\parallel [(u,U)]-[({v}_{h},{V}_{h})]\parallel }_{H/R}\ \mathrm{for}\ \mathrm{any}[({v}_{h},{V}_{h})]\in {H}_{h}/R.\end{eqnarray} \tag{ 25 }$$

The (25) implies (19).

Second, we use a density argument as in the proof of theorem 3.2.3 in (Ciarlet 2002) to show

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{h\to 0}\mathop{\inf }\limits_{[({v}_{h},{V}_{h})]\in {H}_{h}/R}\parallel [(u,U)]-[({v}_{h},{V}_{h})]{\parallel }_{H/R}=0.\end{eqnarray} \tag{ 26 }$$

We start by introducing space

$$\begin{eqnarray}&&{ \mathcal V }\mathop{=}\limits^{\mathrm{def}}{W}^{2,\infty }({\rm{\Omega }})\cap {H}^{1}({\rm{\Omega }}),\end{eqnarray} \tag{ 27 }$$

where ${W}^{2,\infty }({\rm{\Omega }})$ is a Sobolev space based on ${L}^{\infty }({\rm{\Omega }})$. Since ${ \mathcal V }$ is dense in H¹(Ω), given  > 0, there exists a ${v}^{* }\in { \mathcal V }$ such that

$$\begin{eqnarray}&&\parallel u-{v}^{* }{\parallel }_{1}\lt \epsilon /2.\end{eqnarray} \tag{ 28 }$$

For v^*, we claim

$$\begin{eqnarray}&&\parallel {v}^{* }-{{\rm{\Pi }}}_{h}{v}^{* }{\parallel }_{1,{\rm{\Omega }}}\leqslant {C}_{3}h{(\mathrm{meas}({\rm{\Omega }}))}^{1/2}| {v}^{* }{| }_{2,\infty ,{\rm{\Omega }}}.\end{eqnarray} \tag{ 29 }$$

Here ${{\rm{\Pi }}}_{h}{v}^{* }$ is the unique interpolant of v^* associated with finite element space ${H}_{h}^{1}({\rm{\Omega }})$. meas(Ω) is Lebesgue measure of domain Ω. $\parallel \cdot {\parallel }_{1,{\rm{\Omega }}}$ is the natural norm of space H¹(Ω), and $| \cdot {| }_{2,\infty ,{\rm{\Omega }}}$ the semi-norm of space ${W}^{2,\infty }({\rm{\Omega }})$.

To prove (29), assume $(\hat{K},{P}_{\hat{K}},{{\rm{\Sigma }}}_{\hat{K}})$ as defined in (Ciarlet 2002) is the reference element which all the finite elements are affine-equivalent to. On $\hat{K}$, we have continuous embeddings:

$$\begin{eqnarray}&&{W}^{2,\infty }(\hat{K})\,\hookrightarrow \,{C}^{s}(\hat{K}),s=0\ \mathrm{or}\ 1\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{W}^{2,\infty }(\hat{K})\,\hookrightarrow \,{H}^{1}(\hat{K}),\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&{P}_{1}(\hat{K})=\hat{P}\subset {H}^{1}(\hat{K})\end{eqnarray} \tag{ 32 }$$

when linear triangular or tetrahedron elements are employed. According to (Ciarlet 2002, Theorem 3.1.6), there exists constant C₃ depending only the reference element such that

$$\begin{eqnarray}&&\parallel v-{{\rm{\Pi }}}_{K}v{\parallel }_{1,K}\leqslant {C}_{3}{(\mathrm{meas}(K))}^{1/2}{h}_{K}| v{| }_{2,\infty ,K},\qquad \mathrm{for}\ \mathrm{any}\ v\in { \mathcal V }.\end{eqnarray} \tag{ 33 }$$

Here Π_K v is the P_K interpolate of v as defined in (Ciarlet 2002). h_K is the diameter of element K. The (33) together with

$$\begin{eqnarray}&&\parallel v-{{\rm{\Pi }}}_{h}v{\parallel }_{1,{\rm{\Omega }}}={\left(\displaystyle \sum _{K\in {{ \mathcal T }}_{h}}\parallel v-{{\rm{\Pi }}}_{K}v{\parallel }_{1,K}^{2}\right)}^{1/2}\qquad \mathrm{for}\ \mathrm{any}\ v\in { \mathcal V },\end{eqnarray} \tag{ 34 }$$

gives (29).

The (29) indicates that given  > 0 , we can find a h such that

$$\begin{eqnarray}&&\parallel {v}^{* }-{{\rm{\Pi }}}_{h}{v}^{* }{\parallel }_{1}\lt \epsilon /2.\end{eqnarray} \tag{ 35 }$$

From (28) and (35), we have

$$\begin{eqnarray}&&\parallel u-{{\rm{\Pi }}}_{h}{v}^{* }{\parallel }_{1}\leqslant \parallel u-{v}^{* }{\parallel }_{1}+\parallel {v}^{* }-{{\rm{\Pi }}}_{h}{v}^{* }{\parallel }_{1}\lt \epsilon .\end{eqnarray} \tag{ 36 }$$

From (36), we have

$$\begin{eqnarray}&&\parallel [(u,U)]-[({{\rm{\Pi }}}_{h}{v}^{* },U)]{\parallel }_{H/R}=\mathop{\inf }\limits_{c\in R}\parallel u-{{\rm{\Pi }}}_{h}{v}^{* }+c{\parallel }_{1}\leqslant \parallel u-{{\rm{\Pi }}}_{h}{v}^{* }{\parallel }_{1}\lt \epsilon .\end{eqnarray} \tag{ 37 }$$

To sum up, given any  > 0, we can find a h and $[({{\rm{\Pi }}}_{h}{v}^{* },U)]\in {H}_{h}/R$ such that (37) holds. So (26) is proved.

Third, (19) and (26) leads to

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{h\to 0}{\parallel [(u,U)]-[({u}_{h},{U}_{h})]\parallel }_{H/R}=0\end{eqnarray} \tag{ 38 }$$

With definition (20), (38) implies (13).□

The (13) is simplified to

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{h\to 0}\parallel U-{U}_{h}{\parallel }_{{R}^{L}}=0,\end{eqnarray} \tag{ 39 }$$

when the sum of voltages on all the electrodes are set to zero in practice. This is because

$$\begin{eqnarray}&&\mathop{\inf }\limits_{c\in R}\parallel U-{U}_{h}+c{\parallel }_{{R}^{L}}=\parallel U-{U}_{h}+{u}_{h,\mathrm{avg}}-{u}_{\mathrm{avg}}{\parallel }_{{R}^{L}}\end{eqnarray} \tag{ 40 }$$

with ${u}_{\mathrm{avg}},{u}_{h,\mathrm{avg}}$ the average of U and U_h respectively. With (40), the (13) becomes (39) when ${u}_{\mathrm{avg}}={u}_{h,\mathrm{avg}}=0$.

The (39) means that the finite element approximation for the voltages on the electrodes converge point-wisely to the true ones as the sizes of elements in the underlying mesh all tend to zero. Here convergence is only discussed with respect to the voltages on the electrodes. This is because it is these voltages together with current through electrodes that are normally measured in practice. Finally, rather than Ma (2017), the convergence result (13) here also holds true for EIT over a three-dimensional domain.

Simulations were also run to test the convergence of finite element approximation on a toy 2-D EIT model which is depicted in figure 1. The figure shows a disk of radius R that is made of material obeying the ohm's law. The disk has a constant conductivity distribution that is denoted by σ. On the boundary of this disk, two identical electrodes are placed. Current of amount I is injected into the region through the bottom electrode and extracted out of the region through the top one. Contact impedance denoted by z exists underneath each electrode. They are assumed to be a constant and is the same for each electrode. This EIT system is assumed to be perfectly modelled by the complete electrode model (1a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When voltages were measured, grounding was made at the rightmost point on the boundary of this disk (See figure 1). By doing this, voltages on the electrode and the potential distribution can be easily found with the Fourier method as shown later on.

In simulation, the BVP for this 2-D EIT model was solved by FEM, an Fourier method and the method in (Demidenko 2011). The finite element method was employed with linear triangle element. The computation were done on a set of nine unstructured meshes with varying densities. All these meshes were generated by the 'distmesh' mesh generator (Persson and Strang 2004). The number of elements in each mesh is contained in table 1. Two of the meshes are shown in figure 2. In all the nine meshes, elements that are close to the ends of any electrode are of small size. This non-uniformity in the sizes of elements is to accommodate the rapid change of current density around the ends of the electrodes.

Zoom In Zoom Out Reset image size

The Fourier method employed here is similar to (Paulson et al 1992), (Demidenko 2011). This method results in solving a system of linear equations with infinite number of equations each of which has infinite number of terms. When the first N terms were kept in the Fourier representation of the current density, the infinite system becomes the following finite system of linear equations

$$\begin{eqnarray}&&\displaystyle \frac{R}{\sigma }\displaystyle \sum _{n=1}^{N}\displaystyle \frac{a(k,n)}{n}{c}_{n}+\displaystyle \frac{\pi z}{2}{c}_{k}+{U}_{1}b(k)=0,k=1,2,\,\ldots ,\,N,\end{eqnarray} \tag{ 41a }$$

$$\begin{eqnarray}&&\displaystyle \sum _{n=1}^{N}b(n){c}_{n}=I,\end{eqnarray} \tag{ 41b }$$

with

$$\begin{eqnarray}&&a(k,n)={\int }_{\tfrac{\pi }{3}}^{\tfrac{2\pi }{3}}\sin (n\theta )\sin (k\theta ){\rm{d}}\theta \quad \mathrm{and}\quad b(k)={\int }_{\tfrac{\pi }{3}}^{\tfrac{2\pi }{3}}\sin (k\theta ){\rm{d}}\theta .\end{eqnarray} \tag{ 42 }$$

In (41a), the unknown are c_k = 1,2, ..., N and U₁ which is the voltage on the bottom electrode. Details about the derivation of this method is in appendix I. When the number of electrode increases and the contact impedance are different for different electrodes, an analytic solution for our 2-D EIT model is still available, e.g. in (Demidenko 2011). This solution can be used to test FEM and compared to the theoretical solution for more complex model.

When the analytic solution described in (Demidenko 2011) was employed for this toy model of EIT, the potential at the point with the polar coordinate (r, θ) is

$$\begin{eqnarray}&&u(r,\theta )=\displaystyle \sum _{i=1}^{\infty }{b}_{n}\sin (n\theta )\end{eqnarray} \tag{ 43 }$$

with infinite dimension vector ${\bf{b}}={({b}_{1},{b}_{2},\ldots )}^{T}$ satisfying

$$\begin{eqnarray}&&(2{{\bf{M}}}_{221}+\pi {\bf{n}}){\bf{b}}=2{U}_{1}{{\bf{r}}}_{21}.\end{eqnarray} \tag{ 44 }$$

In (44), M₂₂₁ is the infinite dimension matrix whose (k,n)th element is

$$\begin{eqnarray}{{\bf{M}}}_{221{kn}}=\left\{\begin{array}{ll}\displaystyle \frac{\pi }{6}-\displaystyle \frac{\cos (-k\pi )\sin (k\pi /3)}{2k}, & \mathrm{if}\ n=k\\ \displaystyle \frac{\cos ((n-k)\pi /2)\sin ((n-k)\pi /6)}{n-k}-\displaystyle \frac{\cos ((n+k)\pi /2)\sin ((n+k)\pi /6)}{n+k}, & \mathrm{if}\ n\ne k.\end{array}\right.\end{eqnarray} \tag{ 45 }$$

n =  diag{ 1, 2, ...}, and r₂₁ is the infinite dimension vector with components ${r}_{21n}=(2/n)\sin (\pi n/6)\sin (-\pi n/2)$.

The computation results from the three methods were compared in terms of the difference of voltages on the two electrodes which is denoted by ΔU. The ΔU computed by the three methods were shown in figure 3. These results were obtained when σ, R and z in the set-up of EIT model all took unit value in consistent units of measurement. The amount of current through the region was I ≈ 0.999 7. This value was considered to be the true amount of current when supplied voltages were U₁ = −U₂ = 1.629 5. This value was obtained with the analytic solution (43) where truncation was made after 13000 terms. In our analytic method, estimations of ΔU were recorded when truncation was made right after N = 500, 1000, 3000, 5000, 7000, 9000, 11000 terms respectively.

As seen in figure 3, the ΔU estimated by FEM and our analytic method tend to converge to its true value, as the underlying mesh get finer in FEM and more terms in the Fourier expansion are kept in the analytic method. The discrepancy between the two estimations is largely due to the error in the discretization of the disk region. Recall the theorem on FEM convergence holds under the assumption that the region is polygonal and discretized without any error. In our simulation, however, the region is a disk and was discretized into polygons which were all within it. Because of this, the estimated voltage difference by FEM is consistently smaller, when the amount of current passing through the region between the two electrodes is fixed.

The relative error

$$\begin{eqnarray}&&{\delta }_{{\rm{\Delta }}U}=| {\rm{\Delta }}{U}_{{FEM}}-2{U}_{1}| /(2{U}_{1})\end{eqnarray} \tag{ 46 }$$

of FEM approximation against the number of elements was shown in figure 4.

Zoom In Zoom Out Reset image size

As seen, the relative error can be made to be well within 0.1% with sufficiently fine meshes.

In the computation of the forward problem of EIT with the complete electrode model, the voltages on the electrodes that are estimated by the conforming finite element method with linear triangle or tetrahedron elements, was proved to converge point-wisely to the true ones up to an additive constant, as the sizes of elements in the underlying mesh all approach zero. Simulation on a toy EIT model for the convergence of the finite element approximation also produced results that were consistent with the theoretical prediction, when the error in the discretization of the region is considered.

Our result implies that theoretically we can make our finite element approximation of voltages on the electrodes as much accurate as we want by refining the underlying mesh which discretizes the region without any error. This may interest some practitioners. Because in practice high accuracy in the computation of the estimated voltage by FEM is desirable, when this quantity is compared with the measured voltages in some iterative reconstruction algorithms of EIT.

In the future, work could tackle the optimal mesh design for EIT because it plays the central role in the FEM approximation. A good design may provide a better approximation even with a smaller number of elements.

When the conductivity distribution σ over the region is a constant, the number of electrodes are L = 2 and the two electrodes are placed on the boundary in a way as shown in figure 1, the BVP (1a) for the toy EIT model becomes

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}u=0\end{eqnarray} \tag{ 47a }$$

$$\begin{eqnarray}&&-\sigma \displaystyle \frac{\partial u}{\partial n}=f(\theta )\end{eqnarray} \tag{ 47b }$$

$$\begin{eqnarray}&&u(R,0)=0,\end{eqnarray} \tag{ 47c }$$

with

$$\begin{eqnarray}f(\theta )=\left\{\begin{array}{ll}\displaystyle \frac{-u-{U}_{1}}{z} & \quad \mathrm{for}\ \pi /3\leqslant \theta \leqslant 2\pi /3\\ \displaystyle \frac{-u-{U}_{2}}{z} & \quad \mathrm{for}\ -2\pi /3\leqslant \theta \leqslant -\pi /3\\ 0 & \quad \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 48 }$$

and

$$\begin{eqnarray}&&{\int }_{\pi /3}^{2\pi /3}f(\theta ){\rm{d}}\theta =I.\end{eqnarray} \tag{ 49 }$$

Here f(θ) is the current density in the direction of out normal and is at the point on the boundary that has polar coordinate (R, θ). With (47c), it is easy to see that U₁ = −U₂ and f is an odd function of θ. So we assume

$$\begin{eqnarray}&&f(\theta )=\displaystyle \sum _{k=1}^{\infty }{c}_{k}\sin (k\theta ).\end{eqnarray} \tag{ 50 }$$

with coefficients c_k, k = 1, 2, ... to be determined.

Note that the BVP

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}u=0\end{eqnarray} \tag{ 51a }$$

$$\begin{eqnarray}&&-\displaystyle \frac{\partial u}{\partial n}=\sin (k\theta )\end{eqnarray} \tag{ 51b }$$

$$\begin{eqnarray}&&u(R,0)=0\end{eqnarray} \tag{ 51c }$$

has the solution

$$\begin{eqnarray}&&u(r,\theta )=-\displaystyle \frac{{r}^{k}\sin (k\theta )}{\sigma {{kR}}^{k-1}}.\end{eqnarray} \tag{ 52 }$$

With (50) and the principle of superposition, BVP (47a) has the solution

$$\begin{eqnarray}&&u(r,\theta )=-\displaystyle \sum _{k=1}^{\infty }{c}_{k}\displaystyle \frac{{r}^{k}\sin (k\theta )}{\sigma {{kR}}^{k-1}}.\end{eqnarray} \tag{ 53 }$$

The (53) leads to

$$\begin{eqnarray}&&u(R,\theta )=-\displaystyle \sum _{k=1}^{\infty }\displaystyle \frac{{c}_{k}R\sin (k\theta )}{\sigma k}=-\displaystyle \frac{R}{\sigma }\displaystyle \sum _{k=1}^{\infty }\displaystyle \frac{{c}_{k}}{k}\sin (k\theta ).\end{eqnarray} \tag{ 54 }$$

With (54), the (48) and the oddness of f(θ) indicates that

$$\begin{eqnarray}&&{c}_{k}={\int }_{0}^{\pi }f(\theta )\sin (k\theta ){\rm{d}}\theta =\displaystyle \frac{2}{\pi }{\int }_{\tfrac{\pi }{3}}^{\tfrac{2\pi }{3}}\displaystyle \frac{1}{z}\left(-\displaystyle \frac{R}{\sigma }\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{c}_{n}}{n}\sin (n\theta )-{U}_{1}\right)\sin (k\theta ){\rm{d}}\theta ,k=1,2,\ldots \end{eqnarray} \tag{ 55 }$$

On the other hand, the constraint (49) together with (50) requires

$$\begin{eqnarray}&&{\int }_{\pi /3}^{2\pi /3}\displaystyle \sum _{k=1}^{\infty }{c}_{k}\sin (k\theta ){\rm{d}}\theta =I.\end{eqnarray} \tag{ 56 }$$

The (55) and (56) forms a system of infinite number of linear equations with unknown c_k, k = 1, 2, ... and U₁. This system cannot be solved exactly. For approximation, truncation was made in the Fourier representation (50). Suppose it was truncated at the N-th term, then derivation similar to the above gives the system of linear equations (41a). The solution of this finite system are approximations to the value of c_k, k = 1, 2, ..., N and U₁.