Carlos Quesada

G Xu

David González

Iciar Alfaro

Adrien Leygue

Michel Visonneau

Elías Cueto

Francisco Chinesta

Un método de descomposición propia generalizada para operadores diferenciales de alto orden

Keywords:

Introducción

Los ingenieros han acudido tradicionalmente a la reducción de los modelos como única metodología disponible para solucionar problemas complejos en épocas en las que la capacidad de cálculo disponible no era mucha. Así, por ejemplo, el modelo de viga de Euler, Bernoulli y Navier desarrolla unas ecuaciones monodimensionales para la resolución práctica del campo de desplazamientos, deformaciones y tensiones en piezas que son realmente tridimensionales. De la misma manera, el método de los modos normales para resolver problemas elastodinámicos lineales supone que la solución del problema puede escribirse como una combinación lineal de unos pocos modos que son realmente las frecuencias propias de vibración de la estructura a analizar.

Más recientemente (ya en el siglo XX, de hecho), esta metodología de reducción de la complejidad de los modelos ha sido sistematizada y aplicada con éxito en multitud de campos de las ciencias aplicadas y la ingeniería, tanto para problemas lineales [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF][START_REF] Loève | Probability theory[END_REF][START_REF] Ryckelynck | On the a priori Model Reduction: Overview and recent developments[END_REF] como no lineales [START_REF] Tenenbaum | A global framework for nonlinear dimensionality reduction[END_REF][START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF]. En general, estas técnicas, que reciben nombres distintos en distintos campos de la ciencia, y que han sido redescubiertas en multitud de ocasiones, como Principal Component Analysis (PCA), Proper Orthogonal Decomposition (POD), empirical eigenvectors, o muchos otros nombres similares, buscan funciones que permitan una buena aproximación de la solución empleando un número mínimo de ellas, es decir, de grados de libertad. Para ello, generalmente se han de resolver ciertos problemas similares al dado, de los cuales se extraen, mediante un proceso de minimización estadística, las funciones que contienen más información (energía) de la solución, en la esperanza de que constituyan, de la misma manera, una buena elección para el problema dado, que solo difiere ligeramente de aquellos de los cuales se conoce la solución completa.

En problemas definidos en coordenadas espacio-temporales, se trataría, en suma, de encontrar aproximaciones del tipo:

u(x, t) ≈ N i=1 i (x) • i (t),
donde i (t) representa una colección de coeficientes temporales y i (x) una serie de funciones dependientes únicamente de variables espaciales y que son las que mejor representan, desde un punto de vista estadístico, la solución del problema. Si consideramos los resultados de una simulación de elementos finitos, de la que dispondríamos de los valores nodales para cada incremento de tiempo, y los agrupáramos todos en una matriz A cuyas columnas fueran los valores nodales u(x j , t k ), j = 1, . . ., M, k = 1, . . ., T, entonces los N vectores i (x) se corresponderían con los autovectores asociados a los N mayores valores singulares de la matriz A.

En general, este tipo de técnicas ha gozado de un amplio predicamento en muchas ramas de la ciencia, pero posee severas limitaciones. La primera, el hecho de que sean necesarias las soluciones de ciertos problemas completos similares para determinar la base i con la que resolver el problema dado, en la esperanza de que ésta sea también óptima para él. La segunda, que su complejidad aumenta exponencialmente si el problema a considerar es paramétrico, es decir, del tipo: u = u(x, t, p 1 , p 2 , . . ., p P ), donde los distintos p m , m = 1, . . ., P, representan parámetros del problema. Intentar obtener soluciones completas mallando el hiperdominio que resultaría es simplemente imposible, por el crecimiento exponencial del número de grados de libertad de la malla resultante con el número de dimensiones del espacio de fase. Es lo que se ha dado en llamar la maldición de la dimensionalidad [START_REF] Laughlin | The theory of everything[END_REF].

Para intentar solucionar la primera de las dificultades antes mencionadas, P. Ladeveze desarrolló un método a priori, es decir, que no requería la solución de problemas similares a posteriori, obteniendo así una aproximación espacio-tiempo que denominó radial loading approximation dentro de su método de LArge Time INcrements (LATIN) [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF][START_REF] Ladeveze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF].

