Alessandro José Queiroz Sarnaglia

Anselmo Valdério

Pascal Reisen

Céline Bondon

Lévy-Leduc

Valdério Anselmo Reisen
email: valderio.reisen@ufes.br

Pascal Bondon
email: pascal.bondon@l2s.centralesupelec.fr

Céline Lévy-Leduc
email: celine.levy-leduc@agroparistech.fr

M-regression spectral estimator for periodic ARMA models. An empirical investigation

Keywords: Periodic stationarity, PARMA models, Robust estimation, Outliers, Whittle estimation, CO concentration

Introduction

Stochastic processes exhibiting periodic correlation are frequently named as periodically correlated, periodically stationary (PS) or cyclostationary. [START_REF] Tiao | Hidden periodic autoregressive-moving average models in time series[END_REF] point out that periodic correlation may be neglected and misspecified as stationary seasonality when the standard time series tools are used. Since the introduction of PS processes in the literature by [START_REF] Gladyshev | Periodically correlated random sequences[END_REF], many authors have identified the periodic correlation phenomenon in time series of different areas, see e.g. [START_REF] Gardner | Characterization of cyclostationary random signal processes[END_REF] and [START_REF] Bloomfield | Periodic correlation in stratospheric ozone data[END_REF]. Recent reviews on PS processes can be found, for instance, in [START_REF] Gardner | Cyclostationarity: Half a century of research[END_REF] and [START_REF] Hurd | Periodically correlated random sequences: Spectral theory and practice[END_REF].

Cyclostationary counterparts of standard parametric stationary models, e.g., ARMA processes, are obtained by authorizing the parameters to vary periodically in time. In this context, the PARMA framework is a natural candidate for parsimoniously fitting PS time series. Estimation methods for PARMA models have been investigated in the literature. For example, [START_REF] Lund | Recursive prediction and likelihood evaluation for periodic ARMA models[END_REF] have considered the Gaussian maximum likelihood estimator (GMLE), [START_REF] Basawa | Large sample properties of parameter estimates for periodic ARMA models[END_REF] have studied the least square method, and [START_REF] Sarnaglia | Periodic ARMA models: Application to particulate matter concentrations[END_REF] have proposed a Whittle likelihood estimator (WLE). These authors have shown that GMLE and WLE display similar performances in finite sample sizes, while WLE is faster to calculate than GMLE, mainly in large dimensional parameter spaces.

The methodology proposed in this paper is not only motivated by its analytical and empirical features, but also from a real problem in the health and environmental fields. The concern about air pollution problems has increased considerably in the last 50 years. Especially in developing countries, the air quality has been degraded as a result of industrialization, population growth, high rates of urbanization, and inadequate or non-existent policies to control air pollution. The problems caused by air pollution produce local, regional and global impacts. Among different environmental pollution problems, air pollution is reported to cause the greatest damage to health and loss of quality of life from environmental causes.

Air quality control is accomplished by monitoring stations, which constitute the air quality monitoring network, where each station may monitor different pollutants according to the needs of the region where it is installed. Once data are collected from a monitoring system, they must be stored in data management systems and databases. Subsequently, the data must be retrieved and analyzed to see what they reveal about the effectiveness of regulatory standards, the accuracy of modeling, impacts on health endpoints, among others.

The records over time of particulate matter (PM), sulfur dioxide (SO 2 ) and CO, among other atmospheric pollutants, present some interesting features to investigate by time series tools, such as stationarity and nonstationarity (local trend, unit root), periodicity, asymmetry, weak and strong (long-memory) dependence. For example, He and Lu (2012) have studied the periodic variation of the concentrations of CO, carbon dioxide (CO 2 ) and PM at a typical traffic intersection in Hong Kong. In addition, the pollution series may present observations with high levels of concentrations producing sample densities with heavy tails. Besides the consequence of the atypical concentrations of the air pollutant on the environment and on the public health, these high levels of concentrations are not sequentially repetitive events, and can be identified as outliers which affect the statistical properties of sample functions, such as the standard mean and covariance. Since the estimation of time series models is connected with these sample functions, the estimated model can be strongly affected by the large peaks of the concentration. Therefore, it is important to develop statistical tools for a time series with these characteristics in order to describe the dynamics of the pollutants as well as to forecast future concentrations.

There are several types of outliers which cause different effects on the estimates. In general, the following three types are considered [START_REF] Denby | Robust estimation of the first order autoregressive parameter[END_REF]: innovation outliers, which affect all subsequent observations; AO and replacement outliers, which have no effect on subsequent observations. AO affect the parameter estimates more than innovation outliers, and they have the same effect as replacement outliers [START_REF] Ma | Highly robust estimation of the autocovariance function[END_REF]. In the case of PS processes, the effect of AO was discussed by Sarnaglia et al. (2010). These authors have proposed a robust autocovariance function for PS processes which is used in the periodic Yule-Walker equations to provide robust estimates for PAR models. [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF] has also suggested a robust estimation method for PAR models. [START_REF] Solci | Empirical study of robust estimation methods for PAR models with application to PM 10 data[END_REF] have presented a comparative simulation study considering the robust methods proposed by [START_REF] Ma | Highly robust estimation of the autocovariance function[END_REF], [START_REF] Shao | Robust estimation for periodic autoregressive time series[END_REF] and Sarnaglia et al. (2010) in the scenario of PAR processes with asymmetric errors and atypical observations. Recent studies in the environmental area related to the use of robust methods to estimate and forecast time series data are, for example, Ursu and Pereau (2016), [START_REF] Sguera | Functional outlier detection by a local depth with application to no x levels[END_REF][START_REF] Cabana | Robust regression based on shrinkage with application to living Environment Deprivation[END_REF], and references therein.

In the frequency domain, the periodogram is a powerful tool in the context of spectral analysis. The periodogram is a function of a time series sample and is strongly affected by outliers. One way to circumvent the outlier's effect is to use robust methods to compute alternative spectral estimators and these functions have been recently introduced as alternatives to the classical periodogram for stationary linear time series. It is well known that the periodogram is related to the least square estimator of the coefficients of a linear regression model with sine and cosine regressors, see e.g., Priestley (1981). Alternatively, several authors have defined M -periodogram by using the non-linear method of Mregression, see e.g. [START_REF] Li | Laplace periodogram for time series analysis[END_REF] and [START_REF] Li | A nonlinear method for robust spectral analysis[END_REF]. In [START_REF] Reisen | An Mestimator for the long-memory parameter[END_REF] and [START_REF] Fajardo | M-periodogram for the analysis of long-rangedependent time series[END_REF], the authors have studied the M -periodogram for long-memory processes based on the M -regression approach discussed by [START_REF] Koul | M-estimators in linear models with long range dependent errors[END_REF]. [START_REF] Reisen | Time Series and Cyclostationary process, Springer series Sarnaglia AJQ, Reisen VA, Lévy-Leduc C[END_REF] have discussed this spectral estimator for the short-memory case such as ARMA processes. Therefore, based on the above references, the asymptotic theory related to M -regression spectral estimation method in the standard stationary univariate linear time series context is now fully established.

Here, we extend the M -regression method to PARMA models by introducing a multivariate M -periodogram spectral estimator in the Whittle likelihood function, denoted as WLE M . The empirical performance of WLE M is evaluated through Monte Carlo simulations. The results show that WLE M behaves similarly to GMLE and WLE for uncontaminated time series. On the other hand, for contaminated time series with AO, GMLE and WLE are seriously affected, while WLE M gives close estimates to the ones in the uncontaminated data case. Although the asymptotic properties of WLE M method are not derived here due to their level of complexities, the finite sample size investigation supports its use in practical problems.

As an application, we consider CO daily maximum concentrations. PARMA models are fitted to the data with WLE and WLE M , respectively, and the models are used to forecast the data. In addition, we investigate whether possible atypical observations cause effects on the identification of the model and the estimation of its parameters. To the best of our knowledge, the Mregression estimation method for PARMA model and its application have not been discussed in the literature yet.

The rest of the paper is structured as follows: Section 2 describes the PARMA model with AO; Section3 introduces the M -Whittle estimation method; The finite sample performance of the robust estimator is investigated through a Monte Carlo study in Section 4; An application of the methodology to CO daily mean concentrations is presented in Section 5.

PARMA model with AO

Let Z be the set of integer numbers and (Z t ), t ∈ Z, be a real valued stochastic process satisfying E(Z 2 t ) < ∞ for all t ∈ Z. Let µ Z,t = E(Z t ) and γ Z,t (τ ) = Cov(Z t , Z t-τ ). We say that (Z t ) is PS with period S (PS S ) if for every (t, τ ) ∈ Z 2 , µ Z,t+S = µ Z,t and γ Z,t+S (τ ) = γ Z,t (τ ),

(1) and there are no smaller values of S > 0 for which (1) holds. This definition implies that µ Z,t and γ Z,t (τ ) are periodic functions in t and need to be known only for t = 1, . . . , S. If (Z t ) is PS 1 then it is weakly stationary in the usual sense. In the following, we assume without loss of generality that µ Z,t = 0 for all t ∈ Z, and we use the notation t = (r -1)S + ν where r ∈ Z and the season ν = 1, . . . , S.

One of the most popular PS S process is the PARMA model which generalizes the ARMA model, see e.g. Vecchia (1985). (Z t ) is said to be a PARMA model if it satisfies the difference equation

pν j=0 φ ν,j Z (r-1)S+ν-j = qν k=0 θ ν,k ε (r-1)S+ν-k , (2)
where, for each season ν, p ν and q ν are the AR and MA orders, respectively, φ ν,1 , . . . , φ ν,pν and θ ν,1 , . . . , θ ν,qν are the AR and MA coefficients, respectively, and φ ν,0 = θ ν,0 = 1. The innovations (ε t ) satisfy ε (r-1)S+ν = σ ν (r-1)S+ν , where ( t ) is a sequence of zero-mean independent and identically distributed random variables with unit variance. In the following, we set p = max ν p ν , q = max ν q ν , φ ν,j = 0 for j > p ν , θ ν,k = 0 for k > q ν , and we refer to (2) as the PARMA(p, q) S model.

Let (Z r ), r ∈ Z, be the S-variate time series defined by Z r = [Z (r-1)S+1 , . . . , Z (r-1)S+S ], where Z r denotes the transpose of Z r . It is well known that (Z t ) is PS S if and only if (Z r ) is weakly stationary. The covariance matrix function of (Z r ) is Γ Z (τ ) = Cov(Z r , Z r-τ ) and is related to γ Z,t (τ ) by [Γ Z (τ )] l,m = γ Z,l (τ S + l -m) for every l, m = 1, . . . , S. Now, ( 2) is equivalent to the vector ARMA (VARMA) difference equation

P j=0 Φ j Z r-j = Q k=0 Θ k ε r-k , (3)
where P = p/S , Q = q/S and x denotes the smallest integer greater than or equal to x. The entries of matrix

Φ j are [Φ j ] l,m = φ l,jS+l-m with the conven- tion that [Φ 0 ] l,m = 0 when l < m. The definition of Θ k is similar. The white noise vector process (ε r ) is defined by ε r = [ε (r-1)S+1 , . . . , ε (r-1)S+S ] and has the covariance matrix Σ ε = diag(σ 2 1 , . . . , σ 2 S ). For all z ∈ C, let Φ(z) = P j=0 Φ j z j and Θ(z) = Q k=0 Θ k z k ,
and assume that det Φ(z)Θ(z) = 0 for |z| ≤ 1. Therefore, (Z r ) is causal and invertible and the spectral density matrix of (Z r ) is

f Z (ω) = 1 2π Φ -1 (e -iω )Θ(e -iω )Σ ε Θ (e iω )Φ -1 (e iω ),
where ω ∈ (-π, π]. The causality and invertibility do not ensure that Σ ε , Φ(z) and Θ(z) are uniquely determined by the covariance matrix function Γ Z (τ ), or equivalently the spectral matrix f Z (ω), see e.g. (Brockwell and Davis, 2006, page 431). Additional restrictions have to be imposed in order to obtain identifiable models, see [START_REF] Dunsmuir | Vector Linear Time Series Models[END_REF] and [START_REF] Deistler | Vector linear time series models: Corrections and extensions[END_REF] for general conditions for VARMA models. Sarnaglia et al. ( 2020) have developed specific identifiability conditions for PARMA models. In the following, we assume that model (3) is identifiable.

Example 1 When p ν = q ν = 1 for all ν = 1, . . . , S, we have p = q = 1,

P = Q = 1, Φ 0 =        1 0 • • • • • • 0 φ 2,1 1 • • • • • • 0 0 φ 3,1 1 • • • 0 . . . . . . . . . . . . . . . 0 • • • • • • φ S,1 1        , Φ 1 =         0 • • • • • • • • • φ 1,1 0 • • • • • • • • • 0 . . . • • • • • • • • • . . . . . . • • • • • • • • • . . . 0 • • • • • • • • • 0         and det Φ(z) = det(Φ 0 + Φ 1 z) = 1 + (-1) S+1 S ν=1 φ ν,1 z.
Similar results can be derived for Θ 0 , Θ 1 and det Θ(z).

The causality condition det Φ(z) = 0 for |z| ≤ 1, and the invertibility condition det Θ(z) = 0 for |z| ≤ 1 are equivalent respectively to

ϑ φ = S ν=1 |φ ν,1 | < 1 and ϑ θ = S ν=1 |θ ν,1 | < 1. ( 4
)
As discussed in Section 1, we shall consider AO since they cause the more harmful effect in the inference of time series. Sarnaglia et al. (2010) have shown theoretically and empirically, that AO can induce a spurious memory loss by increasing the variance of the process. Let (Y t ) be defined by

