A posteriori verification of the positivity of solutions to elliptic boundary value problems

The purpose of this paper is to develop a unified a posteriori method for verifying the positivity of solutions of elliptic boundary value problems by assuming neither $H^2$-regularity nor $ L^{\infty} $-error estimation, but only $ H^1_0 $-error estimation. In [J. Comput. Appl. Math, Vol. 370, (2020) 112647], we proposed two approaches to verify the positivity of solutions of several semilinear elliptic boundary value problems. However, some cases require $ L^{\infty} $-error estimation and, therefore, narrow applicability. In this paper, we extend one of the approaches and combine it with a priori error bounds for Laplacian eigenvalues to obtain a unified method that has wide application. We describe how to evaluate some constants required to verify the positivity of desired solutions. We apply our method to several problems, including those to which the previous method is not applicable.


Introduction
In recent decades, numerical verification (also known as computer-assisted proof, validated numerics, or verified numerical computation) has been developed and applied to various partial differential equations, including those where purely analytical methods have failed (see, for example, [28,5,29,23,26,24,31,25,33] and the references therein). One such successful application is to the semilinear elliptic boundary value problem where Ω ⊂ R N (N = 2, 3, · · · ) is a given bounded domain, ∆ is the Laplacian, and f : R → R is a given nonlinear map. Further regularity assumptions for Ω and f will be shown later for our setting. Positive solutions of (1) have attracted significant attention [19,10,18,4,11,6]. For example, positive solutions of problem (1) with f (t) = λt + t|t| p−1 , λ ∈ [0, λ 1 (Ω)), p ∈ (1, p * ) have been investigated from various points of view such as uniqueness, multiplicity, nondegeneracy, and symmetry, [10,18,4,11,6], where p * = ∞ when N = 2 and p * = (N + 2)/(N − 2) when N ≥ 3; λ 1 (Ω) is the first eigenvalue of −∆ with the homogeneous Dirichlet boundary value condition in the weak sense. Another important nonlinearity is f (t) = λ(t−t 3 ), λ > 0. This corresponds to the stationary problem of the Allen-Cahn equation [1]. The problem (1) with this nonlinearity may have a positive solution when λ ≥ λ 1 (Ω). However, when λ < λ 1 (Ω), no positive solution is admitted; this can be confirmed by multiplying −∆u = λ(u − u 3 ) with the first eigenfunction of −∆ and integrating both sides. The Allen-Cahn equation is a special case of the Nagumo equation [22] with the nonlinearity f (t) = λt(1 − t)(t − a), λ > 0 and 0 < a < 1, and both of these equations have been investigated by many researchers. We are moreover interested in the related case in which f (t) = λ(t + At 2 − Bt 3 ) with A, B > 0. The bifurcation of problem (1) with this nonlinearity was analyzed in [19]. This problem has two positive solutions when λ * < λ < λ 1 (Ω) for some λ * > 0. Despite these results, quantitative information about the positive solutions, such as their shape, has not been clarified analytically. Throughout this paper, H k (Ω) denotes the kth order L 2 Sobolev space. We define H 1 0 (Ω) := {u ∈ H 1 (Ω) : u = 0 on ∂Ω}, with the inner product (u, v) H 1 0 := (∇u, ∇v) L 2 and norm u H 1 0 := (u, u) H 1 0 . Numerical verification methods enable us to obtain an explicit ball containing exact solutions of (1). More precisely, for a numerical approximation u ∈ H 1 0 (Ω) that satisfies the assumption required by such methods, they prove the existence of an exact solution u ∈ H 1 0 (Ω) of (1) that satisfies with an explicit error bound ρ > 0. Under an appropriate condition, we can obtain an L ∞ -estimation with bound σ > 0. For instance, when u,û ∈ H 2 (Ω), we can evaluate the L ∞bound σ > 0 by considering the embedding H 2 (Ω) → L ∞ (Ω) [28, Theorem 1 and Corollary 1]. Thus, these approaches have the advantage that quantitative information about a target solution is provided accurately in a strict mathematical sense. We can identify the approximate shape of solutions from the error estimates. Despite these advantages, information about the positivity of solutions is not guaranteed without further considerations, irrespective of how small the error bound (ρ or σ) is. In the homogeneous Dirichlet case (1), it is possible for a solution that is verified by such methods to be negative near the boundary ∂Ω. Therefore, we developed methods of verifying the positivity of solutions of (1) in previous studies [35,37,36,33] and applied these result to sing-changing solutions [34]. These methods succeeded in verifying the existence of positive solutions by checking simple conditions. In [35,37,36], we proposed methods for verifying the positivity of solutions of (1) by assuming both error estimates (2) and (3). Subsequently, in [33], we extended our method to a union of two different approaches [33, Theorems 2.1 and 3.2] under certain conditions for nonlinearity f . Table 1 summarizes the error-estimate types that are required by these theorems when f is a subcritical polynomial Theorem 2.1 in [33] can be applied to cases in which λ < λ 1 (Ω). This theorem is based on the constructive norm estimation for the minus part u − := max{−u, 0} of a solution u, and does not assume an L ∞ -estimation (3) but only requires an H 1 0 -error estimation (2). When a i ≤ 0 for all i and λ ≥ λ 1 (Ω), we used a completely different approach [33,Theorem 3.2] that was based on the Newton iteration that retains nonnegativity. This theorem needs no L ∞estimation but requires an explicit evaluation of the minimal eigenvalue of a certain linearized operator around approximationû. Actually, this eigenvalue evaluation itself is not trivial to obtain. In [33], the eigenvalue was estimated using the existing method [33,20] based on Galerkin approximations. Theorems 2.1 and 3.2 in [33] were applied numerically to problem (1) with the above-mentioned nonlinearities: f (t) = λt + t|t| p−1 and f (t) = λ(t − t 3 ). However, a problem still remains in the sense that we need the L ∞ -error estimation (3) to prove positivity when a i a j < 0 for some i, j and λ ≥ λ 1 (Ω). For example, the nonlinearity f (t) = λ(t + At 2 − Bt 3 ) in which we are interested requires such an estimation. These requirements may narrow the applicability of the methods because the existing approach provided in [28, Theorem 1 and Corollary 1] requires u,û ∈ H 2 (Ω), evaluating the bound for the embedding H 2 (Ω) → L ∞ (Ω), to obtain σ from ρ. The H 2 -regularity of u may fall when Ω is a nonconvex polygonal domain.
The purpose of this paper is to develop a unified method for verifying the positivity of solutions u of (1) by assuming neither H 2 -regularity nor (3), but only H 1 0 -error estimation (2), which can be applied to all the cases in Table  Table 1 Error estimates required in [33] when f is the subcritical polynomial (4). 1 for arbitrary bounded domains Ω. Our method is based on a posteriori constructive norm estimation for the minus part u − and can be regarded as an extension of [33,Theorem 2.1]. In short, we confirm that the norm of u − vanishes by checking certain inequalities while assuming (2) (see Lemma 6.1). One of the key points is to estimate lower bounds of eigenvalue λ 1 (supp u − ) explicitly because the inequality λ 1 (supp u − ) ≥ λ has to be confirmed for the success of our positivity proof (again, see Lemma 6.1). Here, supp u − := {x ∈ Ω : u − (x) = 0} is the support of u − , and λ 1 (supp u − ) is understood as the first eigenvalue on the interior of supp u − . When the interior of supp u − is empty, we interpret λ 1 (supp u − ) = ∞, and all real numbers are lower bounds of λ 1 (supp u − ). The difficulty is that we cannot identify the location and shape of supp u − from (2) even whenû is nonnegative. If a polygon or polyhedron S enclosing supp u − is obtained concretely, we can apply the Liu-Oishi method [21,20] based on finite element methods to obtain a lower bound for λ 1 (S). Then, the inequality λ 1 (supp u − ) ≥ λ 1 (S) gives the desired lower bound of λ 1 (supp u − ). Such a supremum set S over supp u − can be obtained when we have an L ∞ -estimation (3). By setting Ω + := {x ∈ Ω :û − σ ≥ 0} where u ≥ 0 therein, we can construct such a domain S as a supremum set of Ω\Ω + . Again, we cannot determine such a supremum set S only from H 1 0 -error estimation (2). The Temple-Lehmann-Goerisch method can help us to evaluate λ 1 (supp u − ) more accurately (see, for example, [25,Theorem 10.31]).
To estimate the lower bound of λ 1 (supp u − ), we rely on the following argument: Fact 1.1 For a bounded domain Ω ⊂ R N (N = 2, 3, · · · ), there exists a constant A k,N independent of Ω such that where λ k (Ω) denotes the k-th eigenvalue of −∆ with the homogeneous Dirichlet boundary condition.
We also refer to [17] for an easy-to-estimate formula for A k,N for all k ≥ 1 (see Remark 4.3). Section 4 provides explicit lower bounds of A k,N based on these results. To estimate a lower bound of λ 1 (supp u − ) using the inequality (5), we focus on estimating upper bounds of |supp u − | while assuming only the H 1 0 -error estimation (2) without knowing the specific shape and location of supp u − . Suppose that (2) is proved for a positive approximationû with sufficient accuracy. Then, the upper bound for |supp u − | can be estimated very small using Lemma 3.2 provided later, and therefore, λ 1 (supp u − ) ≥ λ can be confirmed using (5) for a moderately large λ. The established estimation for |supp u − | is used to evaluate not only λ 1 (supp u − ) but also some Sobolev embedding constants on |supp u − |, which play an essential role for our positivity proof (again, see Lemma 6.1) The remainder of this paper is organized as follows. Section 2 introduces required notation and definitions. In Section 3, we evaluate an upper bound for the volume |supp u − | assuming the H 1 0 -estimation (2) for a continuous or piecewise continuous approximationû ∈ H 1 0 (Ω). Subsequently, we use the bound for |supp u − | to evaluate lower bounds for λ 1 (supp u − ) in Section 4. In Section 5, required Sobolev embedding constants on bounded domains are evaluated. In Section 6, we extend the previous formula [33, Theorem 2.1] and combine it with the estimates derived from Sections 4 and 5, thereby designing a unified method for proving positivity. Finally, Section 8 presents numerical examples where the proposed method is applied to problem (1) with several nonlinearities, including those to which the previous method is not applicable without an L ∞ -error estimation.