Más recientemente, F. Chinesta generalizó de manera independiente esta aproximación, orientándola hacia el problema multidimensional, dando así lugar a la PGD, por oposición a la descomposición propia ortogonal (POD, por sus siglas en inglés) [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. part II: transient simulation using space-time separated representations[END_REF][START_REF] Gonzalez | Recent advances on the use of separated representations[END_REF][START_REF] Chinesta | Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF]. Desde entonces la PGD ha tenido un tremendo desarrollo en multitud de aplicaciones, como lo prueba la larga lista de referencias posibles que, sin ánimo de ser exhaustivos, se incluyen: [START_REF] Ladeveze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Ammar | Reduction of the chemical master equation for gene regulatory networks using proper generalized decompositions[END_REF][START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF][START_REF] Ghnatios | Proper generalized decomposition based dynamic data-driven control of thermal processes[END_REF][START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF][START_REF] Ladeveze | The latin multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF][START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the Proper Generalized Decomposition[END_REF][START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF]. Las «reduced bases» [START_REF] Maday | Global a priori convergence theory for reduced-basis approximations of single-parameter symmetric coercive elliptic partial differential equations[END_REF] han tratado de obtener la base reducida haciendo un muestreo óptimo del espacio paramétrico. Sin embargo, este muestreo requiere conocer estimadores de error a priori cuya obtención resulta delicada en muchas circunstancias. Por ello, en lo que sigue el trabajo se centra en la PGD que resulta ser una técnica mas sistemática.

En la sección 2 se realiza un breve resumen de las características más importantes de la técnica PGD en su versión más estándar, aquella basada en la aproximación por polinomios a trozos lineales en un método de Galerkin.

Un repaso a la descomposición propia generalizada

La descomposición propia generalizada (PGD) es una metodología de reduccion de modelos a priori basada en el uso de representaciones separadas. Como se ha comentado anteriormente, inicialmente se desarrolló para la resolución de problemas estructurales no lineales en espacio-tiempo [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF], pero pronto evolucionó hacia su aplicación en modelos definidos en espacios con un alto número de dimensiones y, posteriormente, se extendió a modelos generales en mecánica computacional.

La característica principal de la PGD es que aproxima la solución exacta del problema imponiendo separación de variables. Así, si u(x), x ∈ R d , es la solución del problema, esta se aproxima en la PGD como:

u(x) ≈ N i=1 F i 1 (x 1 ) • F i 2 (x 2 ) • . . . • F i d (x d ), (1)
donde d es el número de dimensiones del problema. A diferencia de otras metodologías (como la POD), las funciones F s n (x n ) son funciones desconocidas a priori. Nótese que no es estrictamente necesario realizar una descomposición completa del espacio de fase, pudiéndose agrupar ciertas dimensiones y otras no. Por ejemplo, en [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3d solutions with 2d computational complexity[END_REF] se utiliza de manera especialmente eficiente la PGD para obtener soluciones tridimensionales en geometrías degeneradas como placas y láminas mediante una descomposición del tipo:

u(x, y, z) ≈ i F i 1 (x, y) • F i 2 (z),
donde la complejidad resultante del método es esencialmente la de la determinación de las funciones F i 1 (x, y) y, por tanto, aproximadamente la de un problema bidimensional, siendo la determinación de las funciones F i 2 (z) de coste despreciable, al ser problemas unidimensionales.

La estrategia de determinación de las funciones incógnita F i j sigue un algoritmo voraz (greedy), por lo que en cada determinación de cada sumando se selecciona un óptimo local con la perspectiva de encontrar el óptimo global. Pero cada sumando está compuesto por un producto de funciones desconocidas, que se aproximan en general mediante elementos finitos, cuyos valores nodales deben determinarse. Esto da lugar, como puede comprenderse fácilmente, a un problema no lineal. El método con el que se calculan los nuevos términos de la suma puede ser cualquiera de los muchos existentes (Newton, cuasi-Newton...), pero los autores prefieren el de punto fijo, siempre que se pueda, por su sencillez conceptual y de programación. En cada paso del algoritmo se obtiene una de las funciones separadas, que se actualizan sucesivamente para minimizar el residuo.

En este artículo se desarrolla la aplicación de la PGD a la resolución de problemas de alto orden, n > 2. Se presentan y comparan 2 estrategias alternativas para resolver el problema biarmónico (n = 4): la primera hace uso del método espectral de colocación mediante polinomios de Chebyshev, mientras que la segunda combina el método de Galerkin con polinomios de Hermite para las funciones de forma. Por supuesto, no son las únicas posibles. Otras, como el uso de funciones de aproximación suaves como MLS [START_REF] Belytschko | Element-Free Galerkin Methods[END_REF], máxima entropía local [START_REF] González | A higher order method based on local maximum entropy approximation[END_REF], etc., serían igualmente posibles.

PGD para problemas de alto orden

El problema biarmónico, además de ser un problema de referencia para métodos numéricos y analíticos, aparece en muchas aplicaciones prácticas. Por ejemplo, modela la flexión de placas delgadas o el flujo de un fluido viscoso en 2D cuando se utiliza la llamada función de corriente, por ejemplo. Por ello, se ha tomado como problema de referencia sobre el que estudiar el comportamiento de las soluciones desarrolladas en este trabajo: la flexión de placas delgadas y el problema del flujo en una cavidad forzado «por la tapa»(lid-driven cavity flow).

El problema biarmónico bidimensional se describe mediante la siguiente ecuación:

2 u = f (x, y) en , (2)
y se suponen las siguientes condiciones de contorno:

u = 0 y ∂u ∂n = 0 en ∂ , (3)
para el problema con condiciones de contorno homogéneas, donde

2 = ∂ 4 ∂x 4 + 2 ∂ 4 ∂x 2 ∂y 2 + ∂ 4 ∂y 4 .
En el caso de la PGD, a diferencia de las técnicas de elementos finitos tradicionales, la imposición de condiciones de contorno no homogéneas juega un papel preponderante en el desarrollo del método [START_REF] Gonzalez | Recent advances on the use of separated representations[END_REF]. En el siguiente epígrafe se desarrollan estos aspectos en detalle.

Un método PGD basado en la colocación de polinomios de Chebyshev

El método de colocación espectral mediante polinomios de Chebyshev se puede describir de la siguiente manera. Por un lado, se aproxima cada una de las funciones F i j de la ecuación (1) mediante el uso de polinomios de Chebyshev. Por otro, se fuerza la anulación del residuo de las ecuaciones que las gobiernan en los nodos de la discretización. Las condiciones de contorno, en este caso asociadas a u y ∂u ∂n , se deben satisfacer también en los puntos de colocación del contorno.

En este contexto, el objetivo del método consiste en hallar N parejas de funciones {X i (x), Y i (y)}, i = 1, . . ., N tales que la solución bidimensional para el problema dado por la ecuación (2) tendrá la forma:

u(x, y) = N i=1 X i (x) • Y i (y).
A pesar de que posteriormente se trabajará en forma fuerte, como corresponde a un método de colocación, el desarrollo de la PGD parte de la forma débil del problema dado por la ecuación (2), que se obtiene a partir de:

x y u * (x, y)( 2 u(x, y) -f (x, y)) dx • dy = 0,
ecuación que, como es sabido, debe cumplirse para todas las funciones test u * (x, y) seleccionadas en un espacio funcional adecuado.

A continuación se calculan las funciones que forman parte de la representación separada. Se supone conocido, porque se ha calculado previamente, el conjunto de parejas funcionales {X i (x), Y i (y)}, i = 1, . . ., n con 1 ≤ n < N. En la iteración actual se busca la pareja de enriquecimiento {R(x), S(y)} mediante la aplicación de un algoritmo de punto fijo en direcciones alternadas que, tras la convergencia, constituirá la siguiente pareja funcional {X n+1 , Y n+1 }. Por lo tanto, en la iteración n + 1 suponemos la representación separada como:

u(x, y) ≈ n i=1 X i (x) • Y i (y) + R(x) • S(y)
La función test u * (x, y) se supone entonces como:

u * (x, y) = R * (x) • S(y) + R(x) • S * (y) (4)
Introduciendo las funciones de prueba y ponderación en la formulación débil se obtiene:

[R * (x)S(y) + R(x)S * (y)] ∂ 4 R(x) ∂x 4 S(y) + 2 ∂ 2 R(x) ∂x 2 ∂ 2 S(y) ∂y 2 + R(x) ∂ 4 S(y) ∂y 4 dxdy = [R * (x)S(y) + R(x)S * (y)] f (x, y) - n i=1 ( ∂ 4 X i (x) ∂x 4 Y i (y) + 2 ∂ 2 X i (x) ∂x 2 ∂ 2 Y i (y) ∂y 2 + X i (x) ∂ 4 Y i (y)
∂y 4 ) dxdy.

(5)

Primero, se supone que R(x) es conocida, lo que implica que R * (x) = 0. De esta forma, la ecuación (5) resulta:

y S * (y) ˛Rx S(y) + 2ˇR x ∂ 2 S(y) ∂y 2 + Rx ∂ 4 S(y) ∂y 4 dy = y S * (y) Á Rx (y) - n i=1 (˛i Rx Y i (y) + 2ˇi Rx ∂ 2 Y i (y) ∂y 2 + i Rx ∂ 4 Y i (y) ∂y 4 ) dy, donde ˛Rx = x R(x) ∂ 4 R(x) ∂x 4 dx ˇRx = x R(x) ∂ 2 R(x) ∂x 2 dx Rx = x R 2 (x)dx ˛i Rx = x R(x) ∂X 4 i (x) ∂x 4 dx ˇi Rx = x R(x) ∂X 2 i (x) ∂x 2 dx i Rx = x R(x) X i (x)dx Á Rx (y) = x R(x) f (x, y)dx.
Puesto que la formulación débil se satisface para todo S * (y), es posible regresar a su forma fuerte asociada:

˛Rx S(y) + 2ˇR x ∂ 2 S(y) ∂y 2 + Rx ∂ 4 S(y) ∂y 4 = Á Rx (y) - n i=1 (˛i Rx Y i (y) + 2ˇi Rx ∂ 2 Y i (y) ∂y 2 + i Rx ∂ 4 Y i (y) ∂y 4 ).
Es esta ecuación diferencial de cuarto orden la que se resolverá usando el método de colocación de Chebyshev.

Una vez obtenida la función S(y), se busca R(x). En este caso, la función conocida es S(y) (y por lo tanto, S * (y) es nula), por lo que la ecuación (5) se expresará como:

x R * (x) ˛Sy ∂ 4 R(x) ∂x 4 + 2ˇS y ∂ 2 R(x) ∂x 2 + Sy R(x) dx = x R * (x)Á Sy (x)dx - x R * (x) n i=1 (˛i Sy ∂ 4 X i (x) ∂x 4 +2ˇi Sy ∂ 2 X i (x) ∂x 2 + i Sy X i (x)) dx,
∂x 4 + 2ˇS y ∂ 2 R(x) ∂x 2 + Sy R(x) = Á Sy (x) - n i=1 (˛i Sy ∂ 4 X i (x) ∂x 4 + 2ˇi Sy ∂ 2 X i (x) ∂x 2 + i Sy X i (x)),
que, de nuevo, se calculará mediante el método de colocación de Chebyshev. Estos 2 pasos se repetirán hasta alcanzar el punto fijo. Si se denota la función R(x) en las iteraciones actual y previa como R p (x) y R p-1 (x), respectivamente, y de la misma forma para la función S(y), S p (y) y S p-1 (y), se puede definir el error en la iteración actual como:

e = x × y (R p (x)S p (y) -R p-1 (x)S p-1 (y))
2 dxdy ≤ ε, donde ε es un valor suficientemente peque ño. Tras la convergencia, se puede dar por calculada la siguiente pareja funcional, haciendo la identificación

X n+1 = R y Y n+1 = S.
El procedimiento de enriquecimiento debe continuar hasta que se alcance la convergencia, que se puede evaluar a partir del error E:

E = 2 u -f (x, y) f (x, y) ≤ ˜ (6)
siendo ˜ otro parámetro lo suficientemente peque ño. Para la resolución de los problemas unidimensionales de cuarto orden que aparecen en la formulación anterior se supone en primer lugar una ecuación diferencial ordinaria (ODE) de cuarto orden de tipo general:

a d 4 u dx 4 + b d 2 u dx 2 + cu = g(x). (7)
La función incógnita u(x) se aproxima en el dominio x = (-1, 1) mediante la expresión:

u(x) = i=M i=1 ˛i • T i (x),
donde M denota el número de nodos considerados en x . Las coordenadas de estos nodos vienen determinadas por:

x i = cos (i -1) • M -1 , i = 1, . . ., M.
Los interpolantes T i (x) basados sobre los polinomios de Chebyshev verifican la propiedad de la delta de Kroenecker, es decir, T i (x k ) = ı ik . Para cada nodo k, con 3 ≤ k ≤ M -2 (los cuatro nodos restantes se utilizarán para la imposición de las condiciones de contorno), la ecuación discreta [START_REF] Laughlin | The theory of everything[END_REF] se puede reescribir como:

a • i=M i=1 ˛i • dT 4 i dx 4 | x k + b • i=M i=1 ˛i • dT 2 i dx 2 | x k + c • ˛k = f (x k ).
Las funciones R(x) y S(y) están sujetas a las condiciones de contorno de Dirichlet y Neumann. Se debe cumplir entonces que u(x 1 ) = u(x M ) = 0 y du dx | x 1 = du dx | x M = 0. La imposición de las condiciones da como resultado:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˛1 = 0 i=M i=1 ˛i • dT i (x) dx | x 1 = 0 ˛M = 0 i=M i=1 ˛i • dT i (x) dx | x M = 0.

Método de Galerkin con polinomios de Hermite

En esta sección se utilizará otra estrategia para hallar una solución aproximada al problema biarmónico. Se empleará el método de Galerkin para la discretización del problema y se usarán los polinomios de Hermite como funciones de forma (fig. 1). Los polinomios de Hermite permiten obtener de manera sencilla una solución con continuidad C 1 , al precio de tener que introducir los giros nodales (primera derivada de la variable esencial) como incógnitas del problema.

El procedimiento parte nuevamente de la forma débil del problema, ecuación (2), resultante de la aplicación de integración por partes y de aplicar las condiciones de contorno: La función de ponderación u * (x, y) tiene la misma expresión que la ecuación (4), de tal forma que, al introducirla junto a la función de aproximación en la formulación débil, se obtiene:

n i=1 x d 2 X i (x) dx 2 • d 2 R * (x) dx 2 dx • y Y i (y) • S(y)dy+ + x d 2 R(x) dx 2 • d 2 R * (x) dx 2 dx • y S 2 (y)dy+ + n i=1 x d 2 X i (x) dx 2 • d 2 R(x) dx 2 dx • y Y i (y) • S * (y)dy+ + x d 2 R(x) dx 2 2 dx • y S(y) • S * (y)dy+ + n i=1 x X i (x) • R * (x)dx • y d 2 Y i (y) dy 2 • d 2 S(y) dy 2 dy+ + x R(x) • R * (x)dx • y d 2 S(y) dy 2 2 dy+ + n i=1 x X i (x) • R(x)dx • y d 2 Y i (y) dy 2 • d 2 S * (y) dy 2 dy+ + x R 2 (x)dx • y d 2 S(y) dy 2 • d 2 S * (y) dy 2 dy+ +2 n i=1 x dX i (x) dx • dR * (x) dx dx • y dY i (y) dy • dS(y) dy dy+ +2 x dR(x) dx • dR * (x) dx dx • y dS(y) dy 2 dy+ +2 n i=1 x dX i (x) dx • dR(x) dx dx • y dY i (y) dy • dS * (y) dy dy+ +2 x dR(x) dx 2 dx • y dS(y) dy • dS * (y) dy dy+ - x R * (x) • f (x)dx • y S(y) • f (y)dy+ - x R(x) • f (x)dx • y S * (y) • f (y)dy = 0.
En la etapa de enriquecimiento de la PGD, se utilizará el método iterativo del punto fijo para despejar S(y), suponiendo que R(x) es conocida (y por lo tanto, R * (x) = 0). Una vez obtenido S(y), se seguirá el mismo procedimiento para obtener R(x). Estos 2 pasos se repetirán hasta que el error en la iteración actual sea menor que un determinado valor de tolerancia.

La etapa de enriquecimiento continuará hasta que se alcance la convergencia, evaluada a partir del error E, como en la ecuación [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF].

Desde este punto de vista, la formulación PGD del problema biarmónico empleando polinomios de Hermite parecería no introducir ninguna dificultad adicional sobre la formulación estándar con operadores de orden 2. Sin embargo, la imposición de condiciones de contorno esenciales no homogéneas, definidas ahora tanto en la variable esencial como en su primera derivada, merece un comentario especial.

Imposición de condiciones de contorno esenciales no homogéneas

Como es fácil comprobar, el método PGD proporciona siempre soluciones al problema homogéneo. Para la resolución, como es el caso, de problemas con condiciones de contorno Dirichlet no homogéneas, habitualmente se acude a una técnica de cambio de variable. Así, se supone que es posible encontrar una función , continua en , que satisface las condiciones de contorno no homogéneas. Entonces, la solución al problema se puede obtener como: u = + z, de forma que el problema se traslada a uno homogéneo en la variable z:

2 z = f -2 en , (8)
z = ∂z ∂n = 0 en u ∈ ∂ . (9)
Se puede encontrar información más detallada sobre cómo cambiar la condición de contorno inhomogénea a la homogénea en el contexto de la PGD en [START_REF] Gonzalez | Recent advances on the use of separated representations[END_REF].

En lo que sigue se proponen 2 alternativas:

1. Empleo de R-funciones La determinación eficiente de esta función se vuelve entonces un ingrediente fundamental de la técnica desarrollada. En algunos trabajos previos de los autores [START_REF] Gonzalez | Recent advances on the use of separated representations[END_REF], la función se determina mediante el uso de R-funciones [START_REF] Rvachev | Transfinite interpolation over implicitly defined sets[END_REF][START_REF] Rvachev | On completeness of rfm solution structures[END_REF]. Es importante hacer notar que, en el caso de problemas de alto orden, las condiciones de contorno están impuestas tanto en la variable del problema como en sus derivadas. Por tanto, obtener una función que las verifique todas y hacerlo de manera eficiente se convierte en un problema de primera magnitud. Sin embargo, la técnica de R-funciones permite una solución especialmente eficiente para este problema.

En efecto, en [START_REF] Rvachev | On completeness of rfm solution structures[END_REF] se propone una construcción de la función ϕ como sigue. Supóngase en primer lugar que el contorno se define como un Planar Straight Line Graph (PSLG), es decir, un conjunto de líneas rectas ∂ i que lo discretizan. Sea entonces ω una función normalizada, es decir, que se comporta como una pseudo-distancia:

ω| ∂ = 0, ∂ω ∂n ∂ = 1,
P = f * 0 + m k=1 1 k! f * k ω k .
Nótese que hay un polinomio P i para cada uno de los s segmentos del contorno sobre los que se define una condición esencial no homogénea. Se recomienda al lector interesado en la técnica de R-funciones acudir a la referencia [START_REF] Rvachev | On completeness of rfm solution structures[END_REF] para una explicación más detallada.

En la figura 2 se muestra un ejemplo de construcción de una R-función para un contorno esencial arbitrario con condiciones esenciales no homogéneas de primer orden. En la sección 4 se incluyen ejemplos de R-funciones con condiciones en la primera derivada de la variable esencial.

Formulación penalizada

Otra posibilidad, si no se quiere obtener la función verificando las 2 condiciones sobre la variable esencial y su primera derivada, es obtener una que solo verifique la condición u(x ∈ ) = f 0 e introducir la condición en la primera derivada en la formulación diferencial con un coeficiente de penalización, es decir En lo que sigue se aplicará la formulación penalizada combinada con la colocación espectral y el uso de las R-funciones junto con una discretización Hermite.

Resultados numéricos

En esta sección se comparan los resultados obtenidos con las 2 estrategias anteriormente presentadas para 2 problemas que se pueden modelar con la ecuación biarmónica: la flexión de una placa de Kirchhoff rectangular y la cavidad forzada (lid-driven cavity, en inglés).

Flexión de una placa de Kirchhoff rectangular

El problema propuesto en las ecuaciones (2) y (3) es equivalente al cálculo de la deformación en una placa de Kirchhoff a la que se le aplica una carga distribuida definida por la función f(x, y) (véase la fig. 3). En este caso, sin embargo, se considera una placa simplemente apoyada, de forma que las condiciones de contorno esenciales en la primera derivada del desplazamiento no están fijadas, pasando a estarlo las segundas derivadas (condiciones de Neumann, proporcionales al momento flector en el borde que, al considerarse ahora libre el giro, serán nulas).

Se ha tomado para este caso, f (x, y) = -8 4 [cos(2 x)sin 2 ( y)

+sin 2 ( x) cos(2 y) -cos(2 x) cos(2 y)],
de forma que la solución analítica del problema es u(x, y) = sin 2 ( x)sin 2 ( y).

En la figura 4 se representan los modos de la solución (un solo producto X 1 (x) • Y 1 (y)) en el enfoque colocación-Chebyshev. Como corresponde a la solución analítica del problema, ecuación [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF], la solución numérica debe obtenerse con un único producto funcional. Lo mismo sucede con el enfoque Galerkin-Hermite, que proporciona un único modo. La solución así construida se muestra en la figura 5.

Problema con condiciones de contorno no homogéneas: la cavidad forzada

Se considera ahora el problema clásico de la cavidad forzada o lid-driven cavity flow. En este problema, un fluido se encuentra u=0, ∂u/∂x=0 y x u=0, ∂u/∂y=1-x 2 u=0, ∂u/∂y=0 u=0, ∂u/∂x=0 Figura 6. Condiciones de contorno en el problema de la cavidad forzada. contenido en un depósito generalmente cuadrado, con condiciones de contorno de velocidad conocida en todo el contorno (velocidad nula en 3 de sus lados y forzada en el cuarto). Este último es precisamente el que fuerza el movimiento del fluido en el seno del depósito. En este caso, la solución de la ecuación biarmónica representa la función de la corriente, a partir de la cual puede calcularse la velocidad del fluido como:

⎧ ⎪ ⎨ ⎪ ⎩ v x = ∂u ∂y v y = - ∂u ∂x .
La geometría y las condiciones de contorno del problema de la cavidad forzada se muestran en la figura 6.

El problema se considera como biarmónico bidimensional, en el dominio = (-1, 1) 2 y sujeto a las condiciones de contorno:

u = f 0 = 0 (11a) ∂u ∂n = f 1 (11b)
donde f 0 y f 1 son funciones dadas. La forma de f 1 para este problema se muestra en la figura 6.

Enfoque colocación-Chebyshev

Se ha obtenido la función (x, y) que verifica en u(x ∈ ) = f 0 , a la que se le ha aplicado la técnica de descomposición en valores singulares (SVD) para obtener los modos previos para el método PGD. En este caso, esta operación es trivial, al ser f 0 = 0. A continuación se aplica el método de Chebyshev combinado con la PGD anteriormente expuesto, sobre la formulación penalizada, para resolver el problema con condiciones de contorno homogéneas.

La solución numérica que se obtiene aplicando el método de Chebyshev se muestra en la figura 7. Las líneas de corriente se muestran en la figura 8 y los modos necesarios en la figura 9. La evolución del error se muestra en la figura 10. La solución en este caso se ha obtenido con un número total de 8 modos.

Enfoque Galerkin-Hermite

En este enfoque, como se ha comentado anteriormente, la solución del problema biarmónico se redefine, mediante un simple cambio de variable para la imposición de condiciones de contorno, como u = z + , donde es una función cuya traza es igual a las condiciones de contorno esenciales (definidas para el operador bilaplaciano sobre la función y su primera derivada) no homogéneas impuestas al problema. Una vez obtenida la función mediante la aplicación de la técnica de R-funciones mencionada anteriormente, la forma débil para el nuevo problema en z pasa a ser, para la i-ésima iteración de la etapa de enriquecimiento y supuesto conocido R (R * = 0): contorno. Para obtenerla, se ha recurrido al método de R-funciones detallado anteriormente. Una vez obtenida la superficie, se le ha aplicado descomposición en valores singulares (SVD) para la obtención de las funciones Y j (y) y X j (x).

(Figura 1 .

 1 Figura 1. Polinomios de Hermite

=

  k = 2, . . ., m, con n la normal al contorno esencial. Entonces, se denomina normalizador de una función f(x) a la funciónf * (x) = f (x -ω∇ω).Entonces, por construcción, f * verifica quef * (x)| ∂ = f (x)| ∂ ,decir, f * se comporta esencialmente como una constante en la dirección normal al contorno. Supóngase entonces que cada segmento ∂ i , i = 1, . . ., s, del PSLG se define mediante una ecuación implícita ω i = 0. De esta forma, si las condiciones de contorno definidas sobre ∂ son del tipo u(x)| ∂ = f 0 (x), f k (x), k = 1, 2, . . ., m, entonces, si se han definido en el contorno condiciones de orden m i -1, se puede construir la función buscada de la forma P i se corresponde con los polinomios de Taylor generalizados

  x y u * (x, y)( 2 u(x, y)f (x, y)) dx dy+ u * ∂u ∂n f 1 dx dy = 0.

Figura 2 .

 2 Figura 2. Ejemplo de aplicación del método de las R-funciones a la imposición de condiciones de contorno esenciales no homogéneas. Arriba, PSLG que define el contorno esencial. Abajo, R-función resultante que, como se puede observar, se anula en el contorno esencial.

Figura 3 .

 3 Figura 3. Definición del problema de una placa de Kirchhoff en flexión, solicitada por una carga distribuida f(x, y).

Figura 4 .

 4 Figura 4. Modos obtenidos con el enfoque colocación-Chebyshev para la placa de Kirchhoff rectangular.

Figura 5 .

 5 Figura 5. Solución en desplazamientos (izquierda) y giros (derecha) para el problema de flexión de placas se la Sección 4.1. Las solución está compuesta por un único modo.

Figura 7 .Figura 8 .Figura 9 .

 789 Figura 7. Solución del problema de la cavidad forzada obtenida mediante colocación-Chebyshev.

Figura 10 .Figura 11 .

 1011 Figura 10. Solución del problema de la cavidad forzada obtenida mediante colocación-Chebyshev. Convergencia del error con el número de funciones empleadas en la aproximación.

Figura 12 .

 12 Figura 12. R-función para el problema de la cavidad forzada de la Sec ćion 4.2.

2 • d 2 SR

 22 ) • f (x)dx • y S * (y) • f (y)dy = 0.Por su parte, cuando la variable conocida es S (S * = 0): x) • R * (x)dx • y d 2 Y i (y) dy * (x) • f (x)dx • y S(y) • f (y)dy = 0. donde = m j=1 X j (x)• Y j (y). La función bidimensional representada por la matriz que se obtiene al multiplicar Y j (y) • X j (x) debe cumplir las condiciones de

Figura 13 .

 13 Figura 13. Solución para el problema de la cavidad forzada de la Sec ćion 4.2.

En la figura 11 se muestra la solución obtenida para la variable homogénea z. Por el contrario, en la figura 12 se muestra la función que permite la imposición de las condiciones esenciales, es decir, la R-función obtenida. La forma de la solución obtenida se muestran en la figura 13.5. ConclusionesEn este trabajo se ha presentado una discusión referente a la aplicación del método de la descomposición propia generalizada a problemas de alto orden, con un énfasis especial en problema de tipo bilaplaciano. Se han estudiado, de entre la pléyade de posibles formulaciones, 2 que se han considerado como representativas. La primera, basada en el uso de polinomios de Chebyshev y métodos de colocación. La segunda, basada en un método de Galerkin y polinomios de Hermite.Mientras que puede afirmarse con rotundidad que ambas aproximaciones al problema proporcionan resultados satisfactorios, la principal diferencia práctica que se ha encontrado entre ambas es la relativa a la imposición de condiciones de contorno esenciales. En el método de colocación se ha hecho uso de la penalización para imponer condiciones de contorno esenciales, mientras que en el caso de Galerkin se ha visto que la técnica de R-funciones, aunque formalmente compleja, proporciona una forma muy práctica de imponer tales condiciones. La penalización, sin embargo, actúa en todos los modos que se calculan, sumando los errores de cada modo en la solución final. Sin embargo, la imposición mediante R-funciones elimina esos problemas y proporciona formulaciones robustas, aunque para ello sea necesario calcular la descomposición en valores singulares de la R-función (que, por otra parte, tiene un coste computacional muy reducido).AgradecimientosEl trabajo recogido en el presente artículo ha sido financiado parcialmente por el Ministerio de Economía e Innovación de Espa ña, a través del proyecto CICYT-DPI2011-27778-C02-01.


Dominique Leguillon
email: aleguillo@lmm.jussieu.fr

Sébastien Murer
email: bmurer@lmm.jussieu.fr

A criterion for crack kinking out of an interface

Keywords: Brittle fracture, interfaces, crack kinking, matched asymptotics

. The mixed criterion for fracture which involves both an energetic and a stress condition (Leguillon 2002) allows defining a positive threshold of the T-stress below which no branching can occur (Selvarathinam and Goree 1998). This analysis enters within a more general mixed-mode analysis (I+II+T-stress). Despite the complex terms and the oscillations, results extend to interfacial cracks.

The assumption of a crack jump as a consequence of the energy balance allows getting rid of the problem brought by the oscillations due to these complex terms. This approach brings a new insight on the prediction of crack kinking out of a bimaterial interface.

Introduction

Within the plane strain elasticity framework, the displacement field in the vicinity of a crack tip at the interface between two homogeneous materials can be described by the so-called Williams' series made of power terms [START_REF] Cotterell | Slightly curved or kinked cracks[END_REF]. The leading terms involve a complex singularity exponent 1/ 2 i λ ε = + associated with a complex mode u and the conjugates. The first non singular term called "T-stress" corresponds to a tension acting in a direction parallel to the crack. The related power is 1; the associated mode is denoted t and T holds for the corresponding intensity factor, i.e. the remote tension

0 1/ 2 i 1/ 2 i 0 1/ 2 i 1/ 2 i ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) r R K r u K r u T r t … U r K r s K r s T … θ θ τ θ + - -+ - -  , = + + + +   , = + + +   (1)
where r and θ are the polar coordinates, with the crack tip O chosen as origin. The meaning of index 0 will be explained in the next section. The constant vector R is for consistency only. The dots refer to modes with higher order real parts (3/2, 2, 5/2 ...). Oscillations are due to the terms i r ε ± and (1) 2 bring into evidence the singular behaviour in the vicinity of the crack tip.

