Convergence study of 2D forward problem of electrical impedance tomography with high-order finite elements

A convergence study of the forward problem of electrical impedance tomography is performed using triangular high-order piecewise polynomial finite-element methods (p-FEM) on a square domain. The computation of p-FEM for the complete electrode model (CEM) is outlined and a novel analytic solution to the CEM on a square domain is presented. Errors as a function of mesh-refinement and computational time, as well as convergence rates as a function of contact impedance, are computed numerically for different polynomial approximation orders. It is demonstrated that p-FEM can generate more accurate forward solutions in less computational time, which implies more accurate simulated interior potentials, electrode voltages and conductivity Jacobians.


Introduction
Electrical impedance tomography (EIT) is an emerging imaging modality that aims to reconstruct the interior conductivity distribution of an object from electrical measurements obtained on the periphery. The use of non-invasive conductivity imaging has exciting biomedical applications, such as breast cancer detection, respiratory and brain function monitoring, as well as geophysical and chemical engineering applications (see [1] for a review of EIT application areas).
Mathematically, the objective of EIT is to recover the conductivity coefficient of a second-order linear elliptic partial differential equation (PDE) from knowledge of electrical measurements on the boundary. Even assuming full knowledge of boundary data, the recovery of the conductivity coefficient is severely ill-posed. Hence, the vast majority of research in the EIT community has focused on applying regularization techniques for stable recovery of the conductivity coefficient (see [2,3] for review articles on EIT reconstruction algorithms).
Many EIT reconstruction algorithms require the simulation of electrical boundary measurements from some conductivity sufficiently close to the true conductivity, σ T , known as the forward problem. Typically, a set of low frequency electric currents is applied to a set of electrodes on the boundary and for each current the resulting electric potential is measured on the electrodes. The forward problem is typically solved using piecewise linear finite elements, with a mesh sufficiently fine to ensure inaccuracies in the forward problem are significantly less than the noise in the electrical measurements, and a mesh generated with increased mesh density in a neighbourhood of each electrode where the gradient of the electric potential is largest. [4,5] The complete electrode model (CEM) is the most commonly used set of boundary conditions in EIT, which takes into account the formation of contact impedances at the electrode-object interface. [6] The CEM boundary conditions are non-local, since the input current is an integral over the electrode of a current density, and this has resulted in a difficulty in construction of analytic solutions. Semi-analytic solutions for the CEM on the unit disc were originally discussed in [7,8], and extended in [9] for homogeneous conductivities on the unit disc and more recently in [10] for circular inhomogeneities on the unit disc. In this paper, a novel semi-analytic solution to the CEM on a 2D square domain is computed, and subsequently used to perform a convergence study of p-FEM on this domain. It is further demonstrated that one can solve analytically for a basis of input electrode voltages for the interior potentials and electrode currents, and use these to synthesize electric potentials for an arbitrary input electrode current.
It is instructive to understand what norms to measure convergence in to be most relevant to improving EIT reconstruction. Given a vector of measured voltages U(σ T ) ∈ R m at σ T , a typical EIT reconstruction algorithm seeks to find the conductivity as the minimizer of a data misfit and regularization functional where R is a regularization functional used to stabilize the mapping from measured voltages to conductivity due to the inherent ill-posed nature of EIT. V ( · ) is a representation of the forward model and is a (non-linear) map from the best current guess of the conductivity, σ r , to a simulated voltage data vector, typically computed using FEM. If the voltage measurements are known to be Gaussian distributed, then the most appropriate misfit functional is a (weighted) 2-norm, and hence it is natural to consider convergence of simulated voltages in the 2-norm. A typical EIT measurement consists of a tetrapolar measurement in which current is driven between a pair of electrodes and a voltage measurement is recorded between a pair of electrodes. Voltage measurements on current active electrodes are subject to larger errors than current passive electrodes, and thus convergence is considered on both current active and current passive electrodes in the 2-norm. It is well known in EIT that voltage measurements are Fréchet differentiable with respect to L ∞ conductivity perturbations. [11,12] In iterative reconstruction algorithms, a discretized Fréchet derivative, the Jacobian matrix, is iteratively updated at the current estimate of the conductivity, σ . In integral form, the kernel of the Fréchet derivative is given by the inner product of gradients of electric potentials at σ , and hence it is natural to consider convergence of the interior potential in the H 1 norm. It has been shown in [13] that an accurate Jacobian is less important than accurate CEM electrode potentials during iterative algorithms, but convergence in the H 1 norm is still considered here for completeness.
The present author is aware that more advanced numerical finite-element methods, such as hp-FEM, are available to solve the forward problem. In hp-FEM, the solution can be updated locally both through a spatial refinement of the elements (h-refinement), when discontinuities are present in the solution, or increasing the approximating polynomial degree (p-refinement), where the solution has higher regularity. hp-FEM can achieve exponential convergence as a function of the number of degrees of freedom of the approximation (see [14,15] for an application of hp-FEM in 2D EIT for the continuum model and CEM, respectively, as well as [16] for an application of p-FEM in 2D EIT for the CEM). However, this is the first study to use a novel analytic solution to the CEM as a testbed for p-FEM convergence.
This paper is organized as follows. In Section 2, the governing PDE, electrode models and weak formulations of EIT forward problems are introduced. In Section 3, the computation of the CEM with p-FEM is outlined. In Section 4, a derivation of the novel semi-analytic solution to the CEM is performed (Section 4.2) followed by a p-FEM convergence study with different electrode models (Section 4.3).

