Rima Gandlin

Faculty of Mathematics and Computer Science

The Weizmann Institute of Science

POB 26, Rehovot 76100, ISRAEL

Achi Brandt

**Abstract **

A *direct* partial differential problem involves an
interior differential equation and a set of initial and/or boundary conditions
which stably determines a unique solution. An* inverse* problem is one in
which the differential equation and/or initial and/or boundary conditions are
not fully given and instead the results of a set of solution observations
(measurements) are known. They may contain errors, and even without errors the
problem is usually ill-posed: the known data may be approximated by widely
different solutions, but the ill-posedness of an inverse problem does not
necessarily imply more expensive solution process. On the contrary: once
the nature of the ill-posedness has been generally understood, to solve an
inverse problem may even be much less expensive than to solve a corresponding
direct problem. I am going to demonstrate this methodological point on an
example of inverse problem for Electrical Impedance Tomography (EIT).

An EIT
device for medical use consists of a set of electrodes attached to the chest of
a patient. A small *known* current is
passed between two driver electrodes. In each measurement a current is passed
through a different electrode pair, while the voltage drops at *all *the
electrodes are recorded. The collected data are used in order to compute the
conductivity distribution in a part of the patient's chest and then to display
it on a screen in order to detect anomalies, such as tumors. The electrical potential satisfies the equation
Ñª(sÑ*u*)=0, where
s is an electrical conductivity. The
set of measurements gives ideally (in the limit of many small electrodes and as
many measurements) the *Neumann to Dirichlet* mapping: the Dirichlet *
u*
boundary condition resulting from any Neumann
¶*u*/¶**n** condition. The inverse problem is to evaluate
s from this mapping.
The conductivity depends on the EIT data in a very weak way. Therefore the
inverse problem of EIT is ill-posed.

There exist some works on numerical methods for the
relevant problems, but their number is rather sparse and even those papers do
not consider the question of numerical efficiency, despite its importance for
applications. In this specific inverse EIT problem, employing local Fourier
decompositions, we have shown that all conductivity-function components of
wavelenght l are ill-posed at distances
*r*>>l from the
boundary. Hence there is no need to use at such distances fine solution grids,
since all we can know about the solution can be calculated with grids whose meshsizes
increase proportionality to *r*.

Two multigrid solvers for this problem are presented:

The **first multigrid solver** uses the Tikhonov
regularization technique. The inverse problem is reformulated as a variational
minimization problem. The resulting Euler equations form a PDE system (for *
u*, s and Lagrange multipliers), which makes the problem suitable
in principle for an effective numerical solution by multigrid methods. The
multigrid solver starts with large and then procceds to progressively smaller
regularization parameters. Many features of classical multigrid were adapted to
this particular problem (such as intergrid communications, boundary condition
treatment and coarse grid solution).

In the case of a large regularization parameter,
numerical results demonstrate a good convergence of the developed solver, but
the obtained solution is too smeared and doesn't approximate the real
conductivity too well. At small regularization the final approximation is much
better, especially near the boundary. However, in this case the system is no
longer elliptic and for efficiency the multigrid solver must use more
sophisticated relaxation scheme, which effectively decomposes the system into
its scalar factors. Although the multigrid cycles assymptotically slow down, the
final approximation to the *conductivity* is obtained by just one multigrid
cycle per grid refinement even for the smallest regularization for which
solution still exists.

In the **second multigrid solver
**a regularization in the
classical sense is replaced by *a careful choice of grids*.
In order to avoid introducing ill-posed
components into the current approximation, changes to the solution far from the
measurements are calculated only on coarse enough grids (with meshsizes at least
comparable to the distance from the measurements), which only define large-scale
*averages *of the solution. On the other hand, near the boundary of
measurements we introduce both high and low oscillatory components to the
solution. For this we use semi-coarsening in the direction parallel to the
boundary of measurements and full-coarsening in the perpendicular direction.

The algorithm is
based on three main principles: 1) due to the
ill-posedness, the computational work on the
grid of meshsize *h*
should only be done near the boundary, at
distances *O(h)*
from it; 2) similar changes introduced into
s
at the different *O(h)* distances from the
boundary of measurements cause *very different* changes in *
u* on
this boundary Þ
for smooth Fourier components in *y* a full coarsening (both in *x*
and *y*) does not resolve these differences, hence we must use *a semi
coarsening* in *y *direction (i.e. without coarsening the meshsize in
the *x* direction); 3) it is necessary to start always from the *
better-posed* changes Þ
in each cycle (FMG, V-cycle, W-cycle, semi-cycle) we introduce changes to the
current approximation first on the coarsest (or semi-coarsest) grid, where those
changes are well-posed, and only after that changes from the finer grids.

Numerical **results** show that the
second algorithm approximates s *
better* than the regularized method with its many
artificial parameters (the regularization coefficients, which should change over
the domain), for less work
(no Lagrange multipliers). The accuracy of both algorithm, in
particular near the boundary of measurements, is nearly the same. But the
second algorithm does not use any regularization in a classical
sense, and as a result significant space and time saving is achieved. In this
algorithm there is no need to solve a large system of discretized differential
equations and boundary conditions; also, no Lagrange-multiplier equations need
to be treated.