Preliminaries
We begin by introducing required notation. We denote by H −1 the topological dual of H 1 0 (Ω). When describing norms and inner products, we may omit the domain of a function space unless there is a risk of misunderstanding. For example, we simply write · L p = · L p (Ω) if no confusion arises. For two Banach spaces X and Y , the set of bounded linear operators from X to Y is denoted by L(X, Y ) with the usual supremum norm T L(X,Y ) := sup{ T u Y / u X : u ∈ X \ {0}} for T ∈ L(X, Y ). The norm bound for the embedding H 1 0 (Ω) → L p+1 (Ω) is denoted by C p+1 (= C p+1 (Ω)); that is, C p+1 is a positive number that satisfies where p ∈ [1, ∞) when N = 2 and p ∈ [1, p * ] when N ≥ 3. If no confusion arises, we use the notation C p+1 to represent the embedding constant on the entire domain Ω, whereas, in some parts of this paper, we need to consider an embedding constant on some subdomain Ω ⊂ Ω. This is denoted by C p+1 (Ω ) to avoid confusion. Moreover, λ 1 (Ω) denotes the first eigenvalue of −∆ imposed on the homogeneous Dirichlet boundary condition. This is characterized by Throughout this paper, we assume that f is a C 1 function that satisfies for some a 0 , a 1 , b 0 , b 1 ≥ 0 and p ∈ [1, p * ). We define the operator F as Moreover, we define another operator F : where The Fréchet derivatives of F and F at ϕ ∈ H 1 0 (Ω), denoted by F ϕ and F ϕ , respectively, are given by Under the notation and assumptions, we look for solutions u ∈ H 1 0 (Ω) of which corresponds to the weak form of (1). We assume that some verification method succeeds in proving the existence of a solution u ∈ H 1 0 (Ω) of (13) satisfying inequality (2) givenû ∈ H 1 0 (Ω) and ρ > 0. Although the regularity assumption forû (to be in H 1 0 (Ω)) is sufficient to obtain the error bound (2) in theory, we further assume thatû is continuous or piecewise continuous throughout this paper. This assumption impairs little of the flexibility of actual numerical computation methods. We recall u − = max{−u, 0}, and define u + := max{u, 0}.