Forward problem in EIT
Assuming time harmonic electric and magnetic fields, a second-order linear elliptic PDE can be derived from Maxwell's equations (see [17, p.37-39]) that governs the relationship between the interior voltage and conductivity distribution where ⊂ R n for n = 2, 3 is the domain of interest, u : → C is the electric potential and σ , : → R n×n are the conductivity and permittivity tensors. The frequency ω is assumed sufficiently small so that the permittivity can be neglected. σ is a symmetric positive definite matrix almost everywhere with σ ij ∈ L ∞ ( ), and in the isotropic case (σ ij (x) = s(x)δ ij ), this reduces to 0 < c ≤ s(x) ≤ C almost everywhere for constants c, C ∈ R, denoted σ ∈ L ∞ + ( ). Mathematically, the objective of EIT is to infer σ from (partial, noisy) knowledge of the Dirichlet-Neumann map on the boundary ∂ .

Boundary conditions: modelling electrodes
Electrode models have been studied extensively in [6,8] with the most prominent being the continuum, point electrode and complete electrode models. Denote the lth electrode as E l ⊂ ∂ , and let E T := l E l and E T := ∂ \ E T , respectively. Let ν denote the outward pointing unit normal at x ∈ ∂ . It is assumed that (σ ∇u) · ν| E T = 0 because the surrounding air is an insulator to a good level of approximation.

Continuum model
In this model, a current flux density, f, is prescribed as a Neumann condition on ∂ to the conductivity Equation (2) σ ∂u ∂n where f : ∂ → R is a current injection function, with ∂ f = 0 to conserve current. This model assumes all of ∂ is accessible to input current and so there is effectively an infinite number of electrodes, which is clearly unrealistic.

Point electrode model
The point electrode model (PEM) is a limiting case of the continuum model in which the current injection function f is of the form where δ is the Dirac delta distribution. This models L electrodes, all of infinitesimal size, where x l ∈ ∂ is the spatial coordinate of the lth electrode and L i=1 I l = 0 to conserve current. The PEM is unphysical because the electrodes are of finite size in practice.

Complete electrode model
When electrodes are placed in contact with an object, a contact impedance layer exists between the electrode and object. [6] If electrode l is a perfect conductor, the upper voltage on the electrode, U l , is constant and there is a voltage drop across the contact impedance layer. The contact impedance, z l , is assumed to be constant for each electrode and positive, z l > 0. This model is the most widely used in practical EIT systems (see [6] for experimental validation).
The CEM consists of the conductivity Equation (2) along with the boundary conditions along with the consistency condition L l=1 I L = 0. The first two equations above can be combined, to give the following equivalent form where λ(x) = L l=1 χ l , μ(x) = 1, and χ l is the characteristic function of the lth electrode.

Weak formulation
The weak formulation of the forward problems is derived by multiplying the conductivity Equation (2) by a test function, v, integrating over and invoking the divergence theorem (see e.g. [18,Chap. 6] for a general text on weak formulation for elliptic PDE).

Continuum and PEM
The weak formulation for the continuum model is: where It is well known there exists a unique solution (up to a constant) for this problem (see e.g. [18,Chap. 6]). Note that for the PEM, the current injection function in (4) lies just short of H − 1 2 (∂ ) (see Section 4.1.2), however, it is still common in the literature to consider the PEM having the same weak formulation as the continuum model.