Y t = Z t + ηW t , (5)
where η > 0 is the magnitude of the outlier and (W t ) is a sequence of independent and identically distributed random variables taking the values -1, 0, 1 with probabilities P(W t = 0) = 1 -ζ and P(W t = -1) = P(W t = 1) = ζ/2, where ζ ∈ (0, 1) is the probability of occurrence of an outlier. We assume that (W t ) and (Z t ) are independent processes. The probabilities for positive and negative outliers are the same. If ζ = 0 or η = 0, then Y t = Z t and (Y t ) is uncontaminated. It is worth to point out that in model ( 5), the position (in time) of the outliers is random, which seems to be appropriate in real applications, since the position and the occurrence of the outliers in the sample is usually unknown. There are other ways to describe atypical observations. For example, one can consider heavy-tailed distributions for the white noise process (ε t ), see e.g. Katkovnik (1998). However, in this kind of outlier generating mechanism, there is no explicit definition for the magnitude of the outliers and their investigation can not be performed directly.

The autocovariance function γ Y,t (τ ) = Cov(Y t , Y t-τ ) of the contaminated process (Y t ) in ( 5) is given by

γ Y,t (τ ) = γ Z,t (0) + η 2 ζ, τ = 0, γ Z,t (τ ), τ = 0, while E(Y t ) = µ Y,t = µ Z,t . Therefore, (Y t ) is also a PS S process but with larger variance than (Z t ). Let Y r = [Y (r-1)S+1 , . . . , Y (r-1)S+S ] and Γ Y (τ ) = Cov(Y r , Y r-τ ). Then Γ Y (τ ) = Γ Z (τ ) + D 1 {τ =0} , ( 6
)
where D is the (S × S) diagonal matrix with diagonal entries (D) l,l = η 2 ζ, l = 1, . . . , S. Therefore, the spectral density matrix of the contaminated vector process (Y r ) is given by

f Y (ω) = f Z (ω) + 1 2π D, ω ∈ (-π, π].
Note that, letting η → ∞ makes the diagonal matrix D to dominate f Y , which becomes close to the spectral density of a vector white noise process. This is the frequency domain counterpart of the memory loss property of processes with AO.

3 The Whittle M -estimator

We now introduce the estimation method proposed in this paper. Firstly, we define a robust alternative to the Fourier transform based on the non-linear M -regression approach. Next, the robust Whittle-type method is presented.

M -Fourier transform

Let (X r ), r ∈ Z, be any S-dimensional vector process and X = [X 1 , . . . , X N ] be a sample of size N observed from (X r ). The discrete Fourier transform of X, at the frequency ω ∈ (-π, π], is defined by

W X (ω) = (2πN ) -1/2 N r=1 X r e -irω . (7)
Expression ( 7) can be fastly computed at the Fourier frequencies ω j = 2πj/N , j = 1, . . . , N , where N = (N -1)/2 is the greatest integer smaller than or equal to (N -1)/2.

For any fixed frequency ω j , the S-dimensional vector W X (ω j ) can be viewed as a linear regression vector as follows. Let X r,ν be the νth component of vector X r , r = 1, . . . , N , ν = 1, . . . , S. Define the vector of covariates C r,j = [cos(rω j ), sin(rω j )] and consider the linear model

X r,ν = C r,j β ν (ω j ) + ξ r,ν ,
where ξ r,ν is a random error term and the coefficient vector β ν (ω j ) = (β ν,1 (ω j ), β ν,2 (ω j )) can be seen as describing the impact of the jth harmonic on the time series X 1,ν , . . . , X N,ν . The classical least square estimator of the vector β ν (ω j ) is given by βν (ω j ) = argmin

βν (ωj )∈R 2 N r=1 X r,ν -C r,j β ν (ω j ) 2 . (8) Now, define the vector d(ω j ) = [d 1 (ω j ), . . . , d S (ω j )] ,
where d ν (ω j ) = βν,1 (ω j ) -i βν,2 (ω j ). It can be shown that

W X (ω j ) = N/8π d(ω j ), j = 1, . . . , N .
It is well known that βν (ω j ) is not robust to atypical observations and heavy-tailed distributions. To improve robustness, one idea is to replace βν (ω j ) by a non-linear M -regression estimator in W X (ω j ). This will lead to the robust periodogram for PS processes proposed here.

The key idea is to replace the quadratic loss function in ( 8) by an alternative function ρ(•), which gives

argmin βν (ωj )∈R 2 N r=1 ρ(X r,ν -C r,j β ν (ω j )) .
Equivalently, one can define the M -estimator βν,ψ (ω j ) of β ν (ω j ) as the solution of

N r=1 C r,j ψ X r,ν -C r,j βν,ψ (ω j ) = 0,
where ψ is the first derivative of the Huber function given by

ρ(x) = x 2 /2, |x| ≤ δ, δ(|x| -δ/2), |x| > δ, (9)
see [START_REF] Huber | Robust Estimation of a Location Parameter[END_REF]. The choice of the tuning parameter δ > 0 is quite important and provides the compromise between robustness and efficiency of the M -estimator. Larger values of δ would lead to the Huber loss closer to the square loss, which would lead to closer performances between WLE and WLE M . On the other hand, small value of δ would improve robustness, but at the price of efficiency loss in the uncontaminated Gaussian processes. This conclusion is corroborated by other authors, see e.g. [START_REF] Maronna | Robust statistics: theory and methods (with R)[END_REF]. Following the usual choice given in the literature, the tuning parameter is set to 1.345, which ensures efficiency of 95% of the regression estimator in the Gaussian case. Finally, by defining d ψ (ω j ) similarly to d(ω j ) with βν (ω j ) replaced by βν,ψ (ω j ), the robust alternative to W X (ω j ) proposed here is given by W X,ψ (ω j ) = N/8π d ψ (ω j ), j = 1, . . . , N .

(10)

Whittle M -estimator of PARMA parameters

Assume that Z 1 , . . . , Z n is a sample from a PARMA process with known orders p and q. For simplicity, suppose that n = N S, such that every season ν = 1, . . . , S is observed N times. As previously, define the vector Z r = [Z (r-1)S+1 , . . . , Z (r-1)S+S ] corresponding to the rth cycle, the full sample being given by Z = [Z 1 , . . . , Z N ] .

The parameter vector of model ( 2) is ϕ = (ϕ φ , ϕ θ , ϕ σ ) where ϕ φ and ϕ θ contain all the AR and MA parameters, respectively, and ϕ σ = (σ 2 1 , . . . , σ 2 S ) . We define the parameter space P ⊂ R (p+q)S as the set of points (ϕ φ , ϕ θ ) for which model ( 2) is identifiable in the sense of [START_REF] Deistler | Vector linear time series models: Corrections and extensions[END_REF]. We denote by R >0 the set of positive real numbers. For any ϕ ∈ P × R S >0 , let Γ N (ϕ) be the N S × N S matrix with Γ Z (m -l) in the (l, m)th block of S × S elements, 1 ≤ l, m ≤ N . The Gaussian log-likelihood multiplied with -2N -1 is given by

LN (ϕ) = N -1 log det Γ N (ϕ) + N -1 Z Γ -1
N (ϕ)Z. We denote by ϕ 0 the true parameter vector ϕ from which the sample Z 1 , . . . , Z n is generated. We assume that (ϕ φ0 , ϕ θ0 ) ∈ P, and we have Γ N (ϕ 0 ) = Cov(Z, Z). The GMLE of ϕ 0 is φN = argmin ϕ∈P×R S >0 LN (ϕ). To obtain φN , an optimization algorithm is used and can demand high computational effort due to the fact that Γ N (ϕ) has to be inverted. To circumvent this difficulty, we use the multivariate version of Whittle's methodology to approximate LN (ϕ), see [START_REF] Dunsmuir | Vector Linear Time Series Models[END_REF] . For a PARMA process, it was shown by [START_REF] Sarnaglia | Periodic ARMA models: Application to particulate matter concentrations[END_REF] that the corresponding Whittle likelihood estimator (WLE) of ϕ 0 is φN = ( φ φ N , φ θ N , φ σ N ) where

( φ φ N , φ θ N ) = argmin (ϕ φ ,ϕ θ ) ∈P S ν=1 log σ2 N,ν (ϕ φ , ϕ θ ), σ2 N,ν (ϕ φ , ϕ θ ) = 1 2πN N j=1 Θ -1 (e -iωj )Φ(e -iωj )W Z (ω j ) ν 2 , ν = 1, . . . , S (11)
and W Z (ω j ) is given by (7) in which X r is replaced by Z r , and the νth component of φσ N is σ2 N,ν ( φφ N , φθ N ) for ν = 1, . . . , S.

The Whittle M -estimator (WLE M ) of ϕ 0 , φN,ψ , is obtained by replacing W Z (ω j ) in (11) by W Z,ψ (ω j ) defined in (10) in which X r is replaced by Z r .

It was pointed out by [START_REF] Sarnaglia | Periodic ARMA models: Application to particulate matter concentrations[END_REF] that (ϕ φ , ϕ θ ) involves (p + q)S parameters whereas the dimension of ϕ is (p + q + 1)S. Then, φN is easier to calculate and is obtained faster than φN . The same remark applies to the calculation of φN,ψ . However, the computation time of φN,ψ may be larger than φN because a numerical optimization method is needed to obtain W Z,ψ (ω j ) since this function does not have a closed form expression.

Numerical experiments

In this section, we compare the finite sample behaviors of WLE M , GMLE and WLE. Samples of the processes (Z t ) and (Y t ) were generated according to (2) and ( 5), respectively, with p = 1, q ∈ {0, 1}, S = 2, ζ = 0.01 and η = 10. The parameter values are given in Table 1. The sample sizes are n = N S = 300, 800, and the Huber function ( 9) is used with δ = 1.345. For each model, 200 sample paths were generated. Other PARMA models were also considered and the results lead to similar conclusions.

We evaluate the finite sample performances of the estimators by computing their sample root mean square errors (RMSE), when the random variables ( t ) are Gaussian and non-Gaussian in (2). The results are displayed in Tables 2, 3, 4 and 5 when ( t ) are Gaussian; in Tables 6, 7, 8 and 9 when ( t ) follow a standardized Student's t distribution with 5 degrees of freedom; in Tables 10, 11, 12 and 13 when ( t ) follow a Laplace distribution with location 0 and scale √ 2. The values with " * " refer to the sample RMSE for the contaminated series.

In the case of uncontaminated data, the three estimators present very similar sample RMSE in all cases. This quantity decreases as the sample size increases, and increases with the orders of the models, which are expected results. Relating to the estimation of the variance of the innovations, GMLE and WLE are more precise than WLE M when the data are Gaussian (Tables 2 to 5), and this is also an expected result. But, when the data are not Gaussian (Tables 6 to 13), the RMSE of the estimation of the variance of the innovations are similar for the three estimators.

In the case of contaminated data with AO, the sample RMSE of GMLE and WLE increase substantially, which means that these estimators are totally corrupted by the atypical observations, while WLE M displays similar sample RMSE values as in non-contaminated scenario. WLE M presents generally accurate estimates, even for a large number of outliers. Its superiority over GMLE and WLE appears clearly in the tables for Gaussian and non-Gaussian innovations. This empirical study indicates that GMLE and WLE should be avoided when the series contains AO, and WLE M is an alternative estimator to deal with a periodic time series with and without possible AO. These results are similar to those obtained with the M -spectral estimator in the context of standard stationary time series with short and long-range dependence discussed by [START_REF] Reisen | Time Series and Cyclostationary process, Springer series Sarnaglia AJQ, Reisen VA, Lévy-Leduc C[END_REF], [START_REF] Reisen | An Mestimator for the long-memory parameter[END_REF] and [START_REF] Fajardo | M-periodogram for the analysis of long-rangedependent time series[END_REF].

Analysis of air pollution data

We analyze the daily maximum concentrations of CO collected at the station located in Enseada do Suá area which belongs to the automatic air quality monitoring network of the Greater Vitória Region, ES, Brazil. This monitoring network has eight monitoring stations which are located at urban sites of four cities and monitors PM 10 concentrations, total suspended particles (TSP), ozone (O 3 ), nitrogen oxides (NO x ), CO, hydrocarbons (HC) and meteorological variables.

The data consist of n = N S = 728 observations measured from January 3, 2010 to December 31, 2011, and 31 observations from January 1, 2012 to January 31, 2012. The first 728 data were used to fit the model, and the last 31 records were used to perform the forecast study. The sample mean and standard deviation of the first 728 data are 2501.35 µg/m 3 and 869.41 µg/m 3 , respectively. The plot of the data over time and the boxplots for each day of the week are displayed in Figures 1 and2, respectively. These figures show that the series displays a weekly pattern (seasonality) without shifts or "abrupt" changes in the mean and long-term trends.