Evaluation of the volume of supp u −
To estimate an upper bound for |supp u − | from the information of the inclusion (2), we define D(v) := {x ∈ Ω :û(x) ≤ v(x)} and consider the maximization problem for fixed q ∈ (1, ∞) and c > 0. This maximization takes place over the set of all functions v ∈ L q (Ω) satisfying v L q (Ω) = c. When û + L q (Ω) ≤ c, the maximal value of the problem (14) is |Ω|. Therefore, we consider the case where û + L q (Ω) > c. In the following, we denote D(l) where D is a set that satisfies û + L q (D) = c and for some l ∈ R. The maximal value of the problem (14) is |D|.
Let v c be nonnegative. Ifû−v c is strictly negative in some part Ω ⊂ supp v c satisfying |Ω | = 0 and vanishes in (supp v c )\Ω , another v c ∈ {v ∈ L q (Ω) : v L q (Ω) = c} with the same c > 0 can be constructed to obtain larger D(v c ) as follows.
Finally, we consider a subset D (= supp v c ) ⊂ Ω with the largest volume satisfying û + L q (D) = c. In the following, we prove that the volume of such D is maximized when (16) holds for some l ∈ R. Sinceû is continuous or piecewise continuous on Ω, there exist D and l ∈ R that satisfy û + L q (D) = c and (16).
Then, we have thatû and therefore, |D| ≥ |D |. Thus, the assertion of this lemma is proved.

