\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{amssymb,graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2018 (2018), No. 146, pp. 1--18.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2018 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2018/146\hfil Shape derivative]
{A shape-derivative approach to some PDE model in image restoration}

\author[C. Baroncini, J. Fern\'andez Bonder \hfil EJDE-2018/146\hfilneg]
{Carla Baroncini, Juli\'an Fern\'andez Bonder}

\address{Carla Baroncini \newline
IMAS - CONICET and Departamento de Matem\'atica,
FCEyN - Universidad de Buenos Aires,
Ciudad Universitaria, Pabell\'on I  (1428)
Buenos Aires, Argentina}
\email{cbaronci@dm.uba.ar}

\address{Juli\'an Fern\'andez Bonder \newline
IMAS - CONICET and Departamento de Matem\'atica,
FCEyN - Universidad de Buenos Aires,
Ciudad Universitaria, Pabell\'on I  (1428)
Buenos Aires, Argentina}
\email{jfbonder@dm.uba.ar}
\urladdr{http://mate.dm.uba.ar/~jfbonder}


\dedicatory{Communicated by Claudianor O. Alves}

\thanks{Submitted February 26, 2018. Published July 25, 2018.}
\subjclass[2010]{49Q10, 49J45}
\keywords{Shape optimization; sensitivity analysis; nonstandard growth}

\begin{abstract}
 In this paper we analyze the shape derivative of a cost functional appearing
 in image restoration. The main feature of this cost functional is the
 appearance of a variable exponent.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

Shape derivative (or Hadamard derivative) has been proved to be a valuable
tool in order to study shape optimization problems. The main ideas go back
to Hadamard's original paper \cite{Hadamard} and has been further developed since.
See for instance the books \cite{Banichuk, Pironneau, Sokolowski-Zolesio}.

In this article we are devoted to the analysis of shape derivative of certain
functionals arising in image restoration, whose main feature is that it
involves a variable exponent.

Let us begin by discussing the model where these functionals appear.
The goal in image restoration is to obtain an image which is modeled by
a function $u\colon \Omega \to \mathbb R $, where
$\Omega=(0,1)\times(0,1)\subset \mathbb R ^2$, given that one has a distorted image
$I\colon \Omega \to \mathbb R $.

It is customary to assume that the introduced error, $e=u-I$, is small and
the objective is to recover $u$ from $I$ without making any further assumptions
on the error $e$.
A classical PDE model introduced by Chambolle and Lions in \cite{Chambolle-Lions}
in 1997, propose to obtain $u$ as a minimizer of the functional
$$
\min \frac{1}{2\beta} \Big(\int_{\{|\nabla v|\le \beta\}} |\nabla v|^2\,dx
+ \int_{\{|\nabla v|>\beta\}} |\nabla v|\,dx\Big)
+ \frac{\beta}{2} \int_{\Omega} (v-I)^2\,dx,
$$
where $\beta>0$ is a parameter that needs to be adjusted by the operator of
the method for each image.
The idea behind this method is that the real image must be smooth in regions
where there are no boundaries (which are interpreted as regions where the
derivatives are not big) and, in the ones which contains boundaries,
the solution must admit discontinuities.
This method can be re-written as follows
$$
\min \frac{1}{2\beta} \int_\Omega |\nabla v|^{p(|\nabla v|)}\,dx
 + \frac{\beta}{2}\int_\Omega (v-I)^2\,dx,
$$
where the exponent $p$ is defined as
$$
p(t) = \begin{cases}
2 & \text{if } t\le \beta\\
1 & \text{if } t>\beta.
\end{cases}
$$
This method is extremely difficult to study rigorously since the space where
the functional is defined is not a good functional space.
That is why, in 2006, Chen, Levine and Rao introduced in \cite{CLR} a
modification by which the exponent $p$ is computed from $I$ but it is fixed.
In this second model,
$$
p(x) = 1 + \frac{1}{1+k|\nabla G_\sigma * I|^2},
$$
where $G_\sigma(x) = \frac{1}{\sigma} \exp(-|x|^2/4\sigma^2)$ is the Gaussian
filter, with $k, \sigma>0$ parameters. Therefore, $p\sim 1$ where $I$ is
discontinuous and $p\sim 2$ where $I$ is smooth.

Then, the problem to be minimized is
$$
\min \frac{1}{2\beta} \int_\Omega |\nabla v|^{p(x)}\,dx
+ \frac{\beta}{2}\int_\Omega (v-I)^2\,dx.
$$
By considering a fixed regular exponent, the authors can use the Sobolev
and Lebesgue spaces with variable exponent, thoroughly studied since the sixties.
See \cite{Diening}.

Here we consider a variant of these methods, that can be thought of being
in between these two,  that approximates the one created by Chambolle and
Lions preserving the good functional properties given by the one presented
by Chen, Levine and Rao.

We start by dividing the region $\Omega$ into two sub regions $D_1$ and $D_2$
such that for $i=1,2$,
\begin{equation}\label{prop.particion}
D_i \subset \Omega \text{ is open},\ \mathring{\overline{D_i}}
= D_i,\  D_1\cap D_2 = \emptyset,\  \text{and}\ \overline{\Omega}
= \overline{D_1}\cup \overline{D_2}.
\end{equation}

By this partition, we make sure that $D_1$ contains the regions with boundaries
of the image and $D_2$ its complement. One way of creating this partition
is the following:
$$
D_1 = \{x\in\Omega \colon |\nabla G_\sigma * I|> \beta\},\quad
D_2 =\{x\in \Omega \colon |\nabla G_\sigma * I|<\beta\}.
$$

We define an exponent  $p\colon \Omega\to \mathbb R $ given by
$$
p(x) = \begin{cases}
1+\epsilon & \text{if } x\in D_1\\
2 & \text{if } x\in D_2.
\end{cases}
$$
Then we compute $u$ by minimizing the functional
$$
J(v) = \frac{1}{2\beta}\int_\Omega |\nabla v|^{p(x)}\,dx
 + \frac{\beta}{2} \int_\Omega (v-I)^2\,dx.
$$
To improve the image found, we then may apply an iterative
\textit{steepest descent type method} by following the shape derivative of
the functional.
So the main objective of this paper is to compute this shape derivative.


Let us recall that  a related minimization problem was studied in
\cite{Acerbi-Fusco}. In that article it is shown that minimizers
are H\"older-continuous across the interphase.
Hence we are left with the problem of computing the shape derivative of $J(u)$
with respect to $D_i$, which we describe now.
Given $V\colon \mathbb R ^N \to \mathbb R ^N$ a Lipschitz deformation field,
the associated flow $\{\Phi_t\}_{t\in\mathbb R }$ is defined by
\begin{equation}\label{flux}
\begin{gathered}
\frac{d}{dt}\Phi_{t}(x) = V(\Phi_{t}(x)),\quad t\in\mathbb R ,\; x\in \mathbb R ^N,\\
\Phi_{0}(x) = x \quad x\in\mathbb R ^N.
\end{gathered}
\end{equation}
Let us observe that $\Phi_t\colon \mathbb R ^N \to \mathbb R ^N$ is a group of diffeomorfisms.
That is, $\Phi_t\circ \Phi_s = \Phi_{t+s}$ and $\Phi_t^{-1} = \Phi_{-t}$.

We will assume that $\operatorname{suppt}(V)\subset \Omega$, so that $\Phi_t(\Omega) = \Omega$
for every $t\in \mathbb R $.
Then, the regions $D_i$ are deformed by $\Phi_{t}$ and we obtain a family of
 partitions $D_i^t = \Phi_t(D_i),\quad i=1,2$ that verify \eqref{prop.particion}
and we define
$$
p_t(x) = \begin{cases}
1+\epsilon & \text{if } x\in D_1^t\\
2 & \text{if } x\in D_2^t.
\end{cases}
$$
Observe that $p_t = p\circ \Phi_{-t}$.

Then, for each $t\in\mathbb R $ we define the  functional
$$
J_t(v) = \frac{1}{2\beta}\int_\Omega |\nabla v|^{p_t(x)}\,dx
 + \frac{\beta}{2} \int_\Omega (v-I)^2\,dx,
$$
Let $u_t$ be the minimizer of $J_t$. We can consider the function
$j\colon \mathbb R \to\mathbb R $ given by $j(t) = J_t(u_t)$.
The shape derivative consists then in computing $j'(0)$.

Then, by finding a good expression for such derivative, it will be possible
to compute the deformations field $V$ which makes it as negative as possible
and so choose the optimal deformation field to then iterate
$$
D_i^{\Delta t} \simeq (id + \Delta t V)(D_i).
$$

This article is organized as follows:
In section 2 we collect the preliminaries con variable exponent Lebesgue
and Sobolev spaces that will be used in the article. With the exception
of Proposition \ref{prop.24} there is no new material so the expert reader
can safely skip this section and move on to the rest of the article.
In section 3 we prove one of the mail results of the paper, namely the
differentiability of the cost functional $j(t) = J_t(u_t)$ (Theorem \ref{teo.diff}).
In section 4, we prove the second mail result of the paper
(Theorem \ref{improvement.formula}), that is the improvement of the formula
for the derivative of the functional $j'(0)$.
Finally, we close this paper with an appendix on some $\Gamma$-convergence
results that are needed in some parts of our arguments.

\section{Preliminaries}

Because of the nature of our problem, which deals with piecewise constant
exponents, we are unable to assume any regularity on the variable exponent $p$.
Therefore, since most of the known results for variable exponent Sobolev
spaces assume that the exponent is at least log-H\"older continuous,
we need to review the results that are needed here and prove the missing
parts in the case of piecewise constant exponents.

\subsection{Definitions and well-known results}

Given $\Omega\subset \mathbb R ^N$ a bounded open set, we consider the class of exponents
$\mathcal P(\Omega)$ given by
$$
\mathcal P(\Omega) := \{p\colon \Omega\to [1,\infty)\colon p
\text{ is measurable and bounded}\}.
$$
The variable exponent Lebesgue space $L^{p(x)}(\Omega)$ is defined by
$$
L^{p(x)}(\Omega):= \{f\in L^1_{\rm loc}(\Omega)\colon \rho_{p(x)}(f)<\infty\},
$$
where the modular $\rho_{p(x)}$ is
$$
\rho_{p(x)}(f) := \int_{\Omega} |f|^{p(x)}\,dx.
$$
This space is endowed with the Luxemburg norm
$$
\|f\|_{L^{p(x)}(\Omega)} = \|f\|_{p(x),\Omega} = \|f\|_{p(x)}
:= \inf\{\lambda>0\colon \rho_{p(x)}\big(\frac{f}{\lambda}\big)<1\}.
$$

The infimum and the supremum of the exponent $p$ play an important role in
the estimates as the next elementary proposition shows.
For further references, the following notation will be imposed
$$
1\le p_-:= \inf_{\Omega}p \le \sup_{\Omega} p =: p_+<\infty.
$$
The proof of the following proposition can be found in
\cite[Theorem 1.3, p.p. 427]{FanyZhao}.

\begin{proposition}\label{propdesigualdades}
Let $f\in L^{p(x)}(\Omega)$, then
$$
\min\big\{\|f\|_{p(x)}^{p_-}, \|f\|_{p(x)}^{p_+}\big\}
\leq \rho_{p(x)}(f)
\leq \max\big\{\|f\|_{p(x)}^{p_-}, \|f\|_{p(x)}^{p_{+}}\big\}.
$$
\end{proposition}

\begin{remark}\label{minmax} \rm
Proposition \ref{propdesigualdades}, is equivalent to
$$
\min\big\{\rho_{p(x)}(f)^{\frac{1}{p_-}}, \rho_{p(x)}(f)^{\frac{1}{p_{+}}} \big\}
\leq \|f\|_{p(x)}
\leq \max\big\{\rho_{p(x)}(f)^{\frac{1}{p_-}}, \rho_{p(x)}(f)^{\frac{1}{p_{+}}} \big\}.
$$
\end{remark}

We will use the following form of H\"older's inequality for variable exponents.
The proof, which is an easy consequence of Young's inequality, can be
found in \cite[Lemma 3.2.20]{Diening}.


\begin{proposition}[H\"older's inequality]\label{propholder}
Assume $p_->1$. Let $u\in L^{p(x)}(\Omega)$ and $v\in L^{p'(x)}(\Omega)$, then
$$
\int_{\Omega} |u v|\,dx\leq 2\|u\|_{p(x)} \|v\|_{p'(x)},
$$
where $p'(x)$ is, as usual, the conjugate exponent, i.e. $p'(x):= p(x)/(p(x)-1)$.
\end{proposition}

The variable exponent Sobolev space $W^{1,p(x)}$ is defined by
$$
W^{1,p(x)}(\Omega):=\big\{u\in W^{1,1}_{\rm loc}(\Omega)
\colon u\in L^{p(x)}(\Omega) \text{ and } \partial_i u\in L^{p(x)}(\Omega)\;
 i=1,\dots,N\big\},
$$
where $\partial_i u$ stands for the $i-$th partial weak derivative of $u$.

This space posses a natural modular given by
$$
\rho_{1,p(x)}(u) := \int_\Omega |u|^{p(x)} + |\nabla u|^{p(x)}\,dx,
$$
so $u\in W^{1,p(x)}(\Omega)$ if and only if $\rho_{1,p(x)}(u)<\infty$.

The corresponding Luxemburg norm associated to this modular is
$$
\|u\|_{W^{1,p(x)}(\Omega)} = \|u\|_{1,p(x),\Omega}
= \|u\|_{1,p(x)}
:= \inf\big\{\lambda>0\colon \rho_{1,p(x)}\big(\frac{u}{\lambda}\big)<1\big\}.
$$
Observe that this norm turns out to be equivalent to
$\|u\|:= \|u\|_{p(x)} + \|\nabla u\|_{p(x)}$.
Now we state and prove a simple proposition that characterizes the Sobolev
space when the variable exponent is piecewise constant.

\begin{proposition}\label{prop.24}
Let $\Omega\subset \mathbb R ^N$ be an open of finite measure and let
$D_1, D_2\subset \Omega$ be a partition verifying \eqref{prop.particion}.
Let $1\le p_1,p_2<\infty$ and let $p\in \mathcal P(\Omega)$ be such that
$p=p_1\chi_{D_1} + p_2\chi_{D_2}$.
Then $u \in W^{1,p(x)}(\Omega)$ if and only if $u \in W^{1,p_1}(D_1)$,
$u \in W^{1,p_2}(D_2)$ and $u \in W^{1,\min \{p_1,p_2\}}(\Omega)$.
\end{proposition}

\begin{proof}
Observe that \eqref{prop.particion} implies that
 $|\Omega \setminus (D_1\cup D_2)|=0$. Then
\begin{gather*}
\int_{\Omega} |u|^{p(x)}\,dx = \int_{D_1} |u|^{p_1}\,dx
+ \int_{D_2} |u|^{p_2}\,dx, \\
\int_{\Omega} |\nabla u|^{p(x)}\,dx
= \int_{D_1} |\nabla u|^{p_1}\,dx + \int_{D_2} |\nabla u|^{p_2}\,dx.
\end{gather*}
Moreover, assuming that $p_1<p_2$ by H\"older's inequality, we have
\begin{align*}
\int_{\Omega} |\nabla u|^{p_1}\,dx
&=  \int_{D_1} |\nabla u|^{p_1}\,dx + \int_{D_2} |\nabla u|^{p_1}\,dx\\
&\le \int_{D_1} |\nabla u|^{p_1}\,dx
+ |D_2|^{\frac{p_2-p_1}{p_2}}\Big(\int_{D_2} |\nabla u|^{p_2}\,dx\Big)^{p_1/p_2}
<\infty.
\end{align*}
Analogously, $u\in L^{p_1}(\Omega)$.

For the converse, we just observe that since
$u\in W^{1, \min\{p_1,p_2\}}(\Omega)$, then $\nabla u$ is defined in the whole
of $\Omega$. Then is easy to see that $\nabla u\in L^{p(x)}(\Omega)$ by
the same arguments as before.
\end{proof}

\section{Differentiability}

Let $V$ be a Lipschitz vector field with support in $\Omega$ and let
$\{\Phi_t\}_{t\in\mathbb R }$ its associated flux given by  \eqref{flux}.
Let us begin with the following observation.

\begin{remark}\label{asintoticas} \rm
By Taylor expansion, we have
$$
\Phi_t(x)=x+V(x)t+o(t)
$$
and so we have the  asymptotic formulas
\begin{gather*}
D\Phi_t(x)=Id+tDV(x)+o(t)=Id+O(t), \\
J\Phi_t(x)=1+t \operatorname{div} V (x)+o(t)=1+O(t),
\end{gather*}
for all $x \in \mathbb R ^N$, where $J\Phi_t$ is the Jacobian of $\Phi_t$.
\end{remark}

The following proposition, though elementary, will be useful in the sequel
and shows that any diffeomorphism $\Phi\colon \mathbb R ^N\to\mathbb R ^N$, induces a
bounded linear isomorphism between Sobolev spaces.

\begin{proposition}\label{prop.iso}
Let $\Phi \colon \Omega_1 \to \Omega_2$ be a diffeomorphism and
$p \in \mathcal P(\Omega_1)$ be a bounded exponent.

Then, $\Phi$ induces a bounded linear isomorphism
$$
\mathcal F \colon W^{1,p}(\Omega_1) \to W^{1,q}(\Omega_2),
$$
where $q \colon \Omega_2 \to [1, +\infty)$ is given by $q(x):=p(\Phi^{-1}(x))$,
by the expression
$$
\mathcal F(u):=u\circ \Phi^{-1}.
$$
\end{proposition}

\begin{proof}	
We first observe that $\mathcal F$ is clearly a linear isomorphism with inverse given by
$$
\mathcal F^{-1} \colon W^{1,q}(\Omega_2) \to W^{1,p}(\Omega_1), \quad
 \mathcal F^{-1}(v):=v\circ \Phi.
$$
Let us now show that it is also bounded.
Let us consider $\lambda>0$ and, for simplicity, let us denote $v=\mathcal F(u)$.
Then, by changing variable $y=\Phi^{-1}(x)$,
\begin{align*}
\int_{\Omega_2}\big|\frac{v(x)}{\lambda}\big|^{q(x)}\,dx
&=\int_{\Omega_2}\big|\frac{u(\Phi^{-1}(x))}{\lambda}\big|^{p(\Phi^{-1}(x))}\,dx\\
&=\int_{\Omega_1}\big|\frac{u(y)}{\lambda}\big|^{p(y)}J\Phi(y)\,dy\\
&\leq \|J \Phi\|_{\infty}\int_{\Omega_1}\big|\frac{u(y)}{\lambda}\big|^{p(y)}dy
\end{align*}
Let us observe that, if $C:=\|J \Phi\|_{\infty}\leq 1$, clearly we have
\begin{align*}
\|u\|_{p,\Omega_1}
&=\inf \big\{\lambda>0 \colon \int_{\Omega_1}
\big|\frac{u(y)}{\lambda}\big|^{p(y)}dy\leq 1\big\}
\geq \inf \big\{\lambda>0 \colon \int_{\Omega_2}\big|\frac{v(y)}{\lambda}
 \big|^{q(y)}dy\leq 1\big\} \\
&=\|v\|_{q,\Omega_2}
\end{align*}
Let us now assume that $C>1$. Then, since
\begin{align*}
\big\{\lambda>0 \colon \int_{\Omega_1}\big|\frac{C^{\frac{1}{p_-}}u(y)}{\lambda}
\big|^{p(y)}\,dy \leq 1\big\}
 & \subset \big\{\lambda>0 \colon \int_{\Omega_1}\big|\frac{u(y)}{\lambda}
\big|^{p(y)}\,dy \leq \frac{1}{C}\big\}\\
&\subset \big\{\lambda>0 \colon \int_{\Omega_2}\big|\frac{v(x)}{\lambda}
\big|^{q(x)}\,dx\leq 1\big\},
\end{align*}
taking the infimum, we conclude that
\begin{align*}
C^{\frac{1}{p_-}}\|u\|_{p,\Omega_1}
&=\|C^{\frac{1}{p_-}}u\|_{p,\Omega_1}
 \geq\inf\big\{\lambda>0 \colon \int_{\Omega_1}\big|\frac{u(y)}{\lambda}
 \big|^{p(y)}\,dy \leq \frac{1}{C}\big\} \\
&\geq \|v\|_{q,\Omega_2}=\|\mathcal F(u)\|_{q, \Omega_2}.
\end{align*}
Analogously,
\begin{align*}
\int_{\Omega_2}\big|\frac{\nabla v(x)}{\lambda}\big|^{q(x)}\,dx
&= \int_{\Omega_2}\big|\frac{\nabla (u\circ \Phi^{-1})(x)}{\lambda}\big|^{q(x)}\,dx\\
&= \int_{\Omega_2}\big|\frac{\nabla u(\Phi^{-1}(x)) D\Phi^{-1}(x)}{\lambda}
 \big|^{p(\Phi^{-1}(x))}\,dx\\
&= \int_{\Omega_1}\big|\frac{\nabla u(y) D\Phi^{-1}(\Phi(y))}{\lambda}
 \big|^{p(y)}J \Phi(y)\,dy\\
&\leq \max\{1, \|D \Phi^{-1}\|_{\infty}\}^{p_+}
 \|J \Phi\|_{\infty}\int_{\Omega_1}\big|\frac{\nabla u(y)}{\lambda}\big|^{p(y)}\,dy.
\end{align*}
Therefore, $\|\nabla \mathcal F(u)\|_{q, \Omega_2}\leq C \|\nabla u\|_{p, \Omega_1}$,
which completes the proof.
\end{proof}

\begin{remark} \rm
In the previous proof, given $A \colon \Omega \to \mathbb{R}^{N\times N}$,
we considered the norm $\|A\|_{\infty}:=\sup_{x\in \Omega}\|A(x)\|$ and,
given $B \in \mathbb{R}^{N\times N}$, we considered the norm
$\|B\|:=\sup_{\xi\neq 0} |B\xi|/|\xi|$.
\end{remark}

Observe that, since $\operatorname{suppt}(V)\subset\subset \Omega$, it follows that
$\Phi_t(\Omega)=\Omega$ for every $t\in\mathbb R $ and that if
 $p = p_1\chi_{D_1} + p_2\chi_{D_2}$ then
$p_t:=p\circ \Phi_{-t} = p_1 \chi_{D_1^t} + p_2 \chi_{D_2^t}$, where
$D_i^t = \Phi_t(D_i)$, $i=1,2$.

Therefore, in view of Proposition \ref{prop.iso}, we have that
$$
\mathcal F_t\colon W^{1,p}(\Omega)\to W^{1,p_t}(\Omega), \quad u\mapsto u\circ \Phi_{-t}
$$
is a bounded linear isomorphism.

Let us consider the space $X_t:=W^{1,p_t}(\Omega)\cap L^2(\Omega)$ equipped
with the norm
$$
\|\cdot\|_{X_t}:=\|\cdot\|_{W^{1,p_t}(\Omega)}+\|\cdot\|_{L^2(\Omega)}
$$
and the space $X:=W^{1,p}(\Omega)\cap L^2(\Omega)$ equipped with the norm
$$
\|\cdot\|_X:=\|\cdot\|_{W^{1,p}(\Omega)}+\|\cdot\|_{L^2(\Omega)}.
$$
It is clear that $\mathcal F_t\colon X\to X_t$ is still a bounded linear isomorphism.

Given $f \in L^2(\Omega)$, we define the quantity
$$
\tilde s(t):= \inf_{v \in X_t} \int_{\Omega}\frac{|\nabla v|^{p_t}}{p_t}\,dx
+\int_{\Omega}\frac{|v-f|^2}{2}\,dx
$$
which is clearly equivalent to
\begin{equation}\label{st}
s(t):= \inf_{v \in X_t} \int_{\Omega}\frac{|\nabla v|^{p_t}}{p_t}\,dx
+\int_{\Omega} \frac{|v|^2}{2} \,dx- \int_{\Omega}v f\,dx.
\end{equation}
In fact, $\tilde s(t) = s(t) + \|f\|_2^2$.

Observe that, since $\mathcal F_t$ is an isomorphism, one actually has
$$
s(t) = \inf_{u\in X} \int_{\Omega} \frac{|\nabla (u\circ \Phi_{-t})|^{p_t}}{p_t}\,dx
 + \int_{\Omega} \frac{|u\circ \Phi_{-t}|^2}{2}\,dx
- \int_{\Omega} (u\circ \Phi_{-t}) f\,dx.
$$
So, in view of our previous discussions, our primary goal is to find
an expression for $\frac{ds}{dt}(0)$.

\begin{remark}\label{rem.equivalencia} \rm
Let us observe that, by changing variables $y=\Phi_{-t}(x)$,
$$
s(t)= \inf_{u \in X} \int_{\Omega}\frac{|\nabla u D\Phi_{-t}\circ
\Phi_{t} |^p}{p} J\Phi_t\,dy +\int_{\Omega} \frac{|u|^2}{2} J\Phi_t\,dy
-\int_{\Omega} u f\circ \Phi_t J\Phi_t\,dy.
$$
\end{remark}
Let
\begin{gather*}
\mathcal J_t u := \int_{\Omega}\frac{|\nabla u D\Phi_{-t}\circ \Phi_{t} |^p}{p}
J\Phi_t\,dy +\int_{\Omega} \frac{|u|^2}{2} J\Phi_t\,dy
- \int_{\Omega} u f\circ \Phi_t J\Phi_t\,dy, \\
\mathcal J u := \int_{\Omega}\frac{|\nabla u|^p}{p}\,dy
+ \int_{\Omega} \frac{|u|^2}{2} \,dy - \int_{\Omega}u f\,dy.
\end{gather*}

\begin{lemma}\label{equicoercividad}
There exists $\delta>0$ such that the functionals $\{\mathcal J_t\}_{|t|<\delta}$
are uniformly coercive with respect to the weak topology of $X$.
That is, for any $\lambda\in\mathbb R $, there exists a weakly compact
set $K\subset X$ such that
$$
\{\mathcal J_t \le \lambda\} \subset K, \quad \text{for every } |t|<\delta.
$$
\end{lemma}

\begin{proof}
Take $\delta>0$ such that $1/2 \leq J\Phi_t \leq 2$. Therefore,
\begin{equation}\label{cotaJ}
\mathcal J_t u \geq \frac{1}{2}\int_{\Omega}\frac{|\nabla u D\Phi_{-t}\circ \Phi_{t} |^p}{p}
 \,dy +\frac{1}{2}\int_{\Omega} \frac{|u|^2}{2} \,dy
 - 2 \int_{\Omega} |f| |u|\,dy.
\end{equation}
By Young's inequality with $\epsilon=1/8$,
\begin{equation}\label{young4}
2\int_{\Omega} |f| |u|\,dy
\leq \frac{1}{8}\int_{\Omega} |u|^2\,dy+8\int_\Omega |f|^2\,dy.
\end{equation}

As $D\Phi_{-t}\rightrightarrows Id$ uniformly on $\Omega$, it follows that
$\|D\Phi_t\|_\infty$ is bounded away from zero and infinity for every $|t|< \delta$,
so
\begin{equation}\label{D-Id}
\int_{\Omega}\frac{|\nabla u D\Phi_{-t}\circ \Phi_{t} |^p}{p} \,dy
\geq c \int_{\Omega}|\nabla u|^p \,dy.
\end{equation}
So, combining \eqref{cotaJ}, \eqref{young4} and \eqref{D-Id}, we obtain
$$
\mathcal J_t u \geq c\int_\Omega |\nabla u|^p\,dy + \frac18 \int_\Omega |u|^2\,dy
- 8\|f\|_2^2.
$$
By Proposition \ref{minmax} we easily conclude that there exists a radius
$R=R(\lambda)$ such that $\{\mathcal J_t\le \lambda\}\subset B_X(0, R)$.
Therefore, if we denote $K:=\{\|u\|_{X}<R\}$, satisfies our requirements.
This completes the proof.
\end{proof}


\begin{lemma}\label{lemmaextremales}
There exists a unique extremal for $s(t)$ and $s(0)$.
\end{lemma}

\begin{proof}
The proof is an immediate consequence of the fact that both $\mathcal J_t$ and $\mathcal J$
are strictly convex and sequentially weakly lower semicontinuous
on $W^{1,p}(\Omega)$.
\end{proof}

Our first result shows that $s(t)$ is continuous with respect to $t$ at $t=0$.


\begin{theorem}\label{teo.continuidad}
With the previous notation,
\begin{equation}\label{st.to.s}
\lim_{t \to 0^+}s(t)=s(0).
\end{equation}
Moreover, if $u_t$ and $u$ are the extremals associated to $s(t)$ and $s(0)$
respectively, then $u_t\rightharpoonup u$ weakly in $W^{1,p}(\Omega)$.
Finally, if $p^*:=\frac{pN}{N-p}>2$ then $u_t\to u$ strongly in $W^{1,p}(\Omega)$.
\end{theorem}

\begin{remark} \rm
The hypothesis $p^*>2$ is needed in order to secure the compact embedding
$W^{1,p}(\Omega)\subset L^2(\Omega)$ for any dimension $N$.
For the case $N=2$, one has $p^*>2$ for any $p>1$ so no extra hypothesis is needed.
\end{remark}

\begin{proof}[Proof of Theorem \ref{teo.continuidad}]
Since, by Lemma \ref{equicoercividad}, we know that the functionals $\mathcal J_t$
are uniformly coercive, the proof of \eqref{st.to.s} will follow from
Remark \ref{obs gamma convergencia} if we show that $\mathcal J_t \rightrightarrows \mathcal J$
uniformly on bounded sets of $X$. Observe that since the minimizers are unique,
we will then have that the whole sequence of minimizers is weakly convergent.

Let us consider now $B \subset X$ a bounded subset and $u \in B$.
By Remark \ref{asintoticas},
\begin{align*}
\mathcal J_t u&=\int_{\Omega}\frac{|\nabla u (Id+O(t)) |^p}{p} (1+O(t))\,dy
  +\int_{\Omega} \frac{|u|^2}{2} (1+O(t))\,dy \\
 &\quad - \int_{\Omega} u (f\circ \Phi_{t})(1+O(t))\,dy\\
&=(1+O(t))\Big\{\int_{\Omega}\frac{|\nabla u(Id +O(t))|^p}{p} \,dy
 +\int_{\Omega} \frac{|u|^2}{2}\,dy - \int_{\Omega} u (f\circ \Phi_t)\,dy \Big\}.
\end{align*}
Again by Remark \ref{asintoticas}, and by Taylor expansion formula, we obtain
$$
\int_{\Omega}\frac{|\nabla u(Id +O(t))|^p}{p} \,dy
=\int_{\Omega}\frac{|\nabla u|^p}{p}\,dy + O(t),
$$
uniformly in $B$.

Assume for a moment that $f$ is a continuous function with compact support.
Then, since $\Phi_{t}\to id$ uniformly as $t \to 0$, we have that
$f\circ\Phi_t\to f$ uniformly as $t \to 0$ and therefore,
$$
\|f\circ\Phi_t-f\|^2_{2}=\int_{\Omega}|f\circ\Phi_t-f|^2\,dx
\leq \|f\circ\Phi_t-f\|^2_{\infty}|\Omega|\to 0,\quad (t\to 0).
$$
So we have that $\|f\circ\Phi_t-f\|_{2}\to 0$, $(t\to 0)$.
Now, by a standard density argument, it is easy to see that the same result
holds for any $f \in L^2(\Omega)$.

Then, by H\"older inequality and since $u \in B$, there is a constant $C$,
independent of $u$, such that
$$
\big|\int_{\Omega}u(f\circ \Phi_t-f)\big|\leq C\|f\circ \Phi_t - f\|_2\to 0
$$
as $t\to 0^+$.

Assume now that $p^*>2$. It remains to see the strong convergence of $u_t$
to $u$ in $W^{1,p}(\Omega)$.
Let us observe that  to see the strong convergence it is enough to
show the convergence of the modulars (see \cite{Diening}).
Let us now recall that
\begin{align*}
\int_{\Omega}\frac{|\nabla u_t|^p}{p}\,dy
+\int_{\Omega} \frac{|u_t|^2}{2}\,dy
&= s(t) + \int_{\Omega}\frac{|\nabla u_t|^p}{p}\,dy
 -\int_{\Omega}\frac{|\nabla u_t D\Phi_{-t}\circ \Phi_{t} |^p}{p} J\Phi_t\,dy\\
&\quad +\int_{\Omega} \frac{|u_t|^2}{2}(1-J\Phi_t)\,dy
 +\int_{\Omega} u_t (f\circ \Phi_t) J\Phi_t\,dy.
\end{align*}
By Remark \ref{asintoticas},
\[
\int_{\Omega}\frac{|\nabla u_t D\Phi_{-t}\circ \Phi_{t} |^p}{p} J\Phi_t\,dy
=\int_{\Omega} \frac{|\nabla u_t - t\nabla u_t DV
 + o(t)|^p}{p} (1 + t\operatorname{div} V + o(t))\,dy.
\]
Using the Taylor expansion
$$
|\nabla u_t - t \nabla u_t DV + o(t)|^p
= |\nabla u_t|^p- p t |\nabla u_t|^{p-2} \nabla u_t\cdot\nabla u_t DV + o(t),
$$
we find that
\begin{align*}
&\int_{\Omega}\frac{|\nabla u_t D\Phi_{-t}\circ \Phi_{t}|^p}{p} J\Phi_t\,dy \\
&= \int_{\Omega}\frac{|\nabla u_t|^p+t(|\nabla u_t|^p
\operatorname{div} V-p|\nabla u_t|^{p-2}\nabla u_t\cdot\nabla u_t DV)}{p}\,dy+o(t).
\end{align*}
So we have
\begin{align*}
&\int_{\Omega}\frac{|\nabla u_t|^p}{p}\,dy
 -\int_{\Omega}\frac{|\nabla u_t D\Phi_{-t}\circ \Phi_{t} |^p}{p} J\Phi_t\,dy \\
&=-\int_{\Omega}\frac{t(|\nabla u_t|^p\operatorname{div} V-p|\nabla u_t|^{p-2}
 \nabla u_t\cdot\nabla u_t DV)}{p}\,dy+o(t)
\end{align*}

Now, for our fourth term, we only need to observe that $\frac{|u_t|^2}{2}$
is bounded and $1-J\Phi_t \to 0$ uniformly.
Then, since $s(t)\to s(0)$ and
$$
\int_{\Omega} u_t (f\circ \Phi_t) J\Phi_t\,dy\to \int_{\Omega} u f,
$$
we can conclude that
$$
\int_{\Omega}\frac{|\nabla u_t|^p}{p}\,dy+\int_{\Omega} \frac{|u_t|^2}{2}\,dy
 \to \int_{\Omega}\frac{|\nabla u|^p}{p}\,dy+\int_{\Omega} \frac{|u|^2}{2}\,dy,
$$
which completes the proof.
\end{proof}

Now we prove the main result of the section, namely the differentiability
of the cost functional $s(t)$. For this result we will need the function $f$
to be of class $C^1$.

\begin{theorem}\label{teo.diff}
$s(t)$ is differentiable at $t=0$ and
$$
\frac{ds}{dt}(0)=R(u) - \int_\Omega uf \operatorname{div} V\,dy
- \int_\Omega u \nabla f\cdot V\,dy,
$$
where
$$
R(u):=\int_{\Omega}\frac{|\nabla u|^p}{p}\operatorname{div} V-|\nabla u|^{p-2}
\nabla u\cdot\nabla uDV+ \operatorname{div} V\frac{|u|^2}{2}\,dy
$$
and $u$ is the extremal of $s(0)$.
\end{theorem}

\begin{proof}
By Lemma \ref{lemmaextremales}, we can consider $u$ the extremal of $s(0)$.
Then, by Remark \ref{rem.equivalencia},
\begin{align*}
s(t)&=\inf_{X} \mathcal J_{t} \leq \mathcal J_{t}(u) \\
&= \int_{\Omega}\frac{|\nabla u D\Phi_{-t}\circ \Phi_{t} |^p}{p} J\Phi_t\,dy
 +\int_{\Omega} \frac{|u|^2}{2} J\Phi_t\,dy
  - \int_{\Omega} u f\circ \Phi_t J\Phi_t\,dy.
\end{align*}
Now, by Remark \ref{asintoticas}, as in the proof of Theorem \ref{teo.continuidad}
we find that
\begin{align*}
&\int_{\Omega}\frac{|\nabla u D\Phi_{-t}\circ \Phi_{t}|^p}{p} J\Phi_t\,dy \\
&= \int_{\Omega}\frac{|\nabla u|^p+t(|\nabla u|^p\operatorname{div} V-p|\nabla u|^{p-2}
 \nabla u\cdot\nabla uDV)}{p}\,dy+o(t).
\end{align*}
On the other hand, again by Remark \ref{asintoticas},
\begin{align*}
\int_{\Omega} \frac{|u|^2}{2} J\Phi_t\,dy
&= \int_{\Omega} \frac{|u|^2}{2} (1+t\operatorname{div} V+o(t))\,dy\\
&= \int_{\Omega} \frac{|u|^2}{2}\,dy+t\int_{\Omega} \operatorname{div} V\frac{|u|^2}{2}\,dy+o(t).
\end{align*}
Therefore, setting
$$
R(u):=\int_{\Omega}\frac{|\nabla u|^p}{p}\operatorname{div} V-|\nabla u|^{p-2}
\nabla u\cdot\nabla uDV+ \operatorname{div} V\frac{|u|^2}{2}\,dy,
$$
we conclude that
$$
s(t)\leq \int_{\Omega}\frac{|\nabla u|^p}{p}\,dy
 +\int_{\Omega} \frac{|u|^2}{2}\,dy+tR(u)+o(t)
 -\int_{\Omega} u (f\circ \Phi_t) (1 + t\operatorname{div} V + o(t))\,dy.
$$
Recall that
$$
s(0) = \int_{\Omega}\frac{|\nabla u|^p}{p}\,dy+\int_{\Omega} \frac{u^2}{2}\,dy
- \int_\Omega u f\,dy.
$$
Therefore,
$$
\frac{s(t)-s(0)}{t}\leq R(u)+\frac{o(t)}{t}
 -\int_{\Omega}u (f\circ \Phi_t) \operatorname{div} V\,dy
 - \int_\Omega u \frac{(f\circ \Phi_t) - f}{t}\,dy.
$$
Taking the limit as $t\to 0^{+}$, we obtain
$$
\limsup_{t\to 0+} \frac{s(t)-s(0)}{t} \leq R(u)
- \int_\Omega uf \operatorname{div} V\,dy - \int_\Omega u \nabla f\cdot V\,dy,
$$
where we have used  that $\Phi_0=id$ and $\dot \Phi_t = V\circ \Phi_t$.

Let us consider now $\{t_n\}_{n \in \mathbb{N}}$ such that $t_{n}\to 0^+$ and
$$
\liminf_{t \to 0^+}\frac{s(t)-s(0)}{t}
=\lim_{n \to \infty}\frac{s(t_n)-s(0)}{t_n}.
$$
Let $u_n:=u_{t_n}\in X_{t_n}$ be the extremal associated to $s(t_n)$.
By Remark \ref{rem.equivalencia},
$$
s(t_n)= \int_{\Omega}\frac{|\nabla u_n D\Phi_{-t_n}\circ \Phi_{t_n} |^p}{p}
J\Phi_{t_n}\,dy +\int_{\Omega} \frac{|u_n|^2}{2} J\Phi_{t_n}\,dy
-\int_{\Omega} u_{t_n} f\circ \Phi_{t_n} J\Phi_{t_n}\,dy
$$
Arguing as in the previous case, we have that
\begin{align*}
\frac{s(t_n)-s(0)}{t_n}
\geq & \int_{\Omega}\frac{|\nabla u_n|^p}{p}\operatorname{div} V-|\nabla u_n|^{p-2}
 \nabla u_n\cdot\nabla u_{n}DV+ \operatorname{div} V\frac{|u_n|^2}{2}\,dy\\
&+\frac{o(t_n)}{t_n}-\int_{\Omega}u_{n} (f\circ \Phi_{t_n}) \operatorname{div} V\,dy
 - \int_\Omega u_n \frac{(f\circ \Phi_{t_n}) - f}{t_n}\,dy\\
=& R(u_n) + \frac{o(t_n)}{t_n}-\int_{\Omega}u_{n} (f\circ \Phi_{t_n}) \operatorname{div} V\,dy
 - \int_\Omega u_n \frac{(f\circ \Phi_{t_n}) - f}{t_n}\,dy.
\end{align*}
Since $R(u_n) \to R(u)$ when $n \to \infty$ (just observe that $R$ is continuous
with respect to the strong topology and $u_n\to u$ in $W^{1,p}(\Omega)$ by
Theorem \ref{teo.continuidad}), we have
$$
\liminf_{t \to 0^+}\frac{s(t)-s(0)}{t}\geq R(u)
- \int_\Omega uf \operatorname{div} V\,dy - \int_\Omega u \nabla f\cdot V\,dy.
$$
So we  conclude that $s(t)$ is differentiable at $t=0$ and
$$
\frac{ds}{dt}(0)=R(u) - \int_\Omega uf \operatorname{div} V\,dy
- \int_\Omega u \nabla f\cdot V\,dy,
$$
where $u\in X$ is the extremal of $s(0)$. This completes the proof.
\end{proof}

\section{Improvement of the formula}

Now we try to find a more explicit formula for $s'(0)$.
In the following study, we will need the solution $u$ to the equation
\begin{equation}
\begin{gathered}
-\Delta_{p(x)}u+u=f \quad \text{in } \Omega,\\
u=0 \quad \text{on } \partial\Omega,
\end{gathered}
\end{equation}
to be $C_{\rm loc}^2(D_1)\cap C^2_{\rm loc}(D_2)$ in order for our
computations to work. However, this is not true since the optimal regularity
is known to be $C_{\rm loc}^{1,\alpha}(D_1)\cap C^{1,\alpha}_{\rm loc}(D_2)$.
See \cite{Tolksdorf}.
To overcome such difficulty, we will proceed as follows.

\subsection{Domain regularization}

Let us first define $D_i(t):=\Phi_t(D_i)$.
Now given a fixed $\delta>0$, we define the  sets
$$
D^{\delta}_i:=\{x\in D_i \colon \operatorname{dist}(x,D_j)>\delta\},\ i\neq j
$$
and consider $D^{\delta}_i(t):=\Phi_t(D^{\delta}_i)$. Now consider the sets
$$
\Gamma^{\delta}_i(t):= \partial D^{\delta}_i(t) \cap \Omega.
$$
Let us observe that, in each $D^{\delta}_i$, the exponent $p(x)=p_i$
is constant so we can apply the classic regularity results.
See for instance \cite{Tolksdorf}.

Now we define the sets $A_\delta:= \Omega \setminus (D_1^\delta \cup D_2^\delta)$
and observe that
$$
\partial A_\delta \cap \Omega=\Gamma^{\delta}_1\cup \Gamma^{\delta}_2,
\quad \Omega=D^{\delta}_1\cup D^{\delta}_2\cup A^{\delta}.
$$
See Figure \ref{fig1}.

\begin{figure}
\begin{center}
  \includegraphics[width=0.6\textwidth]{fig1}
\end{center}
\caption{Partition of $\Omega$} \label{fig1}
\end{figure}

\subsection{Operator regularization}
Now, for $\epsilon \geq 0$, we consider the regularized problems
\begin{equation}\label{eqregul}
\begin{gathered}
-\operatorname{div}((|\nabla v|^2+\epsilon^2)^{\frac{p_i-2}{2}}\nabla v)+v=f^{\epsilon} \quad
 \text{in } D^{\delta}_i(t),\\
v=0 \quad \text{on } \partial\Omega \cap (\overline{D^{\delta}_1(t)}
 \cup \overline{D^{\delta}_2(t)}),\\
v=u(0)\circ \Phi^{-1}_t \quad \text{on } \Gamma^{\delta}_i(t),
\end{gathered}
\end{equation}
with $f^{\epsilon}\in C^{\infty}$ such that $f^{\epsilon}\to f$ in $L^{p'}$.

\begin{remark} \rm
Applying classical estimates (see for instance \cite{GT} it is possible to
see that the solution of \eqref{eqregul} is
$C^{2,\alpha}_{\rm loc}(D^{\delta}_i)\cap C^{1}(\overline{D^{\delta}_i})$ if
$\epsilon > 0$, since $u(0)$ is $C^1(\overline{D^{\delta}_i})$ and $\Phi_t$
is the identity map in a neighborhood of $\partial\Omega$.
See also \cite{Tolksdorf} for regularity estimates in Sobolev spaces.
\end{remark}

Let us define the following sets
\begin{gather*}
X^{\delta}_i:=\{v \in W^{1,p_i}(D^{\delta}_i) \text{ such that } v=0 \text{ in }
\partial\Omega \cap \overline{D^{\delta}_i} \text{ and }
v=u  \text{ in } \Gamma^{\delta}_i\}, \\
\begin{aligned}
X^{\delta}_i(t):=\big\{&v \in W^{1,p_i}(D^{\delta}_i(t)) \text{ such that } v=0
 \text{ in } \partial\Omega \cap \overline{D^{\delta}_i(t)} \text{ and }\\
& v=u(0)\circ \Phi^{-1}_t  \text{ in } \Gamma^{\delta}_i(t)\}.
\end{aligned}
\end{gather*}
Let us also consider the functionals
$\tilde{\mathcal J}^{\epsilon,\delta}_{t,i}\colon X^{\delta}_i(t) \to \mathbb R $ defined by
$$
\tilde{\mathcal J}^{\epsilon,\delta}_{t,i}(v):=\int_{D^{\delta}_i(t)}\frac{\big(|\nabla v|^2
 + \epsilon^2\big)^{\frac{p_i}{2}}}{p_i}\,dy +\int_{D^{\delta}_i(t)} \frac{|v|^2}{2}\,dy -\int_{D^{\delta}_i(t)} v f^{\epsilon}\,dy.
$$

\begin{remark} \rm
$X^{\delta}_i(t)$ is strongly closed and convex, therefore it is weakly closed.

The solutions of \eqref{eqregul} are the minimuma of the functionals
$\tilde{\mathcal J}^{\epsilon,\delta}_{t,i}$ in $X^{\delta}_i(t)$.

Since the functional $\tilde{\mathcal J}^{\epsilon,\delta}_{t,i}$ is continuous
for the strong topology, strictly convex and coercive, it has a unique
minimum in $X^{\delta}_i(t)$ and, therefore, \eqref{eqregul} has a
unique weak solution.
We will denote $\tilde{u}^{\epsilon,\delta}_i(t)$ as the function where
the minimum is attained.
\end{remark}

\begin{remark} \rm
Observe that $\psi_t \colon X^{\delta}_i \to X^{\delta}_{i}(t)$ defined by
 $v \mapsto v \circ \Phi_{t}^{-1}$ is a biyection between $X^{\delta}_i$
and $X^{\delta}_{i}(t)$ and the following equality holds
$$
\mathcal J^{\epsilon,\delta}_{t,i}=\tilde{\mathcal J}^{\epsilon,\delta}_{t,i}\circ \psi_{t}^{-1}.
$$
\end{remark}

By a change of variables as in the previous section we obtain the functional
$\mathcal J^{\epsilon,\delta}_{t,i}\colon X^{\delta}_i \to \mathbb R $ given by
\begin{align*}
\mathcal J^{\epsilon,\delta}_{t,i}(v)
:=&\int_{D^{\delta}_i}\frac{\big(|\nabla v D\Phi_{-t}\circ \Phi_{t} |^2
 + \epsilon^2\big)^{\frac{p_i}{2}}}{p_i} J\Phi_t\,dy
  +\int_{D^{\delta}_i} \frac{|v|^2}{2} J\Phi_t\,dy \\
&-\int_{D^{\delta}_i} v f^{\epsilon}\circ \Phi_t J\Phi_t\,dy.
\end{align*}
and define
$$
s^{\epsilon,\delta}_i(t)= \inf_{v \in X^{\delta}_i} \mathcal J^{\epsilon,\delta}_{t,i}(v)
= \inf_{v \in X^{\delta}_i(t)} \tilde{\mathcal J}^{\epsilon,\delta}_{t,i}(v).
$$
We will denote $u^{\epsilon,\delta}_i(t)\in X_i^\delta$ as the function where
the minimum of $\mathcal J^{\epsilon,\delta}_{t,i}$ is attained.
Observe that $u^{\epsilon,\delta}_i(t)(x)=\tilde{u}^{\epsilon,\delta}_i(t)(\Phi_t(x))$.

To make the notation lighter, we will focus on the needed parameter in each step.
First, $u_i$, then $u^{\epsilon}$ and finally, $u^{\delta}$.
Let us now define  $s^{\epsilon,\delta}:=s^{\epsilon,\delta}_1+s^{\epsilon,\delta}_2$.

\begin{proposition}\label{convergenciasolu}
If $2<p_i^*$, then $u^{\epsilon,\delta}_i(0)$ converges to
$u^{0,\delta}_i(0)(=u^{\delta}_i)$ strongly in $W^{1,p_i}(D^{\delta}_i)$
and $s^{\epsilon,\delta}_i(0)$ converges to $s^{0,\delta}_i(0)(=s^{\delta}_i(0))$
when $\epsilon \to 0$.
\end{proposition}

\begin{proof}
Let us begin by observing that
$$
\mathcal J^{\epsilon,\delta}_{0,i}(v):=\int_{D^{\delta}_i}\frac{\big(|\nabla v |^2
 + \epsilon^2\big)^{\frac{p_i}{2}}}{p_i}\,dy
+\int_{D^{\delta}_i} \frac{|v|^2}{2} \,dy -\int_{D^{\delta}_i} v f^{\epsilon}\,dy.
$$
Now let us denote
$$
\mathcal J^{\delta}_i(v):=\int_{D^{\delta}_i}\frac{|\nabla v |^2}{p_i}\,dy
+\int_{D^{\delta}_i} \frac{|v|^2}{2} \,dy -\int_{D^{\delta}_i} v f\,dy.
$$
Observe that $\mathcal J^{\epsilon,\delta}_{0,i}$, $\mathcal J^{\delta}_i(v):X^{\delta}_i \to \mathbb R $.
By Theorem \ref{convergencia sobre acotados}, it is enough to prove that
$\mathcal J^{\epsilon,\delta}_{0,i}$ $\Gamma$- converges to $\mathcal J^{\delta}_i$ in
$W^{1,p_i}(D^{\delta}_i)$ for the weak topology.

First, let $v^{\epsilon}\rightharpoonup v$ weakly in $W^{1,p_i}(D^{\delta}_i)$.
 Let us observe that $v \in X^{\delta}_i$ since $X^{\delta}_i$ is weakly closed.

Observe that the first and second terms in $\mathcal J^{\delta}_{i}$ are convex and
strongly continuous, therefore weakly lower semicontinuous. And the third term
is linear and continuous, therefore weakly continuous.
Therefore,
$$
\mathcal J^{\delta}_{i}(v)\leq \lim \mathcal J^{\delta}_{i}(v^{\epsilon})
\leq \liminf\mathcal J^{\epsilon,\delta}_{0,i}(v^{\epsilon})
+\int_{D^{\delta}_i} v^{\epsilon}(f-f^{\epsilon}).
$$
Applying H\"older's inequality for the last term above, we have
$$
\int_{D^{\delta}_i} v^{\epsilon}(f-f^{\epsilon})
\leq ||v^{\epsilon}||_{p_i}||f-f^{\epsilon}||_{p'_i}.
$$
Since $||v^{\epsilon}||_{p_i}$ is bounded (because of the weak convergence)
and $f^{\epsilon} \to f$ in $L^{p'_i}$, the last term goes to 0.
Now, taking $\{v^{\epsilon}\}=v$ as recovery sequence, we have that
$\mathcal J^{\epsilon,\delta}_{0,i}(v)\to \mathcal J^{\epsilon,\delta}(v)$, which completes
the proof.
\end{proof}

Performing analogous computations as in the previous section, we can see that
$s^{\epsilon,\delta}_i(t)$ is differentiable at $t=0$ and
$$
\frac{ds^{\epsilon,\delta}_i(0)}{dt}=R^{\epsilon,\delta}_{i}(u^{\epsilon, \delta}_i)
- \int_{D^{\delta}_i} u^{\epsilon, \delta}_i f \operatorname{div} V\,dy
- \int_{D^{\delta}_i} u^{\epsilon, \delta}_i \nabla f^{\epsilon}\cdot V\,dy
$$
where
$$
R^{\epsilon,\delta}_{i}(v)
:=\int_{D^{\delta}_i}\frac{(|\nabla v|^2+\epsilon^2)^\frac{p_i}{2}}{p_i}\operatorname{div} V
-(|\nabla v|^2+\epsilon^2)^{\frac{p_i}{2}-1} \nabla v\cdot\nabla v DV
+ \operatorname{div} V\frac{|v|^2}{2}\,dy.
$$
Since the expression of $\frac{ds^{\epsilon,\delta}_{i}(0)}{dt}$ given above
only involves first derivatives, we can conclude the following result
from Corollary \ref{convergenciasolu}.

\begin{proposition}\label{convergenciaderivada}
$\frac{ds^{\epsilon,\delta}_{i}(0)}{dt}$ converges to
$\frac{ds^{0,\delta}_{i}(0)}{dt}$ when $\epsilon \to 0$.
\end{proposition}

Observe that, by Propositions \ref{convergenciasolu} and \ref{convergenciaderivada},
 if we find an expression for the shape derivative of the regularized operator,
we will have found one for the original operator.

\subsection{Improvement of the formula for the regularized operator}
Our main concern in this part of our work will be to find a formula for
the shape derivative that does not involve second order derivatives.
Therefore, we will be able to pass to the limit when $\epsilon$ goes to $0$.
And so, by Proposition \ref{convergenciaderivada}, we will have found an
expression for the shape derivative of the original operator.

We start with some preliminaries computations in which we will see the
need to have $C^2$ regularity for our solutions.
Since
\begin{align*}
&\operatorname{div}((|\nabla u^{\epsilon}|^2+\epsilon^2)^{\frac{p_i}{2}}\cdot V)\\
&=\frac{p_i}{2}(|\nabla u^{\epsilon}|^2+\epsilon^2)^{\frac{p_i}{2}-1}
 D(|\nabla u^{\epsilon}|^2+\epsilon^2) \cdot V+(|\nabla u^{\epsilon}|^2
 +\epsilon^2)^{\frac{p_i}{2}}\operatorname{div} V\\
&=p_i(|\nabla u^{\epsilon}|^2+\epsilon^2)^{\frac{p_i}{2}-1}
 \nabla u^{\epsilon} D^2u^{\epsilon} \cdot V
 +(|\nabla u^{\epsilon}|^2+\epsilon^2)^{\frac{p_i}{2}}\operatorname{div} V,
\end{align*}
we have that
\begin{align*}
&\frac{1}{p_i}\int_{D^{\delta}_i}(|\nabla u^{\epsilon}|^2+\epsilon^2)^{\frac{p_i}{2}}
 \operatorname{div} V \\
&=\frac{1}{p_i}\int_{D^{\delta}_i}\operatorname{div}((|\nabla u^{\epsilon}|^2
 +\epsilon^2)^{\frac{p_i}{2}}V)- \int_{D^{\delta}_i}
 (|\nabla u^{\epsilon}|^2+\epsilon^2)^{\frac{p_i}{2}-1}
 \nabla u^{\epsilon} D^2u^{\epsilon} \cdot V
\end{align*}
Therefore,
\begin{align*}
\frac{ds^{\epsilon,\delta}_i}{dt}(0)
&=\frac{1}{p_i}\int_{D^{\delta}_i}\operatorname{div}((|\nabla u^{\epsilon}|^2+\epsilon^2
 )^{\frac{p_i}{2}}V)- \int_{D^{\delta}_i}(|\nabla u^{\epsilon}|^2
 +\epsilon^2)^{\frac{p_i}{2}-1} \nabla u^{\epsilon} D^2u^{\epsilon} \cdot V \\
&\quad -\int_{D^{\delta}_i}(|\nabla u^{\epsilon}|^2+\epsilon^2)^{\frac{p_i}{2}-1}
  \nabla u^{\epsilon} \cdot \nabla u^{\epsilon} DV\,dy+\frac{1}{2}\int_{D^{\delta}_i}
  \operatorname{div} (|u^{\epsilon}|^2V)\,dy \\
&\quad -\int_{D^{\delta}_i}u^{\epsilon}\nabla u^{\epsilon} \cdot V \,dy
 -\int_{D^{\delta}_i} u^{\epsilon} f^{\epsilon} \operatorname{div} V\,dy
 - \int_{D^{\delta}_i} u^{\epsilon} \nabla f^{\epsilon}\cdot V\,dy.
\end{align*}
Let us call $\nu^{\delta}_i$ the exterior unit normal vector to
$\partial D^{\delta}_i$ and observe that, since $\operatorname{suppt} V \subset\subset \Omega$,
$$
\int_{D^{\delta}_i} \operatorname{div} (|u^{\epsilon}|^2V)\,dy
 =\int_{\Gamma^{\delta}_i} |u^{\epsilon}|^2V \cdot \nu^{\delta}_i \,dS.
$$
Since $u^{\epsilon}$ is a weak solution of our equation, for every test function
$\varphi$ we have
\begin{align*}
&\int_{D^{\delta}_i}(|\nabla u^{\epsilon}|^2+{\epsilon}^2)^{\frac{p_i-2}{2}}
 \nabla u^{\epsilon}\nabla \varphi+\int_{D^{\delta}_i} u^{\epsilon} \varphi \\
&=\int_{D^{\delta}_i} f^{\epsilon} \varphi.
\end{align*}
Let us consider $\varphi=\nabla u^{\epsilon} \cdot V$ as a test function.
Since $\nabla (\nabla u^{\epsilon} \cdot V)=D^2u^{\epsilon} \cdot V^t
+\nabla u^{\epsilon} DV$, we obtain
\begin{align*}
&\int_{D^{\delta}_i}(|\nabla u^{\epsilon}|^2+{\epsilon}^2)^{\frac{p_i-2}{2}}
 \nabla u^{\epsilon}(D^2u^{\epsilon} \cdot V^t+\nabla u^{\epsilon} DV) \\
&=\int_{\Gamma^{\delta}_i}(|\nabla u^{\epsilon}|^2
 +{\epsilon}^2)^{\frac{p_i-2}{2}}\nabla u^{\epsilon} \cdot \eta
  \nabla u^{\epsilon} \cdot V+\int_{D^{\delta}_i}
  (f^{\epsilon}-u^{\epsilon}) \nabla u^{\epsilon} \cdot V.
\end{align*}
Since $V$ has compact support in $\Omega$, we arrive at
$$
\int_{\partial D^{\delta}_i}(|\nabla u^{\epsilon}|^2
+{\epsilon}^2)^{\frac{p_i-2}{2}}\nabla u^{\epsilon} \cdot
\eta \nabla u^{\epsilon} \cdot V
=\int_{\Gamma^{\delta}_i}(|\nabla u^{\epsilon}|^2
+{\epsilon}^2)^{\frac{p_i-2}{2}}\nabla u^{\epsilon} \cdot \nu_i^{\delta}
 \nabla u^{\epsilon} \cdot V.
$$
Therefore, taking into account that
$\nabla u^{\epsilon} D^2 u^{\epsilon} \cdot V =\nabla u^{\epsilon}
\cdot D^2 u^{\epsilon} V^T$, we have
\begin{align*}
\frac{ds^{\epsilon, \delta}_i}{dt}(0)
&= \frac{1}{p_i}
\int_{\Gamma_i^\delta}\Big(|\nabla u^{\epsilon}|^2
 + \epsilon^2\Big)^\frac{p_i}{2}V\nu_i^\delta
- \int_{\Gamma_i^\delta}\Big(|\nabla u^{\epsilon}|^2
+ \epsilon^2\Big)^{\frac{p_i}{2}-1} \nabla u^{\epsilon} \nu_i^\delta
\nabla u^{\epsilon} V \\
&\quad +\frac{1}{2}\int_{\Gamma_i^\delta}|u^{\epsilon}|^2V \nu_i^\delta
 - \int_{D_i^\delta} \Big(\underbrace{f^{\epsilon}\nabla u^{\epsilon} \cdot V
 + u^{\epsilon} f^{\epsilon} \operatorname{div} V + u^{\epsilon}\nabla f^{\epsilon}
 \cdot V}_{\operatorname{div}(u^{\epsilon}f^{\epsilon}V)}\Big)\,dy.
\end{align*}
Again since $V$ has compact support in $\Omega$, we have
$$
\int_{D^{\delta}_i}\operatorname{div} (u^{\epsilon}f^{\epsilon}V)\,dy
=\int_{\partial D^{\delta}_i}u^{\epsilon}f^{\epsilon}V \cdot \nu^{\delta}_i \,dS
=\int_{\Gamma^{\delta}_i}u^{\epsilon}f^{\epsilon}V \cdot \nu^{\delta}_i \,dS.
$$
Observe that we arrive at an expression for the shape derivative that does
not involve second order derivatives of $u^{\epsilon}$:
\begin{align*}
\frac{ds^{\epsilon, \delta}_i}{dt}(0)
&= \frac{1}{p_i} \int_{\Gamma_i^\delta}\Big(|\nabla u^{\epsilon}|^2
+ \epsilon^2\Big)^\frac{p_i}{2}V\nu_i^\delta
- \int_{\Gamma_i^\delta}\Big(|\nabla u^{\epsilon}|^2
+ \epsilon^2\Big)^{\frac{p_i}{2}-1} \nabla u^{\epsilon}
\nu_i^\delta \nabla u^{\epsilon} V\\
&\quad +\frac{1}{2}\int_{\Gamma_i^\delta}|u^{\epsilon}|^2V \nu_i^\delta
 -\int_{\Gamma^{\delta}_i}u^{\epsilon}f^{\epsilon}V \cdot \nu^{\delta}_i \,dS.
\end{align*}

\subsection{Back to the original operator: the limit when $\epsilon$ approaches $0$}
Now we able to apply Tolksdorf's regularity estimates (see \cite{Tolksdorf}).
These estimates give us uniform bounds for $||u^{\epsilon}||_{C^{1,\alpha}}$
so we have $u^{\epsilon}\to u$ in $C^{1}$. And so we can pass to the limit
when $\epsilon$ goes to $0$. Therefore,
\begin{align*}
\frac{ds^{0, \delta}_i}{dt}(0)
&= \frac{1}{p_i} \int_{\Gamma_i^\delta}|\nabla u|^{p_i}V\nu_i^\delta
 - \int_{\Gamma_i^\delta}|\nabla u|^{p_i-2} \nabla u \nu_i^\delta \nabla u V
 +\frac{1}{2}\int_{\Gamma_i^\delta}|u|^2V \nu_i^\delta \\
&\quad - \int_{\Gamma^{\delta}_i}ufV \cdot \nu^{\delta}_i \,dS.
\end{align*}
In conclusion we arrive at
\begin{align*}
\frac{ds^{0,\delta}}{dt}(0)
&=\frac{ds^{0,\delta}_1}{dt}(0)+\frac{ds^{0,\delta}_2}{dt}(0) \\
&=\frac{1}{p_1} \int_{\Gamma_1^\delta}|\nabla u|^{p_1}V\nu_1^\delta
 +\frac{1}{p_2} \int_{\Gamma_2^\delta}|\nabla u|^{p_2}V\nu_2^\delta
 - \int_{\Gamma_1^\delta}|\nabla u|^{p_1-2} \nabla u \nu_1^\delta \nabla u V \\
&\quad - \int_{\Gamma_2^\delta}|\nabla u|^{p_2-2} \nabla u \nu_2^\delta \nabla u V
 +\frac{1}{2}\int_{\Gamma_1^\delta}|u|^2V \nu_1^\delta
 - \int_{\Gamma^{\delta}_1}ufV \cdot \nu^{\delta}_1 \,dS \\
&\quad +\frac{1}{2}\int_{\Gamma_2^\delta}|u|^2V \nu_2^\delta
 - \int_{\Gamma^{\delta}_2}ufV \cdot \nu^{\delta}_2 \,dS.
\end{align*}
Let us now observe that $\nu_1^\delta \to \nu_1$ and
$\nu_2^\delta \to \nu_2=-\nu_1$ when $\delta \to 0$.
Therefore, taking limit when $\delta \to 0$, the last four terms in the
expression above vanish and so we have proved the following.

\begin{theorem}\label{improvement.formula}
Let $\Omega\subset \mathbb R ^N$ be open and bounded. Let $D_1, D_2\subset \Omega$
be such that \eqref{prop.particion} is satisfied, let
$p=p_1\chi_{D_1} + p_2\chi_{D_2}$, where $1<p_1<p_2$ and
$\Gamma=\bar{D_1}\cap\bar{D_2}$.
Let $V\colon \mathbb R ^N\to\mathbb R ^N$ be a Lipschitz deformation field, such that
$\operatorname{suppt}(V)\subset\subset\Omega$ and let $s(t)$ be defined
by \eqref{st}. Then, the following formula for the derivative $s'(0)$ holds:
$$
\frac{ds}{dt}(0)= \int_{\Gamma} \Big[\frac{|\nabla u|^p}{p}\Big]V \cdot \nu \,dS
- \int_{\Gamma}[|\nabla u|^{p-2}] (\nabla u\cdot \nu) (\nabla u\cdot V)\,dS,
$$
where
$$
\int_\Gamma [f] G\cdot \nu\,dS := \lim_{\delta\to 0}
\Big(\int_{\Gamma_1^\delta} f G\cdot\nu_1\,dS
 - \int_{\Gamma_2^\delta} f G\cdot \nu_2\,dS\Big).
$$
\end{theorem}


\section{Appendix: Gamma convergence results}

In this appendix we recall some basic concepts of $\Gamma-$convergence that
are needed in the present paper. Although these results are well-known,
we decide to include this appendix in order to make the paper self contained.
Also, the results presented here are not stated in the most general form,
 but in a for that will be enough for our work. For a complete presentation
of the theory of $\Gamma$-convergence, see the book of Dal Maso \cite{DalMaso}.

Let $\psi_n$ and $\psi$ defined in a topological space $X_\tau$ with $T^2$ topology.
 For our applications, $X_\tau$ will be a Banach space and we will consider
the weak topology. Then, a family of functionals $\psi_n$ $\Gamma$-converges
to $\psi$ if
\begin{itemize}
	\item (liminf inequality) $x_n \to_{\tau} x$ implies that
$\psi(x)\leq \liminf \psi_n(x_n)$
and
  \item (limsup inequality) there exists $y_n \rightharpoonup x$ such that
$\psi(x)\geq \limsup \psi_n(y_n)$.
\end{itemize}

\begin{theorem}\label{convergencia sobre acotados}
Let $X$ be a Banach space, and $C \subset X$ be closed and convex.
Let $\psi_n, \psi\colon C\to [-\infty,\infty]$ be weakly lower semicontinuous,
strictly convex and uniformly coercive functionals (i.e. for every $\lambda$,
the set $\{x \in C : \psi_n(x)\leq \lambda\}\subset B_r$ for every $n$)
such that $\psi_n$ $\Gamma$-converges to $\psi$.
Then $\inf_{C}\psi_n=\min_{C}\psi_n \to \inf_{C}\psi=\min_{C}\psi$.

Moreover, if $x_n \in C$ is such that $\psi_n(x_n)=\min_{C}\psi_n$, then
$	\{x_n\}_{n\in\mathbb N }$ is precompact and $\psi(x_0)=\min_{C}\psi$ where $x_0$
is any accumulation point of the sequence $\{x_n\}_{n\in\mathbb N }$.
\end{theorem}

\begin{proof}
Let us start by observing that, since $\psi_n$ weakly lower semicontinuous,
strictly convex and uniformly coercive functionals, for every $n$ there is
 a unique $x_n$ such that $\psi_n(x_n)=\inf_{C}\psi_n$ and $(x_n)$ is bounded
if $\psi_n(x_n)$ is bounded.
Let us consider now the following recovery function: $x \in C$ such that
$y_n \rightharpoonup x$. Therefore,
$$
\psi_n(x_n)=\inf_{C}\psi_n\leq \psi_n(y_n).
$$
So for every $x$ we have
$$
\limsup \psi_n(x_n)\leq \limsup \psi_n(y_n)\leq \psi(x).
$$
Therefore,
$$
\limsup \psi_n(x_n)\leq \inf_{C}\psi<\infty
$$
and we can conclude that $x_n \in \{x \in C : \psi_n(x)\leq \lambda\}\subset B_r$
for every $n\geq n_0$ taking $\lambda=\inf_{C}\psi+1$. So $(x_n)$ is bounded and,
via subsequences if necessary, $x_n\rightharpoonup x_0 \in C$
(remember that $C$ is convex and closed, therefore weakly closed)).

Finally, observing that
$$
\inf_{C}\psi\leq \psi(x_0)\leq \liminf \psi_n(x_n)\leq \liminf (\inf_{C}\psi_n),
$$
the proof is complete.
\end{proof}

\begin{remark}\label{obs gamma convergencia} \rm
If $\psi_n \to \psi$ point-wise, the inequality of the inferior limit always
holds (it is enough to take $y_n$ equal to $x$ for every $n$).
Therefore, to obtain the convergence of the functionals it would only be
necessary to check the superior limit inequality.
\end{remark}


\subsection*{Acknowledgements}

This research was partially supported by grants UBACyT 20020130100283BA,
CONICET PIP 11220150100032CO and ANPCyT PICT 2012-0153.
J. Fern\'andez Bonder is a member of CONICET.


\begin{thebibliography}{00}

\bibitem{Acerbi-Fusco} E.~Acerbi, N.~Fusco;
 \emph{A transmission problem in the calculus of
  variations}, Calc. Var. Partial Differential Equations \textbf{2} (1994),
  no.~1, 1--16. \MR{1384391}

\bibitem{Banichuk} N.~V. Banichuk;
 \emph{Introduction to optimization of structures},
  Springer-Verlag, New York, 1990, Translated from the Russian by Vadim Komkov.
  \MR{1101808}

\bibitem{Chambolle-Lions} Antonin Chambolle, Pierre-Louis Lions;
 \emph{Image recovery via total   variation minimization and related problems},
 Numer. Math. \textbf{76}   (1997), no.~2, 167--188. \MR{1440119}

\bibitem{CLR} Yunmei Chen, Stacey Levine,  Murali Rao;
 \emph{Variable exponent, linear   growth functionals in image restoration},
SIAM J. Appl. Math. \textbf{66}  (2006), no.~4, 1383--1406 (electronic).
\MR{2246061 (2007d:94004)}

\bibitem{DalMaso} Gianni Dal~Maso;
\emph{An introduction to {$\Gamma$}-convergence}, Progress in
  Nonlinear Differential Equations and their Applications, 8, Birkh\"auser
  Boston, Inc., Boston, MA, 1993. \MR{1201152}

\bibitem{Diening} Lars Diening, Petteri Harjulehto, Peter H\"ast\"o,  Michael
  R$\mathaccent"7017 {\rm u}$\v{z}i\v{c}ka;
 \emph{Lebesgue and {S}obolev spaces with   variable exponents},
Lecture Notes in Mathematics, vol. 2017, Springer,
  Heidelberg, 2011. \MR{2790542}

\bibitem{FanyZhao} Xianling Fan, Dun Zhao;
 \emph{On the spaces {$L^{p(x)}(\Omega)$} and   {$W^{m,p(x)}(\Omega)$}},
J. Math. Anal. Appl. \textbf{263} (2001), no.~2,
  424--446. \MR{1866056 (2003a:46051)}

\bibitem{GT} David Gilbarg, Neil~S. Trudinger;
 \emph{Elliptic partial differential  equations of second order},
 Classics in Mathematics, Springer-Verlag, Berlin,
  2001, Reprint of the 1998 edition. \MR{1814364 (2001k:35004)}

\bibitem{Hadamard} Jaques Hadamard;
 \emph{M\'emoire sur le probl\`eme d'analyse relatif \`a
  l'equilibre des plaques \'elastiques encastr\'ees}, Bull. Soc. Math. France,
  Nabu Press, 1907.

\bibitem{Pironneau} Olivier Pironneau;
 \emph{Optimal shape design for elliptic systems}, Springer
  Series in Computational Physics, Springer-Verlag, New York, 1984. \MR{725856
  (86e:49003)}

\bibitem{Sokolowski-Zolesio} Jan Soko\l~owski, Jean-Paul Zol\'esio;
 \emph{Introduction to shape  optimization}, Springer Series in Computational
Mathematics, vol.~16,  Springer-Verlag, Berlin, 1992, Shape sensitivity analysis.
 \MR{1215733}

\bibitem{Tolksdorf} Peter Tolksdorf;
 \emph{Regularity for a more general class of quasilinear   elliptic equations},
J. Differential Equations \textbf{51} (1984), no.~1,   126--150. \MR{727034}

\end{thebibliography}

\end{document}