It is worth noting that many authors ([3], [START_REF] Larsson | Influence of non-singular stress terms and specimen geometry on small-scale yielding at crack tips in elastic-plastic materials[END_REF], [START_REF] Leevers | Inherent stress biaxiality in various fracture specimen geometries[END_REF]) emphasize on the difficulty to extract T from a finite element approximation of (1), even in the real homogeneous case. They propose different procedures leading to quite large inaccuracies. Herein, the path-independent integral H ([4], [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF]) is used and proved to be a simple and accurate procedure to determine either the real intensity factor T or the complex one K . Introducing dimensionless mix mode parameters

2i 1/ 2 i ( ) ; ( ) K T m r r M r r K K ε ε - - = = (2)
equation ( 1) can be rewritten as follows

0 1/ 2 i 0 1/ 2 i ( ) ( ( ) ( ) ( ) ( ) ( ) ) ( ) ( ( ) ( ) ( ) ( ) ( ) ) U r R K r u m r u M r t … r Kr s m r s M r
See [START_REF] He | Kinking of a crack out of an interface[END_REF] and [START_REF] Rice | Fracture mechanics concepts for interfacial cracks[END_REF] for the r dependency of m . The parameter M depends on r as well and is a generalisation of the so-called "stress biaxiality ratio" B used by Leevers and Radon [START_REF] Leevers | Inherent stress biaxiality in various fracture specimen geometries[END_REF]. The aim of this paper is to extend Leguillon's criterion ( [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF], [START_REF] Yosibash | A failure criterion for brittle elastic materials under mixed mode loading[END_REF]) in order to include the T-stress in the analysis. This criterion is a combination of energy and stress conditions which avoids any arbitrary choice of a critical length and has proved to work well to predict crack initiation at Vnotches under symmetric [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF] and complex loadings [START_REF] Yosibash | A failure criterion for brittle elastic materials under mixed mode loading[END_REF] in homogeneous materials. Williams' expansion ( 1) is used to derive the stress field and a matched asymptotic expansions procedure provides the energy release rate due to the onset of a short crack increment. The main difference with previous works is twofold, a complex singularity holds at the crack tip and the T-stress is taken into account in addition to singular modes. The crucial point lies in the comparison between two states: one relies on the initial structure, prior to crack initiation; the other refers to the same structure after onset of a short crack increment (length l ) making an angle α with the primary crack (figure 1). The corresponding displacement fields are denoted respectively 0 ( )

