# A Penalization Method for the Elliptic Bilateral Obstacle Problem

Cornel Marius Murea, Dan Tiba

► To cite this version:

Cornel Marius Murea, Dan Tiba. A Penalization Method for the Elliptic Bilateral Obstacle Problem. 26th Conference on System Modeling and Optimization (CSMO), Sep 2013, Klagenfurt, Austria. pp.189-198, 10.1007/978-3-662-45504-3_18 . hal-01286413

HAL Id: hal-01286413

https://inria.hal.science/hal-01286413

Submitted on 10 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A penalization method for the elliptic bilateral obstacle problem

Cornel Marius Murea<fnref n="1" /> and Dan Tiba<fnref n="2" />,<fnref n="3" />*

**Abstract.** In this paper we propose a new algorithm for the well-known elliptic bilateral obstacle problem. Our approach enters the category of fixed domain methods and solves just linear elliptic equations at each iteration. The approximating coincidence set is explicitly computed. In the numerical examples, the algorithm has a fast convergence.

**Keywords:** obstacle problem, free boundary problems, penalization

## 1 Introduction

The obstacle problem may be formulated as an elliptic variational inequality. Detailed theoretical discussions of various variational inequalities may be found in [4, 15, 23]. Applications, including optimal control problems are investigated in the books [1, 5, 6, 24]. From the point of view of the numerical approximation, we quote just the monographs [6, 7, 21].

In this paper we propose an algorithm for the elliptic bilateral obstacle problem which is of fixed domain type in the sense that the finite element discretization is given in the whole domain, independently of the position of the unknown free boundary. In each iteration a linear elliptic equation has to be solved in the whole domain and the corresponding stiffness matrix is common for all iterations. This is a clear advantage from the point of view of the implementation and the approximating coincidence set is explicitly computed in each iteration and it converges in the Hausdorff-Pompeiu sense [20] to the searched geometry. Moreover, we need just a scalar penalization parameter in our method. A similar strategy was employed in [18] for the elliptic unilateral obstacle problem and for parabolic variational inequalities.

<footnote n="1">¹ Laboratoire de Mathématiques, Informatique et Applications,
Université de Haute Alsace
4-6, rue des Frères Lumière, 68093 Mulhouse Cedex, France
cornel.murea@uha.fr
http://www.edp.lmia.uha.fr/murea/</footnote>

<footnote n="2">² Institute of Mathematics of the Romanian Academy
P.O. Box 1-764, 014700 Bucharest, Romania,</footnote>

<footnote n="3">³ Academy of Romanian Scientists, Bucharest, Romania
dan.tiba@imar.ro</footnote>

<footnote>* Supported by Grant 145/2011 CNCS, Romania.</footnote>
Our approach is inspired from shape optimization techniques, but no shape optimization problem is used here although this is a known method in free boundary problems, [2]. One may compare the present approach to the recent works [8, 19, 22]. An efficient Lagrangian method together with a primal-dual active set strategy with regularization is studied in [13]. But our approach and arguments are certainly different. We also quote the multi grid method employed in [12], the path-following method for semi-smooth Newton schemes [10] and a duality-type method [16, 17].

## 2 Formulation of the problem and the algorithm

Let $D$ be a smooth domain in $\mathbb{R}^d$, $d \in \mathbb{N}^*$ and $f \in L^2(D)$ be given. We denote the obstacles by $\psi_1, \psi_2: D \to \mathbb{R}$, such that $\psi_1, \psi_2 \in H^2(D)$, $\psi_1 \le \psi_2$ in $D$, $\psi_{1|\partial D} \le 0$, $\psi_{2|\partial D} \ge 0$. The admissible set

$$K = \{v \in H_0^1(D); \psi_1(x) \le v(x) \le \psi_2(x) \text{ a.e. in } D\}$$

is a nonvoid closed convex subset of $H_0^1(D)$.

To $K$, the following variational inequality, may be associated:

$$\int_D \nabla y \cdot (\nabla y - \nabla v) dx \le \int_D f(y-v) dx, \quad \forall v \in K. \quad (1)$$

The existence of a unique solution $y \in K$ is well-known.

We introduce $\beta \subset \mathbb{R} \times \mathbb{R}$ the maximal monotone graph given by

$$\beta(r) = \begin{cases} \emptyset, & r < 0, \\ ]-\infty, 0], & r = 0, \\ 0, & r > 0, \end{cases} \quad (2)$$

the maximal monotone graph $\gamma \subset \mathbb{R} \times \mathbb{R}$ given by

