SOLUTION METHOD FOR DISCRETE DOUBLE OBSTACLE PROBLEMS BASED ON A POWER PENALTY APPROACH

. We develop a power penalty approach to a ﬁnite-dimensional double obstacle problem. This problem is ﬁrst approximated by a system of non- linear equations containing two penalty terms. We show that the solution to this penalized equation converges to that of the original obstacle problem at an exponential rate when the coeﬃcient matrices are M -matrices. Numerical examples are presented to conﬁrm the theoretical ﬁndings and illustrate the eﬃciency and eﬀectiveness of the new method.

iteration method is another popular method to solve the discrete double obstacle problems [2]. Owing to its equivalence to the generalized Newton method, the policy iteration method performs well in most cases. However, as stated in [7], there exist some cases where the policy iteration method loses its efficiency. For instance, in our numerical example 2 the policy iteration method fails to compute a solution. Active set method is proposed to solve Problem 1 as well in [8], where the method is verified to be very fast, when the scale of the problem is large enough. Despite of its fastness, the active set method requires much more computer memory comparing with other methods, which makes it less attractive.
It is easy to see that Problem 1 can also be stated as the following equivalent complementarity form: The above complementarity form inspires us to apply the power penalty method to solve Problem 1, since the method has been demonstrated to be a very successful means to solve general complementarity problems, see, for example [12,14,15,17]. The advantages of the power penalty method lie in several aspects. First, it is easy to implement and robust to the scale of the problem. Second, it can be applied to more complicated problems, such as nonlinear complementarity problem [13], Hamilton-Jacobi-Bellman problem [16], etc. Third, the accuracy of the penalty method can be simply controlled by the penalty parameters with an exponential convergence rate.
Despite of its popularity in approximating standard complementarity problems, the power penalty approach to double obstacle problems has drawn little attention.
There are only few works are contributed to this area in the literature. In [10], a higher order penalty method is proposed to solve the continuous double obstacle problems arising from the convertible bond pricing. Zhou et al. [19] present a power penalty approximation for a semilinear parabolic double obstacle problem in an infinite dimensional space. In the finite dimensional space, Wang et al. [13] propose a power penalty approach to a discrete bounded variational inequalities. Under some conditions a discrete double obstacle problem can be equivalently transformed into a bounded variational inequality. Hence, the power penalty method can be indirectly applied to solve the double obstacle problem. However, this transformation makes the power penalization less intuitive and inconvenient to implementation.
This work mainly contributes to the direct application of the power penalty method to the discrete double obstacle problem. By penalizing the upper obstacle g and the lower obstacle f simultaneously, we develop a power penalty method, especially the lower order penalty method, to the discrete double obstacle problem (1). More importantly, an exponentially convergence rate of the power penalty method with respect to (w.r.t.) the penalty parameter is established. Another contribution of this paper is to design an efficient computational method to solve the discrete double obstacle problem, based on the lower order penalty approach. The computational efficiency of the new method is verified by numerical comparison with the widely-used policy iteration method.
We comment that in [16] we propose a power penalty approach to the discrete Hamilton-Jacobi-Bellman equation. However, the discrete double obstacle problem (1) is an Hamilton-Jacobi-Bellman-Issac equation, or an Hamilton-Jacobi-Bellman complementarity problem, rather than an Hamilton-Jacobi-Bellman equation. Thus, the results developed in [16] do not apply to Problem 1 any longer. As can be seen below, both the formulation and analysis for the discrete double obstacle problem (1) are different from that in [16]. We also comment that the interior penalty approach to the discrete mixed complementarity problem, developed in [18], is a promising method to solve (1). However, to guarantee its convergence, this method usually requires much stronger conditions than the power penalty method. This limits the application of the interior penalty method to more general double obstacle problems.
The reminder of this paper is organized as follows. We will first develop a power penalty approach to Problem 1 in Section 2. The solvability of the penalized problem is also given therein. In Section 3, the exponential convergence rate of the power penalty method w.r.t. the penalty parameter is established. Section 4 gives a solution algorithm for the penalized problem. In Section 5, three numerical examples are presented to illustrate the convergence rate, efficiency and effectiveness of the new method. Finally, conclusions are given in Section 6.
Before further discussion, we first make the following assumption on the matrix A: The above assumption means that a ij ≤ 0, a ii ≥ i =j |a ij | and for at least one i, a ii > i =j |a ij |. This implies that A is an M -matrix. It is also worth noting that the above assumption is normally guaranteed by a proper discretization method such as the upwind finite/difference/finite element or a fitted finite volume method for the differential equations in the continuous problem associated with Problem 1 (see, for example, [17]).
In the rest of our discussion, we assume that the assumption (A) is fulfilled. Under this assumption, the unique solvability of Problem 1 is a direct consequence of Theorem 5.1 of [2].
Problem 2 is the penalization corresponding to Problem 1.1, in which the first penalty term penalizes the violation of lower obstacle constraint while the second term penalizes the violation of the upper obstacle constraint. The essence of this penalization is to force all of the constraints to be satisfied when λ → ∞. We expect that any solution x λ to Problem 2 converges to that of Problem 1 when λ → ∞. Clearly, the rate of convergence depends on the parameter in the penalty equation. In the following sections we will establish the convergence rate for the penalty approach (2) w.r.t. the penalty parameter λ.
Before considering the convergence of this power penalty approach, we first show that there exists a solution to Problem 2. We start this discussion with the following lemma.
Lemma 2.1. Let x λ be the solution to Problem 2. Then there exists a positive constant C, independent of λ, such that where · ∞ denotes the usual ∞ norm on a finite dimensional space.
Proof. Rearranging (2), we get 1/k + . By defining three disjoint nonempty index subsets of I := {1, 2, · · · , N }: we consider the components of the above equation in the following three cases.
Now, we construct an N × N matrix A * 1 and an N × 1 matrix b * 1 such that their jth rows satisfy respectively, for any j ∈ I, where I denotes the N × N identity matrix. From the above three cases we have the following inequality.
Similarly, we introduce matrix A * 2 and b * 2 such that their jth rows satisfy Obviously, we have From the constructions of A * 1 and A * 2 , it follows that these two new matrices are also irreducibly diagonally dominant M -matrices. Hence, (A * 1 ) −1 > 0 and (5) and (6) we get Thus, x λ is bounded both below and above, and so there exists a positive constant C, independent of λ, and x λ , such that (3) holds.
With the above lemma, we are now ready to establish the unique solvability results of Problem 2 in the following Proposition.
Proposition 1. For any λ > 0, there exists a unique solution x λ to Problem 2.
Proof. For notation-simplicity, we omit the subscript λ of x λ in this proof. We show that Problem 2 has a solution in the bounded region S := {x ∈ R N : −ε −1 e < x < δ −1 e}, where e = (1, ..., 1) ∈ R N and ε and δ are (sufficiently small) positive constants. Let To prove this theorem, it suffices to verify the conditions of Miranda's theorem 1 . We first show that F (x) = 0 for x on the boundary ∂S of S. More specifically, we will show that 0 / ∈ F (∂S) when both ε > 0 and δ > 0 are sufficiently small. To prove this, we assume that 0 ∈ F (∂S), that is, there exists an x ∈ ∂S such that F (x) = 0. Then, we show this is not possible when both δ and ε are sufficiently small in the following two cases: From Assumption (A), we know that A is an irreducibly diagonally dominant Mmatrix. Hence, a ll > 0 and a l,j ≤ 0 for l = j. Thus, combining x λ ∞ ≤ C, we have x l ≤ 0 as δ → +0. This is contradicted to the assumption x l = δ −1 > 0 as δ → +0, and thus we conclude that when δ > 0 is sufficiently small, 0 / ∈ F (∂S) with x l = δ −1 for a feasible index l.
Case 2. We now consider the case x λ on ∂S such that at least one component of x λ is equal to O(−ε −1 ), say x l = −ε −1 for a feasible index l. In this case we have Using the similar argument as in Case 1, we also have that x l ≥ 0 as ε → +0, which is contradicted to the assumption x l = −ε −1 ≤ 0 as ε → +0, Combining the above two cases we see that when ε > 0 and δ > 0 are both sufficiently small, 0 / ∈ F (∂S). Now, we will check whether the conditions (3). In fact, since A is an irreducibly diagonally dominant Mmatrix, we have a ii > 0, a l,j ≤ 0 for l = j, and a ii − i =j |a ij | > 0, for i, j ∈ I. 1 Let G = {x ∈ R n : |x i | < L, for 1 ≤ i ≤ n} and suppose the mapping F = (f 1 , · · · , fn) :Ḡ → Rn is continuous on the closureḠ of G such that F (x) = 0 for x on the boundary ∂G of G, and Then, F (x) = 0 has a solution in G. See [1].
Hence, combing this property and x ∞ ≤ C, we have for 1 ≤ i ≤ n when C is sufficiently large. In the same way, we also have that when C is sufficiently large, for 1 ≤ i ≤ n From the above analysis, we see that all the conditions of Miranda's theorem are satisfied. Hence, the existence of the solution to the penalized Problem 2 is proved.