Figure 2 shows that the data are generally symmetric for each day of the week. We observe a high level of CO on Tuesday, September 13, 2011 (619th observation), which is approximately 5 standard deviations from the sample mean. Although, high levels of pollutants are an important information that should be considered in the context of air pollution and health impact, these observations may be identi- To verify if the 619th observation can be identified as an AO or not, a new series was generated from the original data by replacing this observation by the Hu- (ACFs) of the original and modified series were computed for each season (day of the week), and are displayed in Table 14. These values reveal that the 619th observation has some impact on the ACFs, especially at lag 1 for the seasonal periods ν = 3, 4. In these cases, the original series has smaller ACF values than the modified series. This means that an AO generates a memory-loss, see [START_REF] Fajardo | Robust estimation in long-memory processes under additive outliers[END_REF] and Sarnaglia et al. (2010). The ACF at the other periods are the same. This simple empirical investigation supports the idea that the largest observation may be considered as an AO. Now we investigate whether this observation has an impact on the identification and estimation of the model. The WLE was calculated from the centered original data obtained by subtracting the sample mean of each period. The WLE M was calculated from the centered original data obtained by subtracting the Huber location M -estimator of each period, and taking δ = 1.345 in (9). The modified WLE (MWLE) was obtained by calculating the WLE with the centered modified data got by subtracting the modified sample mean of each period. As part of the modeling strategy, the values of the Bayesian information criterion (BIC) were computed using each estimator. The GMLE was not considered in this study since GMLE and WLE have similar behaviors, see Section 4.

ν = 1 ν = 2 Eq. (4) Model φ 1,1 0 θ 1,1 0 σ 2 1 0 φ 2,1 0 θ 2,1 0 σ 2 2 0 ϑ φ ϑ θ 1 -0.2 0.0 1.0 -0.5 0.0 1.0 0.1 0.0 2 -0.2 -0.5 1.0 -0.5 -0.2 1.0 0.1 0.1 3 -1.0 0.0 1.0 -0.5 0.0 1.0 0.5 0.0 4 -1.0 -0.5 1.0 -0.5 -0.2 1.0 0.5 0.1
For WLE, BIC was computed by

BIC = N S ν=1 log(σ 2 N,ν ) + log(N ) S ν=1 (p ν + q ν ), (12)
where, for each ν = 1, ..., S, σ2 N,ν is the residual variance. For WLE M , BIC was obtained by replacing σ2 N,ν by σ2

N,ν,ψ in ( 12) and is denoted by BIC M . In the case of MWLE, BIC was computed from ( 12), where σ2 N,ν was replaced by the residual variance obtained with the modified series, and is denoted by MBIC. For the three estimators, the initial values for AR and MA estimates were set as the AR and MA coefficients estimated from a stationary ARMA(p, q) model with p = max ν p ν and q = max ν q ν , respectively, and for white noise variances set as S 2 X,ν , which is the data sample variance of the period ν. Identifiable models were fitted using the conditions proposed by [START_REF] Sarnaglia | Identifiability and Whittle Estimator of periodic ARMA models[END_REF]. BIC selected a PAR model with orders (p 1 , . . . , p S ) = (0, 0, 0, 0, 1, 1, 1), while BIC M and MBIC selected a PAR model with orders (p 1 , . . . , p S ) = (0, 0, 1, 1, 1, 1, 1). The parameters of these models and the (in sample) RMSE and median absolute deviation (MAD) between the predicted values calculated from the model and the actual values, are displayed in Table 15.

First, we observe that BIC suggested a smaller model order than MBIC. This is in agreement with the results in Table 14, and was also noticed by Sarnaglia et al. (2010) who discussed the reduction of the order of PAR models with AO when using Akaike information criterion and BIC. Since BIC has selected a different model than BIC M and MBIC, the WLE is not expected to have similar values to WLE M and MWLE, as seen in Table 15, especially for ν = 3, the season in which the atypical observation is present, and for ν = 4. Notice that even though only one observation is suspected to be atypical, it has an impact on the estimation of the parameters. Finally, WLE M and MWLE gave similar results, which corroborates the robustness of WLE M . The next step is to investigate whether the aberrant observation has impact on the forecasts. Figure 3 displays the one-step-ahead forecasts and their respective 90% forecast intervals, and Table 16 presents the RMSE and MAD for h-step-ahead forecasts for the models fitted by WLE and WLE M . The one-step-ahead 90% forecast intervals were calculated as Ẑ(r-1)S+ν-1 (1) ± 1.64σ ν , where Ẑt-1 (1) = E(Z t |Z t-1 , Z t-2 , . . .) denotes the one-step-ahead point forecast of Z t . Deviations from normality could be addressed by forecast intervals calculated with the methods addressed in [START_REF] Sarnaglia | Modeling and forecasting daily maximum hourly ozone concentrations using the regar model with skewed and heavy-tailed innovations[END_REF]. Similar conclusions related to the impact of the outlier in the confidence intervals are expected. We observe that the model fitted by WLE M provides more accurate forecasts across the h-step-ahead, h = 1, . . . , 7, since both RMSE and MAD are smaller using this robust approach. has the same finite sample performance as GMLE and WLE in the absence of outliers. Therefore, when the presence of outliers is suspected, the robust spectral estimator should be preferred. As an application to air pollution data, daily maximum CO concentrations are analyzed. The results strongly suggest that the high levels of CO concentration lead to similar features as AO do. Fitting a PARMA model to these data using the robust estimator provides also more accurate multistep predictors.

data

Fig. 1

 1 Fig. 1 Daily maximum concentrations of CO concentrations in Vitória, ES, Brazil.

Fig. 2

 2 Fig. 2 Daily maximum concentrations of CO by day of the week.

)Fig. 3

 3 Fig.3Actual observations (solid black line), one-step-ahead point forecasts (solid line), and 90% forecast intervals (dashed lines) of the CO concentrations of January 2012.

Table 1

 1 Parameters.

Table 2

 2 RMSE of the estimates for Model 1 with Gaussian innovations. RMSE of the estimates for Model 2 with Gaussian innovations. RMSE of the estimates for Model 3 with Gaussian innovations. RMSE of the estimates for Model 4 with Gaussian innovations. 300 0.272; 0.296 * 0.288; 0.311 * 0.153; 0.207 * 0.134; 0.140 * 0.148; 0.152 * 0.141; 0.183 * WLE M 800 0.150; 0.149 * 0.155; 0.154 * 0.112; 0.169 * 0.082; 0.085 * 0.094; 0.097 * 0.107; 0.145 RMSE of Model 1 with Student's t innovations.

	Method	n 300 0.067; 0.121 * 0.117; 1.366 * 0.079; 0.252 * 0.111; 1.363 * φ 1,1 σ 2 1 φ 2,1 σ 2 1
	GMLE	800 0.048; 0.101 * 0.079; 1.122 * 0.046; 0.239 * 0.074; 1.253 *
		300 0.068; 0.121 * 0.117; 1.368 * 0.079; 0.252 * 0.111; 1.364 *
	WLE	800 0.048; 0.101 * 0.079; 1.122 * 0.046; 0.239 * 0.074; 1.253 *
		300 0.067; 0.067 * 0.147; 0.179 * 0.083; 0.089 * 0.147; 0.189 *
	WLE M 800 0.051; 0.054 * 0.118; 0.149 * 0.051; 0.058 * 0.108; 0.152 Method n φ 1,1 θ 1,1 σ 2 1 φ 2,1 300 0.364; 1.393 Method n φ 1,1 σ 2 1 φ 2,1 σ 2 1 300 0.064; 0.394 * 0.107; 1.770 * 0.047; 0.159 * 0.107; 1.327 *	θ 1,1	σ 2 1
	GMLE	800 0.036; 0.373 * 0.069; 1.642 * 0.030; 0.142 * 0.070; 1.279 *
		300 0.066; 0.397 * 0.116; 1.780 * 0.047; 0.159 * 0.107; 1.328 *
	WLE	800 0.037; 0.374 * 0.073; 1.646 * 0.030; 0.142 * 0.070; 1.279 *
		300 0.080; 0.107 * 0.193; 0.334 * 0.058; 0.067 * 0.156; 0.227 *
	WLE M 800 0.053; 0.077 Method n φ 1,1 300 0.233; 1.539 Method n φ 1,1 300 0.0730; 0.1242 * 0.176; 1.296 * 0.0816; 0.2505 * 0.215; 1.493 * θ 1,1 σ 2 1 φ 2,1 σ 2 1 φ 2,1 σ 2 2	θ 1,1	σ 2 1
	GMLE	800 0.0444; 0.0960 * 0.124; 1.110 * 0.0491; 0.2480 * 0.125; 1.197 *
		300 0.0727; 0.1235 * 0.176; 1.297 * 0.0815; 0.2505 * 0.215; 1.494 *
	WLE	800 0.0445; 0.0961 * 0.124; 1.111 * 0.0495; 0.2481 * 0.125; 1.197 *
		300 0.0729; 0.0732 * 0.206; 0.185 * 0.0896; 0.0898 * 0.191; 0.183 *
	WLE M	800 0.0450; 0.0454 * 0.198; 0.171 * 0.0573; 0.0570 * 0.145; 0.115

*

Table 3 * 0.371; 1.398 * 0.120; 1.433 * 0.638; 2.219 * 0.649; 2.249 * 0.114; 1.257 * GMLE 800 0.171; 0.492 * 0.184; 0.500 * 0.065; 1.118 * 0.167; 1.007 * 0.171; 1.031 * 0.064; 1.102 * 300 0.369; 0.981 * 0.377; 0.985 * 0.119; 1.433 * 0.545; 1.309 * 0.559; 1.340 * 0.114; 1.255 * WLE 800 0.171; 0.439 * 0.184; 0.448 * 0.065; 1.119 * 0.167; 0.466 * 0.171; 0.504 * 0.064; 1.102 * 300 0.357; 0.344 * 0.371; 0.363 * 0.135; 0.164 * 0.747; 0.704 * 0.760; 0.714 * 0.150; 0.180 * WLE M 800 0.178; 0.186 * 0.193; 0.201 * 0.101; 0.132 * 0.189; 0.200 * 0.194; 0.206 * 0.101; 0.131 * Table 4 * 0.168; 0.317 * 0.038; 0.047 * 0.134; 0.215 * Table 5 * 0.252; 1.594 * 0.118; 1.307 * 0.124; 0.319 * 0.150; 0.350 * 0.107; 1.337 * GMLE 800 0.140; 0.278 * 0.146; 0.376 * 0.068; 1.230 * 0.077; 0.134 * 0.089; 0.181 * 0.071; 1.141 * 300 0.236; 0.792 * 0.255; 0.864 * 0.121; 1.313 * 0.125; 0.238 * 0.150; 0.280 * 0.107; 1.336 * WLE 800 0.141; 0.276 * 0.147; 0.372 * 0.068; 1.233 * 0.077; 0.133 * 0.089; 0.180 * 0.071; 1.141 * *

Table 6 * fied, from a statistical point of view, as being AO if they are not sequentially repetitive events.

Table 7

 7 RMSE of Model 2 with Student's t innovations. 0.305; 1.431 * 0.181; 1.256 * 0.617; 3.456 * 0.625; 3.502 * 0.204; 1.437 * GMLE 800 0.154; 0.426 * 0.164; 0.461 * 0.135; 1.202 * 0.158; 1.431 * 0.163; 1.443 * 0.139; 1.118 * RMSE of Model 3 with Student's t innovations. RMSE of Model 4 with Student's t innovations. GMLE 800 0.133; 0.263 * 0.139; 0.366 * 0.143; 1.197 * 0.073; 0.119 * 0.088; 0.158 * 0.167; 1.151 * 300 0.216; 0.216 * 0.218; 0.269 * 0.196; 1.251 * 0.125; 0.161 * 0.154; 0.215 * 0.194; 1.411 * WLE 800 0.125; 0.140 * 0.130; 0.221 * 0.143; 1.198 * 0.072; 0.093 * 0.088; 0.153 * 0.167; 1.152 * 300 0.217; 0.221 * 0.227; 0.228 * 0.181; 0.165 * 0.120; 0.127 * 0.145; 0.151 * 0.198; 0.174 * WLE M 800 0.127; 0.132 * 0.133; 0.138 * 0.129; 0.101 * 0.071; 0.073 * 0.087; 0.087 * 0.160; 0.129 * Table 10 RMSE of Model 1 with Laplace innovations.Table 11 RMSE of Model 2 with Laplace innovations.Table 12 RMSE of Model 3 with Laplace innovations.

	Method	n 300 0.291; 1.425 300 0.225; 0.256 * 0.240; 0.286 * 0.181; 1.263 * 0.262; 0.241 * 0.280; 0.294 * 0.203; 1.453 * φ 1,1 θ 1,1 σ 2 1 φ 2,1 θ 2,1 σ 2 2
	WLE	800 0.133; 0.150 * 0.141; 0.192 * 0.135; 1.205 * 0.137; 0.170 * 0.144; 0.184 * 0.139; 1.121 *
		300 0.223; 0.231 * 0.237; 0.248 * 0.204; 0.181 * 0.255; 0.256 * 0.270; 0.270 * 0.217; 0.187 *
	WLE M 800 0.145; 0.148 * 0.154; 0.155 * 0.180; 0.150 * 0.145; 0.155 * 0.149; 0.159 * 0.183; 0.154 Method n φ 1,1 σ 2 1 φ 2,1 σ 2 2 300 0.0640; 0.3882 * 0.2213; 1.8123 * 0.0502; 0.1705 * 0.194; 1.386 *
	GMLE	800 0.0395; 0.3895 * 0.1160; 1.6595 * 0.0322; 0.1398 * 0.127; 1.314 *
		300 0.0659; 0.3915 * 0.2208; 1.8239 * 0.0501; 0.1705 * 0.194; 1.382 *
	WLE	800 0.0406; 0.3911 * 0.1188; 1.6650 * 0.0321; 0.1398 * 0.127; 1.314 *
		300 0.0812; 0.1013 * 0.1717; 0.2800 * 0.0582; 0.0697 * 0.170; 0.176 *
	WLE M 800 0.0495; 0.0702 * 0.0932; 0.1882 * 0.0358; 0.0421 * 0.118; 0.110 Method n φ 1,1 θ 1,1 σ 2 1 φ 2,1 300 0.227; 1.102 Method n φ 1,1 σ 2 1 φ 2,1 σ 2 2 300 0.0725; 0.1149 * 0.200; 1.334 * 0.0746; 0.2470 * 0.182; 1.309 *	θ 2,1	σ 2 2
	GMLE	800 0.0431; 0.0940 * 0.114; 1.151 * 0.0535; 0.2493 * 0.100; 1.237 *
		300 0.0728; 0.1155 * 0.200; 1.335 * 0.0747; 0.2468 * 0.182; 1.309 *
	WLE	800 0.0432; 0.0940 * 0.115; 1.151 * 0.0535; 0.2494 * 0.100; 1.237 *
		300 0.0715; 0.0701 * 0.294; 0.266 * 0.0894; 0.0888 * 0.222; 0.199 *
	WLE M	800 0.0418; 0.0398 * 0.272; 0.241 * 0.0647; 0.0599 * 0.190; 0.155 *
	Method	n 300 0.319; 1.311 * 0.327; 1.297 * 0.190; 1.350 * 2.087; 3.162 * 2.088; 3.192 * 0.166; 1.377 * φ 1,1 θ 1,1 σ 2 1 φ 2,1 θ 2,1 σ 2 2
	GMLE	800 0.167; 0.699 * 0.177; 0.740 * 0.113; 1.173 * 0.162; 1.409 * 0.161; 1.401 * 0.110; 1.125 *
		300 0.247; 0.224 * 0.252; 0.263 * 0.190; 1.353 * 0.247; 0.260 * 0.271; 0.307 * 0.165; 1.397 *
	WLE	800 0.145; 0.153 * 0.154; 0.201 * 0.113; 1.175 * 0.145; 0.154 * 0.146; 0.178 * 0.110; 1.128 *
		300 0.246; 0.247 * 0.259; 0.260 * 0.288; 0.264 * 0.268; 0.272 * 0.283; 0.285 * 0.263; 0.233 *
	WLE M	800 0.137; 0.139 * 0.146; 0.149 * 0.244; 0.211 * 0.148; 0.155 * 0.147; 0.155 * 0.232; 0.200 *
	Method	n 300 0.0615; 0.3880 * 0.178; 1.696 * 0.0499; 0.1661 * 0.172; 1.405 * φ 1,1 σ 2 1 φ 2,1 σ 2 2
	GMLE	800 0.0392; 0.3823 * 0.104; 1.626 * 0.0304; 0.1423 * 0.102; 1.286 *
		300 0.0642; 0.3913 * 0.187; 1.715 * 0.0498; 0.1664 * 0.172; 1.403 *
	WLE	800 0.0391; 0.3838 * 0.105; 1.632 * 0.0305; 0.1423 * 0.103; 1.286 *
		300 0.0793; 0.1003 * 0.189; 0.297 * 0.0518; 0.0622 * 0.191; 0.174 *
	WLE M	800 0.0477; 0.0705 * 0.104; 0.200 * 0.0350; 0.0430 * 0.138; 0.106 *
	ber location M -estimator of the corresponding season	(ν = 3, Tuesday). The sample autocovariance functions

