Accurate and efficient numerical solutions for elliptic obstacle problems

Elliptic obstacle problems are formulated to find either superharmonic solutions or minimal surfaces that lie on or over the obstacles, by incorporating inequality constraints. In order to solve such problems effectively using finite difference (FD) methods, the article investigates simple iterative algorithms based on the successive over-relaxation (SOR) method. It introduces subgrid FD methods to reduce the accuracy deterioration occurring near the free boundary when the mesh grid does not match with the free boundary. For nonlinear obstacle problems, a method of gradient-weighting is introduced to solve the problem more conveniently and efficiently. The iterative algorithm is analyzed for convergence for both linear and nonlinear obstacle problems. An effective strategy is also suggested to find the optimal relaxation parameter. It has been numerically verified that the resulting obstacle SOR iteration with the optimal parameter converges about one order faster than state-of-the-art methods and the subgrid FD methods reduce numerical errors by one order of magnitude, for most cases. Various numerical examples are given to verify the claim.


Introduction
Variational inequalities have been extensively studied as one of key issues in calculus of variations and in the applied sciences. The basic prototype of such inequalities is represented by the so-called obstacle problem, in which a minimization problem is often solved. The obstacle problem is, for example, to find the equilibrium position u of an elastic membrane whose boundary is held fixed, with an added constraint that the membrane lies above a given obstacle ϕ in the interior of the domain Ω ⊂ R d : where Γ = ∂Ω denotes the boundary of Ω and f is the fixed value of u on the boundary. The problem is deeply related to the study of minimal surfaces and the capacity of a set in potential theory as well. Other classical applications of the obstacle problem include the study of fluid filtration in porous media, constrained heating, elasto-plasticity, optimal control, financial mathematics, and surface reconstruction [-].
The problem in (.) can be linearized in the case of small perturbations by expanding the energy functional in terms of its Taylor series and taking the first term, in which case the energy to be minimized is the standard Dirichlet energy A variational argument [] shows that, away from the contact set {x|u(x) = ϕ(x)}, the solution to the obstacle problem (.) is harmonic. A similar argument (which restricts itself to variations that are positive) shows that the solution is superharmonic on the contact set. Thus both arguments imply that the solution is a superharmonic function. As a matter of fact, it follows from an application of the maximum principle that the solution to the obstacle problem (.) is the least superharmonic function in the set of admissible functions. The Euler-Lagrange equation for (.) reads In modern computational mathematics and engineering, the obstacle problems are not extremely difficult to solve numerically any more, as shown in numerous publications; see [-], for example. However, most of those known methods are either computationally expensive or yet to be improved for higher accuracy and efficiency of the numerical solution. In this article, we consider accuracy-efficiency issues and their remedies for the numerical solution of elliptic obstacle problems. This article makes the following contributions.
-Accuracy improvement through subgrid finite differencing of the free boundary: It can be verified either numerically or theoretically that the numerical solution easily involve a large error near the free boundary (the edges of obstacles), particularly when the grid mesh does not match with the obstacle edges. We suggest a post-processing algorithm which can reduce the error (by about a digit) by detecting accurate free boundary in subgrid level and introducing nonuniform finite difference (FD) method. The main goal of the subgrid FD algorithm is to produce a numerical solution of a higher accuracy u h , which guarantees u h (x) ≥ ϕ(x) for all points x ∈ Ω. -Obstacle SOR: The iterative algorithm for solving the linear system of the obstacle problem is implemented based on one of simplest iterative algorithms, the successive over-relaxation (SOR) method. Convergence of the obstacle SOR method is analyzed and compared with modern sophisticated methods. We also suggest an effective way to set the optimal relaxation parameter ω. Our simple obstacle SOR method with the optimal parameter performs better than state-of-the-art methods in both accuracy and efficiency. -Effective numerical methods for nonlinear problems: For the nonlinear obstacle problem (.), a method of gradient-weighting is introduced to solve the problem more conveniently and efficiently. In particular, the suggested numerical schemes for the gradient-weighting problem produce an algebraic system of a symmetric and diagonally dominant M-matrix of which the main diagonal entries are all the same positive constant. Thus the resulting system is easy to implement and presumably converges fast; as one can see from Section , the obstacle SOR algorithm for nonlinear problems converges in a similar number of iterations as for linear problems. The article is organized as follows. The next section presents a brief review for stateof-the-art methods for elliptic obstacle problems focusing the one in []. Also, accuracy deterioration of the numerical solution (underestimation) is discussed by exemplifying an obstacle problem in D where the mesh grid does not match with edges of the free boundary. In Section , the SOR is applied for both linear and nonlinear problems and analyzed for convergence; the limits of iterates are proved to satisfy discrete obstacle problems. A method of gradient-weighting and second-order FD schemes are introduced for nonlinear problems. An effective strategy is suggested to find the optimal relaxation parameter. Section  introduces subgrid FD schemes near the free boundary in order to reduce accuracy deterioration of the numerical solution. In Section , various numerical examples are included to verify the claims we just made. Section  concludes the article summarizing our experiments and findings.