Complete electrode model
The weak formulation is: where whereḢ( ) = (H 1 ( ) ⊕ R L )/R, reflecting the fact that the potential is only defined up to a constant. In [19], it is shown that this problem is uniquely solvable for (u, U) ∈Ḣ( ).

Methods
As discussed in the introduction, different forward problems can only be solved analytically for simple domain geometries, conductivities and boundary conditions, and in practice, a numerical method is almost always deployed. Different methods are available including the finite difference method (FDM), boundary element method (BEM) or the FEM. The FEM is the most widely used numerical method and has the advantage that it can be applied to domains with curved boundaries, inhomogeneous conductivities and results in a sparse system matrix (see e.g. [20] for general text on FEM). The implementation of high-order polynomial finite elements for the potential and anisotropic piecewise constant conductivities in 2D is now outlined (for full details on implementation see [21,Chap. 2]). These elements have been implemented in the open source reconstruction software EIDORS 3.8. [22,23] In the FE method, the domain ⊂ R 2 is decomposed into a set of E elements chosen here to be triangles. Within each element, each corner is called a vertex and there are V vertices in total over . On each element, a number of nodes are defined, of which the vertices are a subset, with N distinct nodes over . One looks for an approximate potential of the form where ψ i are known as the shape functions, satisfying the interpolation property where x i is the spatial coordinate of node i and δ ij is the Kronecker delta. From (11) the vector α = (α 1 , . . . , α N ) ∈ R N represents the discretized approximation to the potential, u. A piecewise constant conductivity approximation over elements is used, σ = E i=1 σ i C i , where C i and σ i ∈ R n×n are the characteristic function and restriction of the conductivity to the ith finite element, respectively.

Continuum and PEM FE approximation
In the Galerkin finite-element method, a finite dimensional subspace, S 1 h ( ) ⊂ H 1 ( ), is chosen to represent the approximate potential u h , and the test function v h chosen to lie in the same space. Substituting the expression for the approximate potential (11) into the weak formulation (8) and . . , N leads to the following system of equations for α where the matrix A and vector f have entries respectively. The RHS integral can be computed numerically for a given current injection f. In the special case of the PEM boundary condition (4), f has value I l at the index corresponding to the node of the lth electrode. The potential u is defined up to a constant resulting in a 1-dimensional null space of A. This is resolved by choosing, say, the gth interior node, with coordinate x g , to be at zero potential i.e. u(x g ) = 0. The gth row and column of A and the gth row of f and u are removed to generateÃ,f andα, respectively. The N − 1 dimensional systemÃα =f is now uniquely solvable forα.

CEM FE approximation
In the CEM, an approximate solution for the interior potential is sought as well as an approximation, U h , to the potentials on the electrodes where β ∈ R L . As first described by Vauhkonen [17, p.44], (α, β) is given by the solution of the system of equations where the matrix A is given by (14) and the matrices B, C and D by where i, j = 1, . . . , N and l = 1, . . . , L. This can be written compactly as Similarly to the continuum model, since the potential is only defined up to a constant, an N + L − 1 dimensional linear systemSb =Ĩ is solved forb, whereS,Ĩ andb are generated by removing the gth row and column of S and gth row of I and b, respectively.