* *

Table 8 * Table 9 * 0.230; 1.154 * 0.196; 1.247 * 0.124; 0.208 * 0.151; 0.244 * 0.194; 1.407 *

Table 13

 13 RMSE of Model 4 with Laplace innovations.

	Method	n 300 0.203; 2.041 * 0.215; 2.080 * 0.170; 1.468 * 0.124; 0.634 * 0.147; 0.644 * 0.168; 1.273 * φ 1,1 θ 1,1 σ 2 1 φ 2,1 θ 2,1 σ 2 2
	GMLE	800 0.131; 0.286 * 0.133; 0.381 * 0.102; 1.280 * 0.072; 0.123 * 0.085; 0.176 * 0.116; 1.164 *
			300 0.192; 0.205 * 0.200; 0.260 * 0.171; 1.475 * 0.123; 0.186 * 0.147; 0.230 * 0.168; 1.281 *
		WLE	800 0.125; 0.142 * 0.126; 0.221 * 0.102; 1.283 * 0.071; 0.098 * 0.085; 0.173 * 0.116; 1.165 *
			300 0.201; 0.204 * 0.215; 0.218 * 0.196; 0.191 * 0.124; 0.130 * 0.147; 0.151 * 0.250; 0.217 *
	WLE M	800 0.127; 0.131 * 0.134; 0.132 * 0.154; 0.121 * 0.069; 0.069 * 0.082; 0.081 * 0.221; 0.187 *
	CO (µg m 3 )	1000 3000 5000 7000
			Sun	Mon	Tue	Wed	Thu	Fri	Sat
					Day of the week

Table 14

 14 ACFs of the original and modified series.

	ν	lag 1	Original lag 2	lag 3	lag 1	Modified lag 2	lag 3
	1 0.130 0.141 -0.048	0.130 0.141 -0.048
	2 0.045 0.178	0.181	0.045 0.178	0.181
	3 0.206 0.138	0.105	0.313 0.096	0.153
	4 0.205 0.110	0.052	0.219 0.110	0.052
	5 0.238 0.127	0.029	0.238 0.046	0.029
	6 0.554 0.267	0.151	0.554 0.267	0.108
	7 0.314 0.184	0.160	0.314 0.184	0.160

Table 15

 15 WLE, WLE M and MWLE for each season ν.

	ν	φ ν,1	WLE	σ ν	WLE M φ ν,1 σ ν	MWLE φ ν,1 σ ν	S X,ν
	1	-	475.43	-	503.20	-	475.43 475.43
	2	-	835.60	-	830.62	-	835.60 835.61
	3	-	901.63	-0.313 767.41	-0.293 742.58 901.63
	4	-	781.63	-0.246 770.44	-0.219 762.73 781.63
	5	-0.244 775.66	-0.296 794.00	-0.243 775.66 798.65
	6	-0.581 696.79	-0.550 769.66	-0.581 696.79 837.09
	7	-0.278 701.43	-0.274 695.97	-0.278 701.43 738.88
	RMSE	749.03			744.71		720.77	869.41
	MAD	482.31			471.95		469.32	583.98

