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 : 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 [1-7].

The problem in (1.1) 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 [2] shows that, away from the contact set , the solution to the obstacle problem (1.2) 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 (1.2) is the least superharmonic function in the set of admissible functions. The Euler-Lagrange equation for (1.2) 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 [8-14], 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 , which guarantees for all points .

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 (1.1), 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 5, 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 state-of-the-art methods for elliptic obstacle problems focusing the one in [14]. Also, accuracy deterioration of the numerical solution (underestimation) is discussed by exemplifying an obstacle problem in 1D where the mesh grid does not match with edges of the free boundary. In Section 3, 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 4 introduces subgrid FD schemes near the free boundary in order to reduce accuracy deterioration of the numerical solution. In Section 5, various numerical examples are included to verify the claims we just made. Section 6 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 -like penalty term (PDL1P) studied by Zosso et al. [14]. Primal-dual splitting methods have a great deal of attention, particularly in the context of total variation (TV) minimization and -type problems in image processing [15-19].

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. [13] proposed the following minimization problem of a -like penalty term: where μ is a Lagrange multiplier and . It is shown that, for sufficiently large but finite μ, the minimizer of the unconstrained problem (2.1) is also the minimizer of the original, constrained problem (1.2).

The PDL1P [14] is a hybrid method which combines primal-dual splitting algorithm and the -like penalty method in (2.1); it can be summarized as follows. where denotes the numerical approximation of the gradient ∇, associated with the mesh size h, and are constants to be determined, is the dual variable representing the gradient of the primal variable (), and is an intermediate solution. Here is an obstacle projection defined by