$$\gamma(r) = \begin{cases} 0, & r < 0, \\ [0, +\infty[, & r = 0, \\ \emptyset, & r > 0 \end{cases} \quad (3)$$

and denote by $\beta_\epsilon$, $\gamma_\epsilon$, $\epsilon > 0$, their Yosida approximations. We have

$$\beta_{\epsilon}(r) = \begin{cases} \frac{1}{\epsilon}r, & r \le 0, \\ 0, & r > 0, \end{cases} \quad \gamma_{\epsilon}(r) = \begin{cases} 0, & r < 0, \\ \frac{1}{\epsilon}r, & r \ge 0. \end{cases}$$

Notice that $\beta_\epsilon$ is a concave and $\gamma_\epsilon$ is a convex function in $\mathbb{R}$.

In the case when $f \in L^2(D)$, $\psi_1, \psi_2 \in H^2(D)$ with the compatibility condition $\psi_1 \le \psi_2$ in $D$, $\psi_{1|\partial D} \le 0$, $\psi_{2|\partial D} \ge 0$, it is known that the solution of (1) satisfies the regularity property $y \in H^2(D)$. Moreover, in this case, the obstacle problem may be written as a multivalued equation

$$-\Delta y + \beta(y - \psi_1) + \gamma(y - \psi_2) \ni f \text{ in } D. \quad (4)$$
One can define two coincidence sets, corresponding to the two obstacles:

$$D_1 = \{x \in D; y(x) = \psi_1(x)\}$$

$$D_2 = \{x \in D; y(x) = \psi_2(x)\}$$

and associated to (1).

We state now our algorithm.

**Algorithm**

1) Choose $n = 0, \epsilon_0 > 0, \Omega_1^0 \subset D, \Omega_2^0 \subset D$ open subsets such that $(D \setminus \Omega_1^0) \cap (D \setminus \Omega_2^0) = \emptyset, \tilde{y}_{-1} = 0$;

2) Compute $y_n \in H_0^1(D)$ as solution of the linear elliptic equation

$$-\Delta y_n + \frac{1}{\epsilon_n} \chi_{D \setminus \Omega_1^n} (y_n - \psi_1) + \frac{1}{\epsilon_n} \chi_{D \setminus \Omega_2^n} (y_n - \psi_2) = f \quad \text{in } D \quad (5)$$

3) Compute $\mathbf{y}_n = \min\{\psi_2, \max\{y_n, \psi_1\}\}, \Omega_1^{n+1} = \{x \in D; \mathbf{y}_n(x) > \psi_1(x)\}, \Omega_2^{n+1} = \{x \in D; \mathbf{y}_n(x) < \psi_2(x)\} \epsilon_{n+1} = \frac{\epsilon_n}{2}$;

4) If $\|\mathbf{y}_n - \mathbf{y}_{n-1}\|_{H^1(D)} < \text{tol}$ then STOP else n=n+1 GO TO step 2.

*Remark 1.* By the classical result of [3], the elastic-plastic torsion problem is equivalent with a variational inequality of obstacle type and our algorithm may be applied as well.

We convene to extend the value $\frac{1}{\epsilon}$ for $\beta'_\epsilon$ and $\gamma'_\epsilon$ in the origin and we can rewrite the step 2 of the Algorithm as

$$-\Delta y_n + (\beta'_{\epsilon_n}(y_{n-1} - \psi_1))(y_n - \psi_1) + (\gamma'_{\epsilon_n}(y_{n-1} - \psi_2))(y_n - \psi_2) = f. \quad (6)$$

Recall that the usual approximation by regularization of the variational inequality (1) is

$$-\Delta \tilde{y}_n + \beta_{\epsilon_n}(\tilde{y}_n - \psi_1) + \gamma_{\epsilon_n}(\tilde{y}_n - \psi_2) = f \quad \text{in } D, \quad (7)$$

plus homogeneous boundary conditions on $\partial D$.

Notice that $\beta_\epsilon(r) = \beta'_\epsilon(r)r$ and $\gamma_\epsilon(r) = \gamma'_\epsilon(r)r$, under the above convention, which shows that (6) and (7) have very similar structure. Clearly, (7) is a non-linear elliptic equation, while the decoupling operated in (6) allows to use linear elliptic equations.

## 3 Stability

We present in this section a stability result in $L^2(D)$ for the algorithm introduced above applied to the bilateral obstacle problem.