Table 16

 16 RMSE and MAD for h-step-ahead forecasts of the CO concentrations of January 2012.This paper proposes to use the M -regression method to estimate PARMA models in the frequency domain. Contrarily to the usual GMLE and WLE, the proposed spectral estimator is robust against AO. Moreover, it

					CO (µg m 3	1000 2000 3000 4000	Actual WLE WLE M
					2012-01-01	2012-01-31
	h	RMSE WLE WLE M	MAD WLE WLE M
	1 660.29	630.70	559.26	442.78
	2 661.70	655.85	569.76	547.25
	3 662.41	652.47	559.26	539.53
	4 638.23	628.78	526.07	492.95
	5 624.54	616.44	492.88	457.57
	6 636.13	628.26	526.07	498.55
	7 631.67	625.05	492.88	457.57
	6 Conclusions

  . Biometrika 67:365-73 Ursu E, Pereau JC (2016) Application of periodic autoregressive process to the modeling of the garonne river flows. Stochastic Environmental Research and Risk Assessment 30(7):1785-1795 Vecchia AV (1985) Periodic autoregressive-moving average (PARMA) modeling with applications to water resources. Water Resour Bull 21:721-30

Acknowledgements The authors would like to thank the agencies CNPq, CAPES and FAVES, Brazil and CNRS and CentraleSupélec, France for their financial support. This research was also partially supported by DATAIA Convergence Institute as part of the "Programme d'Investissement d'Avenir, (ANR17-CONV-0003) operated by CentraleSupélec", and by the iCODE Institute, research project of the IDEX Paris-Saclay, and by the Hadamard Mathematics LabEx (LMH) through the grant number ANR-11-LABX-0056-LMH in the "Programme d'Investissement d'Avenir". The authors also thank the referees and the editor for their valuable suggestions which have contributed significantly for the improvement of the paper.


Mohamed Iguernane
email: mohamed.iquernane@iecn.u-nancy.fr

Serguei A Nazarov
email: srgnazarov@yahoo.co.uk

Jean-Rodolphe Roche
email: roche@iecn.u-nancy.fr

Jan Sokolowski
email: jan.sokolowski@iecn.u-nancy.fr

Katarzyna Szulc
email: katarzyna.szulc@iecn.u-nancy.fr

TOPOLOGICAL DERIVATIVES FOR SEMILINEAR ELLIPTIC EQUATIONS

Keywords: Shape optimization, topological derivative, levelset method, variational inequality, asymptotic analysis

The form of topological derivatives for integral shape functional is derived for a class of semilinear elliptic equations. The convergence of finite element approximation for the topological derivatives is shown and the error estimates in the L ∞ norm are obtained. Results of numerical experiments which confirm the theoretical convergence rate are presented.

Introduction 1.Topological derivatives in shape optimization

Topological derivatives are introduced for linear problems in [START_REF] Sokolowski | Modelling of topological derivatives for contact problems[END_REF], and for variational inequalities in [START_REF] Sokolowski | Introduction to shape optimization[END_REF]. The mathematical theory of asymptotic analysis is applied in [START_REF] Nazarov | Asymptotic analysis of shape functional[END_REF], [START_REF] Nazarov | Self-adjoint extensions for the Neumann Laplace and applications[END_REF] for derivation of topological derivatives in shape optimization of elliptic boundary values problems. The numerical solutions of the shape optimization problems for variational inequalities obtained by the level set method combined with the topological derivatives are presented in [START_REF] Fulmanski | A level set method in shape and topology optimization for variational inequalities[END_REF] In the paper we present the topological derivatives for semilinear elliptic boundary value problems. In the first part, the asymptotic analysis is performed for a class of the second order semilinear elliptic equations. In the second part the convergence of the finite element approximation for the topological derivatives is established, and results of numerical experiments are presented.

The topological sensitivity analysis aims to provide an asymptotic expansion of a shape functional with respect to the size of a small hole created inside the domain. For the sake of simplicity we assume that the hole is located at the origin, but the same results hold for the center of the hole being an arbitrary point in Ω. For a criterion j(Ω) = J (u Ω ; Ω)

where Ω ⊂ | N (N = 2 or 3) and u Ω is a solution of a of partial differential equation defined over Ω, this expansion can be generally written in the form:

j(Ω \ (O + ω ε )) -j(Ω) = f (ε)T Ω (O, ω) + o(f (ε)). (1)
Here ε and O denote respectively the diameter and the center of the hole, ω is a fixed domain containing the origin O and f (ε) is a positive function tending to zero with ε → 0. The coefficient T Ω is commonly called the topological derivative.

Semilinear elliptic equation

Let Ω and ω be bounded domains in | 3 with the smooth boundaries ∂Ω and ∂ω and the compact closures Ω and ω, respectively. The origin O of the coordinate system is assumed to belong to the domains Ω and ω. The following sets are introduced:

ω ε = x ∈ | 3 : ξ := ε -1 x ∈ ω , Ω(ε) := Ω \ ω ε , (2)
where x = (x 1 , x 2 , x 3 ) are Cartesian coordinates in the domain Ω and ε > 0 is a small parameter. The upper bound ε 0 > 0 is chosen in such a way that for ε ∈ (0, ε 0 ] the set ω ε belongs to the domain Ω. We can diminish the value of ε 0 > 0 in the sequel, if necessary, however, the notation for the bound ε 0 remains unchanged. The set ω ε is called hole, or opening, in the domain Ω(ε).

In this paper, we consider a nonlinear elliptic problem in the singularly perturbed domain Ω(ε) :

-∆ x u ε (x) = F (x, u ε (x)), x ∈ Ω(ε),

u ε (x) = 0, x ∈ ∂Ω(ε). (3)
Here F ∈ C 0,α (Ω × |) and f ∈ L 2 (Ω) are given functions, independent of the parameter ε. The asymptotic analysis in the linear case is well known (see monographs [START_REF] Il | in: Matching of asymptotic expansions of solutions of boundary value problems[END_REF], [START_REF] Mazja | Asymptotische Theorie elliptischer Randwertaufgaben in singulär gestörten[END_REF]), e.g. for the Dirichlet boundary value problem for the Poisson equation:

-∆ x u ε (x) = f (x), x ∈ Ω(ε), u ε (x) = 0, x ∈ ∂Ω(ε), (4)
According to the method of compound asymptotic expansions [START_REF] Mazja | Asymptotische Theorie elliptischer Randwertaufgaben in singulär gestörten[END_REF], in asymptotic analysis of (4) there appear two limit problems. The first one is obtained by formally taking ε = 0, e.g. filling the hole ω ε :

-∆ x u(x) = f (x), x ∈ Ω, u(x) = 0, x ∈ ∂Ω, (5)
and the second one is the boundary value problem which furnishes the leading boundary layers term:

-∆ ξ w(ξ) = 0, ξ ∈ | 3 \ ω, w(ξ) = -u(O), ξ ∈ ∂ω, (6)
where u(O) is the value at the origin of the solution of (5).

As in [START_REF] Mazja | On the asymptotic behavior of solutions of elliptic boundary value problems with irregular perturbations of the domain[END_REF] (see also [START_REF] Mazja | Asymptotische Theorie elliptischer Randwertaufgaben in singulär gestörten[END_REF]; §5.7), for the nonlinear problem (3) we obtain also two limit problems, the first limit problem is nonlinear

-∆ x v(x) = F (x, v(x)), x ∈ Ω, v(x) = 0, x ∈ ∂Ω, (7)
and the second limit problem is the linear exterior problem [START_REF] Demlow | Sharply localized pointwise and W -1 ∞ estimates for finite element methods for quasilinear problems[END_REF] with u(O) := v(O) given by the solution to [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF]. Our aim in this paper is the construction of asymptotic approximations for solutions to (3) in such a way that we are able to obtain an expansion of a given shape functional

J (u ε ; Ω(ε)) = Ω(ε) J(x, u ε (x))dx, (8)
of the first order with respect to ε, namely,

J (u ε ; Ω(ε)) = J (v; Ω) + εT Ω (O) + o(ε), (9)
(cf. ( 1)), where

J (v; Ω) = Ω J(x, v(x))dx, (10)
and T Ω is the topological derivative of the functional J .

Beside that we need the linearized problem [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF], which gives us the regular terms in the asymptotic approximation,

-∆ x V (x) -F ′ v (x, v(x))V (x) = F(x), x ∈ Ω, V (x) = g(x), x ∈ ∂Ω (11)
The solution V is but the so-called adjoint state. The adjoint state is introduced in order to simplify the expression for the topological derivative. Appropriate function spaces are employed to analyze the solvability of all boundary value problems introduced above.

The weighted Hölder spaces Λ l,α β (Ω) are defined [START_REF] Mazja | Estimates in Lp and Hölder classes and the Miranda-Agmon maximum principle for solutions of elliptic boundary value problems in domains with singular points on the boundary[END_REF] as the closure of C ∞ c (Ω \ O) (smooth functions vanishing in the vicinity of O) in the norm

Z; Λ l,α β (Ω) = l k=0 sup x∈Ω |x| β-l-α+k |∇ k x Z(x)| + sup x,y∈Ω,|x-y|< |x| 2 |x| β |x -y| -α |∇ l x Z(x) -∇ l y Z(y)|.
The standard norm in the Hölder space C l,α (Ω) looks as follows:

Z; C l,α (Ω) = l k=0 sup x∈Ω |∇ k x Z(x)| + sup x,y∈Ω,|x-y|< |x| 2 |x -y| -α |∇ l x Z(x) -∇ l y Z(y)|.
Here l ∈ {0, 1, ...}, α ∈ (0, 1) and β ∈ |.

Now we introduce several assumptions which are required to define the topological derivatives:

(H1) The limit problem [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF] has a solution v ∈ C 2,α (Ω) and F ∈ C 0,1 (Ω × |) with a certain α ∈ (0, 1).

(H2) The linear problem [START_REF] Il | in: Matching of asymptotic expansions of solutions of boundary value problems[END_REF] with

F ∈ C 0,α (Ω), g ∈ C 2,α (∂Ω) has a unique solution V ∈ C 2,α (Ω), V ; C 2,α (Ω) ≤ c( F; C 0,α (Ω) + g; C 2,α (∂Ω) ). ( 12
)
Here and in the sequel c stand for a positive constant that may change from place to place but never depends on ε.

(H3) F ′ v ∈ C 0,α (Ω × |).
If (H3) holds true and F ′ v (x, v(x)) ≤ 0 for x ∈ Ω, (H2) is also satisfied. Hypothesis (H2) means the existence and uniqueness of classical solutions to the linearized problem in Hölder spaces C 2,α (Ω) with the a priori estimate [START_REF] Jackowska-Strumillo | On numerical solution of shape inverse problems[END_REF]. It turns out that the linear mapping for problem [START_REF] Il | in: Matching of asymptotic expansions of solutions of boundary value problems[END_REF]

S : {F, g} -→ V, (13)
is an isomorphism in the Hölder spaces

C 2,α (Ω) → C 0,α (Ω) × C 2,α (∂Ω)
. By a general result in [START_REF] Mazja | Estimates in Lp and Hölder classes and the Miranda-Agmon maximum principle for solutions of elliptic boundary value problems in domains with singular points on the boundary[END_REF], (see also [START_REF] Nazarov | Elliptic Problems in Domains with Piecewise Smooth Boundaries[END_REF]; Ch. 3,4) the operator remains to be an isomorphism in weighted Hölder spaces under the proper choice of indices.

Theorem 1. Under assumptions (H2) and (H3), the mapping [START_REF] Kondratiev | Boundary problems for elliptic equations in domains with conical or angular points[END_REF] considered in the weighted Hölder spaces

S : Λ 0,α β (Ω) × C 2,α (∂Ω) -→ Λ 2,α β (Ω) is an isomorphism if and only if β -α ∈ (2, 3).
The following result on asymptotics is due to [START_REF] Kondratiev | Boundary problems for elliptic equations in domains with conical or angular points[END_REF][START_REF] Mazja | Estimates in Lp and Hölder classes and the Miranda-Agmon maximum principle for solutions of elliptic boundary value problems in domains with singular points on the boundary[END_REF] (see also [START_REF] Mazja | On the behavior of solutions to quasilinear elliptic boundary-value problems in a neighborhood of a conical point[END_REF] and, e.g., [START_REF] Nazarov | Elliptic Problems in Domains with Piecewise Smooth Boundaries[END_REF]; Ch. 3, 4). Theorem 2. If the right hand side in [START_REF] Il | in: Matching of asymptotic expansions of solutions of boundary value problems[END_REF] F ∈ Λ 0,α γ (Ω) and γ -α ∈ (1, 2), then the solution V to [START_REF] Il | in: Matching of asymptotic expansions of solutions of boundary value problems[END_REF] can be decomposed V (x) = V (x) + V (O) and the following estimate holds

|V (O)| + V ; Λ 2,α γ (Ω) ≤ c( F; Λ 0,α γ (Ω) + g; C 2,α (∂Ω) ). (14)
An assertion, similar to Theorem 1, is valid for the perforated domain Ω(ε) as well. The following result is due to [START_REF] Mazja | On the asymptotic behavior of solutions of elliptic boundary value problems with irregular perturbations of the domain[END_REF] (see also [START_REF] Mazja | Asymptotische Theorie elliptischer Randwertaufgaben in singulär gestörten[END_REF], Ch. 2.4, [START_REF] Nazarov | Elliptic Problems in Domains with Piecewise Smooth Boundaries[END_REF], Ch. 6).

Theorem 3. Under assumptions (H2) and (H3), the linearized problem

-∆ x v ε (x) -F ′ v (x, v(x))v ε (x) = F ε (x), x ∈ Ω(ε), v ε (x) = g ε (x), x ∈ ∂Ω(ε) ( 15
)
is uniquely solvable and the solution operator

S ε : {F ε , g ε } -→ v ε ( 16
)
is bounded in the weighted Hölder spaces

S ε : Λ 0,α β (Ω(ε)) × Λ 2,α β (∂Ω(ε)) -→ Λ 2,α β (Ω(ε)) . Moreover, in the case β -α ∈ (2, 3) the estimate v ε ; Λ 2,α β (Ω(ε)) ≤ c β ( F ε ; Λ 0,α β (Ω(ε)) + g ε ; Λ 2,α β (∂Ω(ε)) ) ( 17
)
is valid, where the constant c β is independent of ε ∈ (0, ε 0 ].

Remark Since |x| ≥ cε > 0 in Ω(ε), the weighted norm •; Λ 2,α β (Ω(ε)) is equivalent to the usual norm •; C 2,α (Ω(ε)) , however, the equivalence constants depend on ε. Thus Λ 2,α β (Ω(ε)) and C 2,α (Ω(ε)) coincide algebraically and topologically but are normed in different way. The norm of operator S ε is uniformly bounded for ε ∈ (0, ε 0 ] for any β, although the constant c β depends on ε provided β / ∈ (2, 3), that is, the norm of the inverse operator is uniformly bounded in ε ∈ (0, ε 1 ] only in the case of β ∈ (2, 3).

For the nonlinear problem (3), we shall use the classical solutions to boundary value-problem [START_REF] Casas | Uniform convergence of the FEM. Applications to state constrained control problems[END_REF], which means that for given F ∈ C 0,α (Ω × |), α ∈ (0, 1), the solutions lives in C 2,α (Ω), we refer to [START_REF] Ladyzhenskaya | tseva: Linear and quasilinear elliptic equations[END_REF] and [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF] for a result on the existence and uniqueness of solutions to semilinear elliptic boundary value-problems. It means, in particular, that problem (3) admits the unique solution u ε ∈ C 2,α (Ω(ε)) for some 0 < α < 1 and for all ε ∈ [0, ε 0 ).

Topological derivative for semilinear problems in 3D

We present here the complete analysis of the semilinear elliptic problem in three spatial dimensions. Such an analysis is interesting on its own, since in the existing literature there is no elementary derivation of the form of topological derivatives for nonlinear problems besides [START_REF] Mazja | On the asymptotic behavior of solutions of elliptic boundary value problems with irregular perturbations of the domain[END_REF], (see also [START_REF] Mazja | Asymptotische Theorie elliptischer Randwertaufgaben in singulär gestörten[END_REF],Ch. 5.7) i.e., using the asymptotic approximations of solutions to non linear PDE's. There are some results on topological derivatives of shape functional for non linear problems, see e.g., [START_REF] Amstutz | Topological sensitivity analysis for some nonlinear PDE system[END_REF], however such results are given in terms of the one term exterior approximation of the solutions and without asymptotically exact estimate.

Formal asymptotic analysis

Referring to [START_REF] Mazja | Asymptotische Theorie elliptischer Randwertaufgaben in singulär gestörten[END_REF], we set

u ε (x) = v(x) + w(ε -1 x) + εv ′ (x) + ..., (18)
where v, v ′ and w are components of regular and boundary layer types, respectively. Thus,

-∆ x v(x) -ε -2 ∆ ξ w(ξ) -ε∆ x v ′ (x) + ... = F (x, v(x) + w(ε -1 x) + εv ′ (x) + ...) = F (x, v(x)) + (w(ε -1 x) + εv ′ (x))F ′ v (x, v(x)) + ... (19)
In view of [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF], the first terms on the left and the right are cancelled and, moreover, w satisfies the problem [START_REF] Demlow | Sharply localized pointwise and W -1 ∞ estimates for finite element methods for quasilinear problems[END_REF] with

u(O) = v(O), -∆ ξ w(ξ) = 0, ξ ∈ | 3 \ ω, w(ξ) = -v(O), ξ ∈ ∂ω, (20)
while the boundary datum comes from the relation

v(x) + w(ε -1 x) + εv ′ (x) = v(O) + w(ε -1 x) + O(ε), x ∈ ∂ω ε .
We have

w(ξ) = -v(O)P (ξ) ( 21
)
where P is the capacitary potential [START_REF] Landkof | Fundamentals of modern potential theory. Izdat[END_REF][START_REF] Pólya | Isoperimetric inequalities in mathematical phisics[END_REF], e.g., a harmonic function in | 3 \ ω such that P (ξ) = 1 on ∂ω and

P (ξ) = |ξ| -1 cap(ω) + O(|ξ| -2 ), (22)
where cap(ω) is the capacity of the set ω. Since

w(ε -1 x) = -|x| -1 εv(O)cap(ω) + O(ε 2 |x| -2 ), (23)
we collect coefficients on ε in ( 19) and obtain

-∆ x v ′ (x) -v ′ (x)F ′ v (x, v(x)) = -aΦ(x)F ′ v (x, v(x)), x ∈ Ω, v ′ (x) = aΦ(x), x ∈ ∂Ω, (24)
where a = 4πv(O)cap(w) and Φ(x) = (4π|x|) -1 is the fundamental solution of the Laplace equation in | 3 .

Since a direct calculation yields

F ′ (•, v)Φ ∈ Λ 0,α γ (Ω) with any γ > 1 + α, we obtain the solution v ′ ∈ Λ 2,α β (Ω) of problem (24) such that v ′ -v ′ (O) ∈ Λ 2,α γ (Ω) where β -α ∈ (2, 3
) and γ -α ∈ (1, 2) can be taken arbitrarily.

Justification of asymptotic

We search for a solution of problem (3) in the form

u ε (x) = v(x) + w(ε -1 x) + εv ′ (x) + ûε (x) (25)
where ûε is a small remainder, satisfying the problem

   -∆ x ûε (x) = Fε (x; û), x ∈ Ω(ε), ûε (x) = ĝε Ω (x), x ∈ ∂Ω, ûε (x) = ĝε ω (x), x ∈ ∂ω(ε). (26)
Here

Fε (x; û) = F (x, v(x) + w(ε -1 x) + εv ′ (x) + ûε (x)) -F (x, v(x)) -ε(v ′ (x) -aΦ(x))F ′ v (x, v(x)), ĝε Ω (x) = -w(ε -1 x) -aεΦ(x), ĝε ω (x) = -v(x) + v(O) -εv ′ (x). (27)
We are going to employ the Banach contraction principle and, thus, we need to estimate the norms of [START_REF] Sokolowski | On topological derivative in shape optimization[END_REF].

Owing to ( 21), ( 22), the function x → w(ε -1 x) + aεΦ(x) is smooth on the surface ∂Ω, where |x| ≥ c > 0, and

|w(ε -1 x) + aεΦ(x)| ≤ |v(O)||P (ξ) -cap(ω)|ξ| -1 | ≤ cε 2 |x| -2 ≤ cε 2 , |∇ k x w(ε -1 x) + aε∇ k x Φ(x)| ≤ ε -k |v(O)||∇ k ξ P (ξ) -cap(ω)∇ k ξ |ξ| -1 | ≤ cε -k |ξ| -2-k = cε 2 |x| -2-k ≤ cε 2 . ( 28
)
Hence, by the above inequalities for the function x → w(ε -1 x) + aεΦ(x), we obtain the following estimates of the norm of ĝε Ω in the weighted Hölder space :

ĝε Ω ; Λ 2,α β (∂Ω) ≤ c ĝε Ω ; C 2,α (∂Ω) ≤ c ĝε Ω ; C 3 (∂Ω) ≤ cε 2 . (29)
Moreover, for β -β ′ > 0, we have

ĝε ω ; Λ 2,α β (∂Ω) ≤ c( sup x∈∂ωε 2 k=0 |x| β-2-α+k (|∇ k x (v(x) -v(O))| + ε|∇ k x v ′ (x)| + sup x,y∈∂ωε |x| β |x -y| -α (|∇ 2 x v(x) -∇ 2 y v(y)| + ε|∇ 2 x v ′ (x) -∇ 2 y v ′ (y)|)) ≤ c(ε β-1-α v; C 2,α (Ω) + ε 1+β-β ′ v ′ ; Λ 2,α β ′ (Ω) ). (30)
Notice that v ′ ∈ Λ 2,α β ′ (∂Ω) with arbitrary β ′ ∈ (2 + α, 3 + α), we shall further select the indices β and β ′ in an appropriate way. Let us denote

F(x, V (x)) = F (x, v(x) + V (x)) -F (x, v(x)) -V (x)F ′ v (x, v(x)). ( 31
) so that Fε (x; ûε ) = F(x, w(ε -1 x) + εv ′ (x) + ûε (x)) + (w(ε -1 x) + εaΦ(x) + ûε (x))F ′ v (x, v(x)). (32) Since (x → F ′ v (x, v(x))) ∈ C 0,α (Ω)
, by (H3), we take into account representation [START_REF] Nazarov | Asymptotic analysis of shape functional[END_REF] together with the inequality β -α > 2 and, as a result, we obtain

(w + εaΦ)F ′ v ; Λ 0,α β (Ω(ε)) ≤ c sup x∈Ω(ε) |x| β-α |x| ε -2 + sup x,y∈Ω(ε),|x-y|<|x|/2 |x| β |x -y| -α x ε - y ε |x| ε -3 ≤ cε 2 sup x∈Ω(ε) (|x| β-α |x| -2 + |x| β+1-α |x| -3 ) ≤ cε 2 .
(33) To estimate the first term on the right of (32), we need the following assumption on the mapping F: (H4) With a certain κ ∈ (0, 1), the inequality |F(x, V (x))| ≤ c|V (x)| 1+κ and the following relations are valid:

|F(x, V 1 (x)) -F(y, V 2 (y))| ≤ c(|x -y| α (|V 1 (x)| + |V 2 (y)|) 1+κ + |V 1 (x) -V 2 (y)|(|V 1 (x)| κ + |V 2 (y)| κ )), |F(x, V 1 (x)) -F(x, V 2 (x)) -(F(y, V 1 (y)) -F(y, V 2 (y)))| ≤ c|V 1 (x) -V 2 (x) -(V 1 (y) -V 2 (y))|(|V 1 (x)| κ + |V 2 (x)| κ + |V 1 (y)| κ + |V 2 (y)| κ ),
In other words, the mapping F enjoys the Hölder condition in both arguments and has a power-law growth in the second one; moreover, the second order difference satisfies estimate [START_REF] Kondratiev | Boundary problems for elliptic equations in domains with conical or angular points[END_REF].

Lemma 1. 1) Let V ∈ Λ 2,α
β (Ω(ε)) and β -α ∈ (2, 3), α ∈ (0, 1), κ ∈ (0, 1). Then, for x ∈ Ω(ε) and |x -y| < |x|/2, the estimates

|x| β-α |V (x)| 1+κ ≤ c V ; Λ 2,α β (Ω(ε)) 1+κ , |x| β-(β-α) κ 1+κ |x -y| -α |V (x) -V (y)| ≤ c V ; Λ 2,α β (Ω(ε)) are valid.
2) Under the same restrictions on α, β, κ and x, y as above

