A Monotone Discretization for the Fractional Obstacle Problem and Its Improved Policy Iteration

Abstract

In recent years, the fractional Laplacian has attracted the attention of many researchers, the corresponding fractional obstacle problems have been widely applied in mathematical finance, particle systems, and elastic theory. Furthermore, the monotonicity of the numerical scheme is beneficial for numerical stability. The purpose of this work is to introduce a monotone discretization method for addressing obstacle problems involving the integral fractional Laplacian with homogeneous Dirichlet boundary conditions over bounded Lipschitz domains. Through successful monotone discretization of the fractional Laplacian, the monotonicity is preserved for the fractional obstacle problem and the uniform boundedness, existence, and uniqueness of the numerical solutions of the fractional obstacle problems are proved. A policy iteration is adopted to solve the discrete nonlinear problems, and the convergence after finite iterations can be proved through the monotonicity of the scheme. Our improved policy iteration, adapted to solution regularity, demonstrates superior performance by modifying discretization across different regions. Numerical examples underscore the efficacy of the proposed method.

1. Introduction

In this work, we consider the obstacle problem associated with the integral form of the fractional Laplacian, referred to as the fractional obstacle problem. For $n\ge 1$, let $\Omega \subset {\mathbb{R}}^n$ be a bounded Lipschitz domain satisfying the exterior ball condition, with the boundary denoted as $\partial \Omega$. Additionally, ${\Omega }^c$ represents the complement of $\Omega$. The specific form of the fractional obstacle problem is as follows: given two prescribed functions $f:\Omega \to \mathbb{R}$ and the obstacle function $\psi :{\mathbb{R}}^n\to \mathbb{R}$ with the nondegeneracy condition ${\psi |}_{{\Omega }^c}<0$ and $s\in (0,1)$, we seek a function $u:{\mathbb{R}}^n\to \mathbb{R}$ that satisfies the following nonlinear equation with homogeneous Dirichlet boundary conditions:

