Rigorous numerical inclusion of the blow-up time for the Fujita-type equation

Multiple studies have addressed the blow-up time of the Fujita-type equation. However, an explicit and sharp inclusion method that tackles this problem is still missing due to several challenging issues. In this paper, we propose a method for obtaining a computable and mathematically rigorous inclusion of the L2(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\varOmega )$$\end{document} blow-up time of a solution to the Fujita-type equation subject to initial and Dirichlet boundary conditions using a numerical verification method. More specifically, we develop a computer-assisted method, by using the numerically verified solution for nonlinear parabolic equations and its estimation of the energy functional, which proves that the concerned solution blows up in the L2(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\varOmega )$$\end{document} sense in finite time with a rigorous estimation of this time. To illustrate how our method actually works, we consider the Fujita-type equation with Dirichlet boundary conditions and the initial function u(0,x)=1925x(x-1)(x2-x-1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u(0,x)=\frac{192}{5}x(x-1)(x^2-x-1)$$\end{document} in a one-dimensional domain Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega$$\end{document} and demonstrate its efficiency in predicting L2(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\varOmega )$$\end{document} blow-up time. The existing theory cannot prove that the solution of the equation blows up in L2(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\varOmega )$$\end{document}. However, our proposed method shows that the solution is the L2(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\varOmega )$$\end{document} blow-up solution and the L2(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\varOmega )$$\end{document} blow-up time is in the interval (0.3068, 0.317713].


Introduction
Let ⊂ ℝ N (N ∈ ℕ) be a bounded Lipschitz domain. In this paper, we propose a numerical verification algorithm for obtaining an explicit and computable inclusion (an interval containing the exact value) of the blow-up time of solutions to the Fujita-type equation: For 1 ≤ q < ∞ , L q ( ) is the set of the q-th power Lebesgue integrable real-valued functions in endowed with the norm ‖u‖ L q ( ) ∶= � ∫ �u(x)� q dx � 1 q . The L 2 ( ) inner product is defined as (u, v) L 2 ( ) ∶= ∫ u(x)v(x)dx for u, v ∈ L 2 ( ) . The t v(t) ∶= ( t v)(t, ⋅) , where t denotes the weak derivative of t ∈ J . For a real Banach space X, the function space C(J; X) is defined by the set of continuous functions as J ∋ t ↦ v(t) ∈ X . Let C 1 (J;X) be the set of C 1 (J) functions as J ∋ t ↦ v(t) ∈ X . For a time interval (a, b) (0 ≤ a < b < ∞) , we define the function space: Assume that u 0 ∈ H 1 0 ( ) and 1 < p < N+2 N−2 , where we replaced N+2 N−2 by ∞ for N = 1, 2 . Recall that we consider a blow-up solution of the Fujita-type equation: Fujita pioneered in studies of a solution for (2) with = ℝ N in 1966 [6]. This seminal work has inspired many researchers to study solutions of variable timeevolution partial differential equations including (2) (see e.g., [13,26,27] and references therein). Currently, the existence of solutions of (2) in a bounded time (1) ⟨Au, v⟩ = (u, v) H 1 0 ( ) ∀v ∈ H 1 0 ( ).

3
Rigorous numerical inclusion of the blow-up time for the Fujita-type equation interval [0, T) (T < ∞) has been proven. For instance, under some assumptions, we can prove the existence of the solution in Z (0,T) (see e.g., Theorem 3.1 in [1]). For the existence of a blow-up solution of (2), the sufficient conditions for the initial function u 0 should be satisfied. It can be achieved using the first eigenvalue of the operator A defined in (1) [9], the energy functional (see e.g., [1,25]), and an unstable stationary solution of (2) (see e.g., [12]), etc. In the numerical analysis, many approximations of the blow-up time have been suggested and proved to converge to the exact blow-up time (see, e.g., [3,4,17]). Computer-assisted methods for finding the exact value of the blow-up time of the solution for ordinary differential equations have also been established (see, e.g., [14,15,23]). However, the same problem for the partial differential equations has not been properly addressed so far. We define L 2 ( ) blow-up solution as a solution of (2) satisfying for some t * max > 0 . The t * max is also defined as L 2 ( ) blow-up time. We provide an explicit and computable inclusion of the L 2 ( ) blow-up time for the solution u of (2) in Z (0,t * max ) . We propose a time-stepping algorithm using numerical verification methods for obtaining a sharp inclusion of the L 2 ( ) blowup time. Our proposed algorithm is initialized at t 0 = 0 , i = 0 and includes the following steps: Step 1: Estimate an upper bound t max of the L 2 ( ) blow-up time using the initial value u(t i ) . Then, note that the true L 2 ( ) blow-up time t * max is included in the interval (t i , t max ] . We proceed to Step 2 even if t max is not obtained. Step 2: Prove the existence of the solution in the time interval [t i , t i+1 ] using a computer-assisted proof and obtain the solution u(t i+1 ) at the time t i+1 .
Step 3: Replace i with i + 1 and return to Step 1.
Note that in Step 3, we update the initial value of the Fujita-type equation from u(t i ) to u(t i+1 ) . Hereinafter, i ∈ ℕ ∪ {0} in Steps 1, 2, and 3 are called the "step number". In Step 2, we prove that the solution exists in time interval [t i , t i+1 ] by computer-assisted proofs of the existence for solutions of parabolic equations. Since the first derivation by Nakao [18], several methods of the computer-assisted proofs of the existence for solutions of parabolic equations have been proposed. These methods have been improved with time and can now verify the existence of solutions for parabolic equations in a long-time interval J ⊂ ℝ [7,24]. In this paper, we use the method in [7] (see Appendix B for details) in Step 2.
Here, we focus on Step 1. In Step 1, we introduce the L 2 ( ) blow-up time estimate using the energy functional E ∶ H 1 0 ( ) ∩ L p+1 ( ) → ℝ defined below for obtaining the upper bound t max : Several well-known theorems provide some sufficient conditions concerned with the value of E(u 0 ) , where u 0 is an initial function of (2), for proving that the solution u of (2) is the L 2 ( ) blow-up solution. However, the exact value of the upper bound of the L 2 ( ) blow-up time is not obtained in the well-known theorems (see e,g., [1]). In Theorem 1, we formulate a theorem for giving an explicit and computable inclusion of the L 2 ( ) blow-up time as well as verifying the existence of L 2 ( ) blow-up solution by means of numerical verification methods: Note that we can obtain an inclusion of the L ∞ ( ) blow-up time if the inequality (3) holds under some assumptions (see Appendix A for details).
In Step 1, we first derive the inclusion of ‖u(t i )‖ L 2 ( ) and E(u(t i )) using a numerical verification method (see Appendix C for details). Moreover, if the sufficient condition (3) is satisfied, we can compute the upper bound t max of the L 2 ( ) blow-up time using (4).
We can verify the existence of the solution in time interval [0, t i+1 ] and demonstrate it blowing up for t > t i+1 using iterative procedure. Therefore, our method can verify whether the solution blows up even if we cannot determine that the blowingup occurs from the initial function u 0 . We provide the detailed algorithm description in Algorithm 1.
This paper is organized as follows. In Sect. 2, we prove Theorem 1. In Sect. 3, we introduce several numerical examples obtained by the Algorithm 1.

Proof of Theorem 1
In the below, we first present some auxiliary lemmas and a corollary, which are used to prove Theorem 1. Corollary 1 is well-known with a smooth domain (see e.g., [20,Lemma 17.5]). We show Lemma 1 for proving Corollary 1 for the case, when the bounded domain has Lipschitz boundary.
We introduce the well-known Lemma 3, which leads an upper bound of the blow-up time in variable topologies for solutions of several differential equations (see e.g., [11,16,28]). We provide Lemma 2 for proving Lemma 3:

positive and continuous functions. Assume that b is non-decreasing and functions u,v:
Then we have the following lemma: Furthermore, we assume that , and x 0 = x 1 = y(t M ) , respectively in Lemma 2, Lemma 2 implies that the function y satisfies (9). Now, let us prove Theorem 1.

Proof of Theorem 1
Provided that the solution u ∈ Z (0,t i ) of (2) exists, we consider the Fujita-type equation (2) with the initial function . First, assuming that the solution u ∈ Z J ∞ exists, we show a contradiction. Let t ∈ J ∞ . Note that the solution u ∈ Z J ∞ satisfies (see, e.g., [5, Theorem 3 in Section 5.9]) and t can be treated as the usual derivative and integrating in , we obtain From (10) and (11) it follows that

3
Here we accounted for (12). (3). Therefore, the solution u does not exists in Z J ∞ . This is a contradiction. Next, we derive an upper bound of the L 2 ( ) blow-up time t * max . Setting s = t −2 , we obtain It follows from (13) that the upper bound t max of the t * max can be derived from

Numerical results
We consider the Fujita-type equation with p = 2: , and u 0 > 0 means that u 0 (x) > 0 a.e. x ∈ . We verify that the solution u of (14) blows up in the L 2 ( ) sense and obtain explicit numerical inclusions of the L 2 ( ) blow-up time using Algorithm 1. Let V h and V k be the sets depending on h > 0 and k > 0 of the piecewise Hermite spline ( C 1 -class with 5-degree) functions in the domain and the piecewise quadratic ( C 0 -class) functions in J, respectively. We define V h,k ∶= V h ⊗ V k and construct a full-discrete approximation u h,k ∈ V h,k of (14). The method in [7] verifies the existence of the solution u ∈ Z J of (14) (see Appendix B for details). If the method in [7] succeeds to verify, we can compute explicit error bounds L 2 and H 1 0 satisfying that , we can also compute an explicit error bound E err satisfying that We have an upper bound of the energy functional E(u(t i )) and a lower bound of the L 2 ( ) norm of the forms: Then, we can check that the condition (3) holds or not, where we write the details of the procedure of checking the condition of (3) in Appendix C.
In this section, we put h = 1∕64 and k = (t i+1 − t i )∕128 for the step numbers i and set the minimal time-step width min = 10 −4 in Algorithm 1. All computations were carried out on the Dell Precision 7920 Intel Xeon Gold 6134 CPU 3.20GHz with MATLAB (var. R2018a) and Symbolic Math Toolbox. Furthermore, we added INTLAB(var. 10.1), which calculated estimates including all rounding errors [22].
We verified the existence of the solution u, which blows up, and led an explicit inclusion of the L 2 ( ) blow-up time using Algorithm 1. For the step numbers i, the detailed estimates of (15), t max − t i , and (t max − t i )∕t i , which were obtained by Algorithm 1, are summarized in Table 1, where the "Failed1" means the failure of checking the sufficient condition (3) in Step 1 and the "Failed2" means that we failed to verify the solution u in time [t i , t i+1 ] when ≥ min in Step 2 in Algorithm 1, respectively. When we started from t = 0 and the step number 0 ≤ i ≤ 8 , it could not be proved that the solution u blows up using the procedure of Step 1 in Algorithm 1; that is, we could not check that the sufficient condition (3) in Theorem 1 holds. The procedure of Step 2 in Algorithm 1 with the computer-assisted method given by [7] started to prove that the solution u exists

3
Rigorous numerical inclusion of the blow-up time for the Fujita-type equation Table 1 The detailed estimates for step numbers i in Sect. 3 in a neighborhood centered at u h,k for the step numbers i. When the step number i ≥ 9 , the procedure of Step 1 in Algorithm 1 proved that the condition (3) in Theorem 1 is satisfied. More specifically, we verified that the solution blows up and obtain the upper bound t max of the L 2 ( ) blow-up time. When i = 52 , Algorithm 1 were aborted because the time-step width satisfied min > in the procedures of Step 2. We wrote the lower bounds t i and the upper bounds t max obtained by Algorithm 1 in Table 2 and plotted them in Fig. 3. It demonstrated that the absolute errors t max − t i became smaller for i ≤ 46 compared with the corresponding values in the previous steps. However, for i ≥ 47 , because of the full-discrete approximate solutions u h,k with low accuracy, the absolute errors have been bigger compared with the corresponding values in the previous steps (Table 1). We have that the sharpest inclusion of the L 2 ( ) blow-up time is in the interval (0.3068, 0.317713] from results corresponding with i = 52 and i = 46 in Table 2, respectively. The inclusion result is mathematically interesting because we cannot show that the solution u in this numerical example blows up in the L 2 ( ) sense using the existing results (see e.g., [20,Theorem 17.6]).

Numerical example 2
Let = (0, 1) . We take the initial function u 0 as (see Fig. 4) The full-discrete approximations u h,k of (14) with the initial function u 0 in time intervals The maximal value of approximation u h,k became smaller in the time interval [0, 0.0188]. However, the value grew for the time t > 0.0188 . We proved that the solution u exists and blows up using Algorithm 1. Moreover, we led an explicit

3
Rigorous numerical inclusion of the blow-up time for the Fujita-type equation  (15), t max − t i , and (t max − t i )∕t i , which were provided by Algorithm 1, are written in Table 3, where the "Failed1" and "Failed2" are defined in Sect. 3.1.
The procedure of Step 1 in Algorithm 1 could not check that the condition (3) in Theorem 1 holds at the time t = 0 . The procedure of Step 2 in Algorithm 1 with the computer-assisted method by [7] started to verify the existence of the solution u in a neighborhood centered at u h,k for the step numbers i. When the step number i ≥ 8 , we could prove that the solution u blows up and compute the upper bound t max of the L 2 ( ) blow-up time using the procedure of Step 1 in Algorithm 1. However, when i = 43 , Algorithm 1 was aborted because the method by [7] could not prove the existence of the solution u for the time t > t i even if we set = min in the procedures of Step 2. For the step numbers i, the lower bounds t i and the upper bounds t max obtained by the method in Algorithm 1 were shown in Table 4 and plotted them in Fig. 6, respectively. It yields that the absolute errors t max − t i become smaller for i ≤ 36 and bigger for i ≥ 37 compared with the corresponding values in the previous steps. For i ≥ 37 , the approximate solutions u h,k with low accuracy prevent from deriving the sharp inclusion of the L 2 ( ) blow-up time ( Table 3). The sharpest inclusion of the L 2 ( ) blow-up time is in the interval (0.2285, 0.257074] from the results with i = 43 and i = 36 in Table 4, respectively.

Conclusion
We have proposed a time-stepping algorithm using numerical verification methods for obtaining explicit inclusion of the L 2 ( ) blow-up time of a solution to the Fujitatype equation. The algorithm consists of three iteration steps.
Step 1 provides a criterion whether the solution u blows up at t > t i as well as a rigorous upper bound of t max . Moreover, we have formulated Theorem 1 for constructing the detailed procedure of Step 1. In Step 2, we have verified the existence of the solution in each time interval [t i , t i+1 ] using a computer-assisted existential proof of solutions for parabolic equations. In Step 3, we have updated the step number i to i + 1 and attempted to verify the L 2 ( ) blow-up solution for time t > t i+1 by Step 1. Our method not only verifies that the solution blows up but also derives the explicit inclusion of the L 2 ( ) blow-up time by iterating the procedures of Steps 1, 2, and 3. Our method remains valid even when we cannot determine whether the solution blows up in finite time from the initial function u 0 . We believe that our method will find applications in other problems involving various partial differential equations. Rigorous numerical inclusion of the blow-up time for the Fujita-type equation  Table 3 The detailed estimates for step numbers i in Sect. 3.  (see e.g., [8] and Theorem 8.2 in [19]), where N ∈ ℕ is the dimension of the domain .
We consider an inclusion of t ∞ . Since u ∈ Z (0,t i ) ⊂ C((0, t i );D(A)) ⊂ C((0, t i ); where t * max is the L 2 ( ) blow-up time, we have a contradiction because H ölder's inequality implies that Therefore, the estimate t ∞ ≤ t * max holds and we obtain that the L ∞ ( ) blow-up time

Appendix B: The verification method in [7]
We briefly present a numerical verification method of solutions for the nonlinear problem (14) with g(u) ∶= |u|u . Note that g(u) = u 2 since u(t) ≥ 0 for all t > 0 on by the assumption u 0 > 0 in (14). Dividing interval J into l subintervals, we use a time evolving algorithm in the below.
In order to get an approximate solution for the problem (14) on × J , we use a finite element subspace (14) is equivalent to the following residual equation: h,k is a residual function for i = 1, … , l. We define the operators L i ∶ Z J i → L 2 J i ;L 2 ( ) by where g � [u (i) h,k ] = 2u (i) h,k denote the Fréchet derivative of g at u (i) h,k . Using the operators L i , a solution ū of the problem (16) can be decomposed as ū = v + w by using solutions v and w of and respectively, where Since the solution v of (17) can be determined independently of w in (18), if the solution v of the linear problem (17) is numerically verified, then the problem (16) can be reduced to find a solution w to the nonlinear problem (18). First, we present a numerical verification method of solutions for nonlinear problems (18). Note that the problem (18) is rewritten as the following fixed-point equation of the map L −1 i : Here, the map L −1 i ∶ L 2 J i ;L 2 ( ) → Z J i is considered as the solution operator for the linear parabolic problem with homogeneous initial condition.
For any positive constants i and i , we define the candidate set W i , i of solutions in (18) as where constants M (i) 1 , M (i) t i and C (i) are numerically determined which presented by Lemma 4 and Theorem 5 in [7]. Moreover, by Lemma 6 in [7], the solution v of (17) is estimated as follows: where (T i ) ∶= exp(− 2 T i ) and (T i ) ∶= √ 1 2 2 (1 − (2T i )) , if = (0, 1). Therefore, defining the function G( i , i ) of i and i satisfying the sufficient condition of (20) is given as follows: Next, we describe several remarks on the verification step from the interval J i to J i+1 . Let * i and * i be two positive numbers satisfying the condition (23). Then there exists a solution w * ∈ W * i , * i ( , J i ) of the problem (18) and the following estimates hold When denoting v * as a solution of the problem (17), the solution u * of the nonlinear problem (14) on × J i can be written by u * = u (i) h,k + v * + w * . Note that the initial condition of the next time-step problem on × J i+1 is given by u * (t i ) = u (i) h,k (t i ) + v * (t i ) + w * (t i ) . Since we take the approximate solution , an initial function i+1 of the problem (16) on × J i+1 is given by Therefore, we can obtain the following estimations:

3
Rigorous numerical inclusion of the blow-up time for the Fujita-type equation

Appendix C: The procedure for checking the condition (3)
For checking the condition (3) with p = 2 , we show the procedure for the upper bound of E(u(t i )) and the lower bound ‖u(t i )‖ 3 L 3 ( ) in the below. Let u ∈ Z j and u h,k ∈ V h,k ( , J i ) be an exact and an approximate solutions of (14) on × J i , respectively. Then we have u(t i ) = u h,k (t i ) + i+1 , where i+1 is defined in (24). Since u(t) ≥ 0 for all t ∈ J i on , it follows that Note that ‖u(t i )‖ H 1 0 ( ) ≤ ‖u h,k (t i )‖ H 1 0 ( ) + H 1 0 by (25). Thus, we can obtain the following estimation.