Lemma 3.1 ensures that this maximal value is realized when
where D is a set that satisfies û + L q (D) = c and (16) for some l ∈ R. The maximal value is |D|, which is not greater than |D(l)| due to (16). Inequality (19) ensures that Suppose that D(m) ⊂ D(l) in the strict sense. Then, we have m < l, and thus, D(m) ⊆D(l) ⊆ D due to (16). This contradicts (21). Therefore, we have D(l) ⊆ D(m) and conclude (20).
The choice of q does not greatly affect realizing (19), and we can usually set q = 2. Meanwhile, appropriately setting m is important for confirming (19) as discussed below: Let q, ρ, andû be fixed. Then, û + L q (D(m)) monotonically decreases as m decreases, and û + L q (D(m)) ↓ 0 as m ↓ 0. Therefore, although smaller m gives a good upper bound of |supp u − | as in (20), too small m > 0 leads to failure in ensuring (19). Some concrete choices of m can be found in our numerical experiments in Section 8.

Lower bound for the minimal eigenvalue
The purpose of this section is to estimate a lower bound of the minimal eigenvalue λ 1 (supp u − ) while we assume an H 1 0 -estimation (2) only; therefore, supp u − cannot be identified explicitly. To this end, we use the following Rayleigh-Faber-Krahn constant.
where B N = π N/2 /Γ (N/2 + 1) denotes the volume of the unit N -ball with the usual gamma function Γ , and j N 2 −1,1 is the first positive zero of the Bessel function of order N 2 − 1. The equality in (5) is attained if and only if Ω is a ball in R N .
In general, evaluating A 1,N in explicit decimal form using (22) is not trivial. However, one can find in Table 2 rigorous enclosures of B N , j N 2 −1,1 , and A 1,N for several dimensions N . These were derived by strictly estimating all numerical errors; therefore, the correctness is mathematically guaranteed in the sense that correct values are included in the corresponding closed intervals. The enclosures were obtained using the kv library [12], a C++ based package for rigorous computations. The kv library includes four interval arithmetic operations and a function for rigorously calculating the gamma functions needed to derive B N . However, no function for enclosing the Bessel function j N 2 −1,1 is built therein. Accordingly, we present a rigorous algorithm for calculating j N 2 −1,1 in Appendix A based on the bisection method. Combining Lemma 3.2 and Theorem 4.1 for N = 2, 3 and using the lower bound in Table 2, we immediately have the following lower bounds for the minimal eigenvalue on supp u − .
This estimation is somewhat rough compared with that in Corollary 4.2 but stands alone in the sense that the lower bound can be calculated by hand as long as we know B N .