Synthesizing applied current solutions from applied potential solutions
The typical forward problem to solve in EIT is to simulate a pair (u, U) given I. Due to non-local boundary condition for I i.e. the last equation in (5), a direct analytic solution for this problem is difficult to determine. Instead, a set of problems can be solved and subsequently used to synthesize a solution (ũ,Ũ) for an arbitraryĨ. The semi-analytic solution for a single pair (u i , I i ) given U i is postponed until Section 4.2.3 for a square domain, and the general synthesis procedure is now outlined. Assume a basis of electrode potentials {U i } L i=1 has been applied, and for each the interior voltage/current pair (u i , I i ) computed. Define the transfer admittance, Y σ ,z , as Each computed current is assembled to form the columns of the admittance To synthesise the solution for a current,Ĩ, this is equated to a linear combination of the computed currents so that for some coefficients c := (c 1 , . . . , c L ) ∈ R L .Ỹ σ ,z has a 1-dimensional null space since U is defined up to a constant. To compute the inverse i.e.Z σ ,z :=Ỹ −1 σ ,z , the singular value decomposition is computed numerically and subsequently the Moore-Penrose pseudoinverse,Ỹ + σ ,z , computed to give c =Ỹ + σ ,zĨ , and thus (ũ, where γ ∈ R L , then the solution of (α, γ ) given U is given by where 1 L is the L × L identity matrix. Thus, an approximate FEM solution (ũ h ,Ũ h ) for a driven currentĨ can also be synthesized from the FEM solution pairs given a basis of input potentials {U i } L i=1 , making this synthesis method possible if one does not have access to analytic forms for the currents or one only has access to a FEM solver that solves local Robin problems.

Reference elements and boundaries
The computation of the stiffness matrix in (14) using reference element transformations is now outlined (for more details, see e.g. [24,Chap. 1]). This well-known transformation is outlined since analogous transformations are required to compute integrals appearing in the CEM matrices. Let d be the dimension of the set of polynomials of degree p in n dimensions (see Section 3.5). The element shape functions are defined as the restriction of the d global shape functions ψ i in (12). These form a basis set

Stiffness matrix
Using a piecewise-constant conductivity basis, (14) can be written and hence A is assembled from elemental stiffness matrices Define a linear map F k : R 2 → R 2 from a reference element, R, with local coordinates ( , η) chosen as the unit triangle, to the kth element with local coordinates (x, y). The simplest map is where i is the linear basis function associated with the ith reference element vertex. The basis functions i satisfy (12) and are uniquely determined as This linear map sends a straight sided reference element to a straight sided global element. Equipped with this transformation, the derivatives of the shape functions with respect to global coordinates (x, y) are transformed to derivatives with respect to local coordinates ( , η) through the Jacobian matrix, J = ∂(x,y) ∂( ,η) . The gradient operator in global coordinates, , is written in terms of the local gradient operator, where |J k | is the absolute value of the determinant of the Jacobian matrix, R is the reference element and

CEM matrices
The forward problem for the CEM requires the additional computation of the matrices B, C, and D in Equation (17). The matrix D is independent of the choice of basis functions and is trivial to compute, but B and C are line integrals for a 2D problem. During FE meshing, it is convenient for electrodes to be the union of some triangular element edges in 2D so that integrals over electrodes are the sum of integrals over these edges. Let n l be the number of edges that the lth electrode contains, and l(b p ) denotes the pth edge of the lth electrode. The matrix entries of B and C in (17) can be written respectively. Assume the pth edge of the lth electrode is an edge of the kth element. Let (24), be the parametrization of the pth edge. The partial derivative of G k and magnitude are given by The line integrals are given by The only contributions are from basis functions ψ i whose associated node x i lies on l(b p ).
The line integrals appearing in the matrix entries of B and C in (27) are written respectively. After the integrals have been evaluated, the matrices B and C are assembled in an analogous manner to the stiffness matrix.

Lagrange elements
Denote P n p as the set of polynomials of degree p in n dimensions. The dimension of P n p is d n p = n+p n , which is the number of coefficients to represent such a polynomial, and when n = 2, d 2 p = 1 2 (p + 2)(p + 1). The kth element has d 2 p local degrees of freedom, corresponding to the number of nodes used in the element to interpolate the potential, and there are d basis functions not identically 0 on the element.
For p = 1, the coefficients are uniquely determined from the interpolation property at the nodes (12), and thus shape functions are taken as the mapping basis functions for the potential on the element, that is ψ (24). For p > 1, there is freedom to choose nodal positions within the reference element satisfying the interpolation property (12). Lagrange elements, where nodes are placed nodes equispaced within an element, are implemented (see e.g. [25, p.72-74]).

Numerical quadrature
Gauss quadrature is used to evaluate the integrals (26), (29) and (30). An N point quadrature rule for h ∈ P n p (R) where R ⊂ R n , consists of a set of weights, w i , and Gauss points, g i ∈ R. The smallest integer k n such that is the quadrature order. 1D quadrature rules for polynomials of degree 2p and p are required for B and C in (29) and (30), respectively, and hence for p = 1, 2 and 3 rules for degree 2, 4 and 6 polynomials are sufficient. 2D quadrature rules for polynomial of degree 2(p − 1) are required for A in (26), and hence for p = 1, 2 and 3 rules for degree 2, 4 and 6 polynomials are sufficient. For tabulated Gauss points and weights, see e.g. [26, p.357-361].

Convergence study
In this section, a convergence study under uniform refinement of a square domain for the continuum model, PEM and CEM is performed. First solution regularity and finiteelement approximation error estimation results are discussed.

Finite-element error and regularity estimates
There are classical error estimates for high-order FE approximations for elliptic problems.
Let be a polygonal domain, and assume the solution u ∈ H 1 ( ) has increased regularity u ∈ H p+1 ( ), with p ≥ 1 an integer, and let u p h be the pth order finite-element approximation to u. Then where h is largest edge in the FE mesh and C 1 , C 2 > 0 are constants independent of h (see [20, p.139] and [20, p.134] respectively), and the semi-norm | · | p+1 is defined as Different EIT forward problems do not necessarily have such high regularity as is discussed on an individual basis shortly. Convergence results have been extended to solutions with fractional Sobolev regularity [27, p.14]. In particular, let u ∈ H s ( ) and s ∈ R with p + 1 ≥ s ≥ 1, then where C 3 , C 4 > 0 are constants independent of h. Thus if the regularity of u ∈ H s ( ) is such that p + 1 ≥ s, convergence rates independent of the degree p are expected.
Although the solutions of EIT forward problems are not necessarily smooth globally, there is a notion of elliptic regularity for solutions of elliptic PDE. In particular, the solution operator for an elliptic PDE preserves singular support, where the singular support is the complement of the smallest open set on which the function is smooth [28, p.246]. Hence, away from ∂ and in regions when the conductivity is smooth (C ∞ ) then the solution u is smooth (C ∞ ).

Continuum model
If ∂ is smooth and the conductivity coefficient and boundary condition are smooth functions, then u ∈ C ∞ ( ) and exponential convergence is guaranteed in the H 1 and L 2 norm for increasing p at a fixed h by (32). If the conductivity coefficient is discontinuous, u is continuous (C 0 ) along the discontinuity, but due to elliptic regularity, u is smooth (C ∞ ) away from the discontinuity. Hence if a finite-element triangulation conforms to a discontinuity, exponential convergence for increasing p at a fixed h is expected.

Point electrode model
For the PEM, a Dirac delta distribution, δ (n) , is prescribed as Neumann condition to the conductivity equation. δ (n) belongs to the dual space of H (n−1)/2+ (∂ ), δ (n) ∈ H −(n−1)/2− (∂ ) for all > 0, and for n ≥ 2 a careful analysis demonstrates that u ∈ H min{(4−n)/2− ,1} ( ) for all > 0. [29] Hence, in 2D, u ∈ H 1− ( ) for all > 0, and is not necessarily bounded in the energy norm, highlighting the non-physicality of the PEM. Global convergence is not expected in either the L 2 or H 1 norm. However, due to elliptic regularity, and assuming the conductivity is smooth, the solution will be smooth away from current driven electrodes, and hence exponential convergence for increasing p at a fixed h away from the electrodes in both the L 2 and H 1 norms is expected.

Complete electrode model
In [30, p.5, Remark 1] an analysis of the CEM for σ ∈ C ∞ ( ) shows that u ∈ H 2− ( ) for all > 0. Hence from (34) optimal convergence rates of O(h 1− ) and O(h 2− ) as h → 0 in the H 1 and L 2 norms respectively for all approximation orders are predicted. Further theoretical results [31] give tighter bounds on the regularity of u. Let the vector of contact impedances have entries that are constant still denoted by z. The authors demonstrate that, for any s ∈ [0, 1 2 ), u ∈ H 3 2 +s ( ) and ||u|| and for any ∈ (0, 1 − s), with C 5 , C 6 > 0 are independent of h. These theoretical results demonstrate that as the contact impedance decreases, the solution becomes unbounded in the H 3 2 +s -norm, and thus smaller convergence rates than O(h 1− ) and O(h 2− ) as h → 0 in the H 1 and L 2 norms respectively are likely to be obtained.

Analytic solutions on square domain
We consider the conductivity equation on = [−π, π] 2 for the continuum model, PEM and CEM. The conductivity is given a jump discontinuity at y = y and is of the form σ = σ r y < y In all cases, Neumann zero boundary conditions on the y = π, x = 0 and x = π boundaries of are considered. After separation of variables, an analytic expression for u is given by (n(π − y)) + N n sinh (n(π − y))) cos (nx) y ≤ π 6 .
(38) The constants M 0 and M 0 are determined from the choice of ground point and there are three sets of equations coupling the coefficients M n , M n and N n for n ≥ 1. The first two sets of equations are given by enforcing continuity of u and ∂u ∂n along y = π 6 leading to M n = M n + γ n N n , γ n M n = σ r (γ n M n + N n ), n = 1, 2, 3, . . .
where γ n = tanh (n π 6 ). The last equation coupling the coefficients comes from the particular boundary condition of a given electrode model on the y = 0 boundary. Analytic solutions for specific continuum, PEM and CEM boundary conditions are discussed in Sections 4.2.1 and 4.2.3, followed by a p-FEM convergence study with these solutions in Section 4.3.