|x| β-α |w(ε -1 x)| 1+κ ≤ cε 1+κ , |x| β-(β-α) κ 1+κ |x -y| -α |w(ε -1 x) -w(ε -1 y)| 1+κ ≤ cε.
Proof. First, we readily show the first assertion :

|x| β-α |V (x)| 1+κ ≤ |x| β-α |x| -(1+κ)(β-2-α) (|x| β-2-α |V (x)|) 1+κ ≤ |x| 2-κ(β-2-α) V ; Λ 2,α β (Ω(ε)) 1+κ
The second inequality has followed from the relation

2 -κ(β -2 -α) ≥ 2 -1(3 -2 -α) > 1 > 0. Since 1 2 |x| < |y| < 3 2
|x|, in view of |x -y| < |x|/2 and using the Newton-Leibnitz formula, we conclude that

|x| β-(β-α) κ 1+κ |x -y| -α |V (x) -V (y)| 1+κ ≤ c|x| β-(β-α) κ 1+κ |x -y| -α |x| -β+1+α |x -y| sup x∈Ω(ε)) (|x| β-1-α |∇ x V (x)|) ≤ c|x| β-(β-α) κ 1+κ |x| 1-α |x| -β+1+α V ; Λ 2,α β (Ω(ε)) ,
while applying the inequalities

β -(β -α) κ 1 + κ + 1 -α -β + 1 + α = 2 -(β -α) κ 1 + κ ≥ 2 -(β -α -2)κ 1 + κ > 0.
Based on the assumptions β -α > 2 and 1 + κ < 2, we prove the second assertion. We have

|x| β-α |w(ε -1 x)| 1+κ ≤ c|x| β-α (1 + |x| ε ) -1-κ = cε 1+κ |x| β-α (ε + |x|) 1+κ ≤ cε 1+κ .
Owing to the estimate |P (ξ)| ≤ c(1 + |ξ|) -1 for the capacity potential and the boundary condition [START_REF] Nazarov | Self-adjoint extensions for the Neumann Laplace and applications[END_REF], it follows that

|x| β-(β-α) κ 1+κ |x -y| -α |w(ε -1 x) -w(ε -1 y)| ≤ c|x| β-(β-α) κ 1+κ |x -y| -α | x ε - y ε |(1 + |x| ε ) -2 sup ξ∈| 3 \ω (1 + |ξ|) 2 |∇ ξ w(ξ)|) ≤ cε|x| β-(β-α) κ 1+κ |x| 1-α (ε + |x|) -2 ≤ cε
Indeed, in the first inequality we have again applied the Newton-Leibnitz formula, in the second inequality we have used that

|∇ ξ P (ξ)| ≤ c(1 + |ξ|) -2 and β -α -(β -α) κ 1 + κ = β -α 1 + κ ≥ 1.
We now list the necessary estimates based on Lemma 1 and (H5). We start with the boundedness of the first term in (32) multiplied by a weight. We obtain

|x| β-α |F(x, w(ε -1 x) + εv ′ (x) + ûε (x))| ≤ c|x| β-α (|w(ε -1 x)| 1+κ + ε 1+κ |v ′ (x)| 1+κ + |û ε (x))| 1+κ ) ≤ c(ε 1+κ + ûε ; Λ 2,α β (Ω(ε)) 1+κ ). (34)
Second, we verify the boundedness of the weighted difference, namely,