Embedding constant
Explicitly estimating the embedding constant C p is important for our method. We use [36, Corollary A.2] to obtain an explicit value of C p for bounded domains based on the best constant in the classical Sobolev inequality provided in [2,32].
Here, T p,N is defined by where Γ is the gamma function and q = N p/(N + p).
Remark 5.2 Another formula to estimate the embedding constant C p can be found in [25,Lemma 7.10], which is applicable not only to bounded domains but also to unbounded domains with a more generalized norm in H 1 0 (Ω). Moreover, in [36], one can find very sharp estimations of the best values of C p for p = 3, 4, 5, 6, 7 on Ω = (0, 1) 2 . Table 3 shows the strict upper bounds of the embedding constants for several cases. These were evaluated using MATLAB 2019a with INTLAB version 11 [30] with rounding errors strictly estimated; the required gamma functions were strictly computed via the function "gamma" packaged in INTLAB.  In the case of p = 2, to which Theorem 5.1 is inapplicable, the following best evaluation can be used instead of Theorem 5.1: For example, when Ω = (0, 1) N (N = 1, 2, 3, · · · ), we have the exact eigenvalue λ 1 (Ω) = N π 2 . Even otherwise, lower bounds of λ 1 (Ω) (and therefore upper bounds for the embedding constant) can be estimated for bounded domain Ω using the formulas in Section 4.

A posteriori verification for positivity
In this section, we design a unified method for proving positivity. We first clarify the assumption imposed on nonlinearity f with some explicit parameters.
for some λ ∈ R, nonnegative coefficients a 1 , a 2 , · · · , a n , and subcritical exponents p 1 , p 2 , · · · , p n ∈ (1, p * ). This assumption does not break the generality of f satisfying (8). Our algorithm for verifying positivity is based on the following argument, a generalization of [33, Thorem 2.1]. In the following, λ 1 (supp u − ) and C p+1 (supp u − ) are interpreted as λ 1 (int (supp u − )) and C p+1 (int (supp u − )), respectively, where int (supp u − ) denotes the interior of supp u − . We denote by λ ∈ R a lower bound of λ 1 (supp u − ). When the interior of supp u − is empty, we can set λ to an arbitrarily large value.
Lemma 6.1 Suppose that a solution u ∈ H 1 0 (Ω) of (13) exists in B(û, ρ), given some approximationû ∈ H 1 0 (Ω) and radius ρ > 0. If then u is nonnegative. Note that when supp u − is disconnected, (28) is understood as the set of inequalities for all connected components of supp u − .