Preliminaries
As preliminaries, we first present a brief review for state-of-the-art methods for elliptic obstacle problems and certain accuracy issues related to the free boundary.

State-of-the-art methods for elliptic obstacle problems
This subsection summarizes state-of-the-art methods for elliptic obstacle problems focusing on the primal-dual method incorporating L  -like penalty term (PDLP) studied by Zosso et al. []. Primal-dual splitting methods have a great deal of attention, particularly in the context of total variation (TV) minimization and L  -type problems in image processing [-].
In the literature of optimization problems, one of common practices is to reformulate a constrained optimization problem for a unconstrained problem by incorporating the constraint as a penalty term. Recently, Tran et al. [] proposed the following minimization problem of a L  -like penalty term: where μ is a Lagrange multiplier and (·) + = max(·, ). It is shown that, for sufficiently large but finite μ, the minimizer of the unconstrained problem (.) is also the minimizer of the original, constrained problem (.). The PDLP [] is a hybrid method which combines primal-dual splitting algorithm and the L  -like penalty method in (.); it can be summarized as follows.
where ∇ h denotes the numerical approximation of the gradient ∇, associated with the mesh size h, r  and r  are constants to be determined, p n is the dual variable representing the gradient of the primal variable (u n ), and u * is an intermediate solution.
Here P ϕ is an obstacle projection defined by (.) The above algorithm can be implemented effectively. It follows from (.)(a) that where h is the discrete Laplacian. Thus S n+ ≡ ∇ h · p n+ can be considered as a variable and updated in each iteration, averaging its previous iterate S n and h u n as in (.). As analyzed in [], PDLP (away from the obstacle) can be compared to either the forward Euler (explicit) scheme for discrete heat equation or a three-level time stepping method for a damped acoustic wave equation, where r  r  plays the role of the time-step size. PDLP converges when where h is the operator/induced norm of the discrete Laplacian h (= , when the mesh size h = ). The authors in [] claimed that '[Their] results achieve state-of-the-art precision in much shorter time; the speed up is one-two orders of magnitude with respect to the method in [], and even larger compared to older methods [-]. ' Thus, in this article our suggested method would be compared mainly with PDLP (the best-known method), in order to show its superiority.

