{'fraction_non_alphanumeric': 0.048101653398652786, 'fraction_numerical': 0.022535211267605635, 'mean_word_length': 4.108869075551384, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Fourier analysis for discontinuous conductivities: computational results

Two reconstruction methods of Electrical Impedance Tomography (EIT) are numerically compared for nonsmooth conductivities in the plane based on the use of complex geometrical optics (CGO) solutions to D-bar equations involving the global uniqueness proofs for Calder\'on problem exposed in [Nachman; Annals of Mathematics 143, 1996] and [Astala and P\"aiv\"arinta; Annals of Mathematics 163, 2006]: the Astala-P\"aiv\"arinta theory-based"low-pass transport matrix method"implemented in [Astala et al.; Inverse Problems and Imaging 5, 2011] and the"shortcut method"which considers ingredients of both theories. The latter method is formally similar to the Nachman theory-based regularized EIT reconstruction algorithm studied in [Knudsen, Lassas, Mueller and Siltanen; Inverse Problems and Imaging 3, 2009] and several references from there. New numerical results are presented using parallel computation with size parameters larger than ever, leading mainly to two conclusions as follows. First, both methods can approximate piecewise constant conductivities better and better as the cutoff frequency increases, and there seems to be a Gibbs-like phenomenon producing ringing artifacts. Second, the transport matrix method loses accuracy away from a (freely chosen) pivot point located outside of the object to be studied, whereas the shortcut method produces reconstructions with more uniform quality.

New numerical results are presented using parallel computation with size parameters larger than ever, leading mainly to two conclusions as follows. First, both methods can approximate piecewise constant conductivities better and better as the cutoff frequency increases, and there seems to be a Gibbs-like phenomenon producing ringing artifacts. Second, the transport matrix method loses accuracy away from a (freely chosen) pivot point located outside of the object to be studied, whereas the shortcut method produces reconstructions with more uniform quality.

Introduction
We study a widely applicable nonlinear Fourier transform in dimension two. We perform numerical tests related to the nonlinear Gibbs phenomenon with much larger cutoff frequencies than before. Furthermore, we compare two computational inverse transformations, called low-pass transport matrix method and shortcut method in terms of accuracy.
The inverse conductivity problem of Calderón [10] is the main source of applications of the nonlinear Fourier transform we consider. Let Ω ⊂ R 2 be the unit disc and let σ : Ω → (0, ∞) be an essentially bounded measurable function satisfying σ(x) ≥ c > 0 for almost every x ∈ Ω. Let u ∈ H 1 (Ω) be the unique solution to the following elliptic Dirichlet problem: The inverse conductivity problem consists on recovering σ from the Dirichlet-to-Neumann (DN) map or voltage-to-current map defined by Here ν is the unit outer normal to the boundary. Note that the map Λ : σ → Λ σ is nonlinear. The inverse conductivity problem is related to many practical applications, including the medical imaging technique called electrical impedance tomography (EIT). There one attaches electrodes to the skin of a patient, feeds electric currents into the body and measures the resulting voltages at the electrodes. Repeating the measurement with several current patterns yields a current-to-voltage data matrix that can be used to compute an approximation Λ δ σ to Λ σ . Since different organs and tissues have different conductivities, recovering σ computationally from Λ δ σ amounts to creating an image of the inner structure of the patient. See [41,12] for more information on EIT and its applications.
Recovering σ from Λ δ σ is a nonlinear and ill-posed inverse problem, whose computational solution requires regularization. Several categories of solution methods have been suggested and tested in the literature; in this work we focus on so-called D-bar methods based on complex geometrical optics (CGO) solutions. There are three main flavors of D-bar methods for EIT: • Schrödinger equation approach for twice differentiable σ. Introduced by Nachman in 1996 [43], implemented numerically in [45,40,31,30,34,41]. • First-order system approach for once differentiable σ. Introduced by Brown and Uhlmann in 1997 [7], implemented numerically in [32,36,19].
Using these approaches, a number of conditional stability results have been studied for the Calderón problem. The most recent results in the plane were obtained by Clop, Faraco and Ruiz in [13], where stability in L 2 -norm was proven for conductivities on Lipschitz domains in the fractional Sobolev spaces W α,p with α > 0, 1 < p < ∞, and in [15], where the Lipschitz condition on the boundary of the domain was removed. In dimension d ≥ 3, conditional stability in Hölder norm for just C 1+ε conductivities on bounded Lipschitz domains was proved by Caro, García and the third author in [11] using the method presented in [18].
The three aforementioned D-bar methods for EIT in the two-dimensional case are based on the use of nonlinear Fourier transforms specially adapted to the inverse conductivity problem. Schematically, the idea looks like this: The main point above is that the nonlinear Fourier transform can be calculated from the infinite-precision data Λ σ , typically via solving a second-kind Fredholm boundary integral equation for the traces of the CGO solutions on ∂Ω.
In practice one is not given the infinite-precision data Λ σ , but rather the noisy and finite-dimensional approximation Λ δ σ . Typically all we know about Λ δ σ is that Λ σ − Λ δ σ Y < δ for some (known) noise level δ > 0 measured in an appropriate norm · Y . Most CGO-based EIT methods need to be regularized by a truncation |k| < R in the nonlinear frequency-domain as illustrated in Figure 1.
The regularization step results in a smooth reconstruction. This smoothing property of the nonlinear low-pass filter resembles qualitatively the effect of linear low-pass filtering of images. In particular, the smaller R is, the blurrier the reconstruction becomes.
The cut-off frequency R is determined by the noise amplitude δ, and typically R cannot exceed 7 in practical situations. However, it is interesting to understand how the conductivity is represented by its nonlinear Fourier transform. To study this numerically, we compute nonlinear Fourier transforms in large discs such as  Figure 1: Schematic illustration of the nonlinear low-pass filtering approach to regularized EIT. This image corresponds to the so-called shortcut method defined below. The simulated heart-and-lungs phantom σ(z) (bottom left) gives rise to a finite voltage-to-current matrix Λ δ σ (orange square), which can be used to determine the nonlinear Fourier transform. Measurement noise causes numerical instabilities in the transform (irregular white patches), leading to a bad reconstruction (a). However, multiplying the transform by the characteristic function of the disc |k| < 4 yields a lowpass-filtered transform, which in turn gives a noise-robust approximate reconstruction (b).
|k| < 60 and compute low-pass filtered reconstructions using inverse nonlinear Fourier transform and observe the results.
We compare two reconstruction methods based on the use of the CGO solutions (2.1): the low-pass transport matrix method implemented in [3], and a shortcut method based on solving a D-bar equation. The latter method does not have rigorous analysis available yet, but it is formally similar to the regularized EIT reconstruction algorithm studied in [43,45,40,31,30,33,34].
Our new computational findings can be roughly summarized by the following two points. First, both methods can approximate piecewise constant conductivities better and better as the cutoff frequency R increases, and there seems to be a Gibbs-like phenomenon producing ringing artifacts. Second, the transport matrix method loses accuracy away from a (freely chosen) pivot point located R = 20 R = 40 R = 50 True Figure 2: Convergence of nonlinearly low-pass filtered conductivity to the true discontinuous conductivity. Note the ringing artefacts reminiscent of the Gibbs effect of linear low-pass filtering. This picture illustrates the so-called shortcut method, and it is the first-ever computation of this kind with cutoff frequencies larger than 20.
outside of Ω, whereas the shortcut method produces reconstructions with more uniform quality. Figure 2 shows some of our numerical results via the shortcut method for a nonsymmetric conductivity distribution.
The rotationally symmetric examples presented in Section 4.1 provide numerical evidence of the fact that discontinuous conductivities can be reconstructed with the shortcut method more and more accurately when R tends to infinity.
In the computational experiments we show in this paper, the truncated scattering transform is computed with a very large cutoff frequency through the Beltrami equation solver. In practice, it is not possible to obtain such scattering data from boundary measurements unless we had unrealistically high-precision EIT measurements available.
Notwithstanding the aforesaid, our new results are relevant for applications at least in these two ways: (1) We provide strong new evidence suggesting that the shortcut method can recover organ boundaries at the low-noise limit. Even with practical noise levels we can recover low-pass filtered jumps in the cases of nested inclusions or checkerboard-type discontinuity curves. While there are many dedicated methods for inclusion detection in EIT, for example [16,9,28,8,29,25,27,39,21,17,46,22,26,23], none of them can locate nested inclusions. (2) The nonlinear Fourier transform studied here can be used to solve the nonlinear Noviko-Veselov (NV) equation. It is a nonlinear evolution equation generalizing the celebrated Korteweg-de Vries (KdV) equation into dimension (2+1). There has been significant recent progress in linearizing the NV equation using inverse scattering methods, see [38,37,44,42,14].
In the potential applications of the NV equation there should be no restriction to as low frequencies as in EIT. We remark that all our numerical experiments have been computed using either a 192 Gigabyte computer with Linux or parallel computation provided by Techila. 1

The nonlinear Fourier transform
In dimension two Calderón's problem [10] was solved using complex geometrical optics (cgo) solutions by Astala and Päivärinta [6]. In the case of L ∞conductivities σ the cgo solutions need to be constructed via the Beltrami equation Here z = z 1 + iz 2 ∈ C, ∂ z = (∂/∂z 1 + i∂/∂z 2 )/2 and k is a complex parameter. More specifically, for σ ∈ L ∞ (C) with c −1 ≤ σ(z) ≤ c almost every z ∈ C and σ ≡ 1 outside the unit disc, and for κ < 1 with |µ(z)| ≤ κ a.e. z ∈ C, Theorem 4.2 in [6] establishes that for any 2 < p < 1 + 1/κ and complex parameter k, there exists a unique solution f µ (·, k) in the Sobolev space W 1,p loc (C) such that holds with the asymptotic formula In addition, f µ (z, 0) ≡ 1. Write f −µ for the solution to the Beltrami equation with coefficient −µ. The scattering transform τ of Astala and Päivärinta (AP) is given by the formula Defining the functions as follows for k, z ∈ C: it turns out that u 1 , u 2 are the unique solutions to the following elliptic problems: The complex parameter k above can be understood as a nonlinear frequencydomain variable. As well as the function τ (k), defined in (2.4), plays the role of the nonlinear Fourier transform in the shortcut method described below, it is played by the following function in the low-pass transport matrix method presented in Section 3.1, where z 0 is a pivot point outside Ω. A practical algorithm was introduced in [5, 6, 3] for recovering σ approximately from a noisy measurement of Λ σ . This algorithm involves a truncation operation (nonlinear low-pass filtering) in the frequency-domain.
The scattering transform used in [43] is defined as where q = σ −1/2 ∆σ 1/2 and σ is a strictly positive-real valued function in W 2,p (R 2 ) for some p > 1 with σ ≡ 1 on the unit disc near the boundary and outside the unit disc, and m solves the equation Here * denotes convolution over the plane and g k the Faddeev fundamental solution of −∆ − 4ik∂. For smooth enough σ one can prove that In the rest of the paper, for non-smooth conductivities t denotes the function defined through formula (2.10). Why are t and τ called nonlinear Fourier transforms? There is admittedly some abuse of terminology involved, but the main reason is this. If one replaces m in (2.9) with 1, the result is the linear Fourier transform of q at (−2k 1 , 2k 2 ) ∈ R 2 . But m is only asymptotically close to 1 and actually depends on q.
The following theorem holds: is real-valued, radially symmetric, bounded away from zero and σ ≡ 1 outside the unit disc. Then the function t is real-valued and rotationally symmetric.

Inverse transforms and reconstruction algorithms
3.1. The low-pass transport matrix method. Given noisy EIT data Λ δ σ , we can evaluate the traces of M µ ( · , k) = 1 + ω( · , k) at the boundary ∂Ω for |k| < R with some R > 0 depending on the noise level δ. This is achieved by the numerical method described in Section 4.1 of [3] based on solving the boundary integral equations studied in [5].
Once M µ is known on the boundary, so is f µ ( · , k)| ∂Ω , which can be used to expand f µ as a power series outside Ω. Choose a point z 0 ∈ R 2 \ Ω. Set We next solve the truncated Beltrami equations with solutions represented in the form where ε(k)/k → 0 and ε(k)/k → 0 as k → ∞. Requiring To compute α (R) we first find solutions η 1 and η 2 to the equation For computational construction of the frequency-domain cgo solutions η j satisfying (3.7-3.9), see [3, Section 5.2] and [24].
Fix any nonzero k 0 ∈ C and choose any point z inside the unit disc. Denote We can now use the truncated transport matrix We know the approximate solutions u (R) j (z, k 0 ) for z ∈ Ω and one fixed k 0 . We use formulas u Finally we reconstruct the conductivity σ approximately as 3.2. The shortcut method. Fix R > 0 corresponding to a certain noise level δ in the EIT data. For all |k| < R, we compute the traces of M µ (·, k) and M −µ (·, k) by solving the boundary integral equation aforementioned in Section 3.1. The analyticity of M µ (·, k), M −µ (·, k) outside the unit disc allows to develop these factors as follows for |z| > 1.
In Section 5 of [6] it is proved that the scattering transform τ satisfies Notice that this formula is consistent with (2.4). Let τ δ R (k) be a numerical approximation to τ (k) such that τ δ R (k) = 0 for |k| > R. For all z ∈ Ω, we solve the D-bar equation This method is proven to work for C 2 conductivities σ by computing the approximation τ δ R (k) in a different way, namely the methods described in [43,45,40,31,30,33,34]. Nevertheless, there is no theoretical understanding of what happens when σ is nonsmooth and only essentially bounded.

3.3.
Periodization of the nonlinear inverse transform. The algorithm we use to solve equation (3.14) is presented in [35] and is a modification of the method introduced by Vainikko in [47] for the Lippmann-Schwinger equation related to the Helmholtz equation. Such algorithm is based in reducing the equation to a periodic one. Let us see how this reduction is done.
Definition 3.1. We say a function f defined on the plane is 2s-periodic if f (k) = f (k + 2s(j 1 + ij 2 )) for any k ∈ C, j 1 , j 2 ∈ Z. That is to say, we mean that f is 2s-periodic in each coordinate.
For any measurable subset A of R 2 we denote by χ A both the characteristic function of A and the multiplier operator with such function in k.
We write E 0 to denote the extension by zero outside the torus Q defined as for any f (z, ·) ∈ L 1 (Q).
Define a 2s-periodic approximation to the principal value 1/(πk) by setting β(k) = η(k)/(πk) for k ∈ Q and extending periodically: Define a periodic approximate solid Cauchy transform in k by where * denotes convolution on the torus. We shall use another periodization of the Cauchy transform given by where E denotes the periodic extension operator defined as Define the multiplication operator if k ∈ C and its periodization as Let us denote the complex conjugation operator by c, i.e. c(f ) = f . We periodize the complex conjugation operator c in the same way.
We shall need some boundedness properties of the so-called Cauchy transform. The Cauchy transform, we denote it by ∂ −1 , is defined as the following weakly singular integral operator: where k = k 1 + ik 2 and dk 1 dk 2 denotes the Lebesgue measure on R 2 . Theorem 4.3.8 in [2] states that for some constant C, 1 < q < 2 and q * = 2q/(2 − q).
For any 1 < q < ∞ and 1/q + 1/q = 1 the space L q (C) ∩ L q (C) with the norm · L q + · L q is a Banach space. Theorem 4.3.11 in [2] claims the following: Theorem 3.1. When 1 < q < 2 the Cauchy transform maps continuously the space L q (C)∩L q (C) into the space of continuous functions on the plane vanishing at infinity with the sup norm and satisfies the estimate It is easy to see that if 1 < q 1 < 2 < q 2 < ∞ but not necessarily q −1 1 + q −1 2 = 1 then and ∂ −1 φ is a continuous function assuming the right hand side bounded, where the implicit constant just depends on q 1 , q 2 . We conclude this section with the Theorem 3.2. Its proof uses the aforementioned properties of the Cauchy transform and is quite similar to the arguments presented in Sections 14.3.3 and 15.4.1 of the book [41]. Some of the techniques for Theorem 4.1 in [43] are involved.
Theorem 3.2. Assume c −1 ≤ σ ≤ c a.e. in the plane and σ ≡ 1 outside the unit disc. Let z ∈ C. There exists a unique 2s-periodic (in k) solution m R (z, k) to with m R (z, ·) ∈ L r (Q) for some 2 < r < ∞. Furthermore, the solutions of (3.15) and (3.18) agree on the disc of the k-plane with radius R: m R (z, k) = m R (z, k) for |k| < R.

Computational results
Once scattering data with big cutoff frequencies are generated, we test the nonlinear inversion procedure consisting of computing the approximate conductivity from the truncated scattering transform by solving The unique solution m R (z, k) to (4.1) with the normalization condition (4.2) can be obtained by solving (3.15). Finally, the approximate reconstruction σ R (z) is obtained by doing σ R (z) := m R (z, 0) 2 for |z| < 1.

Nonlinear Gibbs phenomenon in radial cases.
In this section we test the shortcut method for two radially symmetric conductivities as follows: Firstly, we solve numerically the Beltrami equation and obtain reliable scattering data over a mesh for a real interval of the form [0, R] with R > 0. Secondly, from such scattering data the approximate reconstruction is computed by solving numerically the equation (3.18) on a mesh for the real interval [0, 1]. It suffices to confine such computations to real positive numbers since the results on the whole disc can be obtained by simple rotation with respect to the origin by symmetry. The pictures for such reconstructions corresponding to the largest cutoff frequencies exhibit certain artifacts near the jump discontinuities of the actual conductivity quite similar to the Gibbs phenomenon for the truncated linear Fourier transform. 4.1.1. Example σ 1 . We consider the rotationally symmetric conductivity σ 1 on the plane defined as σ 1 (z) = 2 if |z| < 0.5 and σ 1 (z) = 1 if |z| ≥ 0.5. Figure 3 shows its representation.

Computing the scattering transform
The numerical evaluation of the AP scattering transform τ (k) was computed for this example conductivity σ 1 . The used algorithm generates for each fixed k cgo solutions for the Beltrami equation in the z-plane using the Huhtanen-Perämäki preconditioned solver ( [24]) by executing standard C-linear GMRES. To do so a z-grid is required. Such a z-grid is determined by a size parameter m z so that 2 mz × 2 mz equidistributed points are taken in the z-plane within the square [−s z , s z ) × [−s z , s z ) with s z = 2.1 and step h z = s z /2 mz−1 . Let us denote the numerical computation of τ with a z-grid size parameter m z by τ mz .
We know that for radial, real-valued conductivities the function t is rotationally symmetric too (Theorem 2.1). Therefore, we can restrict the numerical computation of t to positive values of the k-plane. The chosen k-grid was all the points in the interval [0, 150] with step h = 0.1. Thus, the scattering transform is computed farther away from the origin than any previous computation.
The vector τ mz corresponding to the aforementioned k-grid was computed for both m z = 10 and m z = 11. Next, the function t mz defined as t mz (k) = −4πi k τ mz (k), was evaluated on the same k-grid in [0, 150] for m z = 10,11 from the AP scattering transform using identity (2.10).
One can consider the approximate scattering transform τ 11 (with m z = 11) reliable since its accuracy is high in a double sense.
On one hand, Figure 4 shows the lack of significant differences between t 10 and t 11 for k close to the origin. As a matter of fact, using the notation  On the other hand, |Re( τ 11 )| is very small. By Theorem 2.1 and formula (2.10), for a rotationally symmetric, real-valued conductivity the scattering transform τ restricted to real numbers has to be pure-imaginary valued. We obtain a nonzero real part Re( τ mz ) for both m z values due to precision errors. In particular, |Re( τ 11 )| ≤ 3.7196 × 10 −9 .
Note that one can extend t 11 to the disc centered at the origin with radius 150 in the k-plane by simple rotation since t 11 (k 1 ) = t 11 (k 2 ) if |k 1 | = |k 2 |.

Solving the D-bar equation
We compute a number of reconstructions of σ 1 using the D-bar method truncating the approximate scattering transform τ 11 . This way the computation is optimized by restricting the degrees of freedom to the values of τ 11 (k) at grid points satisfying |k| < R, for some positive number R called the cutoff frequency.
The program constructs the set of evaluation points in the k-plane for which the D-bar equation is solved for every fixed z, as a k-grid of 2 m k × 2 m k points within the square [−s k , s k ) × [−s k , s k ) with s k = 2.3 × R and step h k = s k /2 m k −1 . Let σ 1 (z; m k , R) denote the shortcut reconstruction of σ 1 with size parameter m k and truncation radius R.
Again we can restrict the computation (of σ 1 (z; m k , R) in this case) to positive values in the z-domain and afterwards extend to the plane by rotation. We evaluated σ 1 (z; m k , R) for each point z in a grid of a certain number N x of equispaced points in the real interval [0,1].
The vector σ 1 (z; m k , R) was computed for R = 5, 10, 15, 20, 25, 30, 35, 40, 50,..., 150. For R = 5, 10 we took m k = 8, 9 and for the rest of R values we chose m k = 12. These are by far the largest nonlinear cut-off frequencies studied numerically so far. In Figure 5 the reader can see three of these reconstruction examples.
Since the shortcut reconstruction for σ 1 is real-valued, the obtained imaginary parts are due to precision errors. Such errors are displayed multiplied by 10 8 in the following tables: where error 1 = max |z|≤1 |Im( σ 1 (z; m k , R))|.

Computing the scattering transform
For conductivity σ 2 the scattering transform τ 11 (with m z = 11) turns out to be very accurate too. Indeed, t 10 and t 11 corresponding to σ 2 are quite similar for |k| small in view of Figure 7 and the following facts (see notation

4.2.
Comparison of the two methods in nonsymmetric cases. We apply both reconstruction methods to some discontinuous non-radial "checkerboardstyle" conductivities without assuming EIT data. We choose the conductivity examples σ 3 , σ 4 shown in Figure 9. Both take value 1 near the boundary of the unit disc. Notice that the contrasts of these examples σ 3 , σ 4 (defined as the difference between the maximum and the minimum) are 1.5 and 2.8, respectively. We denote such errors as follows: where σ denotes the actual conductivity and σ the approximate reconstruction.

Low-pass transport matrix method
In order to test the low-pass transport matrix method for each conductivity, we obtain the CGO solutions f µ (z 0 , k) , f −µ (z 0 , k) directly from the known conductivities without simulating EIT data. The point z 0 is outside the unit disc and we take the k-grid of 2 m k × 2 m k equispaced points in the square [−R, R) 2 with step h k = R/2 m k −1 and cut-off frequencies R ≤ 50. Next, transportation of f µ (z 0 , k) , f −µ (z 0 , k) to the unit disc and the final reconstruction σ (R) are performed following the steps described in Section 3.1. To this end, in addition to the k-grid, a grid in the z-variable of 2 mz × 2 mz equidistributed points in the square [−s z , s z ) 2 , with step h z = s z /2 mz−1 , is required.

Shortcut method
Concerning the shortcut method, firstly the approximate scattering transform, denote it by τ SC j for σ j (j = 3, 4), is computed directly from the known conductivity through the Beltrami equation solver using formula (2.4). The code evaluates the approximation over the points inside the disc of radius R (cut-off frequency) from a k-grid of 2 m k × 2 m k equispaced points in the square [−R, R] 2 . Another z-grid with 2 mz × 2 mz points in [−s z , s z ) 2 (s z > 2) is involved here.
To compute τ SC 3 we take R = 50, m k = 8, m z = 10, and for τ SC 4 we choose R = 50, m k = 7, m z = 7. Let us remark that this choice of R enables reconstructions corresponding to cut-off frequencies less than or equal to 50.
Secondly, Figures 10, 11 show plots of τ SC j and t SC j , where j = 3, 4 and t SC j is computed from τ SC j by formula (2.10). Finally, we evaluate the approximate reconstruction by solving the D-bar equation with the truncated scattering transform τ SC j as a coefficient. The reconstructions are computed over the points of the same kind of z-grid as above belonging to the unit disc with size parameter m z . The code computes the solution to the D-bar equation on every point in the z-grid independently and uses a number of k-grids with size parameter m k and cut-off frequency R.
These reconstructions were evaluated with m z = 7 and some R values with R ≤ 50. In addition, with respect to σ 3 we took m k = 9 and for σ 4 , m k = 7.
Notice that the case R = 50 was correctly computed for both conductivities σ 3 , σ 4 using the shortcut method, but for σ 4 the highest R-value through the low-pass transport matrix method was R = 20.       Figure 13: Comparison of reconstructions for σ 4 by the low-pass transport matrix method (left) and the shortcut method (right) using different cut-off frequencies R. The chosen point z 0 appears in black. The true conductivity is represented at the bottom row.

4.3.
Numerical evidence of the transport matrix efficiency. In addition, further numerical experiments were made aimed at comparing the actual complex geometric optics solution to the Beltrami equation and the transported solution f µ (z, k 0 ) to the unit disc from f µ (z 0 , k 0 ) and f −µ (z 0 , k 0 ) at certain z 0 ∈ C \ Ω. Figure 14 shows some pictures for f µ (z, k 0 ) and f µ (z, k 0 ) corresponding to the conductivity example σ 4 . Here |z| < 1, k 0 ≈ 1 and the pivot point outside the unit disc is z 0 ≈ (0, 1.266). To compute f µ the transport matrix was generated with a cut-off frequency R = 20 and we took the size parameters as follows: m k = 7, m z = 7, s z = 1.5.
The relative errors with sup and l 2 norms Real part Imaginary part Figure 14: Comparison of f µ (z, k 0 ) and f µ (z, k 0 ) for the conductivity σ 4 with R = 20. In the first row we show real and imaginary parts of both functions. In a) the real part of the actual solution f µ (z, k 0 ) is depicted on the left and the real part of the transported f µ (z, k 0 ) on the right (where the point z 0 is represented in black). In b) the same information is showed for the imaginary part of f µ (z, k 0 ) and f µ (z, k 0 ). The last row c) shows the difference f µ (z, k 0 ) − f µ (z, k 0 ). The real part is exhibited on the left and the imaginary part on the right. Again, the pivot point z 0 appears in black. All the pictures were generated under the same colormap using Matlab.