Matched asymptotics

r U θ , and ( ) r U θ , l
. The asymptotic procedure is carried out with respect to the small crack increment length and matching conditions with (1) lead to

0 0 1/ 2 i 1/ 2 i 1/ 2 i 1/ 2 i 1/ 2 i 1/ 2 i 1/ 2 i 1/ 2 i ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ( ) ( )] [ ( ) ( )] [ ( ) ( )] r R K u K u T t … U U r R K u V K u V U U T t W … ε ε ε ε
where X is a coefficient depending on the crack initiation angle α . It derives directly from the functions V and W and can be computed using the path-independent integral H ([4], [9]). Stress condition. The second condition is based on the maximum tension that a material can sustain before failure. It states that fracture can occur only if the opening stress along the expected crack path (defined above by the angle α and the length l ) exceeds the material strength c σ . Due to the oscillations, the forthcoming reasoning can be carried out using only an analogy with the approach used in the real case [START_REF] Leguillon | Crack deflection in a biaxial stress state[END_REF]. The stress condition where the index θθ stands for the hoop component of the tensors and s τ .

Mix criterion. The compatibility between ( 5) and ( 6) gives an equation for the crack initiation length c l as a function of α 2 ( ( ), ( )) ( ( ), ( ))

c c c c c c c X m M G Y m M α σ , = l l l l l (7)
Finally eq. ( 5) with c = l l gives a condition on K for crack initiation in the direction α : , Min ( ( ), ( ))

c f c c G K K K K X m M α α α α ≥ = = , l l (8)
The critical value K α depends on α and the actual kink angle f α maximizes the denominator, i.e.

minimizes K α giving f K (i.e. K at failure). ). The stiffer material looks roughly like Alumina while the other looks either like a porous Alumina (weak contrast) or a polymer (strong contrast). Results are in table 1. The presence of complex terms makes the interpretation quite entangled. A complete sensibility analysis is now in progress for a better understanding of the influence of the various parameters mentioned above. As already observed in the real case [START_REF] Leguillon | Crack deflection in a biaxial stress state[END_REF], f K defined in [START_REF] Leguillon | Crack deflection in a biaxial stress state[END_REF] decreases as T increases (not illustrated herein). The particular case 1 R = (no contrast) is given for comparison taking i

Discussion

I II K K K = +
. It is worth noting that the sign of ε is a convention since both + and -occur in eq. ( 1).

Conclusion

The T-stress together with the mix criterion defines an initiation length that allows getting rid of the oscillation terms in the prediction of the critical load and the deflection angle of a crack kinking out of an interface. For a complete analysis the competition between the deflection mechanism and the delamination growth must be analysed. It requires in addition the knowledge of the interface toughness, which determination is difficult, since it depends in particular on the modes mix.

Figure 1 .

 1 Figure 1. Crack kinking out of an interface.

Figure 2 . 3 -

 23 Figure 2. 3-point bending test, / 0.25 F L = , / 0.2 H L = (1/2 specimen).

Table 1 .

 1

			Im( ) / Re( ) K K	T	/ Re( ) K	f α (°)	f K (MPa.m 1/2 )	c l (µm)
	1	0	0.93		1.29	60	2.27	17
	2	0.03	1.32		2.79	50	1.39	16
	0.5 -0.03	0.89		1.33	60	0.88	18
	50 0.09	0.41		1.20	90	0.48	51
	0.02 -0.09	1.58		0.38	70	0.81	68

R

ε

where / r ρ = l holds for the dimensionless radial coordinate, and where the complex function V and the real one W are solutions to well-posed problems and vanish at infinity. They depend neither on the global geometry of the structure nor on the loading, and can be computed once and for all as functions of the crack initiation angleα .

The fracture criterion

Using the previously described matched asymptotic analysis (section 2), it has been shown [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF] that the initiation of a crack at a V-notch in a homogeneous material can be accurately predicted using two conditions involving both stress and energy. Furthermore, the proposed criterion coincides with Griffith's one for a pure crack. As a consequence of exponents greater than 1 2 / in the Williams' expansion in the general case, this initiation process is shown to be unstable; the crack jumps a short length. The presence of the T-stress term leads to a similar reasoning which is carried out herein. Energy condition. The first condition results from an energy balance between two states of the structure prior and after the onset of a short crack increment. It states that the incremental energy release rate


Habib Ammari
email: habib.ammari@math.ethz.ch

Faouzi Triki
email: faouzi.triki@univ-grenoble-alpes.fr

Chun-Hsiang Tsou
email: chun-hsiang.tsou@univ-grenoble-alpes.fr

NUMERICAL DETERMINATION OF ANOMALIES IN MULTIFREQUENCY ELECTRICAL IMPEDANCE TOMOGRAPHY

Keywords: Mathematics Subject Classification. Primary: 35R30. Key words and phrases. Inverse problems, multifrequency electric impedance tomography, anomalies reconstruction. [Ru]

The multifrequency electrical impedance tomography consists in retrieving the conductivity distribution of a sample by injecting a finite number of currents with multiple frequencies. In this paper we consider the case where the conductivity distribution is piecewise constant, takes a constant value outside a single smooth anomaly, and a frequency dependent function inside the anomaly itself. Using an original spectral decomposition of the solution of the forward conductivity problem in terms of Poincaré variational eigenelements, we retrieve the Cauchy data corresponding to the extreme case of a perfect conductor, and the conductivity profile. We then reconstruct the anomaly from the Cauchy data. The numerical experiments are conducted using gradient descent optimization algorithms.

The mfEIT Mathematical Model

