0 in Ω, (1.1) where Ω is a bounded open set in R2 with smooth boundary ∂Ω and the matrix γ(x

ANISOTROPIC EQUATIONS AND APPLICATIONS

HIDEKI TAKUWA, GUNTHER UHLMANN, AND JENN-NAN WANG

Abstract. In this article, we construct complex geometrical op- tics solutions with general phase functions for the second order elliptic equation in two dimensions. We then use these special so- lutions to treat the inverse problem of reconstructing embedded inclusions by boundary measurements.

1. Introduction

In previous works [21] and [22], special complex geometrical optics (CGO) solutions for certain isotropic systems in the plane were con- structed and their applications to the object identification problem were also demonstrated theoretically and numerically. Those systems include the conductivity equation, the elasticity equations, and the Stokes equations, all with isotropic inhomogeneous coefficients. In this paper we will extend the method in [21] and [22] to scalar equations with anisotropic inhomogeneous coefficients in the plane. We shall fo- cus on the conductivity equation:

L_γu =:

2

X

i,j=1

∂_x^j(γ^ij(x)∂_xⁱu) = 0 in Ω, (1.1) where Ω is a bounded open set in R² with smooth boundary ∂Ω and the matrix γ(x) = [γ^ij(x)] is symmetric and there exists C₀ > 0 such that

γ(x) ≥ C₀I for all x ∈ Ω.

The precise regularity of γ will be specified later on. Our method can treat the equation with L_γ as the leading order without essential modification. The main aim here is to construct solutions to (1.1) with

The first author was partially supported by Grant-in-Aid for Young Scientists (B) 18740073, Japan.

The second author was partially supported by NSF and a Walker Family En- dowed Professorship.

The third author was supported in part by the National Science Council of Taiwan (NSC 96-2628-M-002-009).

1

general complex phases. These types of solutions are very useful in the reconstruction of objects embedded in Ω by boundary measurements.

The key step in treating (1.1) is to transform it to an isotropic equa- tion using isothermal coordinates, which are closely related to quasi- conformal mappings [1]. This idea is crucial in studying the Calder´on problem in two dimensions (see, for example, [2] and [20]). After reduc- ing the anisotropic equation (1.1) to an isotropic one, we can apply the method in [21] to construct CGO solutions for the isotropic equation and then transform the solutions back to the original coordinates.

A useful application of these CGO solutions is to study the object identification problem by boundary measurements when the known background medium is anisotropic. Using the complex geometrical op- tics solutions in the reconstruction of embedded objects from boundary measurements was first introduced by Ikehata [11], [12]. He called this method the enclosure method. Our method shares the same spirit as the enclosure method. For brevity, we will not give a detailed Ikehata’s enclosure method here. We refer the reader to a nice survey article by Ikehata [13] for more details and related applications.

In the implementation of the enclosure type method, how much in- formation of the embedded object one can reconstruct depends on the phase functions of the CGO solutions used in the method. For exam- ple, one can reconstruct the convex hull of the object with Calder´on type solutions [11], [12]. With complex spherical waves, one may be able to reconstruct some nonconvex parts of the object [8], [19]. In the two-dimensional case, we are able to get much more information of the object by using Mittag-Leffler functions [14], [15] or the spe- cial solutions constructed in [21]. Note that Mittag-Leffler functions can be only applied to the case where the background equation is the Laplacian. Like the results in [21], the method of this paper works for any general second order elliptic equations with coefficients having appropriate finite smoothness.

In this paper, we provide theoretical grounds for our method. In order to find the coordinates transformation, we need to solve the Bel- trami equation. In Section 4, we will describe a numerical scheme to solve the Beltrami equation. This algorithm is a combination of the fast algorithm developed by Daripa et al. [3], [4], [5] (also see [6]) and the NUFFT (nonuniform fast Fourier transform) [7]. The rest of the paper is organized as follows. In Section 2, we will show how to construct complex geometrical optics solutions with general phases for (1.1). The key step is to reduce the equation to the isotropic equation with the

help of a quasiconformal mapping. In Section 3, we provide theoret- ical backgrounds of reconstructing the embedded object by boundary measurements with these special solutions.

2. Construction of CGO solutions

In this section we will give a framework of constructing complex geo- metrical optics solutions for L_γ. We first recall a fundamental fact: let F : Ω → eΩ, y = F (x), be a C¹ bijective map (orientation-preserving), then u(x) solves (1.1) if and only if v(y) = u ◦ F⁻¹(y) solves

Leγv = 0 in Ω,e (2.1)

where eγ = F∗γ defined by (F∗γ)^k`(y) = 1

det(DF )

2

X

i,j=1

∂_xⁱF^k∂_x^jF^`γ^ij

x=F⁻¹(y).

Our aim here is to find a suitable F so that F∗γ is isotropic. In the plane, the existence of such map F is well-known and it is a quasicon- formal map. In the following, we will review some essential properties we need from the theory of quasiconformal mapping. For more details, we refer to [1], [16], and [18].

Let us denote R² as the complex plane C and z = x¹+ ıx² from now on. Assume that γ ∈ C¹(Ω). We now extend γ to C, still denoted by γ, such that γ ∈ C¹(C) and γ = I for |z| > R for some R > 0.

Theorem 2.1. There exists a C¹ bijective map F : C → C such that F∗γ = (detγ ◦ F⁻¹)^1/2I. (2.2) Proof. The proof of this theorem can be found in [20] for C³ or C^1,1 conductivities and in [2] for L^∞ conductivities. For the purpose of numerical computation in the following sections, we sketch the proof here. The relation (2.2) holds if and only if F is a quasiconformal map with complex dilatation

µ_γ = γ²²− γ¹¹− 2iγ¹² γ¹¹+ γ²²+ 2√

detγ, i.e.,

∂F = µ¯ _γ∂F, (2.3)

where

∂ =¯ 1

2(∂_x¹ + i∂_x²) and ∂ = 1

2(∂_x¹ − i∂_x²).

Equation (2.3) is the well-known Beltrami equation. Note that here µγ

is supported in a disc of radius R and kµ_γk_L^∞ ≤ K < 1.

Following Ahlfors’ approach [1, Chapter V], there are two integral operators which play key roles in solving (2.3):

P f (ξ) = ı 2π

Z

C

1 z − ξ − 1

z
f (z)d¯z ∧ dz and

T f (ξ) = lim

ε→0

ı 2π

Z

|z−ξ|>ε

f (z)

(z − ξ)²d¯z ∧ dz.

The existence of F to (2.3) can be achieved in two steps [1, page 92].

Firstly, we solve the integral equation

h = T (µ_γh) + T µ_γ. (2.4) Let h be a solution of (2.4) then

F = P (µ_γ(h + 1)) + z (2.5) solves (2.3). Since ∂µ_γ ∈ L^p for any p > 2, it follows from [1, Lemma 3, page 94] that F ∈ C¹ and is bijective.

Before ending the proof, we say a few words on the existence of (2.4).

The celebrated Calder´on-Zygmund inequality indicates that for p > 1 kT f k_L^p ≤ C_pkf k_L^p

with the constant C_p only depending on p. Also, one has C_p → 1 as p → 2. Now by taking KC_p < 1, the existence of (2.4) can be obtained using the Neumann series. On the other hand, the operator is well-defined on L^p functions with p > 2.  Now we discuss how to construct complex geometrical optics solu- tions for L_γ. Suppose that the original conductivity function γ, defined on Ω, was suitably extended and a quasiconformal map F described in Theorem 2.1 was found. Let us define ˜F = F |_Ω. Then ˜F : Ω → ˜Ω with Ω = ˜˜ F (Ω) is a C¹ diffeomorphism and

˜

γ(y) = ˜F∗γ = (detγ ◦ ˜F⁻¹)^1/2(y) ∈ C¹( ˜Ω).

As before, we know that v(y) solves 1

˜

γL_˜γv = ∆v + ∇˜γ

˜

γ · ∇v = 0 in Ω˜ (2.6) if and only if u(x) = v ◦ F (x) solves

L_γu = 0 in Ω.

Therefore, it suffices to construct CGO solutions for (2.6). We recall that some CGO solutions with general phase functions for (2.6) have been established in [21]. We now demonstrate how to adapt the ar- guments in [21] to our setting here. Pick a y₀ ∈ C. Without loss of generality, we take y₀ = 0. For N ∈ N, let

φ˜_N(y) = c_N(y − y₀)^N = c_Ny^N = ˜ϕ_N(y) + ı ˜ψ_N(y),

where y = y1 + ıy2 ∈ C, ˜ϕN = Re(cNy^N), ˜ψN = Im(cNy^N), cN ∈ C with |cN| = 1. Some level curves of ˜ϕN were shown below (also see [21]).

As in [21], let Γ_N be the open cone with vertex at 0 and opening angle π/N in which ˜ϕ_N > 0 in Γ_N and ˜ϕ_N = 0 on the boundary of Γ_N. We assume that Γ_N ∩ ˜Ω 6= ∅. According to [21], we can construct solutions of the form

˜

v_N,h= e^φ˜N/h(˜a_N(y) + ˜r_N(y, h)) to (2.6), where ˜a_N(y) never vanishes in Γ_N ∩ ˜Ω and

k∂_y^βr˜Nk_L2(ΓN∩ ˜Ω) ≤ Ch^1−|β| ∀ |β| ≤ 1.

We now extend ˜v_N,h defined in Γ_N ∩ ˜Ω to the whole ˜Ω by introducing a suitable cut-off function.

For s > 0, let ˜`_s = {y ∈ Γ_N : ˜ϕ_N = s⁻¹}. This is the level curve of

˜

ϕ_N in Γ_N. Let 0 < t < t₀ such that

(∪_s∈(0,t)`˜_s) ∩ ˜Ω 6= ∅