**Theorem 1.** i) The sequence $\{y_n\}$ is bounded in $L^2(D)$.
ii) There is $C > 0$, independent of $n$, such that:

$$ \int_{D \setminus \Omega_2^n} (y_n - \psi_2)_+^2 dx + \int_{D \setminus \Omega_1^n} (y_n - \psi_1)_-^2 dx \leq C \epsilon_n. \quad (8) $$

*Proof.* Using $\beta_{\epsilon}(r) = \beta'_{\epsilon}(r)r$, the concavity of $\beta_{\epsilon}(\cdot)$ and the definition of the subdifferential of concave mapping, we obtain

$$ \begin{aligned} & (\beta'_{\epsilon_n}(y_{n-1} - \psi))(y_n - \psi) = (\beta'_{\epsilon_n}(y_{n-1} - \psi))(y_{n-1} - \psi) \\ & + (\beta'_{\epsilon_n}(y_{n-1} - \psi))(y_n - \psi - y_{n-1} + \psi) \geq \beta_{\epsilon_n}(y_{n-1} - \psi) \\ & + \beta_{\epsilon_n}(y_n - \psi) - \beta_{\epsilon_n}(y_{n-1} - \psi) = \beta_{\epsilon_n}(y_n - \psi). \end{aligned} $$

We use the above inequality in the equation (5). We get

$$ -\Delta y_n + \beta_{\epsilon_n}(y_n - \psi_1) + \frac{1}{\epsilon_n} \chi_{D \setminus \Omega_2^n} (y_n - \psi_2) \leq f, \quad (9) $$

where $\beta$ is given by (2) and $\beta_{\epsilon_n}$ is its regularization.

We multiply (9) by $(y_n - \psi_2)_+$ and we use that

$$ \beta_{\epsilon_n}(y_n - \psi_1)(y_n - \psi_2)_+ = 0. \quad (10) $$

While $\beta_{\epsilon_n}(y_n - \psi_1)$ may take negative values, this happens for $y_n \le \psi_1$, that is $y_n - \psi_2 \le 0$ (since $\psi_1 \le \psi_2$). Then $(y_n - \psi_2)_+ = 0$ and (10) follows. We infer

$$ \begin{aligned} & \int_D |\nabla(y_n - \psi_2)_+|^2 + \frac{1}{\epsilon_n} \int_{D \setminus \Omega_2^n} (y_n - \psi_2)_+^2 \quad (11) \\ & \le \int_D f(y_n - \psi_2)_+ + \int_D \nabla\psi_2 \cdot \nabla(y_n - \psi_2)_+. \end{aligned} $$

By the conditions $\psi_2|_{\partial D} \ge 0$ and $y_n|_{\partial D} = 0$, we have $(y_n - \psi_2)_+ = 0$ on $\partial D$ and the Poincaré inequality shows that $\{(y_n - \psi_2)_+\}$ is bounded in $H_0^1(D)$, by (11).

Equation (5) may be rewritten in the form

$$ -\Delta y_n + \frac{1}{\epsilon_n} \chi_{D \setminus \Omega_1^n} (y_n - \psi_1) + \gamma'_{\epsilon_n} (y_{n-1} - \psi_2)(y_n - \psi_2) = f. \quad (12) $$

We compute

$$ \begin{aligned} & \gamma'_{\epsilon_n}(y_{n-1} - \psi_2)(y_n - \psi_2) = \gamma'_{\epsilon_n}(y_{n-1} - \psi_2)(y_{n-1} - \psi_2) \quad (13) \\ & + \gamma'_{\epsilon_n}(y_{n-1} - \psi_2)(y_n - \psi_2 - y_{n-1} + \psi_2) = \gamma_{\epsilon_n}(y_{n-1} - \psi_2) \\ & + \gamma'_{\epsilon_n}(y_{n-1} - \psi_2)(y_n - \psi_2 - y_{n-1} + \psi_2) \le \gamma_{\epsilon_n}(y_{n-1} - \psi_2) \\ & + \gamma_{\epsilon_n}(y_n - \psi_2) - \gamma_{\epsilon_n}(y_{n-1} - \psi_2) = \gamma_{\epsilon_n}(y_n - \psi_2) \end{aligned} $$

using the subdifferential of convex mappings.
By (12), (13), we obtain

$$ -\Delta y_n + \frac{1}{\epsilon_n} \chi_{D \setminus \Omega_1^n} (y_n - \psi_1) + \gamma_{\epsilon_n} (y_n - \psi_2) \geq f. \quad (14) $$

Multiply (14) by $-(y_n - \psi_1)_- \in H_0^1(D)$, due to $y_n = 0$ on $\partial D$, $\psi_1 \le 0$ on $\partial D$:

$$ \int_D |\nabla(y_n - \psi_1)_-|^2 + \frac{1}{\epsilon_n} \int_{D \setminus \Omega_1^n} (y_n - \psi_1)_-^2 - \gamma_{\epsilon_n} (y_n - \psi_2)(y_n - \psi_1)_- \quad (15) \\ \le -\int_D f(y_n - \psi_1)_- + \int_D \nabla \psi_1 \cdot \nabla (y_n - \psi_1)_-. $$

Notice that

$$ -\gamma_{\epsilon_n} (y_n - \psi_2)(y_n - \psi_1)_- = 0 \quad (16) $$

since $-\gamma_{\epsilon_n}(y_n - \psi_2)$ may take negative values just for $y_n \ge \psi_2 \ge \psi_1$ and in this case $(y_n - \psi_1)_- = 0$. By (15), (16) we obtain:

$$ \int_D |\nabla(y_n - \psi_1)_-|^2 + \frac{1}{\epsilon_n} \int_{D \setminus \Omega_1^n} (y_n - \psi_1)_-^2 \quad (17) \\ \le -\int_D f(y_n - \psi_1)_- + \int_D \nabla \psi_1 \cdot \nabla (y_n - \psi_1)_-. $$

Relation (17) shows that $\{(y_n - \psi_1)_-\}$ is bounded in $H_0^1(D)$, by the Poincaré inequality.

We use the inequality

$$ (x - b)_+ \le (x - a)_+ + (a - b)_+ $$

and we have

$$ (y_n - \psi_1)_+ \le (y_n - \psi_2)_+ + (\psi_2 - \psi_1)_+ = (y_n - \psi_2)_+ + \psi_2 - \psi_1. $$

Relation (11) shows that $\{(y_n - \psi_1)_+\}$ is bounded in $L^2(D)$. In combination with (17), it yields $\{y_n\}$ bounded in $L^2(D)$.

Relation (8) follows by adding (11) and (17) and using the already established boundedness of all the terms except the penalization term. This ends the proof.

*Remark 2.* In fact, the above proof shows that $\{y_n\}$ bounded in $L^p(D)$, $p > 2$ depending on the dimension of $D$. Relation (8) says that the sequence $\{y_n\}$ does not overpass the obstacles $\psi_1, \psi_2$, in the limit. The proof also provides partial information on $\{\nabla y_n\}$, but it is unclear whether $\{y_n\}$ is bounded in $H_0^1(D)$.

# 4 Numerical tests

We have used the software FreeFem++ v 3.19, [9]. For all tests, we use the same initial guess for the coincidence set $D \setminus \Omega_1^0 = \emptyset$ and $D \setminus \Omega_2^0 = \emptyset$.
## Test 1.

We consider the torsion of an elastic-plastic prism studied in [7] p. 133 and [25]. The cross-section of the prism is $D = [0, 1] \times [0, 1]$. We solve the problem (1) with $f(x) = -8$ where $K = \{v \in H_0^1(D); -1 \le \nabla v(x) \le 1 \text{ a.e. in } D\}$. For $v \in K$, the set $\{x \in D; |\nabla v(x)| < 1\}$ is the elastic zone and $\{x \in D; |\nabla v(x)| = 1\}$ is the plastic zone, [5, p. 264]. By the result of [3], the elastic-plastic torsion problem is equivalent with a variational inequality of obstacle type $K = \{v \in H_0^1(D); \psi_1(x) \le v(x) \le \psi_2(x) \text{ a.e. in } D\}$ where $\psi_1(x) = -\text{dist}(x, \partial D)$ and $\psi_2(x) = \text{dist}(x, \partial D)$. If $f < 0$, then $y \le 0$, consequently the top obstacle will be inactive.

We use a mesh of 39158 triangles, 19836 vertices and size $h = \frac{1}{128}$. The tolerance for the stopping test is $tol = 10^{-3}$ and the penalization parameter is $\epsilon_n = 0.003$. The coincidence set of the solution presented in Figure 1 at the right is similar to the above references.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Test 1. The computed coincidence set of the plastic zone for the bottom obstacle at the first (left), second (middle) and last (right) iteration.</figcaption></figure>

Our algorithm stops after 6 iterations and the relative error in the $H^1$ norm at the last iteration is $\|y_n - y_{n-1}\|_{H^1(D)} = 1.5 \times 10^{-5}$.

In [18], we have tested numerically with positive results the stability of a similar algorithm when $f$ the right-hand side in (1) is perturbed.

## Test 2.

We solve the problem (1) where

$$K = \{v \in H_0^1(D); \psi_1(x) \le v(x) \le \psi_2(x) \text{ a.e. in } D\},$$

$D = [0, 1] \times [0, 1]$, $\psi_1(x) = -\text{dist}(x, \partial D)$, $\psi_2(x) = \text{dist}(x, \partial D)$ and $f(x) = 11(x + y - 1)$. Now both obstacles are active.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Test 2. Coincidence sets of the plastic zone for the bottom obstacle (left) and for the top obstacle (right).</figcaption></figure>

We use a mesh of 39158 triangles, 19836 vertices, the size $h = \frac{1}{128}$, the tolerance for the stopping test tol = $10^{-3}$ and the penalization parameter is $\epsilon_n = 0.003$. The algorithm stops in 4 iterations and the relative error in the $H^1$ norm at the last iteration is $\|y_n - y_{n-1}\|_{H^1(D)} = 0.000209$.

The coincidence sets are presented in Figure 2 and the computed solution in Figure 3.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Test 2. Computed solution.</figcaption></figure>

We solved the problem on different meshes, see Table 1. We denote by $u_i$, the solution obtained using the mesh no. $i$.

In [11], for the semi-smooth Newton method, it is proved a mesh-independence result: the continuous and the discrete process, converge q-linearly with the same rate.
<table><caption>**Table 1.** Test 2. Mesh parameters.</caption><thead><tr><th>mesh no.</th><th>mesh size h</th><th>triangles</th><th>vertices</th><th>||u<sub>i</sub> - u<sub>i-1</sub>||<sub>H<sup>1</sup>(D)</sub></th></tr></thead><tbody><tr><td>1</td><td>1/64</td><td>9720</td><td>4989</td><td>-</td></tr><tr><td>2</td><td>1/128</td><td>39158</td><td>19836</td><td>0.019005</td></tr><tr><td>3</td><td>1/256</td><td>154050</td><td>77538</td><td>0.008826</td></tr><tr><td>4</td><td>1/512</td><td>630326</td><td>316188</td><td>0.005768</td></tr></tbody></table>

**Test 3.**

Now we test the algorithm for the torsion of the elastic-plastic prism discussed in [14] and [26]. We can put this problem in the form (1). Let $D = [0, 1] \times [0, 1]$, $\psi_1(x, y) = -\text{dist}((x, y), \partial D)$, $\psi_2(x, y) = 0.2$ for all $(x, y) \in D$ and set

$$g(x) = \begin{cases} 6x, & 0 < x \le 1/6, \\ 2(1 - 3x), & 1/6 < x \le 1/3, \\ 6(x - 1/3), & 1/3 < x \le 1/2, \\ 2(1 - 3(x - 1/3)), & 1/2 < x \le 2/3, \\ 6(x - 2/3), & 2/3 < x \le 5/6, \\ 2(1 - 3(x - 2/3)), & 5/6 < x \le 1 \end{cases}$$

and

$$f(x, y) = \begin{cases} 300, & (x, y) \in S = \{(x, y) \in D; |x - y| \le 0.1 \ \& \ x \le 0.3\}, \\ -70 \exp(y)g(x), & x \le 1 - y \text{ and } (x, y) \notin S, \\ 15 \exp(y)g(x), & x > 1 - y \text{ and } (x, y) \notin S. \end{cases}$$

We use a mesh of 39158 triangles, 19836 vertices, the size $h = \frac{1}{128}$, the tolerance for the stopping test $\text{tol} = 10^{-3}$ and the penalization parameter is $\epsilon_n = 0.03$. The computed solution after 6 iterations is presented in Figure 4 and the corresponding coincidence sets in Figure 5. The relative error in the $H^1$ norm at the last iteration is $||y_n - y_{n-1}||_{H^1(D)} = 2.7 \times 10^{-5}$.

In [14], an augmented lagrangian active set strategy is employed. At each iteration, a reduced linear system associated with the inactive set is solved. In [26], at each iteration, linear systems associated to the complementary of the coincidence sets are solved.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Test 3. Computed solution.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Test 3. Coincidence sets for the bottom obstacle (left) and for the top obstacle (right).</figcaption></figure>

References

1. Barbu, V.: Optimal control of variational inequalities. Research Notes in Mathematics, 100. Pitman, Boston (1984)

2. Burger, M., Matevosyan, N., Wolfram, M.T.: A level set based shape optimization method for an elliptic obstacle problem. Math. Models Methods Appl. Sci. 21, no. 4, 619–649 (2011)

3. Brézis, H., Sibony, M.: Équivalence de deux inéquations variationnelles et applications. (French) Arch. Rational Mech. Anal. 41, 254–265 (1971)

4. Brézis, H.: Problèmes unilatéraux, J. Math. Pures Appl. (9) 51, 1–168 (1972)

5. Duvaut, G., Lions, J.-L.: Les inéquations en mécanique et en physique. Travaux et Recherches Mathématiques, No. 21. Dunod, Paris (1972)
6. Elliott, C. M., Ockendon, J. R.: Weak and variational methods for moving boundary problems. Research Notes in Mathematics, 59. Pitman, London (1982)

7. Glowinski, R., Lions, J.-L., Trémolières, R.: Numerical analysis of variational inequalities. Studies in Mathematics and its Applications, 8. North-Holland Publishing Co., Amsterdam-New York (1981)

8. Halanay, A., Murea, C.M., Tiba, D.: Existence and approximation for a steady fluid-structure interaction problem using fictitious domain approach with penalization, Mathematics and its Applications 5, no. 1-2, 120-147 (2013)

9. Hecht, F.: http://www.freefem.org

10. Hintermuller, M., Kunisch, K.: Path-following methods for a class of constrained minimization problems in function space. SIAM J. Optim. 17, no. 1, 159-187, (2006)

11. Hintermuller, M., Ulbrich, M.: A mesh-independence result for semismooth Newton methods, Math. Program., Ser. B 101, 151-184 (2004)

12. Hoppe, R. H. W., Kornhuber, R.: Adaptive multilevel methods for obstacle problems. SIAM J. Numer. Anal. 31, no. 2, 301-323, (1994)

13. Ito, K., Kunisch, K.: Semi-smooth Newton methods for variational inequalities of the first kind. M2AN Math. Model. Numer. Anal. 37, no. 1, 41-62 (2003)

14. Karkkainen, T., Kunisch, K., Tarvainen, P.: Augmented Lagrangian active set methods for obstacle problems. J. Optim. Theory Appl. 119, no. 3, 499-533 (2003)

15. Kinderlehrer, D., Stampacchia, G.: An introduction to variational inequalities and their applications. Reprint of the 1980 original. Classics in Applied Mathematics, 31, SIAM (2000)

16. Merlusca, D.: A duality algorithm for the obstacle problem. Annals of the Academy of Romanian Scientists, Vol. 5, No. 1-2, 209-215 (2013)

17. Merlusca, D.: A duality-type method for the fourth order obstacle problem, U.P.B. Sci. Bull., Series A, Vol. 76, No. 2, 147-158 (2014)

18. Murea, C.M., Tiba, D.: A direct algorithm in some free boundary problems, submitted to Journal of Numerical Mathematics (2014)

19. Neittaanmaki, P., Pennanen, A., Tiba, D.: Fixed domain approaches in shape optimization problems with Dirichlet boundary conditions. Inverse Problems, 25, no. 5, 1-18 (2009)

20. Neittaanmaki, P., Sprekels, J., Tiba, D.: Optimization of elliptic systems. Theory and applications. Springer Monographs in Mathematics. Springer, New York (2006)

21. Neittaanmaki P., Tiba D.: Optimal Control of Nonlinear Parabolic Systems. Theory, Algorithms, and Applications, Monographs and Textbooks in Pure and Applied Mathematics, 179. Marcel Dekker, Inc., New York, (1994)

22. Neittaanmaki, P., Tiba, D.: Fixed domain approaches in shape optimization problems, Inverse Problems, 28, no. 9, 1-35 (2012)

23. Rodrigues, J.-F.: Obstacle Problems in Mathematical Physics, North-Holland Publishing Co., Amsterdam, (1987)

24. Tiba, D.: Optimal Control of Nonsmooth Distributed Parameter Systems, Springer, Berlin, (1990)

25. Xue, L., Cheng, X.-L.: An algorithm for solving the obstacle problems. Comput. Math. Appl. 48, no. 10-11, 1651-1657 (2004)

26. Wang, F., Cheng, X.-L.: An algorithm for solving the double obstacle problems. Appl. Math. Comput. 201, no. 1-2, 221-228 (2008)