Experimental research has found that the conductivity of many biological tissues varies strongly with respect to the frequency of the applied electric current within certain frequency ranges [GPG]. In [AGGJS], using homogenization techniques, the authors analytically exhibit the fundamental mechanisms underlying the fact that effective biological tissue electrical properties and their frequency dependence reflect the tissue composition and physiology. The multifrequency electrical impedance tomography (mfEIT) is a diffusive imaging modality that recovers the conductivity distribution of the tissue by using electrodes to measure the resulting voltage on its boundary, induced by two known injected currents and for many frequency values. The principal idea behind the (mfEIT) is that the dependance of the effective conductivity of the tissue with respect to the frequency of the electric current is extremely related to its state. In fact, its frequency dependence changes with its composition, membrane characteristics, intra-and extra-cellular fluids and other factors [AGGJS]. Therefore, the frequency dependence of the conductivity of the tissue can provide some information about the tissue microscopic structure and its physiological and pathological conditions. In other words, the frequency dependence of the conductivity of the tissue can help to determine if it is healthy or cancerous. The advantages of the (mfEIT) is canceling out errors due to boundary shape, the electrode positions, and other systematic errors that appear in -the more conventional imaging modality-electric impedance tomography (EIT) [Bor].

In the following we introduce the mathematical model of the (mfEIT). Let Ω be the open bounded smooth domain in R 2 , occupied by the sample under investigation and denote by ∂Ω its boundary. The mfEIT forward problem is to determine the potential u(•, ω) ∈ H 1 (Ω) := {v ∈ L 2 (Ω) : ∇v ∈ L 2 (Ω)}, solution to    -∇ • (σ(x, ω)∇u(x, ω)) = 0 in Ω, σ(x, ω)∂ ν Ω u(x, ω)(x) = f (x) on ∂Ω, ∂Ω u(x, ω)ds = 0,

(1)

where ω denotes the frequency, ν Ω (x) is the outward normal vector to ∂Ω, σ(x, ω) is the conductivity distribution, and f ∈ H -1 2 (∂Ω) := {g ∈ H -1 2 (∂Ω) : ∂Ω g ds = 0} is the input current.

In this work we are interested in the case where the frequency dependent conductivity distribution takes the form

σ(x, ω) = k 0 + (k(ω) -k 0 )χ D (x) (2)
with χ D (x) being the characteristic function of a C 2 domain D in Ω (D ⊂ Ω), k 0 being a fixed strictly positive constant, and k(ω) : R + → C \ R -, being a continuous complex-valued function.

Here k 0 represents the conductivity of the background medium, is known, and k(ω) is the conductivity of the biological tissue, given by the empirical model

k(ω) := κ 1 - κ 2 ω 2 + iωκ 3 , (3)
where κ p > 0, p = 1, 2, 3, are constants that only depend on the biological tissue properties (see for instance [AGGJS] ). The frequency profile k(ω) is somehow a meromorphic approximation with a single pole of the graph of experimental measurements for a given biological tissue [AGGJS]. It also appears as a homogenized model for periodically distributed biological cells in the dilute limit [AGGJS], and is similar to Drude models that describes the frequency dependence of the electric permittivity of a real metal within the visible frequency range [MFZ].

The mfEIT inverse problem is to determine the anomaly D and the characteristics κ p , p = 1, 2, 3, of the biological tissue from measurements of the boundary voltages u(x, ω) on ∂Ω, for ω ∈ (ω, ω), 0 ≤ ω < ω.

There have been several numerical approaches on multifrequency electrical impedance tomography. Most of them are dealing with small-volume anomalies [ABG, ABGW, GH] and frequencydifference imaging [JS, MSHA]. In small-volume-volume imaging, only the location and the multifrequency polarization tensor can be reconstructed from boundary measurements. In frequencydifference imaging, the main idea is to compare the images for different frequencies, and consider only the frequency dependent part. It was numerically shown that the approach can accommodate geometrical errors, including imperfectly known boundary. This approach which seems more natural since it aims to identify the anomaly by only focusing on the changes in the resulting images for different frequencies, is somehow simultaneously complementary and opposite to our analysis in this paper. Taking the difference between two images associated to different frequencies will remove the frequency independent part which is the keystone of the identification path pursued in this paper. Our strategy is based on the plasmonic spectral decomposition derived in [AT] which splits the electric potential on the boundary ∂Ω into two parts u = k -1 0 u 0 + u f , and separate between the frequency dependent and independent parts; see also [ADM, AMRZ]. In fact the part k -1 0 u 0 corresponds to the response of the same anomaly filled with a perfect conductor, that is k -1 0 u 0 is the limit of u when k tends to infinity. Precisely, in [AT], it was proven that the convergence of u to k -1 0 u 0 is linear in 1/k. We first process algebraically the data on the boundary and recover the frequency dependent part in order to acquire the Cauchy data of frequency independent part u 0 . Then, based on known results and approaches on the reconstruction of zero level set of harmonic functions from Cauchy data we determine the anomaly itself. The other approach is based on perturbation techniques, suppose that the contrast k(ω)/k 0 is close to 1, and linearize the problem around the harmonic function in the whole domain that shares the same flux on the boundary [AAJS]. Then, the sparsity issue comes into play and help to speed up the convergence of the iterative algorithm. In our approach, the notion of sparsity appears naturally in post-processing the data on the boundary to recover the frequency independent part. Precisely, it seems that only a finite number of eigenfunctions intervene in the linear inversion, and this can be completely mastered by the shape of the anomaly and its distance to the boundary where the measurement are taken. Finally, the perturbation approach is also complementary to our analysis since the higher the frequency is the better the recovery of the frequency independent part on the boundary is.

The paper is organized as follows. In section 2 we provide the spectral decomposition derived in [AT]. The linearization of the frequency independent part with respect to the shape of the anomaly which is necessary in our identification approach is studied in section 3. Section 4 is devoted to the retrieval of the frequency independent part on the boundary. Here, we will not follow the theoretical approach developed in [AT] based on the unique continuation of meromorphic functions. We solve the problem using algebraic tools under simplification assumptions inspired by the sparsity properties of the conductivity distribution and the behavior of the eigenvalues near the unique accumulation point 1/2. We also determine in the sequel the profile constants κ i , i = 1, 2, 3. In section 5, one we have the Cauchy data of the frequency independent part u 0 on the measurement boundary we use a conventional optimization technique to recover the anomaly. Several numerical examples are presented in section 6. Comments on the obtained results and future directions are given in the conclusion section 7.

Spectral decomposition of u(x, ω)

We first introduce an operator whose spectral decomposition will be later the corner stone of the identification of the anomaly D. Let H 1 (Ω) be the space of functions v in H 1 (Ω) satisfying ∂Ω vds = 0.

For u ∈ H 1 (Ω), we infer from the Riesz theorem that there exists a unique function T u ∈ H 1 (Ω) such that for all v ∈ H 1 (Ω),

Ω ∇T u • ∇vdx = D ∇u • ∇vdx. (4)
The variational Poincaré operator T : H 1 (Ω) → H 1 (Ω) is easily seen to be self-adjoint and bounded with norm T ≤ 1.

The spectral problem for T reads as:

Find (λ, w) ∈ R × H 1 (Ω), w = 0 such that ∀v ∈ H 1 (Ω), λ Ω ∇w • ∇vdx = D ∇w • ∇vdx.
Integrating by parts, one immediately obtains that any eigenfunction w is harmonic in D and in D = Ω \ D, and satisfies the transmission and boundary conditions

w| + ∂D = w| - ∂D , ∂ ν D w| + ∂D = (1 - 1 λ )∂ ν D w| - ∂D , ∂ ν Ω w = 0,
where w| ± ∂D (x) = lim t→0 w(x ± tν D (x)) for x ∈ ∂D. In other words, w is a solution to (1) for k = k 0 (1 -1 λ ) and f = 0.

Let H the space of harmonic functions in D and D , with zero mean ∂Ω uds(x) = 0, and zero normal derivative ∂ ν Ω u = 0 on ∂Ω, and with finite energy semi-norm

u H = Ω |∇u| 2 dx.
Since the functions in H are harmonic in D , the H is a closed subspace of H 1 (Ω). Later on, we will give a new characterization of the space H in terms of the single layer potential on ∂D associated with the Neumann function of Ω.

We remark that T u = 0 for all u in H 1 0 (D ), and T u = u for all u in H 1 0 (D) (the set of functions in H 1 (D) with trace zero).

We also remark that T H ⊂ H and hence the restriction of T to H defines a linear bounded operator. Since we are interested in harmonic functions in D and D = Ω \ D, we only consider the action of T on the closed space H . We further keep the notation T for the restriction of T to H . We will prove later that T has only isolated eigenvalues with an accumulation point 1/2. We denote by (λ - n ) n≥1 the eigenvalues of T repeated according to their multiplicity, and ordered as follows

0 < λ - 1 ≤ λ - 2 ≤ • • • < 1 2 ,
in (0, 1/2] and, similarly,

1 > λ + 1 ≥ λ + 2 ≥ • • • > 1 2 .
the eigenvalues in [1/2, 1). The eigenvalue λ ∞ = 1/2 is the unique accumulation point of the spectrum. To ease the notation we further denote the orthogonal spectral projector on the eigenspace associated to 1/2, by ∂Ω •w ± ∞ (z)ds(z)w ± ∞ (x). Next, we will characterize the spectrum of T via the mini-max principle.

Proposition 2.1. The variational Poincaré operator has the following decomposition

T = 1 2 I + K, ( 5
)
where K is a compact self-adjoint operator. Let w ± n , n ≥ 1 be the eigenfunctions associated to the eigenvalues

(λ ± n ) n≥0 . Then λ - 1 = min 0 =w∈H D |∇w(x)| 2 dx Ω |∇w(x)| 2 dx , λ - n = min 0 =w∈H ,w⊥w 1 ,••• ,w n-1 D |∇w(x)| 2 dx Ω |∇w(x)| 2 dx , = max Fn⊂H , dim(Fn)=n-1 min w∈Fn D |∇w(x)| 2 dx Ω |∇w(x)| 2 dx
, and similarly

λ + 1 = max 0 =w∈H D |∇w(x)| 2 dx Ω |∇w(x)| 2 dx , λ + n = max 0 =w∈H ,w⊥w 1 ,••• ,w n-1 D |∇w(x)| 2 dx Ω |∇w(x)| 2 dx , = min Fn⊂H , dim(Fn)=n-1 max w∈Fn D |∇w(x)| 2 dx Ω |∇w(x)| 2 dx .
We have the following decomposition of u(x, ω) in the basis of the eigenfunctions of the variational Poincaré operator T .

Theorem 2.1. [AT] Let u(x, ω) be the unique solution to the system (1).

Then, the following decomposition holds:

u(x, ω) = k -1 0 u 0 (x) + ∞ n=1 ∂Ω f (z)w ± n (z)ds(z) k 0 + λ ± n (k(ω) -k 0 ) w ± n (x), x ∈ Ω, ( 6
)
where u 0 (x) ∈ H 1 (Ω) depends only on f and D, and is the unique solution to

   ∆v = 0 in D , ∇v = 0 in D, ∂ ν Ω v = f on ∂Ω. (7)