Remark 6.2
The maximum principle ensures the positivity of u (i.e., u(x) > 0 inside Ω) from its nonnegativity for a wide class of nonlinearities f . See, for example, [8] for a generalized maximum principle applicable for weak solutions.
By applying the estimations in Sections 4 and 5 to Lemma 6.1, we have the following theorem. The constants in this theorem can be computed explicitly using formulas provided before.
Theorem 6.5 Suppose that a solution u ∈ H 1 0 (Ω) of (13) exists in B(û, ρ), given some approximationû ∈ H 1 0 (Ω) and radius ρ > 0. Moreover, we assume (19), given q ∈ [2, p * + 1) and m > 0. Then, we define C 1 = C 1 (N, f,û, ρ, m) and C 2 = C 2 (N, f,û, m) as where T i := T pi+1,N is defined by (24) and A 1,N is the constant of (22). If C 1 < C 2 , u is nonnegative. Table 4 shows calculation examples of C 1 and C 2 for some concrete nonlinearities f , where we can use the estimations of T p,N in Tables 2 and 3. For the first nonlinearity f (t) = λt + |t| p−1 t with a subcritical exponent p ∈ (1, p * ), positive solutions are admitted only when λ < λ 1 (Ω). Therefore, calculating |D(m)| can be avoided if we can evaluate λ 1 (Ω) with sufficient accuracy. Note that λ 1 (Ω) can be obtained analytically, such as when Ω is a hyperrectangle. When λ = 0, |D(m)| does not need to be calculated to derive C 2 . In cases where we do not require |D(m)| to estimate C 2 , it is reasonable to replace |D(m)| with |Ω| in calculating C 1 to reduce calculation cost especially when ρ is sufficiently small. Table 4 Calculation examples of C 1 and C 2 for several nonlinearities f .
The values of the gamma functions should be estimated explicitly to compute T p+1,N . There are several toolboxes that enable us to calculate the gamma functions rigorously. Indeed, as mentioned in Sections 4 and 5, kv library [12] and INTLAB [30] have such a rigorous function. Therefore, we are left to estimate upper bounds for û − L p and |D(m)|, and a lower bound for û + L p (D(m)) . We describe the ways to calculate these bounds in the following subsections. For sufficient accuracy, we divide Ω into a union of small We can obtain such a division "for free" such as when using finite element methods to computeû. Otherwise, we should create a mesh {K i } N K i=1 of Ω that satisfies the above property. When Ω is polygonal, a convenient way is to divide Ω into a rectangular or triangular mesh. In preparation, we estimate a lower bound To calculate an upper bound of û − L p , we use the inequality Upper bounds for max x∈supp u−û − (x) and |supp u − | can be obtained as follows: Thus, we can obtain the desired lower bound according to (32).

Upper bound of |D(m)|
Recall that D(m) 6.3 Lower bound of û + L p (D(m)) for p ∈ (1, ∞) It should be noted that û + L p (D(m)) needs to be estimated from below, whereas upper bounds are required for the other constants. More precise estimation is required to satisfy (19) compared with û − L p . Therefore, we use the following estimation: where i∈Λm K i ⊂ D(m). By applying min x∈Kiû + (x) ≥ max{0, m i } to each i ∈ Λ m , estimation of û + L p (D(m)) from below is completed.