$$\left\{\begin{array}{cc}\hfill \mathcal{G}\left[u\right]:=min\left\{{(-\Delta )}^{s}u-f,u-\psi \right\}& =0,x\in \Omega ,\hfill \\ \hfill u& =0,x\in {\Omega }^c.\hfill \end{array}\right.$$

(1)

Here, the integral fractional Laplacian of order $s\in (0,1)$ is defined by

$${(-\Delta )}^{s}u\left(x\right):=C_{n,s}\mathrm{P}.\mathrm{V}.{\int }_{{\mathbb{R}}^n}\frac{u\left(x\right)-u\left(y\right)}{{|x-y|}^{n+2s}}\mathrm{d}y.$$

(2)

The normalization constant is provided by $C_{n,s}:=\frac{2^{2s}s\Gamma (s+\frac{n}{2})}{{\pi }^{n/2}\Gamma (1-s)}$. According to Animasaun et al., the fractional Laplacian captures nonlocal effects and is used to model a wide range of phenomena [1]. Consequently, the fractional obstacle problem finds applications in various fields, such as mathematical finance theory, systems of particles, and elasticity problems [2,3,4,5,6]. For a comprehensive overview of obstacle problems, including traditional obstacle problems and those related to integral-differential operators (of which the fractional obstacle problem is a particular case), we refer readers to the survey article in [7].

For fractional operators, it has been shown that the maximum principle holds [8]. In terms of numerical stability, it is desirable to preserve this property at the discrete level, known as the monotonicity of the scheme. Constructing monotone schemes is a natural advantage of finite difference methods. For the Laplace operator $(-\Delta )$, the monotone difference scheme has been attended to by researchers for a long time, [9,10]. In the context of the finite element method, achieving monotonicity in the discretization process depends on specific grid conditions [11]. Furthermore, considering that the regularity of the solution is essential for the convergence of numerical schemes, the fractional obstacle problem has received significant attention in mathematics, particularly in the field of PDE regularity of solutions and the free boundary. In this regard, a series of works [12,13,14] have established the Hölder regularity of $C^{1,s}\left({\mathbb{R}}^n\right)$ for solutions to the s-order obstacle problem. For the case of a bounded domain, the interior regularity is expected to exhibit the same $C^{1,s}$ Hölder regularity. However, lower regularity may arise near the boundary of the domain. For the fractional Laplace equation, a well-known result is that when the problem domain satisfies the exterior ball condition and the right-hand side belongs to $L^{\infty }$, the solution exhibits singularity near the boundary $\partial \Omega$ in terms of $\mathrm{dist}{(x,\partial \Omega )}^s$, resulting in global $C^s$ Hölder regularity. This characteristic is observed widely in a large class of nonlocal elliptic equations, as discussed in the survey article [8]. Therefore, the global regularity (boundary regularity) of the fractional obstacle problem on bounded domains requires careful investigation.

In recent years, numerous numerical algorithms have been developed for solving fractional-order partial differential equations (PDEs). For the fractional Laplace equation, the finite element method can be employed, utilizing the variational formulation and accounting for the low regularity near the boundary [15]. Similarly, the fractional obstacle problem can be formulated as a variational inequality problem, where the solution u is sought in order to minimize a functional subject to the constraint $u\ge \psi$ [12]. To address this variational inequality, [16] utilized the finite element method, establishing interior and boundary regularity results for bounded fractional obstacle problems and providing error estimates based on these regularity results. However, for fractional operators, due to their non-local nature, it is challenging to establish general conditions that guarantee the monotonicity of finite element discretization. This challenge is further amplified when dealing with the finite element discretization of nonlinear problems.On the other hand, several finite difference discretizations that satisfy monotonicity have been proposed [17,18,19] for the fractional operator ${(-\Delta )}^s$. However, most of these works focus on problems solved on structured grids, and tend to make overly strong assumptions about the Hölder regularity of the true solution. In our recent work [20], we proposed a monotone difference scheme based on quadrature formulas, which can be applied to solve the fractional Laplace equation on general bounded Lipschitz domains. This work provides a comprehensive analysis of the scheme’s consistency, taking into account the actual regularity of the problem. Furthermore, we obtain rigorous pointwise error estimates by utilizing discrete barrier functions.

The objective of this study is to develop a monotone scheme for solving the problem defined by Equation (1) on a bounded Lipschitz domain and to devise an efficient solution solver. Building upon our previous work [20], we extend this discretization to the fractional obstacle problem. A crucial aspect of our work is the introduction of the enhanced discrete comparison principle, as demonstrated in Lemma 5, for the discrete fractional Laplacian operator. By utilizing this result, we establish that the discretization scheme maintains monotonicity and satisfies the discrete maximum principle, guaranteeing the uniqueness of the discrete solution. Additionally, we employ the discrete Perron method to establish the existence and uniform boundedness of the discrete solution. In addition, this study explores the policy iteration method for solving nonlinear discrete problems. Leveraging the enhanced discrete comparison principle, we establish the finite convergence of the standard policy iteration method. Furthermore, the paper provides an in-depth discussion on the relationship between the discretization scheme for the fractional Laplace equation and the regularity results of the problem. Based on this discussion, we propose an improved policy iteration method that considers the low regularity near the contact set of the problem and selects different discrete scales for different regions. Numerical results demonstrate the superior performance of the improved method compared to the standard policy iteration method.

The remaining parts of this paper are organized as follows. In Section 2, we introduce necessary preliminary results. In Section 3, we present the monotone discretization for the fractional obstacle problem and discuss the properties of the scheme. In Section 4, we apply the policy iteration method for solving nonlinear discrete problems and prove the convergence of the iteration. Additionally, taking into account the specific regularity of the problem, we propose an improved policy iteration method that exhibits better numerical performance in practical computations. Numerical examples are presented in Section 5 to support our theoretical results. Furthermore, for the boundary Hölder regularity of the problem, a simple proof is provided in Appendix A.

2. Preliminary Results

In this section, preliminary knowledge and results are introduced, including notation, definitions of function spaces, and regularity results. Considering an open set $U\subset {\mathbb{R}}^n$ with $\partial U\ne ⌀$, the $C^{\beta }\left(U\right)$ Hölder seminorm (with $\beta >0$) is denoted by ${|·|}_{C^{\beta }\left(U\right)}$. More precisely, for $\beta =k+t$ with k integer and $t\in (0,1]$, we define

$$\begin{array}{cc}\hfill {\left|w\right|}_{C^{\beta }\left(U\right)}={\left|w\right|}_{C^{k,t}\left(U\right)}& :=\underset{x,y\in U,x\ne y}{sup}\frac{|D^{k}w\left(x\right)-D^{k}w\left(y\right)|}{{|x-y|}^t},\hfill \\ \hfill {\parallel w\parallel }_{C^{\beta }\left(U\right)}& :=\sum _{l=0}^k\left(\underset{x\in U}{sup}\left|D^{l}w\left(x\right)\right|\right)+{\left|w\right|}_{C^{\beta }\left(U\right)}.\hfill \end{array}$$

As usual, the nonessential constant denoted by C may vary from line to line; $X\lesssim Y$ means that there exists a constant $C>0$ such that $X\le CY$, and $X\eqsim Y$ means that $X\lesssim Y$ and $Y\lesssim X$.

2.1. The Regularity of the Solution for the Fractional Laplacian Problem

Here, we present a number of well-known conclusions regarding the regularity of the integral fractional Laplacian:

$${(-\Delta )}^s{w}_f=f\mathrm{in}\Omega ,w_f=0\mathrm{in}{\Omega }^c.$$

(3)

The first lemma pertains to the interior regularity of fractional harmonic functions. This lemma, together with the forthcoming Section 2.2, establishes that the solution of Equation (1) attains smoothness within the region not in contact with the obstacle.

Lemma 1

(balayage). Let $w\in L^{\infty }\left({\mathbb{R}}^n\right)$ be such that ${(-\Delta )}^{s}w=0$ in $B_R$. Then, $w\in C^{\infty }\left(B_{R/2}\right)$.

For the global regularity result, however, it can be proven that the solution belongs to $C^s\left({\mathbb{R}}^n\right)$ and this result is sharp [21]. This is due to the limited regularity at the boundary, and suggests the need for a more comprehensive discussion using weighted Hölder spaces, as discussed in [15,21].

2.2. The Regularity of the Solution for the Fractional Obstacle Problem

Next, regularity results for the fractional obstacle problem are presented. It is worth noting that when discussing regularity the assumption that the forcing term f is zero can be made without loss of generality. This is because the general solution can be decomposed as follows:

$$u=u_1+w_f,$$

where $u_1$ solves Equation (1) with zero forcing ($f=0$) and obstacle ${\psi }_1:=\psi -w_f$. Here, $w_f$ is the solution to the linear problem Equation (3).

The contact set is defined as follows:

$$\Lambda =\{x\in \Omega :u(x)=\psi (x\left)\right\}.$$

The first result on the interior Hölder regularity of the fractional obstacle problem was established by Silvester in ${\mathbb{R}}^n$ [12]. When considering the fractional obstacle problem in bounded domains, Nochetto, Borthagaray, and Salgado in [16] derived a corresponding result under the assumption that

$$\mathrm{dist}(\Lambda ,\partial \Omega )=r_0>0.$$

(4)

The main technique employed in their work is the use of truncation functions. The result is presented as follows.

Proposition 1

(interior Hölder regularity). Let Ω be a bounded Lipschitz domain and let $\psi \in C\left(\overline{\Omega }\right)\cap C^{2,1}(\Omega )$ with tthe nondegeneracy condition ${\psi |}_{{\Omega }^c}<0$. Then, the solution u of Equation (1) with $f=0$ satisfying $u\in C^{1,s}\left(D\right)$ for all $D\subset \Omega$ is compact.

Similarly, utilizing truncation function techniques, [16] provided the following boundary Hölder regularity result:

Proposition 2

(boundary Hölder regularity). Let Ω be a bounded Lipschitz domain satisfying the exterior ball condition, and let $\psi \in C\left(\overline{\Omega }\right)\cap C^{2,1}(\Omega )$ with the nondegeneracy condition ${\psi |}_{{\Omega }^c}<0$. Let u solve Equation (1); then,

$${\parallel u\parallel }_{C^s\left({\mathbb{R}}^n\right)}\le C(\Omega ,f,\psi ).$$

In fact, in the boundary Hölder regularity estimate, the smoothness assumption on $\psi$ can be relaxed to $\psi \in L^{\infty }(\Omega )$. The main idea is to revert back to the fundamental theoretical framework of boundary regularity [21] and explore its extension to the fractional obstacle problem. Additionally, compared to [16], the proof is more direct. The detailed proof of this new approach is included in Appendix A.

Before concluding this section, we provide an intuitive summary of the implications of the regularity results for the design of numerical schemes:

• Lemma 1 (balayage) indicates that the solutions are smooth within both $\Lambda$ and $\Omega \setminus \Lambda$ if $\psi$ is smooth. Therefore, the standard discretization methods would be suitable for these regions.
• Combining the interior regularity discussed in Proposition 1 suggests that the smoothness across the contact boundary decreases to $C^{1,s}$. As a result, discretization methods that rely on higher-order derivatives would be inappropriate for this region.
• The boundary regularity result, Proposition 2, is similar to that of linear problems. Thus, the techniques used for linear problems, such as the regularity under weighted norms, remain applicable. From a numerical perspective, this phenomenon motivates the use of graded grids to enhance the convergence order, as explored in the work by [15,20]. This issue is revisited in Section 5.

3. A Monotone Discretization

For simplicity, it is assumed that $\Omega$ is a polygon in 2D and a polyhedron in 3D. Let ${\mathcal{T}}_h$ be a triangulation of the computation domain $\Omega$, i.e., ${\cup }_{T\in {\mathcal{T}}_h}\overline{T}=\overline{\Omega }$. For $T\in {\mathcal{T}}_h$, let $h_T$ denote the diameter of element T and let ${\rho }_T$ denote the radius of the largest inscribed ball contained in T. The triangulation is referred to as local quasi-uniform if there exists a constant ${\lambda }_1$ such that

$$h_T\le {\lambda }_1{h}_{T^{\prime }}\forall \overline{T}\cap \overline{T^{\prime }}\ne ⌀.$$

Meanwhile, the triangulation ${\mathcal{T}}_h$ is called shape regular if there is ${\lambda }_2>0$ such that $h_T\le {\lambda }_2{\rho }_T$. For $n\ge 2$, the shape regularity property implies the local quasi-uniform properties; see [22]. Let ${\mathcal{N}}_h$ denote the node set of ${\mathcal{T}}_h$ and let ${\mathcal{N}}_h^b:=\left\{x_i\in {\mathcal{N}}_h:x_i\in \partial \Omega \right\}$ be the collection of boundary nodes; then, the interior node set is denoted by ${\mathcal{N}}_h^0:={\mathcal{N}}_h\setminus {\mathcal{N}}_h^b$.

Let ${\mathbb{V}}_h:=\left\{u\in C\left({\mathbb{R}}^n\right):u{{|}_T\in {\mathcal{P}}_1\left(T\right),u|}_{{\Omega }^c}=0\right\}$, where ${\mathcal{P}}_1\left(T\right)$ denotes the linear function collection defined on element T. In [20], a monotone discretization for the integral fractional Laplacian was proposed and a corresponding pointwise error estimate was conducted under a realistic Hölder regularity assumption. In this work, the monotone discretization for the integral fractional Laplacian is applied, denoted by ${(-\Delta )}_h^s$:

$${(-\Delta )}_h^s\left[v_h\right]\left(x_i\right):=\underset{:={\mathcal{L}}_h^S\left[v_h\right]\left(x_i\right)}{\underbrace{-{\kappa }_{n,s,i}\frac{{\Delta }_{FD}{v}_h(x_i;H_i)}{H_i^{2s}}}}+\underset{:={\mathcal{L}}_h^T\left[v_h\right]\left(x_i\right)}{\underbrace{{\int }_{{\Omega }_i^c}\frac{v_h\left(x_i\right)-v_h\left(y\right)}{|x_i{-y|}^{n+2s}}\mathrm{d}y}}\forall v_h\in {\mathbb{V}}_h.$$

(5)

Here, ${\Omega }_i\subset \Omega$ is a star-shaped domain centered at $x_i$ satisfying some symmetrical condition ([20] Equations (3.5) and (3.6)) with a typical scaling:

$$H_i:=h_i^{{\alpha }_i}{\delta }_i^{1-{\alpha }_i}\forall x_i\in {\mathcal{N}}_h^0,$$

(6)

where $h_i$ is the mesh size around $x_i$ and ${\delta }_i:=\mathrm{dist}(x_i,\partial \Omega )$. A concrete example that satisfies the condition is the n-cubic domain ${\Omega }_i=x_i+{[-H_i,H_i]}^n$, which is utilized in our numerical experiments. The parameter ${\alpha }_i$ is carefully selected to address the regularity concern of the integral fractional Laplacian [20,21].

The singular integral within ${\Omega }_i$ is approximated using the finite difference method, denoted as ${\Delta }_{FD}u(x_i;H_i):={\sum }_{j=1}^{n}u(x_i+H_i{e}_j)-2u\left(x_i\right)+u(x_i-H_i{e}_j)$, where $e_j$ represents the unit vector of the jth coordinate. The coefficient ${\kappa }_{n,s,i}$ is a known positive constant, and is provided in ([20], Equation 3.10).

3.1. Properties of ${(-\Delta )}_h^s$

In this subsection, several fundamental properties of the ${(-\Delta )}_h^s$ operator are revisited and established. These properties are utilized in subsequent discussions. One of the key features is its monotonicity, as proven in ([20], Lemma 3.4), which is closely related to the matrix discretization being an M-matrix. An even stronger property of ${(-\Delta )}_h^s$ is explored here, namely, that the discretization matrix exhibits strong diagonal dominance. Additionally, the consistency of ${(-\Delta )}_h^s$ is discussed.

To begin, let us revisit the discrete barrier function and monotone property provided in ([20], Lemmas 3.4 and 3.5).

Lemma 2

(discrete barrier function). Let $b_h\in {\mathbb{V}}_h$ satisfy $b_h\left(x_i\right):=1$ for all $x_i\in {\mathcal{N}}_h^0$. It follows that

$${(-\Delta )}_h^s\left[b_h\right]\left(x_i\right)\ge C{\delta }_i^{-2s}>0\forall x_i\in {\mathcal{N}}_h^0,$$

(7)

where the constant C depends only on s and Ω.

Lemma 3

(monotonicity of ${(-\Delta )}_h^s$). Let $v_h,w_h\in {\mathbb{V}}_h$. If $v_h-w_h$ attains a non-negative maximum at an interior node $x_i\in {\mathcal{N}}_h^0$, then ${(-\Delta )}_h^s\left[v_h\right]\left(x_i\right)\ge {(-\Delta )}_h^s\left[w_h\right]\left(x_i\right)$.

The above monotonicity result is insufficient for the obstacle problem, as the fractional order operator equation in the obstacle problem holds only in certain regions of the domain. Let us denote the resulting discrete matrix of ${(-\Delta )}_h^s$ as $\mathbf{L}\in {\mathbb{R}}^{N\times N}$, where N is the number of interior nodes. Below, it is shown that $\mathbf{L}$ is a strongly diagonal dominant M-matrix.

Lemma 4

(strongly diagonal dominant M-matrix). The discrete matrix $\mathbf{L}\in {\mathbb{R}}^{N\times N}$ satisfies

$$\begin{array}{}\hfill \mathrm{(8a)}& {\mathbf{L}}_{ii}>0\forall i;{\mathbf{L}}_{ij}\le 0\forall i,jandi\ne j;\hfill \hfill \mathrm{(8b)}& {\mathbf{L}}_{ii}>\sum _{i\ne j}^N\left|{\mathbf{L}}_{ij}\right|\forall i.\hfill \end{array}$$

Proof. 

The property in Equation (8a) follows directly from the definition of ${(-\Delta )}_h^s$. In fact, for Equation (5) at $x_i$, it is straightforward to observe that ${\mathbf{L}}_{ii}=2n{\kappa }_{n,s,i}{H}_i^{-2s}+{\int }_{{\Omega }_i^c}{|x_i-y|}^{-(n+2s)}\mathrm{d}y$. Furthermore, all the other contributions to the off-diagonal terms are non-positive.

Utilizing Lemma 2, the discrete barrier function can be expressed as the vector $B\in {\mathbb{R}}^N$, where all entries are equal to one. Therefore, Equation (7) implies that

$${\left(\mathbf{L}B\right)}_i={\mathbf{L}}_{ii}+\sum _{i\ne j}^N{\mathbf{L}}_{ij}={\mathbf{L}}_{ii}-\sum _{i\ne j}^N\left|{\mathbf{L}}_{ij}\right|\ge C{\delta }_i^{-2s}>0,$$

which leads to Equation (8b).    □

Next, the enhanced discrete comparison principle applicable to the obstacle problem is presented. This principle plays a crucial role in the subsequent analysis.

Lemma 5

(enhanced discrete comparison principle for ${(-\Delta )}_h^s$). Let $v_h,w_h\in {\mathbb{V}}_h$ be such that

$${(-\Delta )}_h^s\left[v_h\right]\left(x_i\right)\ge {(-\Delta )}_h^s\left[w_h\right]\left(x_i\right)\forall x_i\in {\Lambda }_h,$$

(9a)

$$v_h\left(x_i\right)\ge w_h\left(x_i\right)\forall x_i\in {\mathcal{N}}_h^0\setminus {\Lambda }_h,$$

(9b)

where ${\Lambda }_h\subset {\mathcal{N}}_h^0$ is an arbitrary subset. Then, $v_h\ge w_h$ in Ω.

Proof. 

Because $v_h,w_h$ are piecewise linear functions, it suffices to prove $v_h\left(x_i\right)\ge w_h\left(x_i\right)$ for all $x_i\in {\mathcal{N}}_h^0$. Let $z_h:=v_h-w_h$ and let $Z\in {\mathbb{R}}^N$ denote its coefficient under the basis function. Our objective is to demonstrate that $Z\ge 0$. Note that the resulting discrete matrix of ${(-\Delta )}_h^s$, denoted as $\mathbf{L}$, exhibits strong diagonal dominance. Then, Equation (9) has the matrix representation

$$\mathbf{L}Z=\left(\begin{array}{cc}{\mathbf{L}}_{11}& {\mathbf{L}}_{12}\\ 0& \mathbf{I}\end{array}\right)\left(\begin{array}{c}{Z}^1\\ Z^2\end{array}\right)\ge 0.$$

Here, the indices of all interior points in ${\Lambda }_h$ are combined into the first group. The remaining indices form the second group. This implies that ${\mathbf{L}}_{11}{Z}^1+{\mathbf{L}}_{12}{Z}^2\ge 0$ and $Z^2\ge 0$.

Per the above lemma (strongly diagonal dominant M-matrix), the entries in the off-diagonal block ${\mathbf{L}}_{12}$ are non-positive, which yields

$${\mathbf{L}}_{11}{Z}^1\ge -{\mathbf{L}}_{12}{Z}^2\ge 0.$$

Because $\mathbf{L}$ is a strongly diagonal dominant M-matrix, so is ${\mathbf{L}}_{11}$, that is ${\left({\mathbf{L}}_{11}^{-1}\right)}_{ij}\ge 0$. Therefore, it can be deduced that $Z^1\ge 0$ and $Z\ge 0$, which completes the proof.    □

3.2. Numerical Scheme

Now, the numerical scheme for solving the fractional obstacle problem is proposed. Find $u_h\in {\mathbb{V}}_h$ such that

$$\begin{array}{c}\hfill {\mathcal{G}}_h\left[u_h\right]\left(x_i\right):=min\left\{{(-\Delta )}_h^s\left[u_h\right]\left(x_i\right)-f\left(x_i\right),u_h\left(x_i\right)-\psi \left(x_i\right)\right\}=0\forall x_i\in {\mathcal{N}}_h^0.\end{array}$$

(10)

It should be noted that, although not explicitly stated, there is a parameter $\alpha$ in the discretization of the integral fractional Laplacian that is influenced by the regularity of the solution u. The discrete comparison principle for ${\mathcal{G}}_h$ is demonstrated, from which the uniqueness of Equation (10) follows directly.

Lemma 6

(discrete comparison principle for ${\mathcal{G}}_h$). Let $v_h,w_h\in {\mathbb{V}}_h$ be such that

$${\mathcal{G}}_h\left[v_h\right]\left(x_i\right)\ge {\mathcal{G}}_h\left[w_h\right]\left(x_i\right)\forall x_i\in {\mathcal{N}}_h^0.$$

Then, $v_h\ge w_h$ in Ω.

Proof. 

Because $v_h,w_h\in {\mathbb{V}}_h$, it suffices to prove that $v_h\left(x_i\right)\ge w_h\left(x_i\right)$ for all $x_i\in {\mathcal{N}}_h^0$. The interior points are classified into two sets based on the values of ${\mathcal{G}}_h\left[w_h\right]$:

$${\Lambda }_h^w:=\{x_i\in {\mathcal{N}}_h^0:{(-\Delta )}_h^s\left[w_h\right]\left(x_i\right)-f\left(x_i\right)\ge w_h\left(x_i\right)-\varphi \left(x_i\right)\}.$$

For any $x_i\in {\Lambda }_h^w$, the following holds:

$$v_h\left(x_i\right)-\psi \left(x_i\right)\ge {\mathcal{G}}_h\left[v_h\right]\left(x_i\right)\ge {\mathcal{G}}_h\left[w_h\right]\left(x_i\right)=w_h\left(x_i\right)-\psi \left(x_i\right),$$

Furthermore, for any $x_i\in {\mathcal{N}}_h^0\setminus {\Lambda }_h^w$, the following inequalities are true:

$${(-\Delta )}_h^s\left[v_h\right]\left(x_i\right)-f\left(x_i\right)\ge {\mathcal{G}}_h\left[v_h\right]\left(x_i\right)\ge {\mathcal{G}}_h\left[w_h\right]\left(x_i\right)={(-\Delta )}_h^s\left[w_h\right]\left(x_i\right)-f\left(x_i\right).$$

This leads to the desired result by taking ${\Lambda }_h={\Lambda }_h^w$ in Lemma 5 (enhanced discrete comparison principle for ${(-\Delta )}_h^s$).    □

Corollary 1

(uniqueness). For the problem Equation (10), the discrete solution is unique.

3.3. Stability and Existence

The existence and stability of Equation (10) is shown in this section.

Theorem 1

(existence and stability). There exists a unique $u_h\in {\mathbb{V}}_h$ that solves Equation (10). The solution $u_h$ is stable in the sense that $\parallel u_h{\parallel }_{L^{\infty }(\Omega )}\le C$, where the constant C is independent of h.

Proof. 

To establish existence, a monotone sequence of discrete functions ${\left\{u_h^k\right\}}_{k=0}^{\infty }$ is constructed starting with an initial guess $u_h^0\in {\mathbb{V}}_h$ that satisfies the following condition:

$${\mathcal{G}}_h\left[u_h^0\right]\left(x_i\right)\ge 0\forall x_i\in {\mathcal{N}}_h^0.$$

(11)

Step 1 (existence of $u_h^0$). Let $u_h^0:=E{b}_h\in {\mathbb{V}}_h$, where the discrete barrier function $b_h$ is defined in Equation (7); the constant $E>0$ is specified later. Per Lemma 2 (discrete barrier function), the following relation holds:

$$\begin{array}{cc}\hfill {(-\Delta )}_h^s\left[E{b}_h\right]\left(x_i\right)-f\left(x_i\right)& \ge CE{\delta }_i^{-2s}-f\left(x_i\right)\ge CE\mathrm{diam}{(\Omega )}^{-2s}-{\parallel f\parallel }_{L^{\infty }(\Omega )}\hfill \\ \hfill E{b}_h\left(x_i\right)-\psi \left(x_i\right)& \ge E-{\parallel \psi \parallel }_{L^{\infty }(\Omega )}.\hfill \end{array}$$

By setting $E=max\{\parallel \psi {\parallel }_{L^{\infty }(\Omega )},C^{-1}\mathrm{diam}{(\Omega )}^{2s}\parallel f{\parallel }_{L^{\infty }(\Omega )}\}$, Equation (11) is ensured.

Step 2 (Perron construction). Here, induction is employed. Suppose that there already exists a discrete function $u_h^k\in {\mathbb{V}}_h$ that satisfies

$${\mathcal{G}}_h\left[u_h^k\right]\left(x_i\right)\ge 0\forall x_i\in {\mathcal{N}}_h^0.$$

(12)

The construction of $u_h^{k+1}\in {\mathbb{V}}_h$ with the properties $u_h^{k+1}\le u_h^k$ while satisfying Equation (12) is as follows. All interior nodes are considered sequentially, and auxiliary functions $u_h^{k,i-1}\in {\mathbb{V}}_h$ are constructed using the first $i-1$ nodes, starting with $u_h^{k,0}:=u_h^k$. At $x_i\in {\mathcal{N}}_h^0$, whether or not ${\mathcal{G}}_h\left[u_h^{k,i-1}\right]>0$ is checked. If this is the case, the value of $u_h^{k,i-1}\left(x_i\right)$ is decreased, resulting in the function $u_h^{k,i}$, until

$${\mathcal{G}}_h\left[u_h^{k,i}\right]\left(x_i\right)=0.$$

This is achievable because the discrete matrix of ${(-\Delta )}_h^s$ satisfies ${\mathbf{L}}_{ii}>0$, as stated in Equation (8a). Moreover, at other nodes we have $x_j\ne x_i$, ${\mathbf{L}}_{ij}\le 0$, which implies that

$${\mathcal{G}}_h\left[u_h^{k,i}\right]\left(x_j\right)\ge {\mathcal{G}}_h\left[u_h^{k,i-1}\right]\left(x_j\right)\ge 0\forall x_j\ne x_i.$$

This process is repeated with the remaining nodes $x_j$ for $i<j\le N$, with $u_h^{k+1}:=u_h^{k,N}$ set as the last intermediate function. By construction, we obtain the following:

$${\mathcal{G}}_h\left[u_h^{k+1}\right]\left(x_i\right)\ge 0,u_h^{k+1}\left(x_i\right)\ge u_h^k\left(x_i\right)\forall x_i\in {\mathcal{N}}_h^0.$$

Step 3 (convergence). The sequence ${\left\{u_h^k\left(x_i\right)\right\}}_{k=1}^{\infty }$ is monotonically decreasing and clearly bounded from below by $\psi \left(x_i\right)$. Hence, the sequence converges, and the limit

$$u_h\left(x_i\right):=\underset{k\to \infty }{lim}{u}_h^k\left(x_i\right)\forall x_i\in {\mathcal{N}}_h^0$$

defines $u_h\in {\mathbb{V}}_h$, which satisfies the desired equality

$${\mathcal{G}}_h\left[u_h\right]\left(x_i\right)=0\forall x_i\in {\mathcal{N}}_h^0,$$

because ${\mathcal{G}}_h\left[u_h\right]\left(x_i\right)=\underset{k\to \infty }{lim}{\mathcal{G}}_h\left[u_h^k\right]\left(x_i\right)\ge 0$. If the last inequality were strict, Step 2 could be applied to further improve $u_h$. This demonstrates the existence of a discrete solution to Equation (10). Furthermore, the above proof implies

$${\parallel u_h\parallel }_{L^{\infty }(\Omega )}\le max\{{\parallel \psi \parallel }_{L^{\infty }(\Omega )},C^{-1}\mathrm{diam}{(\Omega )}^{2s}{\parallel f\parallel }_{L^{\infty }(\Omega )}\},$$

which is the uniform bound, as was asserted.    □

Remark 1

(discussion on consistency and convergence). An important concept in numerical methods is consistency, which involves investigating the property of $|\mathcal{G}\left[v\right]\left(x\right)-{\mathcal{G}}_h\left[{\mathcal{I}}_{h}v\right]\left(x_i^h\right)|$ tending to zero as the mesh size decreases. In fact, for the integral fractional Laplacian operator, the detailed analysis provided in [20] reveals that the discretization achieves consistency at interior points that are a constant number of mesh sizes away from the boundary. The inconsistency near boundary points can be addressed by incorporating the discrete barrier function in Equation (7). This analysis technique is commonly used in the convergence analysis of semi-Lagrangian or two-scale methods [23,24], and can be viewed as an extension of the Barles–Souganidis [25] analytical framework for nonlinear problems. Due to space limitations, elaboration of the corresponding results is not provided here.

4. An Efficient Solver for the Discrete Nonlinear Problem

In this section, the policy iteration (Howard’s algorithm) is employed to solve the discrete fractional obstacle problem Equation (10). It is worth noting that policy iteration is a well-established and extensively studied technique in dynamic control problems. For the problem Equation (10), the convergence of the policy iteration is established by leveraging the monotonicity of the discrete operator. Furthermore, an improved policy iteration is proposed by incorporating prior knowledge, such as the regularity of the solution to the fractional obstacle problem discussed in Section 2.

4.1. Policy Iteration

Policy iteration is utilized to solve the discrete problem Equation (10). The algorithm is outlined as Algorithm 1.    

Algorithm 1: Policy iteration for the fractional obstacle problem Equation (10)

Input: right hand side f, obstacle function $\psi$. Output: solution $u_h$. Initialize $u_h^{\left(0\right)}\left(x_i\right)=\psi \left(x_i\right)\forall x_i\in {\mathcal{N}}_h^0$; Iterate for $k\ge 0$: 1.  Update discrete contact set ${\mathcal{C}}_h^{(k+1)}$:   ${\mathcal{C}}_h^{(k+1)}:=\{x_i:{(-\Delta )}_h^s\left[u_h^{\left(k\right)}\right]\left(x_i\right)-f\left(x_i\right)\ge u_h^{\left(k\right)}\left(x_i\right)-\psi \left(x_i\right)\}$. 2.  Update solution $u_h^{(k+1)}$ such that $$\begin{array}{cc}\hfill u_h^{(k+1)}\left(x_i\right)-\psi \left(x_i\right)& =0\forall x_i\in {\mathcal{C}}_h^{(k+1)},\hfill \\ \hfill {(-\Delta )}_h^s\left[u_h^{(k+1)}\right]\left(x_i\right)-f\left(x_i\right)& =0\forall x_i\in {\mathcal{N}}_h^0\setminus {\mathcal{C}}_h^{(k+1)}.\hfill \end{array}$$ (13) 3.  If $u_h^{(k+1)}=u_h^{\left(k\right)}$, then stop; Otherwise, continue to Step 1.

There exists a convergence result for policy iteration when solving the nonlinear problem Equation (10), which was initially proven in [26]. Here, a brief overview of the proof using our notation is provided for clarity. It is worth noting that the monotonicity property plays a crucial role in the convergence analysis.

Theorem 2

(convergence of policy iteration). The sequence ${\left\{u_h^{\left(k\right)}\right\}}_{k=0}^{\infty }$ provided by the above algorithm satisfies the following:

1.

$u_h^{\left(k\right)}\left(x_i\right)\le u_h^{(k+1)}\left(x_i\right)$ for all $k\ge 0$ and $x_i\in {\mathcal{N}}_h^0$.

2.

The algorithm converges in at most N iterations.

Proof. 

With $u_h^{\left(0\right)}\left(x_i\right)=\psi \left(x_i\right)$ for all $x_i\in {\mathcal{N}}_h^0$, ${\mathcal{C}}_h^{\left(0\right)}$ can be defined as ${\mathcal{N}}_h^0$ to ensure that Equation (13) holds for the initialization.

Step 1 (increasing sequence). Equation (13) implies that, for $k\ge 0$,

$${\mathcal{G}}_h\left[u_h^{\left(k\right)}\right]\left(x_i\right)=min\{{(-\Delta )}_h^s\left[u_h^{\left(k\right)}\right]\left(x_i\right)-f\left(x_i\right),u_h^{\left(k\right)}\left(x_i\right)-\psi \left(x_i\right)\}\le 0\forall x_i\in {\mathcal{N}}_h^0.$$

Then, the definition of ${\mathcal{C}}_h^{(k+1)}$ in Step 1 implies

$$\begin{array}{cc}\hfill u_h^{\left(k\right)}\left(x_i\right)-\psi \left(x_i\right)& \le 0\forall x_i\in {\mathcal{C}}_h^{(k+1)},\hfill \\ \hfill {(-\Delta )}_h^s\left[u_h^{\left(k\right)}\right]\left(x_i\right)-f\left(x_i\right)& \le 0\forall x_i\in {\mathcal{N}}_h^0\setminus {\mathcal{C}}_h^{(k+1)}.\hfill \end{array}$$

Next, using the update condition in Equation (13), the following relation holds:

$$\begin{array}{cc}\hfill u_h^{\left(k\right)}\left(x_i\right)-\psi \left(x_i\right)& \le u_h^{(k+1)}\left(x_i\right)-\psi \left(x_i\right)\forall x_i\in {\mathcal{C}}_h^{(k+1)},\hfill \\ \hfill {(-\Delta )}_h^s\left[u_h^{\left(k\right)}\right]\left(x_i\right)-f\left(x_i\right)& \le {(-\Delta )}_h^s\left[u_h^{\left(k\right)}\right]\left(x_i\right)-f\left(x_i\right)\forall x_i\in {\mathcal{N}}_h^0\setminus {\mathcal{C}}_h^{(k+1)}.\hfill \end{array}$$

By combining Lemma 5 (enhanced discrete comparison principle for ${(-\Delta )}_h^s$), the increasing sequence is shown as follows:

$$u_h^{\left(k\right)}\left(x_i\right)\le u_h^{(k+1)}\left(x_i\right)\forall x_i\in {\mathcal{N}}_h^0.$$

Step 2 (strictly decreasing discrete contact set). First, the increasing sequence ${\left\{u_h^{\left(k\right)}\right\}}_{k=0}^{\infty }$ implies

$$u_h^{\left(k\right)}\left(x_i\right)-\psi \left(x_i\right)\ge 0,\forall x_i\in {\mathcal{N}}_h,\forall k\ge 0.$$

Now, for any $x_j\in {\mathcal{N}}_h^0\setminus {\mathcal{C}}_h^{\left(k\right)}$, the update condition Equation (13) and the above property imply

$${(-\Delta )}_h^s\left[u_h^{\left(k\right)}\right]\left(x_j\right)-f\left(x_j\right)=0\le u_h^{\left(k\right)}\left(x_j\right)-\psi \left(x_j\right),$$

whence we have $x_j\in {\mathcal{N}}_h^0\setminus {\mathcal{C}}_h^{(k+1)}$ for the definition of ${\mathcal{C}}_h^{(k+1)}$ in Step 1. This means that

$${\mathcal{N}}_h^0\setminus {\mathcal{C}}_h^{\left(k\right)}\subset {\mathcal{N}}_h^0\setminus C_h^{(k+1)},\mathrm{or}{\mathcal{C}}_h^{(k+1)}\subset {\mathcal{C}}_h^{\left(k\right)}\forall k\ge 0.$$

Further, if ${\mathcal{C}}_h^{(k+1)}={\mathcal{C}}_h^{\left(k\right)}$, then the iteration clearly stops at the $(k+1)$-th step due to the update condition in Equation (13); that is, the discrete contact set is strictly decreasing unless the iteration stops, which makes for at most N iterations.    □

Remark 2

(convergence profile of Algorithm 1). The N-step iteration is proven to be sharp in the general case, as demonstrated by a specific example presented in [27]. However, in numerical experiments convergence is often observed within only a few iterations.

4.2. An Improved Policy Iteration

Next, improvements are introduced to the discrete operator in the policy iteration algorithm. It is observed that the scaling $H_i$ of the singular integral part of ${(-\Delta )}_h^s$ over the region ${\Omega }_i$ can be adjusted, as shown by the expression in Equation (6). In the policy iteration, modifications are made to the selection of $H_i$, primarily for the following reasons:

• During the updating of $u_h$ in Step 3, if the discrete contact set ${\mathcal{C}}_h$ has been updated, Equation (13) can be seen as the discrete fractional Laplacian on the non-contact points ${\mathcal{N}}_h^0\setminus {\mathcal{C}}_h$. Instead, the values at the contact points ${\mathcal{C}}_h$ can be viewed as the Dirichlet “boundary” data.
• In the discretization of the fractional Laplacian, the singular part is approximated by a scaled Laplace operator [20]. However, it is important to note that the regularity estimate for the obstacle problem, Proposition 1, indicates that the true solution has only $C^{1,s}$ continuity over the boundary of the contact set. This limitation restricts the accuracy of such approximation.

Based on the above, when discretizing the fractional Laplacian on ${\mathcal{N}}_h^0\setminus {\mathcal{C}}_h$, it is necessary to restrict the size of $H_i$ to ensure that ${\Omega }_i$ does not intersect with the region relating to the discrete contact set ${\mathcal{C}}_h$. Therefore, the following improved strategy is proposed:

$$H_i=min\{h_i^{{\alpha }_i}{\delta }_i^{1-{\alpha }_i},\theta \mathrm{dist}(x_i,{\mathcal{C}}_h)\}\forall x_i\in {\mathcal{N}}_h^0,$$

(14)

where $\theta$ is a constant that depends on the shape regularity of the mesh. In the numerical experiments, $\theta$ is set to $1/4$.

Because the discrete contact set in policy iteration is continuously updated, the improvement in the algorithm primarily lies in updating the discretization of the fractional Laplacian through Equation (14) after updating the contact set. The improved algorithm is outlined as follows.

Remark 3

(convergence profile of Algorithm 2). In the improved algorithm, it is important to note that the solutions in the consecutive steps correspond to different discretized fractional Laplacian operators, as seen in Equation (15). Therefore, it is not possible to directly apply Theorem 2 to provide a rigorous convergence analysis. However, it is worth noting that although the discretized fractional Laplacian operators vary across steps, they are all appropriate discretizations of the continuous operator and maintain monotonicity. From this perspective, the convergence behavior of the improved policy iteration algorithm should be similar to the original version. Indeed, this observation has been confirmed in numerical experiments; the improved algorithm exhibits better performance, particularly near the contact interface.

Algorithm 2: Improved policy iteration for the fractional obstacle problem Equation (10)

Input: right hand side f, obstacle function $\psi$. Output: solution $u_h$. Initialize $u_h^{\left(0\right)}\left(x_i\right)=\psi \left(x_i\right)\forall x_i\in {\mathcal{N}}_h^0$; Initialize discrete fractional Laplacian: ${(-\Delta )}_h^{s,\left(0\right)}={(-\Delta )}_h^s$. Iterate for $k\ge 0$: 1.  Update discrete contact set ${\mathcal{C}}_h^{(k+1)}$:   ${\mathcal{C}}_h^{(k+1)}:=\{x_i:{(-\Delta )}_h^{s,\left(0\right)}\left[u_h^{\left(k\right)}\right]\left(x_i\right)-f\left(x_i\right)\ge u_h^{\left(k\right)}\left(x_i\right)-\psi \left(x_i\right)\}$. 2.  Update discrete fractional Laplacian:   ${(-\Delta )}_h^{s,(k+1)}$ by choosing $H_i=min\{h_i^{{\alpha }_i}{\delta }_i^{1-{\alpha }_i},\theta \mathrm{dist}(x_i,{\mathcal{C}}_h^{(k+1)})\}$. 3.  Update solution $u_h^{(k+1)}$ such that $$\begin{array}{cc}\hfill u_h^{(k+1)}\left(x_i\right)-\psi \left(x_i\right)& =0\forall x_i\in {\mathcal{C}}_h^{(k+1)},\hfill \\ \hfill {(-\Delta )}_h^{s,(k+1)}\left[u_h^{(k+1)}\right]\left(x_i\right)-f\left(x_i\right)& =0\forall x_i\in {\mathcal{N}}_h^0\setminus {\mathcal{C}}_h^{(k+1)}.\hfill \end{array}$$ (15) 4.  If $u_h^{(k+1)}=u_h^{\left(k\right)}$, then stop; Otherwise, continue to Step 1.

5. Numerical Experiments

In this section, several numerical experiments are conducted to test the stability and effectiveness of the methods along with the convergence of the improved policy iteration.

According to the discussion in Section 2.2, the fractional obstacle problem exhibits low regularity at the domain boundary, which is consistent with the fractional linear problem. Various approaches have been proposed to address this issue, including the use of graded grids. In this context, as discussed in [15,16,20], the concept of graded grids is introduced with a parameter h. Let $\mu \ge 1$ be a constant such that for any $T\in {\mathcal{T}}_h$,

$$h_T\eqsim \left\{\begin{array}{cc}{h}^{\mu }\hfill & \hfill \mathrm{if}T\cap \partial \Omega \ne ⌀,\\ h\mathrm{dist}{(T,\partial \Omega )}^{\frac{\mu -1}{\mu }}\hfill & \hfill \mathrm{if}T\cap \partial \Omega =⌀.\end{array}\right.$$

(16)

It is worth noting that the quasi-uniform grid corresponds to the case where $\mu =1$. Additionally, the choice of $\mu$ has an impact on both the grid quality and the overall computational complexity. The value of $\mu$ is specified in each experiment.

In the discretization of the fractional Laplacian operator, the scale $H_i$ of the singular approximation part depends on both the grid and on ${\alpha }_i$, as shown in Equation (6). The optimal choice of $\alpha$ depends on the local smoothness of the solution or the smoothness of the forcing term f, as discussed in [20]. In the tests conducted in this section, a smooth f is assumed; based on the findings in ([20], Theorem 6.1), the optimal choice of $\alpha$ is used, i.e., $\alpha =\frac{1}{2}$.

5.1. Convergence Order Test

In this 1D test, the convergence rate is tested on both uniform and graded grids. The domain considered is $\Omega =(-1,1)$, and an explicit solution for problem Equation (1) is constructed as follows:

$$u\left(x\right)=\frac{2^{-2s}\Gamma (n/2)}{\Gamma (n/2+s)\Gamma (1+s)}{(1-|x|}^2{)}_{+}^s.$$

A direct calculation shows that

$${(-\Delta )}^{s}u=1x\in \Omega ,u=0x\in {\Omega }^c.$$

Next, the following forcing and obstacle terms are considered:

$$f\left(x\right)=1-5(\frac{1}{2}-{\left|x\right|)}_{+}\mathrm{and}\psi =u\left(x\right)-\frac{1}{2}{(x^2-\frac{1}{4})}_{+},$$

such that ${(-\Delta )}^{s}u-f>0$ in $B_{\frac{1}{2}}\left(0\right)$ and ${(-\Delta )}^{s}u-f=0$ outside of this set. Consequently, u represents the solution of problem Equation (1), and the contact set is $B_{\frac{1}{2}}\left(0\right)$.

Figure 1 illustrates the convergence rates for both uniform and graded grids. According to ([20], Section 7.1), the optimal choice for the pointwise estimate of the linear fractional Laplacian is $\mu =\frac{2-s}{s}$, which is used in this test as well. The convergence rate over uniform grids is observed to be approximately s, which is consistent with the behavior observed in the quasi-uniform grids for the linear case. Moreover, the graded grid yields a significantly improved convergence rate. It is worth noting that the convergence order under the optimal choice is $2-s$ for the linear case, which is observed for certain nonlinear cases as well (e.g., $s=0.6$). However, for other cases, the strategy of choosing $H_i$ in Equation (14) based on the local distance to the contact set may lead to different convergence rates compared to those for the linear case.

5.2. Quantitative Behavior and Comparison of Different Policy Iterations

In this section, the problem Equation (1) is investigated over the domain $\Omega =(-1,1)$. The force is set as $f=0$, and the obstacle function is provided by

$$\psi =1-4|x-1/4|.$$

In Figure 2, the numerical solutions $u_h$ obtained for different values of s are presented below. A clear qualitative difference is observed between the solutions for different choices of s. When $s=0.9$, the discrete solution resembles the expected solution of the classical obstacle problem, where the operator ${(-\Delta )}^s$ is effectively replaced by $-\Delta$. As s decreases, the solution approaches ${\psi }_{+}$. Throughout the experiments, a consistent observation is made that as s increases, the contact set decreases, and always contains the point $x_0$.

The difference between the improved method and the original method is remarkable. Observations reveal that when using the original method with larger values of s, oscillations occur near the free boundary due to the handling of the singular part across it. On the other hand, the improved method yields a numerically smoother solution without oscillations. As s decreases, the solutions obtained by the two methods become closer to each other, and the differences diminish.

5.3. Convergence History of Improved Policy Iteration

Next, the convergence history of the improved policy iteration (Algorithm 2) is examined. The experiment is conducted with $f=1$, and the obstacle function is defined as follows:

$$\psi \left(x\right)=3-6|x-1/4|.$$

As mentioned in Remark 3, it is highly likely that the improved algorithm possesses the property of an increasing solution sequence (or decreasing discrete contact set), as stated in Theorem 2 for the standard algorithm. In this experiment, as demonstrated in Figure 3, the convergence profile is observed through the values of $|{\mathcal{C}}_h^{\left(k\right)}|$.

Furthermore, the actual number of iterations obtained in the experiments is presented in

Table 1. Number of iterations for $s=0.3,0.6,0.9$ using graded grids ($\mu =\frac{2-s}{s}$) with a different number of degrees of freedom (DOFs).

	$\mathit{s}=0.3$				$\mathit{s}=0.6$				$\mathit{s}=0.9$
N	126	254	510	1022	118	237	476	953	77	155	311	623
iterations	4	5	5	6	5	6	7	8	6	7	9	11

. It can be observed that the actual number of iterations required is significantly smaller than the theoretical upper bound N. Additionally, as the grid is refined, the algorithm’s performance remains stable. This demonstrates the efficiency of our algorithm.

5.4. Two-Dimensional Test

Finally, the problem Equation (1) is considered in the domain $B_1(0)\subset {\mathbb{R}}^2$ with $f=0$ and the obstacle function provided by

$$\psi =\frac{1}{2}-|x-x_0|,where{x}_0=(\frac{1}{4},\frac{1}{4}).$$

For $s=0.1$, the numerical solution requires 13,395 degrees of freedom and the number of improved policy iterations is 4. Conversely, for $s=0.9$ the numerical solution involves 4253 degrees of freedom and the number of improved policy iterations is 10. The number of iterations increases as s grows, though significantly less than N, as observed in the 1D case. As expected, the 2D numerical results (Figure 4) demonstrate similar qualitative characteristics as observed in the 1D case.

6. Conclusions

In this paper, a discrete scheme for the fractional obstacle problem is proposed based on the monotone discretization for the fractional Laplace operator. Utilizing the distinctive structure of the problem Equation (10) and conducting a comprehensive study of the discrete operator, this study reveals that the nonlinear discrete operator ${\mathcal{G}}_h$ upholds the discrete comparison principle. Based on this property, the existence, uniqueness, and uniform boundedness of numerical solutions are established.

Due to the limited regularity of the true solution near the domain boundary, a graded grid is introduced to capture this behavior. The policy iteration method is used to solve nonlinear problems. Benefiting from the monotonicity of the discrete scheme of the fractional Laplace operator, the policy iterative process can converge in finite iterations. Moreover, considering the reduced regularity of the actual solution near the contact set boundary, the discretization of the fractional Laplacian is adaptively refined by iteratively updating contact nodes. This refinement leads to an improved policy iteration approach. In contrast to the conventional policy iteration, this improved method demonstrates superior numerical performance across a range of numerical experiments. One-dimensional and two-dimensional examples are provided to support the theoretical results.