Proof. In order to have a self-contained document we give the proof of the theorem.

We first observe that frequency dependent part

u f = u -k -1
0 u 0 , lies in H . Since the eigenfunctions w ± (x) form an orthonormal basis of H , the frequency part u f posses the following spectral decomposition:

u f (x) = ∞ n=1 Ω ∇u f (z)∇w ± n (z)dzw ± n (x), x ∈ Ω.
A forward computation leads to

Ω ∇u f (z)∇w ± n (z)dz = Ω ∇u(z)∇w ± n (z)dz.
On the other hand, since u ∈ H 1 (Ω), we obtain

Ω ∇u(z)∇w ± n (z)dz = λ ± n D ∇u(z)∇w ± n (z)dz = k 0 k(ω) λ ± n ∂D ∂ ν D u(z)| + w ± n (z)ds(z) = k 0 k(ω) λ ± n D ∇u(z)∇w ± n (z)dz - k 0 k(ω) λ ± n ∂Ω f (z)w ± n (z)ds(z).
Using the simple fact that

Ω ∇u(z)∇w ± n (z)dz = D ∇u(z)∇w ± n (z)dz + D ∇u(z)∇w ± n (z)dz,
we obtain the desired decomposition.

In [AT], assuming that the profile k(ω) is given, the spectral decomposition (6) of the solution of the forward conductivity problem has been used to retrieve the Cauchy data corresponding to the extreme case of perfect conductor u 0 . Based on unique continuation techniques, the uniqueness of the mfEIT problem, and rigorous stability estimates have been obtained from the knowledge of u 0 | ∂Ω , in the case where the anomaly is within a class of star shaped domains.

Assume that X 0 ∈ Ω, and let d 1 = dist(X 0 , ∂Ω) and let d 0 < d 1 . For δ > 0 small enough, and m > 0 large enough, define the set of anomalies:

D := D = X 0 + Υ(θ) cos θ sin θ , θ ∈ [0; 2π) ; Υ ∈ Π ,
where

Π := {d 0 < Υ(θ) < d 1 -δ; Υ(2π) = Υ(0); Υ C β ≤ m, β ≥ 2} .
In this paper we consider the reconstruction of the profile function k(ω) defined in (2), and anomalies D within the set D. At first glance, the numerical reconstruction of the inclusion in the mfEIT problem does not need to follow the path of the theoretical results derived in [AT], that is, to determine first u 0 | ∂Ω , and then find the high conductor anomaly D that produces the recovered Cauchy data of u 0 . In fact preliminary numerical calculations show that a blind minimization approach that searches the anomaly D and the profile k(ω) using boundary multifrequency data does not converge in most cases, and if it happens to converge the rate turns out to be very slow. These difficulties are well known in inverse conductivity problem, usually it is very hard to distinguish between the conductivity value and the size of the anomaly [AK]. On the other hand the numerical identification of a high conductor anomaly is a well known inverse problem, and many works have been done on it (see for instance [KS, LL]). We can cite for example quasi-reversibility type based methods [KS, LL]. Here, we will consider the parameterization type based methods [CK, ACLZ, Ru]. Since the problem is ill-posed we will use a cut-off regularization approach that consists on taking into account in the computation only the important Fourier modes of the parametrization function Υ(θ). Then, the identification is transformed into an optimization problem with a finite number of degree of freedom. The problem is still strongly nonlinear we propose here to solve it using the gradient adjoint method.

The algorithm we propose in this paper for identifying numerically the anomaly and the frequency conductivity profile is inspired by the theoretical approach developed in [AT], and it can be summarized as follows:

(i)To recover u 0 (x)| ∂Ω and κ p , p = 1, 2, 3, from the knowledge of u(x, ω)| ∂Ω , ω ∈ (ω, ω). Here, we will use an linear algebraic approach based on the understanding of the behavior of the spectrum of Neumann-Poincaré operator near its unique accumulation point, and the sparsity of the considered conductivity distribution.

(ii) To identify the anomaly D from the Cauchy data (u 0 (x), f (x)) on the boundary ∂Ω using a cut-off parameterization/Fourier approach. Here to further stabilize and speed up the convergence of the iterative gradient based method we will use two linearly independent boundary currents (f 1 , f 2 ).

The linearized map

We shall use gradient methods to identify the anomaly from the Cauchy data (u 0 (x), f (x)) on the boundary ∂Ω. In this section, we determine u h , the derivative of u 0 with respect to a shape perturbation in the direction h(x)ν D (x), where h(x) is a scalar function defined on ∂D. We will follow the analysis in [AKLZ] for a non-degenerate conductivity inside the anomaly, based on integral equations techniques.

Let D ∈ D be a given anomaly. We define

X(t) : [a, b] → R 2 to be a smooth clockwise parametrization of ∂D, where a, b ∈ R, a < b. We assume that X ∈ C β ([a, b]) and |X (t)| = 1 for all t ∈ [a, b]. Then (8) ∂D = {x = X(t), t ∈ [a, b]}.
Let h ∈ C 2 (∂D), and define the boundary of the perturbed anomaly D ε by ( 9)

∂D ε = { x = X(t) := X(t) + εh(X(t))ν D (X(t)), t ∈ [a, b]}.
Define u ε to be the unique solution to the system (7) associated to the perturbed anomaly D ε , that is

(10)        v = 0 in Ω \ D ε , ∇v = 0 in D ε , ∂ ν v = f on ∂Ω,
∂Ω vdσ = 0, where ν is the outward normal vector on ∂D ε .

The objective of this section is to derive a linear correction u h of u ε , such that

u ε = u 0 + εu h + O(ε 2 ), as ε → 0.
The main result of this section is the following.

Theorem 3.1. Let h ∈ C 2 (∂D) be fixed. Then, u h is the unique solution to the system

(11)            v = 0 in Ω \ D, ∇v = 0 in D, u| + -u| -= -h∂ ν D u 0 | + on ∂D, ∂ ν v = 0
on ∂Ω, ∂Ω vdσ = 0. Proof. We first derive an integral equation representation of the field.

Let G(x, z) = 1 2π log(|x -y|) be the Green function for the Laplacian in R 2 , and define the single layer potentials respectively on ∂D, and ∂Ω by

S D : H -1 2 (∂D) → H 1 2 (∂D), S D ϕ(x) = ∂D G(x, z)ϕ(z)ds(z),
and

S Ω : H -1 2 (∂Ω) → H 1 2 (∂Ω), S Ω ψ(x) = ∂Ω G(x, z)ψ(z)ds(z),
They satisfy the following jump relations through the boundary of respectively D and Ω [AK]

∂ ν D S D ϕ(x)| ± = ± 1 2 ϕ(x) + K * D ϕ(x), for x ∈ ∂D, ∂ ν Ω S Ω ϕ(x)| ± = ± 1 2 ϕ(x) + K * Ω ϕ(x), for x ∈ ∂Ω,
where

K * D : H -1 2 (∂D) → H -1 2 (∂D) K * D ϕ(x) = ∂D ∂ ν D (x) G(x, z)ϕ(z)ds(z),
and

K * Ω : H -1 2 (∂Ω) → H -1 2 (∂Ω) K * Ω ψ(x) = ∂Ω ∂ ν Ω (x) G(x, z)ψ(z)ds(z),
are compact operators.

Since G(x, z) is harmonic in R 2 \{z}, S Ω ψ has a unique harmonic extension in R 2 \∂Ω. Similarly, S D ϕ has a unique harmonic extension in R 2 \ ∂D. In addition, we have

S D ϕ| + = S D ϕ| - on ∂D, S Ω ψ| + = S Ω ψ| - on ∂Ω.
We also define the double layer potential D D

D D : H -1 2 (∂D) → H 1 2 loc (R 2 \ ∂D), D D ϕ(x) = ∂D ∂ D G(x, z)ϕ(z)ds(z).
It satisfies the following jump relations

∂ ν D D D ϕ(x)| + = ∂ ν D D D ϕ(x)| - for x ∈ ∂D, D D ϕ(x)| ± = ± 1 2 ϕ(x) + K D ϕ(x), for x ∈ ∂D,
where

K D : H -1 2 (∂D) → H -1 2 (∂D) K D ϕ(x) = ∂D ∂ ν D (z) G(x, z)ϕ(z)ds(z), is the L 2 (∂D)-adjoint of K * D .
The solution u 0 can be written as

(12) u 0 (x) = S D φ(x) + S Ω ψ(x) x ∈ Ω, where φ = ∂ ν D u 0 | + ∈ H -1 2 (∂D) and ψ ∈ H -1 2 (∂Ω).
Similarly, we have

(13) u ε (x) = S Dε φ ε (x) + S Ω ψ ε (x), x ∈ Ω,
where

φ ε = ∂ ν u ε | + ∈ H -1 2 (∂D ε ) and ψ ε ∈ H -1 2 (∂Ω).
Using the jump relations and the facts that ∂ ν D u 0 | -= 0 on ∂D, and ∂ ν Ω u 0 | -= f on ∂Ω, the densities φ and ψ satisfy the following system (-

1 2 I + K * D )φ(x) + ∂ ν D S Ω ψ(x) = 0 on ∂D, (14)
∂ ν Ω S D φ(x) + (- 1 2 I + K * Ω )ψ(x) = f on ∂Ω. ( 15
)
This system can be also represented in a matrix form

M φ ψ := -1 2 I + K * D ∂ ν D S Ω ∂ ν Ω S D -1 2 I + K * Ω φ ψ = 0 f . ( 16
)
The same analysis leads to the system

(17) M ε φ ε ψ ε := -1 2 I + K * Dε ∂ ν S Ω ∂ ν Ω S Dε -1 2 I + K * Ω φ ε ψ ε = 0 f .
From the parameterization of ∂D, we deduce that the outward unit normal vector ν D (x) is given by ν

D (X(t)) = R -π 2 T (X(t))
, where R -π 2 is the rotation with the angle -π 2 , and T (X(t)) = X (t) is the tangential normal vector. Let γ(X(t)) be the curvature, it satisfies (18) X (t) = γ(X(t))ν(X(t)).

Using the parameterization of ∂D ε , we deduce the following asymptotic expansion

(19) ν( x) = ν(x) -εh (t)T (x) + O(ε 2 ), for x ∈ ∂D ε ,
where h (t) = d dt h(X(t)) (we also use h (x) to denote this quantity). In the same way obtain the asymptotic expansion of the length element (20)

ds ε ( z) = ds(z)(1 -εγ(y)h(z) + O(ε 2 )).
Let Ψ ε be the diffeomorphism from ∂D onto ∂D ε given by Ψ ε (x) = x + εh(x)ν(x). From [AKLZ], we deduce the asymptotic expansion of