The above algorithm can be implemented effectively. It follows from (2.2)(a) that where is the discrete Laplacian. Thus can be considered as a variable and updated in each iteration, averaging its previous iterate and as in (2.4). As analyzed in [14], PDL1P (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 plays the role of the time-step size. PDL1P converges when where is the operator/induced norm of the discrete Laplacian (=8, when the mesh size ). The authors in [14] 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 [13], and even larger compared to older methods [20-22].’ Thus, in this article our suggested method would be compared mainly with PDL1P (the best-known method), in order to show its superiority.

Accuracy issues

The solution of obstacle problems must lie on or over the obstacle (), which is also one of requirements for numerical solutions. For FD methods and finite element (FE) methods for the obstacle problem (1.3), for example, this requirement can easily be violated when edges of the free boundary does not match with mesh grids. See Figure 1, where the shaded rectangle indicates the obstacle defined on one-dimensional (1D) interval : which is not matching with the mesh grids . The figure shows the true solution u (red solid curve) and a numerical solution (blue dashed curve) of the linear obstacle problem (1.3) in 1D. The numerical solution is clearly underestimated and the magnitude of the error is maximized at : which is .

Figure 1

A non-matching grid: The true solution (red solid curve) and the numerical solution on the non-matching grid (blue dashed curve). Here the obstacle is the shaded region, , and .

Let denote the numerical contact set: where 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 . Then, for the example in Figure 1, and .

The accuracy of the numerical solution can be improved by applying a post-processing in which a subgrid FD method is applied at grid points in the neighboring set. For example, at , can be approximated by employing nonuniform FD schemes over the grid points , given as where , and therefore numerical solution of at must satisfy As r is approaching 0 (i.e., becomes smaller proportionally), the obstacle value is more weighted. On the other hand, when , and the scheme in (2.9) becomes the standard second-order FD scheme. Let the numerical solution ũ be obtained from where and . Then it is not difficult to prove that ũ is exactly the same as the true solution u at all grid points (except numerical rounding error), regardless of the grid size .

The above example has motivated the authors to develop an effective numerical algorithm for elliptic obstacle problems in 2D 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 for all (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 4.

Obstacle relaxation methods

This section introduces and analyzes effective relaxation methods for solving (1.3) and its nonlinear problem as shown in (3.10) below.

The linear obstacle problem

We will begin with second-order approximation schemes for . For simplicity, we consider a rectangular domain in , . Then the following second-order FD scheme can be formulated on the grid points: where, for some positive integers and , Let . Then, at each of the interior points , the five-point FD approximation of reads Multiply both sides of (3.2) by to have where and at boundary grid points .

Now, consider the following Jacobi iteration for simplicity. Given an initialization , find iteratively as follows.

Algorithm

where at boundary grid points .

Note that Algorithm produces a solution u of which the function value at a point is a simple average of four neighboring values, satisfying the constraint .

Theorem 1

Let û be the limit of the iterates of Algorithm . Then û satisfies the FD discretization of (1.3). That is, where denotes the set of interior grid points and is the set of boundary grid points.

Proof

It is clear to see from Algorithm that Let at an interior point . Then it follows from (3.4)(b) that which implies that On the other hand, let at . Then, since , we must have which implies that . This completes the proof. □

One can easily prove that the algebraic system obtained from (3.3) 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 [23]. 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 û and therefore satisfy (3.5). In this article, variants of Algorithm for the GS and the SOR would be denoted, respectively, by and , where ω is an over-relaxation parameter for the SOR, . For example, is formulated as

where at boundary grid points .

Note that the right side of (3.9)(a) involves updated values wherever available. When , Algorithm becomes Algorithm ; that is, .

The nonlinear obstacle problem

Applying the same arguments for the linear problem (1.3), the Euler-Lagrange equation for the nonlinear minimization problem (1.1) can be formulated as where Thus the solution to the nonlinear problem (3.10) can be considered as a minimal surface satisfying the constraint given by the obstacle function φ.

Since , the nonlinear obstacle problem (3.10) can equivalently be formulated as where

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 , we first rewrite as where both and are the same as ; 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 th pixel, we first compute second-order FD approximations of at , , , and : Then the directional-derivative terms at the pixel point can be approximated by Now, we discretize the surface element as follows: where the right-hand sides are harmonic averages of FD approximations of in x- and y-coordinate directions, respectively. Then it follows from (3.14), (3.16), and (3.17) that where Note that . As for the linear problem, it is easy to prove that the algebraic system obtained from (3.18) is an M-matrix.

Given FD schemes for as in (3.18), the nonlinear obstacle problem (3.12) can be solved iteratively by the Jacobi iteration.

where at boundary grid points .

The superscript on the coefficients , , indicate that they are obtained using the last iterate . Algorithm produces a solution u of which the function value at a point is a weighted average of four neighboring values, satisfying the constraint . One can prove the following corollary, using the same arguments introduced in the proof of Theorem 1.

Corollary 1

Let û be the limit of the iterates of Algorithm . Then û satisfies the FD discretization of (3.12). That is, where denotes the FD scheme of as defined in (3.18) with .

Variants of Algorithm for the GS and the SOR can be formulated similarly as for the linear obstacle problem; they would be denoted respectively by and . In practice, such symmetric coercive optimization problems, the SOR methods are much more efficient than the Jacobi and Gauss-Seidel methods. We will exploit and 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: , 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 [23], Section 4.3, where is the spectral radius of the iteration matrix of the Jacobi method . The iteration matrix can be explicitly presented as a block tridiagonal matrix where is the m-dimensional identity matrix and

For such a matrix , it is well known that Thus it follows from (3.24) and (3.26) that the optimal SOR parameter corresponding to the mesh size h, , can be expressed as where . Hence, for general mesh size h, the corresponding optimal SOR parameter can be found as follows. It is often the case that the calibration (3.28)(a)-(3.28)(b) can be carried out with a small problem, i.e., with 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 .

Let û be the numerical solution of an obstacle problem. Then it would satisfy the discrete obstacle problem (3.5), particularly at all (interior) grid points . However, when the mesh grid is not matching with the free boundary, the obstacle constraint may not be satisfied at all points . This implies that when the mesh is not fine enough, the numerical solution can be underestimated near the free boundary, as shown in Figure 1 in Section 2.2. Note that the error introduced by non-matching grids is in , while the numerical truncation error is in for second-order FD schemes. That is, the underestimation is in , 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 post-processing 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 û and φ be the numerical solution and the lower obstacle, respectively. Then, for example, the characteristic set of contact points can be determined as follows.

As defined in Section 2.2, 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 .

Subgrid determination of the free boundary

Let be a neighboring point with two of its adjacent points are contact points (), as in Figure 2. 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.

Figure 2

Contact points (red solid circle) and neighboring points (blue open circle). The red dashed curve indicates a possible free boundary.

We first focus on the horizontal line segment connecting and in the east (E) direction. Define Then the corresponding linear interpolation between and over the line segment is formulated as Let Since and are a neighboring point and a contact point, respectively, we have If the free boundary passes between and , then there must exist such that . Let be such that represents the intersection between the line segment and the free boundary. Then it can be approximated as follows.

The maximization problem in (4.6) 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 4.3. The quantity will be used as the Dirichlet value on the free boundary.

For other directions D (=W, S, or N), one can define corresponding difference functions as shown in (4.2)-(4.4) for . Then can be obtained by applying (4.7) with being replaced with . When the adjacent point is not a contact point, you may simply set . Thus each neighboring point produces an array of four values and free boundary values for the directions D where .

Assuming that (4.6) has a unique solution and the obstacle is given as a smooth function, the maximum error for the detection of the free boundary using (4.7) is where h is mesh size. It has been numerically verified that the choice is enough for an accurate detection of the free boundary, for which the upper bound of the error becomes .

Nonuniform FD schemes on the neighboring set

Let be a neighboring point. Then would be available for each ; is also available for . Thus the FD scheme for can be formulated over three points as follows. where Similarly, the FD scheme for can be formulated over three points in the y-direction : where

Thus, the post-processing algorithm of the obstacle SOR (3.9), , can be formulated by replacing the two terms in the right side of (3.2) with the right sides of (4.9) and (4.10), and computing in (3.9.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 2-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 i5-3450S 2.80 GHz processor. The optimal relaxation parameter is calibrated with the lowest resolution to find a constant (3.28) and the constant is used for all other cases. For a comparison purpose, we implemented a state-of-the-art method, PDL1P [14], and its parameters ( and in (2.2)) are found heuristically for cases where the parameters are not suggested in [14]. The iterations are stopped when the maximum difference of consecutive iterates becomes smaller than the tolerance ε: where mostly; Section 5.3 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 with , defined by We solve the linear obstacle problem varying resolutions. The tolerance is set hereafter except for examples in Section 5.3. Table 1 presents the number of iterations and CPU (the elapsed time, measured in second) for the linear problem of the non-smooth obstacle (5.3). One can see from the table that our suggested method requires less iterations and converges about one order faster in the computation time than the PDL1P, a state-of-the-art method. We have also implemented the primal-dual hybrid gradient (PDHD) algorithm in [15, 18, 19] for obstacle problems. The PDL1P turns out to be a simple adaptation of the PDHD and their performances are about the same, particularly when μ is set large. For the resolution , Figure 3 depicts the numerical solutions of the PDL1P and the obstacle SOR and their contour lines. For this example, both the PDL1P and the obstacle SOR resulted in almost identical solutions.

Table 1

The number of iterations and CPU for the linear problem of the non-smooth obstacle ( )

(ε=10^−6)	PDL1P		Obstacle SOR
Resolution	Iter	CPU	Iter	CPU
32 × 32	1,284	0.13	84	0.005
64 × 64	1,744	0.71	148	0.03
128 × 128	2,111	3.58	268	0.22
256 × 256	2,097	14.09	525	1.70

For PDL1P [14], set μ = 108, , and .

Figure 3

Solutions to the linear problem for the obstacle ( ) at resolution . (a) The numerical solution by the PDL1P, (b) its contour plot, (c) the numerical solution by the obstacle SOR, and (d) its contour plot.

As the second example, we consider the radially symmetric obstacle with defined by where  , 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 . See Figure 4, in which we give plots of and the true solution .

Figure 4

The true solution to the linear problem for the obstacle ( ) at resolution . (a) The obstacle and (b) the true solution (5.6)

In Table 2, we compare performances of the PDL1P and the obstacle SOR applied for the linear obstacle problem with (5.4). The PDL1P uses the parameters suggested in [14] (, , ). As one can see from the table, our suggested method takes about one order less CPU time than the PDL1P for the computation of the numerical solution. In Figure 5, we show the numerical solutions and the errors produced by the PDL1P 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.

Table 2

-errors, the number of iterations, and the CPU for linear obstacle problem with ( )

(ε=10^−6)	PDL1P		Obstacle SOR
Resolution	L^∞ -error	Iter (CPU)	L^∞ -error	Iter (CPU)
32 × 32	8.91⋅10−3	715 (0.09)	8.69⋅10−3	62 (0.02)
64 × 64	3.01⋅10−3	1,340 (0.60)	3.05⋅10−3	122 (0.04)
128 × 128	7.66⋅10−4	1,971 (3.51)	7.64⋅10−4	244 (0.22)
256 × 256	1.86⋅10−4	2,072 (14.85)	1.88⋅10−4	489 (1.63)

The PDL1P uses the parameters suggested in [14] (μ = 0.1, , ).

Figure 5

Numerical solutions and errors at the resolution. (a)-(b) by the PDL1P and (c)-(d) by the obstacle SOR.

As a more general obstacle problem, we consider the elastic-plastic torsion problem in [22]. 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 , and the force defined by , and where See Figure 6, where the obstacles and the force are depicted.

Figure 6

Elastic-plastic torsion problem. (a) The obstacles ( and ) and (b) the force v.

In Table 3, we present performances of the PDL1P and the obstacle SOR applied for the elastic-plastic torsion problem (5.7). For the PDL1P, we use the parameters suggested in [14] (, , ). As one can see from the table, our suggested method again resulted in the numerical solution about one order faster than the PDL1P measured in the CPU time. In Figure 7, we illustrate the simulated membranes in the equilibrium satisfying (5.7) and their contact sets at resolution . In Figures 7(b) and 7(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.

Table 3

The number of iterations and the CPU time for the elastic-plastic torsion problem ( )

(ε=10^−6)	PDL1P		Obstacle SOR
Resolution	Iter	CPU	Iter	CPU
32 × 32	887	0.13	47	0.02
64 × 64	1,287	0.68	98	0.04
128 × 128	1,609	3.43	193	0.22
256 × 256	1,866	17.27	368	1.58

For the PDL1P, we use the parameters suggested in [14] (μ = 0.1, , ).

Figure 7

The numerical solutions and the contact sets for the elastic-plastic torsion problem at the resolution. (a)-(b) by the PDL1P and (c)-(d) by the obstacle SOR.

Nonlinear obstacle problems

The obstacle SOR is implemented for nonlinear obstacle problems as described in Section 3.2.

In Table 4, we present experiments for which the obstacle SOR is applied for nonlinear obstacle problems with , . From a comparison with linear cases presented in Tables 1, 2, and 3, 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 (3.19). For , the nonlinear solution is plotted in Figure 8. Compared with the linear solutions in Figure 3, 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 (3.19), 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 9.

Table 4

The performance of the obstacle SOR applied for nonlinear obstacle problems with ,

(ε=10^−6)	φ_1		φ_2		φ_3
Resolution	Iter	CPU	Iter	CPU	Iter	CPU
32 × 32	79	0.02	62	0.04	47	0.04
64 × 64	133	0.16	121	0.17	98	0.16
128 × 128	257	1.25	239	1.16	192	1.09
256 × 256	513	10.19	477	9.20	368	8.24

Figure 8

Nonlinear obstacle problem with at resolution . (a) The nonlinear numerical solution and (b) its contour plot.

Figure 9

The difference between the linear solution and the nonlinear solution, at the resolution, for the obstacle problem with . (a) and (b) its density plot.

Post-processing algorithm

In Figure 5, 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 4 as a post-processing (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 and subgrid FD schemes are applied at neighboring grid points as in Section 4, and (c) another round of iterations is applied to satisfy the tolerance .

First, we consider a step function for an one-dimensional (1D) obstacle, as in Section 2.2. Let and defined by The analytic solution to the linear problem is given as

Figure 10 shows the numerical solutions to the linear problem associated to (5.10) 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: . The error is reduced to after the post-process.

Figure 10

The post-processing algorithm applied for an obstacle problem in 1D at resolution 16. (a) The computed solutions without the post-process (blue solid curve) and with the post-process (red dotted curve with × marks) and (b) their errors. Here the subscript pp indicates the post-process.

The post-processing algorithm is applied to the linear obstacle problem in 2-D involving . Table 5 contains efficiency results that compare performances of the PDL1P, the obstacle SOR (without post-process), and the obstacle SOR with the post-process (Obstacle SOR+PP) at various resolutions; while Table 6 presents an accuracy comparison for those methods. According to Table 5, the post-processed solution requires about 40% more iterations than the non-post-processed one; the incorporation of the post-process makes the iterative algorithm as twice expensive measured in CPU time as the original iteration. However, one can see from Table 6 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 PDL1P requires 10.47 seconds and the obstacle SOR completes the task in 0.84 seconds; when the obstacle SOR+PP takes only 0.1 second.

Table 5

CPU time and iteration comparisons for the suggested post-process, applied to the linear problem with

(ε=10^−7)	PDL1P		Obstacle SOR		Obstacle SOR+PP
Resolution	Iter	CPU	Iter	CPU	Iter	CPU
25 × 25	814	0.07	56	0.02	92	0.07
50 × 50	1,486	0.40	108	0.03	147	0.10
100 × 100	2,092	2.26	213	0.13	299	0.26
200 × 200	2482	10.47	416	0.84	570	1.34

Table 6

and error comparisons for the suggested post-process, applied to the linear problem with

(ε=10^−7)	PDL1P		Obstacle SOR		Obstacle SOR+PP
Resolution	L^∞ error	L^2 error	L^∞ error	L^2 error	L^∞ error	L^2 error
25 × 25	1.94⋅10−2	4.35⋅10−3	1.94⋅10−2	4.38⋅10−3	8.44⋅10−4	2.80⋅10−4
50 × 50	4.38⋅10−3	8.10⋅10−4	4.39⋅10−3	8.41⋅10−4	2.01⋅10−4	6.93⋅10−5
100 × 100	1.25⋅10−3	2.73⋅10−4	1.25⋅10−3	2.87⋅10−4	5.16⋅10−5	1.79⋅10−5
200 × 200	5.45⋅10−4	7.29⋅10−5	5.46⋅10−4	7.76⋅10−5	1.40⋅10−5	4.74⋅10−6

Figure 11 includes plots of the error () at the resolution for the linear obstacle problem with , produced by the PDL1P, the obstacle SOR, and the obstacle SOR+PP. The numerical solutions of the PDL1P 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 11(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.

Figure 11

Plots of the error ( ) at the resolution for the linear obstacle problem with . (a) by the PDL1P, (b) by the obstacle SOR, and (c) by the obstacle SOR+PP.

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 (3.28), we first select . Then by using a trial-by-error method, we found the calibrated optimal relaxation parameter , which results in the following calibrated constant: Thus it follows from (3.27) that the calibrated optimal relaxation parameter reads which is used for the results of the obstacle SOR included in Table 3.

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 3. For and , the line search algorithm returned the curves as shown in Figure 12 with Note that when the calibrated parameters are used, the iteration counts of the obstacle SOR presented in Table 3 are 47, 98, and 193, respectively, for , , and . Thus the calibrated optimal parameters in (5.13) are quite accurate for the optimal convergence.

Figure 12

The relaxation parameter (horizontal axis) vs. the number of iterations (vertical axis) for solving the linear obstacle problem with by the obstacle SOR. (a) when and (b) when .

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.