- Email: [email protected]
Numerical verification of solutions for obstacle problems
Journal of Computational and Applied Mathematics 161 (2003) 405 – 416 www.elsevier.com/locate/cam
Numerical veri$cation of solutions for obstacle problems Cheon Seoung Ryooa;∗ , Mitsuhiro T. Nakaob a
Department of Mathematics, Hannam University, Daejeon 306-791, South Korea b Faculty of Mathematics, Kyushu University 33, Fukuoka 812-8581, Japan Received 19 August 2002; received in revised form 6 May 2003
Abstract This article is an extension of the previous paper (Numer. Math. 81 (1998) 305) by the same authors. We propose a method to prove the existence of solutions for more general obstacle problems by numerical computations. The main task in this paper consists of the numerical determination of some constants which appear in a priori error estimations for the $nite element approximation of a simple variational inequality. We present a numerical example for veri$cation. c 2003 Elsevier B.V. All rights reserved. Keywords: Numerical veri$cation; Obstacle problems; Error estimates; Quadratic programming
1. Introduction In recent years, several numerical methods have been proposed to verify the solutions of the differential equations. In the case of partial di:erential equations, several techniques have been studied in [6–11]. However, for some problems governed by the variational inequality, there are very few approaches. In the previous paper [14], we described a method to verify the existence of solutions for simple obstacle problems with convex set K ={v ∈ H01 () : v ¿ 0} based upon $nite element approximations, explicit a priori error estimates and Schauder’s $xed point theorem. In the present paper, we extend the method to more generalized obstacle problems with convex set K = {v ∈ H01 () ; v ¿ } and is the height of the obstacle. The main improvement compared with the previous paper [14] is that we newly established the constructive a priori error estimates of the $nite element approximation for the nonhomogeneous, which means is not identically 0, obstacle problems. It is readily seen that the proof technique used in [14] can no longer be applied to the present case. Actually, we
This work was supported by Korea Research Foundation Grant, KRF-2001-003-D00011. Corresponding author. Tel.: +82-42-629-7525; fax: +82-42-629-7894. E-mail addresses: [email protected] (C.S. Ryoo), [email protected] (M.T. Nakao).
∗
c 2003 Elsevier B.V. All rights reserved. 0377-0427/$ - see front matter doi:10.1016/j.cam.2003.05.006
406
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
have used quite a di:erent approach to get the desired error estimations. Also note that, in general, it is not possible to transform the concerned problem to the case that ≡ 0, which is actually a di:erent situation from the case of the elliptic equations. Getting constructive error estimates plays an essential role in our veri$cation method. The remainder of the paper is organized as follows. In Section 2, we describe the basic concept of numerical veri$cation methods which have been developed by the authors. Using a $xed point formulation, we can verify the existence of solutions via Schauder’s $xed point theorem. More speci$cally, we can verify a set which includes a solution. Such a set is constituted by the sum of two subsets in H01 (). One is the computation of a projection into a closed convex subset of some $nite-dimensional subspace, the other is the estimation of the error for the projection. Thus, the values of constants appearing in error estimations play an essential role in the veri$cation. We describe a method to estimate such constants numerically in Section 3. In this setting we describe, in Section 4, procedures verifying the existence of solutions by numerical enclosing. Finally, a numerical example for a generalized obstacle problem is presented. 2. The problem and the methods of verication Let be a bounded convex domain in Rn , 1 6 n 6 2, with piecewise smooth boundary 9. We set H01 () = {v ∈ H 1 (); v|9 = 0} and a(u; v) = ∇u · ∇v d x;
where ∇u · ∇v = 99xu1 99xv1 + 99xu2 99xv2 . Let A be the operator de$ned by Au = −Mu. Then, for u; v ∈ H01 (), we have a(u; v) = (Au; v), where (·; ·) denotes the L2 -inner product on . Denote · L2 as usual the L2 norm on . Next, we de$ne K = {v ∈ H01 (); v ¿ a:e: on }, where, is a given function in H 2 () such that 6 0 on 9. First, we note that, by the well-known result [4], for any g ∈ L2 , the problem a(u; v − u) ¿ (g; v − u);
∀v ∈ K; u ∈ K
(2.1)
has a unique solution u ∈ H01 () ∩ H 2 (), and the estimate |u|H 2 () 6 max{ g L2 ; }
(2.2)
holds. Here is the solution of the equation = max(A ; 0) L2 ( + g L2 )=( − g L2 ) and |u|H 2 implies the semi-norm of u on H 2 () de$ned by 2 92 u 2 2 |u|H 2 () ≡ 9xi 9xj 2 : L i; j=1 We adopt (∇; ∇ ) as the inner product on H01 (). Hence, the associated norm is de$ned by
H01 () = ∇ L2 . Now, let us consider the following obstacle problem: Find u ∈ K such that a(u; v − u) ¿ (f(u); v − u);
∀v ∈ K:
(2.3)
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
407
Here, assume that f satis$es the following hypotheses: A1. f is a continuous map from H01 () to L2 (). A2. For each bounded subset W ⊂ H01 (), f(W ) is also a bounded set in L2 (). To verify the existence of a solution of (2.3) on a computer, we use a $xed point formulation. Since a(·; ·) is a continuous bilinear form on H01 () × H01 (), for each u ∈ H01 (), from the Riesz representation theorem, there exists a unique element F(u) ∈ H01 () such that a(F(u); v) = (f(u); v) for all v ∈ H01 (). That is, ∃ F(u) ∈ H01 () such that − MF(u) = f(u) in ; F(u) = 0 on 9: Then the map F: H01 () → H01 () is a compact operator by the above assumptions on f. As in the preceding paper [14], de$ning the projection PK : H01 () → K, problem (2.3) is equivalent to that of $nding u ∈ H01 () such that u = PK F(u):
(2.4)
Now we describe a numerical veri$cation method to verify the existence of solutions of (2.3). We determine a set V ≡ PK F(U ) for a bounded, convex and closed subset U ⊂ H01 () as V = {v ∈ H01 (); v = PK F(u); ∀u ∈ U }:
(2.5)
From Schauder’s $xed point theorem, if V ⊂ U holds, then there exists a solution of (2.3) in the set U . Our goal is to $nd a set U which includes V . A procedure to verify V ⊂ U using a computer is similar to that in [14] as below. We now take a $nite element subspace Sh of H01 () for 0 ¡ h ¡ 1 with an appropriate mesh. Letting N denote the set of nodes associated with the space Sh , we then de$ne Kh , an approximation of K, by Kh = {vh ∈ Sh ; vh (p) ¿ (p); ∀p ∈ N}: Note that usually Kh ⊂ K, namely, we use an outer approximation. It is well known [3] that, for an arbitrary solution u ∈ H01 () of (2.1) and its $nite element approximation uh ∈ Kh de$ned by a(uh ; vh − uh ) ¿ (g; vh − uh )
∀vh ∈ Kh ;
(2.6)
there exists a C(g; ) such that
u − uh H01 () 6 C(g; )h:
(2.7)
We describe the method to estimate C(g; ) numerically in Section 3. For any u ∈ H01 (), we de$ne the rounding R(PK F(u)) ∈ Kh as the solution of the following problem: a(R(PK F(u)); vh − R(PK F(u))) ¿ (f(u); vh − R(PK F(u)))
∀vh ∈ Kh :
We de$ne the rounding R(V ) ⊂ Kh for the set V de$ned by (2.5) as R(V ) = {vh ∈ Kh ; vh = R(PK F(u)); u ∈ U }:
408
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
Also, we de$ne the rounding error RE(V ) ⊂ H01 () as 1 RE(V ) = v ∈ H0 (); v H01 () 6 sup C(f(u); )h : u ∈U
(2.8)
From the de$nition, we have V ⊂ R(V ) + RE(V ). Then it is suNcient to $nd U satisfying R(V ) + RE(V ) ⊂ U . Next let us introduce the procedure for $nding such a set U , which is referred as a candidate set, using computers. In order to $nd a set U satisfying the above condition, we use a simple iterative method as in [14]. First, we obtain an approximate solution uh(0) ∈ Kh to (2.3) by some appropriate method. Set U (0) = {uh(0) } and 0 = 0. Then, we iterate the following procedure: For i ¿ 0, and a candidate set U (i) = Uh(i) + [ i ], where i ∈ R+ and [ i ] = {v ∈ H01 (); v H01 () 6 i }, set V (i) = {v(i) ∈ K; v(i) = PK F(u(i) ); u(i) ∈ U (i) }: In order to enclose V (i) , let us de$ne R(V (i) ) as follows. R(V (i) ) is de$ned by the subset of Kh which consists of all elements vh(i) ∈ Kh such that a(vh(i) ; ! − vh(i) ) ¿ (f(u(i) ); ! − vh(i) );
∀! ∈ Kh
holds for some u(i) ∈ U (i) . Note that R(V (i) ) can be enclosed by R(V (i) ) ⊂ nj=1 Aj j where Aj = [Aj ; Aj ] are intervals, and {}nj=1 is a basis of Sh . Secondly, RE(V (i) ) is de$ned by RE(V (i) ) =
v ∈ H01 (); v H01 () 6 sup C(f(u(i) ); )h : u(i) ∈U (i)
Then the veri$cation condition is written as V (i) ⊂ R(V (i) ) + RE(V (i) ) ⊂ U (i) ;
(2.9)
which is implied by R(V (i) ) ⊂ Uh(i) and RE(V (i) ) ⊂ [ i ]. If the veri$cation condition is satis$ed, then U (i) is the desired set, and a solution to (2.3) exists in V (i) , and hence in U (i) . When the veri$cation condition is not satis$ed, we continue the simple iteration by using "-in6ation, i.e., let " be a certain positive constant given beforehand, and take i+1
[
= sup C(f(u(i) ); )h + ";
i+1 ]
Uh(i+1)
u(i) ∈U (i)
= {v ∈ H01 (); v H01 () 6 =
n
[ Aj − ";
i+1 };
Aj + " ]j ;
j=1
U (i+1) = Uh(i+1) + [
i+1 ]:
Then we continue the above iteration process (see [14]).
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
409
3. Computation of the value C (g; ) In this section, we consider the one-dimensional case. In order to verify solutions numerically, it is necessary to determine the value C(g; ) that appears in (2.7). Let = (0; 1) and g ∈ L2 . Now, let n be a positive integer and let h = 1=n + 1. We de$ne xi := ih for i = 0; 1; 2; : : : ; n + 1 (that is, a uniform partition of ) and ei := (xi−1 ; xi ); i = 1; 2; : : : ; n + 1. We then approximate H01 () by Sh = {vh ∈ C 0 (); vh (0) = vh (1) = 0; vh |ei ∈ P1 ; i = 1; 2; 3; : : : ; n + 1} with, as usual, P1 representing the space of polynomials of degree 6 1, thus dim Sh = n, and Kh = {vh ∈ Sh ; vh (xi ) ¿ (xi ); i = 0; 1; 2; : : : ; n; n + 1}. Let {j }j=1:::n be a basis of Sh such that j (x) ¿ 0; ∀x ∈ , which satis$es 1; i = j; j (xi ) = 0; i = j: Regarding the approximation error uh − u H01 () , we then have Theorem 3.1. Let u and uh be the solutions of (2.1) and (2.6), respectively. If g ∈ L2 (), then we have
uh − u H01 () 6 C(g; )h; where
(3.1)
1 2 |u|H 2 () + 2(||g||L2 + |u|H 2 () )(|u|H 2 () + | |H 2 () ): $ is estimated by g L2 and using (2.2).
C(g; ) 6 Here, |u|H 2 ()
Proof. We have a(uh − u; uh − u) 6 a(u; u) + a(uh ; uh ) − a(u; uh ) − a(uh ; u) and, using (2.1) and (2.6), a(u; u) 6 a(u; v) + (g; u − v)
∀v ∈ K;
a(uh ; uh ) 6 a(uh ; vh ) + (g; uh − vh )
∀vh ∈ Kh :
Hence, we deduce that, for all v ∈ K and all vh ∈ Kh , a(uh − u; uh − u) 6 a(u; v − uh ) + a(uh ; vh − u) + (g; u − v) + (g; uh − vh ) = a(u; v − uh ) − (g; v − uh ) + a(u; vh − u) − (g; vh − u) + a(uh − u; vh − u) = (g − Au; uh − v) + (g − Au; u − vh ) + a(u − uh ; u − vh ): Since g − Au ∈ L2 (), we obtain, for all v ∈ K and vh ∈ Kh , that
u − uh 2H 1 () 6 g − Au L2 ( u − vh L2 + uh − v L2 ) 0
+ u − uh H01 () u − vh H01 () :
410
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
Since
u − uh H01 () u − vh H01 () 6 12 ( u − uh 2H 1 () + u − vh 2H 1 () ) 0
0
we have, by combining the two previous inequalities, 1
u 2
− uh 2H 1 () 6 g − Au L2 ( u − vh L2 + uh − v L2 ) 0
+ 12 u − vh 2H 1 () :
(3.2)
0
Next, for v ∈ K, we de$ne the linear interpolation rh v by r h v ∈ Sh ;
(rh v)(xi ) = v(xi );
i = 0; 1; : : : ; n; n + 1:
(3.3)
Note that rh v ∈ Kh . Then replacing vh by rh u in (3.2), we have 1
uh 2
− u 2H 1 () 6 g − Au L2 ( u − rh u L2 + uh − v L2 ) 0
+ 12 rh u − u 2H 1 () :
(3.4)
0
Therefore, by (3.2) and standard results of approximation theory [16], we have 1
rh u − u H01 () 6 h|u|H 2 () $ and 1
rh u − u L2 6 2 h2 |u|H 2 () : $ Next, in order to estimate uh − v L2 , it is convenient to introduce the function uh∗ = max{uh ; that uh∗ ¿ hold in . For both functions uh and being in the space H 1 (), it follows that also in H 1 (). Hence, the condition 6 0 on 9 implies that uh∗ ∈ H01 (). Thus the function an element of the set K. Let & = {x ∈ ; uh ¡ }. Since uh − uh∗ = 0 on − &, we have ∗ 2
uh − uh L2 = |uh − |2 d x:
(3.5)
(3.6) }, so uh∗ is uh∗ is
&
Now, we consider the following space such that Sh ⊂ Sh∗ by Sh∗ = {v ∈ H 1 (); v|ei ∈ P1 ; i = 1; 2; : : : ; n + 1}: For
∈ H 2 (), we de$ne the linear interpolation 'h 'h ∈ Sh∗ ;
(xi ) = 'h (xi );
by
i = 0; 1; 2; : : : ; n + 1:
Since, uh (xi ) ¿ (xi ) = 'h (xi ); i = 0; 1; 2; : : : ; n; n + 1, it follows that u h − 'h ¿ 0
in :
Consequently, ∀x ∈ &; 0 ¡ |( − uh )(x)| = ( − uh )(x) 6 ( − 'h )(x) = |( − 'h )(x)| and thus, ∗ 2 2
uh − uh L2 = |uh − | d x 6 | − 'h |2 d x 6 − 'h 2L2 : &
Since
'h − L2 6
1 2 h | |H 2 () $2
&
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
411
we obtain 1 2 h | |H 2 () : $2 Therefore, from (3.4)–(3.7) and replacing v by uh∗ in (3.4), we obtain
1 2 1 2 2
uh − u H 1 () 6 2( g L2 + |u|H 2 () ) h |u|H 2 () + 2 h | |H 2 () 0 $2 $
uh − uh∗ L2 6
1 2 2 h |u|H 2 () : $2 From the above arguments and (2.2), we can estimate C(g; ) in (3.1) by 1 2 |u|H 2 () + 2(||g||L2 + |u|H 2 () )(|u|H 2 () + | |H 2 () ): C(g; ) 6 $ +
(3.7)
(3.8)
(3.9)
We now consider the L2 estimates of uh − u via a generalization of the Aubin–Nitsche method. First, we note the following lemma in [14]. Lemma 3.2. Assume that u ∈ H 2 (); u ¿ 0 on . Then, there exists u˜ h ∈ Sh satisfying 0 6 u˜ h 6 u; such that
u − u˜ h H01 () 6
√ 1 4 2 h Mu L2 : + $ 3
The following result is given by arguments similar to those in [5]. Theorem 3.3. Let u and uh be the solutions of (2.1) and (2.6), respectively. If g ∈ L2 (), then we have
uh − u L2 6 C1 h uh − u H01 () ; where
C1 =
√ 1 4 2 + : $ 3
Proof. We omit the proof of this theorem because it is almost the same as in [14]. Therefore, we present the following theorem. Theorem 3.4. Let u and uh be the solutions of (2.1) and (2.6), respectively. If g ∈ L2 (), then √ 1 4 2 + h2 :
uh − u L2 6 C(g; ) $ 3 Here, C(g; ) is the same as in (3.1).
412
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
4. Computing procedures for verication We propose a computer algorithm to obtain a set U (i) which satis$es the veri$cation condition (2.9). Now, we consider the following auxiliary problem associated with (2.3), concerning any g ∈ L2 (): a(u; v − u) ¿ (g; v − u)
∀v ∈ K; u ∈ K:
(4.1)
We then de$ne the approximate problem corresponding to (4.1) as a(uh ; vh − uh ) ¿ (g; vh − uh )
∀vh ∈ Kh ; uh ∈ Kh :
(4.2)
Since the bilinear form a(·; ·) is symmetric, (4.2) is actually equivalent to the quadratic programming problem min [ 12 a(v; v) − (g; v)]:
(4.3)
v ∈ Kh
Let z = (zj ) ∈ Rn be the coeNcient vector for {j } corresponding to the function v in (4.3), and de$ne ˆ := ( j ) ∈ Rn , where j = (xj ); j = 1; 2; : : : ; n. Therefore, we can represent the above quadratic programming problem (4.3) in the form min [ 12 z T Dz − P T z]:
(4.4)
z¿ ˆ
Here, D = (dij ), with dij = (∇i ; ∇j ) and 1 6 i; j 6 n;. Further, P ≡ ((g; j )) is an n-dimensional vector. By the Kuhn–Tucker theorem [15], a vector z − ˆ = (zj ) − ( j ) ∈ Rn with z − ˆ ¿ 0 is an optimal solution to (4.4) if and only if there exists w = (wj ) ∈ Rn such that w − Dz = −(g; j ); w(z − ˆ ) = 0;
1 6 j 6 n;
w ¿ 0; z − ˆ ¿ 0:
Thus we can proceed in the following manner: Let R+ denote the set of all nonnegative real numbers. For
(4.5) ∈ R+ we de$ne
[ ] ≡ { ∈ H01 (); H01 () 6 ; L2 6 C1 h }:
(4.6)
This set corresponds n to the rounding error de$ned in Section 2. Now, let Aj (1 6 j 6 n) be intervals 1 on R and let j=1 Aj j be a linear combination of {j }, i.e., an element of the power set 2Sh in the following sense: n n Aj j = aj j ; aj ∈ Aj ; 1 6 j 6 n : j=1
Then, setting U =
j=1
n
j=1
Aj j + [ ] and g = f(U ) in (4.2), we consider the nonlinear system
w − Dz = −(f(U ); j ); w(z − ˆ ) = 0;
1 6 j 6 n;
w ¿ 0; z − ˆ ¿ 0:
(4.7)
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
413
Eq. (4.7) is in fact a bilinear system of equations whose right-hand side consists of intervals with constrained conditions w ¿ 0 and z − ˆ ¿ 0. To solve the nonlinear system (4.7) with automatic veri$cation of correctness of the result, veri$cation method for nonsmooth equations by a generalized Krawczyk operator as in [1] could be used. While, in order to enclose a solution for (4.7), we can use the same procedure as in [12]. Here we adopt another method as given below. Setting x = (w; z) ∈ R2n , (4.7) is written without constraint as a non-linear system of equations "(x) = 0:
(4.8)
2 ; : : : ; w n ; z1 ; z2 ; : : : ; zn ) be an approximate solution of (4.5). Then, note that wi ≈ 0 or Let x˜ := ( w1 ; w zi − i ≈ 0 for each 1 6 i 6 n. Problem (4.7) can also be reformulated as nonsmooth equations to use other methods, e.g., [1,2]. However, (4.8) is continuous and di:erentiable. Hence, to enclose solutions for (4.8), we use the following theorem proposed in [12]. Theorem 4.1. Let " : R2n → R2n be a function with continuous
−M · "(x) ˜ + {I − M · " (x˜ + X )} · X ⊆ X ˆ = 0. then there exists an xˆ ∈ x˜ + X ◦ with "(x)
Let X = (W; Z) be an enclosure of a solution of the nonlinear system (4.8) by using Theorem 4.1, where W := (W1 ; W2 ; : : : ; Wn ) ∈ IRn and Z := (Z1 ; Z2 ; : : : ; Zn ) ∈ IRn . Then we set Wi := 0 or Zi − i := 0 for each 1 6 i 6 n provided that wi ≈ 0 or zi − i ≈ 0, respectively. If, for all i; {Wi = 0 and inf (Zi − i ) ¿ 0} and {inf (Wi ) ¿ 0 and Zi − i := 0} hold, then it implies that problem (4.7) has an optimal solution x ∈ X (cf. [12]). As one can see, for the case that w˜ i and z˜i − i are both close to zero, this algorithm would not work. Fortunately, we have never encountered such a diNculty 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 [1]. We now consider the fully automatic computer generation of a sequence of sets {U (i) }; i=0; 1; : : : ; which consists of subsets of H01 (), in Section 2. We present an iterative procedure for generating {U (i) }i=0; 1; ::: . For i = 0, we choose appropriate initial values uh(0) ∈ Kh and 0 ∈ R+ , and de$ne U (0) ⊂ V by U (0) = uh(0) + [ 0 ]: Usually, uh(0) is determined as a(uh(0) ; vh − uh(0) ) ¿ (f(uh(0) ); vh − uh(0) )
∀vh ∈ Kh ; uh(0) ∈ Kh :
This corresponds to the Galerkin approximation for (2.3).
(4.9)
414
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
For Uh(i) = nj=1 A(i) j j and and i+1 ∈ R+ according to
i
∈ R+ , we set U (i) = Uh(i) + [ i ], i ¿ 1. Then, we de$ne Uh(i+1) ⊂ Kh
w − Dz = −(f(U (i) ); j ); w(z − ˆ ) = 0; i+1
1 6 j 6 n;
w ¿ 0; z − ˆ ¿ 0:
(4.10)
= sup C(f(u); )h:
(4.11)
u∈U (i)
Here, Uh(i+1) is determined as the solution set of (4.10), as described above. Thus, we de$ne n [−"; +"] j + [ i+1 + "] U i+1 := Uh(i+1) + j=1
and then we go back to the iteration scheme in Section 2. 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 PK F in (2.4) is retractive around the solution u and provided that the mesh size h is suNciently small (cf. [6–8]). We will consider it in the forthcoming paper about general case. 5. Numerical example We provide numerical examples of veri$cation in the one-dimensional case following the procedure described in the previous section. Let = (0; 1). We consider the case f(u) = Qu − sin 2$x, where Q is a constant, = sin $x and use the same $nite element subspace as in Section 3. We now choose the basis {i }ni=1 of Sh as the usual hat functions Fig. 1: Execution conditions Q = 3: dim Sh = 100: Extension parameter: " = 10−5 : Initial values: uh(0) = Galerkin approximation (4:9):
0
= 0:
The form of uh(0) is displayed in Fig: 1: Results Iteration numbers: N = 3: L2 -error bound := 0:01851: Maximum width of coeNcient intervals in {Aj(N ) } = 0:0000107858095: Remark 5.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,
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
415
Fig. 1. Approximate solution uh(0) .
PASCAL-XSC, FORTRAN-XSC, C-XSC, PROFIL, etc.). Therefore, it means that we neglected the round-o: 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−10 times larger than the round-o: errors. That is, there will be in general, some rounding errors at each step. Therefore, it is almost negligible compared with the truncation error which amount 10−3 ∼ 10−2 . Of course, we have to use those veri$cation software systems (e.g. PROFIL, INTLAB) in case that we need the rigorous mathematical proof. PROFIL is a portable C++ class fast interval library that supports an interval linear system solvers proposed by Rump. INTLAB is a Matlab toolbox supporting real and complex interval scalars, vectors, and matrices, as well as sparse real and complex interval matrices, coded in [13]. Acknowledgements The authors express their hearty thanks to the referees for their careful reading of the manuscript and for giving elaborate and useful comments to revise it, as well as to continue further study. References [1] X. Chen, A veri$cation method for solutions of nonsmooth equations, Computing 58 (1997) 281–294. [2] X. Chen, T. Yamamoto, On the convergence of some quasi-Newton methods for nonlinear equations with nondi:erentiable operators, Computing 48 (1992) 87–94. [3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [4] R.S. Falk, Error estimates for the approximation of a class of variational inequalities, Math. Comput. 28 (1974) 963–971. [5] U. Mosco, Error estimates for some variational inequalities, in: Lecture Notes in Mathematics, Vol. 606, Springer, Berlin, 1977. [6] M.T. Nakao, A numerical approach to the proof of existence of solutions for elliptic problems, Japan J. Ind. Appl. Math. 5 (1988) 313–332.
416
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
[7] M.T. Nakao, A numerical approach to the proof of existence of solutions for elliptic problems II, Japan J. Ind. Appl. Math. 7 (1990) 477–488. [8] M.T. Nakao, A numerical veri$cation method for the existence of weak solutions for nonlinear boundary value problems, J. Math. Anal. Appl. 164 (1992) 489–507. [9] M. Plum, Computer-assisted existence proofs for two-point boundary value problems, Computing 46 (1991) 19–34. [10] M. Plum, Existence and enclosure results for continua of solutions of parameter-dependent nonlinear boundary value problems, J. Comput. Appl. Math. 60 (1995) 187–200. [11] M. Plum, Computer assisted enclosure methods for elliptic di:erential equations, J. Linear Algebra Appl. 324 (2001) 147–187. [12] S.M. Rump, Solving algebraic problems with high accuracy, in: U.W. Kulish, W.L. Miranker (Eds.), A New Approach to Scienti$c Computation, Academic Press, New York, 1983, pp. 51–120. [13] S.M. Rump, INTLAB-INTerval LABoratory, in: T. Csendes (Ed.), Developments in Reliable Computing, Kluwer Academic Publishers, Dordrecht, 1999, pp. 77–104. [14] C.S. Ryoo, M.T. Nakao, Numerical veri$cation of solutions for variational inequalities, Numer. Math. 81 (1998) 305–320. [15] R.W.H. Sargent, An eNcient implementation of the Lemke algorithm and its extension to deal with upper and lower bounds, Math. Programming Study 7 (1978) 36–54. [16] G. Strang, G. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cli:s, NJ, 1973.
Numerical veri$cation of solutions for obstacle problems Cheon Seoung Ryooa;∗ , Mitsuhiro T. Nakaob a
Department of Mathematics, Hannam University, Daejeon 306-791, South Korea b Faculty of Mathematics, Kyushu University 33, Fukuoka 812-8581, Japan Received 19 August 2002; received in revised form 6 May 2003
Abstract This article is an extension of the previous paper (Numer. Math. 81 (1998) 305) by the same authors. We propose a method to prove the existence of solutions for more general obstacle problems by numerical computations. The main task in this paper consists of the numerical determination of some constants which appear in a priori error estimations for the $nite element approximation of a simple variational inequality. We present a numerical example for veri$cation. c 2003 Elsevier B.V. All rights reserved. Keywords: Numerical veri$cation; Obstacle problems; Error estimates; Quadratic programming
1. Introduction In recent years, several numerical methods have been proposed to verify the solutions of the differential equations. In the case of partial di:erential equations, several techniques have been studied in [6–11]. However, for some problems governed by the variational inequality, there are very few approaches. In the previous paper [14], we described a method to verify the existence of solutions for simple obstacle problems with convex set K ={v ∈ H01 () : v ¿ 0} based upon $nite element approximations, explicit a priori error estimates and Schauder’s $xed point theorem. In the present paper, we extend the method to more generalized obstacle problems with convex set K = {v ∈ H01 () ; v ¿ } and is the height of the obstacle. The main improvement compared with the previous paper [14] is that we newly established the constructive a priori error estimates of the $nite element approximation for the nonhomogeneous, which means is not identically 0, obstacle problems. It is readily seen that the proof technique used in [14] can no longer be applied to the present case. Actually, we
This work was supported by Korea Research Foundation Grant, KRF-2001-003-D00011. Corresponding author. Tel.: +82-42-629-7525; fax: +82-42-629-7894. E-mail addresses: [email protected] (C.S. Ryoo), [email protected] (M.T. Nakao).
∗
c 2003 Elsevier B.V. All rights reserved. 0377-0427/$ - see front matter doi:10.1016/j.cam.2003.05.006
406
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
have used quite a di:erent approach to get the desired error estimations. Also note that, in general, it is not possible to transform the concerned problem to the case that ≡ 0, which is actually a di:erent situation from the case of the elliptic equations. Getting constructive error estimates plays an essential role in our veri$cation method. The remainder of the paper is organized as follows. In Section 2, we describe the basic concept of numerical veri$cation methods which have been developed by the authors. Using a $xed point formulation, we can verify the existence of solutions via Schauder’s $xed point theorem. More speci$cally, we can verify a set which includes a solution. Such a set is constituted by the sum of two subsets in H01 (). One is the computation of a projection into a closed convex subset of some $nite-dimensional subspace, the other is the estimation of the error for the projection. Thus, the values of constants appearing in error estimations play an essential role in the veri$cation. We describe a method to estimate such constants numerically in Section 3. In this setting we describe, in Section 4, procedures verifying the existence of solutions by numerical enclosing. Finally, a numerical example for a generalized obstacle problem is presented. 2. The problem and the methods of verication Let be a bounded convex domain in Rn , 1 6 n 6 2, with piecewise smooth boundary 9. We set H01 () = {v ∈ H 1 (); v|9 = 0} and a(u; v) = ∇u · ∇v d x;
where ∇u · ∇v = 99xu1 99xv1 + 99xu2 99xv2 . Let A be the operator de$ned by Au = −Mu. Then, for u; v ∈ H01 (), we have a(u; v) = (Au; v), where (·; ·) denotes the L2 -inner product on . Denote · L2 as usual the L2 norm on . Next, we de$ne K = {v ∈ H01 (); v ¿ a:e: on }, where, is a given function in H 2 () such that 6 0 on 9. First, we note that, by the well-known result [4], for any g ∈ L2 , the problem a(u; v − u) ¿ (g; v − u);
∀v ∈ K; u ∈ K
(2.1)
has a unique solution u ∈ H01 () ∩ H 2 (), and the estimate |u|H 2 () 6 max{ g L2 ; }
(2.2)
holds. Here is the solution of the equation = max(A ; 0) L2 ( + g L2 )=( − g L2 ) and |u|H 2 implies the semi-norm of u on H 2 () de$ned by 2 92 u 2 2 |u|H 2 () ≡ 9xi 9xj 2 : L i; j=1 We adopt (∇; ∇ ) as the inner product on H01 (). Hence, the associated norm is de$ned by
H01 () = ∇ L2 . Now, let us consider the following obstacle problem: Find u ∈ K such that a(u; v − u) ¿ (f(u); v − u);
∀v ∈ K:
(2.3)
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
407
Here, assume that f satis$es the following hypotheses: A1. f is a continuous map from H01 () to L2 (). A2. For each bounded subset W ⊂ H01 (), f(W ) is also a bounded set in L2 (). To verify the existence of a solution of (2.3) on a computer, we use a $xed point formulation. Since a(·; ·) is a continuous bilinear form on H01 () × H01 (), for each u ∈ H01 (), from the Riesz representation theorem, there exists a unique element F(u) ∈ H01 () such that a(F(u); v) = (f(u); v) for all v ∈ H01 (). That is, ∃ F(u) ∈ H01 () such that − MF(u) = f(u) in ; F(u) = 0 on 9: Then the map F: H01 () → H01 () is a compact operator by the above assumptions on f. As in the preceding paper [14], de$ning the projection PK : H01 () → K, problem (2.3) is equivalent to that of $nding u ∈ H01 () such that u = PK F(u):
(2.4)
Now we describe a numerical veri$cation method to verify the existence of solutions of (2.3). We determine a set V ≡ PK F(U ) for a bounded, convex and closed subset U ⊂ H01 () as V = {v ∈ H01 (); v = PK F(u); ∀u ∈ U }:
(2.5)
From Schauder’s $xed point theorem, if V ⊂ U holds, then there exists a solution of (2.3) in the set U . Our goal is to $nd a set U which includes V . A procedure to verify V ⊂ U using a computer is similar to that in [14] as below. We now take a $nite element subspace Sh of H01 () for 0 ¡ h ¡ 1 with an appropriate mesh. Letting N denote the set of nodes associated with the space Sh , we then de$ne Kh , an approximation of K, by Kh = {vh ∈ Sh ; vh (p) ¿ (p); ∀p ∈ N}: Note that usually Kh ⊂ K, namely, we use an outer approximation. It is well known [3] that, for an arbitrary solution u ∈ H01 () of (2.1) and its $nite element approximation uh ∈ Kh de$ned by a(uh ; vh − uh ) ¿ (g; vh − uh )
∀vh ∈ Kh ;
(2.6)
there exists a C(g; ) such that
u − uh H01 () 6 C(g; )h:
(2.7)
We describe the method to estimate C(g; ) numerically in Section 3. For any u ∈ H01 (), we de$ne the rounding R(PK F(u)) ∈ Kh as the solution of the following problem: a(R(PK F(u)); vh − R(PK F(u))) ¿ (f(u); vh − R(PK F(u)))
∀vh ∈ Kh :
We de$ne the rounding R(V ) ⊂ Kh for the set V de$ned by (2.5) as R(V ) = {vh ∈ Kh ; vh = R(PK F(u)); u ∈ U }:
408
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
Also, we de$ne the rounding error RE(V ) ⊂ H01 () as 1 RE(V ) = v ∈ H0 (); v H01 () 6 sup C(f(u); )h : u ∈U
(2.8)
From the de$nition, we have V ⊂ R(V ) + RE(V ). Then it is suNcient to $nd U satisfying R(V ) + RE(V ) ⊂ U . Next let us introduce the procedure for $nding such a set U , which is referred as a candidate set, using computers. In order to $nd a set U satisfying the above condition, we use a simple iterative method as in [14]. First, we obtain an approximate solution uh(0) ∈ Kh to (2.3) by some appropriate method. Set U (0) = {uh(0) } and 0 = 0. Then, we iterate the following procedure: For i ¿ 0, and a candidate set U (i) = Uh(i) + [ i ], where i ∈ R+ and [ i ] = {v ∈ H01 (); v H01 () 6 i }, set V (i) = {v(i) ∈ K; v(i) = PK F(u(i) ); u(i) ∈ U (i) }: In order to enclose V (i) , let us de$ne R(V (i) ) as follows. R(V (i) ) is de$ned by the subset of Kh which consists of all elements vh(i) ∈ Kh such that a(vh(i) ; ! − vh(i) ) ¿ (f(u(i) ); ! − vh(i) );
∀! ∈ Kh
holds for some u(i) ∈ U (i) . Note that R(V (i) ) can be enclosed by R(V (i) ) ⊂ nj=1 Aj j where Aj = [Aj ; Aj ] are intervals, and {}nj=1 is a basis of Sh . Secondly, RE(V (i) ) is de$ned by RE(V (i) ) =
v ∈ H01 (); v H01 () 6 sup C(f(u(i) ); )h : u(i) ∈U (i)
Then the veri$cation condition is written as V (i) ⊂ R(V (i) ) + RE(V (i) ) ⊂ U (i) ;
(2.9)
which is implied by R(V (i) ) ⊂ Uh(i) and RE(V (i) ) ⊂ [ i ]. If the veri$cation condition is satis$ed, then U (i) is the desired set, and a solution to (2.3) exists in V (i) , and hence in U (i) . When the veri$cation condition is not satis$ed, we continue the simple iteration by using "-in6ation, i.e., let " be a certain positive constant given beforehand, and take i+1
[
= sup C(f(u(i) ); )h + ";
i+1 ]
Uh(i+1)
u(i) ∈U (i)
= {v ∈ H01 (); v H01 () 6 =
n
[ Aj − ";
i+1 };
Aj + " ]j ;
j=1
U (i+1) = Uh(i+1) + [
i+1 ]:
Then we continue the above iteration process (see [14]).
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
409
3. Computation of the value C (g; ) In this section, we consider the one-dimensional case. In order to verify solutions numerically, it is necessary to determine the value C(g; ) that appears in (2.7). Let = (0; 1) and g ∈ L2 . Now, let n be a positive integer and let h = 1=n + 1. We de$ne xi := ih for i = 0; 1; 2; : : : ; n + 1 (that is, a uniform partition of ) and ei := (xi−1 ; xi ); i = 1; 2; : : : ; n + 1. We then approximate H01 () by Sh = {vh ∈ C 0 (); vh (0) = vh (1) = 0; vh |ei ∈ P1 ; i = 1; 2; 3; : : : ; n + 1} with, as usual, P1 representing the space of polynomials of degree 6 1, thus dim Sh = n, and Kh = {vh ∈ Sh ; vh (xi ) ¿ (xi ); i = 0; 1; 2; : : : ; n; n + 1}. Let {j }j=1:::n be a basis of Sh such that j (x) ¿ 0; ∀x ∈ , which satis$es 1; i = j; j (xi ) = 0; i = j: Regarding the approximation error uh − u H01 () , we then have Theorem 3.1. Let u and uh be the solutions of (2.1) and (2.6), respectively. If g ∈ L2 (), then we have
uh − u H01 () 6 C(g; )h; where
(3.1)
1 2 |u|H 2 () + 2(||g||L2 + |u|H 2 () )(|u|H 2 () + | |H 2 () ): $ is estimated by g L2 and using (2.2).
C(g; ) 6 Here, |u|H 2 ()
Proof. We have a(uh − u; uh − u) 6 a(u; u) + a(uh ; uh ) − a(u; uh ) − a(uh ; u) and, using (2.1) and (2.6), a(u; u) 6 a(u; v) + (g; u − v)
∀v ∈ K;
a(uh ; uh ) 6 a(uh ; vh ) + (g; uh − vh )
∀vh ∈ Kh :
Hence, we deduce that, for all v ∈ K and all vh ∈ Kh , a(uh − u; uh − u) 6 a(u; v − uh ) + a(uh ; vh − u) + (g; u − v) + (g; uh − vh ) = a(u; v − uh ) − (g; v − uh ) + a(u; vh − u) − (g; vh − u) + a(uh − u; vh − u) = (g − Au; uh − v) + (g − Au; u − vh ) + a(u − uh ; u − vh ): Since g − Au ∈ L2 (), we obtain, for all v ∈ K and vh ∈ Kh , that
u − uh 2H 1 () 6 g − Au L2 ( u − vh L2 + uh − v L2 ) 0
+ u − uh H01 () u − vh H01 () :
410
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
Since
u − uh H01 () u − vh H01 () 6 12 ( u − uh 2H 1 () + u − vh 2H 1 () ) 0
0
we have, by combining the two previous inequalities, 1
u 2
− uh 2H 1 () 6 g − Au L2 ( u − vh L2 + uh − v L2 ) 0
+ 12 u − vh 2H 1 () :
(3.2)
0
Next, for v ∈ K, we de$ne the linear interpolation rh v by r h v ∈ Sh ;
(rh v)(xi ) = v(xi );
i = 0; 1; : : : ; n; n + 1:
(3.3)
Note that rh v ∈ Kh . Then replacing vh by rh u in (3.2), we have 1
uh 2
− u 2H 1 () 6 g − Au L2 ( u − rh u L2 + uh − v L2 ) 0
+ 12 rh u − u 2H 1 () :
(3.4)
0
Therefore, by (3.2) and standard results of approximation theory [16], we have 1
rh u − u H01 () 6 h|u|H 2 () $ and 1
rh u − u L2 6 2 h2 |u|H 2 () : $ Next, in order to estimate uh − v L2 , it is convenient to introduce the function uh∗ = max{uh ; that uh∗ ¿ hold in . For both functions uh and being in the space H 1 (), it follows that also in H 1 (). Hence, the condition 6 0 on 9 implies that uh∗ ∈ H01 (). Thus the function an element of the set K. Let & = {x ∈ ; uh ¡ }. Since uh − uh∗ = 0 on − &, we have ∗ 2
uh − uh L2 = |uh − |2 d x:
(3.5)
(3.6) }, so uh∗ is uh∗ is
&
Now, we consider the following space such that Sh ⊂ Sh∗ by Sh∗ = {v ∈ H 1 (); v|ei ∈ P1 ; i = 1; 2; : : : ; n + 1}: For
∈ H 2 (), we de$ne the linear interpolation 'h 'h ∈ Sh∗ ;
(xi ) = 'h (xi );
by
i = 0; 1; 2; : : : ; n + 1:
Since, uh (xi ) ¿ (xi ) = 'h (xi ); i = 0; 1; 2; : : : ; n; n + 1, it follows that u h − 'h ¿ 0
in :
Consequently, ∀x ∈ &; 0 ¡ |( − uh )(x)| = ( − uh )(x) 6 ( − 'h )(x) = |( − 'h )(x)| and thus, ∗ 2 2
uh − uh L2 = |uh − | d x 6 | − 'h |2 d x 6 − 'h 2L2 : &
Since
'h − L2 6
1 2 h | |H 2 () $2
&
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
411
we obtain 1 2 h | |H 2 () : $2 Therefore, from (3.4)–(3.7) and replacing v by uh∗ in (3.4), we obtain
1 2 1 2 2
uh − u H 1 () 6 2( g L2 + |u|H 2 () ) h |u|H 2 () + 2 h | |H 2 () 0 $2 $
uh − uh∗ L2 6
1 2 2 h |u|H 2 () : $2 From the above arguments and (2.2), we can estimate C(g; ) in (3.1) by 1 2 |u|H 2 () + 2(||g||L2 + |u|H 2 () )(|u|H 2 () + | |H 2 () ): C(g; ) 6 $ +
(3.7)
(3.8)
(3.9)
We now consider the L2 estimates of uh − u via a generalization of the Aubin–Nitsche method. First, we note the following lemma in [14]. Lemma 3.2. Assume that u ∈ H 2 (); u ¿ 0 on . Then, there exists u˜ h ∈ Sh satisfying 0 6 u˜ h 6 u; such that
u − u˜ h H01 () 6
√ 1 4 2 h Mu L2 : + $ 3
The following result is given by arguments similar to those in [5]. Theorem 3.3. Let u and uh be the solutions of (2.1) and (2.6), respectively. If g ∈ L2 (), then we have
uh − u L2 6 C1 h uh − u H01 () ; where
C1 =
√ 1 4 2 + : $ 3
Proof. We omit the proof of this theorem because it is almost the same as in [14]. Therefore, we present the following theorem. Theorem 3.4. Let u and uh be the solutions of (2.1) and (2.6), respectively. If g ∈ L2 (), then √ 1 4 2 + h2 :
uh − u L2 6 C(g; ) $ 3 Here, C(g; ) is the same as in (3.1).
412
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
4. Computing procedures for verication We propose a computer algorithm to obtain a set U (i) which satis$es the veri$cation condition (2.9). Now, we consider the following auxiliary problem associated with (2.3), concerning any g ∈ L2 (): a(u; v − u) ¿ (g; v − u)
∀v ∈ K; u ∈ K:
(4.1)
We then de$ne the approximate problem corresponding to (4.1) as a(uh ; vh − uh ) ¿ (g; vh − uh )
∀vh ∈ Kh ; uh ∈ Kh :
(4.2)
Since the bilinear form a(·; ·) is symmetric, (4.2) is actually equivalent to the quadratic programming problem min [ 12 a(v; v) − (g; v)]:
(4.3)
v ∈ Kh
Let z = (zj ) ∈ Rn be the coeNcient vector for {j } corresponding to the function v in (4.3), and de$ne ˆ := ( j ) ∈ Rn , where j = (xj ); j = 1; 2; : : : ; n. Therefore, we can represent the above quadratic programming problem (4.3) in the form min [ 12 z T Dz − P T z]:
(4.4)
z¿ ˆ
Here, D = (dij ), with dij = (∇i ; ∇j ) and 1 6 i; j 6 n;. Further, P ≡ ((g; j )) is an n-dimensional vector. By the Kuhn–Tucker theorem [15], a vector z − ˆ = (zj ) − ( j ) ∈ Rn with z − ˆ ¿ 0 is an optimal solution to (4.4) if and only if there exists w = (wj ) ∈ Rn such that w − Dz = −(g; j ); w(z − ˆ ) = 0;
1 6 j 6 n;
w ¿ 0; z − ˆ ¿ 0:
Thus we can proceed in the following manner: Let R+ denote the set of all nonnegative real numbers. For
(4.5) ∈ R+ we de$ne
[ ] ≡ { ∈ H01 (); H01 () 6 ; L2 6 C1 h }:
(4.6)
This set corresponds n to the rounding error de$ned in Section 2. Now, let Aj (1 6 j 6 n) be intervals 1 on R and let j=1 Aj j be a linear combination of {j }, i.e., an element of the power set 2Sh in the following sense: n n Aj j = aj j ; aj ∈ Aj ; 1 6 j 6 n : j=1
Then, setting U =
j=1
n
j=1
Aj j + [ ] and g = f(U ) in (4.2), we consider the nonlinear system
w − Dz = −(f(U ); j ); w(z − ˆ ) = 0;
1 6 j 6 n;
w ¿ 0; z − ˆ ¿ 0:
(4.7)
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
413
Eq. (4.7) is in fact a bilinear system of equations whose right-hand side consists of intervals with constrained conditions w ¿ 0 and z − ˆ ¿ 0. To solve the nonlinear system (4.7) with automatic veri$cation of correctness of the result, veri$cation method for nonsmooth equations by a generalized Krawczyk operator as in [1] could be used. While, in order to enclose a solution for (4.7), we can use the same procedure as in [12]. Here we adopt another method as given below. Setting x = (w; z) ∈ R2n , (4.7) is written without constraint as a non-linear system of equations "(x) = 0:
(4.8)
2 ; : : : ; w n ; z1 ; z2 ; : : : ; zn ) be an approximate solution of (4.5). Then, note that wi ≈ 0 or Let x˜ := ( w1 ; w zi − i ≈ 0 for each 1 6 i 6 n. Problem (4.7) can also be reformulated as nonsmooth equations to use other methods, e.g., [1,2]. However, (4.8) is continuous and di:erentiable. Hence, to enclose solutions for (4.8), we use the following theorem proposed in [12]. Theorem 4.1. Let " : R2n → R2n be a function with continuous
−M · "(x) ˜ + {I − M · " (x˜ + X )} · X ⊆ X ˆ = 0. then there exists an xˆ ∈ x˜ + X ◦ with "(x)
Let X = (W; Z) be an enclosure of a solution of the nonlinear system (4.8) by using Theorem 4.1, where W := (W1 ; W2 ; : : : ; Wn ) ∈ IRn and Z := (Z1 ; Z2 ; : : : ; Zn ) ∈ IRn . Then we set Wi := 0 or Zi − i := 0 for each 1 6 i 6 n provided that wi ≈ 0 or zi − i ≈ 0, respectively. If, for all i; {Wi = 0 and inf (Zi − i ) ¿ 0} and {inf (Wi ) ¿ 0 and Zi − i := 0} hold, then it implies that problem (4.7) has an optimal solution x ∈ X (cf. [12]). As one can see, for the case that w˜ i and z˜i − i are both close to zero, this algorithm would not work. Fortunately, we have never encountered such a diNculty 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 [1]. We now consider the fully automatic computer generation of a sequence of sets {U (i) }; i=0; 1; : : : ; which consists of subsets of H01 (), in Section 2. We present an iterative procedure for generating {U (i) }i=0; 1; ::: . For i = 0, we choose appropriate initial values uh(0) ∈ Kh and 0 ∈ R+ , and de$ne U (0) ⊂ V by U (0) = uh(0) + [ 0 ]: Usually, uh(0) is determined as a(uh(0) ; vh − uh(0) ) ¿ (f(uh(0) ); vh − uh(0) )
∀vh ∈ Kh ; uh(0) ∈ Kh :
This corresponds to the Galerkin approximation for (2.3).
(4.9)
414
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
For Uh(i) = nj=1 A(i) j j and and i+1 ∈ R+ according to
i
∈ R+ , we set U (i) = Uh(i) + [ i ], i ¿ 1. Then, we de$ne Uh(i+1) ⊂ Kh
w − Dz = −(f(U (i) ); j ); w(z − ˆ ) = 0; i+1
1 6 j 6 n;
w ¿ 0; z − ˆ ¿ 0:
(4.10)
= sup C(f(u); )h:
(4.11)
u∈U (i)
Here, Uh(i+1) is determined as the solution set of (4.10), as described above. Thus, we de$ne n [−"; +"] j + [ i+1 + "] U i+1 := Uh(i+1) + j=1
and then we go back to the iteration scheme in Section 2. 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 PK F in (2.4) is retractive around the solution u and provided that the mesh size h is suNciently small (cf. [6–8]). We will consider it in the forthcoming paper about general case. 5. Numerical example We provide numerical examples of veri$cation in the one-dimensional case following the procedure described in the previous section. Let = (0; 1). We consider the case f(u) = Qu − sin 2$x, where Q is a constant, = sin $x and use the same $nite element subspace as in Section 3. We now choose the basis {i }ni=1 of Sh as the usual hat functions Fig. 1: Execution conditions Q = 3: dim Sh = 100: Extension parameter: " = 10−5 : Initial values: uh(0) = Galerkin approximation (4:9):
0
= 0:
The form of uh(0) is displayed in Fig: 1: Results Iteration numbers: N = 3: L2 -error bound := 0:01851: Maximum width of coeNcient intervals in {Aj(N ) } = 0:0000107858095: Remark 5.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,
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
415
Fig. 1. Approximate solution uh(0) .
PASCAL-XSC, FORTRAN-XSC, C-XSC, PROFIL, etc.). Therefore, it means that we neglected the round-o: 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−10 times larger than the round-o: errors. That is, there will be in general, some rounding errors at each step. Therefore, it is almost negligible compared with the truncation error which amount 10−3 ∼ 10−2 . Of course, we have to use those veri$cation software systems (e.g. PROFIL, INTLAB) in case that we need the rigorous mathematical proof. PROFIL is a portable C++ class fast interval library that supports an interval linear system solvers proposed by Rump. INTLAB is a Matlab toolbox supporting real and complex interval scalars, vectors, and matrices, as well as sparse real and complex interval matrices, coded in [13]. Acknowledgements The authors express their hearty thanks to the referees for their careful reading of the manuscript and for giving elaborate and useful comments to revise it, as well as to continue further study. References [1] X. Chen, A veri$cation method for solutions of nonsmooth equations, Computing 58 (1997) 281–294. [2] X. Chen, T. Yamamoto, On the convergence of some quasi-Newton methods for nonlinear equations with nondi:erentiable operators, Computing 48 (1992) 87–94. [3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [4] R.S. Falk, Error estimates for the approximation of a class of variational inequalities, Math. Comput. 28 (1974) 963–971. [5] U. Mosco, Error estimates for some variational inequalities, in: Lecture Notes in Mathematics, Vol. 606, Springer, Berlin, 1977. [6] M.T. Nakao, A numerical approach to the proof of existence of solutions for elliptic problems, Japan J. Ind. Appl. Math. 5 (1988) 313–332.
416
C.S. Ryoo, M.T. Nakao / Journal of Computational and Applied Mathematics 161 (2003) 405 – 416
[7] M.T. Nakao, A numerical approach to the proof of existence of solutions for elliptic problems II, Japan J. Ind. Appl. Math. 7 (1990) 477–488. [8] M.T. Nakao, A numerical veri$cation method for the existence of weak solutions for nonlinear boundary value problems, J. Math. Anal. Appl. 164 (1992) 489–507. [9] M. Plum, Computer-assisted existence proofs for two-point boundary value problems, Computing 46 (1991) 19–34. [10] M. Plum, Existence and enclosure results for continua of solutions of parameter-dependent nonlinear boundary value problems, J. Comput. Appl. Math. 60 (1995) 187–200. [11] M. Plum, Computer assisted enclosure methods for elliptic di:erential equations, J. Linear Algebra Appl. 324 (2001) 147–187. [12] S.M. Rump, Solving algebraic problems with high accuracy, in: U.W. Kulish, W.L. Miranker (Eds.), A New Approach to Scienti$c Computation, Academic Press, New York, 1983, pp. 51–120. [13] S.M. Rump, INTLAB-INTerval LABoratory, in: T. Csendes (Ed.), Developments in Reliable Computing, Kluwer Academic Publishers, Dordrecht, 1999, pp. 77–104. [14] C.S. Ryoo, M.T. Nakao, Numerical veri$cation of solutions for variational inequalities, Numer. Math. 81 (1998) 305–320. [15] R.W.H. Sargent, An eNcient implementation of the Lemke algorithm and its extension to deal with upper and lower bounds, Math. Programming Study 7 (1978) 36–54. [16] G. Strang, G. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cli:s, NJ, 1973.