K * Dε (21) (K * Dε φ) • Ψ ε = K * D φ + εK (1) D φ + O(ε 2 ), where φ = φ • Ψ -1
ε , and the operator K

(1) D is defined by

K (1) D φ(x) = 1 2π ∂D [( 1 |x -y| 2 - 2 x -y, ν D (x) 2 |x -y| 4 )h(x) - x -y, T (x) |x -y| 2 h (x) - ν D (x), ν D (y) |x -y| 2 h(y) + 2 x -y, ν D (x) x -y, ν D (y) |x -y| 4 h(y) - x -y, ν D (x) |x -y| 2 γ(y)h(y)]φ(y)ds(y) (22)
Now, we calculate the asymptotic expansion of the operators ∂ ν S Ω on ∂D ε and ∂ ν Ω S D on ∂Ω.

Let ψ ∈ H -1 2 (∂Ω) be fixed. Using the relation x = x + εh(x)ν(x) ∈ ∂D ε for x ∈ ∂D, we obtain

∂ ν S Ω ψ( x) = 1 2π ∂Ω x -y, ν( x) | x -y| 2 ψ(y)ds(y), = 1 2π ∂Ω x + εh(x)ν D (x) -y, ν D (x) -εh (x)T (x) |x + εh(x)ν(x) -y| 2 ψ(y)ds(y) + O(ε 2 ), = ∂ ν D S Ω ψ(x) + ε(-h (x)∂ T S Ω ψ(x) + h(x)S (1) Ω ψ(x)) + O(ε 2 ),
where ∂ T denotes the tangential derivative, and S

Ω is defined by ( 23) S

(1)

Ω ψ(x) = 1 2π ∂Ω [ 1 |x -y| 2 - 2 x -y, ν(x) 2 |x -y| 4
]ψ(y)ds(y), for x ∈ ∂D.

We further determine the asymptotic expansion of ∂ ν Ω S Dε on ∂Ω. Let φ ∈ H -1 2 (∂D), and x ∈ ∂Ω, we have

∂ ν Ω S Dε φ(x) = 1 2π Dε x -y, ν Ω (x) |x -y| 2 φ( y)ds ε ( y), = 1 2π D x -y -εh(y)ν D (y), ν Ω (x) |x -y -εh(y)ν D (y)| 2 φ(y))(1 -εγ(y)h(y))ds + O(ε 2 ), = ∂ ν Ω S D φ(x) + ε{ 1 2π ∂D [- ν D (y), ν Ω (x) |x -y| 2 + 2 x -y, ν Ω (x) x -y, ν D (y) |x -y| 4 ]h(y)φ(y)ds(y) - 1 2π ∂D x -y, ν Ω (x) |x -y| 2 γ(y)h(y)φ(y)ds(y)} + O(ε 2 ), = ∂ ν Ω S D φ(x) + ε∂ ν Ω (D D (hφ) -S D (γhφ)) (x) + O(ε 2 ). ( 24
)
Consequently

M ε = M + εM h + O(ε 2 ),
where the operator

M h on H -1 2 (∂D) × H -1 2 (∂Ω) is defined by (25) M h := K (1) D -h ∂ T S Ω + hS (1) Ω ∂ ν Ω [D D (h•) -S D (γh•)] 0 .
So, the systems ( 16) and ( 17) imply ( 26)

φ ε ψ ε = φ ψ + ε φ h ψ h + O(ε 2 ), where φ h ψ h is given by (27) φ h ψ h = -M -1 M h φ ψ .
Thus, using the representation formula and following the same calculus, we determine the asymptotic expansion of the solution

u ε | ∂Ω (28) u ε (x) = u 0 (x) + ε(S D φ h (x) + S Ω ψ h (x) + D D (hφ)(x) -S D (γhφ)(x)) + O(ε 2 ).
We further denote

u h = S D φ h + S Ω ψ h , u h = u h + D D (hφ) -S D (γhφ).
We deduce from ( 27)

(29) ∂ ν D u h | -+ K (1) D φ -h ∂ T S Ω ψ + hS (1)
Ω ψ = 0, on ∂D, and ( 30)

∂ ν Ω u h | -+ ∂ ν Ω [D D (hφ) -S D (γhφ)] = 0,
on ∂Ω, and hence

(31) ∂ ν Ω u h | -= 0 on ∂Ω.
The equality ( 12) and the fact that ∂ T u 0 = 0 on ∂D, lead to

0 = ∂ T (h∂ T u 0 (x)) , = ∂ T (h∂ T S D φ(x)) + h (x)∂ T S Ω ψ(x) + h(x)∂ 2 T S Ω ψ(x), = ∂ T (h∂ T S D φ(x)) + h (x)∂ T S Ω ψ(x) + h(x) 1 2π ∂Ω [ -1 |x -y| 2 + 2 x -y, ν D (x) 2 |x -y| 4 + γ(x) x -y, ν D (x) |x -y| 2 ]ψ(y)ds(y), = ∂ T (h∂ T S D φ(x)) + h (x)∂ T S Ω ψ(x) -h(x)S (1) Ω ψ(x) + γ(x)h(x)∂ ν D S Ω ψ(x), which implies -h ∂ T S Ω ψ + hS (1) Ω ψ = ∂ T (h∂ T S D φ) + γh( 1 2 I -K * D )φ.
A similar calculus gives

∂ T (h∂ T S D φ(x)) = 1 2π ∂D [h (x) x -y, T (x) |x -y| 2 + h(x)( -1 |x -y| 2 + 2 x -y, ν D (x) 2 |x -y| 4 + 2γ(x) x -y, ν D (x) |x -y| 2 )]φ(y)ds(y), for x ∈ ∂D. Thus K (1) D φ(x) -h ∂ T S Ω ψ(x) + hS (1) Ω ψ(x) = 1 2π ∂D [- ν D (x), ν D (y) |x -y| 2 + x -y, ν D (x) x -y, ν D (y) |x -y| 4 ]h(y)φ(y)ds(y) - 1 2π ∂D x -y, ν D (x) |x -y| 2 γ(y)h(y)φ(y)ds(y) + 1 2 γ(x)h(x)φ(x). (32)
By the continuity of the normal derivative of double layer potentials and the jump relation, we have, for x ∈ ∂D,

∂ ν D [D D (hφ) -S D (γhφ)](x)| - = 1 2π ∂D [- ν D (x), ν D (y) |x -y| 2 + 2 x -y, ν D (x) x -y, ν D (y) |x -y| 4 ]h(y)φ(y)ds(y) + 1 2 γ(x)h(x)φ(x) - 1 2π ∂D x -y, ν D (x)
|x -y| 2 γ(y)h(y)φ(y)ds(y). ( 33) Using ( 29), ( 32) and ( 33), we have,

∂ ν D u h | -= ∂ ν D u h | -+ ∂ν D [D D (hφ) -S D (γhφ)]| -= 0,
which gives the desired result.

4. Reconstruction of u 0 (x)| ∂Ω and k(ω)

In this section we construct u 0 | ∂Ω from the knowledge of u(x, ω)| ∂Ω , ω ∈ (ω, ω). Here we recall the spectral decomposition ( 6), also valid on the boundary ∂Ω, which is the keystone of our approach.

u(x, ω) = k -1 0 u 0 (x) + ∞ n=1 ∂Ω f (z)w ± n (z)ds(z) k 0 + λ ± n (k(ω) -k 0 ) w ± n (x), x ∈ ∂Ω.
The first observation is that the functions w ± n (x) do not need to be orthogonal on the boundary ∂Ω. Then, varying the frequency, and so the coefficients of the expansion above do not guarantee the complete separation between the frequency and the non frequency parts. The second observation is that the simultaneous determination of the plasmonic resonances λ ± n , n ≥ 1, the frequency profile k(ω), and u 0 | ∂Ω is strongly nonlinear while if we assume that k(ω) and λ ± n , n ≥ 1 are given the problem becomes a linear one.

We further consider M ≥ 2 frequencies ω 1 , • • • , ω M in (ω, ω), and their associated solutions u(x, ω 1 ), • • • , u(x, ω M ). Since 1/2 is the unique accumulation point of the eigenvalues (λ ± n ) n≥1 , we only consider the N f ≥ 0 first eigenvalues as unknown variables, and we approximate the others eigenvalues by the limiting value 1/2. In fact it has been shown in [MS] that if D is C β with β ≥ 2 then for any α > -2β + 3, we have

|λ ± n -1/2| = o(n α ), n → +∞.
Thus the boundary regularity is essential to the decay rate of eigenvalues. Consequently if the boundary is C ∞ smooth, then the plasmonic eigenvalues will decay faster than any power order. Recently H. Kang and his collaborators have proved the exponential convergence of the eigenvalues in the case of analytic anomalies [AKM]. This theoretical work justifies the exponential decay behavior that has been observed numerically [PP], and checked for sample geometries like ellipses. If α > 0 is the modified maximal Grauert radius of ∂D, then

|λ ± n -1/2| = O(e -nα ), n → +∞.
These asymptotic properties of the spectrum of the Neumann-Poincaré operator suggest to consider only a finite number of them in the spectral decomposition (6). We further make the following

approximation for x ∈ Ω, 1 ≤ p ≤ M : u(x, ω p ) ≈ 2 k(ω p ) + k 0 v 1 (x), if N f = 0, (34) u(x, ω p ) ≈ k -1 0 u 0 (x) + N f n=1 1 k 0 + λ ± n (k(ω p ) -k 0 ) v ± n (x) + 2 k(ω p ) + k 0 v N f +1 (x), if N f ≥ 1, where v ± n (x) = ∂Ω f (z)w ± n (z)ds(z)w ± n (x), v N f +1 (x) = n>N f ∂Ω f (z)w ± n (z)ds(z)w ± n (x).
A simple integration by parts, leads for all n ≥ 1 (35

) ∂Ω f (z)w ± n (z)ds(z) = Ω ∇f(x)∇w ± n (x)dx,
where f is the unique solution in H 1 (Ω) to

(36) f = 0 in Ω, ∂ ν f = f on ∂Ω.
Consequently, the function

P 0 f = N f +1 n=1 v ± n (x),
where P 0 is the orthogonal projection onto the space H . On the other hand, u 0 satisfies

Ω ∇u 0 (x)∇w ± n (x)dx = Ω\D ∇u 0 (x)∇w ± n (x)dx = ∂Ω u 0 (x)∂ ν Ω w ± n (x)ds(x) - ∂D u 0 (x)∂ ν D w ± n (x)ds(x) = 0,
for all n ≥ 1.

Since P 0 u 0 = 0, and f -u 0 ∈ H the orthogonal projection of f onto the space