Continuum model
For the continuum model, the coefficients M n and N n , n ≥ 1, are related by where f (x) is the Neumann condition along y = 0. We consider a smooth current injection function f (x) = cos (x) on y = 0, resulting in 3 non-zero Fourier coefficients, M 1 , M 1 , N 1 . From Equation (39), the following linear system ⎡ for the non-zero coefficients M 1 , M 1 and N 1 results.

Point electrode model
For the PEM, we let σ r = 1 which implies M n = M n and N n = 0 for all n ≥ 1. For a current injection function 3 ) the coefficients are given by The solution to the PEM is unbounded in the H 1 norm as discussed in the previous section. To efficiently compute the infinite series, it is instructive to split into the interior and boundary nodes. On the y = 0 boundary As n → ∞, tanh (nπ) → 1 and the series is a slowly converging harmonic-type series. Away from the boundary y > 0 which converges quickly as n → ∞ because of the dominant sinh (nπ) term.

Complete electrode model
As discussed in the introduction, a similar technique to [9,10] is applied to obtain a novel semi-analytic solution to the CEM on a square domain with a discontinuous coefficient. The solution is only semi-analytic because an infinite dimensional system of equations to solve for the Fourier coefficients results. As discussed in Section 3.3, if one can solve for (u, I) given U, then L systems for {(u i , can be used to synthesize an analytic solution (ũ,Ũ) for any desired current input current patternĨ. The semi-analytic solution of a single pair (u i , I i ) given U i is now considered.
Three electrodes are placed on the y = 0 boundary with centres (x 1 , x 2 , x 3 ) = ( 3π 14 , π 2 , 11π 14 ) and half-widths w i = π 7 for i = 1, . . . , 3 (see Figure 2). The interior potential solves Laplace's equation which implies the analytic solution u(x, y) still has the form (38). In general, for the continuum model with a jump discontinuity, a 3 × 3 system (41) must be solved for each Fourier coefficient,f (n) of the current injection function f (x). However, the situation is more complicated for the CEM, since a dense system of equations for the Fourier coefficients on the y = 0 boundary results i.e. the Fourier coefficients are coupled. The jump discontinuity implies that the solution has an exponentially growing and decaying component, and it is informative to write the solution as u(x, y) = M 0 + ∞ n=1 M n cosh (n(π − y)) cos (nx) y ≥ π 6 M 0 + ∞ n=1 (M n e −n(π−y) + N n e n(π−y) ) cos (nx) y ≤ π 6 ,  (1) and T (2) triangulations for the PEM. The green circles on the y = 0 boundary are point electrodes. The bottom two figures are the T (1) and T (2) triangulations for the CEM. The green lines on y = 0 boundary are driven electrodes, the red region is a neighbourhood of driven electrodes and the union of the red and yellow region is where σ = σ r .  instead of with (M n cosh (n(π − y)) + N n sinh (n(π − y))) in (38). Equating the solution u and σ ∂u ∂n along y = y gives M n cosh (n(π − y )) = M n e −n(π−y ) + N n e n(π−y ) M n sinh (n(π − y )) = σ r ( − M n e −n(π−y ) + N n e n(π−y ) ), and combining these two equations M n = c n N n where c n = e 2n(π−y ) σ r − tanh (n(π − y )) σ r + tanh (n(π − y )) .
The CEM boundary conditions are incorporated by substituting (38) into the first equation of (6), and further using (46) gives for x ∈ [0, π]  where e n = (1+c n e −2nπ ) and f n = n(1−c n e −2nπ ). Multiplying by cos (mx) and integrating over [0, π] gives ∞ n=1 π 0 f n (N n e nπ )σ r cos (nx) cos (mx) dx The orthogonality of cos (nx) over [0, π] implies the left hand side is (N n e nπ )σ r f n δ nm π 2 . The right-hand side is a linear combination of N n weighted by integrals of cosine products which can be computed analytically. An infinite system of equations for the coefficients N = (N 0 , N 1 , N 2 for m, n = 0, 1, 2 . . .. The current on the lth electrode is then computed as The linear system of equations are inverted numerically by truncating to a system of (N T + 1) equations for N = (N 0 , N 1 , N 2 , . . . , N N T ) ∈ R N T +1 . The matrix P is ill conditioned due to the exponentially growing factors e nπ along each row. However P can be factorized as P = P Q, where Q is diagonal with entries Q nn = e nπ , and the inverse is given by P −1 = Q −1 P −1 , where the matrix P has improved conditioning. The matrix Q is still ill conditioned, but is diagonal and thus has an explicit inverse of the form (Q −1 ) nn = e −nπ . Effectively, N n e nπ is solved for instead of N n .
To choose a truncation level N T , note that by [30,p.5,Remark 1] if σ ∈ C ∞ ( ) then u ∈ H 2− ( ) for all > 0. By the trace theorem u| ∂ ∈ H for all > 0 and some constant C 7 > 0 independent of n. Since (N n e nπ ) in (45) is the Fourier cosine coefficient of u(x, 0) to leading order, the system PN = R is solvable at a truncation level N T . Clearly, N T must be chosen such that N T π w where w is the smallest electrode half-width. The largest inaccuracies in the electric potential are at the electrode edges (where the current density has a jump discontinuity), and thus to compute a feasible N T a potential U = ( − 1, 1, 0) is applied, and the system P (QN ) = R is solved for a range of N T values. For each solution, QN ∈ R N T +1 , u N (x, 0) is computed through (45). Figure 1 illustrates plots of , for x p = 2π 7 , 3π 7 and 4π 7 . When N T = 3000, the maximum pointwise relative error is approximately ≈ 10 −6 , and this truncation level is used in subsequent computations.

Numerical convergence on square domain
A sequence of FE triangulations, T (i) , indexed by i ≥ 1, of the square domain = [0, π] 2 is generated (see Figure 2). Each finite-element triangulation has a maximum associated element length h i , where h is the maximum element length in the coarsest mesh, T (1) . The global linear, quadratic and cubic finite-element errors are computed as a function

Continuum model
First consider the current injection function f (x) = cos (x) on the y = 0 boundary and let σ r = 1 which implies that solution is smooth (u ∈ C ∞ ( )). The top two figures in Figure 3 demonstrate convergence for each p as a function of h and errors approximately of the form O(h p+1 ) and O(h p ) as h → 0 in the L 2 and H 1 norm respectively are observed. These results demonstrate that global high-order finite-element approximations are more accurate at a given triangulation, for this smooth boundary condition, and are in agreement with the error estimates (32). Secondly, a jump discontinuity σ r = 5 at y = π 6 is considered. The solution u is only C 0 along y = π 6 , but is smooth in each region y < π 6 and y > π 6 due to elliptic regularity. and H 1 norm, respectively. Since the finite-element mesh refinement conforms to the discontinuity exponential convergence is retained due to elliptic regularity.

Point electrode model
We now consider the PEM current injection function 3 ) on the y = 0 boundary, with σ r = 1. The top two figures in Figure 4 are the L 2 and H 1 errors are plotted for each p as a function of h. The PEM has approximately O(h) convergence as h → 0 in the L 2 norm, independent of the polynomial approximation, and there is no convergence in the H 1 norm because the analytic solution is unbounded in this norm. These results are in agreement with the theory in Section 4.1.2 because, in 2D, the solution has regularity u ∈ H 1− ( ) for all > 0, and thus ∇u ∈ L 2 . This example highlights the danger of using the PEM since a globally accurate conductivity Jacobian cannot be computed with this model.
The bottom two figures in Figure 4 further illustrates convergence plots of the errors, neglecting a neighbourhood of elements near the boundary (see Figure 2). In the interior, there is convergence of the form O(h p+1 ) and O(h p ) as h → 0 in both the L 2 and H 1 norm due to elliptic regularity of the solution operator (see Section 4.1). This is an interesting observation for the inverse problem, since if the conductivity is known in a neighbourhood of the boundary, this suggests that the conductivity Jacobian can still be accurately computed in the interior.

Complete electrode model
For the CEM L semi-analytic solution pairs (u i , I i ) are computed from a basis of potentials , where e i is ith canonical basis vector in R L (with e L+1 := −e 1 ). Using the method described in Section 3.3, the solution pairs (u i , I i ) are used to compute a semi-analytic solution (ũ,Ũ) for a pair drive current between E 2 and E 3 i.e.Ĩ = (0, 1, −1). A FEM solution for this pair drive current, (ũ h ,Ũ h ), is also computed through (18). Convergence of the electric potential at measurement electrodes is of additional interest, and so convergence of the electrode potential is measured in the 2-norm for both a current passive electrode, E 1 , and a current active electrode, E 2 , which are denoted as the current-passive 2-norm and current-active 2-norm, respectively. We consider a specific jump discontinuity along y = 2π 7 with σ r = 5. Figures 5 and 6 display errors under uniform h-refinement, for each p as a function of h, with z = 1000 and z = 0.1. Firstly, for a given triangulation, the errors decrease as p increases. If z = 1000, when p = 2, 3 approximate convergence of the form O(h 2 ) and O(h) as h → 0 is observed in the L 2 and H 1 norm respectively but for p = 1 the rates are less in both norms. The convergence rate in the current-active 2-norm is also approximately of the form O(h 2 ) and O(h) as h → 0, whereas as the convergence rate in the current-passive 2-norm is significantly larger. The convergence rate is also less when the contact impedance is smaller, which is in agreement with the theory discussed in Section 4.1.3. Figure 7 displays errors away from driven electrodes for z = 1000. Convergence of the form O(h p+1 ) and O(h p ) as h → 0 is retained in the L 2 and H 1 norm respectively due to elliptic regularity. Figures 8 and 9 display errors under uniform h-refinement for each p as a function of the FEM system matrix assembly time, T A , and subsequent FEM solution time using the MATLAB backslash operator, T S , respectively with z = 1000. As p increases, the same level of accuracy is observed in each norm in less time for both assembly and solution. Figures 10 and 11 display errors, under uniform h-refinement, for each p as a function of the total time, T = T A + T S , with z = 1000 and z = 0.1 respectively. As p increases, the same level of accuracy is observed in each norm in less computational time. Figure 12 displays the convergence rate for each p as a function of z in different norms, with σ r = 1. The convergence rate increases as z increases for all p, which is in agreement with the theoretical discussion in Section 4.1.3, and the rates are in agreement with (35) and (36). The quadratic and cubic rate is significantly larger than the linear rate in the H 1 , L 2 norms and current-active 2-norms, especially when z > 0.01. The convergence rates of cubic approximation are slightly greater than quadratic convergence rates for z < 0.01 in the H 1 , L 2 norms and the current-active 2-norm. However, for z > 0.01, the quadratic rates are greater than cubic rates in all norms, even though the magnitude of the errors for cubic approximation were less than quadratic approximation as was observed in Figures 5 and 6. Significantly larger convergence rates are observed in the current-passive 2-norm than the H 1 , L 2 norms and current-active 2-norm as expected, since the largest FEM errors are situated at the edge of current active electrodes.

Conclusions
A convergence study of different EIT forward problems was performed on a 2D square domain using p-FEM, and a novel semi-analytic solution to the CEM on the square domain was derived. For the PEM global convergence is not obtained in the H 1 norm and hence a globally accurate conductivity Jacobian cannot be computed with this model. For the CEM higher order polynomial approximations resulted in greater convergence rates under uniform h-refinement. For a given finite-element triangulation, the errors were consistently smaller in all norms with increasing polynomial approximation order, and the same level of error was obtained in all norms in less computational time with higher order polynomials. The results suggest that more accurate simulated potentials can be obtained in less computational time with p-FEM. The convergence rate increased as the contact impedance increased, and for z > 0.01, the convergence rate of quadratic and cubic is significantly greater than linear. The same convergence order was observed when discontinuous conductivities were used if the finite-element triangulation conformed to the discontinuity. It is of interest to determine analytic solutions to the CEM in 3D to understand if the results of this convergence study extends to domains with a realistic number of spatial dimensions.