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

Abstract

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 370:112647, 2020), 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.

1 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, [5, 23,24,25,26, 28, 29, 31, 33] and the references therein). One such successful application is to the semilinear elliptic boundary value problem

$$\begin{aligned} \left\{ \begin{array}{ll} -\varDelta u=f(u)&{}\text {in}~\varOmega ,\\ u=0&{}\text {on}~\partial \varOmega , \end{array} \right. \end{aligned}$$

(1)

where $\varOmega \subset \mathbb {R}^{N}$ $(N=2,3,\ldots )$ is a given bounded domain, $\varDelta$ is the Laplacian, and $f: \mathbb {R} \rightarrow \mathbb {R}$ is a given nonlinear map. Further regularity assumptions for $\varOmega$ and f will be shown later for our setting.

Positive solutions of (1) have attracted significant attention [4, 6, 10, 11, 18, 19]. For example, positive solutions of problem (1) with $f(t) = \lambda t + t|t|^{p-1}$, $\lambda \in [0,\lambda _1(\varOmega ))$, $p\in (1,p^*)$ have been investigated from various points of view such as uniqueness, multiplicity, nondegeneracy, and symmetry, [4, 6, 10, 11, 18], where $p^*=\infty$ when $N=2$ and $p^*=(N+2)/(N-2)$ when $N\ge 3$; $\lambda _1(\varOmega )$ is the first eigenvalue of $-\varDelta$ with the homogeneous Dirichlet boundary value condition in the weak sense. Another important nonlinearity is $f(t)=\lambda (t-t^3)$, $\lambda >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 $\lambda \ge \lambda _1(\varOmega )$. However, when $\lambda < \lambda _1(\varOmega )$, no positive solution is admitted; this can be confirmed by multiplying $-\varDelta u = \lambda (u-u^3)$ with the first eigenfunction of $-\varDelta$ and integrating both sides. The Allen–Cahn equation is a special case of the Nagumo equation [22] with the nonlinearity $f(t)=\lambda t(1-t)(t-a)$, $\lambda >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)=\lambda (t+{\mathrm{A}}t^2 - {\mathrm{B}}t^3)$ with ${\mathrm{A}},{\mathrm{B}}>0$. The bifurcation of problem (1) with this nonlinearity was analyzed in [19]. This problem has two positive solutions when $\lambda ^*< \lambda < \lambda _1(\varOmega )$ for some $\lambda ^*>0$. Despite these results, quantitative information about the positive solutions, such as their shape, has not been clarified analytically. Throughout this paper, $H^k(\varOmega )$ denotes the kth order $L^2$ Sobolev space. We define $H^1_0(\varOmega ):= \{ u \in H^1(\varOmega ) : u = 0~\text {on} ~ \partial \varOmega \}$, with the inner product $(u, v)_{H^1_0}:=(\nabla u, \nabla v)_{L^2}$ and norm $\Vert u \Vert _{H^1_0} := \sqrt{(u, u)_{\smash {H^1_0}}}$.

Numerical verification methods enable us to obtain an explicit ball containing exact solutions of (1). More precisely, for a numerical approximation $\hat{u} \in H^1_0(\varOmega )$ that satisfies the assumption required by such methods, they prove the existence of an exact solution $u \in H^1_0(\varOmega )$ of (1) that satisfies

$$\begin{aligned} \left\| u-\hat{u}\right\| _{H_{0}^{1}} \le \rho \end{aligned}$$

(2)

with an explicit error bound $\rho >0$. Under an appropriate condition, we can obtain an $L^{\infty }$-estimation

$$\begin{aligned} \left\| u-\hat{u}\right\| _{L^{\infty }} \le \sigma \end{aligned}$$

(3)

with bound $\sigma >0$. For instance, when $u,\hat{u} \in H^2(\varOmega )$, we can evaluate the $L^{\infty }$-bound $\sigma >0$ by considering the embedding $H^2(\varOmega ) \hookrightarrow L^{\infty }(\varOmega )$ [28, Theorem 1, 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 ($\rho$ or $\sigma$) 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 $\partial \varOmega$.

Therefore, we developed methods of verifying the positivity of solutions of (1) in previous studies [33, 35,36,37] 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,36,37], 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, 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

$$\begin{aligned} f(t) = \lambda t + \sum _{i=2}^{n(<p^*)} a_{i}t|t|^{i-1},~\lambda , a_i \in \mathbb {R},~a_i \ne 0 ~\text {for some}~ i. \end{aligned}$$

(4)

Theorem 2.1 in [33] can be applied to cases in which $\lambda <\lambda _1(\varOmega )$. 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^{\infty }$-estimation (3) but only requires an $H^1_0$-error estimation (2). When $a_i \le 0$ for all i and $\lambda \ge \lambda _1(\varOmega )$, 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^{\infty }$-estimation but requires an explicit evaluation of the minimal eigenvalue of a certain linearized operator around approximation $\hat{u}$. Actually, this eigenvalue evaluation itself is not trivial to obtain. In [33], the eigenvalue was estimated using the existing method [20, 33] 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) = \lambda t + t|t|^{p-1}$ and $f(t) = \lambda (t-t^3)$. However, a problem still remains in the sense that we need the $L^{\infty }$-error estimation (3) to prove positivity when $a_i a_j < 0$ for some i, j and $\lambda \ge \lambda _1(\varOmega )$. For example, the nonlinearity $f(t)=\lambda (t+{\mathrm{A}}t^2 - {\mathrm{B}}t^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, Corollary 1] requires $u,\hat{u} \in H^2(\varOmega )$, evaluating the bound for the embedding $H^2(\varOmega ) \hookrightarrow L^{\infty }(\varOmega )$, to obtain $\sigma$ from $\rho$. The $H^2$-regularity of u may fall when $\varOmega$ is a nonconvex polygonal domain.

Table 1 Error estimates required in [33] when f is the subcritical polynomial (4)

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 1 for arbitrary bounded domains $\varOmega$. 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 $\lambda _{1}({\mathrm{supp}}\,{u_{-}})$ explicitly because the inequality $\lambda _1({\mathrm{supp}}\,{u_{-}}) \ge \lambda$ has to be confirmed for the success of our positivity proof (again, see Lemma 6.1). Here, ${\mathrm{supp}}\,{u_{-}}:= \overline{\{x \in \varOmega : u_{-}(x) \ne 0\}}$ is the support of $u_{-}$, and $\lambda _1({\mathrm{supp}}\,{u_{-}})$ is understood as the first eigenvalue on the interior of ${\mathrm{supp}}\,{u_{-}}$. When the interior of ${\mathrm{supp}}\,{u_{-}}$ is empty, we interpret $\lambda _{1}({\mathrm{supp}}\,{u_{-}})=\infty$, and all real numbers are lower bounds of $\lambda _{1}({\mathrm{supp}}\,{u_{-}})$. The difficulty is that we cannot identify the location and shape of ${\mathrm{supp}}\,{u_{-}}$ from (2) even when $\hat{u}$ is nonnegative.

If a polygon or polyhedron S enclosing ${\mathrm{supp}}\,{u_{-}}$ is obtained concretely, we can apply the Liu–Oishi method [20, 21] based on finite element methods to obtain a lower bound for $\lambda _{1}(S)$. Then, the inequality $\lambda _{1}({\mathrm{supp}}\,{u_{-}}) \ge \lambda _{1}(S)$ gives the desired lower bound of $\lambda _{1}({\mathrm{supp}}\,{u_{-}})$. Such a supremum set S over ${\mathrm{supp}}\,{u_{-}}$ can be obtained when we have an $L^{\infty }$-estimation (3). By setting $\varOmega _{+}:=\{x\in \varOmega : \hat{u}-\sigma \ge 0\}$ where $u \ge 0$ therein, we can construct such a domain S as a supremum set of $\varOmega \backslash \varOmega _{+}$. 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 $\lambda _1({\mathrm{supp}}\,{u_{-}})$ more accurately (see, for example, [25, Theorem 10.31]).

To estimate the lower bound of $\lambda _{1}({\mathrm{supp}}\,{u_{-}})$, we rely on the following argument:

Fact 1.1

For a bounded domain $\varOmega \subset \mathbb {R}^{N}$ $(N=2,3,\ldots )$, there exists a constant $A_{k,N}$ independent of $\varOmega$ such that