|x| β |x -y| -α |F(x, =V (x) w(ε -1 x) + εv ′ (x) + ûε (x)) -F(y, V (y))| ≤ c(|x| β (|V (x)| 1+κ + |x -y| -α |V (x) -V (y)|(|V (x)| κ + |V (y)| κ )) ≤ c(ε 1+κ + ûε ; Λ 2,α β (Ω(ε)) 1+κ + (ε κ + ûε ; Λ 2,α β (Ω(ε)) κ )|x| β-(β-α) κ 1+κ |x -y| -α {|w( x ε ) -w( y ε )| + ε|v ′ (x) -v ′ (y)| + |û ε (x) -ûε (y)}) ≤ c(ε 1+κ + ûε ; Λ 2,α β (Ω(ε)) 1+κ ), (35)
Now, we deduce the local Lipschitz continuity of the first part of mapping (32):

|x| β-α |F(x, =V1(x) w(ε -1 x) + εv ′ (x) + ûε 1 (x)) -F(x, =V2(x) w(ε -1 x) + εv ′ (x) + ûε 2 (x))| ≤ c|x| β-α |û ε 1 (x) -ûε 2 (x)|(|V 1 (x)| κ + |V 2 (x)| κ ) ≤ c ûε 1 -ûε 2 ; Λ 2,α β (Ω(ε)) (ε κ + ûε 1 ; Λ 2,α β (Ω(ε)) + ûε 2 ; Λ 2,α β (Ω(ε)) ), (36)
Finally, we prove the local Lipschitz continuity for the weighted second order differences of the mapping F(., V i (.)), i = 1, 2. Indeed, we have

|x| β |x -y| -α |F(x, V 1 (x)) -F(x, V 2 (x)) -(F(y, V 1 (y)) -F(y, V 2 (y)))| ≤ c|x| β-(β-α) κ 1+κ |x -y| -α |(V 1 (x) -V 2 (x)) -(V 1 (y) -V 2 (y))|(ε κ + ûε 1 ; Λ 2,α β (Ω(ε)) + ûε 2 ; Λ 2,α β (Ω(ε)) ) ≤ c ûε 1 -ûε 2 ; Λ 2,α β (Ω(ε)) (ε κ + ûε 1 ; Λ 2,α β (Ω(ε)) κ ) (37)
The above estimates allow us to apply the Banach fixed point theorem to verify the existence of the remainder ûε . To this end, we rewrite problem [START_REF] Sokolowski | Introduction to shape optimization[END_REF] in the form of an abstract equation in the Banach space R = Λ 2,α β (Ω(ε)), namely,

ûε = Gû ε (38)
where Gû ε = S ǫ ( Fε (.; ûε ), ĝε Ω , ĝε ω ) and S ε denotes isomorphism [START_REF] Mazja | On the asymptotic behavior of solutions of elliptic boundary value problems with irregular perturbations of the domain[END_REF]. Let ûε belong to the ball B ⊂ R of radius Cε 1+κ . We further need to verify two properties. First, the mapping C maps the ball B into itself,

B ∋ ûε ⇒ Gû ε ∈ B, (39)
and, second, the mapping becomes a strict contraction on the ball, i.e.,

Gv -Gw; R ≤ k v -w; R , v, w ∈ R with k < 1. (40)
By ( 29), (30), (33) and (34), (35), we have

Gû ε ; R ≤ c( Fε ; Λ 0,α β (Ω(ε)) + ĝε Ω ; Λ 2,α β (Ω(ε)) + ĝε ω ; Λ 2,α β (Ω(ε)) ) ≤ c(ε 1+κ + ûε ; R 1+κ + ε 2 + ε β-1-α + ε 1+β-β ′ ). (41)
Let us fix β, α and β ′ , κ such that

(1, 2) ∋ β -α -1 ≥ 1 + κ, (42)
β -β ′ ≥ κ. (43)
Recall that β -α and β ′ -α belong to the interval (2, 3). Thus, to satisfy (43), we must put β -α near 3 (satisfying (42) as well) and β ′ -α near 2. This allows to create a gap of any length κ ∈ (0, 1). If (42) and (43) hold true, we obtain

Gû ε ; R ≤ c(4ε 1+κ + ûε ; R 1+κ ) ≤ Cε 1+κ , while the desired inequality C ≥ c(4 + C 1+κ ε (1+κ)κ
) is achieved by a proper choice of the constant C (e.g., C = 5c) and the bound for the parameter ε 0 in the condition ε ∈ (0, ε 0 ]. By virtue of (36) and (37), the estimate

Gv -Gw; R ≤ c(ε κ + 2C κ ε (1+κ)κ ) k v -w; R
is valid. The necessary relation k < 1 can be achieved by diminishing, if necessary, the upper bound ε 0 for ε again.

Theorem 4. Let the indices β, α and κ ∈ (0, 1) satisfy (42) and β -2 > κ, while (H2) and (H4) hold true. Then there exist positive constants C and ε 0 such that for ε ∈ (0, ε 0 ], the non-linear problem [START_REF] Sokolowski | Introduction to shape optimization[END_REF] has a unique small solution ûε , namely,

ûε ; Λ 2,α β (Ω(ε)) ≤ Cε 1+κ . ( 44
)
Consequently, the singularity perturbed problem (3) has at least one solution [START_REF] Sokolowski | Modelling of topological derivatives for contact problems[END_REF].

In the theorem we have proven the existence of a small remainder ûε in (25), i.e;, we have verified that problem (3) has a unique solution in a small ball centred at the approximate asymptotic solution. If the uniqueness of the solution ûε is known, for exemple, F in (3) gives rise to a monotone operator, the remainder is unique without any smallness assumption.

The formal asymptotic of the shape functional

We have

J (u ε ; Ω(ε)) = Ω(ε) J(x, v(x))dx + Ω(ε) (w(ε -1 x) + εv ′ (x))J ′ v (x, v(x))dx + ... = Ω J(x, v(x))dx + ε Ω (v ′ (x) -aΦ(x))J ′ v (x, v(x))dx + . . . (45)
We now introduce the following assumption:

(H5) J ∈ C 0,α (Ω × |), J ′ v ∈ C 0,α (Ω × |) Let p ∈ C 2,α (Ω) be a solution of the problem -∆ x p(x) -F ′ v (x, v(x))p(x) = J ′ v (x, v(x)), x ∈ Ω, p(x) = 0, x ∈ ∂Ω. (46)
Integrating by parts in

Ω \ ❇ δ = {x ∈ Ω : |x| > δ} yields Ω (v ′ (x) -aΦ(x))J ′ v (x, v(x))dx = -lim δ→0 Ω\❇ δ (∆ x p(x) + F ′ v (x, v(x))p(x))(v ′ (x) -aΦ(x))dx = -lim δ→0 Ω\❇ δ p(x)(∆ x + F ′ v (x, v(x)))(v ′ (x) -aΦ(x))dx -lim δ→0 ∂Ω ∂ n p(x)(v ′ (x) -aΦ(x))dx + lim δ→0 ∂❇ δ (∂ |x| p(x)(v ′ (x) -aΦ(x)) -p(x)∂ |x| (v ′ (x) -aΦ(x)))dx.
By [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérives partielles[END_REF], we have v ′ (x) -aΦ(x) = 0 for x ∈ Ω and

(∆ x + F ′ v (x, v(x)))(v ′ (x) -aΦ(x)) = ∆ x v ′ (x) + v ′ (x)F ′ v (x, v(x)) -aΦ(x)F ′ v (x, v(x)) = 0.
On the other hand,

∂ |x| p(x)(v ′ (x) -aΦ(x)) = O(δ -1
) and, hence,

Ω (v ′ (x) -aΦ(x))J ′ v (x, v(x))dx -lim δ→0 ∂❇ δ (∂ |x| p(x)(v ′ (x) -aΦ(x)) -p(x)∂ |x| (v ′ (x) -aΦ(x)))dx = -a lim δ→0 ∂❇ δ p(0)(4π|x| 2 ) -1 ds x = -ap(0) = -4πv(O)p(0)cap(ω).
Thus,

J (u ε ; Ω(ε)) = J (v; Ω) -ε4πv(O)p(0)cap(ω) + ... (47)
Let similarly to the first inequality in (H4) the following assumption be valid: (H6) With σ ∈ (0, 1),

|J(x, v(x) + V (x)) -J(x, v(x)) -V (x)J ′ v (x, v(x))| ≤ c|V (x)| 1+σ .
Using this assumption leads to the relation

|J (u ε ; Ω(ε)) -J (v; Ω(ε)) - Ω(ε) (w(ε -1 x) + εv ′ (x) + ûε (x))J ′ v (x, v(x))dx| ≤ c Ω(ε) |w(ε -1 x) + εv ′ (x) + ûε (x)| 1+σ dx ≤ c Ω(ε) x ε -1-σ + |x| -(1+σ)(β-2-α) × (ε 1+σ v ′ ; Λ 2,α β (Ω) 1+σ + ûε (x); Λ 2,α β (Ω) 1+σ ) dx ≤ cε 1+σ   1 ε r -1-σ r 2 dr + 1 ε r -(1+σ)(β-2-α) r 2 dr(ε 1+σ + ε (1+κ)(1+σ ) ≤ cε 1+σ .
Here we have taken into account that 1 + σ ≤ 2, (1 + σ)(β -2 -α) ≤ 2, and both the integrals, extended on the interval (0, 1), do converge. It suffices to mention the following inequalities:

|J (u ε ; Ω(ε)) -J (v; Ω)| ≤ c mes 3 (ω ε ) ≤ cε 3 , Ω(ε) |w(ε -1 x) + aεΦ(x)||J ′ v (x, v(x))|dx ≤ c 1 ε r ε -2 rdr ≤ cε 2 , Ω(ε) |û ε ||J ′ v (x, v(x))|dx ≤ cε 1+κ Ω(ε) |x| -(β-2-α) dx ≤ cε 1+κ .
This confirms the formal calculations performed in section 2.2. Let us formulate the main result in three dimensions.

Theorem 5. Under the assumptions listed above, we have

|J (u ε ; Ω(ε)) -J (v; Ω) + ε4πv(O)p(0)cap(ω)| ≤ cε 1+min(σ,κ) .
3 Topological derivative for mixed semilinear elliptic problem in two spatial dimensions

The numerical analysis is performed in two spatial dimensions. Therefore, we introduce a mixed semilinear problem and analyze the asymptotic in such a case. Since the proof uses the same arguments as in three spatial dimensions, we provide only the formal analysis and impose the Neumann boundary conditions on the hole boundary ∂ω ε . Note that the Dirichlet condition on ∂ω ε changes crutially the form of asymptotic expansions (cf. [START_REF] Il | in: Matching of asymptotic expansions of solutions of boundary value problems[END_REF] and [START_REF] Mazja | On the asymptotic behavior of solutions of elliptic boundary value problems with irregular perturbations of the domain[END_REF], [ [START_REF] Mazja | Asymptotische Theorie elliptischer Randwertaufgaben in singulär gestörten[END_REF]; Ch. 5.7]).

Formal asymptotic analysis

Let Ω and ω be bounded domains in the plane | 2 . We consider the nonlinear mixed problem in the singularly perturbed domain Ω(ε), defined in (2):

   -∆ x u ε (x) = F (x, u ε (x)), x ∈ Ω(ε), u ε (x) = 0, x ∈ ∂Ω, ∂ n u ε (x) = 0, x ∈ ∂ω ε . ( 48
)
Referring to [START_REF] Il | in: Matching of asymptotic expansions of solutions of boundary value problems[END_REF], [START_REF] Mazja | Asymptotische Theorie elliptischer Randwertaufgaben in singulär gestörten[END_REF], especially to [START_REF] Mazja | On the asymptotic behavior of solutions of elliptic boundary value problems with irregular perturbations of the domain[END_REF] and [[17] § 5.7], we set

u ε (x) = v(x) + εw 1 (ε -1 x) + ε 2 w 2 (ε -1 (x) + ε 2 v ′ (x) + ..., (49)
where v, v ′ and w 1 , w 2 are component of regular and boundary layer types, respectively. Precisely, v is a smooth solution of problem [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF] in the two dimensional entire domain Ω. The Taylor formula yields

v(x) = v(O) + x T ∇ x v(O) + 1 2 x T ∇ 2 x v(O)x + O(|x| 3 ).
The second term w 1 in the symmetric ansatz (49) becomes a solution of the exterior problem

-∆ ξ w 1 (ξ) = 0, ξ ∈ | 2 \ ω, ∂ n(ξ) w 1 (ξ) = -∂ n(ξ) ξ T ∇ x v(O), ξ ∈ ∂ω, (50)
Such a solution admits the asymptotic representation

w 1 (ξ) = - 1 2π ξ T |ξ| 2 m(ω)∇ x v(O) + O(|ξ| -2 ), |ξ| → ∞,
where m denotes the virtual mass matrix, see [START_REF] Pólya | Isoperimetric inequalities in mathematical phisics[END_REF]. Then the third term w 2 in (32) satisfies the problem

-∆ ξ w 2 (ξ) = 0, ξ ∈ | 2 \ ω, ∂ n(ξ) w 2 (ξ) = -∂ n(ξ) 1 2 ξ T ∇ x v(O)ξ, ξ ∈ ∂ω, (51)
For such a solution, we write down the classical asymptotic representation

w 2 (ξ) = c 2π ln 1 |ξ| + O( 1 |ξ| ) |ξ| → ∞,
where the constant c in the above equality can be calculated as follows:

∂ω ∂ n(ξ) w 2 (ξ)ds ξ = - ∂❇ R ∂ ∂|ξ| c 2π ln 1 |ξ| dξ = c. (52)
By the Green formula, we compute the left boundary integral

- ∂ω ∂ n(ξ) 1 2 ξ T ∇ 2 x v(O)ξds ξ = ω ∆ ξ 1 2 ξ T ∇ 2 x v(O)ξdξ = mes 2 ω∆ x v(O) = -mes 2 ωF (O; v(O)). (53)
Finally, the fourth term v ′ in (32) is to be found from the Dirichlet problem

       -∆ x v ′ (x) = (- 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x))F ′ v (x, v(x)), x ∈ Ω, v ′ (x) = 1 2π x T |x| 2 m(ω)∇ x v(O) + 1 2π ln ε |x| mes 2 ωF (O; v(O), x ∈ ∂Ω, (54)

The formal asymptotics of the shape functional

We introduce the following hypotheses:

(H7) F ∈ C 0,α (Ω × |), F ′ v ∈ C 0,α (Ω × |) for some α ∈ (0, 1) and F ′ v ≤ 0. (H8) J ∈ C 0,α (Ω × |), J ′ v ∈ C 0,α (Ω × |)
By the monotonicity of F , the Lax-Milgram Lemma and the regularity of J, the problem

-∆ x p(x) -F ′ v (x, v(x))p(x) = J ′ v (x, v(x)), x ∈ Ω, p(x) = 0, x ∈ ∂Ω. ( 55
)
admits a unique solution p ∈ C 2,α (Ω).

We replace the solution u ε by its asymptotic representation (32). As a result we obtain the first asymptotic term of order ε 2 for the shape functional

J (u ε ; Ω(ε)) = Ω(ε) J(x, v(x))dx + Ω(ε) (εw 1 (ε -1 x) + ε 2 w 2 (ε -1 x) + ε 2 v ′ (x))J ′ v (x, v(x))dx + ... = J (v; Ω(ε)) + ε 2 Ω(ε) - 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x))J ′ v (x, v(x) dx + ... = J (v; Ω) -ε 2 mes 2 ωJ(O; v(O)) +ε 2 Ω - 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x) J ′ v (x, v(x))dx + ... ( 56
)
Now we replace the right-hand side of (55 ) according to the equation and twice integrate by parts in the domain

Ω \ ❇ δ = {x ∈ Ω : |x| > δ}. We have Ω (- 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x))J ′ v (x, v(x))dx = -lim δ→0 Ω\❇ δ (∆ x p(x) + F ′ v (x, v(x))p(x))(- 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x))dx = -lim δ→0 Ω\❇ δ p(x)(∆ x + F ′ v (x, v(x)))(- 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x))dx -lim δ→0 ∂Ω ∂ n p(x)(- 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x))dx -lim δ→0 ∂❇ δ (∂ |x| p(x)(- 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x)) -p(x)∂ |x| (- 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x))dx.
On the other hand, the boundary condition (54) implies that

- 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x) = 0. (57)
Furthermore, for the linearized operator ∆ x + F ′ v , the inequality

(∆ x + F ′ v (x, v(x))) + - 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x) = ∆ x v ′ (x) + v ′ (x)F ′ v (x, v(x)) + - 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) F ′ v (x, v(x)) = 0 is valid because the function x → - 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)
is harmonic. Hence, we obtain that

Ω - 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) + v ′ (x))J ′ v (x, v(x) dx -lim δ→0 ∂❇ δ ∂ |x| p(x) - 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O)) -p(x)∂ |x| - 1 2π x T |x| 2 m(ω)∇ x v(O) - 1 2π ln ε |x| mes 2 ωF (O; v(O) dx -lim δ→0 2 i,j,k=1 ∂❇ δ x i |x| ∂p ∂x i (O) - 1 2π x k |x| 2 m kj ∂v ∂x j (O) -(p(O) + x i ∂p ∂x i (O)) - 1 2π x j |x| 3 m jk ∂v ∂x k (O) dx + lim δ→0 ∂❇ δ 2 i=1 (p(O) + x i ∂p ∂x i (O)) × 1 2π 1 |x| mes 2 ωF (O; v(O)) dx -lim δ→0 ∂❇ δ   2 i,j,k=1 x i |x| ∂p ∂x i (O) - 1 π x k |x| 2 m kj ∂v ∂x j (O) -p(O) 1 2π 1 |x| mes 2 ωF (O; v(O)) dx lim δ→0 ∂❇ δ   2 i,j=1 x i |x| ∂p ∂x i (O) 1 π x i |x| 2 m ij ∂v ∂x j (O) +p(O) 1 2π 1 |x| mes 2 ωF (O; v(O)) dx = F (O; v(O))mes 2 ωp(O) + ∇ x p(O) T m(ω)∇ x v(O).
Thus, recalling (56) we conclude that the relation

J (u ε ; Ω(ε)) = J (v; Ω) + ε 2 [-mes 2 ωJ(O; v(O)) + F (O; v(O))mes 2 ωp(O) + ∇ x p(O) T m(ω)∇ x v(O)] + ... (58)
confirms the following result: Theorem 6. Under the assumptions (H1), (H7) and (H8), the following asymptotic expansion is valid :

J (u ε ; Ω(ε)) = J (v; Ω) + ε 2 [-mes 2 ωJ(O; v(O)) + F (O; v(O))mes 2 ωp(O) + ∇ x p(O) T m(ω)∇ x v(O)] + ...

Finite element approximations of topological derivatives

Our aim in this section is to compute a numerical approximation of the topological derivative of the shape functional [START_REF] Fulmanski | A level set method in shape and topology optimization for variational inequalities[END_REF], with u ε the solution of the problem (48) and give L ∞ -estimates of the error.

Family of finite elements

In Ω we consider a family of triangulations {T h } h>0 . With each element T ∈ T h , we associate two parameters ρ(T ) and σ(T ), where ρ(T ) denotes the diameter of the set T , and σ(T ) is the diameter of the largest ball contained in T . We set h = max T ∈T h ρ(T ). We make the following assumptions on the triangulations.

(H10) Regularity assumption : there exists σ > 0 such that ρ(T ) σ(T ) ≤ σ for T ∈ T h and h > 0.

(H11) Inverse assumption : there exists ρ > 0 such that h ρ(T )

≤ ρ for T ∈ T h and h > 0.

(H12) We denote by Ω h = ∪ T ∈T h T the domain obtained by a triangulation, with Ω h its interior and with ∂Ω h its boundary. Then we assume that the vertexes of T h placed on the boundary ∂Ω h belongs also to ∂Ω.

Consider the spaces

V h = v h ∈ C(Ω) : v h | T ∈ P 1 (T ) for T ∈ T h and v h = 0 in Ω \ Ω h } and W h = v h ∈ C(Ω h ) : v h | T ∈ P 1 (T ) for T ∈ T h ,
where P 1 (T ) is the space of polynomials of degree 1 on T . V h is a vector subspace of H 1 0 (Ω) and W h is a subspace of H 1 (Ω).