Convergence analysis.
In what follows, we show that the solution to (2) satisfies Problem 1 as λ → ∞, and establish the rate of the convergence of the power penalty approach.
3.1. Convergence. For the penalty approach (2), we have the following error estimate.
Theorem 3.1. Let x λ be the solution to Problem 2 for λ > 1. There exists a constant C > 0, independent of λ, such that We continue to use the set I 1 , I 2 and I 3 defined in (4). Now, we consider the following three cases.
• ∀i ∈ I 1 , it follows from the penalty equation (2) and the definition I 1 that Hence, • ∀i ∈ I 2 , it follows from the penalty equation (2) and the definition I 2 that Hence, • ∀i ∈ I 3 , it follows from the penalty equation (2) and the definition I 3 that Hence, Combining (7)- (9), we obtain From Lemma 2.1, we see that x λ is bounded for any λ > 0. Thus, the above inequality implies

3.2.
Rate of convergence. With the above error estimation, we are now ready to show that the solution of Problem 2 converges to that of Problem 1 exponentially w.r.t. the penalty parameter. This is given in the following theorem.

Theorem 3.2. Let x *
λ and x * be the solutions to Problems 2 and 1 respectively. When λ is sufficiently large, we have where C is a positive constant, independent of λ.
Proof. Since x * λ and x * are the solutions to Problems 2 and 1, respectively, they satisfy and It follows from (11) and f < g that for each i ∈ I Then for each i ∈ I, we consider i in the following three situations. 1. When i ∈ J 1 := {i ∈ I : (x * λ ) i < f i }, (13) implies that (11) and (12) . Now, we again introduce a matrix, denoted as A * 1 , so that (A * 1 ) i = (I) i when i ∈ J 1 ∪ J 2 , and (A * 1 ) i = (A) i when i ∈ J 3 , where I denotes the N × N identity matrix. From the above three cases we see that where e = (1, . . . , 1) . From its construction and the fact that A is an M -matrix we see that A * 1 is also an M -matrix. Therefore, from the above estimate we have The above inequality defines an upper bound for x * λ −x * . To define a lower bound for it, such as above proof, we divide I into following three sets: . Similar to the proof of (14), we obtain that Combining (14) and (15), we obtain (10) for some constant C > 0 independent of λ.
Remark 1. It is worth noting that a monotonic convergence property for the solution sequence {x λm } w.r.t. the penalty parameter λ m is usually held when the power penalty method is used to approach to the one-side/unilateral obstacle problems, see [12,14,15,17], etc. However, for the double obstacle problems, due to the existence of two-side/bilateral obstacles, the monotonic convergence property is not held any longer. This conclusion is numerically illustrated in Example 1.

