Numerical Verification Method of Solutions for Elliptic Variational Inequalities

Abstract

In this chapter, we propose numerical techniques which enable us to verify the existence of solutions for the free boundary problems governed by two kinds of elliptic variational inequalities. Based upon the finite element approximations and explicit a priori error estimates for some elliptic variational inequalities, we present effective verification procedures that, through numerical computation, generat a set which includes exact solutions. We describe a survey of the previous works as well as show newly obtained results up to now.

1. Introduction

Numerical verification methods of solutions for differential equations have been the subject of extensive study in recent years and much progress have been made both mathematically and computationally [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. However, for some problems governed by the elliptic variational inequalities, there are very few approaches. As far as we know, it is hard to find any applicable methods except for those of Nakao and Ryoo [13, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46].

The authors have studied for several years the numerical verification method of solutions for elliptic variational inequalities using finite element method and the constructive error estimates combining with Schauder’s and Banach’s fixed point theorem. Several results in our research are already published in [13, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. In this chapter, we briefly overview our resent research results including works not yet published.

The outline of this chapter is as follows. In Section 2, the two types of elliptic variational inequalities are considered. In Subsection 2.1, we describe the elliptic variational inequalities and give a fixed point formulation to prove the existence of solutions. In Subsections 2.2 and 2.3, the main tool of the verification method is explained at an abstract level. In Subsection 2.2, we present a simple iteration method for numerical verification of solutions for the elliptic variational inequalities. We construct the concepts of rounding and rounding error for functions and present a computer algorithm to construct the set satisfying the verification conditions. However, it is difficult to apply the method in Subsection 2.2 to a problem in which an associated operator is not retractive in a neighborhood of the solution, because it is based upon a simple iteration method. In Subsection 2.3, we propose another approach to overcome such a difficulty. This method can be applied to general elliptic variational inequalities without any retraction property of the associated operator. We introduce a Newton-like operator and reformulate the problem using it. Particularly, special emphasis is placed on the way to devise the Newton-like operator for a kind of non-differentiable map which defines the original problem. We introduce a computational verification condition. In order to show a concrete usage of the tool, in Section 3, we present an application to some problems governed by the elliptic variational inequalities. Many difficulties remain to be overcome in the construction of general techniques applicable to a broader range of problems. However, the authors have no doubt that investigation along this line will lead to a new approach employing numerical methods in the field of existence theory of solutions for various variational inequalities that appear in mathematical analysis.

2. Elliptic variational inequalities

The theory of elliptic variational inequalities has become a rich source of inspiration in both mathematical and engineering sciences. Elliptic variational inequalities are an effective tool for studying the existence of solutions of constrained problems arising in mechanics, optimization and control, operation research, engineering science, etc. [47, 48, 49, 50, 51, 52]. It is the aim of this chapter to introduce a numerical technique to verify the solutions for elliptic variational inequalities. The basic approach of this technique consists of the fixed point formulation of elliptic variational inequalities and construction of the function set, on computer, satisfying the validation condition of a certain infinite dimensional fixed point theorem. For fixed point formulation, we consider a candidate set which possibly contains a solution. In order to get such a candidate set, we divide the verification procedure into two phases: one is the computation of a projection into a closed convex subset of some finite dimensional subspace (rounding); the other is the estimation of the error for the projection (rounding error). Combining these methods with some iterative technique, the exact solution can be enclosed by sum of rounding parts, which is a subset of finite dimensional space, and the rounding error, which is indicated by a nonnegative real number. These two procedures enable us to treat infinite dimensional problems as finite procedures, thta is, by computer.

Notations

• $V$: real Hilbert space with scalar product $\left(\cdot ,\cdot \right)$ and associated norm $\parallel \cdot \parallel ,$

• $V^{\ast }$: the dual space of $V$,

• $a\left(\cdot ,\cdot \right):V\times V\to \mathbf{R}$ is a bilinear, continuous and $V$-elliptic from on $V\times V$.

A bilinear form $a\left(\cdot ,\cdot \right)$ is said to be $V$-elliptic if there exists a positive constant $\alpha$ such that $a\left(v,v\right)\ge \alpha \parallel v{\parallel }^2,\forall v\in V.$.

In general we do not assume $a\left(\cdot ,\cdot \right)$ to be symmetric, since in some applications nonsymmetric bilinear forms may occur naturally.

• $L:V\to \mathbf{R}$ continuous, linear functional,

• $K$is a closed convex nonempty subset of $V$,

• $j\left(\cdot \right):V\to \mathbf{R}\cup \left\{\infty \right\}$ is a convex lower semicontinuous (l.s.c) and proper functional ($j\left(\cdot \right)$ is proper if $j\left(v\right)>-\infty ,\forall v\in V$ and $j\ne +\infty )$.

The two types of elliptic variational inequalities.

We consider two classes of elliptic variational inequalities.

• Elliptic variational inequalities of the first kind: Find$u\in V$such that$u$is a solution of the problem

$$a\left(u,v-u\right)\ge L\left(v-u\right),\forall v\in K,u\in K.$$

• Elliptic variational inequalities of the second kind: Find$u\in V$such that$u$is a solution of the problem

$$a\left(u,-u\right)+j\left(v\right)-j\left(u\right)\ge L\left(v-u\right),\forall v\in V,u\in V.$$

2.1 The problem and the fixed point formulation

Let us first set a few notations [1, 47, 49, 50, 53, 54, 55, 56, 57, 58, 59, 60, 61]. In what follows we shall make use of the Sobolev spaces $W^{k,p}\left(\mathrm{\Omega }\right)$ of functions which possess generalized derivatives integrable with the $p$th power up to and including the $k$th order. For $p=2,$ we shall write $W^{k,p}\left(\mathrm{\Omega }\right)=H^k\left(\mathrm{\Omega }\right),H^0\left(\mathrm{\Omega }\right)=L^2\left(\mathrm{\Omega }\right).$ Further, we introduce the scalar product in $L^2\left(\mathrm{\Omega }\right)$ by

$$\left(f,g\right)={\int }_{\mathrm{\Omega }}f\left(x\right)g\left(x\right)\mathit{dx}.$$

The norm in $H^k\left(\mathrm{\Omega }\right)$ will be denoted by $\parallel \cdot {\parallel }_{H^k\left(\mathrm{\Omega }\right)}.$ The symbol ${\left|\cdot \right|}_{H^k\left(\mathrm{\Omega }\right)}$ will stand for the seminorm,

$${\left|u\right|}_{H^k\left(\mathrm{\Omega }\right)}={\left(\sum _{\mid \alpha \mid =k}\parallel D^{\alpha }u{\parallel }_{L^2\left(\mathrm{\Omega }\right)}^2\right)}^{\frac{1}{2}},\parallel u{\parallel }_{H^k\left(\mathrm{\Omega }\right)}={\left(\sum _{j=0}^k{\left|u\right|}_{H^j\left(\mathrm{\Omega }\right)}^2\right)}^{\frac{1}{2}}.$$

Let $V$ be a real Hilbert space with a scalar product ${\left(\cdot ,\cdot \right)}_V$ and an associated norm $\parallel \cdot {\parallel }_V,$$V^{\ast }$ its dual space. $K$ denotes a nonempty closed convex subset of $V,a\left(\cdot ,\cdot \right):V\times V\to \mathbf{R}$ is a bilinear, symmetric, continuous and elliptic form of $V$, $a\left(\cdot ,\cdot \right):V\times V\to \mathbf{R}$ is a bilinear, symmetric, continuous and elliptic form of $V\times V$; that is, there exist constants $\alpha >0$, and $\beta >0$ such that$a\left(u,v\right)\le \alpha \parallel u{\parallel }_V\parallel v{\parallel }_V,\forall u,v\in V$ and $a\left(v,v\right)\ge \beta \parallel v{\parallel }_V^2,\forall v\in V$. The pairing between $V$ and $V^{\ast }$ is denoted by $<\cdot ,\cdot >$. Let $\mathrm{\Lambda }$ be a canonical isomorphism from $V^{\ast }$ onto $V$ defined, for $g\in V^{\ast }$, by $<g,v>={\left(\mathrm{\Lambda }g,v\right)}_V,\forall v\in V.$ We can easily see that $\parallel \mathrm{\Lambda }{\parallel }_{V^{\ast }}=\parallel {\mathrm{\Lambda }}^{-1}{\parallel }_V=1$. Now, let us consider the following variational inequality:

$$\text{Find}u\in K\text{such}\text{that}a\left(u,v-u\right)\ge <f\left(u\right),v-u>,\forall v\in K,$$E1

where $f$ is a nonlinear operator such that $f\left(u\right)\in V^{\ast }.$.

In order to obtain a fixed point formulation of variational inequality (1) we need the following standard result.

Lemma 1.Let$K$be a closed convex subset of$V$. Then$u=P_K\omega$, the projection of$\omega$on$K$, if and only if

$$u\in K:{\left(u-\omega ,v-u\right)}_V\ge 0,\forall v\in K.$$E2

For some constant $\rho >0$, let us define a mapping $G:V\to V$ by

$$G\left(u\right)=P_K\mathrm{\Lambda \Phi }\left(u\right),$$E3

where $u\in V$, $\mathrm{\Phi }\left(u\right)\in V^{\ast }$ is defined by

$$<\mathrm{\Phi }\left(u\right),v>={\left(u,v\right)}_V-\mathit{\rho a}\left(u,v\right)+\rho <f\left(u\right),v>,\forall v\in V.$$E4

For some constant $\rho >0$, problem (1) can be written as

$${\left(u,v-u\right)}_V-\left\{{\left(u,v-u\right)}_V-\mathit{\rho a}\left(u,v-u\right)+\rho <f\left(u\right),v-u>\right\}\ge 0,\forall v\in K.$$

Using (4) in the above inequality, problem (1) is equivalent to that of finding $u\in K$ such that

$${\left(u-\mathrm{\Lambda \Phi }\left(u\right),v-u\right)}_V\ge 0,\forall v\in K.$$E5

By (2) and (5), we now have the following fixed point problem for the operator $G$:

$$u=P_K\mathrm{\Lambda \Phi }\left(u\right)=G\left(u\right).$$E6

Under appropriate conditions on the space $V$ and the operator $G:V\to V$(e.g., continuity, compactness), which usually have to be verified by theoretical means, fixed point theorem yields the existence of a solution $u$ of the problem (1) in some suitable set $U\subset V$, provided that

$$G\left(U\right)\subset U.$$E7

In order to compute an explicit inclusion, we must therefore construct $U$ explicitly. For the numerical verification of condition (7), we have to use interval analysis on many levels between basic interval arithmetic and functional analysis. For the appropriate and suitable choice of the operator $f$, the form $a\left(\cdot ,\cdot \right)$, and the convex set $K$; one encounters problems governed by the elliptic variational inequality as special cases from the problem (1) [48, 49, 50, 51, 52]. Inbrief, it is clear that the problem (1) is the most common. Up to now, devising a verification technique for the problem (1) is still an open problem. It is an important and interesting area of future research to find the numerical inclusion methods for the problem (1) by using (6). In this paper, we suppose that $V\subset L^2\left(\mathrm{\Omega }\right)$ and the nonlinear map $f\left(\cdot \right):V\to L^2\left(\mathrm{\Omega }\right)$ satisfies the following assumptions.

A1.$f$ is a continuous map from $V$ to $L^2\left(\mathrm{\Omega }\right).$

A2. For each bounded subset $W\subset V,$$f\left(W\right)$ is also bounded in $L^2\left(\mathrm{\Omega }\right).$

If we restrict the nonlinear map $f$ as above, then it can be shown that the problem (1) can be characterized by a class of variational inequality of the type,

$$\text{find}u\in K\text{such}\text{that}a\left(u,v-u\right)\ge \left(f\left(u\right),v-u\right),\forall v\in K.$$E8

The problem (8) has the restricted condition; even so (8) is an important and very useful class of nonlinear problems arising in mathematical physics, mechanics, engineering sciences, etc. In Section 3, we briefly consider a particular example of interest in applications. Another example is given in [13, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. In the special case in which $K\equiv V$, (8) yields the variational theory of the boundary value problems for partial differential equations. We will discuss existence and inclusion methods for problem (8). These are methods providing the existence of a solution of the problem (8) within explicitly computable bounds. As we have seen before, the transformation of problem (8) into some fixed point formulation (6) can be carried out in the same way. In a conclusion problem (8) is equivalent to the fixed point problem of finding $u\in K$ such that

$$u=S\left(u\right),$$E9

where $S$ denotes a specific operator, not necessarily the same as in (6). In particular for a given problem, we reduced the problem (8) to the fixed point formulation (9) and the continuity and compactness of $S$ is discussed. For this reason, we shall say nothing about this problem for which we refer to [13, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. In order to simplify argument we assume that $S$ is a continuous and compact operator. Since $S$ is continuous and compact, as a result of Schauder’s fixed point theorem, if there exists a nonempty, bounded, convex, and closed subset $U$ such that $S\left(U\right)\subset U$, then there exists a solution of $u=S\left(u\right)$ in $U$. In Sections 2.2 and 2.3, we describe how to construct $U$ explicitly.

2.2 Verification by a simple iteration method

In this subsection, we describe a simple iteration method for numerical verification of solutions for elliptic variational inequalities. In order to treat functions and variational inequalities in the infinite dimensional space $V$ by computer, we introduce two concepts, rounding and rounding error. Now, let $V_h$ be a finite dimensional subspace of $V$ dependent on $h\left(0<h<1\right)$ and let $K_h$ be a nonempty closed convex subset of $V_h$. Usually, $V_h$ is taken to be a finite element subspace with mesh size $h$. For the sake of simplicity, we shall define $K_h$, an approximate subset of $K$, by $K_h=V_h\cap K.$$K_h$ is a closed convex subset of $V_h$. In practical applications, the construction of $K_h$ is one of the difficulties presented by variational inequalities. For a given problem, several approximations are available. For a general study of the approximation of convex sets, we refer the reader to the work of Mosco [51]. We define the projection $P_{K_h}$ from $V$ into $K_h$ [49, 50]. That is, $v_h=P_{K_h}\left(u\right),$ the projection of $u$ into $K_h$, is defined as follows:

$$u=S\left(u\right),v_h\in K_h:{\left(v_h,\zeta -v_h\right)}_V\ge {\left(u,\zeta -v_h\right)}_V,\forall \zeta \in K_h.$$E10

To verify the existence of a solution of (9), we determine a set $W$ for a bounded, convex, and closed subset $U\subset V$ as

$$W=\left\{v\in V:v=S\left(u\right),u\in U\right\}.$$

From Schauder’s fixed point theorem, if $W\subset U$ holds, then there exists a solution of (8) in the set $U$. Our goal is to find a set $U$ which includes $W.$ For any subset $W\subset V$, we define $R\left(W\right)\subset K_h$ by the projection of $V$ to $K_h$, which is called the rounding of $W$. Additionally, we define RE(W), the rounding error of $W$, as a subset of $V$ so that $W\subset R\left(W\right)+\mathit{RE}\left(W\right)$ holds. Using $R\left(W\right)+\mathit{RE}\left(W\right)$ instead of $W$, the verification condition becomes

$$R\left(W\right)+\mathit{RE}\left(W\right)\subset U.$$E11

Let us describe the procedure more concretely. First, we consider the auxiliary problem: given $g\in L^2\left(\mathrm{\Omega }\right)$,

$$\text{find}u\in K\text{such}\text{that}a\left(u,v-u\right)\ge \left(g,v-u\right),\forall v\in K.$$E12

We note that, by well known result [49], there is a unique element $u$ which satisfies (12).

Secondly, we define the approximate problem corresponding to (12) as

$$a\left(u_h,v_h-u_h\right)\ge \left(g,v_h-u_h\right),\forall v_h\in K_h,u_h\in K_h$$E13

and (13) admit one and only one solution [49]. Error estimates for the variational inequalities can be found in [48, 49, 52], etc. Now, using (10), (12), (13) and error estimates, we make the following assumption.

A3. For each $u\in V$, there exists a positive constant $C$, independent of $u$ and $h$, such that

$$\parallel u-P_{K_h}u{\parallel }_V\le \mathit{Ch}\parallel g{\parallel }_{L^2\left(\mathrm{\Omega }\right)}.$$E14

In order to verify the solutions numerically, it is necessary to determine the constant $C$ that appears in a priori error estimations; this constant will be discussed later.

In order to construct the set $U$ satisfying the verification condition (11) in a computer, we use an iterative procedure, that is, the sequential iteration. We propose a computer algorithm to obtain the set $U$ which satisfies the condition (11).

(1) First, we obtain an approximate solution $v_h^{\left(0\right)}\in K_h$ to (8) by an appropriate method. Set $U_h^{\left(0\right)}=\left\{v_h^{\left(0\right)}\right\}$ and ${\alpha }_0=0.$.

(2) Next we will define $R\left(W^{\left(i\right)}\right)$ and $\mathit{RE}\left(W^{\left(i\right)}\right)$ for $i\ge 0,$ where $W^{\left(i\right)}$ is the set defined as follows:

$$W^{\left(i\right)}=\left\{v^{\left(i\right)}\in V:v^{\left(i\right)}=S\left(u^{\left(i\right)}\right),u^{\left(i\right)}\in U^{\left(i\right)}\right\}.$$

$R\left(W^{\left(i\right)}\right)$ is defined by the subset of $K_h$ which consists of all the elements $v_h^{\left(i\right)}\in K_h$ such that

$$a\left(v_h^{\left(i\right)},\psi -v_h^{\left(i\right)}\right)\ge \left(f\left(u^{\left(i\right)}\right),\psi -v_h^{\left(i\right)}\right),\forall \psi \in K_h,$$E15

holds for some $u^{\left(i\right)}\in U^{\left(i\right)}.$ Note that $R\left(W^{\left(i\right)}\right)$ can be enclosed by $R\left(W^{\left(i\right)}\right)\subset {\sum }_{j=1}^M{A}_j{\varphi }_j,$ where $A_j=\left[\underset{¯}{A_j},\overline{A_j}\right]$ are intervals, ${\left\{{\varphi }_j\right\}}_{j=1}^M$ is a basis of $V_h$, and $M=dim{V}_h$. For details of the interval calculation, we refer the reader to Nakao [6, 7, 12]. Next $\mathit{RE}\left(W^{\left(i\right)}\right)$ is defined as

$$\mathit{RE}\left(W^{\left(i\right)}\right)=\left\{v\in V:\parallel v{\parallel }_V\le \mathit{Ch}\underset{u^{\left(i\right)}\in U^{\left(i\right)}}{sup}\parallel f\left(u^{\left(i\right)}\right){\parallel }_{L^2\left(\mathrm{\Omega }\right)}\right\}.$$E16

Here, $C$ is the same constant as in (14). Hence, $W^{\left(i\right)}\subset R\left(W^{\left(i\right)}\right)+\mathit{RE}\left(W^{\left(i\right)}\right)$ holds.

(3) Check the verification condition:

$$R\left(W^{\left(i\right)}\right)+\mathit{RE}\left(W^{\left(i\right)}\right)\subset U^{\left(i\right)}.$$E17

If the condition is satisfied, then $U^{\left(i\right)}$ is the desired set, and a solution to (8) exists in $W^{\left(i\right)}$, and hence in $U^{\left(i\right)}.$

(4) If the condition is not satisfied, we continue the simple iteration by using $\delta -\mathit{\text{inflation}}$; that is, let $\delta$ be a certain positive constant given beforehand, and take

$$\begin{array}{l}{\alpha }_{i+1}=\mathit{Ch}{sup}_{u^{\left(i\right)}\in U^{\left(i\right)}}\parallel f\left(u^{\left(i\right)}\right){\parallel }_{L^2\left(\mathrm{\Omega }\right)}+\delta ,\\ \left[{\alpha }_{i+1}\right]=\left\{v\in V:\parallel v{\parallel }_V\le {\alpha }_{i+1}\right\},\\ U_h^{\left(i+1\right)}={\sum }_{j=1}^M\left[\underset{¯}{A_j}-\delta ,\overline{A_j}+\delta \right]{\varphi }_j,\\ U^{\left(i+1\right)}=U_h^{\left(i+1\right)}+\left[{\alpha }_{i+1}\right],\end{array}$$

and then go back to the second step. The reader may refer to [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46] for the details. If the condition (17) is satisfied, in our inclusion method of solutions for (9), the solution $u$ is enclosed in the set $U^{\left(i\right)}$, which we call ‘a candidate set’ of the form $U^{\left(i\right)}=U_h^{\left(i\right)}+\left[{\alpha }_i\right]$.

2.3 Verification by a Newton-like method

The significance of a Newton-like operator was already pointed out in [29, 43]. Hence we will not discuss it in detail here. In Subsection 2.1, numerical verification of solutions for elliptic variational inequalities using a finite element method have been discussed only for simple iteration method. The method proposed in Subsection 2.2 is such that $\left\{\left(U_h^{\left(i\right)},{\alpha }_i\right)\right\}$ always converges to the limit value $\left\{\left(U_h,\alpha \right)\right\}$ from an arbitrary initial value $\left\{\left(U_h^{\left(0\right)},{\alpha }_0\right)\right\}$ if $S$ in (9) is retractive operator (we refer to Zeidler [59, 60, 61] for the definition of retraction), while no convergence can generally be expected if $S$ is not retractive operator. Briefly, for not retractive operator in the neighborhood of the solution, it is difficult to use the previous scheme proposed in Subsection 2.2. To overcome such a difficulty, in this section, we newly formulate a verification method using the Newton-like method. This approach enables us to remove the restriction in Subsection 2.2 to the retraction property of the operator in the neighborhood of the solution. Namely, this technique can be applied to general variational inequalities without any retraction property of the associated operator $S$. We refer to [29, 43] for a detailed study of the properties of the Newton-like Method.

In this subsection, we use the notation of Section 2.2. We assume that $K_h=V_h\cap K$ is a closed convex cone with vertex at 0 and $K_h^{\ast }$ its dual. We note that $K_h^{\ast }$ is also a closed convex cone with vertex at $0$, which is the only point common to $K_h$ and $K_h^{\ast }$. From (10) it follows that $K_h^{\ast }$ is the set of points whose projections into $K_h$ is $0$. We need some additional lemma.

Lemma 2.Any$u\in V$can be uniquely decomposed into the sum of two orthogonal elements. That is,

$$u=P_{K_h}u\oplus \left(I-P_{K_h}\right)u=P_{K_h}u\oplus P_{K_h^{\ast }}u.$$

Here,$\oplus$denotes the sum of two orthogonal elements in the sense of$V.$.

Note that (9) can be rewritten as the following decomposed form in $K_h$ and $K_h^{\ast }$:

$$\left\{\begin{array}{l}{P}_{K_h}u=P_{K_h}S\left(u\right),\\ \left(I-P_{K_h}\right)u=\left(I-P_{K_h}\right)S\left(u\right).\end{array}\right.$$E18

In order to formulate a Newton-like verification condition for (18), we need a Fréchet derivative of the operator $S$. For most of the variational inequalities, the $S$ in (9) is not Fréchet differentiable at all. Therefore, in order to use a Newton-like type method, a major difficulty in numerically solving the fixed point formulation $u=S\left(u\right)$ is the treatment of the non-differentiable operator $S$. We need a suitable modification of the Fréchet derivative of $S$. Using some techniques, we can devise the approximate Fréchet derivative of $S$. Hence we shall assume that $\tilde{D}S\left(u\right)$ is the approximate Fréchet derivative of the $S\left(u\right)$ at $u$ as the linear operator. Let $\tilde{D}S\left(u\right)$ be designated as the Fréchet-like derivative of $S$ at $u$.

To consider the Newton-like operator for (18), we define the nonlinear operator $N_h:V\to V_h$ as

$$N_h\left(u\right)\equiv P_{K_h}u-{\left[I-\tilde{D}S\left(u_h\right)\right]}_h^{-1}\left(P_{K_h}-P_{K_h}S\right)\left(u\right).$$

Here $I$ is the identity operator and ${\left[I-\tilde{D}S\left(u_h\right)\right]}_h^{-1}$ denotes the inverse on $V_h$ of the restriction operator ${\left.\left[I-\tilde{D}S\left(u_h\right)\right]\right|}_{V_h}.$ Note that we will verify the existence of the inverse operator ${\left[I-\tilde{D}S\left(u_h\right)\right]}_h^{-1}$ from the nonsingularity of the matrix corresponding to ${\left.\left[I-\tilde{D}S\left(u_h\right)\right]\right|}_{V_h}$ in actual calculations.

Next we define the operator $T:V\to V$ as follows:

$$T\left(u\right)\equiv N_h\left(u\right)+\left(I-P_{K_h}\right)S\left(u\right).$$E19

Then $T$ is considered as the Newton-like operator for the former part of (18), but as the simple iterative operator for the latter part. $T$ becomes a compact and continuous map on $V$ by properties of $S$. Using some techniques, for a given problem we can not only define the Newton-like operator, but also devise a Newton-like Method. Furthermore, we obtain the following proposition and theorem.

Proposition 3.Given the assumption that$N_h\left(u\right)\in K_h,$

$$u=S\left(u\right)\iff u=T\left(u\right).$$E20

Theorem 4.If there exists a nonempty, bounded, convex, and closed subset$U\subset K$such that$T\left(U\right)=\left\{T\left(u\right)|u\in U\right\}\subset U$, then by the Schauder fixed point theorem, there exists a solution$u\in U$of$u=S\left(u\right)$.

When we decompose the set $U$ as $U=U_h\oplus U_{\perp }$ in Theorem 8.1, where $U_h\subset K_h$ and $U_{\perp }\subset K_h^{\ast }$, the verification condition can be written by

$$\left\{\begin{array}{l}{N}_h\left(U\right)\subset U_h,\\ \left(I-P_{K_h}\right)S\left(U\right)\subset U_{\perp }.\end{array}\right.$$E21

Here, $U_h$ is represented as the linear combination of the base functions of $V_h$ with interval coefficients, whereas $U_{\perp }$ is the intersection of $K_h^{\ast }$ with a ball in $V$. That is,

$$\begin{array}{c}{U}_h=\left\{{\phi }_h\in K_h:{\phi }_h=\sum _{j=1}^M{A}_j{\varphi }_j\text{with}{a}_j\in \left[\underset{¯}{A_j},\overline{A_j}\right]\right\},\\ U_{\perp }=\left\{\phi \in K_h^{\ast }:\parallel \phi {\parallel }_V\le \alpha \right\},\end{array}$$

respectively.

Note that $N_h\left(U\right)$ can be directly computed from $U_h$ and $U_{\perp }$ with additional information on the a priori error estimates. On the other hand, $\left(I-P_{K_h}\right)S\left(U\right)$ is evaluated using (14), by the following constructive error estimates for the finite approximate solution of variational inequality (8):

$$\parallel \left(I-P_{K_h}\right)S\left(U\right){\parallel }_V\le \mathit{Ch}\underset{u\in U}{sup}\parallel f\left(u\right){\parallel }_{L^2\left(\mathrm{\Omega }\right)}.$$

Therefore, the former condition in (21) is validated as the inclusion relations of corresponding coefficient intervals; the latter part can be checked by comparing two nonnegative real numbers.

Next we show a computer algorithm to construct the set $U$ which satisfies the verification condition (21). In order to realize it, we use the iteration method described in Subsection 2.2. Similarly to that in Subsection 2.2, we now generate the following iteration sequence $\left\{\left(U_h^{\left(n\right)},{\alpha }_n\right)\right\}$ for $n=\mathrm{0,1,2},\cdots .$ For $n\ge 1$, the $\delta$-inflation of $\left(U_h^{\left(n-1\right)},{\alpha }_{n-1}\right)$ is denoted by $\left({\tilde{U}}_h^{\left(n-1\right)},{\tilde{\alpha }}_{n-1}\right)$. Next, for the set ${\tilde{U}}^{\left(n-1\right)}={\tilde{U}}_h^{\left(n-1\right)}\oplus \left[{\tilde{\alpha }}_{n-1}\right]$, define $\left(U_h^{\left(n\right)},{\alpha }_n\right)$ by

$$\left\{\begin{array}{l}{U}_h^{\left(n\right)}\supset N_h\left({\tilde{U}}^{\left(n-1\right)}\right),\\ {\alpha }_n=\mathit{Ch}\underset{u\in {\tilde{U}}^{\left(n-1\right)}}{sup}\parallel f\left(u\right){\parallel }_{L^2\left(\mathrm{\Omega }\right)}.\end{array}\right.$$E22

Finally, the verification condition in a computer is given by the following theorem. The proof of Theorem 4 will be given here for the sake of completeness; it is based on Proposition 3 and Schauder’s fixed point theorem.

Theorem 5.For an integer$N$, if two relationships

$$U_h^{\left(N\right)}\subset {\tilde{U}}_h^{\left(N-1\right)}\mathit{\text{and}}{\alpha }_N<{\tilde{\alpha }}_{N-1}$$E23

hold, then there exists a solution $u$ of (8) in $U_h^{\left(N\right)}\oplus \left[{\alpha }_N\right]$. Here, the first term of (21) means the strict inclusion in the sense of each coefficient interval of $U_h^{\left(N\right)}$ and ${\tilde{U}}_h^{\left(N-1\right)}$.

3. Applications

The study for the numerical verification method for elliptic variational inequalities has been still made less progress than for the differential equation case. The author’s method in the present chapter can be also applied, in principal, to the verification of solutions of the practical problems. Namely, in Section 3.1, we first give, a slightly detailed descroption of the basic principle and formulation of our numerical verification method for the solution of obstacle problems with a homogeneous condition. This should be an appropriate introduction to another applications of our idea. The basic approach of the method consisits of the fixed point formulation of the problems and construction of the function set, in a computer, satisying the validation condition of a certain infinite dimensional fixed point theorem. We also mention that it is possible to extend the method to more general problems with non-homogenerous obstacles. Moreover, in order to apply our method to the problem whose associated operator is not retractive in a neighborhood of the solution, a Newton-like method is introduced. Next, in Section 3.2, we apply our method to another type of free boundary problem with appears in the elasto-plastic deformation theory. This problem causes some properties of non-smoothness in tha associated finite dimensional equations. But, we can also overcome such a difficulty by appling the solution method for non-smooth problems developed by [29, 32, 33]. In the Section 3.3, we briefly remark that our enclosure method can also be applied to the so-called simplified Signorini problem which is a simplified version of a problem accurring in the elasticity theory [43]. Finally, in Section 3.4, we show the way to apply our approach to elliptic variational inequalities of the second kind appearing in the flow problems of a viscos-plastic fulied in a pipe.

3.1 Obstacle problems

We introduce the verification method for solutions of the obstacle problem which is known as a free boundary problem to cahracterize the contacted zone by an obstacle $\psi$ in an elastic membrane region.

3.1.1 Homogeneous case

Here, ‘homogeneous’ stands for the case that obstacle $\psi \equiv 0$ in the whole domain.

3.1.1.1 Basic formulation of verification

Though the basic idea of verification is given in other places [26, 27, 28], in order to keep the paper as self-contained as possible, we describe rather detailed formulation and verification procedure for the present case.

Let $\mathrm{\Omega }$ be a bounded convex domain in ${\mathbb{R}}^n$, $1\le n\le 2$, with piecewise smooth boundary $\mathrm{\partial \Omega }$. We set ${\left.V\equiv H_0^1\left(\mathrm{\Omega }\right)=\{v\in H^1\left(\mathrm{\Omega }\right):v\right|}_{\mathrm{\partial \Omega }}=0\}$ and

$$a\left(u,v\right)=\left(\nabla u,\nabla v\right)$$

which is adopted as the inner product on $V$, where $\left(\cdot ,\cdot \right)$ stands for the inner product on $L^2\left(\mathrm{\Omega }\right)$. We define $K≔\left\{v\in V:v\ge 0a.e.\mathit{on}\mathrm{\Omega }\right\}.$.

First, we note that, by well-known result [49], for any $g\in L^2\left(\mathrm{\Omega }\right)$, the problem:

$$a\left(u,v-u\right)\ge \left(g,v-u\right),\forall v\in K,u\in K,$$E24

has a unique solution $u\in V\cap H^2\left(\mathrm{\Omega }\right)$, and the estimate

$${\left|u\right|}_{H^2\left(\mathrm{\Omega }\right)}\le \parallel g{\parallel }_{L^2\left(\mathrm{\Omega }\right)}$$E25

holds [49], where ${\left|w\right|}_{H^2}$ implies the semi-norm of $w$ in $H^2\left(\mathrm{\Omega }\right)$ defined by

$${\left|w\right|}_{H^2\left(\mathrm{\Omega }\right)}^2\equiv \sum _{i,j=1}^n\parallel \frac{{\partial }^{2}w}{\partial x_i\partial x_j}{\parallel }_{L^2\left(\mathrm{\Omega }\right)}^2.$$

Now consider the following elliptic variational inequalities with nonlinear right-hand side;

$$\left\{\begin{array}{l}\text{Find}w\in K\text{such}\text{that}\\ a\left(w,v-w\right)\ge \left(f\left(w\right),v-w\right),\forall v\in K.\end{array}\right.$$E26

We take an appropriate finite dimensional subspace $V_h$ of $V$ for $0<h<1$. Usually, $V_h$ is taken to be a finite element subspace with mesh size $h$. We then define $K_h$, an approximation of $K$, by

$$K_h=V_h\cap K=\left\{v_h|v_h\in V_h,v_h\ge 0\mathrm{on}\overline{\mathrm{\Omega }}\right\}.$$

We also define the projection $P_K$ from $V$ onto $K$. That is, $v=P_K\left(w\right),$ the projection of $w\in V$ into $K$, is defined as the unique solution of the following problem:

$$v\in K:a\left(v,\zeta -v\right)\ge a\left(w,\zeta -v\right),\forall \zeta \in K.$$E27

And define the projection $P_{K_h}$ from $V$ onto $K_h$. That is, $v_h=P_{K_h}\left(w\right),$ the projection of $w$ into $K_h$, is defined as follows:

$$v_h\in K_h:a\left(v_h,\zeta -v_h\right)\ge a\left(w,\zeta -v_h\right),\forall \zeta \in K_h.$$E28

Now, as one of the approximation properties of $K_h$, assume that.

For each $w\in K\cap H^2\left(\mathrm{\Omega }\right)$, there exists a positive constant $C_1$, independent of $h$, such that

$$\parallel w-P_{K_h}w{\parallel }_V\le C_{1}h{\left|w\right|}_{H^2\left(\mathrm{\Omega }\right)}.$$E29

Here, $C_1$ has to be numerically determined. For example, it is known that we may take $C_1=\frac{\sqrt{5}}{\pi }$ for the linear element in one dimensional case [27]. Furthermore, it will be readily seen that the same constant can be taken for the two dimensional bilinear element from the consideration on the proof of Theorem 5.1 in [27]. To verify the existence of a solution of (26) in a computer, we use the fixed point formulation.

First, note that, for each $w\in V$, there exists a unique $F\left(w\right)\in V$ such that

$$\left(\nabla F\left(w\right),\nabla v\right)=\left(f\left(w\right),v\right),\forall v\in V,$$E30

which also implies that

$$\left\{\begin{array}{l}-\mathrm{\Delta }F\left(w\right)=f\left(w\right)\text{in}\mathrm{\Omega },\\ F\left(w\right)=0\mathrm{on}\mathrm{\partial \Omega }.\end{array}\right.$$E31

Then the map $F:V\to V$ is compact. By (30), the problem (26) is equivalent to finding $w\in V$ such that

$$a\left(w,v-w\right)\ge a\left(F\left(w\right),v-w\right),\forall v\in K.$$E32

Using the definition (27) and (32), we now have the following fixed point problem for the compact operator $P_{K}F$.

$$\text{Find}\exists w\in V\text{such}\text{that}w=P_{K}F\left(w\right).$$E33

3.1.1.2 Verification condition

We introduce two concepts, rounding and rounding error, which enable us to deal with the infinite dimensional problem by finite procedures, that is, in a computer.

Now we define the dual cone of $K_h$ by

$$K_h^{\ast }=\left\{w\in V:a\left(w,v\right)\le 0,\forall v\in K_h\right\},$$

and note that $K_h^{\ast }$ is also closed convex cone in $V$ with vertex at $0$ which is the only point common to $K_h$ and $K_h^{\ast }$. From (28) it follows that $K_h^{\ast }$ is the set of points whose projections into $K_h$ is $0$.

Lemma 6.Any$w\in V$can be uniquely decomposed into the sum of two orthogonal elements. That is,

$$w=P_{K_h}w\oplus \left(I-P_{K_h}\right)w=P_{K_h}w\oplus P_{K_h^{\ast }}w.$$

Here,$\oplus$denotes the sum of two orthogonal elements in the sense of$V$.

For any $w\in V$, we now define the rounding $R\left(P_{K}F\left(w\right)\right)\in K_h$ by the solution of the following problem:

$$a\left(R\left(P_{K}F\left(w\right)\right),v_h-R\left(P_{K}F\left(w\right)\right)\right)\ge \left(f\left(w\right),v_h-R\left(P_{K}F\left(w\right)\right)\right),\forall v_h\in K_h.$$

Next, for any subset $W\subset V$, we define the rounding $R\left(P_K\mathit{FW}\right)\subset K_h$ by

$$R\left(P_K\mathit{FW}\right)=\left\{w_h\in K_h:w_h=R\left(P_{K}F\left(w\right)\right),w\in W\right\}.$$

Usually, $R\left(P_K\mathit{FW}\right)$ is enclosed and represented as a linear conbination of the base functions in $V_h$ with interval coefficients.

Moreover, for $W\subset V$, we define $\mathit{RE}\left(P_K\mathit{FW}\right)$, the rounding error of $P_K\mathit{FW}$, as a subset of $K_h^{\ast }$, that is,

$$\mathit{RE}\left(P_K\mathit{FW}\right)=\left\{v\in K_h^{\ast }:\parallel v{\parallel }_V\le C_{0}h\parallel f\left(W\right){\parallel }_{L^2}\right\},$$E34

where

$$\parallel f\left(W\right){\parallel }_{L^2}\equiv \underset{w\in W}{sup}\parallel f\left(w\right){\parallel }_{L^2}.$$

Here, $C_0\equiv C_1{C}_2$, where $C_1$ is the same positive constant as in (29), and $C_2$ is determined by the following regularity estimate for the solution to (24) of the form

$${\left|u\right|}_{H^2}\le C_2{‖g‖}_{L^2}.$$E35

Thus we may take as $C_2=1$ for the present case from (25). Then, we have

$$P_{K}F\left(w\right)-R\left(P_{K}F\left(w\right)\right)\in \mathit{RE}\left(P_{K}F\left(w\right)\right),\forall w\in W.$$

Therefore, the following verification condition is obtained by Schauder’s fixed point theorem.

Lemma 7.If there exists a nonempty, bounded, convex, and closed subset$W\subset K$such that

$$R\left(P_K\mathit{FW}\right)\oplus \mathit{RE}\left(P_K\mathit{FW}\right)\subset W,$$E36

then there exists a solution of $w=P_{K}F\left(w\right)$ in $W$.

We sometimes refer the above set $W$ as a candidate set, which we generate in computer so that it satisfies the condition (36).

3.1.1.3 Verification procedures

We describe the method to find a set $W$ satisfying (36) in th ebelow.

Consider the following approximate solution $w_h\in K_h$ of (24):

$$a\left(w_h,v_h-w_h\right)\ge \left(g,v_h-w_h\right),\forall v_h\in K_h,w_h\in K_h.$$E37

Since the bilinear form $a\left(\cdot ,\cdot \right)$ is symmetric, (37) is reduced to the quadratic programming problem:

$$\underset{v\in K_h}{min}\left[\frac{1}{2}a\left(v,v\right)-(g,v)\right..$$E38

Let ${\left\{{\varphi }_j\right\}}_{j=1\cdots M}$ be a basis of $V_h$ with usual linear functions such that ${\varphi }_j\left(x\right)\ge 0,\forall x\in \mathrm{\Omega }$ and satisfying

$${\varphi }_j\left(x_i\right)=\left\{\begin{array}{ll}1,& i=j,\\ 0,& i\ne j,\end{array}\right.$$

where $x_i$ is an interior node of the finite element mesh. Then (38) reduces to the following vector form:

$$\underset{w\ge 0}{min}\left[\frac{1}{2}{w}^{\prime }\mathit{Dw}-P^{\prime }w\right],$$E39

where $w\ge 0$ means the componentwise relation. Here, $D≔{\left(d_{\mathit{ij}}\right)}_{1\le i,j\le M}$ with $d_{\mathit{ij}}=\left(\nabla {\varphi }_i,\nabla {\varphi }_j\right),$ and $w$ is the coefficient vector with $\left\{{\varphi }_j\right\}$ of the function $v$ in (38). Also, $P≔{\left(\left(g,{\varphi }_j\right)\right)}_{1\le j\le M}$.

Furthermore, we define for any $\alpha \in R^{+}$, nonnegative real number, we set

$$\left[\alpha \right]\equiv \left\{\varphi \in K_h^{\ast };\parallel \varphi {\parallel }_V\le \alpha \right\}.$$

Then, for a given candidate set $W=W_h\oplus \left[\alpha \right]$ with $W_h\subset K_h$, the computation of the rounding $R\left(P_K\mathit{FW}\right)$ reduces to enclose an interval vector $Z=\left(Z_j\right)$ and $Y=\left(Y_j\right)$ satisfying the following nonlinear system of equations [27]:

$$\left\{\begin{array}{l}Y-\mathit{DZ}=-\left(f\left(W\right),{\varphi }_j\right),1\le j\le M,\\ Y_j{Z}_j=0,1\le j\le M.\end{array}\right.$$E40

Here, $\left(f\left(W\right),{\varphi }_j\right)$ is evaluated as an interval $B_j$ such that $\left\{\left(f\left(w\right),{\varphi }_j\right)|w\in W\right\}\subset B_j.$ In order to solve (40) with guaranteed accuracy, we use some interval approaches for nonlinear system of equations [19, 20]. Thus, using the solution of (40), we can enclose the set $R\left(P_K\mathit{FW}\right)$ in (36). Combining this with (34), we can successfully compute the left-hand side of (36) for any candidate set $W=W_h\oplus \left[\alpha \right]$.

Thus we can present a computational verification condition. In the actual computation, we use an iterative procedure with $\delta$-$\mathit{\text{inflation}}$ technique to find the set $W$ satisfying (36). Several numerical examples for verification are presented in [27] for one dimensional problem using linear finite element.

3.1.2 Non-homogeneous case

In this subsection, we consider the two-dimensional case. In order to verify solutions numerically, it is necessary to determine some constants that appear in the a priori error estimates. For the non-homogeneous case, we define $K≔\left\{v\in V:v\ge \psi \mathrm{a}.\mathrm{e}.\mathrm{on}\mathrm{\Omega }\right\}$, where $\psi$ is a given $H^2\left(\mathrm{\Omega }\right)$ function such that $\psi \le 0$ on $\mathrm{\partial \Omega }$ and is not identically equal to 0. Let $\mathrm{\Omega }$ be a square with side 1 and let ${\mathcal{T}}_h$ be the uniform triangulation of $\mathrm{\Omega }$. We introduce ${\mathrm{\Sigma }}_h=\left\{p;p\in \overline{\mathrm{\Omega }},p\text{is}\mathrm{a}\text{vertex}\text{of}T\in {\mathcal{T}}_h\right\}$ and define the approximate $V_h$ of $H_0^1\left(\mathrm{\Omega }\right)$ by ${\left.V_h=\{v_h;v_h\in H_0^1\left(\mathrm{\Omega }\right)\cap C^0\left(\overline{\mathrm{\Omega }}\right),v_h\right|}_T\in P_1,\forall T\in {\mathcal{T}}_h\}.$ Here, ${\left.v_h\right|}_T$ denotes the restriction of $v_h$ to $T$ and $P_1$ representing the space of polynomials in two variables of degree $\le 1$. It is then quite natural to approximate $K$ by

$$K_h=\left\{v_h\in V_h;v_h\left(p\right)\ge \psi \left(p\right),\forall p\in {\mathrm{\Sigma }}_h\right\}.$$

Note that, in general, $K_h\ne V_h\cap K$. Then, $P_K$ and $P_{K_h}$ are similarly defined as before, and we also have the constructive error estimates of the form, $\forall v_h\in K_h$ and $\forall v\in K$,

$$\parallel u_h-u{\parallel }_{H_0^1\left(\mathrm{\Omega }\right)}\le C\left(g,\psi ,h\right),$$E41

where,

$$C\left(g,\psi ,h\right)\le \underset{g\in L^2\left(\mathrm{\Omega }\right)}{sup}\sqrt{{\left(0.494\right)}^2{h}^2{\left|u\right|}_{H^2}^2+2\left(\parallel g{\parallel }_{L^2}+\parallel \mathit{Au}{\parallel }_{L^2}\right)\left({\left(0.494\right)}^2{h}^2{\left|u\right|}_{H^2}+6h^2{\left|\psi \right|}_{H^2}\right)}.$$

We provide a numerical example of verification in the two-dimensional case according to the procedures described in the previous section. Let $\mathrm{\Omega }=\left(0,1\right)\times \left(0,1\right)$. We consider the case $f\left(u\right)=\mathit{Ku}+sin\mathit{\pi x}sin2\mathit{\pi y}$ and $\psi =sin\mathit{\pi x}sin\mathit{\pi y}$. For simplicity, we only consider the uniform mesh here. First, we divide the domain into small triangles with a uniform mesh size $h$ and choose the basis of $V_h$ as the pyramid functions.

The execution conditions are as follows (Figures 1–3):

Figure 1.

Obstacle function $\psi$.

Figure 2.

Approximate solution $u_h^{\left(0\right)}$.

Figure 3.

Illustration of the contact zone.

$$\begin{array}{l}K=0.1,dim{V}_h=10\\ \text{Obstacle}\text{function}\psi =sin\mathit{\text{π x}}sin\mathit{\text{π y}}\\ \text{the}\text{outline}\text{of}\psi \text{is}\text{shown}\text{in}\text{Figure}1.\\ \text{Initial}\text{value}:u_h^{\left(0\right)}=\text{Galerkin}\text{approximation},{\alpha }_0=0\\ \text{the}\text{outline}\text{of}{u}_h^{\left(0\right)}\text{is}\text{shown}\text{in}\text{Figure}2.\\ \text{Illustration}\text{of}\text{contact}\text{zone}\text{between}\text{obstacle}\\ \text{and}\text{approximate}\text{solution}\text{is}\text{shown}\text{in}\text{Figure}3.\\ \text{Extension}\text{parameters}:\delta ={10}^{-5}.\end{array}$$

Results are as follows:

$$\begin{array}{l}\text{Iteration}\text{numbers}\text{for}\text{verification}:2\\ H_0^1\left(\mathrm{\Omega }\right)-\text{error}\text{bound}:0.15437\\ \text{Maximum}\text{width}\text{of}\text{coefficient}\text{intervals}\text{in}\left\{A_j^{\left(N\right)}\right\}=0.00001.\end{array}$$

Detailed arguments and with numerical examples are presented in [42].

3.1.3 A Newton-type verification method

The idea of the enclosure method for solutions of obstacle problems is based upon simply sequential iterations for the original fixed point operator $P_{K}F$. Therefore, it is difficult to apply the method to the problem of which associated operator is not retractive in a neighborhood of the solution. In order to overcome such a difficulty, we introduce an another formulation using a Newton-like operator. The essential point is the way to devise the Newton-like operator for a kind of non-differentiable map which defines the original problem.

To formulate a Newton-type verification condition, we need a Fréchet derivative of the operator $P_{K}F$. However, $P_{K}F$ is not Fréchet differentiable at all. Therefore, we define the approximate Fréchet-like derivative ${\tilde{D}}_{K}F\left(u_h\right)$ on $V_h$ for some $u_h\in K_h$ instead of the Fréchet derivative. Assume that ${\left\{{\varphi }_j\right\}}_{j=1\cdots M}$ is a basis of $V_h$, where $M=dim{V}_h$, such that ${\varphi }_j\left(x\right)\ge 0$ on $\mathrm{\Omega }$ and satisfying

$${\varphi }_j\left(x_i\right)=\left\{\begin{array}{l}1,i=j,\\ 0,i\ne j,\end{array}\right.$$

where $x_i$ is an interior node of the finite element mesh.

And, for $v_h\in V_h$, we represent it such as

$$v_h=\sum _{j=1}^M{v}_{\mathit{hj}}{\varphi }_j.$$

Here, ${\left(v_{\mathit{hj}}\right)}_{j=1,\cdots ,M}$ is called as the coefficient vector of $v_h$. Now we take a fixed subset $N_0\subset \left\{1,2,\cdots ,M\right\}$, define $V_{h,N_0}$, the closed subspace of $V_h$, by

$$V_{h,N_0}=\left\{v_h|v_h\in V_h,v_{\mathit{hj}}=0\mathit{\text{for}}j\notin N_0\right\}.$$

And let $P_{h,N_0}$ be a $H_0^1$-projection from $V$ onto $V_{h,N_0}$ defined by

$$a\left(u-P_{h,N_0}u,v\right)=0,\forall v\in V_{h,N_0},P_{h,N_0}u\in V_{h,N_0}.$$

In order to define ${\tilde{D}}_{K}F\left(u_h\right):V_h\to V_{h,N_0}$, we differentiate the first equation of (40) in $W$ at $W=u_h$ to get, for arbitrary $\delta \in V_h$,

$$\partial Y^{\ast }-D\partial Z^{\ast }=-{\left\{\left(f^{\text{'}}\left(u_h\right)\delta ,{\varphi }_j\right)\right\}}_{1\le j\le M}.$$E42

Here, $\partial Y^{\ast }={\left({\tilde{Y}}_j^{\ast }\right)}_{1\le j\le M}$ and $\partial Z^{\ast }={\left({\tilde{Z}}_j^{\ast }\right)}_{1\le j\le M}$, where ${\tilde{Y}}_j^{\ast }=0$ for $j\in N_0$ and ${\tilde{Z}}_j^{\ast }=0$ for $j\notin N_0$, respectively.

Then we define the approximate Fréchet-like derivative of $P_{K}F\left(u\right)$ at $u=u_h$, as the linear map ${\tilde{D}}_{K}F\left(u_h\right):V_h\to V_{h,N_0}$ such that, for each $\delta \in V_h$,

$${\tilde{D}}_{K}F\left(u_h\right)\left(\delta \right)≔\sum _{j=1}^M{\widetilde{Z_j}}^{\ast }{\varphi }_j.$$

We now assume that.

A4. The restriction to $V_{h,N_0}$ of the operator $P_{h,N_0}\left[I-{\tilde{D}}_{K}F\left(u_h\right)\right]:V_h\to V_{h,N_0}$ has the inverse operator

$${\left[P_{h,N_0}-{\tilde{D}}_{K}F\left(u_h\right)\right]}_h^{-1}:V_{h,N_0}\to V_{h,N_0}.$$

Here, $I$ means the identity map on $V_h$.

By using the above approximate Fréchet-like derivative, we define the Newton-like operator $N_h:V\to V_h$ by

$$N_h\left(w\right)\equiv P_{K_h}w-{\left[P_{h,N_0}-{\tilde{D}}_{K}F\left(u_h\right)\right]}_h^{-1}{P}_{h,N_0}\left(P_{K_h}-P_{K_h}{P}_{K}F\right)\left(w\right)).$$

Next we define the operator $T:V\to V$ as follows:

$$T\left(w\right)\equiv N_h\left(w\right)+\left(I-P_{K_h}\right)P_{K}F\left(w\right).$$

Then $T$ becomes a compact map on $V$ and it follows the fixed point problem $w=P_K\mathit{Fw}$ is equivalent to $w=T\left(w\right)$. Detailed arguments and with numerical examples are presented in [35].

3.2 Elasto-plastic torsion problems

In this subsection, we consider an enclosure metnod of solutions for elasto-plastic torsion problems governed by an elliptic variational inequalities [25, 32, 33]. The nonlinear elasto-plastic torsion problem is defined as the same type elliptic variational inequalities as (26) with

$$K≔\left\{v\in H_0^1\left(\mathrm{\Omega }\right);|\nabla v|\le 1\mathrm{a}.\mathrm{e}.\mathrm{on}\mathrm{\Omega }\right\}.$$E43

As is well known [56, 58], two sub-domains ${\mathrm{\Omega }}_p$ and ${\mathrm{\Omega }}_e$ defined by

$${\mathrm{\Omega }}_p=\left\{x;x\in \mathrm{\Omega },|\nabla u|=1\right\},$$

and

$${\mathrm{\Omega }}_e=\mathrm{\Omega }\backslash {\mathrm{\Omega }}_p=\left\{x;x\in \mathrm{\Omega },|\nabla u|<1\right\}$$

correspond to the plastic and elastic regions, respectively. The elastic region ${\mathrm{\Omega }}_e$ and the plastic region ${\mathrm{\Omega }}_p$ are not known beforehand and should be determined, therefore $\partial {\mathrm{\Omega }}_e\cap \partial {\mathrm{\Omega }}_p$ is actually the free boundary of the problem (26). The problem (26) has been formulated as the problem of finding $u$ satisfying

$$\left\{\begin{array}{l}-\mathrm{\Delta }u=f\left(u\right)\text{in}{\mathrm{\Omega }}_e,\\ \mid \nabla u\mid =1\text{in}{\mathrm{\Omega }}_p,\\ u=0\mathrm{on}\mathrm{\partial \Omega }.\end{array}\right.$$E44

The finite dimensional convex subset $K_h$ is also defined similarly as before:

$$K_h≔V_h\cap K=\left\{v_h|v_h\in V_h,|\mid \nabla v_h|\le 1\mathrm{a}.\mathrm{e}.\mathrm{on}\mathrm{\Omega }\right\}.$$E45

In order to formulate the verification procedure, we need a verified computational method for solving the finite dimensional part (rounding) and a constructive estimates for infinite dimensional part (rounding error) as in the previous subsection.

Following [49, 56], we define the Lagrangian functional $\mathcal{L}$ associated with (1) by

$$\mathcal{L}\left(v,\mu \right)=\frac{1}{2}{\int }_{\mathrm{\Omega }}{\left|\nabla v\right|}^2\mathit{dx}-\left(g,v\right)+\frac{1}{2}{\int }_{\mathrm{\Omega }}\mu \left({\left|\nabla v\right|}^2-1\right)\mathit{dx}.$$

It follows, from [49, 56], that if $L$ has a saddle point $\left\{u,\lambda \right\}\in H_0^1\left(\mathrm{\Omega }\right)\times L_{+}^{\infty }\left(\mathrm{\Omega }\right),$ then $u$ is a solution of (1), where $L_{+}^{\infty }\left(\mathrm{\Omega }\right)=\left\{q\in L^{\infty }\left(\mathrm{\Omega }\right);q\ge 0a.e.\mathit{\text{in}}\mathrm{\Omega }\right\}.$ We use the Uzawa algorithm to solve (1). Thus we can claculate the rounding $R\left(P_{K}F\left(W\right)\right)$, for a candidate set $W$, by solving the following problem with guaranteed error bounds:

$$\left\{\begin{array}{l}\text{Find}\left\{u_h,{\lambda }_h\right\}\in K_h\times {\mathrm{\Lambda }}_h\text{such}\text{that}\\ {\lambda }_h=max\left[{\lambda }_h+\rho \left({\left|\nabla u_h\right|}^2-1\right),0\right]\text{with}\rho >0.\\ {\int }_{\mathrm{\Omega }}\left(1+{\lambda }_h\right)\nabla u_h\cdot \nabla v_h\mathit{dx}=\left(f\left(W\right),v_h\right),\forall v_h\in V_h,u_h\in V_h,\end{array}\right.$$E46

The problem (46) can be formulated as a system of nonlinear and nonsmooth (nondifferentiable) equations. A verification method for nonsmooth equations by a generalized Krawczyk operator is studied in [1, 55]. We briefly describe the method presented by [55] in the below.

We consider the following equivalent system of nonlinear(and nondifferentiable) equation to (46) for a fixed $w\in W$

$$H\left(x\right)=0.$$E47

Here, we assume that $H:R^n\to R^n$ is locally Lipschitz continuous. The equivalence means that $x^{\ast }$ solves (46) if and only if $x^{\ast }$ solves (47). The method is based on the mean value theorem for local Lipschitz functions of the form

$$H\left(x\right)-H\left(y\right)\in \mathit{co}\partial H\left(\left[x\right]\right)\left(x-y\right),\text{for}\mathrm{all}x,y\in \left[x\right],$$

where [x] stands for an interval vector, “co” denotes the convex hull, and $\partial H$ the generalized Jacobian in Clarke’s sense [57], which is also considered as a slope function, and

$$\mathit{co}\partial H\left(\left[x\right]\right)≔\mathit{co}\left\{V\in \partial H\left(x\right);x\in \left[x\right]\right\}.$$

Let $\left[L_{\left[x\right]}\right]$ be an interval matrix such that $\mathit{co}\partial H\left(\left[x\right]\right)\subseteq \left[L_{\left[x\right]}\right].$ Then for any $x,y\in \left[x\right]\subseteq R^n$ it holds that $H\left(x\right)-H\left(y\right)\in \left[L_{\left[x\right]}\right]\left(x-y\right)$.

Then an interval operator for nonsmooth equations is defined by

$$G\left(x,A,\left[x\right]\right)≔x-A^{-1}H\left(x\right)+\left(I-A^{-1}\left[L_{\left[x\right]}\right]\right)\left(\left[x\right]-x\right).$$E48

The mapping $G\left(x,A,\left[x\right]\right)$ is called a generalized Krawczyk operator. Therefore, the verification condition of solutions for (46) in $\left[x\right]$ is given by

$$G\left(x,A,\left[x\right]\right)\subseteq \left[x\right]\subset D.$$

Thus, we can compute the solution of (46) with guaranteed accuracy. That is, we can enclose the rounding $R\left(P_{K}F\left(U\right)\right)$. On the other hand, in order for the calculation of the rounding error $\mathit{RE}\left(P_{K}F\left(U\right)\right)$, the similar arguements can also be applied for one dimensional problem. Actually, we can prove that the same constant $C_0=\frac{\sqrt{5}}{\pi }$ is also valid for the present problem in one dimensinal case, which implies that we can give a verification procedure besed on the same principle as before [25, 32, 33]. In [33], we extended the approach to the numerical proof of existence of solutions for elasto-plastic torsion problems as well as gave a numerical example for one dimensional case. The verification method in [33] is based on the generalized Krawczyk operator for solving a system of nonsmooth (nondifferentiable) equations. In order to use the generalized Krawczyk operator, we need to calculate the Jacobian. In that case, we need some complicated techniques. However, in many cases, calculating the generalized Jacobian is very difficult. To overcome such difficulties, we proposed a numerical verification method without using the generalized Krawczyk operator. This method is attractive, since calculating the generalized Jacobian is not required in the computational performance. Furthermore, up to know, our verification methods are mainly based on the enclosure of solutions in the sense of $L^2$ or $H^1$ norms. We considered a numerical verification method with guaranteed $L^{\infty }$ error bounds for the solution of elasto-plastic torsion problem.

3.3 Simplified Signorini problems

A simplified Signorini problem is also given by the elliptic variational inequalities of the form (26) with

$$K≔\left\{v\in H_0^1\left(\mathrm{\Omega }\right);|v\ge 0\mathrm{on}\mathrm{\partial \Omega }\right\}$$E49

and

$$a\left(u,v\right)={\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{vdx}+{\int }_{\mathrm{\Omega }}\mathit{\text{uvdx}}.$$E50

where

$$\nabla u\cdot \nabla v=\frac{\partial u}{\partial x_1}\frac{\partial v}{\partial x_1}+\frac{\partial u}{\partial x_2}\frac{\partial v}{\partial x_2}.$$

As well known, the solution $u$ of this elliptic variational inequalities can be characterized as a solution of the following free boundary problem finding $u$ and two subsets ${\mathrm{\Gamma }}_0$ and ${\mathrm{\Gamma }}_{+}$ such that ${\mathrm{\Gamma }}_0\cup {\mathrm{\Gamma }}_{+}=\mathrm{\partial \Omega }$ and ${\mathrm{\Gamma }}_0\cap {\mathrm{\Gamma }}_{+}=\mathrm{Ø}$

$$\left\{\begin{array}{l}-\mathrm{\Delta }u+u=f\left(u\right)\text{in}\mathrm{\Omega },\\ u=0\mathrm{on}{\mathrm{\Gamma }}_0,\frac{\partial u}{\partial n}\ge 0\mathrm{on}{\mathrm{\Gamma }}_0,\\ u>0\mathrm{on}{\mathrm{\Gamma }}_{+},\frac{\partial u}{\partial n}=0\mathrm{on}{\mathrm{\Gamma }}_{+},\end{array}\right.$$E51

where $\frac{\partial }{\partial n}$ the outer normal derivative on $\mathrm{\partial \Omega }$. In the present case, the approximation subspace $K_h$ is taken as

$$K_h≔V_h\cap K=\left\{v_h|v_h\in V_h,v_h\ge 0\mathrm{on}\mathrm{\partial \Omega }\right\}.$$E52

For a candidate set $W$, the computation of rounding $R\left(P_{K}F\left(W\right)\right)$ is also reduced to the quadratic programming problem as in the Section 3.1 [56].

Since the constant $C_2$ in (25) is easily estimated as $C_2=1$, the standard approximation property of the interpolation by $K_h$ gives a constructive error estimates to compute the rounding error $\mathit{RE}\left(P_{K}F\left(W\right)\right)$. For a simplified Signorini problem [43], we constructed a computing algorithm which automatically encloses the solution within guaranteed error bounds. In particular, the method proposed in [43] enables us to verify the free boundary of a simplified Signorini problem, which has been impossible so far. Concerning the numerical verification of solutions for elliptic variational inequalities, we would like to mention that the inclusion method described in this article can be applied to the solution of the elliptic variational inequalities on large space domains.

3.4 Some other problems

In this subsection, we show that our idea of verification method can also be applied to the elliptic variational inequalities of the second kind.

Now, we define the functional $j\left(v\right)={\int }_{\mathrm{\Omega }}\mid \nabla v\mid \mathit{dx}.$ We consider the following problem of the flow of a viscous plastic fluid in a pipe:

$$\left\{\begin{array}{l}\text{Find}u\in H_0^1\left(\mathrm{\Omega }\right)\text{such}\text{that}\\ a\left(u,v-u\right)+j\left(v\right)-j\left(u\right)\ge \left(f\left(u\right),v-u\right),\forall v\in H_0^1\left(\mathrm{\Omega }\right).\end{array}\right.$$E53

As in the previous section, we consider the following auxiliary problem associated with (53) for a given $g\in L^2\left(\mathrm{\Omega }\right):$

$$a\left(u,v-u\right)+j\left(v\right)-j\left(u\right)\ge \left(g,v-u\right),\forall v\in H_0^1\left(\mathrm{\Omega }\right),u\in H_0^1\left(\mathrm{\Omega }\right).$$E54

By the well known result, we have the following lemma.

Lemma 8. There exists a unique solution $u\in H_0^1\left(\mathrm{\Omega }\right)\cap H^2\left(\mathrm{\Omega }\right)$ of (54) for any $g\in L^2$, such that

$$\parallel u{\parallel }_{H^2\left(\mathrm{\Omega }\right)}\le \widehat{C}\parallel g{\parallel }_{L^2\left(\mathrm{\Omega }\right)}.$$

When we denote the solution $u$ of (54) by $u=\mathit{Ag}$ and define the composite map $F$ on $H_0^1\left(\mathrm{\Omega }\right)$ by $F\left(u\right)\equiv \mathit{Af}\left(u\right)$, which is a little bit of different from the previously appeared symbol $F$ in Section 2, we have.

Theorem 9.$F$ is compact on $H_0^1\left(\mathrm{\Omega }\right)$ and the problem (53) is equivalent to the fixed point problem

$$u=F\left(u\right).$$

Proof. First, for a bounded subset $U\subset L^2\left(\mathrm{\Omega }\right)$, we show that $\mathit{AU}\subset H_0^1\left(\mathrm{\Omega }\right)$ is relatively compact. Secondly, prove that $A:L^2\left(\mathrm{\Omega }\right)\to H_0^1\left(\mathrm{\Omega }\right)$ is continuous. By Lemma 3, $\mathit{AU}\subset H^2\left(\mathrm{\Omega }\right)\cap H_0^1\left(\mathrm{\Omega }\right)$ and $\mathit{AU}$ is bounded in $H^2\left(\mathrm{\Omega }\right)$. Since $U$ is bounded in $L^2\left(\mathrm{\Omega }\right),$ by the Sobolev imbedding theorem, we have $\mathit{AU}$ is relatively compact in $H_0^1\left(\mathrm{\Omega }\right).$ Next, for arbitrary $f_1,f_2\in L^2\left(\mathrm{\Omega }\right)$, setting $u_1={\mathit{Af}}_1$ and $u_2={\mathit{Af}}_2$, by using (54), we obtain

$$\begin{array}{l}a\left(u_1,u_2-u_1\right)+j\left(u_2\right)-j\left(u_1\right)\ge \left(f_1,u_2-u_1\right),\\ a\left(u_2,u_1-u_2\right)+j\left(u_1\right)-j\left(u_2\right)\ge \left(f_2,u_1-u_2\right).\end{array}$$

With the above inequalities, we obtain $a\left(u_2-u_1,u_2-u_1\right)=-a\left(u_1,u_2-u_1\right)+a\left(u_2,u_2-u_1\right)\le j\left(u_2\right)-$${\mathcal{T}}_h$ Hence, by the Poincaré inequality, we have

$$\parallel u_2-u_1{\parallel }_{H_0^1\left(\mathrm{\Omega }\right)}^2\le \parallel f_2-f_1{\parallel }_{L^2\left(\mathrm{\Omega }\right)}\parallel u_2-u_1{\parallel }_{L^2\left(\mathrm{\Omega }\right)}\le \overline{C}\parallel f_2-f_1{\parallel }_{L^2}\parallel u_2-u_1{\parallel }_{H_0^1\left(\mathrm{\Omega }\right)}.$$

Therefore, we obtain

$$\parallel u_2-u_1{\parallel }_{H_0^1\left(\mathrm{\Omega }\right)}\le \overline{C}\parallel f_2-f_1{\parallel }_{L^2\left(\mathrm{\Omega }\right)}.$$

That is, $A$ is Lipschitz continuous as a map $L^2\left(\mathrm{\Omega }\right)\to H_0^1\left(\mathrm{\Omega }\right).$ Hence $A$ is compact. The latter half in the theorem is straightforward from the definition of $F$.

We now define the approximate problem corresponding to (54) as

$$a\left(u_h,v_h-u_h\right)+j\left(v_h\right)-j\left(u_h\right)\ge \left(g,v_h-u_h\right),\forall v_h\in V_h,u_h\in V_h.$$E55

In order to apply our verification method to enclose the solutions of (53), we need a guaranteed computation of the exact solution of the problem (55), a rounding procedure, as well as the constructive error estimates between the solution of (54) and (55), rounding error estimates.

A major difficulty in solving the problem (55) numerically is the processing of the nondifferentiable term $j\left(u\right)={\int }_{\mathrm{\Omega }}\mid \nabla u\mid \mathit{dx}.$ One approach is the method of Lagrange multiplier on that term, whose continuous version is as follows [56].

Let us define $\mathrm{\Lambda }=\left\{q|q\in L^2\left(\mathrm{\Omega }\right)\times L^2\left(\mathrm{\Omega }\right),|q\left(x\right)|\le 1a.e.x\in \mathrm{\Omega }\right\}$ with $\mid q\left(x\right)\mid =\sqrt{q_1{\left(x\right)}^2+q_2{\left(x\right)}^2}.$ Then the solution $u$ of (54) is equivalent to the existence of $q$ satisfying

$$\left\{\begin{array}{l}a\left(u,v\right)+{\int }_{\mathrm{\Omega }}q\cdot \nabla v=\left(g,v\right),\forall v\in H_0^1\left(\mathrm{\Omega }\right),u\in H_0^1\left(\mathrm{\Omega }\right),\\ q\cdot \nabla u=\mid \nabla u\mid a.e.,q\in \mathrm{\Lambda }.\end{array}\right.$$E56

Moreover, it is known that (56) is equivalent to the following problem:

$$\left\{\begin{array}{l}a\left(u,v\right)+{\int }_{\mathrm{\Omega }}q\cdot \nabla v=\left(g,v\right),\forall v\in H_0^1\left(\mathrm{\Omega }\right),u\in H_0^1\left(\mathrm{\Omega }\right),\\ q=\frac{q+\rho \nabla u}{sup\left(1,|q+\rho \nabla u|\right)}.\end{array}\right.$$E57

Here $\rho$ is a positive constant. Let ${\mathcal{T}}_h$ be a triangulation of $\mathrm{\Omega }$, and let define $L_h$ and ${\mathrm{\Lambda }}_h$ (approximation of $L^{\infty }\left(\mathrm{\Omega }\right)\times L^{\infty }\left(\mathrm{\Omega }\right)$ and $\mathrm{\Lambda }$, respectively) by

$$L_h=\left\{q_h|q_h=\sum _{\tau \in {\mathcal{T}}_h}{q}_{\tau }{\chi }_{\tau },q_{\tau }\in R^2\right\}\text{and}{\mathrm{\Lambda }}_h=\mathrm{\Lambda }\cap L_h,\text{respectively},$$

where ${\chi }_{\tau }$ is the characteristic function of $\tau$.

Then our first purpose, computing the rounding $\mathit{RF}\left(U\right)$, is to enclose the solution of the following approximation problem of (57):

$$\left\{\begin{array}{l}a\left(u_h,v_h\right)+{\int }_{\mathrm{\Omega }}{q}_h\cdot \nabla v_h=\left(g,v_h\right),\forall v_h\in V_h,u_h\in V_h,\\ q_h=\frac{q_h+\rho \nabla u_h}{sup\left(1,|q_h+\rho \nabla u_h|\right)}.\end{array}\right.$$E58

The Eq. (58) leads to a kind of finite dimensional, nonlinear but nondifferentiable problem. We use a slope function method proposed by Rump [18, 19, 20] to enclose the solutions of (58) with $g=f\left(W\right)$ for a candidate set $W$. On the other hand, the rounding error $\mathit{RE}\left(F\left(U\right)\right)$ can be computed by using the following constructive error estimates:

Theorem 10. Let $u$ and $u_h$ be solutions of (54) and (55), respectively. If $g\in L^2\left(\mathrm{\Omega }\right)$, then there exists a constant $C\left(h\right)$ such that

$$\parallel u_h-u{\parallel }_{H_0^1\left(\mathrm{\Omega }\right)}\le C\left(h\right)\parallel g{\parallel }_{L^2\left(\mathrm{\Omega }\right)}.$$

Here, we may take $C\left(h\right)=\frac{\sqrt{5}}{\pi }h$ for the linear element in one dimensional case, and $C$ is also numerically estimated such that $C\left(h\right)\approx O\left(h^{\frac{1}{2}}\right)$ for the two dimensional linear element. A proof of this theorem is described in Ryoo and Nakao [34]. Thus we can also implement the verification algorithm for the solution of (53) as in the previous section. For details on this subsection, please refer to Ref. [47].

4. Conclusions

We have surveyed numerical verification methods for differential equations, especially around partial differential equations, variational inequalities and the author’s works. But the period of this research is shorter than the history of the numerical methods for differential equations by computer and we can say it is still in the stage of case studies. Indeed, recently, this kind of studies have been referred little by little for practical applications in PDEs and variational inequalities but there are many open problems to be resolve. Therefore, we can make no safe prediction that these approaches will grow into really useful methods for various kinds of equations and variational inequalities in mathematical analysis. Also, since the program description of the verification algorithm is very complicated in general, there is another problem like software technology associated with assurance for the correctness of the verification program itself. Actually, some of the mathematician would not give credit the computer assisted proof in analysis as correct as they believe the theoretical proof, which might cause a kind of seriously emotional problem in the methodology of mathematical sciences. And there is another difficulty from the huge scale of numerical computations which often exceed the capacity of the concurrent computing facilities.

However, in the twenty-first century, the computing environment would make more and more rapid progress, which should be beyond conception in the present state. In any case, a realistic study for partial differential equations and variational inequalities should be the future subject of the numerical computations with guaranteed accuracy. The authors believe that numerical methods with guaranteed accuracy for differential equations and variational inequalities would highly improve the reliability in the numerical simulation of the complicated phenomena in both mathematical and engineering sciences.

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (No. 2017R1A2B4006092).

Mathematics Subject Classification (2000)

65N15; 65G20