and choose a small ε > 0. Define a cut-off function χ_N,t(y) ∈ C^∞(R²) so that χ_N,t(y) = 1 for y ∈ (∪s∈(0,t+ε/2)`˜_s) ∩ ˜Ω and is zero for y ∈

¯˜

Ω \ (∪_s∈(0,t+ε)`˜<sub>s</sub>). We want to remark that ∪s∈(0,t+ε/2)`˜_s ⊂ ∪_s∈(0,t+ε)`˜_s. Now we define

˜

v_N,t,h(y) = χ_N,te^−t−1/h˜v_N,h= χ_N,te^{( ˜ϕN−t−1+ı ˜ψN)/h}(˜a_N + ˜r_N) for y ∈ (∪_s∈(0,t+ε)`˜_s) ∩ ˜Ω. So ˜v_N,t,h can be regarded as a function in ˜Ω which is zero outside of ΓN ∩ ˜Ω. Note that ˜vN,t,h does not solve (2.6) in ˜Ω. We thus take ˜f_N,t,h = ˜v_N,t,h|_{∂ ˜Ω} and consider the boundary value problem:

(∆˜u_N,t,h+^∇˜_˜γ^γ · ∇˜u_N,t,h= 0 in Ω˜

˜

u_N,t,h= ˜f_N,t,h on ∂ ˜Ω. (2.7)

Clearly (2.7) is uniquely solvable. Also it was shown in [21] that there exist ˜C > 0 and ˜ε > 0 such that

k˜u_N,t,h− ˜v_N,t,hk_H1( ˜Ω) ≤ ˜Ce^−˜ε/h

for h  1.

Next we define

u_N,t,h(x) = (˜u_N,t,h◦ F )(x), fN,t,h= ˜fN,t,h◦ F |∂Ω, and

v_N,t,h(x) = (˜v_N,t,h◦ F )(x) = (χ_N,t◦ F )e^{(ϕN−t−1+iψN)/h}(a_N + r_N), where ϕ_N = ˜ϕ_N ◦ F , ψ_N = ˜ψ_N ◦ F , a_N = ˜a_N ◦ F , and r_N = ˜r_N ◦ F . Then, aN is never zero in F⁻¹(ΓN ∩ ˜Ω) and

k∂_x^βr_Nk_L2(F⁻¹(ΓN∩ ˜Ω)) ≤ Ch^1−|β| ∀ |β| ≤ 1.

Moreover, we have that u_N,t,h satisfies (L_γu_N,t,h= 0 in Ω

u_N,t,h= f_N,t,h on ∂Ω and

ku_N,t,h− v_N,t,hk_H¹_(Ω) ≤ Ce^−ε0/h (2.8) for some positive constants C and ε⁰. The solution uN,t,h(x) is a CGO solution of (1.1) and is approximated by v_N,t,h with exponentially small errors. We shall pay attention to the level curves of ϕ_N = ˜ϕ_N ◦ F defined in F⁻¹(Γ_N), which are key to our reconstruction method for the object identification problem.

3. Applications to inverse problems

In this section we demonstrate how to use CGO solutions constructed previously in the object identification problem. To simplify our pre- sentation, we will only discuss the case of identifying inclusions inside of the domain Ω filled with known conductivity. This inverse problem has been extensively studied both theoretically and numerically. We refer to [8] for related references. Let D be an open bounded domain with C¹ boundary such that ¯D ⊂ Ω and Ω \ ¯D is connected. Assume that γ₀(x) ∈ C¹( ¯Ω) is a symmetric matrix satisfying

γ₀(x) ≥ δI

with some δ > 0 for all x ∈ ¯Ω. The conductivity γ(x) is a perturbation of γ₀ described by γ(x) = γ₀ + χ_Dγ₁, where χ_D is the characteristic function of D and γ₁ ∈ C( ¯D) is a positive semi-definite matrix, i.e,

γ₁ ≥ 0 on D.¯ (3.1)

Then we have γ(x) ≥ δI almost everywhere in Ω. Let v be the solution of

(L_γv = 0 in Ω,

v = f on ∂Ω. (3.2)

The meaning of the solution to (3.2) is understood in the following way.

Define

[w]_∂D = tr⁺w − tr⁻w

the jump of the function across ∂D, where tr⁺ and tr⁻ denote re- spectively the trace of w on ∂D from inside and outside of D. For f ∈ H^3/2(∂Ω), we define

V_f = {w ∈ H²(D) ⊕ H²(Ω \ ¯D) : w|_∂Ω= f, [w]_∂D = 0, [γ∂_νw]_∂D = 0}, where ν is unit normal of ∂D and

γ∂νw =

2

X

j,k=1

γ^jkνk∂_x^jw.

We say that v is the solution of (3.2) if v ∈ V_f and L_γv = 0 in D and Ω \ ¯D. The Dirichlet-to-Neumann map is given as

Λ_D : f → γ∂_nv|_∂Ω=

2

X

j,k=1

γ₀^jkn_k∂_x^jv|_∂Ω,

where n is the unit outer normal of ∂Ω. The inverse problem is to determine the inclusion D from Λ_D. Here we are interested in the reconstruction question.

To begin, we would like to derive two important integral inequalities.

Let u(x) ∈ H²(Ω) be the solution of

(L_γ0u = 0 in Ω,

u = f on ∂Ω (3.3)

for f ∈ H^3/2(∂Ω). The Dirichlet-to-Neumann map corresponding to (3.3) is denoted by

Λ₀ : f → γ₀∂_nu|_∂Ω. Then we can prove

Lemma 3.1.

h(Λ_D− Λ₀) ¯f , f i_∂Ω ≤ hγ₁∇u, ∇¯ui_D (3.4) and

h(γ₀⁻¹− γ⁻¹)γ₀∇u, γ₀∇¯ui_D ≤ h(Λ_D− Λ₀) ¯f , f i_∂Ω, (3.5)

where

hA∇u, ∇¯ui_D = Z

D

(∇¯u)^TA∇udx for any matrix A.

Proof. When γ0 and γ1 are isotropic, derivations of (3.4) and (3.5) can be found in [8] (or in [9]). Here we adapt their arguments to treat anisotropic conductivities. Let v and u be solutions of (3.2) and (3.3) with same Dirichlet condition f , respectively. From Green’s formula, we have that

hγ∇v, ∇¯vi_Ω = hγ∇v, ∇¯ui_Ω (3.6) and

hγ₀∇u, ∇¯ui_Ω = hγ₀∇u, ∇¯vi_Ω. (3.7) Next we observe that

hΛ0f , f i¯ ∂Ω= hγ0∇¯u, ∇uiΩ = hγ0∇u, ∇¯uiΩ (3.8) and

hΛ_Df , f i¯ _∂Ω= hγ∇¯v, ∇ui_Ω = hγ∇u, ∇¯vi_Ω. (3.9) Using (3.6), (3.8), (3.9), and some simple computations, we get that

h(Λ_D − Λ₀) ¯f , f i_∂Ω

≤h(Λ_D − Λ₀) ¯f , f i_∂Ω+ hγ∇(v − u), ∇(¯v − ¯u)i_Ω

=hγ∇u, ∇¯vi_Ω− hγ₀∇u, ∇¯ui_Ω+ hγ∇(v − u), ∇(¯v − ¯u)i_Ω

=h(γ − γ₀)∇u, ∇¯ui_Ω

=hγ1∇u, ∇¯uiD, which is exactly (3.4).

By similar computations, we can show that

h(Λ_D − Λ₀) ¯f , f i_∂Ω= hγ₀∇(v − u), ∇(¯v − ¯u)i_Ω+ h(γ − γ₀)∇v, ∇¯vi_Ω. (3.10) On the other hand, it is not hard to check that

hγ₀∇(v − u), ∇(¯v − ¯u)i_Ω+ h(γ − γ₀)∇v, ∇¯vi_Ω

=hγ(∇v − γ⁻¹γ0∇u), (∇¯v − γ⁻¹γ0∇¯u)iΩ+ hγ0(γ₀⁻¹− γ⁻¹)γ0∇u, ∇¯uiΩ. (3.11) Using (3.10), (3.11), and the positive-definiteness of γ immediately

leads to (3.5). 

Now we are ready to consider the reconstruction of D from the mea- surements ΛD. To make sure that the medium has jumps across ∂D, we assume that for any q ∈ ∂D, there exist δq > 0 and εq > 0 such that

γ₁(x) ≥ ε_qI for all x ∈ D ∩ B(q, δ_q). (3.12)

To set up the measurements, we choose the following Dirichlet condi- tion

f_N,t,h= u_N,t,h|_∂Ω= v_N,t,h|_∂Ω= (χ_N,t◦ F )e^{(ϕN−t−1+ıψN)/h}(a_N + r_N)|_∂Ω. Note that the Dirichlet condition f_N,t,h is localized on ∂Ω. With f_N,t,h, we now solve the boundary value problem

(L_γu_N,t,h= 0 in Ω, u_N,t,h= f_N,t,h on ∂Ω and then evaluate the quadratic term

E(N, t, h) := h(Λ_D − Λ₀) ¯f_N,t,h, f_N,t,hi_∂Ω.

It is clear that E(N, t, h) is determined entirely by boundary measure- ments. Before proving the main reconstruction result, we would like to remark that

{x ∈ F⁻¹(Γ_N) : ϕ_N > t⁻¹₁ } ⊂ {x ∈ F⁻¹(Γ_N) : ϕ_N > t⁻¹₂ } for t₁ < t₂. Theorem 3.2. Let t > 0 be given. Then we have:

(i) if ¯D ∩ {x ∈ F⁻¹(Γ_N) : ϕ_N > t⁻¹} = ∅ then there exist C₁ > 0, ε₁ > 0, and h₁ > 0 such that E(N, t, h) ≤ C₁e^−ε1/h for all h ≥ h₁; (ii) if D ∩ {x ∈ F⁻¹(Γ_N) : ϕ_N > t⁻¹} 6= ∅ then there exist C₂ > 0, ε₂ > 0, and h₂ > 0 such that E(N, t, h) ≥ C₂e^ε2/h for all h ≥ h₂. Proof. For statement (i), we use integral inequality (3.4) and (2.8) to obtain that

E(N, t, h) ≤ hγ₁∇u_N,t,h, ∇¯u_N,t,hi_D ≤ C Z

D

|∇u_N,t,h|²dx

≤ C Z

D

|∇v_N,t,h|²dx + Ce^−ε0/h.

(3.13)

Using the structure of v_N,t,h and (3.13) immediately leads to (i).

To prove (ii), we need to investigate the matrix γ₀⁻¹− γ⁻¹ a little bit more. For x ∈ D, we see that

γ₀⁻¹− γ⁻¹ = γ₀⁻¹− (γ₀+ γ₁)⁻¹

= (γ₀+ γ₁)⁻¹{(γ₀+ γ₁)γ₀⁻¹(γ₀ + γ₁) − (γ₀+ γ₁)}(γ₀+ γ₁)⁻¹

= (γ₀+ γ₁)⁻¹(γ₁+ γ₁γ₀⁻¹γ₁)(γ₀+ γ₁)⁻¹

(see similar computations in [10] for anisotropic elastic body). There- fore, γ₀⁻¹− γ⁻¹ is positive semi-definite and from (3.12) we have

(γ₀⁻¹− γ₁⁻¹)(x) ≥ ˜ε_qI ∀ x ∈ D ∩ B(q, δ_q) (3.14)

for some ˜εq > 0. Now assume that

D ∩ B(q, δ_q) ⊂ D ∩ {x ∈ F⁻¹(Γ_N) : ϕ_N > t⁻¹}

for some q ∈ ∂D and δ_q > 0. Using (3.14), (3.5), and (2.8), we can compute

E(N, t, h) ≥ h(γ₀⁻¹− γ⁻¹)γ₀∇u_N,t,h, γ₀∇¯u_N,t,hi_D

≥ C ˜ε_q Z

D∩B(q,δq)

|∇v_N,t,h|²dx − Ce^−ε0/h. (3.15) Combining the exponentially growing behavior of v_N,t,h and (3.15)

yields (ii). 

With the help of Theorem 3.2, we can determine whether the level curve {x ∈ F⁻¹(ΓN) : ϕN = t⁻¹} intersects D or not by looking into the behavior of E(N, t, h) as h  1. In practice, we can start with a smaller t and move the level curve deeper into the domain Ω by choosing a larger t. In addition, we are able to determine more parts of ∂D by increasing N . We refer to [21] for precise arguments.

4. Algorithm of constructing the quasiconformal map In the study of our reconstruction method mentioned above, the forward problem is solved with CGO solutions as the Dirichlet bound- ary condition. Thus, comparing with the reconstruction algorithm in [21] for the isotropic medium case, the only difference is the bound- ary condition used in the forward problem (3.2). In other words, the numerical scheme used in [21] can be adapted to the problem studied here by merely changing the boundary condition. Therefore, we do not intend to repeat computational results in the paper. Nonetheless, since the CGO solutions constructed here involve the quasiconformal map F , we want to discuss how to find this map F numerically. We will leave the implementation of this algorithm to interested readers.

A fast algorithm using FFT for solving the Beltrami equation has been proposed by Daripa et al. [3], [4], [5]. Based on Daripa’s al- gorithm, a slightly modified version was given in [6]. Here we shall follow the description of Daripa’s algorithm given in [6]. Since it is often desirable to refine the grid in the part of boundary where the measurements are taken, we shall replace FFT in Darpia’s algorithm by NUFFT described in [7].

As outlined in the proof of Theorem 2.1, a solution F to the Beltrami equation (2.3) can be constructed by solving integral equations:

h = T (µ_γh) + T (µ_γ) (4.1)

and

F = P (µγ(h + 1)) + z. (4.2) Here the anisotropic conductivity function γ(x) has been extended to C with γ ∈ C¹ and γ = I for x ∈ B(0, R)^c. Viewing the definition of the complex dilatation µ_γ, we get µ_γ ∈ C¹ and µ_γ(x) = 0 for all x ∈ B(0, R)^c.

As shown in the proof of Theorem 2.1, a solution to (4.1) can be constructed by a fix-point type iteration:

hⁿ⁺¹ = T (µ_γhⁿ) + T (µ_γ) (4.3) and hⁿ→ h^∗ in some L^p norm with p > 2. Let g be H¨older continuous and supported in B(0, R). For such g, it is known that T (g) exists in the sense of the Cauchy principal value [1]. Let us express g and T (g) in Fourier series

g(re^ıθ) =

∞

X

−∞

gk(r)e^ıkθ and

T (g)(re^ıθ) =

∞

X

−∞

t_k(r)e^ıkθ.

It was proved in [5] that if T (g) exists as the Cauchy principal value, then the Fourier coefficients of T (g) are given by

t0(0) = −2 lim

ε→0

Z R ε

g₂(ρ)

ρ dρ, and tk(0) = 0, for k 6= 0, (4.4) t_k(r) = A_k

Z r 0

r^k

ρ^k+1g_k+2(ρ)dρ + B_k Z R

r

r^k

ρ^k+1g_k+2(ρ)dρ + g_k+2(r), (4.5) where

A_k =

(0, k ≥ 0,

2(k + 1), k < 0, and B_k =

(−2(k + 1), k ≥ 0, 0, k < 0

(see also [4]). Note that g_k(0) = g(0) if k = 0 and g_k(0) = 0 if k 6= 0.

Having found a solution to (4.1), we can find a solution of the Bel- trami equation by solving (4.2). Let us define

P (˜˜ g)(z) = ı 2π

Z

C

˜ g(z⁰)

z⁰− zd¯z⁰∧ dz⁰, then

P (˜g)(z) = ˜P (˜g)(z) − ˜P (˜g)(0).

An effective algorithm for evaluating the integral transform ˜P was pro- posed in [3]. Precisely, let ˜g and ˜P (˜g) be expressed as

˜

g(re^ıθ) =

∞

X

−∞

˜

g_k(r)e^ıkθ and P (˜˜ g)(re^ıθ) =

∞

X

−∞

p_k(r)e^ıkθ, then

p_k(r) = (2Rr

0(_ρ^r)^kg˜_k+1(ρ)dρ, k < 0,

−2RR

r (^r_ρ)^k˜gk+1(ρ)dρ, k ≥ 0 (4.6) (see [3]). In view of the formulae (4.5) and (4.6), one can easily derive useful relations for t_k(r) and p_k(r) at any pair of nonnegative r and r⁰. Define

d_k(r, r⁰) = 2(k + 1) Z r⁰

r

1 ρ(r⁰

ρ)^kg_k+2(ρ)dρ and

e_k(r, r⁰) = 2 Z r⁰

r

(r

ρ)^k˜g_k+1(ρ)dρ, then we can derive

t_k(r) = (r

r⁰)^k(t_k(r⁰) − g_k+2(r⁰) − d_k(r, r⁰)) + g_k+2(r) (4.7) and

pk(r) = (r

r⁰)^kpk(r⁰) − ek(r, r⁰) (4.8) (see [3], [6]). Armed with recurrence formulae (4.7), (4.8), we can evaluate the Fourier coefficients t_k(r) and p_k(r) more efficiently.

We now describe how to implement the above algorithm in actual computation. In view of (4.2), since supp (µ_γ) ⊂ B(0, R), we only need to construct h in B(0, R). We discretize B(0, R) with a circular M1× M2 grid. In other words, each grid point is represented by (ri, θj), where 0 < r_i ≤ R and θ_j ∈ [0, 2π], for 1 ≤ i ≤ M₁, 0 ≤ j ≤ M₂ − 1.

Note that {θ_j} are not necessarily equi-distributed on [0, 2π]. We may take r₀ = 0 and r_M1 = R. To initiate the iteration (4.3), we can take h⁰ = 0. Then we have h¹ = T (µ_γ). Since in each iteration we need to compute T (µ_γ) and T (µ_γhⁿ), it is helpful to see how to calculate T (g) for a given function g. Suppose that the value of g at each grid point (r_i, θ_j) is given, say g(r_ie^ıθj) = g_i,j. Then implementing the type-1 NUFFT [7], one can compute

g_k(r_i) = 1 M₂

M2−1

X

j=0

g_i,je^−ıkθj (4.9)

for −^N₂ ≤ k ≤ ^N₂ − 1 and 1 ≤ i ≤ M₁. To find the Fourier coefficients of T (g), we start from r_M1 = R. Using (4.5), we see that

t_k(r_M1) = A_k Z r_M1

0

r_M^k₁

ρ^k+1g_k+2(ρ)dρ + g_k+2(r_M1). (4.10) We then approximate the integral in (4.10) with available values of g_k(r) at r_i (0 ≤ i ≤ M₁) using a suitable approximation formula.

Having determined t_k(r_M1), we then determine t_k(r_i) for i ≤ M₁− 1 by the recurrence formula (4.7) with r⁰ = r_M1. Applying (4.3), we can determine the discrete Fourier transform of hⁿ⁺¹. Let h^∗ be an approximate solution of (4.3) with Fourier coefficients h^∗_i,k for 0 ≤ i ≤ M₁ and −^N₂ ≤ k ≤ ^N₂ − 1.

Our next step is to construct a quasiconformal map F from (4.2).

To this end, we need to evaluate µ_γ(h^∗ + 1) =: ˜g. We do so by using the point-wise multiplication over grid points. To get the values of h^∗ at grid points, we simply perform type-2 NUFFT [7] with coefficients h^∗_i,k, i.e.,

h^∗(r_ie^ıθj) =

N 2−1

X

k=−^N₂

h^∗_i,ke^ıkθj.

Then applying the type-1 NUFFT again, we can compute ˜g_k(r_i). By (4.6), we see that for r_M1 = R

p_k(r_M1) =

(2RR

0 (^R_ρ)^kg˜_k+1(ρ)dρ, k < 0, 0, k ≥ 0.

To find p_k(r_i) for i ≤ M₁− 1, we simply apply the recurrence formula (4.8). Having found Fourier coefficients p_k(r_i), we obtain the values of P (˜˜ g) at grid points by using the type-2 NUFFT. Finally, we can get

P (˜g)(r_ie^ıθj) = ˜P (˜g)(r_ie^ıθj) − ˜P (˜g)(0) and

F (rie^ıθj) = P (˜g)(rie^ıθj) + rie^ıθj, where

P (˜˜ g)(0) = p₀(0) = −2 Z R

0

˜

g₁(ρ)dρ.

To get F at other points in B(0, R), we can use appropriate interpola- tion or extrapolation techniques.

References

[1] L.V. Ahlfors, Lectures on Quasiconformal Mappings, Van Norstrand, Prince- ton, NJ, 1966. Reprinted by Wadsworth Inc. Belmont, CA, 1987.

[2] K. Astala, L. P¨aiv¨arinta, and M. Lassas, Calder´on’s inverse problem for anisotropic conductivity in the plane, Comm. PDE, 30 (2005), 207-224.

[3] P. Daripa, A fast algorithm to solve non-homogeneous Cauchy-Riemann equa- tions in the complex plane, SIAM J. Sci. Statist. Comput., 13 (1992), 1418- 1832.

[4] P. Daripa, A fast algorithm to solve the Beltrami equation with applications to quasiconformal mappings, J. Comp. Phys., 106 (1993), 355-365.

[5] P. Daripa and D. Mashat, Singular integral transforms and fast numerical algorithms, Numer. Algor., 18 (1998), 133-157.

[6] D. Gaydashev and D. Khmelev, On numerical algorithms for the solution of a Beltrami equation, SIAM J. Numer. Anal., 46 (2008), 2238-2253.

[7] L. Greengard and J.-Y. Lee, Accelerating the nonuniform fast Fourier trans- form, SIAM Review, 46 (2004), 443-454.

[8] T. Ide, H. Isozaki, S. Nakata, S. Siltanen, and G. Uhlmann, Probing for electrical inclusions with complex spherical waves, to appear in Comm. Pure Appl. Math.

[9] M. Ikehata, Identification of the shape of the inclusion having essentially bounded conductivity, J. Inv. Ill-Posed Problems, 7 (1999), 533-540.

[10] M. Ikehata, G. Nakamura, and K. Tanuma, Identification of the shape of the inclusion in the anisotropic body, Applicable Analysis, 72 (1999), 17-26.

[11] M. Ikehata, Enclosing a polygonal cavity in a two-dimensional bounded do- main from Cauchy data, Inverse Problems, 15 (1999), 1231-1241.

[12] M. Ikehata, How to draw a picture of an unknown inclusion from bound- ary measuremets. Two mathematical inversion algorithms, J. Inv. Ill-Posed Problems, 7 (1999), 255-271.

[13] M. Ikehata, The enclosure method and its applications, Analytic extension formulas and their applications (Fukuoka, 1999/Kyoto, 2000), 87-103, Int.

Soc. Anal. Appl. Comput., 9, Kluwer Acad. Publ., Dordrecht, 2001.

[14] M. Ikehata, Mittag-Leffler’s function and extracting from Cauchy data, In- verse problems and spectral theory, 41-52, Contemp. Math., 348, Amer.

Math. Soc., Providence, RI, 2004.

[15] M. Ikehata and S. Siltanen, Electrical impedance tomography and Mittag- Leffler’s function, Inverse Problems, 20 (2004), 1325-1348.

[16] T. Iwaniec and G. Martin, Geometric Function Theory and Non-linear Anal- ysis, Oxford University Press Inc., New York, 2001.

[17] K. Knudsen, J.L. Mueller, and S. Siltanen, Numerical solution method for the dbar-equation in the plane, Journal of Computational Physics, 198 (2004), 500-517.

[18] O. Lehto and K.I. Virtanen, Quasiconformal Mapping in the Plane, 2nd ed., Springer-Verlag, New York, 1973.

[19] G. Nakamura and K. Yoshida, Identification of a non-convex obstacle for acoustical scattering, J. Inv. Ill-Posed Problems, 15 (2007), 1-14.

[20] J. Sylvester, An anisotropic inverse boundary value problem. Comm. Pure Appl. Math., 43 (1990), 201-232.

[21] G. Uhlmann and J.N. Wang, Reconstructing discontinuities using complex geometrical optics solutions, SIAM J. Appl. Math., 68 (2008), 1026-1044.

[22] G. Uhlmann, J.N. Wang, and C.T. Wu Reconstruction of inclusions in an elastic body, preprint.

Faculty of Science and Engineering, Doshisha University, Kyotan- abe City, 610-0394 Japan

E-mail address: htakuwa@mail.doshisha.ac.jp

Department of Mathematics, University of Washington, Box 354305, Seattle, WA 98195-4350, USA.

E-mail address: gunther@math.washington.edu

Department of Mathematics, Taida Institute for Mathematical Sci- ences, and NCTS (Taipei), National Taiwan University, Taipei 106, Tai- wan.

E-mail address: jnwang@math.ntu.edu.tw