Verification theory
In this section, we prepare verification theory to prove the existence of solutions u of (13) satisfying (2) Moreover, suppose that there exists some β > 0 satisfying where D = B(û, 2α + δ) is an open ball depending on the above value α > 0 for small δ > 0. If Furthermore, F ϕ is invertible for every ϕ ∈ B(û, ρ), and the solution u is unique in B(û, 2α).
We set α and β to upper bounds of respectively, then applying Theorem 7.1 to prove the local existence of solutions. Here, L is a positive number satisfying We estimate the inverse norm F −1 u L(H −1 ,H 1 0 ) using the method described in [38,20] in a finite-dimensional subspace V M ⊂ H 1 0 (Ω) specified later.
We define a finite-dimensional subspace V M (⊂ H 1 0 (Ω)) as the tensor prod- We use [13, Theorem 2.3] to obtain an explicit interpolation-error constant then applying [38,20] where P M is the orthogonal projection from H 1 0 (Ω) to V M defined in (42), and u g ∈ H 1 0 (Ω) is a unique solution of the weak formulation of the Poisson equation given g ∈ L 2 (Ω). The upper bound on F(û) H −1 is evaluated using the Raviart-Thomas mixed finite element method (see, for example, [31]). We estimate the inverse norm F −1 u L(H −1 ,H 1 0 ) using (44) and the method described in [38,20]. Remark 7.2 Solutions of (13) on the L-shaped domain may lose H 2 -regularity due to the re-entrant corner at (x, y) = (0.5, 0.5). Therefore, approximate solutions constructed with only a finite element basis may not lead to sufficiently small residuals. We can obtain smaller residuals by constructing approximates solutions with the sum of finite element basis functions and a singular function of the form r 2 3 sin 2 3 θ , where (r, θ) is a polar coordinate system centered at the re-entrant corner. In [14], a priori error estimates for such a singular function are discussed. We also quote [25,Example 7.7], an example of application to nonlinear elliptic problems. In this paper, we do not use the above singular function, but instead make the size of meshes around (x, y) = (0.5, 0.5) smaller than others to reduce residuals (again, see, Fig. 1). Once an H 1 0 -error estimation of a solution is obtained, even if it is relatively rough, our method can be effective in proving the positivity of the solution.

Numerical Experiments
In this section, we present numerical experiments in which the positivity of solutions of (13) satisfying (2) are verified via the proposed method. All computations were implemented on a computer with 2.90 GHz Intel Xeon Platinum 8380H CPUs × 4, 3 TB RAM, and CentOS 8.2 using MATLAB 2019a with GCC Version 8.3.1. All rounding errors were strictly estimated using the toolboxes INTLAB version 11 [30] and kv library version 0.4.49 [12]. In the tables in this section, we use the following notation: have been studied from various points of view (again, see [10,18,4,11,6]). This equation is covered by Theorem 6.5 (see the first row in Table 4). The Lipschitz constant L satisfying (39) can be estimated as  Table 5 shows the verification result and confirms the positivity of the enclosed solution because C 1 ≤ C 2 .

Allen-Cahn equation and Nagumo equation
In the next example, we consider the stationary problem of the Allen-Cahn equation with λ > 0. This is regarded as a special case of the stationary problem of the Nagumo equation with λ > 0 and 0 < a < 1. By applying u = (v + 1)/2 to (48) with a = 0.5 and adjusting the value of λ, these equations become identical. The Lipschitz constants L satisfying (39) for (47) and (48) were respectively estimated as L ≤ 6λC 3 4 ( û L 4 + C 4 r) , r = 2α + δ for small δ > 0, where we set r to be the next floating-point number after 2α. It should be again noted that, in the previous paper [33], another approach was used for (47) to confirm positivity, requiring the confirmation of the positivity of the minimal eigenvalue of a certain linearized operator around approximationû. Theorem 6.5 can be uniformly applied to the nonlinearities of the form (4) (again, see Table 4). For the square domain Ω = (0, 1) 2 , we constructed approximate solutionŝ u using a Legendre polynomial basis. For (47) (λ = 100, 400, and 1600), we obtained Fig. 3 and the verification results in Table 6. For (48) (λ = 400, a = 1/4), we found multiple solutions displayed in Fig. 4 and the verification results in Table 7.