Solution of the penalty equation.
In this section we will develop an iterative method to numerically solve (2). Clearly, when k ≥ 1, the power penalty equation (2) becomes nonlinear and nonsmooth, which make the classic Newton method inapplicable. To overcome this difficulty, we utilize the smoothing technique, developed in [14], to smoothing out the non-smooth function [z] 1 is a regularization parameter. Likewise, the same technique can be used for the penalty term [z] With this smoothing technique, the equation (2) becomes where P (x) and Q(x) are two vectors defined by for i = 1, · · · , N . We then apply the classic Newton method to solve (16), which yields the following algorithm. For clarity, we omit the subscript λ of x λ in the algorithm.
Step 1.: Choose , ε > 0 sufficiently small; Set l = 0 and an initial guess x 0 such that f ≤ x 0 ≤ g. Step 2.: Solve the following linear system for p l+1 : where J P (x l ) and J Q (x l ) are the Jacobian matrices of P (x) and Q(x), respectively, defined by Step 3.: Set x l+1 = x l + νp l+1 , where 0 < ν < 1 is a damping parameter determined by the Armijo linear search method [5].
Step 4.: The exact solution is x * = (1, 0, 0, 5) . The power penalty approach to this example is 1/k + = 0. To compute the rate of convergence of the power penalty method, we find approximations to the solution to this discrete double obstacle problem by solving the above nonlinear system with k = 1 and 2 using Algorithm 1. To calculate the rates of convergence, we choose λ i = 10 i for i = 2, 3, 4, 5, and the l ∞ -norms of the errors between the numerical solutions and the exact solution are calculated. The computed errors, along with the ratios of errors from two consecutive values of λ, are listed in Table 1. Table 1. Results computed by the power penalty method. Tolerance ε = 10 −6 is chosen for the Newton iteration. The smoothing parameter = 10 −3 . 'Ra' stands for ratio. From the table we see that the numerical rates of convergence are close to O(λ −k ), in consistence with the theoretical result in (10), though some in k = 2 are less than 2. These reduced rates of convergence may due to the relatively large value of . From the error x − x λi ∞ in Table 1 we see that to achieve the same accuracy, the lower order penalty method (k = 2) requires a much smaller value of λ than linear penalty method (k = 1) needs. This verifies the advantage of lower order penalty method. Interestingly, though the system matrix A is not an M -matrix, the rates of convergence are still close to the theoretical one in Theorem 3.2.
As stated in Remark 1, the monotonic convergence property for the whole solution sequence {x λi } w.r.t. the penalty parameter λ is not held any longer. This is clearly observed in Table 1  We then use the following example to explore the advantages of the power penalty method over the policy iteration method.
Example 2. Consider the discrete double-obstacle problem: find (U i ) 1≤i≤N in R N such that, for i = 1, · · · , N, This example is from [2], where it has been shown that the classic policy iteration method failed to solve this problem. Though a modified version of the policy iteration method is proposed to overcome this difficulty therein, the computational cost in terms of total number of iterations is very high. We list all the results computed by the power penalty method (in both cases of k = 1 and k = 2), and the modified policy iteration method in Table 2. The results demonstrate that the power penalty methods can solve this problem effectively and possesses advantages over the modified policy iteration method, since the modified policy iteration method need much more computational costs than the power penalty method. Table 2. Total number of iterations for the linear penalty (k = 1), lower order penalty (k = 2), and modified policy methods with N = 99. Tolerance is set to be 10 −6 . = 10 −3 . λ is set to be 10 6 and 10 3 for k = 1 and k = 2, respectively.
Modified policy iter. Iterations 9 12 88 To illustrate the usefulness of the power penalty method, we plot the computed solution in Figure 1. For reference, the the upper bound (obstacle g) and the lower bound (obstacle f ) are also plotted. From this figure, we see that both the upper and lower bounds are satisfied by our solution.
To further explore the computational efficiency of the power penalty approach, in the following Example 3, we carry out a numerical performance comparison study for the penalty methods, the Gauss-Seidel method and the active set method, since the latter two methods are commonly used to solve the discrete obstacle problem. Example 3. Consider the following two-dimensional bilateral obstacle problem: 6 < x < 1. This example had been carefully examined in the literature (cf. [8]). With the same level of accuracy (ε = 10 −6 ), we list all the average numbers of iterations and computation times on two different mesh grids: 50 × 50 and 60 × 60 in Table 3. The table clearly shows that the computational efficiency of the power penalty methods are rather superior to that of the Gauss-Seidel method, and comparable to that of the active set method. We also note that when the discretization of the continuous problem is not fine enough, i.e. the mesh grid is coarse, the active set method fails. This is a main drawback of the active set method, which has been verified in [8].
Nevertheless, the power penalty method is very robust to the dscretization, which is another advantage of the power penalty method, as we stated in the introduction.
Finally, we plot the solution u(x, y) and active sets in Figures 2 and 3, respectively. These results are computed with the lower penalty method (k = 2) on mesh grid 60 × 60. In Figure 3 the lower and upper coincidence sets are marked with dot '·' and star ' * ' on the mesh grid, respectively. These figures are consistent with those in [8], which again illustrates that the power penalty method is an effective method. Table 3. Comparison of computational efficiency among power penalty methods, Gauss-Seidel iteration method and active set method. Tolerance is set to be 10 −6 . = 10 −3 . λ is set to be 10 6 and 10 3 for k = 1 and k = 2, respectively.  . Solution u(x, y) with the mesh grid 60 × 60, obtained with the lower order penalty method (k = 2). Tolerance is set to be 10 −6 . = 10 −3 . λ = 10 3 .

6.
Conclusions. In this work we developed a penalty-based method to the solution of the discrete double obstacle problem. By penalizing the two obstacles we transform the obstacle problem into a nonlinear system containing two penalty terms. After showing the unique solvability of the nonlinear system, an exponential convergence rate of the power penalty approach w.r.t. the penalty parameter was