An algorithm for solving the double obstacle problems q
This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright
Author's personal copy
Available online at www.sciencedirect.com
Applied Mathematics and Computation 201 (2008) 221–228 www.elsevier.com/locate/amc
An algorithm for solving the double obstacle problems q Fei Wang *, Xiao-Liang Cheng Department of Mathematics, Zhejiang University, Hangzhou 310027, PR China
Abstract In this paper, we extend the algorithm in Xue and Cheng [L. Xue, X.L. Cheng, An algorithm for solving the obstacle problems, Comput. Math. Appl. 48 (2004) 1651–1657] for solving the double obstacle problem. We try to find the approximated region of the contact in the double obstacle problem with less computation time by iteration. Numerical example is given for the double obstacle problem for the elastic–plastic torsion problem to support the algorithm. Ó 2007 Elsevier Inc. All rights reserved. Keywords: The coincidence set; Double obstacle problem; Algorithm
1. Introduction Let X R2 . We consider the double obstacle problem: find u 2 K, such that EðuÞ ¼ min EðvÞ; v2K
ð1Þ
where K ¼ v 2 H 10 ðXÞ; u 6 v 6 w a:e in X ; ð2Þ Z Z 1 2 jrvj dx fvdx; ð3Þ EðvÞ ¼ 2 X X T and the obstacle functions u; w 2 H 1 ðXÞ CðXÞ; u 6 0; w P 0 on oX, which is the boundary of the domain X. Problems (1)–(3) denote the displacement u of elastic membrane lie between u and w under the distributed force f. The admissible set K is convex and v 2 H 10 ðXÞ means v ¼ 0 on oX. Eq. (3) expresses the total potential energy of the deformed membrane. In finding the position of equilibrium, the principle of minimum potential energy reduces the problem to look for, among all functions vðxÞ in admissible set K, the one which minimizes the functional (3).
q *
The projection was supported by National Science Foundation of China (Grant No. 10471129). Corresponding author. E-mail address: [email protected] (F. Wang).
0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.12.015
Author's personal copy
222
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
It is well know that (1)–(3) have the equivalent form: find u 2 K, such that Z Z ru rðv uÞdx P f ðv uÞdx 8v 2 K; X
ð4Þ
X
which is the first kind elliptic variational inequality. The theory and numerical solution T of variational inequalities have been widely studied in the monographs, e.g. [2–7]. If the solution u 2 C 2 ðXÞ H 10 ðXÞ, the variational inequality (4) also can be transformed to the boundary problem 8 Du P f ; X1 ¼ fx 2 XjuðxÞ ¼ uðxÞg; > > > < Du ¼ f ; X0 ¼ fx 2 XjuðxÞ < uðxÞ < wðxÞg; ð5Þ > Du 6 f ; X 2 ¼ fx 2 XjuðxÞ ¼ wðxÞg; > > : u 6 u 6 w in X; and u ¼ 0 on oX: X1 ðX2 Þ is the coincidence set about u with the lower (upper) obstacle. The obstacle problem is not just only introduced for a membrane in elasticity theory, but also for the non-parametric minimal and capillary surfaces as geometrical problems. The elastic–plastic torsion problem and the cavitation problem in the theory of lubrication also can be regarded as obstacle type problems. Since the obstacle problem is highly non-linear, the approximated solution is difficultly computed. The popular method for solving obstacle problems is variable projection method, such as the relaxation method [4], multilevel projection method [8], multigrid method [9– 12] and so on. The authors in [1] propose an algorithm, whose process is similar to simulating of moving obstacle but not same, to find the approximated region of the contact. We give a group of pictures to illustrate the idea of the algorithm for obstacle problem (only one obstacle) in which f ¼ 0. See Fig. 1, the line marked by little circle ‘‘” denotes the displacement of membrane and the line marked by star ‘‘*” is obstacle, they contact when the star ‘‘*” falls in circle ‘‘” completely. First of all, we compute the equation Du ¼ f to get the initial state ð1Þ without obstacle. Then we find X1 where the obstacle first touches the membrane. By setting ð1Þ ð1Þ uðxÞ ¼ uðxÞ in X1 we compute Du ¼ f in X n X1 again to get the next state of the membrane. Repeat this process until no new touch region is found. The intial state
The first step
12
12
10
10
8
8
6
6
4
4
2
2
0
0
−2
0
0.2
0.4
0.6
0.8
1
−2
0
0.2
The second step 12
10
10
8
8
6
6
4
4
2
2
0
0 0
0.2
0.4
0.6
0.6
0.8
1
0.8
1
The third step
12
−2
0.4
0.8
1
−2
0
Fig. 1. The process of the algorithm for
0.2
R X
0.4
0.6
ru rðv uÞdx P 0.
Author's personal copy
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
223
In this paper, we extend this idea for solving double obstacle problem (1)–(3). We present the algorithm in the following section and give one numerical experiment in Section 3. 2. Algorithm The double obstacle problem describes the equilibrium position u of an elastic membrane constrained to lie between two given obstacles u; w under an external force f. To find the approximated region of the contact X1 ; X2 . we first T 1consider the no obstacle case as the initial position of the elastic membrane 2 0 u ðxÞ 2 C ðXÞ H 0 ðXÞ, i.e. Du0 ðxÞ ¼ f
ð6Þ
in X:
We can imagine that the lower obstacle moves from 1 and the upper obstacle moves from þ1 to their ð1Þ final position u; w, respectively. Let lower obstacle move first, and find X1 , the part of lower coincidence set 0 X1 ,where the u ðxÞ first the ‘‘moving lower obstacle”. Then, we obtain the position of the elastic memT touches 2 1 1 brane, u ðxÞ 2 C ðXÞ H 0 ðXÞ, by solving ( ð1Þ Du1 ðxÞ ¼ f in X n X1 ; ð7Þ ð1Þ u1 ðxÞ ¼ uðxÞ in X1 : ð2Þ
ð1Þ
We again T get the contact part X1 , where u1 ðxÞ meets the ‘‘lower moving obstacle” in X n X1 . 2 2 u ðxÞ 2 C ðXÞ H 10 ðXÞ, solve 8 < Du2 ðxÞ ¼ f in X n Xð1Þ S Xð2Þ ; 1 1 ð8Þ : 2 ð1Þ S ð2Þ u ðxÞ ¼ uðxÞ in X1 X1 : Repeat the above procedure until no new touch region is need S ðkÞto be added. ð1Þ S Assume it take k steps of iteration to obtain X1 ¼ X1 X1 and uk . Then let upper obstacle move, we ð1Þ get X2 , the part of region of the contact X2 , where the upper obstacle touches the membrane uðxÞ firstly. By solving S 8 ð1Þ kþ1 > Du ðxÞ ¼ f in X n X ; X > 2 1 < ð9Þ ukþ1 ðxÞ ¼ uðxÞ in X1 ; > > : kþ1 ð1Þ u ðxÞ ¼ wðxÞ in X2 ; T we have the position of the membrane ukþ1 ðxÞ 2 C 2 ðXÞ H 10 ðXÞ. Attentively, the moving of upper obstacle ðkþ1Þ makes deformation of membrane, which brings new touch region X1 between the lower obstacle and the T 2 1 kþ2 membrane. u ðxÞ 2 C ðXÞ H 0 ðXÞ, we solve S 8 ð1Þ S ðkþ1Þ kþ2 > Du ðxÞ ¼ f in X n X ; X X > 2 1 1 < S ðkþ1Þ ð10Þ ukþ2 ðxÞ ¼ uðxÞ in X1 X1 ; > > : kþ1 ð1Þ u ðxÞ ¼ wðxÞ in X2 ; ð2Þ
then obtain the displacement of membrane ukþ2 , and X2 the region of contact with upper obstacle. Again repeat the above procedure until no new coincidence set comes to being. Ultimately we obtain the coincidence S ðkþ1Þ S ð1Þ S ð2Þ S set X1 ¼ X1 X1 and X2 ¼ X2 X2 . During this procedure we first let one obstacle move completely and then let the other obstacle move, so we call this procedure ‘‘One the other”. In what follows we introduce another ‘‘Up down” procedure to get the solution u and the coincidence sets T X1 ; X2 . At the beginning, u0 ðxÞ 2 C 2 ðXÞ H 10 ðXÞ, solve Du0 ðxÞ ¼ f
in X;
ð11Þ
Author's personal copy
224
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
we get u0 the displacement of elastic membrane under the external force f without obstacle as initial state. Let ð1Þ lower obstacle move, X1 is the part of the region of contact between u and u, where lower obstacle first touches membrane. By solving ( ð1Þ Du1 ðxÞ ¼ f in X n X1 ; ð12Þ ð1Þ u1 ðxÞ ¼ uðxÞ in X1 ; T move, we obtain displacement of membrane u1 ðxÞ 2 C 2 ðXÞ H 10 ðXÞ. Now lower obstacle stops, let the upper T ð1Þ after finding X2 the part of touch region of u with w, we can get the next state of u2 ðxÞ 2 C 2 ðXÞ H 10 ðXÞ by solving 8 ð1Þ S ð1Þ 2 > X2 Þ; > < Du ðxÞ ¼ f in X n ðX1 ð1Þ 2 ð13Þ u ðxÞ ¼ uðxÞ in X1 ; > > : 2 ð1Þ u ðxÞ ¼ wðxÞ in X2 ; it is upper obstacle’s turn to stop and lower obstacle moves. RepeatSthis procedure until no new StouchSregion S ðkÞ S isðkÞadded. We obtain the approximated coincidence set ð1Þ ð2Þ S ð1Þ ð2Þ X1 ¼ X1 X1 X1 ; X2 ¼ X2 X2 X2 , and solution u2k we want to have. The above procedure will stop in finite steps for the discreted double obstacle problem. Following, we present the algorithm for solving the discreted obstacle problem. Let V h H 10 ðXÞ be a finite element space. The discrete admissible set is K h ¼ vh 2 V h ;
uðxÞ 6 vh ðxÞ 6 wðxÞ;
ð14Þ
for any node x:
Then, the approximation of problem (4) is to find uh 2 K h , such that Z Z ruh rðvh uh Þdx P f ðvh uh Þdx 8vh 2 K h : X
ð15Þ
X
Denote /1R; /2 ; . . . ; /N , the basis functions at nodesRx1 ; x2 ; . . . ; xN of V h and A be the N N matrix with entries ai;j ¼ X r/i r/j dx. The load vector f is fi ¼ X f r/i dx. Now, we can rewrite (15) in matrix form: find u ¼ ðu1 ; . . . ; uN Þ 2 RN ; uðxi Þ 6 ui 6 wðxi Þ, such that 1 ðAu; v uÞ P ðf; v uÞ; ð16Þ 2 for all v ¼ ðv1 ; . . . ; vN Þ 2 RN ; uðxi Þ 6 vi 6 wðxi Þ; i ¼ 1; . . . ; N . Remark 1. As the author said in [1], their algorithm can be proved to be convergent in finite steps of iteration towards the exact solution of the discrete problem (2.7) in the same way as suggested in [13], the algorithm we proposed here also can be proved to be convergent in the same way. Following, we propose the algorithm based on the idea described above for solving the problem (16). Algorithm 2.1 (One the other type). Solving Au0 ¼ f, we get the initial position of the elastic membrane u0 without obstacle. Then, let k ¼ 0; l ¼ 1: (Lower obstacle move) ðkÞ
ðkÞ
ðkÞ
(i) computing d i ¼ uki uðxi Þ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i ðkÞ ðkÞ (ii) if d min P 0; X1 ¼ £, and go to (v), else n o ðkÞ ðkÞ X1 ¼ xi jd i ðkÞ ¼ d min ; ðkÞ
(iii) for all i 2 X1 , modify matrix A and load vector f by ai;i ¼ 1;
ai;j ¼ 0;
j 6¼ i;
f i ¼ uðxi Þ;
(iv) k ¼ k þ 1, solve Auk ¼ f. Go to (i) S ðnÞ ð0Þ S ð1Þ S (Suppose it takes n steps to go to (v), so we get X1 ¼ X1 X1 X1 and the unþ1 . Then upper obstacle move).
Author's personal copy
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228 ðnþ2l1Þ
ðnþlÞ
225
ðnþlÞ
(v) computing d i ¼ wðxi Þ uinþ2l1 ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i ðnþ2l1Þ ðlÞ (vi) if d min P 0; X2 ¼ £, go to (ix), else n o ðlÞ ðnþlÞ ðnþlÞ ¼ d min ; X2 ¼ xi jd i ðlÞ
(vii) for all i 2 X2 , modify matrix A and load vector f by ai;i ¼ 1; ai;j ¼ 0; j 6¼ i; f i ¼ wðxi Þ; (viii) solve Aunþ2l ¼ f ðnþ2lÞ ðnþ2lÞ ðnþ2lÞ (ix) computing d i ¼ uinþ2l uðxi Þ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i S ðnþ2lÞ ðnþlÞ ðlÞ ðnþlÞ (x) if d min P 0; X1 ¼ £, and if X2 ¼ £, STOP, else go to (v); else X1 n o ðnþlÞ ðnþlþ1Þ ðnþlþ1Þ ¼ xi jd i ¼ d min X1 ; ðnþlÞ
(xi) for all i 2 X1 , modify matrix A and load vector f by ai;i ¼ 1; ai;j ¼ 0; j 6¼ i; f i ¼ uðxi Þ; (xii) l ¼ l þ 1, solve Aunþ2l1 ¼ f, go to (v).Finally, we obtain the solution unþlþl and the lower region S ðnÞ S ðnþ1Þ S S ðnþlÞ ð0Þ S ð1Þ S X1 X1 X1 X1 , the upper coincidence set X2 ¼ of contact X1 ¼ X1 S ðlÞ ð1Þ S ð2Þ S X2 X2 . X2 Algorithm 2.2 (Up down type). Solving Au0 ¼ f, we get the initial position of the elastic membrane u0 . Then, for k ¼ 0; 1; . . . ð2kÞ
ð2kÞ
ð2kÞ
(i) computing d i ¼ u2k i uðxi Þ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i ð2kÞ ð2kÞ (ii) if d min P 0; X1 ¼ £, and go to (v), else n o ð2kÞ ð2kÞ ð2kÞ X1 ¼ xi jd i ¼ d min ; ð2kÞ
(iii) for all i 2 X1 , modify matrix A and load vector f by ai;i ¼ 1; ai;j ¼ 0; j 6¼ i; f i ¼ uðxi Þ; (iv) solve Au2kþ1 ¼ f ð2kþ1Þ ð2kþ1Þ ð2kþ1Þ (v) computing d i ¼ wðxi Þ ui2kþ1 ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i ð2kþ1Þ ð2kþ1Þ ð2kÞ S ð2kþ1Þ (vi) if d min P 0; X2 ¼ £, and if X1 ¼ £, STOP, else go to, (i); else X2 ð2kþ1Þ
X2
ð2kþ1Þ
¼ xi jd i
ð2kþ1Þ
¼ d min ;
ð2kþ1Þ
(vii) for all i 2 X2 , modify matrix A and load vector f by ai;i ¼ 1; ai;j ¼ 0; j 6¼ i; f i ¼ wðxi Þ; (viii) solve Au2kþ2 ¼ f. S ð2kÞ S ð2kþ1Þ ð0Þ S ð1Þ S Ultimately, we get the solution u2k ; X1 ¼ X1 . X2 ; X2 ¼ X2 X2 ðkÞ ðlÞ ðnþlÞ ð2kÞ ð2kþ1Þ In Algorithm 2.1 (Algorithm 2.2), we find very few new contact points in X1 ; X2 ; X1 X1 ; X2 in each iteration and need much iteration. To accelerate the algorithm, we introduce the small parameter e > 0. ðkÞ ðkÞ ðkÞ We think the nodes are contact points if d min 6 d i 6 d min þ e, thus, we obtain useful algorithm through takn o n o ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ing place of X1 ¼ xi jd i ¼ d min by X1 ¼ xi jd min 6 d i 6 d min þ e . In the above algorithm, we need to solve linear system for each iteration. We can use some know methods to solve the linear system, for instance, the fast Poisson solver, incomplete LU decomposition, preconditioned conjugate gradient method, multigrid method, and domain decomposition method. Remark 2. The parameter e in the algorithms determine the pace of ‘‘moving obstacle”, we do not set that it is too small to need much iteration and do not make it large to influence the precision of the algorithm as well.
Author's personal copy
226
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
For obtaining more exact solution as soon as possible, we can set the e larger in the beginning, and make it smaller in the last steps. Remark 3. In above algorithms, the touch region which is found in one step, is regarded as the part of contact region all the times in the following steps. But the algorithms in [13,14] are different, the touch point in one step may become non-contact point in the next steps. 3. Numerical result In this section, one numerical example is given for the elastic–plastic torsion problem. It is seen that the contact region of the obstacle problem is approximated by implementing the new algorithm on the computer, in which we suppose e ¼ 0:01, e.g. suppose X ¼ ½0; 1 ½0; 1 . Let uðx; yÞ ¼ distðx; oXÞ; wðx; yÞ ¼ 0:2, and set 8 if ðx; yÞ 2 S ¼ fðx; yÞ 2 X : jx yj 6 0:1 and x 6 0:3g; > < 300; f ðx; yÞ ¼ 70 expðyÞgðxÞ; if x 6 1 y and ðx; yÞ 62 S; > : 15 expðyÞgðxÞ; if x > 1 y and ðx; yÞ 62 S; where 8 6x; > > > > > 2ð1 3xÞ; > > > < 6ðx 1=3Þ; gðxÞ ¼ > 2ð1 3ðx 1=3ÞÞ; > > > > > 6ðx 2=3Þ; > > : 2ð1 3ðx 2=3ÞÞ;
if if
0 < x 6 1=6; 1=6 < x 6 1=3;
if
1=3 < x 6 1=2;
if if
1=2 < x 6 2=3; 2=3 < x 6 5=6;
if
5=6 < x 6 1:
Here we use the bilinear rectangle finite element space to approach H 10 ðXÞ. It is the 3D picture of the approximated solution u computed by ‘‘One the other” type algorithm in Fig. 2. In Figs. 3 and 4 the lower and upper coincidence sets are marked with star ‘‘*” and dot ‘‘”, respectively. They all have the similar figure with Fig. 3 in Example 5.4 in [14].
0.3 0.2 0.1 0 −0.1 −0.2
1 0.8
1 0.6
0.8 0.6
0.4
0.4
0.2
0.2 0
0
Fig. 2. Up down algorithm.
Author's personal copy
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
227
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.8
0.9
1
Fig. 3. One the other algorithm.
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Fig. 4. Up down algorithm.
Table 1 Comparison of the two algorithm Algorithm type
Consumed time (s)
Compute times of Au ¼ f
One the other Up down
1275.1 916.8
65 46
Author's personal copy
228
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
Comparing the two types of algorithm, we find the consumed times of ‘‘Up down” type is 71.9% of ‘‘One the other” type’s. Detailed information is in Table 1. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
L. Xue, X.L. Cheng, An algorithm for solving the obstacle problems, Comput. Math. Appl. 48 (2004) 1651–1657. G. Duvaut, J.-L. Lions, Inequalities in Mechanics and Physics, Springer-Verlag, Berlin, Germany, 1976. A. Friedman, Variational Principles and Free-boundary Problems, John Wiley, New York, 1982. R. Glowinski, Numerical Methods for Non-linear Variational Problems, Springer-Verlag, New York, 1984. N. Kikuchi, J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, PA, 1988. D. Kinderlehrer, G. Stampacchia, An Introduction to Variational Inequalities and their Applications, Academic Press, New York, 1980. W. Han, B.D. Reddy, Plasticity: Mathematical Theory and Numerical Analysis, Springer-Verlag, New York, 1999. Y. Zhang, Multilevel projection algorithm for solving obstacle problems, Comput. Math. Appl. 41 (12) (2001) 1505–1513. W. Hackbusch, A. Reusken, Analysis of a damped non-linear multilevel method, Numer. Math. 55 (1989) 225–246. R. Kornhuber, Monotone multigrid methods for elliptic variational inequalities I, Numer. Math. 69 (1994) 167–184. R. Kornhuber, Monotone multigrid methods for elliptic variational inequalities II, Numer. Math. 72 (1996) 481–499. B. Imoro, Discretized obstacle problems with penalties on nested grids, Appl. Numer. Math. 32 (2000) 21–34. R. Herbin, A monotinic method for the numerical solution of some free boundary value problems, SIAM J. Numer. Anal. 40 (2003) 2292–2310. T. Ka¨rkka¨inen, K. Kunisch, P. Tarvainen, Augmented Lagrangian active set methods for obstacle problems, J. Optimiz. Theory Appl. 119 (3) (2003) 499–533.
Author's personal copy
Available online at www.sciencedirect.com
Applied Mathematics and Computation 201 (2008) 221–228 www.elsevier.com/locate/amc
An algorithm for solving the double obstacle problems q Fei Wang *, Xiao-Liang Cheng Department of Mathematics, Zhejiang University, Hangzhou 310027, PR China
Abstract In this paper, we extend the algorithm in Xue and Cheng [L. Xue, X.L. Cheng, An algorithm for solving the obstacle problems, Comput. Math. Appl. 48 (2004) 1651–1657] for solving the double obstacle problem. We try to find the approximated region of the contact in the double obstacle problem with less computation time by iteration. Numerical example is given for the double obstacle problem for the elastic–plastic torsion problem to support the algorithm. Ó 2007 Elsevier Inc. All rights reserved. Keywords: The coincidence set; Double obstacle problem; Algorithm
1. Introduction Let X R2 . We consider the double obstacle problem: find u 2 K, such that EðuÞ ¼ min EðvÞ; v2K
ð1Þ
where K ¼ v 2 H 10 ðXÞ; u 6 v 6 w a:e in X ; ð2Þ Z Z 1 2 jrvj dx fvdx; ð3Þ EðvÞ ¼ 2 X X T and the obstacle functions u; w 2 H 1 ðXÞ CðXÞ; u 6 0; w P 0 on oX, which is the boundary of the domain X. Problems (1)–(3) denote the displacement u of elastic membrane lie between u and w under the distributed force f. The admissible set K is convex and v 2 H 10 ðXÞ means v ¼ 0 on oX. Eq. (3) expresses the total potential energy of the deformed membrane. In finding the position of equilibrium, the principle of minimum potential energy reduces the problem to look for, among all functions vðxÞ in admissible set K, the one which minimizes the functional (3).
q *
The projection was supported by National Science Foundation of China (Grant No. 10471129). Corresponding author. E-mail address: [email protected] (F. Wang).
0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.12.015
Author's personal copy
222
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
It is well know that (1)–(3) have the equivalent form: find u 2 K, such that Z Z ru rðv uÞdx P f ðv uÞdx 8v 2 K; X
ð4Þ
X
which is the first kind elliptic variational inequality. The theory and numerical solution T of variational inequalities have been widely studied in the monographs, e.g. [2–7]. If the solution u 2 C 2 ðXÞ H 10 ðXÞ, the variational inequality (4) also can be transformed to the boundary problem 8 Du P f ; X1 ¼ fx 2 XjuðxÞ ¼ uðxÞg; > > > < Du ¼ f ; X0 ¼ fx 2 XjuðxÞ < uðxÞ < wðxÞg; ð5Þ > Du 6 f ; X 2 ¼ fx 2 XjuðxÞ ¼ wðxÞg; > > : u 6 u 6 w in X; and u ¼ 0 on oX: X1 ðX2 Þ is the coincidence set about u with the lower (upper) obstacle. The obstacle problem is not just only introduced for a membrane in elasticity theory, but also for the non-parametric minimal and capillary surfaces as geometrical problems. The elastic–plastic torsion problem and the cavitation problem in the theory of lubrication also can be regarded as obstacle type problems. Since the obstacle problem is highly non-linear, the approximated solution is difficultly computed. The popular method for solving obstacle problems is variable projection method, such as the relaxation method [4], multilevel projection method [8], multigrid method [9– 12] and so on. The authors in [1] propose an algorithm, whose process is similar to simulating of moving obstacle but not same, to find the approximated region of the contact. We give a group of pictures to illustrate the idea of the algorithm for obstacle problem (only one obstacle) in which f ¼ 0. See Fig. 1, the line marked by little circle ‘‘” denotes the displacement of membrane and the line marked by star ‘‘*” is obstacle, they contact when the star ‘‘*” falls in circle ‘‘” completely. First of all, we compute the equation Du ¼ f to get the initial state ð1Þ without obstacle. Then we find X1 where the obstacle first touches the membrane. By setting ð1Þ ð1Þ uðxÞ ¼ uðxÞ in X1 we compute Du ¼ f in X n X1 again to get the next state of the membrane. Repeat this process until no new touch region is found. The intial state
The first step
12
12
10
10
8
8
6
6
4
4
2
2
0
0
−2
0
0.2
0.4
0.6
0.8
1
−2
0
0.2
The second step 12
10
10
8
8
6
6
4
4
2
2
0
0 0
0.2
0.4
0.6
0.6
0.8
1
0.8
1
The third step
12
−2
0.4
0.8
1
−2
0
Fig. 1. The process of the algorithm for
0.2
R X
0.4
0.6
ru rðv uÞdx P 0.
Author's personal copy
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
223
In this paper, we extend this idea for solving double obstacle problem (1)–(3). We present the algorithm in the following section and give one numerical experiment in Section 3. 2. Algorithm The double obstacle problem describes the equilibrium position u of an elastic membrane constrained to lie between two given obstacles u; w under an external force f. To find the approximated region of the contact X1 ; X2 . we first T 1consider the no obstacle case as the initial position of the elastic membrane 2 0 u ðxÞ 2 C ðXÞ H 0 ðXÞ, i.e. Du0 ðxÞ ¼ f
ð6Þ
in X:
We can imagine that the lower obstacle moves from 1 and the upper obstacle moves from þ1 to their ð1Þ final position u; w, respectively. Let lower obstacle move first, and find X1 , the part of lower coincidence set 0 X1 ,where the u ðxÞ first the ‘‘moving lower obstacle”. Then, we obtain the position of the elastic memT touches 2 1 1 brane, u ðxÞ 2 C ðXÞ H 0 ðXÞ, by solving ( ð1Þ Du1 ðxÞ ¼ f in X n X1 ; ð7Þ ð1Þ u1 ðxÞ ¼ uðxÞ in X1 : ð2Þ
ð1Þ
We again T get the contact part X1 , where u1 ðxÞ meets the ‘‘lower moving obstacle” in X n X1 . 2 2 u ðxÞ 2 C ðXÞ H 10 ðXÞ, solve 8 < Du2 ðxÞ ¼ f in X n Xð1Þ S Xð2Þ ; 1 1 ð8Þ : 2 ð1Þ S ð2Þ u ðxÞ ¼ uðxÞ in X1 X1 : Repeat the above procedure until no new touch region is need S ðkÞto be added. ð1Þ S Assume it take k steps of iteration to obtain X1 ¼ X1 X1 and uk . Then let upper obstacle move, we ð1Þ get X2 , the part of region of the contact X2 , where the upper obstacle touches the membrane uðxÞ firstly. By solving S 8 ð1Þ kþ1 > Du ðxÞ ¼ f in X n X ; X > 2 1 < ð9Þ ukþ1 ðxÞ ¼ uðxÞ in X1 ; > > : kþ1 ð1Þ u ðxÞ ¼ wðxÞ in X2 ; T we have the position of the membrane ukþ1 ðxÞ 2 C 2 ðXÞ H 10 ðXÞ. Attentively, the moving of upper obstacle ðkþ1Þ makes deformation of membrane, which brings new touch region X1 between the lower obstacle and the T 2 1 kþ2 membrane. u ðxÞ 2 C ðXÞ H 0 ðXÞ, we solve S 8 ð1Þ S ðkþ1Þ kþ2 > Du ðxÞ ¼ f in X n X ; X X > 2 1 1 < S ðkþ1Þ ð10Þ ukþ2 ðxÞ ¼ uðxÞ in X1 X1 ; > > : kþ1 ð1Þ u ðxÞ ¼ wðxÞ in X2 ; ð2Þ
then obtain the displacement of membrane ukþ2 , and X2 the region of contact with upper obstacle. Again repeat the above procedure until no new coincidence set comes to being. Ultimately we obtain the coincidence S ðkþ1Þ S ð1Þ S ð2Þ S set X1 ¼ X1 X1 and X2 ¼ X2 X2 . During this procedure we first let one obstacle move completely and then let the other obstacle move, so we call this procedure ‘‘One the other”. In what follows we introduce another ‘‘Up down” procedure to get the solution u and the coincidence sets T X1 ; X2 . At the beginning, u0 ðxÞ 2 C 2 ðXÞ H 10 ðXÞ, solve Du0 ðxÞ ¼ f
in X;
ð11Þ
Author's personal copy
224
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
we get u0 the displacement of elastic membrane under the external force f without obstacle as initial state. Let ð1Þ lower obstacle move, X1 is the part of the region of contact between u and u, where lower obstacle first touches membrane. By solving ( ð1Þ Du1 ðxÞ ¼ f in X n X1 ; ð12Þ ð1Þ u1 ðxÞ ¼ uðxÞ in X1 ; T move, we obtain displacement of membrane u1 ðxÞ 2 C 2 ðXÞ H 10 ðXÞ. Now lower obstacle stops, let the upper T ð1Þ after finding X2 the part of touch region of u with w, we can get the next state of u2 ðxÞ 2 C 2 ðXÞ H 10 ðXÞ by solving 8 ð1Þ S ð1Þ 2 > X2 Þ; > < Du ðxÞ ¼ f in X n ðX1 ð1Þ 2 ð13Þ u ðxÞ ¼ uðxÞ in X1 ; > > : 2 ð1Þ u ðxÞ ¼ wðxÞ in X2 ; it is upper obstacle’s turn to stop and lower obstacle moves. RepeatSthis procedure until no new StouchSregion S ðkÞ S isðkÞadded. We obtain the approximated coincidence set ð1Þ ð2Þ S ð1Þ ð2Þ X1 ¼ X1 X1 X1 ; X2 ¼ X2 X2 X2 , and solution u2k we want to have. The above procedure will stop in finite steps for the discreted double obstacle problem. Following, we present the algorithm for solving the discreted obstacle problem. Let V h H 10 ðXÞ be a finite element space. The discrete admissible set is K h ¼ vh 2 V h ;
uðxÞ 6 vh ðxÞ 6 wðxÞ;
ð14Þ
for any node x:
Then, the approximation of problem (4) is to find uh 2 K h , such that Z Z ruh rðvh uh Þdx P f ðvh uh Þdx 8vh 2 K h : X
ð15Þ
X
Denote /1R; /2 ; . . . ; /N , the basis functions at nodesRx1 ; x2 ; . . . ; xN of V h and A be the N N matrix with entries ai;j ¼ X r/i r/j dx. The load vector f is fi ¼ X f r/i dx. Now, we can rewrite (15) in matrix form: find u ¼ ðu1 ; . . . ; uN Þ 2 RN ; uðxi Þ 6 ui 6 wðxi Þ, such that 1 ðAu; v uÞ P ðf; v uÞ; ð16Þ 2 for all v ¼ ðv1 ; . . . ; vN Þ 2 RN ; uðxi Þ 6 vi 6 wðxi Þ; i ¼ 1; . . . ; N . Remark 1. As the author said in [1], their algorithm can be proved to be convergent in finite steps of iteration towards the exact solution of the discrete problem (2.7) in the same way as suggested in [13], the algorithm we proposed here also can be proved to be convergent in the same way. Following, we propose the algorithm based on the idea described above for solving the problem (16). Algorithm 2.1 (One the other type). Solving Au0 ¼ f, we get the initial position of the elastic membrane u0 without obstacle. Then, let k ¼ 0; l ¼ 1: (Lower obstacle move) ðkÞ
ðkÞ
ðkÞ
(i) computing d i ¼ uki uðxi Þ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i ðkÞ ðkÞ (ii) if d min P 0; X1 ¼ £, and go to (v), else n o ðkÞ ðkÞ X1 ¼ xi jd i ðkÞ ¼ d min ; ðkÞ
(iii) for all i 2 X1 , modify matrix A and load vector f by ai;i ¼ 1;
ai;j ¼ 0;
j 6¼ i;
f i ¼ uðxi Þ;
(iv) k ¼ k þ 1, solve Auk ¼ f. Go to (i) S ðnÞ ð0Þ S ð1Þ S (Suppose it takes n steps to go to (v), so we get X1 ¼ X1 X1 X1 and the unþ1 . Then upper obstacle move).
Author's personal copy
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228 ðnþ2l1Þ
ðnþlÞ
225
ðnþlÞ
(v) computing d i ¼ wðxi Þ uinþ2l1 ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i ðnþ2l1Þ ðlÞ (vi) if d min P 0; X2 ¼ £, go to (ix), else n o ðlÞ ðnþlÞ ðnþlÞ ¼ d min ; X2 ¼ xi jd i ðlÞ
(vii) for all i 2 X2 , modify matrix A and load vector f by ai;i ¼ 1; ai;j ¼ 0; j 6¼ i; f i ¼ wðxi Þ; (viii) solve Aunþ2l ¼ f ðnþ2lÞ ðnþ2lÞ ðnþ2lÞ (ix) computing d i ¼ uinþ2l uðxi Þ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i S ðnþ2lÞ ðnþlÞ ðlÞ ðnþlÞ (x) if d min P 0; X1 ¼ £, and if X2 ¼ £, STOP, else go to (v); else X1 n o ðnþlÞ ðnþlþ1Þ ðnþlþ1Þ ¼ xi jd i ¼ d min X1 ; ðnþlÞ
(xi) for all i 2 X1 , modify matrix A and load vector f by ai;i ¼ 1; ai;j ¼ 0; j 6¼ i; f i ¼ uðxi Þ; (xii) l ¼ l þ 1, solve Aunþ2l1 ¼ f, go to (v).Finally, we obtain the solution unþlþl and the lower region S ðnÞ S ðnþ1Þ S S ðnþlÞ ð0Þ S ð1Þ S X1 X1 X1 X1 , the upper coincidence set X2 ¼ of contact X1 ¼ X1 S ðlÞ ð1Þ S ð2Þ S X2 X2 . X2 Algorithm 2.2 (Up down type). Solving Au0 ¼ f, we get the initial position of the elastic membrane u0 . Then, for k ¼ 0; 1; . . . ð2kÞ
ð2kÞ
ð2kÞ
(i) computing d i ¼ u2k i uðxi Þ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i ð2kÞ ð2kÞ (ii) if d min P 0; X1 ¼ £, and go to (v), else n o ð2kÞ ð2kÞ ð2kÞ X1 ¼ xi jd i ¼ d min ; ð2kÞ
(iii) for all i 2 X1 , modify matrix A and load vector f by ai;i ¼ 1; ai;j ¼ 0; j 6¼ i; f i ¼ uðxi Þ; (iv) solve Au2kþ1 ¼ f ð2kþ1Þ ð2kþ1Þ ð2kþ1Þ (v) computing d i ¼ wðxi Þ ui2kþ1 ; i ¼ 1; 2; . . . ; N ; and d min ¼ min16i6N d i ð2kþ1Þ ð2kþ1Þ ð2kÞ S ð2kþ1Þ (vi) if d min P 0; X2 ¼ £, and if X1 ¼ £, STOP, else go to, (i); else X2 ð2kþ1Þ
X2
ð2kþ1Þ
¼ xi jd i
ð2kþ1Þ
¼ d min ;
ð2kþ1Þ
(vii) for all i 2 X2 , modify matrix A and load vector f by ai;i ¼ 1; ai;j ¼ 0; j 6¼ i; f i ¼ wðxi Þ; (viii) solve Au2kþ2 ¼ f. S ð2kÞ S ð2kþ1Þ ð0Þ S ð1Þ S Ultimately, we get the solution u2k ; X1 ¼ X1 . X2 ; X2 ¼ X2 X2 ðkÞ ðlÞ ðnþlÞ ð2kÞ ð2kþ1Þ In Algorithm 2.1 (Algorithm 2.2), we find very few new contact points in X1 ; X2 ; X1 X1 ; X2 in each iteration and need much iteration. To accelerate the algorithm, we introduce the small parameter e > 0. ðkÞ ðkÞ ðkÞ We think the nodes are contact points if d min 6 d i 6 d min þ e, thus, we obtain useful algorithm through takn o n o ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ing place of X1 ¼ xi jd i ¼ d min by X1 ¼ xi jd min 6 d i 6 d min þ e . In the above algorithm, we need to solve linear system for each iteration. We can use some know methods to solve the linear system, for instance, the fast Poisson solver, incomplete LU decomposition, preconditioned conjugate gradient method, multigrid method, and domain decomposition method. Remark 2. The parameter e in the algorithms determine the pace of ‘‘moving obstacle”, we do not set that it is too small to need much iteration and do not make it large to influence the precision of the algorithm as well.
Author's personal copy
226
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
For obtaining more exact solution as soon as possible, we can set the e larger in the beginning, and make it smaller in the last steps. Remark 3. In above algorithms, the touch region which is found in one step, is regarded as the part of contact region all the times in the following steps. But the algorithms in [13,14] are different, the touch point in one step may become non-contact point in the next steps. 3. Numerical result In this section, one numerical example is given for the elastic–plastic torsion problem. It is seen that the contact region of the obstacle problem is approximated by implementing the new algorithm on the computer, in which we suppose e ¼ 0:01, e.g. suppose X ¼ ½0; 1 ½0; 1 . Let uðx; yÞ ¼ distðx; oXÞ; wðx; yÞ ¼ 0:2, and set 8 if ðx; yÞ 2 S ¼ fðx; yÞ 2 X : jx yj 6 0:1 and x 6 0:3g; > < 300; f ðx; yÞ ¼ 70 expðyÞgðxÞ; if x 6 1 y and ðx; yÞ 62 S; > : 15 expðyÞgðxÞ; if x > 1 y and ðx; yÞ 62 S; where 8 6x; > > > > > 2ð1 3xÞ; > > > < 6ðx 1=3Þ; gðxÞ ¼ > 2ð1 3ðx 1=3ÞÞ; > > > > > 6ðx 2=3Þ; > > : 2ð1 3ðx 2=3ÞÞ;
if if
0 < x 6 1=6; 1=6 < x 6 1=3;
if
1=3 < x 6 1=2;
if if
1=2 < x 6 2=3; 2=3 < x 6 5=6;
if
5=6 < x 6 1:
Here we use the bilinear rectangle finite element space to approach H 10 ðXÞ. It is the 3D picture of the approximated solution u computed by ‘‘One the other” type algorithm in Fig. 2. In Figs. 3 and 4 the lower and upper coincidence sets are marked with star ‘‘*” and dot ‘‘”, respectively. They all have the similar figure with Fig. 3 in Example 5.4 in [14].
0.3 0.2 0.1 0 −0.1 −0.2
1 0.8
1 0.6
0.8 0.6
0.4
0.4
0.2
0.2 0
0
Fig. 2. Up down algorithm.
Author's personal copy
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
227
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.8
0.9
1
Fig. 3. One the other algorithm.
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Fig. 4. Up down algorithm.
Table 1 Comparison of the two algorithm Algorithm type
Consumed time (s)
Compute times of Au ¼ f
One the other Up down
1275.1 916.8
65 46
Author's personal copy
228
F. Wang, X.-L. Cheng / Applied Mathematics and Computation 201 (2008) 221–228
Comparing the two types of algorithm, we find the consumed times of ‘‘Up down” type is 71.9% of ‘‘One the other” type’s. Detailed information is in Table 1. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
L. Xue, X.L. Cheng, An algorithm for solving the obstacle problems, Comput. Math. Appl. 48 (2004) 1651–1657. G. Duvaut, J.-L. Lions, Inequalities in Mechanics and Physics, Springer-Verlag, Berlin, Germany, 1976. A. Friedman, Variational Principles and Free-boundary Problems, John Wiley, New York, 1982. R. Glowinski, Numerical Methods for Non-linear Variational Problems, Springer-Verlag, New York, 1984. N. Kikuchi, J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, PA, 1988. D. Kinderlehrer, G. Stampacchia, An Introduction to Variational Inequalities and their Applications, Academic Press, New York, 1980. W. Han, B.D. Reddy, Plasticity: Mathematical Theory and Numerical Analysis, Springer-Verlag, New York, 1999. Y. Zhang, Multilevel projection algorithm for solving obstacle problems, Comput. Math. Appl. 41 (12) (2001) 1505–1513. W. Hackbusch, A. Reusken, Analysis of a damped non-linear multilevel method, Numer. Math. 55 (1989) 225–246. R. Kornhuber, Monotone multigrid methods for elliptic variational inequalities I, Numer. Math. 69 (1994) 167–184. R. Kornhuber, Monotone multigrid methods for elliptic variational inequalities II, Numer. Math. 72 (1996) 481–499. B. Imoro, Discretized obstacle problems with penalties on nested grids, Appl. Numer. Math. 32 (2000) 21–34. R. Herbin, A monotinic method for the numerical solution of some free boundary value problems, SIAM J. Numer. Anal. 40 (2003) 2292–2310. T. Ka¨rkka¨inen, K. Kunisch, P. Tarvainen, Augmented Lagrangian active set methods for obstacle problems, J. Optimiz. Theory Appl. 119 (3) (2003) 499–533.