Accuracy issues
The solution of obstacle problems must lie on or over the obstacle (u ≥ ϕ), which is also one of requirements for numerical solutions. For FD methods and finite element (FE) methods for the obstacle problem (.), for example, this requirement can easily be violated when edges of the free boundary does not match with mesh grids. See Figure , where the shaded rectangle indicates the obstacle defined on one-dimensional (D) interval [x  , x  ]: Let C h denote the numerical contact set: where Ω  h is the set of interior grid points. Define an interior grid point is a neighboring point if it is not in the contact set but one of its adjacent grid points is in the contact set. Let the set of neighboring points be called the neighboring set N h . Then, for the example in Figure  The accuracy of the numerical solution u h can be improved by applying a postprocessing in which a subgrid FD method is applied at grid points in the neighboring set. For example, at x = x  , -u xx can be approximated by employing nonuniform FD schemes over the grid points [x  , x  , p], given as where r = (px  )/h x ∈ (, ], and therefore numerical solution of -u xx =  at x = x  must satisfy As r is approaching  (i.e., (px  ) becomes smaller proportionally), the obstacle value ϕ(p) is more weighted. On the other hand, when r = , ϕ(p) = u  and the scheme in (.) becomes the standard second-order FD scheme. Let the numerical solution u be obtained from where u  =  and u  = . Then it is not difficult to prove that u is exactly the same as the true solution u at all grid points (except numerical rounding error), regardless of the grid size h x . The above example has motivated the authors to develop an effective numerical algorithm for elliptic obstacle problems in D which detects the neighboring set of the free boundary, determines the subgrid proportions (r's), and updates the solution for an improved accuracy using subgrid FD schemes. Here the main goal is to try to guarantee u(x) ≥ ϕ(x) for all x ∈ Ω (whether x is a grid point or not). Since it is often the case that the free boundary is determined only after solving the problem, the algorithm must be a post-process. Details are presented in Section .

Obstacle relaxation methods
This section introduces and analyzes effective relaxation methods for solving (.) and its nonlinear problem as shown in (.) below.

The linear obstacle problem
We will begin with second-order approximation schemes foru. For simplicity, we consider a rectangular domain in R  , Ω = (a x , b x ) × (a y , b y ). Then the following second-order FD scheme can be formulated on the grid points: where, for some positive integers n x and n y , Let u pq = u(x p , y q ). Then, at each of the interior points x pq , the five-point FD approximation ofu reads Multiply both sides of (.) by h  x to have where r xy = h x /h y and u st = f st at boundary grid points (x s , y t ). Now, consider the following Jacobi iteration for simplicity. Given an initialization u  , find u n iteratively as follows.
where u n- st = f st at boundary grid points (x s , y t ).
Note that Algorithm L J produces a solution u of which the function value at a point is a simple average of four neighboring values, satisfying the constraint u ≥ ϕ.
Theorem  Let u be the limit of the iterates u n of Algorithm L J . Then u satisfies the FD discretization of (.). That is, where Ω  h denotes the set of interior grid points and Γ h is the set of boundary grid points.
Proof It is clear to see from Algorithm L J that Let u pq = ϕ pq at an interior point (x p , y q ). Then it follows from (.)(b) that On the other hand, let u pq > ϕ pq at (x p , y q ). Then, since u pq = max( u J,pq , ϕ pq ), we must have which implies thath u pq = . This completes the proof.
One can easily prove that the algebraic system obtained from (.) is irreducibly diagonally dominant and symmetric positive definite. Since its off-diagonal entries are all nonpositive, the matrix must be a Stieltjes matrix and therefore an M-matrix []. Thus relaxation methods of regular splittings (such as the Jacobi, the Gauss-Seidel (GS), and the successive over-relaxation (SOR) iterations) are all convergent and their limits are the same as u and therefore satisfy (.). In this article, variants of Algorithm L J for the GS and the SOR would be denoted, respectively, by L GS and L SOR (ω), where ω is an overrelaxation parameter for the SOR,  < ω < . For example, L SOR (ω) is formulated as Note that the right side of (.)(a) involves updated values wherever available. When ω = , Algorithm L SOR (ω) becomes Algorithm L GS ; that is, L SOR () = L GS .

The nonlinear obstacle problem
Applying the same arguments for the linear problem (.), the Euler-Lagrange equation for the nonlinear minimization problem (.) can be formulated as Thus the solution to the nonlinear problem (.) can be considered as a minimal surface satisfying the constraint given by the obstacle function ϕ.
Since  + |∇u|  ≥ , the nonlinear obstacle problem (.) can equivalently be formulated as Such a method of gradient-weighting will make algebraic systems simpler and better conditioned, as to be seen below. In order to introduce effective FD schemes for M(u), we where both (  + |∇u|  )  and (  + |∇u|  )  are the same as  + |∇u|  ; however, they will be approximated in a slightly different way. The following numerical schemes are of second-order accuracy and specifically designed for the resulting algebraic system to be simpler and better conditioned.
For the FD scheme at the (p, q)th pixel, we first compute second-order FD approximations of  + |∇u|  at x p-/,q (W ), x p+/,q (E), x p,q-/ (S), and x p,q+/ (N): Then the directional-derivative terms at the pixel point x pq can be approximated by Now, we discretize the surface element as follows: where the right-hand sides are harmonic averages of FD approximations of  + |∇u|  in x-and y-coordinate directions, respectively. Then it follows from (.), (.), and (.) that Note that a pq,W + a pq,E = a pq,S + a pq,N = . As for the linear problem, it is easy to prove that the algebraic system obtained from (.) is an M-matrix. Given FD schemes for M(u) as in (.), the nonlinear obstacle problem (.) can be solved iteratively by the Jacobi iteration.
(a n- pq,W u n- p-,q + a n- pq,E u n- p+,q + r  xy a n- pq,S u n- p,q- + r  xy a n- pq,N u n- p,q+ ); (b) u n pq = max(u J,pq , ϕ pq ); end end end where u n- st = f st at boundary grid points (x s , y t ).
The superscript (n -) on the coefficients a pq,D , D = W , E, S, N , indicate that they are obtained using the last iterate u n- . Algorithm N J produces a solution u of which the function value at a point is a weighted average of four neighboring values, satisfying the constraint u ≥ ϕ. One can prove the following corollary, using the same arguments introduced in the proof of Theorem .
Corollary  Let u be the limit of the iterates u n of Algorithm N J . Then u satisfies the FD discretization of (.). That is, where M h ( u) pq denotes the FD scheme of M(u)(x pq ) as defined in (.) with u = u.
Variants of Algorithm N J for the GS and the SOR can be formulated similarly as for the linear obstacle problem; they would be denoted respectively by N GS and N SOR (ω). In practice, such symmetric coercive optimization problems, the SOR methods are much more efficient than the Jacobi and Gauss-Seidel methods. We will exploit L SOR (ω) and N SOR (ω) for numerical comparisons with state-of-the-art methods, by setting the relaxation parameter ω optimal.

The optimal relaxation parameter ω
Consider the standard Poisson equation with a Dirichlet boundary condition for prescribed functions f and g. Let Ω = [, ]  , for simplicity, and apply the second-order FD method for the second derivatives on a uniform grid: h = h x = h y = /(m + ), for some positive integer. The its algebraic system can be written as Then the theoretically optimal relaxation parameter for the SOR method can be determined as [], Section ., where ρ(T J ) is the spectral radius of the iteration matrix of the Jacobi method T J . The iteration matrix T J can be explicitly presented as a block tridiagonal matrix where I m is the m-dimensional identity matrix and For such a matrix T J , it is well known that Thus it follows from (.) and (.) that the optimal SOR parameter corresponding to the mesh size h, ω h , can be expressed as where c  = √ c. Hence, for general mesh size h, the corresponding optimal SOR parameter ω h can be found as follows.
It is often the case that the calibration (.)(a)-(.)(b) can be carried out with a small problem, i.e., with h  of a very low resolution.

Subgrid FD schemes for the free boundary: a post-process
This section describes subgrid FD schemes for the free boundary, focusing on the linear obstacle problem; the arguments to be presented can be applied the same way for nonlinear problems. Again, we assume for simplicity that h = h x = h y .
Let u be the numerical solution of an obstacle problem. Then it would satisfy the discrete obstacle problem (.), particularly u pq ≥ ϕ pq at all (interior) grid points x pq ∈ Ω  h .
However, when the mesh grid is not matching with the free boundary, the obstacle constraint u ≥ ϕ may not be satisfied at all points x ∈ Ω. This implies that when the mesh is not fine enough, the numerical solution can be underestimated near the free boundary, as shown in Figure  in Section .. Note that the error introduced by non-matching grids is in O(h), while the numerical truncation error is in O(h  ) for second-order FD schemes. That is, the underestimation is in O(h), which can be much larger than the truncation error. The strategy below can be considered as a post-processing algorithm designed in order to reduce the underestimation without introducing a mesh refinement. The postprocessing algorithm consists of three steps: (a) finding the numerical contact set and the neighboring set, (b) subgrid determination of the free boundary, and (c) nonuniform FD schemes on the neighboring set.

The contact set and the neighboring set
Finding the numerical contact set is an easy task. Let u and ϕ be the numerical solution and the lower obstacle, respectively. Then, for example, the characteristic set of contact points C h can be determined as follows.
As defined in Section ., an interior grid point is a neighboring point when it is not in the contact set but one of its adjacent grid points is in the contact set. Thus the neighboring points can be found more effectively as follows. Visit each point in the contact set; if any one of its four adjacent points is not in the contact set, then the non-contacting point is a neighboring point. The set of all neighboring points is the neighboring set N h .

Subgrid determination of the free boundary
Let x pq be a neighboring point with two of its adjacent points are contact points (C h (p + , q) = C h (p, q -) = ), as in Figure . Then we may assume that the real free boundary passes somewhere between the contact points and the neighboring points. We will suggest an effective strategy for the determination of the free boundary in subgrid level.
We first focus on the horizontal line segment connecting x pq and x p+,q in the east (E) direction. Define (.) Then the corresponding linear interpolation between u pq and u p+,q over the line segment is formulated as Since x pq and x p+,q are a neighboring point and a contact point, respectively, we have If the free boundary passes between x pq and x p+,q , then there must exist r ∈ (, ) such that F E (r) > . Let r E be such that x E (r E ) represents the intersection between the line segment x E (·) and the free boundary. Then it can be approximated as follows.
The maximization problem in (.) can be solved easily (using the Newton method, for example), when the obstacle is defined as a smooth function. A more robust method can be formulated as a combination of a line search algorithm and the bisection method.

Remarks
-The last evaluation of ϕ (and saving) is necessary for the nonuniform FD schemes on the neighboring set, which will be discussed in Section .. The quantity B E will be used as the Dirichlet value on the free boundary.
where h is mesh size. It has been numerically verified that the choice (k  , k  ) = (, ) is enough for an accurate detection of the free boundary, for which the upper bound of the error becomes h/ = .h.

Nonuniform FD schemes on the neighboring set
Let x pq = (x p , y q ) be a neighboring point. Then r pq,D ∈ (, ] would be available for each D ∈ {W , E, S, N}; B pq,D is also available for r D < . Thus the FD scheme for -u xx (x pq ) can be formulated over three points {(x pr pq,W h x , y q ), (x p , y q ), (x p + r pq,E h x , y q )} as follows.
Similarly, the FD scheme for -u yy (x pq ) can be formulated over three points in the ydirection {(x p , y qr pp,S h y ), (x p , y q ), (x p , y q + r pp,N h y )}: u pq,S r pq,S (r pq,S + r pq,N ) + u pq r pq,S · r pq,N -u pq,N r pq,N (r pq,S + r pq,N ) , Thus, the post-processing algorithm of the obstacle SOR (.), L SOR (ω), can be formulated by replacing the two terms in the right side of (.) with the right sides of (.) and (.), and computing u GS,pq in (..a) correspondingly at all neighboring points.

Numerical experiments
In this section, we apply the obstacle SOR method and the post-processing schemes to various obstacles to verify their effectiveness and accuracy. We mainly concern -D obstacle problems of Dirichlet boundary conditions. The algorithms are implemented, for both one and double obstacles, in Matlab and carried out on a Desktop computer of an Intel i-S . GHz processor. The optimal relaxation parameter is calibrated with the lowest resolution to find a constant c  (.) and the constant is used for all other cases. For a comparison purpose, we implemented a state-of-the-art method, PDLP [], and its parameters (r  and r  in (.)) are found heuristically for cases where the parameters are not suggested in []. The iterations are stopped when the maximum difference of consecutive iterates becomes smaller than the tolerance ε: where ε =  - mostly; Section . uses ε =  - for an accurate estimation of the error.
For all examples, the numerical solution is initialized from ϕ (the lower obstacle) and the boundary condition f . (.)

Linear obstacle problems
We first consider a non-smooth obstacle ϕ  :   As the second example, we consider the radially symmetric obstacle ϕ  : where r * = . . . . , the solution of For the obstacle ϕ  , the analytic solution to the linear obstacle problem can be defined as when the boundary condition is set appropriately using u * . See Figure , in which we give plots of ϕ  and the true solution u * .
In Table , we compare performances of the PDLP and the obstacle SOR applied for the linear obstacle problem with (.). The PDLP uses the parameters suggested in [] (μ = ., r  = ., r  = .). As one can see from the table, our suggested method takes about one order less CPU time than the PDLP for the computation of the numerical solution. In Figure , we show the numerical solutions u h and the errors u hu * produced  The PDL1P uses the parameters suggested in [14] (μ = 0.1, r 1 = 0.008, r 2 = 15.625).
by the PDLP and the obstacle SOR at the  ×  resolution. The solutions are almost identical and the errors are nonpositive. This implies that the numerical solutions of the obstacle problem are underestimated. As a more general obstacle problem, we consider the elastic-plastic torsion problem in []. The problem is to find the equilibrium position of the membrane between two obstacles ϕ, ψ that a force v is acting on: Let Ω = [, ]  and the problem consist of two obstacles ϕ  : Ω → R, ψ  : Ω → R and the force v : (.)  In Table , we present performances of the PDLP and the obstacle SOR applied for the elastic-plastic torsion problem (.). For the PDLP, we use the parameters suggested in [] (μ = ., r  = ., r  = .). As one can see from the table, our suggested method again resulted in the numerical solution about one order faster than the PDLP measured in the CPU time. In Figure , we illustrate the simulated membranes in the equilibrium satisfying (.) and their contact sets at resolution  × . In Figures (b)  For the PDL1P, we use the parameters suggested in [14] (μ = 0.1, r 1 = 0.008, r 2 = 15.625). and (d), the upper and lower contact sets are colored in yellow (brightest in gray scale) and blue (darkest in gray scale), respectively. The results produced by the two methods are apparently the same.

Nonlinear obstacle problems
The obstacle SOR is implemented for nonlinear obstacle problems as described in Section ..
In Table , we present experiments for which the obstacle SOR is applied for nonlinear obstacle problems with ϕ = ϕ i , i = , , . From a comparison with linear cases presented in Tables , , and , we can see for each of the obstacles that the obstacle SOR iteration for the nonlinear problem converges in a similar number of iterations as for the linear problem.   Only the apparent difference is the CPU time; an iteration of the nonlinear solver is about as six time expensive as that of the linear solver, due to the computation of coefficients as in (.). For ϕ = ϕ  , the nonlinear solution is plotted in Figure . Compared with the linear solutions in Figure , the nonlinear solution shows slightly lower function values, which is expected. As the grid point approaches the obstacles, the solution shows an increasing gradient magnitude. This may enlarge weights for far-away grid values as shown in (.), which in return acts as a force to reduce function values. The difference between the linear solution and the nonlinear solution, at the  ×  resolution, is depicted in Figure .

Post-processing algorithm
In Figure , one have seen that the error, the difference between the numerical solution and the analytic solution, shows its highest values near the free boundary. The larger error is due to the result of mismatch between the mesh grid and the obstacle edges. In order to eliminate the error effectively, we apply the subgrid FD schemes in Section  as a postprocessing (PP) algorithm.
For the examples presented in this subsection, the numerical solutions are solved as follows: (a) the problem is solved with ε =  - (pre-processing), (b) the free boundary is estimated with (k  , k  ) = (, ) and subgrid FD schemes are applied at neighboring grid points as in Section , and (c) another round of iterations is applied to satisfy the tolerance ε =  - . First, we consider a step function for an one-dimensional (D) obstacle, as in Section .. Let Ω = [, ] and ϕ  : Ω → R defined by The analytic solution to the linear problem is given as Figure  shows the numerical solutions to the linear problem associated to (.) with and without the post-process, and their errors. The numerical solutions without and with the post-process are obtained iteratively satisfying the tolerance ε =  - . Notice that the solution without post-process is underestimated and shows a relatively high error: uu ,true ∞ = .. The error is reduced to u ppu ,true ∞ = . ×  - after the post-process.
The post-processing algorithm is applied to the linear obstacle problem in -D involving ϕ = ϕ  . Table  contains efficiency results that compare performances of the PDLP, the obstacle SOR (without post-process), and the obstacle SOR with the post-process    (Obstacle SOR+PP) at various resolutions; while Table  presents an accuracy comparison for those methods. According to Table , the post-processed solution requires about % more iterations than the non-post-processed one; the incorporation of the postprocess makes the iterative algorithm as twice expensive measured in CPU time as the original iteration. However, one can see from Table  that the post-process makes the error reduced by a factor of ∼. Thus in order to achieve a three-digit accuracy in the maximum-norm, for example, the PDLP requires . seconds and the obstacle SOR completes the task in . seconds; when the obstacle SOR+PP takes only . second. Figure  includes plots of the error (u hu * ) at the  ×  resolution for the linear obstacle problem with ϕ = ϕ  , produced by the PDLP, the obstacle SOR, and the obstacle SOR+PP. The numerical solutions of the PDLP and the obstacle SOR are almost identical to each other and clearly underestimated, with the maximum discrepancy occurring around the free boundary due to the misfit between the mesh grid and the free boundary. It can be seen from Figure (c) that the post-process can eliminate the misfit error very effectively; the remaining error is the truncation error introduced by the second-order FD schemes.

Parameter choices
Finally, we present experimental results for parameter choices, when the obstacle SOR is applied for the linear problem with ϕ = ϕ  . For an effective calibration of the optimal relaxation parameter as suggested in (.), we first select h  = /. Then by using a trialby-error method, we found the calibrated optimal relaxation parameter ω h  = ., which results in the following calibrated constant: which is used for the results of the obstacle SOR included in Table .
In order to verify effectiveness of the calibration, we implement a line search algorithm to find a relaxation parameter ω that converges in the smallest number of iterations with ε =  - , the same tolerance as for the results in Table  Note that when the calibrated parameters are used, the iteration counts of the obstacle SOR presented in Table  are , , and , respectively, for h = /, h = /, and h = /. Thus the calibrated optimal parameters in (.) are quite accurate for the optimal convergence.

Conclusions
Although various numerical algorithms have been suggested for solving elliptic obstacle problems effectively, most of the algorithms presented in the literature are yet to be improved for both accuracy and efficiency. In this article, the authors have studied obstacle relaxation methods in order to get second-order finite difference (FD) solutions of obstacle problems more accurately and more efficiently. The suggested iterative algorithm is based on one of the simplest relaxation methods, the successive over-relaxation (SOR). The iterative algorithm is incorporated with subgrid FD methods to reduce accuracy deterioration occurring near the free boundary when the mesh grid does not match with the free boundary. For nonlinear obstacle problems, a method of gradient-weighting has been introduced to solve the problem more conveniently and efficiently. The iterative algorithm has been analyzed for convergence for both linear and nonlinear obstacle problems. An effective strategy is also presented to find the optimal relaxation parameter. The resulting obstacle SOR has converged about one order faster than state-of-the-art methods and the subgrid FD methods could reduce the numerical errors by one order of magnitude.