Conclusion
Regarding the radially symmetric examples studied in Section 4.1, the conclusion is clear: numerical evidence suggests that discontinuous conductivities can be reconstructed with the shortcut method more and more accurately when R tends to infinity.
The errors in the reconstructions are very similar to the Gibbs phenomenon observed with truncated linear Fourier transforms.
The numerical study of the nonsymmetric examples comparing both reconstruction types leads us to two observations: -For these discontinuous conductivity examples, the shortcut method generates considerably better reconstructions than the transport matrix method.
-The approximate solution computed by the transport matrix method is reliable on an area within the unit disc close enough to the selected z 0 point outside the unit disc.
Finally, since errors (4.4) are small and plots for both functions f µ (z, k 0 ), f µ (z, k 0 ) are similar in view of Figure 14, we have numerical evidence supporting the fact that the transport matrix method "transports well" and the worse results by the transport matrix method are explained by the final algebraic steps, including numerical differentiation in (3.12), just after the computation of the solutions u (R) 1 (z, k 0 ), u (R) 2 (z, k 0 ) in (3.11) through the transport matrix itself.

Appendix
This Appendix is aimed at presenting two final discussions as follows. On a hand, the noisy scattering transform referred to in Section 3.2 is compared with the free-noise scattering transform computed by solving the Beltrami equation and formula (2.4) for the checkerboard-style conductivity phantom σ 3 from Figure  9. On the other hand, a quantitative discussion on the precision required to the DN map to obtain the accuracy of the scattering transform on the disc of radius 50 is introduced.
Example of noisy scattering transforms For simplicity of notation, let us write τ p for the scattering data obtained on the disc centered at the origin of radius 6 by solving the boundary integral equation mentioned in Section 3.1 from a simulated Dirichlet-to-Neumann map, Λ p , with p % noise. Additionally, denote by τ the aforementioned free-noise scattering transform on the same disc (denoted by τ SC 3 in Section 4.2). The approximate Dirichlet-to-Neumann map Λ p is computed adding p% Gaussian noise to the FEM boundary voltages (for σ 3 ) as it is explained in Section 5.3 of [20].
The k-grid consists of 2 m k × 2 m k equispaced points in the square [−6, 6] 2 with m k = 7. The scattering data are computed on the 12, 851 points of such grid belonging to the disc of radius 6. For the grid in the z-variable used to compute τ , 2 mz × 2 mz equidistributed points are taken in the square [−s z , s z ) 2 , with step h z = s z /2 mz−1 and s z = 2.3.
For different choices of p the following absolute error is computed for r varying over a partition of the interval (0, 6]: where χ D(0,r) denotes the characteristic function of the disc centered at the origin and radius r.
In Figure 15 the profiles of E p (r) for p = 0.1%, 0.5%, 1%, 5% are shown in black, blue, green and red colour, respectively. The maximum value shown on the vertical axis is 1. Accuracy of the DN map to compute t(k) on D(0, 50) How many correct significant digits would one need in the DN map for computing the scattering transform t(k) for |k| < 50? It is well-known [1] that there is a logarithmic relation between measurement accuracy and details in conductivity. In other words, EIT is an exponentially ill-posed inverse problem. In terms of scattering transforms, a fixed noise level p% in the EIT measurement leads to, roughly speaking, the computation of t(k) to be stable and accurate in a disc |k| < R(p) and unstable outside that disc. See [34, Figure 2]. The exponential ill-posedness shows up as a logarithmic dependence of R on p, see [34, Figure  3]. Improving the measurement accuracy so that the DN matrix entries have one more significant digit increases the stability radius from R to R + Γ with Γ > 0 a fixed real number.
Let us make a crude quantitative estimation based on [34, Figure 3]. Two significant digits in the DN map corresponds to R = 3 and six significant digits to R = 7. A quick computation shows that achieving the stability radius R = 50 would require 49 correct significant digits in the elements of the DN matrix (and, actually, using a larger matrix with more oscillatory basis functions involved as well).