Elliptic equation with multiple positive solutions
For the last example, we consider the elliptic boundary value problem given λ, A, B > 0. This problem has two positive solutions when λ * < λ < λ 1 (Ω) for a certain λ * > 0 (see [19]). The Lipschitz constant L for this problem Lower solution Upper solution Fig. 4 Approximations of multiple solutions of (48) on Ω = (0, 1) 2 for λ = 400 and a = 1/4. can be estimated as where we again set r to be the next floating-point number of 2α.
To prove the positivity of solutions of (49), the previous method [33] requires an L ∞ -error estimation (3) and therefore is applicable to (49) only in the special cases where the solution has H 2 -regularity and we can obtain an explicit bound for the embedding H 2 (Ω) → L ∞ (Ω) successfully. However, the proposed method is well applicable even to (49) without assuming H 2regularity; see again the last case in Table 4.
For the square domain Ω = (0, 1) 2 , we constructed approximationsû of multiple positive solutions of (49) (λ = 10, A = 5, B = 1) using a Legendre polynomial basis, obtaining the results in Fig. 6 and Table 9. For the L-shaped Table 7 Verification results for (48) on Ω = (0, 1) 2 for λ = 400, a = 1/4. The values in rows from F −1 u to ρ and in row C 1 represent strict upper bounds in decimal form. The values in row λ 1 (supp u − ) represent strict lower bounds in decimal form.

Solution
Lower Upper  domain Ω = (0, 1) 2 \([0, 0.5] × [0.5, 1]), we constructed multiple approximate solutionsû ∈ V M of (49) (λ = 20, A = 5, B = 1) using a piecewise quadratic basis, obtaining Fig. 7 and Table 10. Since C 1 ≤ C 2 holds for all cases, the positivity of the solutions are confirmed. Because solutions on the L-shaped domain have not H 2 -regularity and problem (49) corresponds to the lower left case in Table 1, the previous method [33] cannot be applied to solutions in Fig. 7. Nevertheless, the proposed method succeeded in proving the positivity of both lower and upper solutions.

Conclusion
We proposed a unified a posteriori method for verifying the positivity of solutions u of elliptic boundary value problem (1) while assuming H 1 0 -error esti-  mation (2) given some numerical approximationû and an explicit error bound ρ. By extending one of the approaches developed in [33], we designed a unified method with wide applicability. We described the way to obtain explicit values of several constants that the proposed method requires. We also presented numerical experiments to show the effectiveness of our method for three types of nonlinearities, including those to which the previous approach is not applicable. Table 9 Verification results for (49) on Ω = (0, 1) 2 when λ = 10, A = 5, B = 1. The values in rows from F −1 u to ρ and in row C 1 represent strict upper bounds in decimal form. The values in rows λ 1 (supp u − ) and C 2 represent strict lower bounds in decimal form.

A Estimates of the first positive zeros of the Bessel functions
This section discusses rigorous estimates of the Bessel functions j n,1 required to obtain the values displayed in Table 2. When n is written in the form n = k + 0.5 with an integer k, explicit formulas for j n,1 can be obtained. Particularly, j 0.5,1 = 2/(πx) sin x; therefore, the first zero of this is π (see, for example, [3, Remark 1.2] and [27,Section 10.16]). Moreover, we have j 1.5,1 = 2/(πx) x −1 sin x − cos x , the zeros of which satisfy x = tan x. The first zero of this function is enclosed by the function "allsol" packaged in the kv library [12], which enables us to obtain all zeros of the function in a given compact interval. We set the initial interval as [π, 1.5π], in which the equation has the first zero. When n = 0, 1, we calculated rigorous values of j n,1 using the bisection method with computer assistance. For preparation, we first prove that there is no positive zero of j 0,1 and j 1,1 in the compact interval [0, 1] as follows: When n is an integer, j n,1 is given by When a ∈ [0, 1] ⊂ [0, π/2), cos a is positive, and therefore, so is j 0,1 . Next, we consider j 1,1 , which satisfies j 1,1 (0) = 0. The first derivative of j 1,1 is given by d dx j 1,1 (x) = 1 π π 0 sin t sin(t − x sin t)dt.