$$\begin{aligned} \lambda _{k}(\varOmega ) \ge A_{k,N} \left( \frac{1}{|\varOmega |}\right) ^{\frac{2}{N}}, \end{aligned}$$

(5)

where $\lambda _{k}(\varOmega )$ denotes the k-th eigenvalue of $-\varDelta$ with the homogeneous Dirichlet boundary condition.

Many articles have investigated this type of inequality in several forms. Among them, we mainly use the Rayleigh–Faber–Krahn inequality, which ensures Fact 1.1 for $k=1$ [9, 15, 16]. This inequality states that if $|\varOmega | = |\varOmega ^*|$ ($\varOmega ^*$ is a ball in $\mathbb {R}^{N}$), then $\lambda _{1}(\varOmega ) \ge \lambda _{1}(\varOmega ^*)$, where the equality holds if and only if $\varOmega =\varOmega ^*$. We also refer to [17] for an easy-to-estimate formula for $A_{k,N}$ for all $k \ge 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 $\lambda _1({\mathrm{supp}}\,{u_{-}})$ using the inequality (5), we focus on estimating upper bounds of $|{\mathrm{supp}}\,{u_{-}}|$ while assuming only the $H^1_0$-error estimation (2) without knowing the specific shape and location of ${\mathrm{supp}}\,{u_{-}}$. Suppose that (2) is proved for a positive approximation $\hat{u}$ with sufficient accuracy. Then, the upper bound for $|{\mathrm{supp}}\,{u_{-}}|$ can be estimated very small using Lemma 3.2 provided later, and therefore, $\lambda _1({\mathrm{supp}}\,{u_{-}}) \ge \lambda$ can be confirmed using (5) for a moderately large $\lambda$. The established estimation for $|{\mathrm{supp}}\,{u_{-}}|$ is used to evaluate not only $\lambda _{1}({\mathrm{supp}}\,{u_{-}})$ but also some Sobolev embedding constants on $|{\mathrm{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 Sect. 3, we evaluate an upper bound for the volume $|{\mathrm{supp}}\,{u_{-}}|$ assuming the $H^1_0$-estimation (2) for a continuous or piecewise continuous approximation $\hat{u} \in H^1_0(\varOmega )$. Subsequently, we use the bound for $|{\mathrm{supp}}\,{u_{-}}|$ to evaluate lower bounds for $\lambda _1({\mathrm{supp}}\,{u_{-}})$ in Sect. 4. In Sect. 5, required Sobolev embedding constants on bounded domains are evaluated. In Sect. 6, we extend the previous formula [33, Theorem 2.1] and combine it with the estimates derived from Sects. 4 and 5, thereby designing a unified method for proving positivity. Finally, Sect. 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^{\infty }$-error estimation.

2 Preliminaries

We begin by introducing required notation. We denote by $H^{-1}$ the topological dual of $H^1_0(\varOmega )$. 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 $\Vert \cdot \Vert _{L^p} = \Vert \cdot \Vert _{L^p(\varOmega )}$ if no confusion arises. For two Banach spaces X and Y, the set of bounded linear operators from X to Y is denoted by ${\mathcal L}(X, Y)$ with the usual supremum norm $\Vert T \Vert _{{\mathcal L}(X, Y)} := \sup \{ \Vert T u \Vert _{Y} / \Vert u \Vert _{X} : {u \in X \setminus \{0\}}\}$ for $T \in {\mathcal L}(X, Y)$. The norm bound for the embedding $H^1_0(\varOmega )\hookrightarrow L^{p+1}\left( \varOmega \right)$ is denoted by $C_{p+1}(=C_{p+1}(\varOmega ))$; that is, $C_{p+1}$ is a positive number that satisfies

$$\begin{aligned} \left\| u\right\| _{L^{p+1}(\varOmega )}\le C_{p+1}\left\| u\right\| _{H^1_0(\varOmega )}~~~\mathrm{for~all}~u\in H^1_0(\varOmega ), \end{aligned}$$

(6)

where $p\in [1,\infty )$ when $N=2$ and $p\in [1,p^*]$ when $N\ge 3$. If no confusion arises, we use the notation $C_{p+1}$ to represent the embedding constant on the entire domain $\varOmega$, whereas, in some parts of this paper, we need to consider an embedding constant on some subdomain $\varOmega '\subset \varOmega$. This is denoted by $C_{p+1}(\varOmega ')$ to avoid confusion. Moreover, $\lambda _1(\varOmega )$ denotes the first eigenvalue of $-\varDelta$ imposed on the homogeneous Dirichlet boundary condition. This is characterized by

$$\begin{aligned} \lambda _{1}(\varOmega ) = \inf _{v\in H^1_0(\varOmega )^\backslash {\{0\}}} \frac{\Vert v\Vert _{H^1_0(\varOmega )}^2}{\Vert v\Vert _{L^2(\varOmega )}^2}. \end{aligned}$$

(7)

Throughout this paper, we assume that f is a $C^1$ function that satisfies

$$\begin{aligned}&|f(t)| \le a_0 |t|^p + b_0 ~~\text {for all}~ t \in \mathbb {R}, \end{aligned}$$

(8)

$$\begin{aligned}&|f'(t)| \le a_1 |t|^{p-1} + b_1 ~~\text {for all}~t \in \mathbb {R} \end{aligned}$$

(9)

for some $a_0,a_1,b_0,b_1\ge 0$ and $p\in [1,p^*)$. We define the operator F as

$$\begin{aligned} F : \left\{ \begin{array}{ccc}{u(\cdot )} &{} {\mapsto } &{} {f(u(\cdot ))}, \\ {H^1_0(\varOmega )} &{} {\rightarrow } &{} {H^{-1}}.\end{array}\right. \end{aligned}$$

Moreover, we define another operator $\mathcal {F} : H^1_0(\varOmega )\rightarrow H^{-1}$ as $\mathcal {F}(u) :=-\varDelta u-F(u)$, which is characterized by

$$\begin{aligned} \left<\mathcal F(u),v\right> = \left( \nabla u,\nabla v\right) _{L^2} - \left<F(u),v\right> ~~\text {for all}~ u,v \in H^1_0(\varOmega ),\end{aligned}$$

(10)

where $\left<F(u),v\right> = \int _{\varOmega } f(u(x)) v(x) dx$. The Fréchet derivatives of F and $\mathcal {F}$ at $\varphi \in H^1_0(\varOmega )$, denoted by ${F'_{\varphi }}$ and ${\mathcal F'_{\varphi }}$, respectively, are given by

$$\begin{aligned}&\langle F'_{\varphi }u,v\rangle = \int _{\varOmega } f'(\varphi (x))u(x) v(x) dx ~~\text {for all}~ u,v \in H^1_0(\varOmega ), \end{aligned}$$

(11)

$$\begin{aligned}&\langle \mathcal F'_{\varphi }u,v \rangle = \left( \nabla u,\nabla v\right) _{L^2} - \langle F'_{\varphi }u,v \rangle ~~\text {for all}~ u,v \in H^1_0(\varOmega ). \end{aligned}$$

(12)

Under the notation and assumptions, we look for solutions $u \in H^1_0(\varOmega )$ of

$$\begin{aligned} \mathcal {F}(u)=0, \end{aligned}$$

(13)

which corresponds to the weak form of (1). We assume that some verification method succeeds in proving the existence of a solution $u \in H^1_0(\varOmega )$ of (13) satisfying inequality (2) given $\hat{u} \in H^1_0(\varOmega )$ and $\rho >0$. Although the regularity assumption for $\hat{u}$ (to be in $H^1_0(\varOmega )$) is sufficient to obtain the error bound (2) in theory, we further assume that $\hat{u}$ 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\}$.

3 Evaluation of the volume of ${\mathrm{supp}}\,{u_{-}}$

To estimate an upper bound for $|{\mathrm{supp}}\,{u_{-}}|$ from the information of the inclusion (2), we define $D(v):=\{ x \in \varOmega : \hat{u}(x)\le v(x)\}$ and consider the maximization problem

$$\begin{aligned} \mathop {\mathrm {maximize}}_{\left\| v\right\| _{L^{q}(\varOmega )}= c} |D(v)| \end{aligned}$$

(14)

for fixed $q \in (1,\infty )$ and $c>0$. This maximization takes place over the set of all functions $v \in L^{q}(\varOmega )$ satisfying $\left\| v\right\| _{L^{q}(\varOmega )}= c$. When $\left\| \hat{u}_{+}\right\| _{L^{q}(\varOmega )} \le c$, the maximal value of the problem (14) is $|\varOmega |$. Therefore, we consider the case where $\left\| \hat{u}_{+}\right\| _{L^{q}(\varOmega )} > c$. In the following, we denote $D(l) := \{ x \in \varOmega : \hat{u}(x)\le l\}$ and $\mathring{D}(l) := \{ x \in \varOmega : \hat{u}(x)< l\}$ for $l \in \mathbb {R}$.

Lemma 3.1

Let $q \in (1,\infty )$ and $c>0$ be fixed. Suppose that $\left\| \hat{u}_{+}\right\| _{L^{q}(\varOmega )} > c$. Then, we have

$$\begin{aligned} \mathop {\mathrm {arg~max}}_{\left\| v\right\| _{L^{q}(\varOmega )}= c} |D(v)| = \left\{ \begin{array}{ll} \hat{u}_+(x),&{}{}x\in D,\\ 0,&{}{} \text{ otherwise }, \end{array} \right. \end{aligned}$$

(15)

where D is a set that satisfies $\left\| \hat{u}_{+}\right\| _{L^{q}(D)}= c$ and

$$\begin{aligned} \mathring{D}(l) \subseteq D \subseteq D(l) \end{aligned}$$

(16)

for some $l \in \mathbb {R}$. The maximal value of the problem (14) is |D|.

Proof

Let us denote by $v_c$ an arbitrary function in $\{v \in L^{q}(\varOmega ): \left\| v\right\| _{L^{q}(\varOmega )}= c\}$. Because $D(v_c)\subseteq D(|v_c|)$ and the equality holds when $v_c \ge 0$ in $\varOmega$, the volume of $D(v_c)$ is maximized for a nonnegative $v_c$.

Let $v_c$ be nonnegative. If $\hat{u}-v_c$ is strictly negative in some part $\varOmega ' \subset {\mathrm{supp}}\,v_c$ satisfying $|\varOmega '| \ne 0$ and vanishes in $({\mathrm{supp}}\,v_c) \backslash \varOmega '$, another $v'_c \in \{v \in L^{q}(\varOmega ): \left\| v\right\| _{L^{q}(\varOmega )}= c\}$ with the same $c>0$ can be constructed to obtain larger $D(v'_c)$ as follows. Since $\left\| \hat{u}_{+}\right\| _{L^{q}(\varOmega )} > c = \left\| v_{c}\right\| _{L^{q}({\mathrm{supp}}\,v_c)}$, we have

$$\begin{aligned} \left\| \hat{u}_{+}\right\| _{L^{q}(\varOmega \backslash {\mathrm{supp}}\,v_c)}^{q}&> \left\| v_{c}\right\| _{L^{q}({\mathrm{supp}}\,v_c)}^{q}-\left\| \hat{u}_{+}\right\| _{L^{q}({\mathrm{supp}}\,v_c)}^{q}\\&= \displaystyle \int _{{\mathrm{supp}}\,v_c}\left( v_{c}(x)^{q}-\hat{u}_{+}(x)^{q}\right) dx\\&=\left\| v_c\right\| ^q_{L^{q}(\varOmega ')} - \left\| \hat{u}_{+}\right\| ^q_{L^{q}(\varOmega ')}. \end{aligned}$$

Therefore, there exists $\varOmega '' \subset \varOmega \backslash ({\mathrm{supp}}\,v_c)$ that satisfies $\left\| \hat{u}_{+}\right\| ^q_{L^{q}(\varOmega '')}=\left\| v_c\right\| ^q_{L^{q}(\varOmega ')} - \left\| \hat{u}_{+}\right\| ^q_{L^{q}(\varOmega ')}$. Defining $v'_c$ as

$$\begin{aligned} v'_c(x)= \left\{ \begin{array}{ll} \hat{u}_+(x),&{}x\in \varOmega ' \cup \varOmega '',\\ v_c(x),&{}\text {otherwise}, \end{array} \right. \end{aligned}$$

we have

$$\begin{aligned} \left\| v'_c\right\| _{L^{q}(\varOmega )}^q&= \left\| \hat{u}_{+}\right\| _{L^{q}(\varOmega ')}^q + \left\| \hat{u}_{+}\right\| _{L^{q}(\varOmega '')}^q + \left\| v_c\right\| _{L^{q}(\varOmega \backslash (\varOmega ' \cup \varOmega ''))}^q\\&= \left\| \hat{u}_{+}\right\| _{L^{q}(\varOmega ')}^q -\left\| \hat{u}_{+}\right\| ^q_{L^{q}(\varOmega ')} + \left\| v_c\right\| ^q_{L^{q}(\varOmega ')} + \left\| v_c\right\| _{L^{q}(\varOmega \backslash (\varOmega ' \cup \varOmega ''))}^q\\&= \left\| v_c\right\| _{L^{q}(\varOmega \backslash \varOmega '')}^q = \left\| v_c\right\| _{L^{q}(\varOmega )}^q \end{aligned}$$

and $D(v_c)\subset D(v'_c)$ in the strict sense because $|\varOmega ''| \ne 0$.

Therefore, when $|D(v_c)|$ is maximized, $\hat{u}-v_c$ vanishes in ${\mathrm{supp}}\,v_c$; that is,

$$\begin{aligned} v_c(x)=\hat{u}_{+}(x),~~x \in {\mathrm{supp}}\,v_c. \end{aligned}$$

(17)

Finally, we consider a subset $D\,(={\mathrm{supp}}\,v_c) \subset \varOmega$ with the largest volume satisfying $\left\| \hat{u}_{+}\right\| _{L^{q}(D)}= c$. In the following, we prove that the volume of such D is maximized when (16) holds for some $l \in \mathbb {R}$.

Since $\hat{u}$ is continuous or piecewise continuous on $\varOmega$, there exist D and $l \in \mathbb {R}$ that satisfy $\left\| \hat{u}_{+}\right\| _{L^{q}(D)}= c$ and (16). Suppose that there exists a different set $D' \subset \varOmega$ that satisfies $\left\| \hat{u}_{+}\right\| _{L^{q}(D')}= c$ so that

$$\begin{aligned} \left\| \hat{u}_{+}\right\| _{L^{q}(D'\backslash (D\cap D'))} = \left\| \hat{u}_{+}\right\| _{L^{q}(D\backslash (D\cap D'))}. \end{aligned}$$

(18)

Then, we have that $\hat{u}_{+}(x)\le l$ for all $x \in D\backslash (D\cap D') \subset D$ and $\hat{u}_{+}(x)\ge l$ for all $x \in D'\backslash (D\cap D') \subset \varOmega \backslash D$ because (16). It follows from (18) that $|D\backslash (D\cap D')|\ge |D'\backslash (D\cap D')|$, and therefore, $|D|\ge |D'|$. Thus, the assertion of this lemma is proved. $\square$

The following lemma provides the desired upper bound for $|{\mathrm{supp}}\,{u_{-}}|$ on the basis of Lemma 3.1.

Lemma 3.2

Given $q \in [2,p^*+1)$ and $m>0$, if we have

$$\begin{aligned} \left\| \hat{u}_{+}\right\| _{L^{q}(D(m))} > C_q \rho , \end{aligned}$$

(19)

then

$$\begin{aligned} |{\mathrm{supp}}\,{u_{-}}|\le |D(m)|. \end{aligned}$$

(20)

Proof

The enclosed solution u can be expressed by $u=\hat{u}-\omega$, $\Vert \omega \Vert _{H^1_0(\varOmega )}\le \rho$. The embedding $H^1_0(\varOmega )\hookrightarrow L^{q}(\varOmega )$ confirms that $\Vert \omega \Vert _{L^{q}(\varOmega )}\le C_q \rho$. We denote $c:=\Vert \omega \Vert _{L^{q}(\varOmega )}$, and then, have

$$\begin{aligned} |{\mathrm{supp}}\,{u_{-}}| \le \max _{\left\| v\right\| _{L^{q}(\varOmega )}= c} |D(v)|. \end{aligned}$$

Lemma 3.1 ensures that this maximal value is realized when

$$\begin{aligned} v(x) = \left\{ \begin{array}{ll} \hat{u}_+(x),&{}x\in D,\\ 0,&{}\text { otherwise}, \end{array} \right. \end{aligned}$$

where D is a set that satisfies $\left\| \hat{u}_{+}\right\| _{L^{q}(D)}= c$ and (16) for some $l \in \mathbb {R}$. The maximal value is |D|, which is not greater than |D(l)| due to (16). Inequality (19) ensures that

$$\begin{aligned} \left\| \hat{u}_{+}\right\| _{L^{q}(D)} = c \le C_q \rho < \left\| \hat{u}_{+}\right\| _{L^{q}(D(m))}. \end{aligned}$$

(21)

Suppose that $D(m) \subset D(l)$ in the strict sense. Then, we have $m<l$, and thus, $D(m) \subseteq \mathring{D}(l) \subseteq D$ due to (16). This contradicts (21). Therefore, we have $D(l) \subseteq D(m)$ and conclude (20). $\square$

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, $\rho$, and $\hat{u}$ be fixed. Then, $\left\| \hat{u}_{+}\right\| _{L^{q}(D(m))}$ monotonically decreases as m decreases, and $\left\| \hat{u}_{+}\right\| _{L^{q}(D(m))} \downarrow 0$ as $m \downarrow 0$. Therefore, although smaller m gives a good upper bound of $|{\mathrm{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 Sect. 8.

4 Lower bound for the minimal eigenvalue

The purpose of this section is to estimate a lower bound of the minimal eigenvalue $\lambda _1({\mathrm{supp}}\,{u_{-}})$ while we assume an $H^1_0$-estimation (2) only; therefore, ${\mathrm{supp}}\,{u_{-}}$ cannot be identified explicitly. To this end, we use the following Rayleigh–Faber–Krahn constant.

Theorem 4.1

([9, 15, 16]) Inequality (5) with $k=1$ holds for

$$\begin{aligned} A_{1,N} = B_N^{\frac{2}{N}} j_{\frac{N}{2}-1,1}^2 \end{aligned}$$

(22)

where $B_N = \pi ^{N/2}/\Gamma (N/2+1)$ denotes the volume of the unit N-ball with the usual gamma function $\Gamma$, and $j_{\frac{N}{2}-1,1}$ is the first positive zero of the Bessel function of order $\frac{N}{2}-1$. The equality in (5) is attained if and only if $\varOmega$ is a ball in $\mathbb {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_{\frac{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_{\frac{N}{2}-1,1}$ is built therein. Accordingly, we present a rigorous algorithm for calculating $j_{\frac{N}{2}-1,1}$ in Appendix A based on the bisection method.

Table 2 Strict enclosures of $B_N$, $j_{\frac{N}{2}-1,1}$, and $A_{1,N}$ for $N=2,3,4,5$

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 ${\mathrm{supp}}\,{u_{-}}$.

Corollary 4.2

If (19) holds given $q \in [2,p^*)$ and $m>0$, then we have

$$\begin{aligned}&\lambda _{1}({\mathrm{supp}}\,{u_{-}}) \ge 18.1684145355 |D(m)|^{-1},&N=2,\\&\lambda _{1}({\mathrm{supp}}\,{u_{-}}) \ge 25.6463452794 |D(m)|^{-\frac{2}{3}},&N=3. \end{aligned}$$

Remark 4.3

Instead of Theorem 4.1, one can use the evaluation provided in [17]: For a bounded domain $\varOmega$, we have

$$\begin{aligned} \lambda _{k}(\varOmega ) \ge \frac{4\pi ^2 N}{N+2} \left( \frac{k}{B_N|\varOmega |}\right) ^{\frac{2}{N}}. \end{aligned}$$

(23)

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$.

5 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].

Theorem 5.1

([36, Corollary A.2]) Let $\varOmega \subset \mathbb {R}^{N}(N \ge 2)$ be a bounded domain, the measure of which is denoted by $| \varOmega |$. Let $p \in (N/(N-1),2N/(N-2) ]$ if $N \ge 3$, $p \in (2, \infty )$ if $N=2$. Then, (6) holds for

$$\begin{aligned} C_p( \varOmega ) = | \varOmega |^{\frac{1}{N}+\frac{1}{p}-\frac{1}{2}}T_{p,N}. \end{aligned}$$

Here, $T_{p,N}$ is defined by

$$\begin{aligned} T_{p,N}=\pi ^{-\frac{1}{2}} N^{-\frac{1}{q}}\left( \frac{q-1}{N-q}\right) ^{1-\frac{1}{q}}\left\{ \frac{\Gamma \left( 1+\frac{N}{2}\right) \Gamma (N)}{\Gamma \left( \frac{N}{q}\right) \Gamma \left( 1+N-\frac{N}{q}\right) }\right\} ^{\frac{1}{N}}, \end{aligned}$$

(24)

where $\Gamma$ is the gamma function and $q=Np/(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(\varOmega )$. Moreover, in [36], one can find very sharp estimations of the best values of $C_p$ for $p=3,4,5,6,7$ on $\varOmega =(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.

Table 3 Upper bounds for embedding constants

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:

$$\begin{aligned} \Vert u\Vert _{L^{2}(\varOmega )} \le \frac{1}{\sqrt{\lambda _{1}(\varOmega )}}\Vert u\Vert _{H^1_0(\varOmega )}. \end{aligned}$$

(25)

For example, when $\varOmega =(0,1)^N$ $(N=1,2,3,\ldots )$, we have the exact eigenvalue $\lambda _{1}(\varOmega )=N\pi ^{2}$. Even otherwise, lower bounds of $\lambda _{1}(\varOmega )$ (and therefore upper bounds for the embedding constant) can be estimated for bounded domain $\varOmega$ using the formulas in Sect. 4.

6 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. Let f satisfy

$$\begin{aligned} -f(-t)\le \lambda t + \displaystyle \sum _{i=1}^{n}a_{i}t^{p_{i}} ~~\text { for all}~t\ge 0 \end{aligned}$$

(26)

for some $\lambda \in \mathbb {R}$, nonnegative coefficients $a_1,a_2,\ldots ,a_n$, and subcritical exponents $p_1,p_2,\ldots ,p_n \in (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, Theorem 2.1]. In the following, $\lambda _1({\mathrm{supp}}\,{u_{-}})$ and $C_{p+1}({\mathrm{supp}}\,{u_{-}})$ are interpreted as $\lambda _1(\text {int }\,({\mathrm{supp}}\,{u_{-}}))$ and $C_{p+1}(\text {int }\,({\mathrm{supp}}\,{u_{-}}))$, respectively, where $\text {int }\, ({\mathrm{supp}}\,{u_{-}})$ denotes the interior of ${\mathrm{supp}}\,{u_{-}}$. We denote by $\underline{\lambda }\in \mathbb {R}$ a lower bound of $\lambda _{1}({\mathrm{supp}}\,{u_{-}})$. When the interior of ${\mathrm{supp}}\,{u_{-}}$ is empty, we can set $\underline{\lambda }$ to an arbitrarily large value.

Lemma 6.1

Suppose that a solution $u \in H^1_0(\varOmega )$ of (13) exists in $\overline{B}(\hat{u},\rho )$, given some approximation $\hat{u} \in H^1_0(\varOmega )$ and radius $\rho >0$. If

$$\begin{aligned} \lambda < \underline{\lambda }\end{aligned}$$

(27)

and

$$\begin{aligned} \displaystyle \sum _{i=1}^{n}a_i C_{p_{i}+1}({\mathrm{supp}}\,{u_{-}})^{2}\left( \left\| \hat{u}_{-}\right\| _{L^{p_{i}+1}}+C_{p_{i}+1}\rho \right) ^{p_{i}-1}<1 - \frac{\lambda }{\underline{\lambda }}, \end{aligned}$$

(28)

then u is nonnegative. Note that when ${\mathrm{supp}}\,{u_{-}}$ is disconnected, (28) is understood as the set of inequalities for all connected components of ${\mathrm{supp}}\,{u_{-}}$.

Proof

We extend the proof of [33, Theorem 2.1] to achieve a more precise evaluation. For $p\in (1,p^*)$, we have

$$\begin{aligned} \left\| u_{-}\right\| _{L^{p+1}} \le \left\| \hat{u}_{-}\right\| _{L^{p+1}}+ C_{p+1}\rho ; \end{aligned}$$

(29)

see the proof of [33, Theorem 2.1]. We then prove that the norm of $u_{-}$ vanishes. Because u satisfies

$$\begin{aligned} \left( \nabla u,\nabla v\right) _{L^2}=\left<F(u),v\right>\mathrm{~~for~all~}v\in H^1_0(\varOmega ), \end{aligned}$$

by fixing $v=u_{-}$, we have from (26) that

$$\begin{aligned} \left\| u_{-}\right\| _{H^1_0(\varOmega )}^{2} \le&\displaystyle \int _{\varOmega }\left\{ \lambda \left( u_{-}(x)\right) ^{2}+\sum _{i=1}^{n}a_{i}\left( u_{-}(x)\right) ^{p_{i}+1}\right\} dx \nonumber \\ =&\displaystyle \lambda \left\| u_{-}\right\| _{L^{2}}^{2}+\sum _{i=1}^{n}a_{i}\left\| u_{-}\right\| _{L^{p_{i}+1}}^{p_{i}+1}\nonumber \\ \le&\left\{ \displaystyle \frac{\lambda }{\underline{\lambda }} +\sum _{i=1}^{n}a_{i}C_{p_{i}+1}({\mathrm{supp}}\,{u_{-}})^{2} \left\| u_{-}\right\| _{L^{p_{i}+1}}^{p_{i}-1} \right\} \left\| u_{-}\right\| _{H^1_0(\varOmega )}^{2}. \end{aligned}$$

(30)

Inequalities (28) and (29) lead to

$$\begin{aligned} \displaystyle \frac{\lambda }{\underline{\lambda }} +\sum _{i=1}^{n}a_{i}C_{p_{i}+1}({\mathrm{supp}}\,{u_{-}})^{2} \left\| u_{-}\right\| _{L^{p_{i}+1}}^{p_{i}-1} < 1, \end{aligned}$$

(31)

which ensures $\Vert u_{-}\Vert _{H^1_0(\varOmega )} = 0$. Thus, the nonnegativity of u is proved. $\square$

Remark 6.2

The maximum principle ensures the positivity of u (i.e., $u(x)>0$ inside $\varOmega$) from its nonnegativity for a wide class of nonlinearities f. See, for example, [8] for a generalized maximum principle applicable for weak solutions.

Remark 6.3

The formula $\left\| \hat{u}_{-}\right\| _{L^{p_{i}+1}}+C_{p_{i}+1}\rho$ in parentheses in (28) indicates that the nonnegativity of u can be confirmed using Lemma 6.1 if the norm of $\hat{u}_{-}$ is sufficiently small and verification succeeds with sufficient accuracy.

Remark 6.4

Lemma 6.1 can be indeed regarded as a generalization of [33, Theorem 2.1 and Corollary A.1] because Lemma 6.1 is formally obtained from [33, Theorem 2.1] via the replacements $\lambda _1(\varOmega ) \rightarrow \underline{\lambda }$ and the left-side $C_{p_{i}+1}(\varOmega ) \rightarrow C_{p_{i}+1}({\mathrm{supp}}\,{u_{-}})$. Actually, the replacement for $\lambda _1(\varOmega )$ was already discussed in [33, Corollary A.1]. However, the embedding constant $C_{p_{i}+1}(\varOmega )$ was not replaced. Because the formula $\left\| \hat{u}_{-}\right\| _{L^{p_{i}+1}}+C_{p_{i}+1}\rho$ in parentheses in (28) can be very small as mentioned in Remark 6.3, replacing $C_{p_{i}+1}({\mathrm{supp}}\,{u_{-}})$ with its rough bound $C_{p_{i}+1}(=C_{p_{i}+1}(\varOmega ))$ is rarely problematic. However, in several cases, roughly estimating the lower bound for $\lambda _1({\mathrm{supp}}\,{u_{-}})$ via $\lambda _1({\mathrm{supp}}\,{u_{-}}) \ge \lambda _1(\varOmega )$ cannot satisfy (28). For example, when $f(t)=\lambda (t - t^3)$, $\lambda$ should satisfy $\lambda \ge \lambda _1(\varOmega )$ to admit a positive solution.

By applying the estimations in Sects. 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 \in H^1_0(\varOmega )$ of (13) exists in $\overline{B}(\hat{u},\rho )$, given some approximation $\hat{u} \in H^1_0(\varOmega )$ and radius $\rho >0$. Moreover, we assume (19), given $q \in [2,p^*+1)$ and $m>0$. Then, we define $\mathcal C_1=\mathcal C_1(N,f,\hat{u},\rho ,m)$ and $\mathcal C_2=\mathcal C_2(N,f,\hat{u},m)$ as

$$\begin{aligned} \mathcal C_1&:=\displaystyle \sum _{i=1}^{n}a_i | D(m) |^{\frac{2}{N}+\frac{2}{p_i+1}-1}T_{i}^2 \left( \left\| \hat{u}_{-}\right\| _{L^{p_{i}+1}}+| \varOmega |^{\frac{1}{N}+\frac{1}{p_i+1}-\frac{1}{2}}T_{i}\rho \right) ^{p_{i}-1},\\ \mathcal C_2&:= 1 - \frac{\lambda }{A_{1,N}} |D(m)|^{\frac{2}{N}}, \end{aligned}$$

where $T_{i} := T_{p_i+1,N}$ is defined by (24) and $A_{1,N}$ is the constant of (22). If $\mathcal C_1 < \mathcal C_2$, u is nonnegative.

Table 4 Calculation examples of $\mathcal C_1$ and $\mathcal C_2$ for several nonlinearities f

Table 4 shows calculation examples of $\mathcal C_1$ and $\mathcal 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)=\lambda t + |t|^{p-1}t$ with a subcritical exponent $p \in (1, p^*)$, positive solutions are admitted only when $\lambda <\lambda _{1}(\varOmega )$. Therefore, calculating |D(m)| can be avoided if we can evaluate $\lambda _{1}(\varOmega )$ with sufficient accuracy. Note that $\lambda _{1}(\varOmega )$ can be obtained analytically, such as when $\varOmega$ is a hyperrectangle. When $\lambda =0$, |D(m)| does not need to be calculated to derive $\mathcal C_2$. In cases where we do not require |D(m)| to estimate $\mathcal C_2$, it is reasonable to replace |D(m)| with $|\varOmega |$ in calculating $\mathcal C_1$ to reduce calculation cost especially when $\rho$ is sufficiently small.

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 Sects. 4 and 5, kv library [12] and INTLAB [30] have such a rigorous function. Therefore, we are left to estimate upper bounds for $\left\| \hat{u}_{-}\right\| _{L^p}$ and |D(m)|, and a lower bound for $\left\| \hat{u}_{+}\right\| _{L^{p}(D(m))}$. We describe the ways to calculate these bounds in the following subsections. For sufficient accuracy, we divide $\overline{\varOmega }$ into a union of small subsets $\left\{ K_{i}\right\} _{i=1}^{N_K}$ that satisfies $\cup _{i}\overline{K_{i}}=\overline{\varOmega }$ and $\mathrm{measure\,} (\overline{K_{i}} \cap \overline{K_{j}})=0$ for $i\ne j$. We can obtain such a division “for free” such as when using finite element methods to compute $\hat{u}$. Otherwise, we should create a mesh $\left\{ K_{i}\right\} _{i=1}^{N_K}$ of $\overline{\varOmega }$ that satisfies the above property. When $\varOmega$ is polygonal, a convenient way is to divide $\overline{\varOmega }$ into a rectangular or triangular mesh. In preparation, we estimate a lower bound $m_i$ and an upper bound $M_i$ for $\hat{u}$ on each mesh $\overline{K_{i}}$, obtaining a closed interval $[m_i, M_i]$ that encloses $[\displaystyle \min _{x \in \overline{K_i}} \hat{u}(x), \displaystyle \max _{x \in \overline{K_i}} \hat{u}(x)]$. When we use a linear finite element approximation, $m_i$ and $M_i$ are the minimal and maximal values at vertexes for each $K_i$, respectively. We denote by $\overline{\Lambda }_m$ and $\underline{\Lambda }_m$ the sets of indices i such that $M_i \le m$ and $m_i \le m$, respectively.

6.1 Upper bound of $\left\| \hat{u}_{-}\right\| _{L^p}$ for $p \in (1,\infty )$

To calculate an upper bound of $\left\| \hat{u}_{-}\right\| _{L^p}$, we use the inequality

$$\begin{aligned} \left\| \hat{u}_{-}\right\| _{L^{p}(\varOmega )} = \left( \int _{{\mathrm{supp}}\,{u_{-}}} \hat{u}_{-}(x)^{p}dx\right) ^{\frac{1}{p}} \displaystyle \le \left( \max _{x\in {\mathrm{supp}}\,{u_{-}}} \hat{u}_{-}(x) \right) \left| {\mathrm{supp}}\,{u_{-}}\right| ^{\frac{1}{p}}. \end{aligned}$$

(32)

Upper bounds for $\displaystyle \max _{x\in {\mathrm{supp}}\,{u_{-}}}\hat{u}_{-}(x)$ and $\left| {\mathrm{supp}}\,{u_{-}}\right|$ can be obtained as follows:

$$\begin{aligned}&\displaystyle \max _{x\in {\mathrm{supp}}\,{u_{-}}} \hat{u}_{-}(x) \le \displaystyle | \min _{i} [ \min \left\{ 0,m_i \right\} ] |,\\&\left| {\mathrm{supp}}\,{u_{-}}\right| \le \displaystyle \sum _{m_i < 0} \left| K_i \right| . \end{aligned}$$

Thus, we can obtain the desired lower bound according to (32).

6.2 Upper bound of |D(m)|

Recall that $D(m):= \{ x \in \varOmega : \hat{u}(x)\le m\}$. When $D(m) \cap \overline{K_i} \ne \phi$, we see that $m_i \le m$.

Therefore, it follows from the definition of $\underline{\Lambda }_m$ that

$$\begin{aligned} D(m)\subset \bigcup _{i\in \underline{\Lambda }_m}K_{i} ~~~\text {and}~~~ |D(m)| \le \displaystyle \sum _{i \in \underline{\Lambda }_m} |K_i|. \end{aligned}$$

6.3 Lower bound of $\left\| \hat{u}_{+}\right\| _{L^{p}(D(m))}$ for $p \in (1,\infty )$

It should be noted that $\left\| \hat{u}_{+}\right\| _{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 $\left\| \hat{u}_{-}\right\| _{L^p}$.

Therefore, we use the following estimation:

$$\begin{aligned}&\left\| \hat{u}_{+}\right\| _{L^{p}(D(m))} = \left( \int _{D(m)} \hat{u}_{+}(x)^{p}dx\right) ^{\frac{1}{p}}\nonumber \\&\ge \left( \sum _{i\in \overline{\Lambda }_{m}}\int _{K_{i}} \hat{u}_{+}(x)^{p}dx\right) ^{\frac{1}{p}} \ge \left( \sum _{i\in \overline{\Lambda }_{m}}|K_{i}|\left( \min _{x\in K_{i}}\hat{u}_{+}(x)\right) ^{p}\right) ^{\frac{1}{p}}, \end{aligned}$$

(33)

where $\displaystyle \bigcup _{i\in \overline{\Lambda }_m}K_{i}\subset D(m)$.

By applying $\displaystyle \min _{x\in K_{i}}\hat{u}_{+}(x) \ge \max \{0,m_i\}$ to each $i \in \overline{\Lambda }_m$, estimation of $\left\| \hat{u}_{+}\right\| _{L^{p}(D(m))}$ from below is completed.

7 Verification theory

In this section, we prepare verification theory to prove the existence of solutions u of (13) satisfying (2) and apply our method to verifying the positivity of u.

We consider a square domain $\varOmega =(0,1)^2$, and an L-shaped domain $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$ where $H^2$-regularity of solutions is lost due to the re-entrant corner at (0.5, 0.5). To obtain $H^1_0$-estimation (2), we use the following affine invariant Newton–Kantorovitch theorem. We omit the notation $(\varOmega )$ when expressing operator norms just to save space. For example, we abbreviate $||\cdot ||_{\mathcal{L}(H^1_0(\varOmega ),H^1_0(\varOmega ))}$ as $||\cdot ||_{\mathcal{L}(H^1_0,H^1_0)}$. We denote $B(\hat{u},r):= \{v\in H^1_0(\varOmega ): \left\| v-\hat{u}\right\| _{H^1_0}<r\}$ and $\overline{B}(\hat{u},r):= \{v\in H^1_0(\varOmega ): \left\| v-\hat{u}\right\| _{H^1_0}\le r\}$ for $r>0$.

Theorem 7.1

([7]) Let $\hat{u} \in H^1_0(\varOmega )$ be some approximate solution of $\mathcal F(u)=0$. Suppose that there exists some $\alpha >0$ satisfying

$$\begin{aligned} ||{\mathcal F'_{\hat{u}}}^{-1}\mathcal F(\hat{u})||_{H^1_0}\le \alpha . \end{aligned}$$

(34)

Moreover, suppose that there exists some $\beta >0$ satisfying

$$\begin{aligned} ||{\mathcal F'_{\hat{u}}}^{-1}(\mathcal F'_v-\mathcal F'_w)||_{\mathcal {L}(H^1_0,H^1_0)}\le \beta ||v-w||_{H^1_0}~~\mathrm {for~all}~ v,w \in D,\end{aligned}$$

(35)

where $D=B(\hat{u},2\alpha +\delta )$ is an open ball depending on the above value $\alpha >0$ for small $\delta >0$. If

$$\begin{aligned} \alpha \beta \le \frac{1}{2}, \end{aligned}$$

(36)

then there exists a solution $u \in H^1_0(\varOmega )$ of $\mathcal F(u)=0$ in $\overline{B}(\hat{u},\rho )$ with

$$\begin{aligned} \rho = \frac{1-\sqrt{1-2\alpha \beta }}{\beta }. \end{aligned}$$

(37)

Furthermore, $\mathcal F'_{\varphi }$ is invertible for every $\varphi \in B(\hat{u},\rho )$, and the solution u is unique in $\overline{B}(\hat{u},2\alpha )$.

We set $\alpha$ and $\beta$ to upper bounds of

$$\begin{aligned} \Vert \mathcal {F}_{\hat{u}}^{\prime -1}\Vert _{\mathcal {L}\left( H^{-1}, H^1_0\right) }\Vert \mathcal {F}(\hat{u})\Vert _{H^{-1}} ~~~\text {and}~~~ \Vert \mathcal {F}_{\hat{u}}^{\prime -1}\Vert _{\mathcal {L}\left( H^{-1}, H^1_0\right) } L,\end{aligned}$$

(38)

respectively, then applying Theorem 7.1 to prove the local existence of solutions. Here, L is a positive number satisfying

$$\begin{aligned} \left\| F_{v}^{\prime }-F_{w}^{\prime }\right\| _{\mathcal {L}(H^1_0, H^{-1})} \le L\Vert v-w\Vert _{H^1_0} ~~\text {for all}~ v, w \in D. \end{aligned}$$

(39)

We estimate the inverse norm $\Vert \mathcal {F}_{\hat{u}}^{\prime -1}\Vert _{\mathcal {L}\left( H^{-1}, H^1_0\right) }$ using the method described in [20, 38] in a finite-dimensional subspace $V_M \subset H^1_0(\varOmega )$ specified later.

7.1 Square domain

For the square domain $\varOmega = (0,1)^2$, we construct approximate solutions $\hat{u}$ with a Legendre polynomial basis. More precisely, we construct $\hat{u}$ as

$$\begin{aligned} \displaystyle \hat{u}(x,y)=\sum _{i=1}^{M}\sum _{j=1}^{M}u_{i,j}\phi _{i}(x)\phi _{j}(y),~~u_{i,j}\in \mathbb {R}, \end{aligned}$$

(40)

where each $\phi _{n}$ ($n=1,2,3,\ldots$) is defined by

$$\begin{aligned} \displaystyle&\phi _{n}(x)=\frac{1}{n(n+1)}x(1-x)\frac{dQ_{n}}{dx}(x)\nonumber \\&\text { with } Q_{n}(x)=\displaystyle \frac{(-1)^{n}}{n!}\left( \frac{d}{dx}\right) ^{n}x^{n}(1-x)^{n},~~n=1,2,3,\ldots . \end{aligned}$$

(41)

The upper bound on $\Vert \mathcal {F}(\hat{u})\Vert _{H^{-1}}$ is evaluated via $C_2 \Vert \mathcal {F}(\hat{u})\Vert _{L^2}$, where $C_2$ is the constant of embedding $L^{2}(\varOmega )\hookrightarrow H^{-1}$, which in fact coincides with the constant of embedding $H^1_0(\varOmega )\hookrightarrow L^{2}(\varOmega )$ (see, for example, [29]). This $L^2$-norm is computed using a numerical integration method with a strict estimation of rounding errors using [12]. For $\varOmega = (0,1)^2$, as mentioned for (25), the embedding constant $C_2$ is calculated as $C_2=(2 \pi ^2)^{-\frac{1}{2}} \approx 0.2251$ with a strict estimation of rounding errors.

We define a finite-dimensional subspace $V_M~(\subset H^1_0(\varOmega ))$ as the tensor product $V_M = \text {span }\,\{\phi _{1},\phi _{2},\ldots ,\phi _{M}\} \otimes \text {span }\,\{\phi _{1},\phi _{2},\ldots ,\phi _{M}\}$, then define the orthogonal projection $P_M$ from $H^1_0(\varOmega )$ to $V_M$ as

$$\begin{aligned} (v- P_M v,v_M)_{H^1_0} = 0 \text { for all } v \in H^1_0(\varOmega )\text { and } v_M \in V_M. \end{aligned}$$

(42)

We use [13, Theorem 2.3] to obtain an explicit interpolation-error constant $C_M$ satisfying

$$\begin{aligned} \left\| v-P_{M}v\right\| _{H^1_0}\le C_{M}\left\| \varDelta v\right\| _{L^{2}} \text { for all } v \in H^1_0(\varOmega )\cap H^2(\varOmega ), \end{aligned}$$

(43)

then applying [20, 38] to estimate the inverse norm $\Vert \mathcal {F}_{\hat{u}}^{\prime -1}\Vert _{\mathcal {L}\left( H^{-1}, H^1_0\right) }$.

7.2 L-shaped domain

For the L-shaped domain $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$, we set $V_M$ to a finite element space of piecewise quadratic basis functions with the non-uniform triangulation displayed in Fig. 1, constructing approximate solutions $\hat{u} \in V_M$.

Using [21, Theorem 3.3], we confirmed that $C_M = 0.011437$ satisfies

$$\begin{aligned} \left\| u_{g}-P_{M}u_{g}\right\| _{H^1_0}\le C_{M}\left\| g\right\| _{L^{2}}~~~\mathrm{for~all}~g\in L^2 (\varOmega ), \end{aligned}$$

(44)

where $P_M$ is the orthogonal projection from $H^1_0(\varOmega )$ to $V_M$ defined in (42), and $u_g\in {H^1_0(\varOmega )}$ is a unique solution of the weak formulation of the Poisson equation

$$\begin{aligned} \left( u_g, v\right) _{H^1_0}=\left( g,v\right) _{L^2}~~~\mathrm{for~all}~v\in H^1_0 (\varOmega ) \end{aligned}$$

(45)

given $g\in L^{2}\left( \varOmega \right)$. The upper bound on $\Vert \mathcal {F}(\hat{u})\Vert _{H^{-1}}$ is evaluated using the Raviart-Thomas mixed finite element method (see, for example, [31]). We estimate the inverse norm $\Vert \mathcal {F}_{\hat{u}}^{\prime -1}\Vert _{\mathcal {L}\left( H^{-1}, H^1_0\right) }$ using (44) and the method described in [20, 38].

Fig. 1

A non-uniform mesh for the L-shaped domain $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$. The number of triangular elements is 3754. The mesh size around the non-convex corner is about a quarter of the other corners

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^{\frac{2}{3}} \sin \left( \frac{2}{3}\theta \right)$, where $(r,\theta )$ 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.

8 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 $\times$ 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:

• $\cdot$ $M_u$: number of Legendre basis functions on $\varOmega = (0,1)^2$ with respect to x and y for constructing approximate solution $\hat{u} \in V_{M_u}$ (see (40))

• $\cdot$ M: number of Legendre basis functions on $\varOmega = (0,1)^2$ with respect to x and y for calculating $\Vert F'^{-1}_{\hat{u}}\Vert _{\mathcal {L}(H^{-1},H^1_0)}$ (see (40))

• $\cdot$ $\Vert F'^{-1}_{\hat{u}}\Vert$: operator norm $\Vert F'^{-1}_{\hat{u}}\Vert _{\mathcal {L}(H^{-1},H^1_0)}$ required in (38)

• $\cdot$ $\Vert F(\hat{u})\Vert$: residual norm $\Vert F(\hat{u})\Vert _{H^{-1}}$ required in (38)

• $\cdot$ L: upper bound for Lipschitz constant satisfying (39)

• $\cdot$ $\alpha$ and $\beta$: constants required in Theorem 7.1

• $\cdot$ $\rho$: error bound $\Vert u- \hat{u} \Vert _{H^1_0}$

• $\cdot$ m: constant that determines D(m); see Lemma 3.2

• $\cdot$ $|{\mathrm{supp}}\,u_{-}|$: volume of support of $u_{-}$

• $\cdot$ $\lambda _{1}({\mathrm{supp}}\,u_{-})$: first eigenvalue of $-\varDelta$ on interior of ${\mathrm{supp}}\,u_{-}$ defined by (7)

• $\cdot$ $\mathcal {C}_1$ and $\mathcal {C}_2$: constants required in Theorem 6.5

8.1 Lane–Emden equation

As mentioned in the Introduction, positive solutions of the Lane–Emden equation with subcritical $p>1$,

$$\begin{aligned} \left\{ \begin{array}{l l} -\varDelta u=u\left| u\right| ^{p-1} &{}\mathrm {in~} \varOmega ,\\ u=0 &{}\mathrm {on~} \partial \varOmega \end{array}\right. \end{aligned}$$

(46)

have been studied from various points of view (again, see [4, 6, 10, 11, 18]). 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

$$\begin{aligned} L \le p(p-1) C_{p+1}^{3}\left( \Vert \hat{u}\Vert _{L^{p+1}}+C_{p+1} r\right) ^{p-2},~~r=2\alpha +\delta ~~\mathrm{for~small}~~\delta >0 \end{aligned}$$

via a simple calculation from the definition, where we set r to be the next floating-point number after $2\alpha$. We refer to [33, Sect. 4] as a numerical experiment for positive solutions of the Lane–Emden equation for $p=3, 5$ on $\varOmega =(0,1)^2$, where the volume of D(m) was roughly estimated as $|D(m)|\le |\varOmega |$.

For the L-shaped domain $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$, we constructed an approximate solution $\hat{u} \in V_M$ of (46) with $p=3$ using a piecewise quadratic basis, obtaining Fig. 2. Table 5 shows the verification result and confirms the positivity of the enclosed solution because $\mathcal {C}_1 \le \mathcal {C}_2$.

Fig. 2

An approximate solution of (46) for $p=3$ on $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$

Table 5 Verification results for (46) for $p=3$ on $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$

8.2 Allen–Cahn equation and Nagumo equation

In the next example, we consider the stationary problem of the Allen–Cahn equation

$$\begin{aligned} \left\{ \begin{array}{l l} -\varDelta u=\lambda (u-u^3) &{}\mathrm {in~} \varOmega ,\\ u=0 &{}\mathrm {on~} \partial \varOmega \end{array}\right. \end{aligned}$$

(47)

with $\lambda >0$. This is regarded as a special case of the stationary problem of the Nagumo equation

$$\begin{aligned} \left\{ \begin{array}{l l} -\varDelta u=\lambda u (1-u) (u-a)= \lambda (-au+(1+a)u^2-u^3) &{}\mathrm {in~} \varOmega ,\\ u=0 &{}\mathrm {on~} \partial \varOmega \end{array}\right. \end{aligned}$$

(48)

with $\lambda >0$ and $0<a<1$. By applying $u=(v+1)/2$ to (48) with $a=0.5$ and adjusting the value of $\lambda$, these equations become identical. The Lipschitz constants L satisfying (39) for (47) and (48) were respectively estimated as

$$\begin{aligned}&L \le 6 \lambda C_{4}^{3}\left( \Vert \hat{u}\Vert _{L^{4}}+C_{4} r\right) ,&~~r=2\alpha +\delta ~~\mathrm{for~small}~~\delta>0,\\&L \le \lambda \left( 2 (1+a) C_{3}^{3} + 6 C_{4}^{3}\left( \Vert \hat{u}\Vert _{L^{4}}+C_{4} r\right) \right) ,&~~r=2\alpha +\delta ~~\mathrm{for~small}~~\delta >0, \end{aligned}$$

where we set r to be the next floating-point number after $2\alpha$. 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 $\hat{u}$.

Fig. 3

Approximate solutions of (47) on $\varOmega =(0,1)^{2}$ for $\lambda =100$, 400, and 1600

Table 6 Verification results for (47) on $\varOmega =(0,1)^{2}$ for $\lambda =100$, 400, and 1600

Theorem 6.5 can be uniformly applied to the nonlinearities of the form (4) (again, see Table 4).

For the square domain $\varOmega =(0,1)^2$, we constructed approximate solutions $\hat{u}$ using a Legendre polynomial basis. For (47) ($\lambda =100$, 400, and 1600), we obtained Fig. 3 and the verification results in Table 6. For (48) ($\lambda =400$, $a=1/4$), we found multiple solutions displayed in Fig. 4 and the verification results in Table 7.

For the L-shaped domain $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$, we constructed approximate solutions $\hat{u} \in V_M$ using a piecewise quadratic basis, obtaining Fig. 5 for (47) ($\lambda =100$) and (48) ($\lambda =180$, $a=1/64$), and the verification results in Table 8. In all cases, we confirmed $\mathcal {C}_1 \le \mathcal {C}_2$ and thus the positivity of the enclosed solutions.

Fig. 4

Approximations of multiple solutions of (48) on $\varOmega =(0,1)^{2}$ for $\lambda =400$ and $a=1/4$

Table 7 Verification results for (48) on $\varOmega =(0,1)^{2}$ for $\lambda =400$, $a=1/4$

Fig. 5

Approximate solutions of (47) and (48) on $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$

Table 8 Verification results for (47) and (48) on $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$

8.3 Elliptic equation with multiple positive solutions

For the last example, we consider the elliptic boundary value problem

$$\begin{aligned} \left\{ \begin{array}{l l} -\varDelta u=\lambda (u+{\mathrm{A}}u^2 - {\mathrm{B}}u^3) &{}\mathrm {in~} \varOmega ,\\ u=0 &{}\mathrm {on~} \partial \varOmega , \end{array}\right. \end{aligned}$$

(49)

given $\lambda , {\mathrm{A}},{\mathrm{B}}>0$. This problem has two positive solutions when $\lambda ^*< \lambda < \lambda _1(\varOmega )$ for a certain $\lambda ^*>0$ (see [19]). The Lipschitz constant L for this problem can be estimated as

$$\begin{aligned} L \le \lambda \left( 2 {\mathrm{A}}C_{3}^{3} + 6 {\mathrm{B}}C_{4}^{3}\left( \Vert \hat{u}\Vert _{L^{4}}+C_{4} r\right) \right) ,~~r=2\alpha +\delta ~~\mathrm{for~small}~~\delta >0, \end{aligned}$$

where we again set r to be the next floating-point number of $2\alpha$.

To prove the positivity of solutions of (49), the previous method [33] requires an $L^{\infty }$-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(\varOmega ) \hookrightarrow L^{\infty }(\varOmega )$ successfully. However, the proposed method is well applicable even to (49) without assuming $H^2$-regularity; see again the last case in Table 4.

For the square domain $\varOmega =(0,1)^2$, we constructed approximations $\hat{u}$ of multiple positive solutions of (49) ($\lambda =10$, ${\mathrm{A}}=5$, ${\mathrm{B}}=1$) using a Legendre polynomial basis, obtaining the results in Fig. 6 and Table 9.

For the L-shaped domain $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$, we constructed multiple approximate solutions $\hat{u} \in V_M$ of (49) ($\lambda =20$, ${\mathrm{A}}=5$, ${\mathrm{B}}=1$) using a piecewise quadratic basis, obtaining Fig. 7 and Table 10. Since $\mathcal {C}_1 \le \mathcal {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.

Fig. 6

Approximations of multiple solutions of (49) on $\varOmega =(0,1)^{2}$ when $\lambda =10$, ${\mathrm{A}}=5$, ${\mathrm{B}}=1$

Table 9 Verification results for (49) on $\varOmega =(0,1)^{2}$ when $\lambda =10$, ${\mathrm{A}}=5$, ${\mathrm{B}}=1$

Fig. 7

Approximations of multiple solutions of (49) on $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$ when $\lambda =20$, ${\mathrm{A}}=5$, ${\mathrm{B}}=1$

Table 10 Verification results for (49) on $\varOmega =(0,1)^2\backslash ([0,0.5]\times [0.5,1])$ when $\lambda =20$, ${\mathrm{A}}=5$, ${\mathrm{B}}=1$

9 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 estimation (2) given some numerical approximation $\hat{u}$ and an explicit error bound $\rho$. 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.

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}=\sqrt{2/(\pi x)} \sin x$; therefore, the first zero of this is $\pi$ (see, for example, [3, Remark 1.2] and [27, Sect. 10.16]). Moreover, we have $j_{1.5,1}=\sqrt{2/(\pi x)} \left( x^{-1} \sin x - \cos x \right)$, 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 $[\pi ,1.5\pi ]$, 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

$$\begin{aligned} j_{n,1}(x) = \frac{1}{\pi } \int _{0}^{\pi } \cos (x\sin t-nt) dt; \end{aligned}$$

see [27, Sect. 10.9]. Let $a\ge 0$ and suppose $0 \le x \le a$ so that $0 \le x \sin t \le a$ ($0\le t \le \pi$). Let us write $A:=[0,a]$. Then, we have

$$\begin{aligned} j_{0,1}(x)&= \frac{1}{\pi } \int _{0}^{\pi } \cos (x\sin t) dt \in \cos A = [\cos a, 1]. \end{aligned}$$

When $a \in [0,1] \subset [0,\pi /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

$$\begin{aligned} \frac{d}{dx}j_{1,1}(x) = \frac{1}{\pi } \int _{0}^{\pi } \sin {t}\sin (t-x\sin t) dt. \end{aligned}$$

Hence, we have, for $a\in [0,\pi ]$,

$$\begin{aligned} \frac{d}{dx}j_{1,1}(x) \in \frac{1}{\pi } \int _{0}^{\pi } \sin {t}\sin (t-A) dt = \frac{1}{4 \pi } \left( 2 \pi \cos A + \sin A - \sin A \right) =:B \end{aligned}$$

For all $t \in B$, we confirm $t \ge \frac{1}{4 \pi } \left( 2 \pi \cos a - \sin a \right)$. This value is positive when $a \in [0,1] \subset [0,\tan ^{-1}(2 \pi ))$. Therefore, $j_{n,1}$ monotonically increases for $x \in [0,1]$.

Using bisection steps, we rigorously computed the first positive zeros of $j_{0,1}$ and $j_{1,1}$. We first found a compact interval that includes the first positive zero of $j_{n,1}$ by searching the first interval $[\underline{x},\overline{x}]:=[1+(i-1)\varepsilon /2,1+(i-1)\varepsilon /2+\varepsilon ]$ ($i=1,2,\ldots$) that satisfies $j_{i,1}(\underline{x})j_{i,1}(\overline{x})<0$. Then, starting from the center of $[\underline{x},\overline{x}]$, we repeated the bisection method until the desired precision was achieved. All rounding errors were strictly estimated using kv library version 0.4.49 [12].