An approach to the numerical verification of solutions for variational inequalities using Schauder fixed point theory

Abstract

In this paper, we describe a numerical method to verify the existence of solutions for a unilateral boundary value problems for second order equation governed by the variational inequalities. It is based on Nakao’s method by using finite element approximation and its explicit error estimates for the problem. Using the Riesz representation theory in Hilbert space, we first transform the iterative procedure of variational inequalities into a fixed point form. Then, using Schauder fixed point theory, we construct a high efficiency numerical verification method that through numerical computation generates a bounded, closed, convex set which includes the approximate solution. Finally, a numerical example is illustrated.

MSC: 65G20, 65G30, 65N15, 65N30.

1 Introduction

A numerical verification method to verify the existence of solutions for mathematical problems is a new approach in the field of existence theory of solutions for mathematical problems that appear in mathematical analysis. Numerical verification methods of solutions for differential equations have been the subject of extensive study in recent years and much progress has been made both mathematically and computationally (see [1]–[11]etc.). These methods are known as new numerical approaches for the problems where it is difficult to prove analytically the existence of solutions for differential equations. However, for some problems governed by the variational inequality, 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. The theory of variational inequalities has become a rich source of inspiration in both mathematical and engineering sciences. So, a high efficiency numerical method for variational inequalities is often beneficial to the relevant subject. It is the aim of this paper to attempt a numerical technique to verify the solutions for elliptic equations of the second order with boundary conditions in the form of inequalities, that is, we construct a computing algorithm which automatically encloses the solution with guaranteed error bounds. In the following section, we describe the elliptic equations of the second order with boundary conditions in the form of inequalities considered and the fixed point formulation to prove the existence of solutions. In Section 3, in order to treat the infinite dimensional operator by computer, we introduce two concepts, rounding and rounding error, and a computational verification condition. In Section 4, we construct a concrete computing algorithm for the verification by computer, which is an efficient computing algorithm from the viewpoint of interval arithmetic. In order to verify solutions numerically, it is necessary to calculate the explicit a priori error estimates for approximate problems. These constants play an important role in the numerical verification method. In Section 5, we determine these constants. Finally, a numerical example is presented. Many difficulties remain to be overcome in the construction of general techniques applicable to a broader range of problems. However, the author has 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. We hope to make progress in this direction in the future.

2 Problem and fixed point formulation

Let us first settle on a few notations. In what follows we shall make use of the Sobolev spaces $W^{k,p}(\mathrm{\Omega })$ 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}(\mathrm{\Omega })=H^k(\mathrm{\Omega })$, $H^0(\mathrm{\Omega })=L^2(\mathrm{\Omega })$. Further, we introduce the scalar product in $L^2(\mathrm{\Omega })$ by

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

Let $\mathrm{\Omega }\subset {\mathbb{R}}^2$ be a bounded domain with piecewise smooth boundary Γ. We consider the unilateral boundary value problem

$$-\mathrm{\Delta }u=g\text{in }\mathrm{\Omega },$$

(2.1)

with boundary conditions of two types:

$$\begin{array}{c}u=0\text{on }{\mathrm{\Gamma }}_0\subset \mathrm{\Gamma },\hfill \\ u\ge 0,\frac{\partial u}{\partial v}\ge 0\text{and}u\frac{\partial u}{\partial v}=0\text{on }{\mathrm{\Gamma }}_{+}=\mathrm{\Gamma }-{\mathrm{\Gamma }}_0.\hfill \end{array}$$

Let us always assume that the right-hand side of (2.1) fulfills $g\in L^2(\mathrm{\Omega })$. We define

$$a(u,v)={\int }_{\mathrm{\Omega }}\mathrm{\nabla }u\cdot \mathrm{\nabla }vdx,\mathrm{\nabla }u\cdot \mathrm{\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}.$$

Set $V=\{v\in H^1(\mathrm{\Omega }):v=0\text{ on }{\mathrm{\Gamma }}_0\}$, and denote the inner product and norm on V, respectively, as follows:

$${(u,v)}_V=(\mathrm{\nabla }u,\mathrm{\nabla }v),{\parallel u\parallel }_V={\parallel \mathrm{\nabla }u\parallel }_{L^2(\mathrm{\Omega })}={|u|}_{H^2(\mathrm{\Omega })}.$$

Problem (2.1) may be formulated as a variational problem. To do this, let us define the set