H is f -u 0 , that is f -u 0 = P 0 (f -u 0 ) = P 0 f. Therefore, the formula (34) becomes u(x, ω p ) ≈ k(ω p ) -k 0 k 0 (k(ω p ) + k 0 ) u 0 (x) + 2 k(ω p ) + k 0 f(x), if N f = 0, ≈ k(ω p ) -k 0 k 0 (k(ω p ) + k 0 ) u 0 (x) + 2 k(ω p ) + k 0 f(x) + N f n=1 ( 1 k 0 + λ ± n (k(ω p ) -k 0 ) - 2 k(ω p ) + k 0 )v ± n (x), if N f ≥ 1. (37)
Next, we reconstruct κ 1 , κ 2 , κ 3 and u 0 (x) by an optimization algorithm. In order to do so, we need an apriori estimation of the eigenvalues λ ± n ∈ [0, 1] for n = 1, • • • , N f . Since the eigenvalues λ ± n are within a relative narrow interval, preliminary calculations showed that the reconstruction of u 0 is indeed not very sensitive to the choice of those eigenvalues. In the rest of this section, we assume that we have we fix λ

± n ∈ [0, 1] for n = 1, • • • , N f .
Let (x j ) 1≤j≤N d ∈ ∂Ω a discretization of the boundary ∂Ω, and define, for n = 1, • • • , N f the scalar functionals

F j (U (j) 0 , V ±(j) 1 , • • • , V ±(j) N f , ω, κ 1 , κ 2 , κ 3 ) := k(ω, κ 1 , κ 2 , κ 3 ) -k 0 k 0 (k(ω, κ 1 , κ 2 , κ 3 ) + k 0 ) U (j) 0 + 2 k(ω, κ 1 , κ 2 , κ 3 ) + k 0 f(x j ) + N f n=1
( 1

k 0 + λ ± n (k(ω, κ 1 , κ 2 , κ 3 ) -k 0 ) - 2 k(ω, κ 1 , κ 2 , κ 3 ) + k 0 )V ±(j) n . ( 38
)
where (U

(j) 0 ) 1≤j≤N d and (V ±(j) n ) 1≤j≤N d are vectors in R N d , that approximate respectively (u 0 (x j )) 1≤j≤N d and (v n (x j )) 1≤j≤N d .
The scheme consists in minimizing the scalar functional

J m (U 0 , V ± 1 , • • • , V ± N f , κ 1 , κ 2 , κ 3 ) := 1 2 M p=1 N d j=1 |u(x j , ω p ) -F j (U (j) 0 , V ±(j) 1 , • • • , V ±(j) N f , ω p , κ 1 , κ 2 , κ 3 )| 2 . (39)
So, we can easily calculate its gradient from ( 38) and ( 3)

, for i = 1, 2, 3, 1 ≤ l ≤ N d and 1 ≤ n ≤ N f , ∂J m ∂κ i = M p=1 N d j=1 (u(x j , ω p ) -F j (•, ω p )) ∂F j ∂κ i (•, ω p ), ( 40
) ∂J m ∂U (l) 0 = M p=1 (u(x j , ω p ) -F l (•, ω p )) ∂F l ∂U (l) 0 (•, ω p ), (41) ∂J m ∂V ±(l) n = M p=1 (u(x j , ω p ) -F l (•, ω p )) ∂F l ∂V ±(l) n (•, ω p ), (42)
we denote here F j (U

(j) 0 , V ±(j) 1 , • • • , V ±(j)
N f , ω p , κ 1 , κ 2 , κ 3 ) by F j (•, ω p ) in order to simplify the notations.

Then, the algorithm follows the standard gradient method for 3 + N d (1 + 2N f ) variables. Once we have reconstructed the conductivity profile, i.e. the approximate values of κ 1 , κ 2 , κ 3 , we can use (37) again to calculate the approximate conductivity k(ω) by ( 3) and the approximate u 0 by the following matrix formula, let x ∈ ∂Ω, (6) The initial estimation of domain D is a centered disk with a radius 1 2 . ( 7) We consider the first 15 Fourier coefficients: N = 15. (8) We use P1 finite elements for the numerical resolution of the PDEs. (9) At each iteration, we remesh the domain to adapt to the new predicted position and shape of the domain. (10) The algorithms stop if J < 10 -5 or the number of iterations exceed 500. All the tests have executed 500 iterations. We present here several numerical simulations of the proposed algorithm. We first give the errors in the reconstruction method of u 0 in Table (1), and the errors in the reconstructed coefficients κ 1 , κ 2 , κ 3 in Table (2). The errors are computed using the L 2 -norm of the difference

(43)      u(x, ω 1 ) u(x, ω 2 ) . . . u(x, ω M )      = U (x,ω 1 ,...,ω M ) ≈        q 0 (ω 1 ) q( λ + 1 , ω 1 ) q( λ - 1 , ω 1 ) • • • q( λ - N f , ω 1 ) q 0 (ω 2 ) q( λ + 1 , ω 2 ) q( λ - 1 , ω 2 ) • • • q( λ - N f , ω 2 ) . . . . . . . . . . . . . . . q 0 (ω M ) q( λ + 1 , ω M ) q( λ - 1 , ω M ) • • • q( λ - N f , ω M )        =L( λ ± 1 ,..., λ ± N f ,ω 1 ,...ω M )        u 0 (x) v + 1 (x) v - 1 (x) . . . v - N f (x)        =V (x)
u 0reconstruct -u 0 = error(u 0reconstruct ) := ∂Ω |u 0reconstruct -u 0 | 2 dx.
We show in the following figures the targets and the reconstruction result. We calculate also the relative symmetric difference |D i D target |/|D target | during the iterations, and we draw the curves of the symmetric difference with respect to log(u i ). We finally give the relative symmetric difference of each shape in Table 3. Finally, we test a reconstruction in domain Ω that has shape different from an ellipse in Figure (5).

Concluding remarks

In this paper, by combining the spectral decomposition derived in [AT] and the linearization of the frequency independent part with respect to the shape of the anomaly, we have provided a new and efficient approach for reconstructing both the shape and conductivity parameter of a conductivity anomaly from multifrequency boundary voltage measurements. The approach and results of this paper can be extended in several directions: (i) to reconstruct multiple anomalies from multifrequency boundary measurements; (ii) to investigate the reconstruction of anisotropic conductivity anomalies from multifrequency boundary measurements, and (iii) to study elastography imaging of visco-elastic anomalies. These new developments will be reported in forthcoming works.

Acknowledgments

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01).

  Relation between J and symmetric differences.

Figure 1 .

 1 Figure 1. Exemple 1: ellipse

Figure 2 .

 2 Figure 2. Exemple 2: a square.

  Symmetric differences as a function of log(J).

Figure 3 .

 3 Figure 3. Exemple 3: a near boundary concave domain.

Figure 4 .

 4 Figure 4. Exemple 4: a centered small domain.

  Symmetric differences as a function of log(J).

Figure 5 .

 5 Figure 5. Exemple 5: A non-ellptical shaped anomaly.

Table 1 .

 1 Errors between u 0reconstruct and u 0 .

	,

Table 2 .

 2 Reconstructed constants in the frequency profile.

	•	ellipse square near-boundary small-central
	|D i D target |/|D target | 0.07055 0.12187	0.24299	0.19471

Table 3 .

 3 Relative symmetric difference.

where u(x, ω) = u(x, ω) -2 k(ω)+k 0 f(x), q 0 (ω) = k(ω)-k 0 k 0 ( k(ω)+k 0 )

, and q( λ, ω) =

Then, the vector V can be obtained by the formula (44)

where (L T L) † is the pseudo-inverse of the matrix L T L. The conditioning of the matrix L T L depends in fact on the distance between the sampling values ω j , j = 1 • • • , M, and the frequency profile (3). The approximate u 0 (x) is then recovered by taking the first coefficient of the vector V (x).

Finally, the algorithm to reconstruct u 0 , can be summarized in the following steps:

(1) Give an apriori estimation

(2) Choose a step length α m > 0 for the gradient descent.

(3) Initialize the vectors

(4) While |∇J m | is larger then a given threshold, we do (a) Calculate the values of the functions F j by ( 38), and ∇J m by ( 40), ( 41), ( 42). (b) Update the parameters

, and

(5) When |∇J m | is smaller then the threshold, we stop the iterations. ( 6) Use ( 44) with the approximate coefficients κ i obtained in the previous step to calculate the approximate value of u 0 (x) for every x ∈ ∂Ω.

Reconstruction of the anomaly from u 0

In this section, we propose a numerical method to identify the anomaly D from a finite number of Cauchy data of (u 0 (f i ), f i ), i = 1, • • • , P, on ∂Ω, where P ≥ 1. We further assume that the anomaly is located within an open subdomain Ω 0 ⊂ Ω with dist(∂Ω 0 , ∂Ω) ≥ δ 0 > 0 The scheme is based on the minimizing of a non convex functional

where u

meas are the measured Dirichlet data corresponding to the i-th Neumann data and where u is the solution to (1) associated to the current domain D ⊂ Ω 0 . In our numerical simulations we take P = 2 with f 1 = e 1 , ν Ω and f 2 = e 2 , ν Ω , where (e 1 , e 2 ) is the canonical base of R 2 .

We further assume that D is within the class D, that is, it is star shaped and its boundary ∂D can be described by the Fourier series:

Using (11) in Theorem 3.1, and integration by parts, we have the expressions of the shape derivative corresponding to each Fourier coefficient c n

and w is the solution of the following equation ( 47)

Formula ( 46) is also valid for the shape derivative corresponding to the displacement of X 0 , in these cases, h = e i , ν D , i = 1, 2. Those expressions are the basis of the following iterative algorithm:

(1) Choose an initial domain D 0 .

(2) For each iteration, i > 0: (a) Calculate the solution u i to (1), associated to the domain D i for which the boundary ∂D i is calculated by ( 45). (b) Calculate the shape derivatives ∂J ∂x 1 , ∂J ∂x 2 and ∂J ∂cn for all -N ≤ n ≤ N . (c) Choose a step length α > 0 for the gradient descent. (d) Update the parameters of the domain

If the updated domain is not entirely in Ω 0 or if R becomes negative, reduce the size of α.

(3) When J(X i , C i ) becomes smaller than a fixed threshold, we stop.

Numerical examples

The numerical tests follow the steps presented here. All the numerical experiments are done using FreeFem++ [FreeFem].

(1) Ω is a centered ellipse defined by the equation:

(2) We use two linearly independent Neumann data: f 1 = e 1 , ν Ω and f 2 = e 2 , ν Ω , where (e 1 , e 2 ) is the canonical base of R 2 . (3) The multifrequence conductivity follows the model (3) with κ 1 = 3, κ 2 = 2, κ 3 = 1 and ω are integers from 1 to 8. (4) Only the first two eigenvalues are taken into consideration, and they are fixed as follows λ + 1 = 3 4 , λ - 1 = 1 4 respectively in all cases. (5) In the algorithm to reconstruct u 0 and the conductivity profile, the initial guess of u 0 is the function f.