We use the Lagrange interpolation operator

Π h : C(Ω) → W h
being Π h z the unique element in W h such that Π h z(x i ) = z(x i ) for every node x i of the triangulation. In the case of a function z vanishing on ∂Ω, we extended Π h z to Ω by zero and we denote this extension by Π h z too. In the last case we have that Π h z ∈ V h .

Numerical solution of semilinear problem

By vitrue of the assumption F ∈ C 0,1 (Ω × |) and by the mean value theorem, we deduce the following local Lipschitz condition: for all M > 0 there exists c M > 0 such that

|F (x, v 1 ) -F (x, v 2 )| ≤ c M |v 1 -v 2 |, x ∈ Ω, |v 1 |, |v 2 | ≤ M. (59)
Using classical arguments we can deduce from the monotonicity of F (x, .) and ( 59), the existence of a unique solution to [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF] in H 1 0 (Ω) ∩ C(Ω); we refer to [START_REF] Stampacchia | Le problème de Dirichlet pour les équations elliptiques du second ordre à coefficients discontinus[END_REF] for the boundedness of the solution. From to the convexity of Ω, we can deduce that the solution is in H 2 (Ω) (see [START_REF] Ladyzhenskaya | tseva: Linear and quasilinear elliptic equations[END_REF]).

The weak formulation of the equation ( 7) is the following

a(v, z) = (F (x, v), z) L 2 (Ω) , z ∈ H 1 0 (Ω), (60)
where

a(v, z) = Ω ∇v(x)∇z(x)dx.
The numerical approximation v h of v is then the solution of the problem:

F ind v h ∈ V h such that a(v h , z h ) = (F (x, v h ), z h ) L 2 (Ω) , for any z h ∈ V h . ( 61
)
The proof of the existence of a solution of (61) is well known (see [START_REF] Stampacchia | Le problème de Dirichlet pour les équations elliptiques du second ordre à coefficients discontinus[END_REF]). It is enough to apply, in a convenient way, Browder's fixed point theorem along with the monotonicity of F (x, .).

We rewrite problem (55 ) in the form

-∆p(x) = F 0 (x, p(x)), x ∈ Ω, p(x) = 0, x ∈ ∂Ω, (62)
where

F 0 (x, p(x)) = F ′ v (x, v(x))p(x)+J ′ (x, v(x)
). F 0 is linear with respect to the second variable, then F 0 ∈ C 0,1 (Ω×

|).

Then the variational formulation of the linear problem (55 ) is the following,

a(p, z) = (F 0 (x, p), z) L 2 (Ω) , z ∈ H 1 0 (Ω). (63)
Thus, that numerical approximation of p is the solution of the variational problem

F ind p h ∈ V h such that a(p h , z h ) = (F 0 (x, p h ), z h ) L 2 (Ω) for any z h ∈ V h . (64)
Due to the hypothesis (H1) and (H7) the linear problem (64) have a solution in space V h . In the next section we give L ∞ -estimates for the finite element approximation of v, ∇v, p, ∇p and the topological derivatives in the case of semilinear elliptic problem.

4.3 Convergence of finite element approximation for semilinear problem [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF] and for adjoint state problem (55)

We are going to use the recent results on the convergence of finite element method in W 1,∞ spaces. The topological derivative is a point wise expression depending on the values of the state equation solution, the adjoint state, as well as the value of the gradient of these functions. In order to derive the error estimates in the L ∞ norm for the topological derivative it is required to have in hand the error estimates in the W 1,∞ norm for the solutions of the semilinear equation as well as the linear adjoint state equation. The results presented below lead to the error estimate for the topological derivative. We refer the reader to [START_REF] Casas | Uniform convergence of the FEM. Applications to state constrained control problems[END_REF], [START_REF] Demlow | Sharply localized pointwise and W -1 ∞ estimates for finite element methods for quasilinear problems[END_REF], [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF], for the proofs of the error estimates for linear and semilinear elliptic equations.

The following L ∞ -estimate for the approximation by finite elements of solutions to problem [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF] was proved in [START_REF] Casas | Uniform convergence of the FEM. Applications to state constrained control problems[END_REF].

Theorem 7. Let v and v h be solutions of the variational problems ( 60) and (61), respectively. Then there exists a constant

C > 0 independent of h such that v -v h ; L ∞ (Ω h ) ≤ Ch v; H 2 (Ω) .
In addition, assuming that Ω h ⊂ Ω, we have the following estimates proved in [START_REF] Demlow | Sharply localized pointwise and W -1 ∞ estimates for finite element methods for quasilinear problems[END_REF] (see also [START_REF] Frehse | Asymptotic L ∞ -error estimates for linear finite element approximations of quasilinear boundary value problems[END_REF])

Theorem 8. Let v and v h be solutions of the variational problems ( 60) and ( 61), respectively. Then there exists a constant C > 0 independent of h such that

∇v -∇v h ; L ∞ (Ω h ) ∞ ≤ Ch v; H 2 (Ω) .
On the other hand, we have the following L ∞ -estimates for the approximation of solutions to the linear problem (55 ): Theorem 9. Let x 0 ∈ Ω and let p and p h be solutions of the variational problems ( 63) and ( 64) respectively. Then there exists a constant C > 0 independent of h such that

p -p h ; L ∞ (Ω h ) + ∇p -∇p h ; L ∞ (Ω h ) ≤ Ch p; H 2 (Ω) .

L ∞ -estimates for the approximation of the topological derivative

We denote by T Ω,h the numerical approximation by finite element method of the topological derivative T Ω . Then

T Ω,h (O) = -mes 2 (ω)J(O, v h (O)) + ∇ x p h (O) T m(ω)∇ x v h (O) + F (O, v h (O))mes 2 (ω)p h (O). (65)
We obtain

|T Ω (O) -T Ω,h (O)| = | -mes 2 (ω)[J(O, v(O)) -J(O, v h (O))] +[∇ x p(O) T m(ω)∇ x v(O) -∇ x p h (O) T m(ω)∇ x v h (O)] +[F (O, v(O))mes 2 (ω)p(O) -F (O, v h (O))mes 2 (ω)p h (O)]| ≤ mes 2 (ω)|J(O, v(O)) -J(O, v h (O))| +|∇ x p(O) T m(ω)∇ x v(O) -∇ x p h (O) T m(ω)∇ x v h (O)| +mes 2 (ω)|F (O, v(O))p(O) -F (O, v h (O))p h (O)|. (66)
We have

|J(O, v(O)) -J(O, v h (O))| ≤ c v -v h ; L ∞ (Ω h ) . (67)
It follows that

|T Ω (O) -T Ω,h (O)| ≤ c v -v h ; L ∞ (Ω h ) + |∇ x p(O) T m(ω)∇ x (v(O) -v h (O))| +|∇ x (p(O) -p h (O)) T m(ω)∇ x v h (O)| + mes 2 (ω)|F (O, v(O))(p(O) -p h (O))| +mes 2 (ω)|(F (O, v(O)) -F (O, v h (O)))p h (O)|. (68)
Thanks to (54) we obtain

|F (O, v(O)) -F (O, v h (O))| ≤ c 2 v -v h ; L ∞ (Ω h ) , (69)
and therefore

|T Ω (O) -T Ω,h (O)| ≤ c 1 v -v h ; L ∞ (Ω h ) + c 2 ∇ x (v(O) -v h (O)); L ∞ (Ω h ) +c 3 ∇ x (p(O) -p h (O)); L ∞ (Ω h ) + c 4 p(O) -p h (O); L ∞ (Ω h ) . (70)
Finally, by the Theorems 8 and 9, we deduce the following result.

Theorem 10. The following error estimate hold for the approximation of the topological derivatives

|T Ω (O) -T Ω,h (O)| ≤ Ch( v; H 2 (Ω) + p; H 2 (Ω) ). (71)

Numerical examples

In this section we presents some numerical examples to show us the behavior of topological derivative approximation w.r.t. the evolution of discretization step size. We derive errors and verifies that the computed error satisfy the Theorem (10) estimate in each case.

For each of the examples we choose the domain Ω as a square (0, 1) × (0, 1) and the following energy functional

J (v; Ω) = 1 2 Ω |∇v(x)| 2 dx + 1 4 Ω v 4 (x)dx - Ω f (x)v(x)dx ( 72
)
where v is the solution of the nonlinear problem

-∆ x v(x) = -v(x) 3 + f (x), x ∈ Ω, v(x) = 0, x ∈ ∂Ω. ( 73
)
We take different right-hand side f in each of the examples. The size of discretization is determined by h, which decrease in each iteration. We compute the error in 20 iteration strarting with h = 0.2 and reduce it by 0.01 in each step.

Example 1. In the first example the function f is given by

f (x) = f (x 1 , x 2 ) = (x 2 1 + x 2
2 ) sin πx 1 sin πx 2 3 + 2 π 2 (x 2 1 + x 2 2 ) -2 sin πx 1 sin πx 2 -4π (x 1 cos πx 1 sin πx 2 + sin πx 1 cos πx 2 )

We calculate the exact solution of (73) The exact and numerical approximation of the topological derivative is presented in Figure 1 . In Figure 2 we plot the relative evolution of the error and the behavior of the error w.r.t. the discretization step size. In this case u(x 1 , x 2 ) = x 1 (x 1 -1)x 2 (x 2 -1) The exact topological derivative and it numerical approximation are presented in Figure 3 and in Figure 4 we show the relative evolution of the error and the behavior of the error w.r.t. the discretization step size.

Example 3. In the last example we take f (x 1 , x 2 ) = 100x 2 1 x 2 2 (x 1 -1) 2 (x 2 -1) 2 3 -200 x 2 1 (x 1 -1) 2 (6x 2 2 -6x 2 + 1) + x 2 2 (x 2 -1) 2 (6x 2 1 -6x 1 + 1) then u(x 1 , x 2 ) = 100x 2 1 x 2 2 (x 1 -1) 2 (x 2 -1) 2 , and p(x 1 , x 2 ) = -50x 2 1 x 2 2 (x 1 -1) 2 (x 2 -1) 2 . In this last example the exact topological derivative and it numerical approximation are presented in Figure 5 and in Figure 6 we show the relative evolution of the error and the behavior of the error w.r.t. the discretization step size.

Conclusion

In the paper the form of topological derivatives of integral shape functional is obtained for semilinear elliptic boundary value problems in two and three spatial dimensions. The finite element method is used to compute an approximation of the topological derivatives. The convergence analysis of the finite element method is performed in two spatial dimensions. Since the application of topological derivatives in shape optimization requires pointwise values we provide L ∞ -estimates for finite element approximations. The result of computations confirm the apriori estimate of the numerical approximation error. The presented results can be used for shape and topology optimization for semilinear elliptic boundary value problems.

u(x 1 x 2 1 + x 2 2

 122 , x 2 ) = x 2 1 + x 2 2 sin πx 1 sin πx 2 and the corresponding adjoint state p(x 1 , x 2 ) = -1 2 sin πx 1 sin πx 2

Figure 1 :

 1 Figure 1: Topological derivative: exact T Ω (left) and approximate T Ω,h (right).

Example 2 .

 2 Let us chosef (x) = f (x 1 , x 2 ) = (x 1 (x 1 -1)x 2 (x 2 -1)) 3 -2 (x 2 (x 2 -1) + x 1 (x 1 -1))

Figure 2 :

 2 Figure 2: Error :|T Ω -T Ω,h | ∞ /|T Ω | ∞ .

Figure 3 :

 3 Figure 3: Topological derivative: exact T Ω (left) and approximate T Ω,h (right).

Figure 4 :

 4 Figure 4: Error : |TΩ -T Ω,h | ∞ /|T Ω | ∞ .

Figure 5 :

 5 Figure 5: Topological derivative: exact T Ω (left) and approximate T Ω,h (right).

Figure 6 :

 6 Figure 6: Error : |TΩ -T Ω,h | ∞ /|T Ω | ∞ .

Authors' addresses: Mohamed Iguernane, I.E.C.N. Université Henri Poincaré, B.P 239, 54506 Vandoeuvre-lès-Nancy,France e-mail: iguernan@iecn.u-nancy.fr. Serguei A. Nazarov, Institute of Mechanical Engineering Problems of Machine, Russian Academy of Sciences, St. Petersburg 199178 e-mail: serna@snark.ipme.ru, srgnazarov@yahoo.co.uk. Jean-Rodolphe Roche, I.E.C.N. Université Henri Poincaré, B.P 239, 54506 Vandoeuvre-lès-Nancy,France e-mail: roche@iecn.u-nancy.fr.

Jan Sokolowski, I.E.C.N. Université Henri Poincaré, B.P 239, 54506 Vandoeuvre-lès-Nancy,France e-mail: Jan.Sokolowski@iecn.u-nancy.f. Katarzyna Szulc, I.E.C.N. Université Henri Poincaré, B.P 239, 54506 Vandoeuvre-lès-Nancy,France e-mail: Katarzyna.Szulc@iecn.u-nancy.fr.