Approximate solutions 
======================

In the section, we will construct approximate solutions for .

Let $R>0$ be a large constant, such that for any $x\in {\Omega}$, ${\Omega}\subset\subset B_R(x)$. Consider the following problem:

$$\label{2.1}
\begin{cases}
-\delta^2\Delta w=( w-a)_+^{p},& \text{in}\; B_R(0),\\
w=0, &\text{on}\;\partial B_R(0),
\end{cases}$$

where $a>0$ is a constant. Then, has a unique solution $W_{{\delta},a}$, which can be written as

$$\label{2.2}
W_{{\delta},a}(x)=
\begin{cases}
a+{\delta}^{2/(p-1)}s_{\delta}^{-2/(p-1)}\phi\bigl(\frac{|x|}{s_{\delta}}\bigr), &  |x|\le s_{\delta},\\
a\ln\frac {|x|} R/\ln \frac {s_{\delta}}R, & s_{\delta}\le |x|\le R,
\end{cases}$$

where $\phi(x)=\phi(|x|)$ is the unique solution of $$-\Delta \phi=\phi^{p},\quad\phi>0,~~\phi\in H_0^1\bigl(B_1(0)\bigr)$$ and $s_{\delta}\in (0,R)$ satisfies $${\delta}^{2/(p-1)}s_{\delta}^{-2/(p-1)}\phi^\prime(1)=\frac{a}{\ln(s_{\delta}/R)},$$ which implies $$\frac{s_{\delta}}{{\delta}|\ln{\delta}|^{(p-1)/2}}\rightarrow\left(\frac{|\phi^\prime(1)|}{a}\right)^{(p-1)/2}>0,\quad\text{as}~~{\delta}\rightarrow0.$$ Moreover, by Pohozaev identity, we can get that $$\int_{B_1(0)}\phi^{p+1}=\frac{\pi(p+1)
}{2}|\phi^\prime(1)|^2~~\text{and}~~\int_{B_1(0)}\phi^{p}=2\pi|\phi^\prime(1)|.$$

For any $z \in {\Omega}$, define $W_{{\delta},z,a}(x)=W_{{\delta},a}(x-z)$. Because $W_{{\delta},z,a}$ does not vanish on $\partial\Omega$, we need to make a projection. Let $PW_{{\delta},z,a}$ be the solution of

$$\begin{cases}
-{\delta}^2 \Delta w=( W_{{\delta},z,a}-a)_+^{p},& \text{in } \; {\Omega},\\
w=0, &\text{on}\; \partial {\Omega}.
\end{cases}$$ Then

$$\label{2.3}
PW_{{\delta},z,a}= W_{{\delta},z,a}-\frac a{\ln \frac{R}{s_{\delta}}} g(x,z),$$

where $g(x,z)$ satisfies

$$\begin{cases}
- \Delta g=0,& \text{in } \; {\Omega},\\
g=\ln\frac{R}{|x-z|}, &\text{on}\; \partial {\Omega}.
\end{cases}$$ It is easy to see that

$$g(x,z)=\ln R +2\pi h(x,z),$$ where $h(x,z)=-H(x,z)$.

We will construct solutions for of the form

$$\sum_{j=1}^m PW_{{\delta},z_{j},a_{{\delta},j}} +\omega_{\delta},$$ where $z_j\in\Omega$ for $j=1,\cdots,m$, $\omega_{\delta}$ is a perturbation term. To obtain 