$$K=\{v\in H^1(\mathrm{\Omega }):u=0\text{ on }{\mathrm{\Gamma }}_0,u\ge 0\text{ on }{\mathrm{\Gamma }}_{+}\}$$

and the potential energy functional

$$J(v)=\frac{1}{2}a(v,v)-(g,v).$$

Then the functional J is continuous, strictly convex, and coercive in the space V. From these properties of J and results of optimization theory [12], it follows that the minimization problem is finding u such that

$$J(u)\le J(v),\mathrm{\forall }v\in K,u\in K.$$

(2.2)

Hence problem (2.2) is equivalent to the problem of finding u such that

$$(J^{\prime }(u),v-u)\ge 0,\mathrm{\forall }v\in K,u\in K,$$

(2.3)

where $J^{\prime }(u)$ is the Gâteaux derivative of J at u. Since $(J^{\prime }(u),v)=a(u,v)-(g,v)$ and $a(u,v)$ is symmetric, problem (2.2) is equivalent to that of finding $u\in K$ such that

$$a(u,v-u)\ge (g,v-u),\mathrm{\forall }v\in K.$$

(2.4)

Now, let us consider the following variational inequality:

$$\text{Find }u\in K\text{ such that }a(u,v-u)\ge (f(u),v-u),\mathrm{\forall }v\in K.$$

(2.5)

Here, we suppose the following conditions for the map f.

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

A2. For each bounded subset $U\in V$, $f(U)$ is also a bounded set in $L^2(\mathrm{\Omega })$.

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

Lemma 2.1

[13]

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

$$u\in K\text{:}a(u,v-u)\ge a(\omega ,v-u),\mathrm{\forall }v\in K.$$

(2.6)

Then, for each $u\in V$, from the Riesz representation theorem, there exists a unique element $F(u)\in V$ such that

$$a(F(u),v)=(f(u),v),\mathrm{\forall }v\in V$$

(2.7)

and the map $F:V\to V$ is a compact operator (see [6]).

By (2.7), problem (2.5) is equivalent to that of finding $u\in V$ such that

$$a(u,v-u)\ge a(F(u),v-u),\mathrm{\forall }v\in K,u\in K.$$

(2.8)

By Lemma 2.1 and (2.8), we have the following fixed point problem for the compact operator $P_{K}F$:

$$\text{Find }u\in V\text{ such that }u=P_{K}F(u).$$

(2.9)

Then, using the fixed point problem (2.9), we can construct the numerical procedure to verify the existence of a solution for the variational inequality (2.5).

3 Rounding and verification conditions

In order to describe the numerical verification procedure, we introduce two concepts, rounding and rounding error. For the sake of simplicity, let us assume that $\mathrm{\Omega }\in {\mathbb{R}}^2$ is a bounded domain with a polygonal boundary Γ. We shall denote by $I_{\mathrm{\Omega }}=\{1,2,3,\dots ,m_0\}$ the set of all indices i associated with the internal nodes $x_i$ of the domain Ω and we shall denote by $I_{\mathrm{\Gamma }}=\{m_0+1,m_0+2,\dots ,m\}$ the set of all node indices i associated with the boundary nodes $x_i$ of the domain Ω and we let $I=I_{\mathrm{\Omega }}\cup I_{\mathrm{\Gamma }}$. Here, for the sake of simplicity, let us assume that $I_{{\mathrm{\Gamma }}_0}=\{m_0+1,m_0+2,\dots ,m_0+j\}$, $I_{{\mathrm{\Gamma }}_{+}}=\{m_0+j+1,\dots ,m\}$, and $I_{\mathrm{\Gamma }}=I_{{\mathrm{\Gamma }}_0}\cup I_{{\mathrm{\Gamma }}_{+}}$. In what follows, we shall consider only a regular system of triangulations. In other words, when refining the partition of $\overline{\mathrm{\Omega }}$, the triangles of the given triangulation do not reduce to segments. Let $\{{\mathcal{T}}_h\}$ be a regular system of triangulations of $\overline{\mathrm{\Omega }}$. The nodes of a triangulation lying on $I_{\mathrm{\Gamma }}$ will be denoted by $p_{m_0+1},p_{m_0+2},\dots ,p_m$. We then approximate V by

$$V_h=\{v_h\in C^0(\overline{\mathrm{\Omega }}):v_h{|}_T\in P_1(T),\mathrm{\forall }T\in {\mathcal{T}}_h\},$$

where $P_k(T)$ denotes the set of all polynomials of degree at most k on the definition domain T. We then define $K_h$, an approximate subset of K, by

