Convergence analysis for a nonlinear system of parabolic variational inequalities

This work aims to provide a comprehensive and unified numerical analysis for a nonlinear system of parabolic variational inequalities (PVIs) subject to Dirichlet boundary condition. This analysis enables us to establish the existence of an exact solution to the considered model and to prove the convergence for the approximate solution and its approximate gradient. Our results are applicable for several conforming and nonconforming numerical schemes.


Introduction
Nonlinear parabolic variational inequalities and PDEs are useful tools to model the coupled biochemical interactions of microbial cells, which are crucial to numerous applications, especially in the medical field and food production [16,19,21]. We consider here a nonlinear parabolic system consisting of PDEs and variational inequalities: Numerical approximation in parabolic systems of inequalities and generalization of inequalities have received considerable attention in the research literature. Wheeler [23] obtains the error estimate of second order in L ∞ (L 2 ) for a linear approximation with respect to space and time, with a strong regularity on the solutions, such as ∂ tB ∈ L 2 (0, T, L 2 ( )).
Johnson [17] analyzes inequality (1.1a), in which F = 0 and D A is constant. The work in [4] considers a model without the barrier and provides O(h) order of convergence in L ∞ (L 2 )norm. An L 2 -error estimate is provided in different studies, such as [5], by using a finite difference in time. Vuik [22] deals with parabolic variational inequalities with a nonlinear source term and derives the convergence rate of the finite element method in space with respect to L ∞ -norm. He also shows that the general finite difference gives O(h) in L ∞ (L 2 )-norm under a strong hypothesis on data. Saker et al. [20] discuss the discrete and continuous forms of a Carlson-type inequality, and [18] introduces Minkowski's inequality by using AB-fractional integral operators.
However, there is a lack of full convergence analyses of numerical schemes for the model (1.1a)-(1.1f) since the coupled nonlinearity of the system and the constraint (the inequality) in the model comprise the primary theoretical challenge. It appears that considerable research is still required, beginning with convergence analysis and testing other varieties of scheme outside conforming methods. Rather than undertaking individual research for every numerical scheme, this work utilizes a gradient discretization method (GDM) to provide a unified and full convergence analysis of numerical methods for (1.1a)-(1.1f) under natural hypotheses on data. The GDM is a generic framework to unify the numerical analysis for diffusion partial differential equations and their corresponding problems. Due to the variety of choice of the discrete elements in the GDM, a series of conforming and nonconforming numerical schemes can be included in the GDM, see [2,3,6,[9][10][11][12][13] for more details.
The outline of this paper is as follows. Section 2 is devoted to writing the model (1.1a)-(1.1f) in an equivalent weak sense. Section 3 defines the discrete space and functions followed by the gradient scheme to our model in the weak sense. Section 4 provides the convergence results, Theorem 4.5, which is proved by following the compactness technique under classical hypothesis on continuous model data. Finally, as an example, we present in Sect. 5 the nonconforming P1 finite element scheme that has not been applied to the nonlinear model (1.1a)-(1.1f), yet.
It is clear that the time-dependent K contains at least the constant in time function t → χ -:= min(0, χ).

Discrete setting
We begin with defining the discrete space and operators. These discrete elements are slightly different from those defined in [2,3], in particular, χ D , I D , and J D are introduced to deal with the nonconstant barrier χ and the initial solutions A ini and B ini .
Remark 3.2 For a general obstacle χ , most of numerical methods fail to approximate the solutionĀ by elements inside the set K D . For an example, in the P1 finite element method, we consider only the values ofĀ at the vertices of the mesh, which only guarantee that these values satisfy the barrier condition (1.1c) only at these vertices, not necessarily at any point in . We define here the set K D based on the approximate barrier ψ D to be able to construct an interpolant that belongs to K D . However, there is no need to use an approximate barrier, if the barrier χ is assumed to be constant.
In order to construct a good approximate scheme, we require four properties: coercivity, consistency, limit-conformity, and compactness. The first three respectively connect to the Poincaré inequality, the interpolation error, and the Stokes formula. The compactness property enables us to deal with the nonlinearity caused by the reaction terms F and G.

Main results
Let the time interval [0, T] be divided into κ intervals of length κ, where κ tends to zero as κ → ∞. Let 1 I i be the characteristic function of I i = [iκ, (i + 1)κ), i = 0, . . . , κ . We define a set of piecewise-constant in time functions by we obtain, for all s ∈ (0, T) and a.e. x x x ∈ ,

An application of the definition of S D m yields
. Using consistency, one obtains S D m (φ i ) → 0 as m → ∞, for any i = 0, . . . , κ , which implies that the second term on the right-hand side vanishes. In the case in which both s, t (n(s)+1) ∈ I i or both s, t (n(s)+1) / ∈ I i , the quantity ξ i m (s) -1 I i (s) equals zero. In the case in which s ∈ I i and t (n(s)+1) / ∈ I i or s / ∈ I i and t (n(s)+1) ∈ I i , one can deduce, This shows that the first term on the right-hand side of (4.4) tends to zero when m → ∞. Hence, (4.2a) is concluded. The proof of (4.2b) is obtained by the same reasoning, replacingw κ by ∇w κ and D m w m by ∇ D m w m .

Lemma 4.2 (Energy estimates) Let Hypothesis
(4.5) Proof We start by taking ϕ := A (n) (it belongs to K D ) and the function ψ := δt (n+ 1 2 (4.7) Applying the fact that (rs) · r ≥ 1 2 |r| 2 -1 2 |s| 2 to the second term on the left-hand side of (4.6) and to the first term on the left-hand side of (4.7), it follows that Summing the above inequalities over n ∈ [0, m -1], where m = 0, . . . , N gives This, together with the Cauchy-Schwarz inequality, implies that Using the Lipschitz condition, we arrive at This, together with Young's inequality, gives, whenever and Thanks to Gronwall inequality [15, Lemma 5.1], inequality (4.11) can be rewritten as where C 3 depends on T, M, and ε 4 . Combining this inequality with (4.10) yields Taking the supremum over m ∈ [0, N] and using the real inequality sup n (r n + s n ) ≤ sup n (r n ) + sup n (s n ), we obtain the desired estimates.
In the following definition, we introduce a dual norm [8], which is defined on the space D (X D,0 ) ⊂ L 2 ( ), to ensure the required compactness results.

Definition 4.3
If D is a gradient discretization, then the dual norm · ,D on D (X D,0 ) is given by Proof Putting ψ = φ in (3.4b), together with the Cauchy-Schwarz inequality and the coercivity property, implieŝ Taking the supremum over φ ∈ X D,0 with ∇ D φ L 2 ( ) d = 1, multiplying by δt (n+1) , summing over n ∈ [0, N -1], and using (4.5) yield the desired estimate. Proof The proof is divided into four stages and its idea is inspired by [1].
Step 1: Existence of approximate solutions. At (n +1), we see that (3.4a) and (3.4b) respectively express a system of nonlinear elliptic variational inequality on A (n+1) and nonlinear equations on B (n+1) . For w = (w 1 , where α := 1 δt (n+ 1 2 ) and the bilinear and linear forms are defined by Stampacchia's theorem implies that there existsĀ ∈ K D satisfying inequality (4.14a). The second equation (4.14b) describes a linear square system. Taking ϕ = B (n+1) in (4.14b), using the similar reasoning as in the proof of Lemma 4.2, and setting G = 0 yield ∇ D B (n+1) L 2 ( ) d = 0. This shows that the matrix corresponding to the linear system is invertible. Consider the continuous mapping T : (A, B) with (A, B) being the solution to (4.14a)-(4.14b). The existence of a solution (A (n+1) , B (n+1) ) to the nonlinear system is a consequence of Brouwer's fixed point theorem.
Let us show the strong convergence of D m B m toB in L ∞ (0, T; L 2 ( )). Indeed, D m B m converges strongly toB in L 2 ( × (0, T)), thanks to estimate (4.13), consistency, limitconformity, and compactness, as well as [8,Theorem 4.14]. We can apply the dominated convergence theorem to show that ( × (0, T)), thanks to the assumptions on G given in Hypothesis 2.1.
Step 4: Proof of the strong convergence of ∇ D m A m and ∇ D m B m . From the weak-strong convergences, we have, for allw κ ∈ L κ , Thanks to the density results, for any ϕ ∈ K, we can find (w κ ) κ>0 that converges to ϕ in L 2 (0, T; H 1 ( )), as κ → 0. Therefore, we infer, for all ϕ ∈ K, Taking ϕ =Ā, the above relation yields Together with (4.20), we obtain showing that ∇ D m A m → ∇Ā strongly in L 2 ( × (0, T)) d . To show the strong convergence of ∇ D m B m , we begin by writinĝ Passing to the limit in (4.22) and using the above inequality, we reach the desired convergence result.

An example of schemes covered by the analysis
Many numerical schemes fit into the analysis provided in this work. We construct here the nonconforming P1 finite element scheme for our model. Let T = (M, E) be the polytopal mesh of defined in [8,Definition 7.2], in which M and E consist of cells K and edges σ , respectively. The elements of gradient discretization D associated with the nonconforming P1 finite element scheme are: • X D,0 = {w = (w σ ) σ ∈E : w σ ∈ R and w σ = 0 for all σ ∈ E ext }.
• For all w ∈ X D,0 and for all K ∈ M, for a.e. x x x ∈ K , where e σ K is a basis function. • For all w ∈ X D,0 and all K ∈ M, for a.e. x x x ∈ K , • The approximate obstacle χ D is defined by • For all ω ∈ W 2,∞ ( ), we can construct the interpolants I D ω = J D ω = (z σ ) σ ∈E with z σ = ω(x σ ). Substituting these elements into scheme (3.4a)-(3.4b) yields the nonconforming P1 finite element scheme for problem (2.1a)-(2.1b) and its convergence is therefore obtained from Theorem 4.5. Droniou et al. [7] show that D given here satisfies the three properties, namely coercivity, limit-conformity, and compactness. Let us discuss the consistency property in the sense of Definition 3.4. It is shown in [7] that S D m (ψ) → 0, for all ψ ∈ H 1 0 ( ), which verifies the second item. Similarly, we can prove that S D m (ϕ) → 0 for all ϕ ∈ C 2 ( ) ∩ K. For any ϕ, let ω = (ω σ ) σ ∈E ∈ X D,0 be the interpolant such that ω σ = ϕ(x σ ), for all σ ∈ E. We clearly deduce D ϕ ≤ χ D in . By the density results established in [14], we see that the first item is fulfilled.
Let ϕ m = (ω σ ) σ ∈E m ∈ K D m and ψ m = (w σ ) σ ∈E m ∈ X D m ,0 be the interpolants such that ϕ m = I D m A ini and ψ m = J D m B ini . Now [8, (B.11) in Lemma B.7] with p = 2 shows that there exist C 5 , C 6 > 0 not depending on m such that Passing to the limit, we see the right-hand sides tend to 0 (thanks to the classical regularity hypothesis on A ini and B ini ), and therefore the third and fourth items of the consistency property are verified. Finally, it is established in [8, proof of Theorem 12.12] that, for ϕ ∈ W 1,p ( ), we can construct a function ω m = (ω σ ) σ ∈E m ∈ K D m and find a C 7 > 0 not depending on m such that Applying this estimate (with p = 2) to ϕ = A ini and ω m = I D m A ini , we deduce that ∇ D m I D m A ini L 2 ( ) d is bounded.