EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6070 ISSN 1307-5543 – ejpam.com Published by New York Business Global Overlapping Domain Decomposition Methods for Noncoercive Hamilton-Jacobi-Bellman Equation Samar Chebbah1,2,∗, Sofiane Madi3, Mohamed Haiour4 1 Department of Mathematics, Institute of Mathematics and Computer Science, Abdelhafid Boussouf University Center–Mila, Mila, Algeria 2 Laboratoire des Mathématiques Appliquées et Didactique (MAD), École Normale Supérieure El Katiba Assia Djebar (ENSC) – Constantine, Constantine, Algeria 3 Department of Mathematics, Faculty of Mathematics and Computer Science, University of Batna 2, Batna, Algeria 4 Department of Mathematics, Faculty of the Sciences, University Badji Mokhtar–Annaba, P.O. 23000 Annaba, Algeria Abstract. This study presents an approach to demonstrate the monotonic and geometric con- vergence of the non-coercive Hamilton-Jacobi-Bellman problem with Dirichlet boundary condi- tions using the finite difference method. The approach combines overlapping domain decom- position with a sequence of lower and upper solutions to characterize the solution in a discrete setting. Numerical experiments are conducted to validate the methodology and illustrate the consistency of the theoretical findings. 2020 Mathematics Subject Classifications: 65N06, 65N12, 49L25, 65N55 Key Words and Phrases: Hamilton-Jacobi-Bellman equation, noncoercive problems, overlap- ping domain decomposition, convergence analysis, lower and upper solutions, finite difference method 1. Introduction We consider the following non-coercive Hamilton-Jacobi-Bellman (HJB) equation, with Dirichlet boundary conditions :max 1≤i≤n (Aiξ −F i) = 0, inΩ, ξ = 0, on ∂Ω, (1) ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6070 Email addresses: s.chebbah@centre-univ-mila.dz (S. Chebbah), s.madi@univ-batna2.dz (S.Madi), mohamed.haiour@univ-annaba.dz (M. Haiour) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 2 of 20 where Ω ⊂ Rd (d ≥ 1) is a bounded domain with a smooth boundary ∂Ω. The operators Ai are uniformly elliptic and take the form: Ai = d∑ t,s=1 aits(x) ∂2 ∂xt∂xs + d∑ s=1 bis(x) ∂ ∂xs + ai0(x), where the coefficients aits(x), b i s(x), a i 0(x) ∈ C2(Ω) satisfy the following conditions: d∑ t,s=1 aits(x)νtνs ≥ δ|ν|2, ∀ν ∈ Rd, δ > 0, ai0(x) ≥ γ > 0. The noncoercive bilinear form for u, v ∈ H1(Ω) is expressed as follows: ai⟨u, v⟩ = ∫ Ω  d∑ t,s=1 aits(x) ∂u ∂xt ∂v ∂xs + d∑ s=1 bis(x) ∂u ∂xs v+ ai0(x)uv  dx. Additionally, the function F i is assumed to be smooth and nonnegative: F i ∈ C2(Ω), F i ≥ 0. This study examines the monotonic and geometric convergence of problem (1) within a discrete framework using the finite difference method. We employ the overlapping domain decomposition approach along with the characterization of sequences of lower and upper solutions. This approach allows us to analyze the convergence properties of the Hamilton-Jacobi-Bellman (HJB) equation, particularly in its non-coercive form. The absence of a coercive condition in the HJB equation presents a significant challenge, making the application of conventional solution methods more complex. The numerical approximation of Hamilton-Jacobi-Bellman (HJB) equations has been an active research area, with various methods proposed to enhance accuracy and compu- tational efficiency. For instance, the study in [1] explores numerical schemes for systems of HJB equations in innovation dynamics, providing insights into their stability and convergence. Mixed finite element techniques have also been developed to handle HJB equations with Cordes coefficients, as discussed in [2]. To address challenges associated with quasi-variational inequalities, domain decomposition strategies such as the Schwarz method have been analyzed in [3]. Furthermore, overlapping domain decomposition ap- proaches have been applied to non-coercive quasi-variational systems related to HJB equations [4], while [5] introduces a contraction-based algorithmic approach to solving these equations. In this work, we adopt an overlapping domain decomposition method that ensures stability through lower and upper solutions, providing a robust numerical framework for solving HJB equations efficiently. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 3 of 20 In this work,we adopt an overlapping domain decomposition approach, leveraging the interaction between lower and upper solutions to ensure numerical stability and convergence. This makes it a reliable method for solving complex problems. This paper is organized as follows: Section 2 provides a detailed analysis of the transformation of the non-coercive problem (1) into the coercive Hamilton-Jacobi-Bellman (HJB) equation, considering both the continuous and discrete cases using the finite difference method. Section 3 establishes the monotonic and geometric convergence of the solutions by applying the overlapping domain decomposition method along with a sequence of lower and upper solutions. The final section presents numerical experiments that confirm the theoretical findings, illustrating the effectiveness of domain decomposition methods in solving HJB equations and evaluating the impact of overlapping on the convergence rate. 2. Coercive Reformulation of the HJB Equation In this section, we reformulate the Hamilton-Jacobi-Bellman (HJB) equation into a coercive form to improve its stability and numerical solvability. Following [6], this is achieved by introducing a regularization parameter µ > 0, which ensures coercivity and enhances the well-posedness of the problem. The modified equation is given by:max 1≤i≤n (Biξ − Gi(ξ)) = 0, inΩ, ξ = 0, on ∂Ω, (2) where the coercive operators B1, · · · ,Bn are defined as: Bi = d∑ t,s=1 aits(x) ∂2 ∂xt∂xs + d∑ s=1 bis(x) ∂ ∂xs + (ai0(x) + µ), and the corresponding equations are defined as: Bi = Ai + µI and Gi(ξ) = F i + µξ, i = 1, n. The bilinear form for the operator Bi is expressed as: bi⟨u, v⟩ = ai⟨u, v⟩+ µ⟨u, v⟩, Previous studies have established the existence, uniqueness, and regularity of the so- lution to the continuous problem (2) (See [7, 8]). To solve this problem numerically, both finite element and finite difference methods have been employed, transforming the continuous problem into a system of algebraic equations through discretization. This approach simplifies the problem and enhances its computational feasibility. Thus, the discrete version of the problem is formulated as: S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 4 of 20 max 1≤i≤n (Biζ − G i(ζ)) = 0, inΩ, ζ = 0, on ∂Ω, (3) here, Bi and G i(ζ) are the discrete approximations of the continuous operators, facili- tating the numerical solution of the HJB equation. These are given by: Bi = A i + µI, G i(ζ) = F i + µζ. The matrices Bi and A i satisfy the discrete maximum principle and are M-matrices[7], ensuring stability and convergence in the numerical approximation of the HJB equation. 3. Decomposition of the Domain into Overlapping Subdomains In this section, we introduce an approach that combines domain decomposition with overlapping subdomains and the sequence of lower and upper solutions to efficiently solve the Hamilton-Jacobi-Bellman (HJB) equation. The computational domain is di- vided into overlapping subdomains, which enhances numerical stability. We establish the geometric and monotonic convergence towards the discrete solution of the problem. 3.1. Definition of Subdomain Decomposition We examine a regular and bounded domain Ω in R2, which is partitioned into two subdomains, Ω1 and Ω2, that overlap along the x-axis. These subdomains are defined by the following relations: Ω = Ω1 ∪ Ω2 and Ω1 ∩ Ω2 ̸= ∅. The boundaries of each subdomain are denoted by Γ1 = ∂Ω1 and Γ2 = ∂Ω2, and the interface where the subdomains meet is represented by: Σ1 = ∂Ω1 ∩ Ω2, Σ2 = ∂Ω2 ∩ Ω1. Figure 1: Illustration of the domain Ω and its subdomains. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 5 of 20 The specific case considered here is where the domain Ω is defined as [0, 1] × [0, 1]. This domain is then subdivided into two overlapping subdomains: Ω1 = [0, β]× [0, 1], Ω2 = [α, 1]× [0, 1], where α and β represent the boundaries that define the partition. The region of overlap between Ω1 and Ω2 is given by the interval: τ = β − α, with τ representing the length of the overlap. Additionally, the geometric length of the overlap region is quantified as: χ(h) = τ × h, where h is a scaling factor. 3.2. Block Decomposition of the Domain Matrix We consider the matrices Bi of dimension (d× d), which are associated with vectors that can be partitioned into block form as follows: ζ = ( ζ1 ζ2 ) , G i(ζ) = ( G i(ζ1) G i(ζ2) ) , Bi = ( Bi 1,1 −Bi 1,2 −Bi 2,1 Bi 2,2 ) . Here, the vectors ζ1 and ζ2 are defined as: ζ1 = ((ζ1)1, . . . , (ζ1)β−1) T , ζ2 = ((ζ2)α+1, . . . , (ζ2)d) T . To construct the matrix Bi using the domain decomposition method, we partition it into two different forms within each subdomain. The decomposition for the first subdomain is defined as follows: Bi = ( Bi 1,1 −Bi 1,2 C i 1 D i 1 ) , where: • Bi 1,1 and D i 1 are square matrices of dimensions (β− 1)× (β− 1) and (d− β+1)× (d− β + 1) respectively. • Bi 1,2 and C i 1 are of dimensions (β − 1) × (d − β + 1) and (d − β + 1) × (β − 1) respectively. The second partition of the matrix Bi is as follows: Bi = ( D i 2 C i 2 −Bi 2,1 Bi 2,2 ) , where: S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 6 of 20 • Bi 2,2 and D i 2 are square matrices of dimensions (d − α) × (d − α) and α × α respectively. • Bi 2,1 and C i 2 are of dimensions (d− α)× α and α× (d− α) respectively. When β = α + 1 (i.e., τ = 1), the matrices Bi 1,1 and D i 2 coincide, and similarly, Bi 2,2 and D i 1 coincide. This results in minimal geometric overlap and complete exclusion of algebraic overlap. Additionally: Bi 1,1 and Bi 2,2 correspond to the discretized operators under homogeneous Dirichlet boundary conditions. To ensure that the enhanced system remains consistent with the original system despite the overlapping region, the non-diagonal blocks Bi 2,1 and Bi 1,2 are extended by adding rows and columns filled with zeros. The resulting extended matrices are defined as follows: B̃i 1,2 = [ 0β−1,τ−1 Bi 1,2 ] , B̃i 2,1 = [ Bi 2,1 0d−α,τ−1 ] . Where: • B̃i 1,2 is a matrix of size (β−1)× (d−α) that connects a vector from Ω2 to Ω1 with zero entries outside of Ω2. • B̃i 2,1 is a matrix of size (d − α) × (β − 1) that connects a vector from Ω1 to Ω2, with zero entries outside of Ω1. The extended system satisfies the following equations:max 1≤i≤n (B̃iζ − G i(ζ)) = 0, inΩ, ζ = 0, on ∂Ω. In the minimal overlap case (τ = 1), the extended system simplifies to the original matrix associated with the domain partition without any overlap. This is expressed by the following relation: B̃i = ( Bi 1,1 −B̃i 1,2 −B̃i 2,1 Bi 2,2 ) = ( Bi 1,1 −Bi 1,2 −Bi 2,1 Bi 2,2 ) = Bi. 3.3. Iterative Algorithm for Domain Decomposition The iterative process begins with selecting an initial solution, either a lower or an upper one, depending on the nature of the problem. We first define the notions of subsolution and supersolution: • A subsolution ζ̌ ∈ Rm satisfies the inequality: max 1≤i≤n { Biζ̌ − G i(ζ̌) } ≤ 0. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 7 of 20 • A supersolution ζ̂ ∈ Rm satisfies the inequality: max 1≤i≤n { Biζ̂ − G i(ζ̂) } ≥ 0. Following these definitions, and based on the approach in [9], the initial subsolution is chosen as: ζ̌0 = (0, 0, . . . , 0), if G i ≥ 0 for all i = 1, . . . , n. To compute the corresponding initial supersolution ζ̂0, we solve the nonlinear system: Biζ̂0 = G i(ζ̂0), for each i = 1, . . . , n. The discrete Hamilton-Jacobi-Bellman (HJB) problem, as expressed in equation (3) can be reformulated using an overlapping domain decomposition strategy, as follows: • Solution for the domain Ω1: To compute ζk+1 1 , we solve the following system: max 1≤i≤n ( Bi 1,1ζ k+1 1 − G i(ζk1 ) ) = 0, in Ω1,( ζk+1 1 ) β = ( ζk2 ) β , on Σ1. • Solution for the domain Ω2: To compute ζk+1 2 , we solve the following system: max 1≤i≤n ( Bi 2,2ζ k+1 2 − G i(ζk2 ) ) = 0, in Ω2,( ζk+1 2 ) α = ( ζk+1 1 ) α , on Σ2. Upon incorporating boundary conditions, the system becomes: max 1≤i≤n { Bi 1,1ζ k+1 1 − B̃i 1,2ζ k 2 − G i(ζk1 ) } = 0, in Ω1, max 1≤i≤n { Bi 2,2ζ k+1 2 − B̃i 2,1ζ k+1 1 − G i(ζk2 ) } = 0, in Ω2. We reformulate the iterative process derived from the overlapping domain decomposition method into the following matrix-based nonlinear system:( Bi 1,1 0 −B̃i 2,1 Bi 2,2 )( ζk+1 1 ζk+1 2 ) = ( B̃i 1,2 0 0 0 )( ζk1 ζk2 ) + ( G i(ζk1 ) G i(ζk2 ) ) . (4) To better understand this formulation, we recast it in the context of the block Gauss- Seidel method. The general form of this iterative scheme is given by: ζ(0) ∈ Rm, Bt,t ζ (k+1) t = − m∑ st Bt,s ζ (k) s + Gt, t = 1, . . . ,m. (5) S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 8 of 20 By comparing systems (4) and (5), it becomes clear that they are structurally equivalent in the nonlinear context, where the source term G is replaced by the nonlinear expression G k at each iteration. The algorithm below outlines the iterative steps of the proposed method. Algorithm Overlapping Domain Decomposition Initialize: Set the initial guess as ζ̌0 = (0, 0, . . . , 0), and set k = 0. Step 1: Compute the updated solution ζ̌k+1 1 by solving the system: max 1≤i≤n { Bi 1,1ζ̌ k+1 1 − B̃i 1,2ζ̌ k 2 − G i(ζ̌k1 ) } = 0. Step 2: Calculate the updated value ζ̌k+1 2 using the most recent value of ζ̌k+1 1 : max 1≤i≤n { −B̃i 2,1ζ̌ k+1 1 + Bi 2,2ζ̌ k+1 2 − G i(ζ̌k2 ) } = 0. Step 3: Verify convergence of the solution. If the solution has not converged, incre- ment k by 1 and repeat from Step 1. 3.4. Monotonic and Geometric Convergence To validate the proposed algorithm, we establish its monotonic and geometric con- vergence using lower and upper solutions. The following theorem presents the proof. Theorem 1. The iterative process described in Algorithm generates two sequences, (ζ̌k) and (ζ̂k), which exhibit monotonic behavior. More precisely, the sequence (ζ̌k) increases, while (ζ̂k) decreases, both converging toward the unique solution ζ of the Hamilton- Jacobi-Bellman system (3). Proof. The proof is structured into five consecutive steps. • Step 1: We establish that the sequence (ζ̌k), produced by the algorithm is mono- tonic and converges to the unique solution of the system (3). Initially, for k = 0, we solve the problem in the first subdomain to compute ζ̌11 as follows: max 1≤i≤n { Bi 1,1ζ̌ 1 1 − B̃i 1,2ζ̌ 0 2 − G i(ζ̌01 ) } = 0. Since ζ̌0 is a lower solution, the following inequality holds: max 1≤i≤n { Bi 1,1ζ̌ 0 1 − B̃i 1,2ζ̌ 0 2 − G i(ζ̌01 ) } ≤ 0. From this, we deduce that: ζ̌01 ≤ ζ̌11 . Next, we proceed to solve the problem in the second subdomain to determine ζ̌12 : max 1≤i≤n { −B̃i 2,1ζ̌ 1 1 + Bi 2,2ζ̌ 1 2 − G i(ζ̌02 ) } = 0. (6) S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 9 of 20 Since ζ̌0 is also a lower solution, we obtain: max 1≤i≤n { −B̃i 2,1ζ̌ 0 1 + Bi 2,2ζ̌ 0 2 − G i(ζ̌02 ) } ≤ 0. (7) Comparing equations (6) and (7), we deduce: ζ̌02 ≤ ζ̌12 . Thus, we conclude that: ζ̌0 ≤ ζ̌1. This confirms that ζ̌1 is also a lower solution. By mathematical induction, assuming that ζ̌k ≤ ζ̌k+1 holds for some k, we extend this property to all iterations: ζ̌0 ≤ ζ̌1 ≤ · · · ≤ ζ̌k ≤ ζ̌k+1. Hence, the sequence (ζ̌k) is non-decreasing. • Step 2: Establishing the convergence of the sequence (ζ̌k) to the solution of system (3). Since ζ̌k is considered a lower solution, the following inequality holds: max 1≤i≤n { Biζ̌k − G i(ζ̌k) } ≤ 0. (8) This directly implies the relation: Biζ̌k − G i(ζ̌k) ≤ 0 ⇒ Biζ̌k − F i − µζ̌k ≤ 0. (9) Now, let us assume the existence of a solution ζ̌∗ satisfying: Biζ̌∗ − F i − µζ̌∗ = 0. (10) From equation (10), we express: F i = Biζ̌∗ − µζ̌∗. By substituting this into equation (9), we obtain: Biζ̌k − Biζ̌∗ + µζ̌∗ − µζ̌k ≤ 0. Rearranging, we get: (Bi − µI)(ζ̌k − ζ̌∗) ≤ 0. Since Bi can be rewritten as A i + µI, it follows that: A i(ζ̌k − ζ̌∗) ≤ 0. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 10 of 20 Given that A i is an M-matrix, we deduce: ζ̌k ≤ ζ̌∗. Thus, the sequence (ζ̌k) is bounded below by ζ̌∗. Since (ζ̌k) is a monotonically increasing and lower-bounded sequence, it necessarily converges. Therefore, we conclude: lim k→∞ ζ̌k = ζ̌∗. (11) • Step 3: The vector ζ̌∗ = (ζ̌∗1 , ζ̌ ∗ 2 ) corresponds to the solution of system (3). From the iterative scheme described in algorithm, we derive the following condi- tions: max 1≤i≤n { Bi 1,1ζ̌ k+1 1 − B̃i 1,2ζ̌ k 2 − G i(ζ̌k1 ) } = 0, max 1≤i≤n { −B̃i 2,1ζ̌ k+1 1 + Bi 2,2ζ̌ k+1 2 − G i(ζ̌k2 ) } = 0. Taking the limit as k → ∞ and employing equation (11), we obtain the following system: max 1≤i≤n { Bi 1,1ζ̌ ∗ 1 − B̃i 1,2ζ̌ ∗ 2 − G i(ζ̌∗1 ) } = 0, max 1≤i≤n { −B̃i 2,1ζ̌ ∗ 1 + Bi 2,2ζ̌ ∗ 2 − G i(ζ̌∗2 ) } = 0. This system can be rewritten in a more compact form as: max 1≤i≤n { Biζ̌∗ − G i(ζ̌∗) } = 0. Consequently, we establish that ζ̌∗ is the unique solution to system (3). • Step 4: Uniqueness of the discrete solution to the Hamilton-Jacobi-Bellman sys- tem (3). Assume that the system (3) admits two possible solutions, denoted by ζ− and ζ+. This leads to the following two conditions: max 1≤i≤n { Biζ− − G i(ζ−) } = 0, max 1≤i≤n { Biζ+ − G i(ζ+) } = 0. From these relations, it follows that: Biζ+ − G i(ζ+) = 0 ⇔ Biζ+ − F i − µζ+ = 0, Biζ− − G i(ζ−) = 0 ⇔ Biζ− − F i − µζ− = 0. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 11 of 20 By subtracting these two equations and considering the decomposition Bi = A i+ µI, we obtain: (Bi − µI)(ζ+ − ζ−) = 0 ⇔ A i(ζ+ − ζ−) = 0. Given that the matrices A i are M-matrices, we conclude that: ζ− − ζ+ = 0 ⇒ ζ− = ζ+. Thus, the solution to the system (3) is uniquely determined. • Step 5: Demonstrating the monotonicity of the upper sequence (ζ̂k) follows an approach similar to that used for the lower sequence (ζ̌k). The goal here is to establish that the two sequences, (ζ̂k) and (ζ̌k), maintain an ordering relationship across all subdomains. Specifically, we aim to show: ζ̌k ≤ ζ̂k, ∀k ∈ N. For the initial iteration k = 0, we assume that ζ̌0 ≤ ζ̂0. In the first subdomain, since ζ̌1 is a lower solution, it satisfies the inequality: max 1≤i≤n { Bi 1,1ζ̌ 1 1 − B̃i 1,2ζ̌ 0 2 − G i(ζ̌01 ) } ≤ 0. (12) On the other hand, since ζ̂1 is an upper solution, it must satisfy: max 1≤i≤n { Bi 1,1ζ̂ 1 1 − B̃i 1,2ζ̂ 0 2 − G i(ζ̂01 ) } ≥ 0. (13) By comparing the inequalities in equations (12) and (13), we deduce that: max 1≤i≤n { Bi 1,1ζ̌ 1 1 − B̃i 1,2ζ̌ 0 2 − G i(ζ̌01 ) } ≤ max 1≤i≤n { Bi 1,1ζ̂ 1 1 − B̃i 1,2ζ̂ 0 2 − G i(ζ̂01 ) } . Thus, we obtain ζ̌02 ≤ ζ̂02 and ζ̌01 ≤ ζ̂01 , leading to: ζ̌11 ≤ ζ̂11 . Next, solving the problem in the second subdomain gives the following inequalities: max 1≤i≤n { −B̃i 2,1ζ̌ 1 1 + Bi 2,2ζ̌ 1 2 − G i(ζ̌02 ) } ≤ 0, (14) and max 1≤i≤n { −B̃i 2,1ζ̂ 1 1 + Bi 2,2ζ̂ 1 2 − G i(ζ̂02 ) } ≥ 0. (15) By comparing equations (14) and (15), we obtain: max 1≤i≤n { −B̃i 2,1ζ̌ 1 1 + Bi 2,2ζ̌ 1 2 − G i(ζ̌02 ) } ≤ max 1≤i≤n { −B̃i 2,1ζ̂ 1 1 + Bi 2,2ζ̂ 1 2 − G i(ζ̂02 ) } . S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 12 of 20 From the assumption ζ̌11 ≤ ζ̂11 and ζ̌02 ≤ ζ̂02 , we conclude that: ζ̌12 ≤ ζ̂12 . Therefore, we can conclude that: ζ̌1 ≤ ζ̂1. By applying induction, we prove that ζ̌k ≤ ζ̂k for all k. Since the sequence (ζ̌k) is increasing and bounded above, while the sequence (ζ̂k) is decreasing and bounded below, and given that these two sequences are interdependent, it can be concluded that both sequences converge to the unique solution of the discrete problem (3). Thus, we obtain: lim k→∞ ζ̂k = lim k→∞ ζ̌k = ζ Theorem 2. For any initial solution ζ0, the iterative sequence (ζk), generated by the algorithm converges geometrically to the unique solution of the system (3), with a positive constant 0 < η < 1 such that: ∥ζk+1 − ζk∥L∞(Ω) ≤ Cηk, k ∈ N. Proof. We begin by examining the formulation of the system and noting the following relations: max 1≤i≤n { Bi 1,1ζ k+1 1 − B̃i 1,2ζ k 2 − G i(ζk1 ) } = 0, and max 1≤i≤n { Bi 1,1ζ k 1 − B̃i 1,2ζ k−1 2 − G i(ζk−1 1 ) } = 0. These conditions correspond to solving the following problems: max 1≤i≤n ( Bi 1,1ζ k+1 1 − G i(ζk1 ) ) = 0, in Ω1,( ζk+1 1 ) β = ( ζk2 ) β , on Σ1, and  max 1≤i≤n ( Bi 1,1ζ k 1 − G i(ζk−1 1 ) ) = 0, in Ω1,( ζk1 ) β = ( ζk−1 2 ) β , on Σ1. By analyzing the boundary conditions of both problems and using the results from [10] applied to the Hamilton-Jacobi-Bellman equation, we deduce the existence of a constant κ1 ∈ (0, 1) such that: κ1 = sup{v(x) : x ∈ Σ1}. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 13 of 20 Thus, we obtain the following inequality: ∥ζk+1 1 − ζk1 ∥L∞(Σ1) ≤ κ1∥ζk2 − ζk−1 2 ∥L∞(Ω2). Similarly, analyzing the second part of the system, we have: max 1≤i≤n { −B̃i 2,1ζ k+1 1 + Bi 2,2ζ k+1 2 − G i(ζk2 ) } = 0, and max 1≤i≤n { −B̃i 2,1ζ k 1 + Bi 2,2ζ k 2 − G i(ζk−1 2 ) } = 0. By applying the results from [10], we conclude the existence of a constant κ2 ∈ (0, 1) such that: κ2 = sup{v(x) : x ∈ Σ2}. This implies the following relationship: ∥ζk+1 2 − ζk2 ∥L∞(Σ2) ≤ κ2∥ζk+1 1 − ζk1 ∥L∞(Ω1). By applying the maximum principle, we obtain: ∥ζk+1 2 − ζk2 ∥L∞(Ω2) ≤ ∥ζk+1 2 − ζk2 ∥L∞(Σ2), ∥ζk+1 2 − ζk2 ∥L∞(Σ2) ≤ κ2∥ζk+1 1 − ζk1 ∥L∞(Ω1), ≤ κ2∥ζk+1 1 − ζk1 ∥L∞(Σ1), ≤ κ1κ2∥ζk2 − ζk−1 2 ∥L∞(Ω2). Similarly, we derive: ∥ζk+1 1 − ζk1 ∥L∞(Ω1) ≤ ∥ζk+1 1 − ζk1 ∥L∞(Σ1), ∥ζk+1 1 − ζk1 ∥L∞(Σ1) ≤ κ1∥ζk2 − ζk−1 2 ∥L∞(Ω2), ≤ κ1∥ζk2 − ζk−1 2 ∥L∞(Σ2), ≤ κ1κ2∥ζk1 − ζk−1 1 ∥L∞(Ω1). Finally, we conclude with the inequality: ∥ζk+1 − ζk∥L∞(Ω) ≤ κ1κ2∥ζk − ζk−1∥L∞(Ω). By induction and noting that η = κ1κ2 ∈ (0, 1), we conclude: ∥ζk+1 − ζk∥L∞(Ω) ≤ ηk∥ζ1 − ζ0∥L∞(Ω). Thus, we obtain: ∥ζk+1 − ζk∥L∞(Ω) ≤ Cηk, k ∈ N. Therefore, the sequence (ζk) converges geometrically to the solution of the system (3). S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 14 of 20 4. Numerical Applications and Computational Results In this section, we formulate the Hamilton-Jacobi-Bellman (HJB) equation using an overlapping domain decomposition method. The computational domain Ω = [0, 1]×[0, 1] is discretized on a uniform grid with step size h based on the finite difference method. The initial solution can be selected as either a lower or upper solution, both ensuring convergence to the discrete solution. The iterative process continues until the stopping criterion ε = 10−6 is satisfied. The system is governed by the following equations:{ max { B1ζ − G 1,B2ζ − G 2 } in Ω, ζ = 0 on ∂Ω. (16) To facilitate the numerical treatment, we introduce the following notations: • B1 = B + µI, G 1 = G + µζ, with µ = 1. • B2 = I, where I is the identity matrix, and G 2 = 0. The operator B, which is non-coercive, is defined as: B = − ∂2 ∂x2 − ∂2 ∂y2 + 0.5x ∂ ∂x + 0.5y ∂ ∂y + 0.045. The function G is given by: G = sin(2πx) sin(2πy). To evaluate the efficiency of the domain decomposition methods, we compare the perfor- mance of the non-overlapping domain decomposition (τ = 1) (See [11]) with the newly introduced overlapping domain decomposition method (τ > 1). This comparative study aims to assess the accuracy and convergence of both methods, verifying their effective- ness against theoretical expectations. The numerical results, illustrating the efficiency of these approaches, are presented in tables and figures generated using MATLAB - R2018a. Table 1: Residual error using DDM for various overlap values. ∥R∥∞ iteration DDM: τ > 1 DDM: τ = 1 1 7.391048136608375e-04 9.591296770784198e-04 4 1.608036639652848e-04 2.363530936250824e-04 8 3.630680843518114e-05 5.836710081515498e-05 Last iteration 2.162254986382806e-05 3.666961301803263e-05 The table (1) presents the infinity norm of the residual error for the approximate solutions obtained using finite difference methods with h = 1 20 under various overlap configurations. The residual is computed as R = max 1≤i≤2 { Biζ̌k − G i ( ζ̌k−1 )} . S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 15 of 20 The results indicate that the residual errors are remarkably small and consistently de- crease across iterations. Notably, configurations with overlap (τ > 1) result in a faster reduction in error compared to the non-overlapping case (τ = 1). This emphasizes the importance of overlap in accelerating convergence and improving the accuracy of the solution approximation. Table 2: Value of ζ̂k at (x, y)T = (0.5, 0.25)T for h = 1 20 . Iteration DDM : τ > 1 DDM : τ = 1 1 -0.003915453822399 -0.003375478664288 4 -0.005478666922059 -0.005035539613455 8 -0.005493072486538 -0.005428621964170 Last iteration -0.005493166542087 -0.005493095359904 Table 3: Value of ζ̌k at (x, y)T = (0.5, 0.25)T for h = 1 20 . Iteration DDM : τ > 1 DDM : τ = 1 1 -0.137907986376806 -1.035248717772950 4 -0.012726797737372 -0.386671348701684 8 -0.005505924360092 -0.098756491226852 Last iteration -0.005493173019551 -0.005493282229003 The numerical results presented in Tables (2) and (3) illustrate the convergence properties of the lower and upper solutions for the Hamilton-Jacobi-Bellman (HJB) equation. As expected, the lower solution ζ̌k exhibits a monotonic increase, while the upper solution ζ̂k follows a monotonic decrease, ensuring a stable convergence towards the exact solution. Moreover, the comparison between the overlapping (τ > 1) and non-overlapping (τ = 1) domain decomposition methods confirms that increasing the overlap parameter significantly enhances the convergence rate, reducing the number of iterations required to reach the desired accuracy. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 16 of 20 Figure 2: Approximated Solutions (DDM). Figure (2) represents the numerical solution of the problem using the overlapping Schwarz domain decomposition method. The red and blue lines indicate the overlapping subdomains along the x-axis. Figure 3: Convergence Behavior of the Iterative Scheme. Figure (3) demonstrates the geometric convergence of the iterative method for h = 1 20 . The graph displays the evolution of the infinity norm of the error ∥ζk+1 − ζk∥∞ as a function of the iteration number. We observe a rapid exponential decrease in the error, indicating that the approximate solution is approaching a stable fixed point. Furthermore, using a larger overlap (blue curve) significantly improves the convergence rate by reducing both the number of iterations and the total computational effort. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 17 of 20 Figure 4: Residuals Evolution for the Overlapping Domain Decomposition . Figure (4) illustrates the residuals R1 and R2, representing the error distribution in the approximate solution while highlighting the active and inactive points at the boundaries. These plots demonstrate the error evolution through iterations within each subdomain, where these values are utilized to adjust and refine the solution at each step. Figure 5: Maximum Residual Over the Computational Domain. Figure (5) illustrates the distribution of the maximum residual over the computa- tional domain. The error is clearly concentrated in the overlapping region between the two subdomains Ω1 and Ω2, where it drops to a minimal value of approximately 10−5. This confirms the efficiency of the iterative method in reducing the residual specifically within the interface area, thereby enhancing the accuracy of the solution and accelerating convergence. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 18 of 20 5. Conclusion In this work, we solved the noncoercive Hamilton-Jacobi-Bellman equation directly, without transforming it into a quasi-variational inequality system. This approach pre- serves the structure of the original problem and simplifies the numerical treatment. The method is based on constructing lower and upper solutions that converge to the discrete solution within the finite difference framework. The overlapping domain decomposition method was applied, which contributed to improving the convergence rate compared to the non-overlapping case. Theoretical anal- ysis confirms that this method ensures monotonic and geometric convergence. Further- more, numerical results demonstrated that the performance achieved by this method is characterized by high accuracy and exhibits significant similarity to methods such as the Two-Level and Multi-Grid approaches, in terms of convergence speed and result precision. In future work, we plan to extend this method to q overlapping subdomains in higher- dimensional spaces, particularly for non-matching grids. Additionally, we plan to com- bine finite element discretizations with different boundary conditions, integrating the domain decomposition method with nonlinear multigrid techniques to further improve efficiency and scalability for large-scale problems. We will also investigate parallel im- plementations to enhance performance and tackle more complex problems. Acknowledgements The authors thank the reviewers and editors for their valuable feedback and support. References [1] L’ubomı́r Baňas, Herbert Dawid, Tsiry Avisoa Randrianasolo, Johannes Storn, and Xingang Wen. Numerical approximation of a system of Hamilton–Jacobi–Bellman equations arising in innovation dynamics. Journal of Scientific Computing, 92:54, 2022. [2] Dietmar Gallistl and Endre Süli. Mixed Finite Element Approximation Of The Hamilton–Jacobi–Bellman Equation With Cordes Coefficients. SIAM Journal on Numerical Analysis, 57(2):592–614, 2019. [3] Mohammed Beggas and Mohamed Haiour. The Maximum Norm Analysis of Schwarz Method for Elliptic Quasi-Variational Inequalities. Kragujevac Journal of Mathematics, 45(4):635–645, 2021. [4] L. Djemaoune and M. Haiour. Overlapping Domain Decomposition Method for a Noncoercive System of Quasi-Variational Inequalities Related to the Hamil- ton–Jacobi–Bellman Equation. Computational Mathematics and Modeling, 27(2), 2016. S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 19 of 20 [5] Madjda Miloudi, Samira Saadi, and Mohamed Haiour. Hamilton-Jacobi-Bellman Equations: An Algorithmic Contraction New Approach. Applied Mathematics & Information Sciences, 15(1):83–88, 2021. [6] Philippe Cortey-Dumont. Sur les inéquations variationnelles à opérateur non coercif. Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, 19(2):195–212, 1985. [7] P. G. Ciarlet and P. A. Raviart. Maximum principle and uniform convergence for the finite element method. Computer Methods in Applied Mechanics and Engineering, 2(1):17–31, 1973. [8] P. L. Lions and J. L. Menaldi. Optimal Control of Stochastic Integrals and Hamilton-Jacobi-Bellman Equations (Part I). SIAM Journal on Control and Opti- mization, 20:58–81, 1982. [9] M. Boulbrachene and M. Haiour. The finite element approximation of Hamilton- Jacobi-Bellman equations. Computers and Mathematics with Applications, 41(7- 8):993–1007, 2001. [10] Mohamed Haiour and Salah Boulaaras. Overlapping domain decomposition meth- ods for elliptic quasi-variational inequalities related to impulse control problem with mixed boundary conditions. Proceedings of the Indian Academy of Sciences (Math- ematical Sciences), 121(4):481–493, 2011. [11] M. Sun. Domain decomposition algorithms for solving hamilton-jacobi-bellman equations. Numerical Functional Analysis and Optimization, 14(1-2), 1993. [12] M. Boulbrachene. A contraction approach for noncoercive hamilton-jacobi-bellman equations. Applied Mathematics Information Sciences, 4(2):261–269, 2010. [13] M. Boulbrachene, M. Haiour, and B. Chentouf. On a Noncoercive System of Quasi- Variational Inequalities Related to Stochastic Control Problems. Journal of In- equalities in Pure and Applied Mathematics, 3(2):Article 30, 14 p, 2002. [14] Salah Boulaaras and Mohamed Haiour. A General Case for the Maximum Norm Analysis of an Overlapping Schwarz Methods of Evolutionary HJB Equation with Nonlinear Source Terms with the Mixed Boundary Conditions. Applied Mathemat- ics & Information Sciences, 9(3):1247–1257, 2015. [15] Mohamed Haiour and Salah Boulaaras. Uniform Convergence of Schwarz Method for Elliptic Quasi-Variational Inequalities Related to Impulse Control Problem. An- Najah University Journal of Research (N. Sc.), 24, 2010. [16] Salah Boulaaras, Mohamed El Amine Bencheikh Le Hocine, and Mohamed Haiour. A new error estimate on uniform norm of a parabolic variational inequality with nonlinear source terms via the subsolution concepts. Journal of Inequalities and Applications, 2020:78, 2020. [17] Ph. Cortey-Dumont. Sur l’analyse Numérique des Équations de Hamilton-Jacobi- Bellman. Mathematical Methods in the Applied Sciences, 9:198–209, 1987. [18] Lawrence C. Evans and Avner Friedman. Optimal Stochastic Switching And The Dirichlet Problem For The Bellman Equations. Transactions of the American Math- ematical Society, 253:365–389, 1979. [19] M. J. Gander and L. Halpren. Méthodes de décomposition de domaines - notions S. Chebbah, S. Madi, M. Haiour / Eur. J. Pure Appl. Math, 18 (2) (2025), 6070 20 of 20 de base. Sciences Fondamentales Mathematics, AF1375(V1), 2012. [20] Abida Harbi, Nasreddine Nemis, and Mohamed Haiour. Optimal error estimates of a class of system of two quasi-variational inequalities. AIMS Mathematics, 6(6):5977– 6002, 2021. [21] Ronald H. W. Hoppe. Multi-grid methods for Hamilton-Jacobi-Bellman equations. Numerische Mathematik, 49:239–254, 1986. [22] P. L. Lions and B. Mercier. Approximation numérique des Équations de Hamilton- Jacobi-Bellman. RAIRO. Analyse Numérique, 14(4):369–393, 1980. [23] Madjda Miloudi, Samira Saadi, and Mohamed Haiour. L∞-error estimates of a fi- nite element method for Hamilton–Jacobi–Bellman equations with nonlinear source terms with mixed boundary condition. Demonstratio Mathematica, 54(1):452–461, 2021. [24] P. Spiteri, J. C. Miellou, and D. El Baz. Asynchronous Schwarz alternating meth- ods with flexible communication for the obstacle problem. Calculateurs Parallèles, Réseaux et Systèmes Répartis, 13(1):47–66, 2001. Introduction Coercive Reformulation of the HJB Equation Decomposition of the Domain into Overlapping Subdomains Definition of Subdomain Decomposition Block Decomposition of the Domain Matrix Iterative Algorithm for Domain Decomposition Monotonic and Geometric Convergence Numerical Applications and Computational Results Conclusion