$$K_h=\{v_h\in V_h:v_h(p_i)=0,\mathrm{\forall }i\in I_{{\mathrm{\Gamma }}_0},v_h(p_i)\ge 0,\mathrm{\forall }i\in I_{{\mathrm{\Gamma }}_{+}}\}.$$

It is easily seen that $K_h$ is a closed, convex, and nonempty subset of $V_h$.

We then define the approximate problem corresponding to (2.4) as

$$a(u_h,v_h-u_h)\ge (g,v_h-u_h),\mathrm{\forall }{v}_h\in K_h.$$

(3.1)

Let u be the solution of (2.4) and $u_h\in K_h$ be the approximate solution of (3.1). Now, as one of the approximation properties of $K_h$, assume the following.

A3. For each $w\in K$, there exists a positive constant $C(h)$ such that

$${\parallel w-P_{K_h}w\parallel }_{L^2(\mathrm{\Omega })}\le C(h){\parallel g\parallel }_{L^2(\mathrm{\Omega })}.$$

(3.2)

Here, $C(h)$ has to be numerically determined.

For any $u\in K$, we now define the rounding $\mathcal{R}(P_{K}F(u))\in K_h$ as the solution of the following variational inequality:

$$a(\mathcal{R}(P_{K}F(u)),v_h-\mathcal{R}(P_{K}F(u)))\ge (f(u),v_h-\mathcal{R}(P_{K}F(u))),\mathrm{\forall }{v}_h\in K_h.$$

For a set $U\subset V$, we define the rounding $\mathcal{R}(P_{K}FU)\subset K_h$ as

$$\mathcal{R}(P_{K}FU)=\{u_h\in K_h:u_h=\mathcal{R}(P_{K}F(u)),u\in U\}.$$

Also, we define for $U\subset V$ the rounding error $\mathcal{RE}(P_{K}FU)\subset V$ as

$$\mathcal{RE}(P_{K}FU)=\{v\in V:{\parallel v\parallel }_{L^2(\mathrm{\Omega })}\le C(h){\parallel f(U)\parallel }_{L^2(\mathrm{\Omega })}\},$$

(3.3)

where

$${\parallel f(U)\parallel }_{L^2(\mathrm{\Omega })}\equiv \underset{u\in U}{sup}{\parallel f(u)\parallel }_{L^2(\mathrm{\Omega })}.$$

The positive constant $C(h)$ appearing here is numerically determined in Section 5 by using the approximation property of $K_h$ expressed by

$$P_{K}F(u)-\mathcal{R}(P_{K}F(u))\in \mathcal{RE}(P_{K}F(u)),\mathrm{\forall }u\in U.$$

With the above, we have the following as a result of the Schauder fixed point theorem.

Theorem 3.1

If there exists a nonempty, bounded, convex, and closed subset$U\subset K$such that$\mathcal{R}(P_{K}FU)+\mathcal{RE}(P_{K}FU)\subset U$, then there exists a solution of$u=P_{K}F(u)$in U.

4 Computing procedures for verification

In this section, we propose a computer algorithm to obtain a set U which satisfies the condition of Theorem 3.1.

Now, we define the approximate problem corresponding to (2.4) as

$$a(u_h,v_h-u_h)\ge (g,v_h-u_h),\mathrm{\forall }{v}_h\in K_h,u_h\in K_h.$$

(4.1)

As parameters to describe a function $v_h\in V_h$ we choose the values $v_h(p_i)$ of $v_h$ at the nodes $p_i$, $i=1,\dots ,m$, of ${\mathcal{T}}_h$. The corresponding basis functions ${\varphi }_j\in V_h$, $j=1,\dots ,m$, are defined by ${\varphi }_j(p_i)={\delta }_{ij}$ (Kronecker’s symbol). A function $v_h\in V_h$ now has the representation

$$v_h(t)=\sum _{j=1}^m{z}_j{\varphi }_j(t),z_j=v_h(p_j)\text{ for }t\in \overline{\mathrm{\Omega }}.$$

By [14], (4.1) is actually equivalent to the following discrete system:

$$\{\begin{array}{l}{D}_{I_{\mathrm{\Omega }}I}{z}_I-P_{I_{\mathrm{\Omega }}}=0,\\ (D_{I_{\mathrm{\Gamma }}I}{z}_I-P_{I_{\mathrm{\Gamma }}})z_{I_{\mathrm{\Gamma }}}=0,\\ z_{I_{\mathrm{\Gamma }}}\ge 0,z_{I_{{\mathrm{\Gamma }}_0}}=0,\\ D_{I_{\mathrm{\Gamma }}I}{z}_I-P_{I_{\mathrm{\Gamma }}}\ge 0.\end{array}$$

(4.2)

Here, $D_{II}\equiv {(a_{ij})}_{i,j\in I}$, with $a_{ij}=(\mathrm{\nabla }{\varphi }_i,\mathrm{\nabla }{\varphi }_j)$ and $z_I$ is the coefficient vector for $\{{\varphi }_i\}$ corresponding to the function $u_h$ in (4.1). Further, $P_I\equiv {((g,{\varphi }_i))}_{i\in I}$ is an m dimensional vector.

Thus we can proceed in the following manner. Let ${\mathbb{R}}^{+}$ denote the set of all nonnegative real numbers. For $\alpha \in {\mathbb{R}}^{+}$ we associate

$$[\alpha ]\equiv \{\varphi \in V:{\parallel \varphi \parallel }_{L^2(\mathrm{\Omega })}\le \alpha \}.$$

(4.3)

Let $A_j$ ($1\le j\le m$) be intervals on ${\mathbb{R}}^1$ and let ${\sum }_{j=1}^m{A}_j{\varphi }_j$ be a linear combination of $\{{\varphi }_j\}$, i.e., an element of the power set $2^{V_h}$ in the following sense:

$$\sum _{j=1}^m{A}_j{\varphi }_j=\{\sum _{j=1}^m{a}_j{\varphi }_j:a_j\in A_j,1\le j\le m\}.$$

Let us denote all the sets of linear combinations of $\{{\varphi }_j\}$ with interval coefficients by

$${\mathbf{D}}_I\equiv \{\sum _{j=1}^m{A}_j{\varphi }_j:A_j;\text{interval in }\mathbb{R},1\le j\le m\}.$$

Then, setting $U={\sum }_{j=1}^m{A}_j{\varphi }_j+[\alpha ]$ and $g=f(U)$ in (4.1), we consider the following nonlinear system:

$$\{\begin{array}{l}{D}_{I_{\mathrm{\Omega }}I}{z}_I-P_{I_{\mathrm{\Omega }}}=0,\\ w_{I_{\mathrm{\Gamma }}}\equiv D_{I_{\mathrm{\Gamma }}I}{z}_I-P_{I_{\mathrm{\Gamma }}}\ge 0,\\ z_{I_{\mathrm{\Gamma }}}\ge 0,z_{I_{{\mathrm{\Gamma }}_0}}=0,\\ w_{I_{\mathrm{\Gamma }}}{z}_{I_{\mathrm{\Gamma }}}=0.\end{array}$$

(4.4)

Here $P_I\equiv {((f(U),{\varphi }_i))}_{i\in I}$.

System (4.4) is in fact a bilinear system of equations whose right-hand side consists of intervals with constraint conditions $z_{I_{\mathrm{\Gamma }}}\ge 0$ and $w_{I_{\mathrm{\Gamma }}}\ge 0$. To solve the nonlinear system (4.4) with automatic verification of the correctness of the result, a verification method for nonsmooth equations by a generalized Krawczyk operator as in [15] could be used. We adopt here another method. Setting $x=(z,w)\in {\mathbb{R}}^{2m-m_0}$, (4.4) without constraint is written as a nonlinear system of equations,

$$\mathcal{F}(x)=0.$$

(4.5)

Let $\tilde{x}:=({\tilde{z}}_1,{\tilde{z}}_2,\dots ,{\tilde{z}}_{m_0},{\tilde{z}}_{m_0+1},\dots ,{\tilde{z}}_m,{\tilde{w}}_{m_0+1},{\tilde{w}}_{m_0+2},\dots {\tilde{w}}_m)$ be an approximate solution of (4.4). Then note that ${\tilde{z}}_i\approx 0$ or ${\tilde{w}}_i\approx 0$ for each $m_0+1\le i\le m$.

Problem (4.5) can also be reformulated by nonsmooth equations using other methods, e.g., [15]. However, (4.5) is continuous and differentiable. Hence, to enclose solutions for (4.5), we use the following theorem proposed by [5].

Theorem 4.1

Let$\mathcal{F}:{\mathbb{R}}^{2m-m_0}\to {\mathbb{R}}^{2m-m_0}$be a function with continuous first derivative and let$M\in {\mathbb{R}}^{(2m-m_0)\times (2m-m_0)}$ (real$(2m-m_0)\times (2m-m_0)$matrix), $\tilde{x}\in {\mathbb{R}}^{2m-m_0}$. Denote the Jacobian matrix of ℱ by${\mathcal{F}}^{\prime }\in {\mathbb{R}}^{(2m-m_0)\times (2m-m_0)}$and for$X\in I{\mathbb{R}}^{2m-m_0}$ (real interval vectors with$2m-m_0$components) define${\mathcal{F}}^{\prime }(X):=\bigcap \{Y\in I{\mathbb{R}}^{m+m_0}:{\mathcal{F}}^{\prime }(x)\in Y\mathit{\text{ for all }}x\in X\}$. If for some$X\in I{\mathbb{R}}^{m+m_0}$with$0\in X$

$$-M\cdot \mathcal{F}(\tilde{x})+\{I-M\cdot {\mathcal{F}}^{\prime }(\tilde{x}+X)\}\cdot X\subseteq \stackrel{\circ }{X},$$

then there exists an$\hat{x}\in \tilde{x}+\stackrel{\circ }{X}$with$\mathcal{F}(\hat{x})=0$.

Let $X=(Z,W)$ be an enclosure of a solution of the nonlinear system (4.5) by using Theorem 4.1, where $Z:=(Z_1,Z_2,\dots ,Z_{m_0},Z_{m_0+1},\dots ,Z_m)\in I{\mathbb{R}}^m$ and $Y:=(W_{m_0+1},\dots ,W_m)\in I{\mathbb{R}}^{m-m_0}$. Then we set $Z_i:=0$ or $W_i:=0$ for each $m_0+1\le i\le m$ provided that ${\tilde{z}}_i\approx 0$ or ${\tilde{w}}_i\approx 0$, respectively. If, for all $i\in I_{\mathrm{\Gamma }}$, $\{Z_i=0\text{ and }inf(W_i)>0\}$ and $\{inf(Z_i)>0\text{ and }{W}_i:=0\}$ hold, then it implies that the problem (4.7) has an optimal solution $x\in X$ (cf.[5]). As one can see, for the case that ${\tilde{z}}_i$ and ${\tilde{w}}_i$ are both close to zero, this algorithm would not work. Fortunately, we have never encountered such a difficulty up to now. But, in order to establish more general applications of our method, it should be necessary to consider the methods for nonsmooth problems such as in [15].

We now consider the fully automatic computer generation of the set U satisfying Theorem 3.1. First, we generate a sequence of sets $\{U^{(i)}\}$, $i=0,1,\dots$ , which consists of subsets of V in the following manner.

We present an iterative procedure for generating ${\{U^{(i)}\}}_{i=0,\dots }$ (cf.[6], [7]). For $i=0$, we choose appropriate initial values $u_h^{(0)}\in K_h$ and ${\alpha }_0\in {\mathbb{R}}^{+}$, and define $U^{(0)}\subset V$ by

$$U^{(0)}=u_h^{(0)}+[{\alpha }_0].$$

Usually, $u_h^{(0)}$ is determined as

$$a(u_h^{(0)},v_h-u_h^{(0)})\ge (f\left(u_h^{(0)}\right),v_h-u_h^{(0)}),\mathrm{\forall }{v}_h\in K_h,u_h^{(0)}\in K_h.$$

(4.6)

This corresponds to the Galerkin approximation for (2.5). The standard selection for ${\alpha }_0$ will be ${\alpha }_0=0$. For $u_h^{(i)}={\sum }_{j=1}^m{A}_j^{(i)}{\varphi }_j$ and ${\alpha }_i\in {\mathbb{R}}^{+}$, we set $U^{(i)}=u_h^{(i)}+[{\alpha }_i]$, $i\ge 1$. Then we define $u_h^{(i+1)}\subset K_h$ and ${\alpha }_{i+1}\in {\mathbb{R}}^{+}$ according to

$$\{\begin{array}{l}{D}_{I_{\mathrm{\Omega }}I}{z}_I-P_{I_{\mathrm{\Omega }}}=0,\\ (D_{I_{\mathrm{\Gamma }}I}{z}_I-P_{I_{\mathrm{\Gamma }}})z_{I_{\mathrm{\Gamma }}}=0,\\ z_{I_{\mathrm{\Gamma }}}\ge 0,\\ D_{I_{\mathrm{\Gamma }}I}{z}_I-P_{I_{\mathrm{\Gamma }}}\ge 0,\\ P_I={((f(U^{(i)}),{\varphi }_i))}_{i\in I},\end{array}$$

(4.7)

$${\alpha }_{i+1}=C(h){\parallel f\left(U^{(i)}\right)\parallel }_{L^2(\mathrm{\Omega })},$$

(4.8)

where $C(h)$ is the same as in (3.3). Here, $u_h^{(i+1)}$ is determined as the solution set of (4.7), as described above. Of course, the solution of (4.7) satisfies the conditions of Theorem 4.1 in an application to the case in which $U=U^{(i)}$. By using (4.7) and (4.8), we define the map $F:{\mathbf{D}}_I\times {\mathbb{R}}^{+}\to {\mathbf{D}}_I\times {\mathbb{R}}^{+}$ by

$$(u_h^{(i)},{\alpha }_i)=F(u_h^{(i-1)},{\alpha }_{i-1})\text{for }i\ge 1,$$

(4.9)

and we can denote the above procedure as

$$U^{(i+1)}=\mathcal{R}\left(P_{K}F{U}^{(i)}\right)+\mathcal{RE}\left(P_{K}F{U}^{(i)}\right),i=0,1,\dots .$$

For $n\ge 1$, first for a given $0<\delta \ll 1$, we define the δ-inflation of $(u_h^{(n-1)},{\alpha }_{n-1})$ by

$$\{\begin{array}{l}{\tilde{u}}_h^{(n-1)}=u_h^{(n-1)}+{\sum }_{j=1}^m[-1,1]\delta {\varphi }_j,\\ {\tilde{\alpha }}_{n-1}={\alpha }_{n-1}+\delta .\end{array}$$

Next, for the set ${\tilde{U}}^{(n-1)}={\tilde{u}}_h^{(n-1)}+[{\tilde{\alpha }}_{n-1}]$, we compute $(u_h^{(n)},{\alpha }_n)$ by

$$(u_h^{(n)},{\alpha }_n)=F({\tilde{u}}^{(n-1)},{\tilde{\alpha }}_{n-1}).$$

(4.10)

Now we have the following verification condition on a computer.

Theorem 4.2

If for an integer N, the two relationships

$$u_h^{(N)}\subset {\tilde{u}}_h^{(N-1)}\mathit{\text{and}}{\alpha }_N<{\tilde{\alpha }}_{N-1}$$

(4.11)

hold, then there exists a solution u of (2.9) in$u_h^{(N)}+[{\alpha }_N]$. Here, the first term of (4.11) means the inclusion in the sense of each coefficient interval of$u_h^{(N)}$and${\tilde{u}}_h^{(N-1)}$.

For a convergence analysis of the iterative method for generating a sequence of set $\{U^{(i)}\}$, we will prove that the concerned sequence converges for the case that the nonlinear operator $P_{K}F$ in (2.9) is retractive around the solution u, and provided that the mesh size h is sufficiently small. We will leave such a general case as a further research topic.

5 Computation of the constants

In this section, we only deal with the one dimensional case. We give a bound of the constant $C(h)$ of (3.3).

Let $\mathrm{\Omega }=(a,b)$ and let $g\in L^2(\mathrm{\Omega })$. Then the basic model problem (2.5) is written as:

$$\text{Find }u\in K\text{ such that }a(u,v-u)\ge (g,v-u),\mathrm{\forall }v\in K.$$

(5.1)

We can represent the above problem (5.1) in the following form:

$$\begin{array}{c}Au=g\text{on }\mathrm{\Omega }\text{ with }Av=-v^{″},\hfill \\ u\ge 0\text{on }\{a,b\},\hfill \\ u(a)=0\text{or}u(b)=0,\hfill \\ u^{\prime }(a)\le 0,u^{\prime }(b)\ge 0,\hfill \\ u{u}^{\prime }=0\text{on }\{a,b\}.\hfill \end{array}$$

(5.2)

Let M be an integer >0 and let $h=\frac{1}{M}$. We consider $x_i=ih$ for $i=0,1,2,\dots ,M$ (that is, a uniform partition of Ω) and $e_i=(x_{i-1},x_i)$, $i=1,2,\dots ,M$. We then approximate $H^1(\mathrm{\Omega })$ by

$$V_h=\{v_h\in C^0(\mathrm{\Omega }):v_h{|}_{e_i}\in P_1,i=1,2,3,\dots ,M\}$$

with, as usual, $P_1$ representing the space of polynomials of degree ≤1, and we approximate K by

$$K_h=\{v_h\in S_h:v_h(a)=0\text{ or }{v}_h(b)=0,v_h\ge 0\text{ on }\{a,b\}\}.$$

The approximate problem is then defined by the following:

$$\text{Find }{u}_h\in K_h\text{ such that }a(u_h,v_h-u_h)\ge (g,v_h-u_h),\mathrm{\forall }{v}_h\in K_h.$$

(5.3)

We now consider the $L^2(\mathrm{\Omega })$ estimates of optimal order (that is, $O(h^2)$) of $u_h-u$ via a generalization of the Aubin-Nitsche method. The following result is given by arguments similar to those in [13], except for obvious modifications. Since the basic notations and results are also the same as that of Natterer [16], we do not discuss it further. The reader may refer to [13] for the details.

Regarding the approximation error ${\parallel u_h-u\parallel }_{L^2(\mathrm{\Omega })}$, we then have the following.

Theorem 5.1

Let u and$u_h$be solutions of problems (5.1) and (5.3), respectively. If$g\in L^2(\mathrm{\Omega })$, then we have

$${\parallel u_h-u\parallel }_{L^2(\mathrm{\Omega })}\le \frac{h^2}{{\pi }^2}{\parallel g\parallel }_{L^2(\mathrm{\Omega })}.$$

Hence, we may take$C(h)=\frac{h^2}{{\pi }^2}$in (3.3).

Proof

Following Natterer [16], we derive that

$$\begin{array}{c}\overline{T}(K,u_h)=\{v\in H^1(\mathrm{\Omega }):v\ge 0\text{ on }{B}_h\},\hfill \\ B_h=\{x\in \{a,b\}:u_h(x)=0\}.\hfill \end{array}$$

Let e be such that $Ae=g$ and $e^{\prime }(a)=e^{\prime }(b)=0$. Hence, for $w\in H^1(\mathrm{\Omega })$,

$$a(e-u,w)=(A(e-u),w)+{(e-u)}^{\prime }w{|}_a^b=-u^{\prime }w{|}_a^b.$$

Then, setting

$$G=\{w\in H^1(\mathrm{\Omega }):w\ge 0\text{ on }{B}_h,u^{\prime }(w+u_h-u){|}_a^b\le 0\},$$

there exists a $z\in G$ such that

$$a(z,w-z)\ge (u-u_h,w-z),\mathrm{\forall }w\in G.$$

(5.4)

Next, we consider $Az=u-u_h$. By using (5.4), we have

$$\begin{array}{rcl}a(z,w-z)& =& {\int }_a^b{z}^{\prime }{(w-z)}^{\prime }\\ \ge & {\int }_a^b(u-u_h)(w-z)\\ =& {\int }_a^b-z^{″}(w-z).\end{array}$$

Hence, we obtain

$$z^{\prime }(w-z){|}_a^b\ge 0,\mathrm{\forall }w\in G.$$

Also, for $z(b)>0$ we have $z^{\prime }(b)\le 0$, and similarly we obtain $z^{\prime }(a)\ge 0$ for $z(a)>0$.

We now have the estimate

$${\parallel z\parallel }_{L^2(\mathrm{\Omega })}^2\le a(z,z)=(Az,z)+z^{\prime }z{|}_a^b\le (u-u_h,z)\le {\parallel u-u_h\parallel }_{L^2(\mathrm{\Omega })}{\parallel z\parallel }_{L^2(\mathrm{\Omega })}.$$

Therefore we have

$${\parallel z\parallel }_{L^2(\mathrm{\Omega })}\le {\parallel u-u_h\parallel }_{L^2(\mathrm{\Omega })}\text{and}{\parallel z^{″}\parallel }_{L^2(\mathrm{\Omega })}\le {\parallel u-u_h\parallel }_{L^2(\mathrm{\Omega })}.$$

(5.5)

Next, for $z\in K$, we define the linear interpolation $r_{h}z$ by

$$r_{h}z\in S_h,(r_{h}z)(x_i)=z(x_i),i=0,1,\dots ,M.$$

Note that $r_{h}z\in K_h$. Therefore, by the standard results of approximation theory [9] and (5.5), we have

$$\begin{array}{r}{\parallel r_{h}z-z\parallel }_{L^2(\mathrm{\Omega })}\le \frac{h^2}{{\pi }^2}{\parallel z^{″}\parallel }_{L^2(\mathrm{\Omega })}\le \frac{h^2}{{\pi }^2}{\parallel u-u_h\parallel }_{L^2(\mathrm{\Omega })},\\ {\parallel r_{h}z-z\parallel }_V\le \frac{h}{\pi }{\parallel z^{″}\parallel }_{L^2(\mathrm{\Omega })}\le \frac{h}{\pi }{\parallel u-u_h\parallel }_{L^2(\mathrm{\Omega })}.\end{array}$$

(5.6)

Also, replacing v by $r_{h}u$ in [13], Theorem 1], we obtain

$$\begin{array}{rcl}{\parallel u_h-u\parallel }_V^2& \le & 2{\parallel Au-g\parallel }_{L^2(\mathrm{\Omega })}{\parallel u-r_{h}u\parallel }_{L^2(\mathrm{\Omega })}+{\parallel u-r_{h}u\parallel }_V^2\\ \le & \frac{h^2}{{\pi }^2}{\parallel g\parallel }_{L^2(\mathrm{\Omega })}^2.\end{array}$$

(5.7)

Hence, replacing y by $r_{h}z$, $y-z=0$ on $\{a,b\}$ in [13], Theorem 2], (5.6), and (5.7), we obtain

$$\begin{array}{rcl}{\parallel u_h-u\parallel }_{L^2(\mathrm{\Omega })}^2& \le & {\parallel Au-g\parallel }_{L^2(\mathrm{\Omega })}{\parallel z-r_{h}z\parallel }_{L^2(\mathrm{\Omega })}+{\parallel u_h-u\parallel }_V{\parallel r_{h}z-z\parallel }_V\\ \le & \frac{h^2}{{\pi }^2}{\parallel g\parallel }_{L^2(\mathrm{\Omega })}{\parallel u_h-u\parallel }_{L^2(\mathrm{\Omega })}.\end{array}$$

Therefore, we deduce that

$${\parallel u_h-u\parallel }_{L^2(\mathrm{\Omega })}\le \frac{h^2}{{\pi }^2}{\parallel g\parallel }_{L^2(\mathrm{\Omega })}.$$

 □

6 Example of numerical verification

In this section, we provide some numerical examples of verification in the one dimensional case according to the procedure described in the previous section.

Let $\mathrm{\Omega }=(0,1)$. We consider the case $f(u)=Ku+(\frac{4K-9{\pi }^2}{4})sin\frac{3\pi }{2}x$ and use a uniform partition of Ω, that is, $x_i=\frac{i}{M}$, $0\le i\le M$. Set $e_i=(x_{i-1},x_i)$; then we have $h=\frac{1}{M}$. We take

$$V_h\equiv \{v_h\in C^0(0,1):v_h{|}_{e_i}\in P_1(e_i),1\le i\le M\},$$

where $P_1(e_i)$ is the space of polynomials of degree ≤1 on $e_i$. We now choose the basis ${\{{\varphi }_i\}}_{i=1}^M$ of $V_h$ as the usual hat functions.

The execution conditions are as follows:

$\text{Numbers of elements}=300$;

$K=1/10$;

Extension parameters: $ϵ={10}^{-3}$;

Initial values: $u_h^{(0)}=\text{Galerkin approximation (4.6)}$; ${\alpha }_0=0$.

The form of $u_h^{(0)}$ is displayed in Figure 1.

Figure 1

Approximation solution $\mathit{y}\mathbf{=}{\mathit{u}}_{\mathit{h}}^{\mathbf{(}\mathbf{0}\mathbf{)}}$ .

The results are as follows:

Iteration numbers for verification: $N=3$;

$L^2$-error bound: 0.00002;

Maximum width of coefficient intervals in $\{A_j^{(N)}\}=0.0000454239$;

Coefficient intervals: as in Table 1.

Table 1 Coefficient intervals

The verification succeeded for h from $1/100$ to $1/300$. In Table 2, we show the values of α and $max|A_j^{(n)}|$, which is the maximum width of the coefficient intervals on the nodes.

Table 2 Maximum width of the coefficient intervals

Remark 6.1

In the above calculations, we carried out all numerical computations using the usual double precision computer arithmetic instead of strict interval computations (e.g., ACRITH-XSC, PASCAL-XSC, FORTRAN-XSC, C-XSC, PROFIL, etc.). Therefore, we neglected the round-off error. The reason is that the main purpose of our numerical experiments is the estimation of the truncation errors which usually, roughly speaking, are over 10⁻¹⁰ times larger than the round-off errors. That is, there will be in general some rounding errors at each step.