Saddle-Type Blow-Up Solutions with Computer-Assisted Proofs: Validation and Extraction of Global Nature

In this paper, blow-up solutions of autonomous ordinary differential equations (ODEs) which are unstable under perturbations of initial points, referred to as saddle-type blow-up solutions, are studied. Combining dynamical systems machinery (e.g., compactifications, time-scale desingularizations of vector fields) with tools from computer-assisted proofs (e.g., rigorous integrators, the parameterization method for invariant manifolds), these blow-up solutions are obtained as trajectories on local stable manifolds of hyperbolic saddle equilibria at infinity. With the help of computer-assisted proofs, global trajectories on stable manifolds, inducing blow-up solutions, provide a global picture organized by global-in-time solutions and blow-up solutions simultaneously. Using the proposed methodology, intrinsic features of saddle-type blow-ups are observed: locally smooth dependence of blow-up times on initial points, level set distribution of blow-up times, and decomposition of the phase space playing a role as separatrixes among solutions, where the magnitude of initial points near those blow-ups does not matter for asymptotic behavior. Finally, singular behavior of blow-up times on initial points belonging to different family of blow-up solutions is addressed.


Introduction
Our concern in the present paper is blow-up solutions of the following initial value problem of an autonomous system of ordinary differential equations (ODEs) in R n : where t ∈ [0, T ) with 0 < T ≤ ∞, f : R n → R n is a C 1 function and y 0 ∈ R n . We call a solution y(t) of the initial value problem (1.1) a blow-up solution if t max def = sup t | a solution y ∈ C 1 ([0,t)) of (1.1) exists < ∞.
The maximal existence time t max is then called the blow-up time of (1.1). Blow-up solutions can be seen in many dynamical systems generated by ODEs, or partial differential equations (PDEs) like nonlinear heat equations or Keller-Segel systems. These dynamical systems are categorized as exhibiting finite-time singularities and have been the center of attention of many researchers, who have studied these phenomena from mathematical, physical, numerical viewpoints and so on (e.g., Fila and Matano 2002;Herrero and Velázquez 1997;Mizoguchi 2016;Winkler 2013 from theoretical viewpoints and, e.g., Anada et al. 2017;Berger and Kohn 1988;Cho 2016;Cho et al. 2007; Zhou and Saito 2017 from numerical viewpoints). Fundamental questions for blow-up solutions are whether or not a solution blows up and, if it does, when, where and how it blows up. In general, blow-up phenomena depend on initial points, and rigorously characterizing them as functions of initial points remains non-trivial. A typical approach for studying and proving existence of blow-up solutions is via energy estimates (see, e.g., Fila and Matano 2002), namely inequalities (involving energy functionals associated with the systems) giving sufficient conditions for existence of blow-up. In such cases, relatively large initial data induce finite-time blow-up. However, in general, these criteria do not provide an answer on how large initial points should be to exhibit blow-up and how solutions behave when these criteria are violated. There are several cases where initial points are divided such that solutions through them either exist globally in time or blow-up by means of bounded stationary solutions (e.g., Fujita 1969). A stationary solution with the above property is referred to as the separatrix, which plays a key role in describing asymptotic behavior of solutions. Despite their importance, results about the existence and explicit description of separatrixes are limited. On the other hand, there are also results about the existence of blow-ups in which the magnitude of initial points does not matter. Alternative approaches to the energy estimates have been introduced to prove such blow-ups, but their dependence on initial points remain unknown in many cases, while arguments based on energy estimates easily yield the continuous dependence of blow-up behavior on initial points by continuity of energy functionals. Furthermore, there are also blow-up solutions whose asymptotic behavior is described not only by divergence, but also by complex behavior like oscillations, some of which are mentioned in Sect. 8.1 (Concluding Remarks). Mathematical and physical importance for studying blow-up behavior follow from such rich nature, but their comprehensive understanding are limited to well-known systems like PDEs mentioned above at present. See, e.g., Fila and Matano (2002) and Galaktionov and Vázquez (2002) for more detailed summaries of blow-up problems including another well-known characterization of blow-up solutions by means of (backward) self-similarity.
Meanwhile, the second author has recently proposed a description of blow-up solutions from the viewpoint of dynamical systems (Matsue 2018). More precisely, compactifications of the phase space R n is applied to mapping the infinity onto points on the boundary E of a compact manifold or their tangent spaces denoted by D with ∂D = E. The boundary E shall be called the horizon in this context. Accordingly the vector field (1.1) is transformed to one on the corresponding manifolds, but the behavior of solutions near the boundary E is still singular reflecting the behavior of the original vector field at infinity. The timescale transformation, which shall be called the timescale desingularization, is then introduced to desingularize the singularity of the vector field around E. Consequently, dynamics at infinity can be characterized through the time-transformed vector field, called the desingularized vector field, on D. Standard arguments in the theory of dynamical systems through compactifications show that divergent solutions of (1.1) correspond to global-in-time solutions of the desingularized vector field converging to invariant sets on E. 1 A significant consequence of the preceding studies is that a solution of (1.1) with bounded initial point is a blow-up solution, namely t max < ∞, if the image of the solution through a compactification mentioned above is on the local stable manifold of a hyperbolic equilibrium on E for the desingularized vector field. 2 Simultaneously, the second and the third authors have developed a computerassisted methodology for proving the existence of blow-up solutions for concretely given dynamical systems with rigorous bounds of their blow-up times t max (Matsue and Takayasu 2020a, b;Takayasu et al. 2017). The basic idea is the combination of compactifications as well as timescale desingularizations mentioned above with rigorous integrator of ODEs based on interval (and affine) arithmetic and topological characterizations of asymptotic behavior such as locally defined Lyapunov functions. Evaluation of t max is one of the most important issues in blow-up studies to estimate upper bounds of the existence of solutions, or the onset of finite-time singularities such as ignition in combustion studies (e.g., Dold 1985), while the study is limited even in numerical studies (e.g., Cho 2016). The proposed methodology provides a rigorous and standard way to obtain both lower and upper bounds of t max through dynamics at infinity.
The methodology works successfully for validating profiles and blow-up times of blow-up solutions generated by hyperbolic stable equilibria at infinity, while blow-up generated by unstable equilibria at infinity is not reported yet due to several technical difficulties. Note that there is another work for characterizing blow-up solutions with computer assistance by the first author and his collaborators based on analytic approach (D'Ambrosio et al. 2015) whose detail is briefly mentioned in Sect. 8.1. On the other hand, from the viewpoint of dynamics at infinity itself, namely when the viewpoint of blow-up characterizations is not considered, asymptotic behavior of unstable invariant sets at infinity is quite natural to study toward description of global bounded dynamics (e.g., Dumortier 2006;Dumortier and Herssens 1999;Dumortier et al. 2006;Giraldo et al. 2020;Kokubu and Roussarie 2004). We then believe that blow-up solutions generated by unstable invariant sets at infinity contribute toward the comprehensive understanding of global dynamics, including characteristics such as criteria for the existence, dependence on initial points and analytic information of blow-up times. Despite many mathematical and numerical studies of blow-ups, characterizations and computations of blow-up solutions which are unstable under perturbations of initial points in a standard way are not realistic, because we have to treat the following difficulties simultaneously: • instability of trajectories exhibiting blow-up solutions under perturbations of initial points and • treatment of infinity.
We shall call blow-up solutions exhibiting instability under perturbation of initial points saddle-type blow-up solutions in the present paper, respecting the structure of equilibria at infinity. This fuzzy nature is difficult to characterize clearly in general, while such behavior can be partially observed in several practical problems as mentioned in Sect. 8.1.
The main aim of the present paper is to reveal a global nature of saddle-type blowup solutions through mathematically rigorous blow-up characterizations mentioned above with computational machineries in dynamical systems, both qualitatively and quantitatively. As any computational method inevitably suffers from numerical errors, due both to rounding and discretizing, one must question the validity of its output. This is especially through when solutions are sensitive to initial conditions, as it is the case for instance for dynamical systems possessing blow-up solutions or exhibiting chaos. In order to address the fundamental issue of reliability of computations, the recent field of computer-assisted proofs in nonlinear analysis emerged at the intersection of scientific computing, functional analysis, approximation theory, numerical analysis and topology. In essence, a computer-assisted proof is the process by which the hypotheses of a theorem are verified rigorously with the help of the computer. In the context of dynamical systems, early pioneering works include the proof of the universality of the Feigenbaum constant (Oscar 1982) and the proof of existence of the strange attractor in the Lorenz system (Tucker 2002). We refer the interested reader to the survey papers (Koch et al. 1996;Nakao 2001;Tucker 2011;van den Berg and Lessard 2015;Gómez-Serrano 2018), as well as the recent book (Nakao et al. 2019). Computerassisted proofs are one way to both characterize and visualize mathematical objects in a mathematically rigorous way. Keeping the success of computer-assisted proofs for various applications to dynamical systems (e.g., Castelli et al. 2018;D'Ambrosio et al. 2015;Matsue and Takayasu 2020a, b;Takayasu et al. 2017) in mind, we believe that studying blow-up solutions with computer-assisted proofs provides rich insights into asymptotic behavior of solutions to differential equations as well as new research directions of global dynamics and finite-time singularities.
To validate saddle-type blow-up solutions, we combine the machinery applied in preceding works, compactifications and timescale desingularizations, with the parameterization method (e.g., see Cabré et al. 2003aCabré et al. , b, 2005. The latter notion is now understood as one of universal machineries in dynamical systems, which aims at characterizing and constructing invariant manifolds, including local (un)stable manifolds of invariant sets such as equilibria and periodic orbits. Moreover, the parameterization method with rigorous ODE integrations has a great compatibility with computerassisted proofs to capture global nature of invariant manifolds in dynamical systems with their explicit enclosures. In particular, globally extended saddle-type blow-up solutions and the corresponding curves of blow-up times can be validated as easily as preceding works (Matsue and Takayasu 2020a, b;Takayasu et al. 2017).
We shall also unravel non-trivial and global nature of saddle-type blow-up solutions with the applicability of our proposed methodology through several examples. The main features of blow-up solutions we shall extract in the present paper are summarized as follows, which are not observed in preceding works or theoretical characterizations of blow-ups: • The blow-up time t max is described by a locally real analytic function of initial points under a generic assumption (Sect. 4.2). • Local foliation structure in level sets of blow-up times which is independent of dynamics at infinity is observed (Sect. 6.2). • Chain of connecting orbits including those corresponding to saddle-type blow-up solutions can separate initial points into several regions possessing significantly different properties, where solutions through these points either exist global-intime or blow up in finite time, no matter how large the magnitude of initial points is (Sect. 7.2). • The above chain of connecting orbits induces discontinuity of blow-up times (Sect. 7.2).
The first feature is one of the biggest benefits of the parameterized method in blowup studies. In preceding works, no explicit expression of local stable manifolds is obtained, yielding at most upper and lower bounds of t max (e.g., Takayasu et al. (2017)).
In the present methodology, the explicit expressions of local stable manifolds as the graphs of locally analytic functions can be applied, and hence, combined with formulae of t max by means of integrals through trajectories, we obtain the explicit formulae of t max as functions of initial points. Through computer-assisted proofs, we obtain explicit distributions of local stable manifolds with their visualizations. We then see an interesting relationship between asymptotic behavior of blow-up solutions and the corresponding t max . As the second feature, we see that the asymptotic dynamics near blow-up do not essentially contribute to determine blow-up times. In other words, only the magnitude of solutions can determine t max . The remaining features are also important and completely different from blow-up solutions possessing persistence of structure under perturbations of initial points. We see that, in the presence of saddle-type blow-up solutions, there is no relationship between the magnitude of initial points and blow-up behavior of solutions through these points. All these features rely on computer-assisted proofs, implying that all results are mathematically rigorous and the methodology toward these results are available to a large class of ODEs without any knowledge of blow-up behavior.
The rest of the present paper is organized as follows. In Sect. 2, we review a methodology for characterizing blow-up solutions from the viewpoint of dynamical systems, which is based on compactifications and timescale desingularizations studied in, e.g., Matsue (2018). Three types of compactifications are shown there: directional, Poincaré-type and parabolic-type ones. The concrete process for characterizing blowup solutions is explained for each compactification for readers' accessibility, while the fundamental idea is identical. Both advantages and disadvantages of each compactification depending the situation are finally mentioned. In Sect. 3, the parameterization method for calculating invariant manifolds is summarized. In the present paper, we restrict our attention to stable manifolds of equilibria. Under an essential assumption called the non-resonance condition of eigenvalues, local stable manifolds can be characterized as zeros of a countable family of nonlinear equations on Banach spaces. Combining with the method of radii polynomials, which is one of standard functional analytic and algebraic machineries for finding zeros of (infinite-dimensional) nonlinear maps, computer-assisted proofs of the existence and characterization of local stable manifolds are provided. Note that the non-resonance condition yields that validated stable manifolds can be given as locally real analytic functions. In Sect. 4, we provide a methodology of computer-assisted proofs of the existence of blow-up solutions. Because the detailed implementations such as the choice of compactifications and timescale transformations is problem-dependent, only the basic idea for validating blow-up solutions is presented therein. We also show that the present methodology enables us to provide an exact and explicit formula of the maximal existence time, equivalently the blow-up time, of solutions as a locally smooth or real analytic function of initial points, provided all our implementations work successfully. The present characterization of the blow-up time provides us with a quantitative feature of blowup solutions such as distributions of blow-up times depending on initial points which are not provided in preceding works (Matsue and Takayasu 2020a, b;Takayasu et al. 2017). As we shall see, the combination of compactifications with the parameterization method provide a universal concept of blow-up validations and characterizations both qualitatively and quantitatively, no matter how stable equilibria on the horizon (for desingularized vector fields) are.
The applicability of the present methodology and global nature of saddle-type blow-up solutions is shown in successive sections. In Sect. 5, a two-dimensional ODE possessing saddle-type blow-up solutions is considered. A locally defined (i.e., directional) compactification is applied, and a saddle-type blow-up solution, as well as the blow-up time as a function of initial points, is validated to check the applicability of our methodology to locally distributed blow-up solutions. In particular, the blowup profile as well as its blow-up time as a function of initial points is successfully validated, extended and visualized. In Sect. 6, we consider a three-dimensional system. The Poincaré-type compactification is applied, and one-and two-dimensional stable manifolds of saddle equilibria on the horizon are validated. The aim is to show the applicability of our methodology to saddle-type blow-up solutions distributed on multi-dimensional stable manifolds of unstable invariant sets on the horizon. Furthermore, distribution of blow-up times as functions of initial points on two-dimensional stable manifolds is validated, which shows a relationship of blow-up times to the structure of stable manifolds around the horizon. Finally, global extension of local stable manifolds is demonstrated to visualize the distribution of blow-up nature. In Sect. 7, a two-dimensional ODE which is quasi-homogeneous in an asymptotic sense is considered. The system possesses both stable and unstable equilibria on the horizon. The parabolic-type compactification is applied and a saddle-type blow-up solution is firstly validated, while validations of blow-up solutions asymptotic to stable equilibria on the horizon are already demonstrated in a preceding work (Matsue and Takayasu 2020a). The main aim of this section is to study global nature of solution families near saddle-type blow-up solutions. We see that saddle-type blow-up solutions can play the role of the separatrix decomposing initial points into collections of blow-up solutions and global-in-time solutions. In other words, saddle-type blow-up solutions can divide initial points into those with globally bounded nature and blow-up nature, no matter how large magnitudes of initial points are. This separation cannot be seen in blowup solutions induced by solutions asymptotic to stable equilibria on the horizon for desingularized vector fields. Moreover, it is also seen that blow-up times can behave in a singular manner across the saddle-type blow-up solutions. Remark that such a singular nature has not been provided only by the local theory, because the global dynamical information requires to unravel it, while many theoretical characterizations of solution structures are stated only in the local sense. We emphasize that computerassisted proofs enable us to clarify the global nature, even in dynamically singular one, with appropriately chosen machineries. All the codes for generating results with computer-assisted proofs in Sects. 5, 6 and 7 are available at Lessard et al. (2021).

Preliminary 1: Characterization of Blow-Up Solutions
In this section, we briefly review a characterization of blow-up solutions for autonomous, finite-dimensional systems of ODEs from the viewpoint of dynamical systems. In particular, we pay attention to several concrete cases which are applied in examples later, while details of the present methodology are already provided in Matsue (2018) and Matsue and Takayasu (2020a, b).
Consider the initial value problem of an autonomous system of ODEs where t ∈ [0, T ) with 0 < T ≤ ∞, f : R n → R n is a C 1 function and y 0 ∈ R n .

Asymptotically Quasi-Homogeneous Vector Fields
First of all, we review a class of vector fields in our present discussions. Matsue 2018) Let f 0 : R n → R be a smooth (i.e., C r with r ≥ 1) function. Let α 1 , . . . , α n , k ≥ 1 be natural numbers. We say that f 0 is a quasi-homogeneous function of type α = (α 1 , . . . , α n ) and order k if Next, let X = n j=1 f j (x) ∂ ∂ x j be a smooth vector field on R n . We say that X , or simply f = ( f 1 , . . . , f n ) is a quasi-homogeneous vector field of type α = (α 1 , . . . , α n ) and order k + 1 if each component f j is a quasi-homogeneous function of type α and order k + α j .
Finally, we say that X = n j=1 f j (x) ∂ ∂ x j , or simply f is an asymptotically quasihomogeneous vector field of type α = (α 1 , . . . , α n ) and order k + 1 at infinity if there is a quasi-homogeneous vector field f α,k = ( f j;α,k ) n j=1 of type α and order k + 1 such that . Throughout successive sections, consider the (autonomous) vector field (2.1), where f : R n → R n is an asymptotically quasi-homogeneous smooth vector field of type α = (α 1 , . . . , α n ) and order k + 1 at infinity.

Compactifications, Dynamics at Infinity and Blow-Up Criteria
Here we summarize the basic strategy used throughout the successive sections. The main idea is application of compactifications: the embedding of the original phase space into compact manifolds or their tangent spaces with boundaries. The boundaries then correspond to the infinity. There are mainly two different types of compactifications: the locally defined one and globally defined one. The local one is simple and applied to many preceding works involving dynamics at infinity, while the global one enables us to treat dynamics including infinity in one chart. After introducing compactifications, we derive vector fields which we mainly concern, and provide the characterization of blow-up solutions by means of dynamical systems. The concrete process for the characterization of blow-up solutions is provided for each compactification which we introduce.

A Basic Strategy
The basic strategy for characterizing blow-up solutions is summarized as follows, which is independent of the choice of compactifications introduced below.
1. For given vector field f provided by (2.1), determine its type α and order k + 1. 2. Choose an appropriate compactification of the same type α (mentioned below) as f . 3. Transform (2.1) into the corresponding one through the compactification. 4. Introduce a timescale transformation to desingularize the vector field determined by the order k +1 of f . The resulting vector field shall be called the desingularized vector field. Dynamics at infinity then makes sense through the desingularized vector field. 5. Validate hyperbolic invariant sets on the special geometric object corresponding to infinity, the horizon and their local stable manifolds for the desingularized vector field.
Once invariant sets, such as equilibria and periodic orbits, on the horizon with their hyperbolicity are validated, their local stable manifolds characterize the collection of blow-up solutions of (2.1) near blow-ups, which is the essence of our proposing methodology. In the successive parts, the blow-up characterization is shown for each compactification. An important point here is a suitable choice of "appropriate" compactifications so that our blow-up problem can be reduced to standard issues in dynamical systems.
Below are examples of such suitable compactifications, which possess both advantages and disadvantages, and hence, these compactifications have to be used according to our needs. Several characteristics of compactifications are summarized in Sect. 2.3.

Directional Compactifications
First a locally defined compactification is introduced, which shall be called a directional compactification.
Definition 2.2 (Directional compactification, cf . Dumortier et al. 2006;Matsue 2018) A directional compactification 3 of type α = (α 1 , . . . , α n ) is defined as with given direction i 0 ∈ {1, . . . , n} and the signature ±. This compactification is bijective in R n ∩ {±y i 0 > 0}, in which sense directional compactifications are local ones. In particular, this compactification is available when we are interested in trajectories of (2.1) such that the i 0 th component has the identical sign during time evolution. The image of T d is 3) The set E = {s = 0} corresponds to the infinity in the original coordinate, which shall be called the horizon.
Other geometric interpretations are mentioned in Sect. 2.3. For simplicity, fix i 0 = 1 in (2.2) in the following arguments. Next transform (2.1) via (2.2), which is straightforward: The resulting vector field is still singular near the horizon, but it turns out that the order of divergence of vector field as s → +0 is O(s −k ), and hence, the following timescale transformation is available.
where τ 0 and t 0 denote the correspondence of initial times, and s(τ d ) is the solution trajectory s(t) under the parameter τ d . We shall call (2.5) the time variable desingularization (of order k + 1).
The componentwise expression is This vector field is as smooth as f including s = 0, and hence, dynamics at infinity makes sense through dynamics generated by (2.7) around the horizon E = {s = 0}. Once the desingularized vector field (2.7) is provided, blow-up solutions can be characterized as follows.
Theorem 2.4 (Stationary blow-up: the directional version, Matsue 2018) Assume that the desingularized vector field (2.7) associated with (2.1) has an equilibrium on the horizon x * = (0, x * ) ∈ E. Also suppose that x * is hyperbolic with n s > 0 (resp. n u = n−n s ) eigenvalues of the Jacobian matrix Dg d (x * ) with negative (resp. positive) real parts. If there is a solution y(t) of (2.1) with a bounded initial point y(0) whose image x = T d (y) is on the local stable manifold W s loc (x * ; g d ), then t max < ∞ holds; namely, y(t) is a blow-up solution. Moreover, where c i is a constant with the same sign as y i (t) as t → t max .

Poincaré-Type Compactifications
The remaining compactifications we introduce here are global ones in the sense that they are embeddings of the whole phase space R n into compact manifolds with boundaries. A suitable class of global type compactifications for characterizing dynamics at infinity for asymptotically quasi-homogeneous vector fields is discussed in Matsue and Takayasu (2020a), where such a class of compactifications are called admissible global compactifications. Among such compactifications, two representative compactifications are reviewed.
Definition 2.6 (Poincaré-type compactification. cf . Matsue 2018) The Poincaré-type compactification (of type α = (α 1 , . . . , α n )) is defined as the mapping T q P : R n → R n as Its geometric interpretation is mentioned in Sect. 2.3. Note from Matsue (2018) that κ = κ q P (y) has an equivalent expression by means of x: Similar to the directional ones, for given vector field f of type α, we apply the Poincarétype compactification of the same type α. Then we have which is the alternate object off j 's in (2.4), κ = κ q P (y), and It is shown in Matsue (2018) that the above vector field is still singular on the horizon E, but the order of divergence is O(κ k ) as p(y) → +∞, equivalently p(x) → 1, which is independent of components. Therefore, a common timescale transformation can be introduced.
Definition 2.7 (Time variable desingularization: the Poincaré-type version) Define the new time variable τ d by dτ q P = κ q P (y(t)) k dt (2.14) equivalently, where τ 0 and t 0 denote the corresponding initial times, and y(τ q P ) is the solution y(t) under the timescale τ q P . We shall call (2.14) the time variable desingularization (of order k + 1).
Using this timescale, we obtaiṅ Theorem 2.8 (Stationary blow-up: the Poincaré-type version, Matsue 2018) Consider the desingularized vector field g q P associated with (2.1) given by (2.16). Assume that g q P is C 1 in a neighborhood of the horizon E and that g q P has an equilibrium on the horizon x * ∈ E. Suppose that x * is hyperbolic with n s > 0 (resp. n u = n − n s ) eigenvalues of Dg q P (x * ) with negative (resp. positive) real parts. If there is a solution y(t) of (2.1) with a bounded initial point y(0) whose image x = T q P (y) is on the local stable manifold W s loc (x * ; g q P ), then t max < ∞ holds; namely, y(t) is a blow-up solution. Moreover, where c i is a constant with the same sign as y i (t) as t → t max .

Parabolic-Type Compactifications
An alternative admissible global compactification, which shall be called the parabolictype compactification, is introduced here. Compactifications of the present type were originally introduced in Gingold (2004) and generalized in Matsue and Takayasu (2020a). Similar to the Poincaré-type compactifications, define a set D ⊂ R n by (2.11). For any x ∈ D, correspond y ∈ R n to x ∈ D by . (2.17) This equality indicates that p(y) = p(S(x)) <κ α (x) holds for all x ∈ D.
Similar to directional and the Poincaré-type ones, we apply the parabolic-type compactification of the type α which is the same as that of f to transforming (2.1). The resulting vector field is . . ,f n ) is (2.12) replacing κ by κ para , in which case Similar to the Poincaré-type case, all components of the transformed vector field are O(κ k ) as p(y) → ∞, equivalently as x approaches to E, and hence, the uniform timescale transformation can be introduced to desingularize the vector field on E.
Definition 2.11 (Time variable desingularization: the parabolic-type version) Define the new time variable τ para by where τ 0 and t 0 denote the correspondence of initial times. We shall call (2.18) the time variable desingularization (of order k + 1).
The change of coordinate and the above desingularization yield the following vector field g para , which is continuous on D = {p(x) ≤ 1}: (2.20) The desingularized vector field g para has the very similar form to g q P . On the other hand, the algebraic structure of κ is quite different from each other. In particular, κ = κ para does not include radicals in x, and hence, the smoothness of f and the asymptotic quasi-homogeneity guarantee the smoothness of the right-hand side g para of (2.20) including the horizon E. See Matsue and Takayasu (2020a) for details. This property yields a relaxation of conditions for characterizing blow-ups.
Theorem 2.12 (Stationary blow-up: the parabolic-type version, cf. Matsue 2018; Matsue and Takayasu 2020a) Consider the desingularized vector field g para associated with (2.1) given by (2.20). Assume that g para has an equilibrium on the horizon x * ∈ E. Also, suppose that x * is hyperbolic with n s > 0 (resp. n u = n − n s ) eigenvalues of Dg para (x * ) with negative (resp. positive) real parts. If there is a solution y(t) of (2.1) with a bounded initial point y(0) whose image x = T para (y) is on the local stable manifold W s loc (x * ; g para ), then t max < ∞ holds; namely, y(t) is a blow-up solution. Moreover, where c i is a constant with the same sign as y i (t) as t → t max .
The proof is essentially the same as Theorem 2.8. Indeed, only the admissible nature (discussed in Matsue and Takayasu 2020a) of T para is used to prove t max < ∞, which is the same as T q P .
The key point of our characterization of blow-ups (Theorems 2.4, 2.8 and 2.12) is that blow-up solutions for (2.1) are characterized as trajectories on (local) stable manifolds of invariant sets 6 on the horizon E for desingularized vector fields. Computations of blow-up solutions are therefore reduced to those of local stable manifolds of invariant sets, such as (hyperbolic) equilibria, for the associated vector field. Although the above theorems only characterizes the existence and local dynamical nature of blow-up solutions, combinations of our characterization with numerical computations and computer-assisted proofs provide global nature of blow-up solutions in the phase space.

Remark on Appropriate Choice of Compactifications
We have introduced three compactifications in this section. Each compactification has its own set of advantages and disadvantages, which depend on our requirements. Here we remark the choice of compactifications in case that the original vector field f is polynomial. 7 In our examples (Sects. 5, 6 and 7), all these compactifications are applied. It is worth mentioning several features of each compactification toward effective choice and applications of our machineries to practical and advanced problems.

Geometric Interpretations of Compactifications
First the geometric interpretation of each compactification is briefly summarized. Directional compactifications are not actually compactifications in the topological sense, while these are still called "compactifications" because these are inclusively discussed in the context of compactifications for applications. In fact, images of directional compactifications are interpreted as the tangent space of the Poincaré's hemisphere considered in the Poincaré-type compactifications at points on the horizon, as shown in Fig. 1a.
Global (Poincaré type and parabolic type) compactifications are geometrically simple in the homogeneous case α = (1, . . . , 1), in which case p(y) = y and we can choose β = (1, . . . , 1) and c = 1. Therefore, κ q P (y) = (1 + y 2 ) 1/2 , which is a well-known (global, but homogeneous) compactification, 8 and the resulting mapping T q R is the embedding of R n into the Poincaré hemisphere A homogeneous compactification of this kind is shown in Fig. 1b.
The geometric nature of the parabolic-type compactification with α = (1, . . . , 1) is also understood in a simple way, in which case T para is defined as See Elias and Gingold (2006) and Gingold (2004) for the homogeneous case, which is called the parabolic compactification. In particular, the parabolic compactification is the embedding of R n onto the bounded parabola in R n+1 with the focus point (x 1 , . . . , x n , x n+1 ) = (0, . . . , 0, 1). The map T para is actually defined as the composition of this embedding and the projection onto the first n components, which is shown in Fig. 1c.
Our compactificaitons introduced here are quasi-homogeneous counterparts of the above homogeneous compactifications. Geometric pictures of quasi-homogeneous Poincaré-type and parabolic-type compactifications are shown in Matsue (2018) and Matsue and Takayasu (2020a), respectively.

Remark 2.13
Any compactifications we have introduced are analytic at any point in D by using the binomial theorem in the standard calculus and the inverse function theorem for analytic mappings (e.g., Dieudonné 1960).

Directional Compactifications: Advantages and Disadvantages
A typical way to study dynamics at infinity is the application of directional compactifications introduced in Sect. 2.2.2, which is simple in the sense that the magnitude of points in the original coordinate can be measured by an independent variable s. Heuristically, associated desingularized vector fields are as complex as the original vector fields because the new variablex i in (2.2) depends only on the original variable y i and the scaling variable s. Moreover,x i is proportional to y i . Characterization of blow-up times is also simple, because they are characterized only by the asymptotic behavior of s = s(τ ). On the other hand, directional compactifications are defined only locally. Poincaré compactification with type α = (1, 1). The image T (y) of the original point y ∈ R 2 is defined as the projection of the intersection point P(M) ∈ H, given by the line segment connecting M = (y, 1) ∈ R 3 and the origin O ∈ R 3 , onto the original phase space R 2 . The horizon is identified with ∂H. The precise definition is its projection onto R 2 × {1}. c Parabolic compactification with type α = (1, 1). The image x of the original point y ∈ R 2 is defined as the projection of the intersection point P(M) ∈ H determined by the paraboloid x 2 1 + x 2 2 = x 3 in R 3 and the line segment connecting M = (y, 0) ∈ R 3 and the focus point (0, 0, 1) ∈ R 3 , onto the original phase space R 2 . The horizon is identified with the circle {x 2 1 + x 2 2 = 1} on the parabola. The precise definition is its projection onto R 2 ×{0} (Color figure online) If our interested blow-up solutions have sign-changing structure, multiple charts of compactifications can be necessary for complete descriptions of blow-up solutions. From the numerical viewpoint, change of coordinates may cause additional computation costs and errors. If one already knows from preceding mathematical or numerical arguments that targeting blow-up solutions have identical signs during time evolutions for a certain component, directional compactifications with appropriate choice of the constant sign components are efficient.

Global Compactifications: Advantages and Disadvantages
If we study blow-up solutions with sign-changing structure, or one does not have sufficient knowledge of solutions near infinity, globally defined compactifications like the Poincaré type and the parabolic type are more appropriate than directional ones, because one does not suffer from violation of integrations of differential equations due to the change of signs, or change of local charts. Because the horizon, topologically sphere-shaped boundary of the compactified space, is invariant under associated desingularized vector fields (cf. Matsue 2018), computed trajectories through points in D for desingularized vector fields are always inside D, unless unrealistic or mathematically inappropriate choice of numerical parameters. On the other hand, application of such global compactifications generally increases the degree of associated desingularized vector fields as polynomial ones, which cause complication of arguments. For example, in the case of the vector field shown in Sect. 7, we have to study (desingularized) polynomial vector fields with degree over 10, while the original one before compactification has degree at most 2 or 3. Without systematic implementations of vector fields or their derivatives like automatic differentiations, applications to concrete systems require lengthy calculations.

Poincaré Type or Parabolic Type?
Among globally defined compactifications, more than one compactifications are introduced here, the Poincaré type and the parabolic type. The simplest one in the class of admissible compactifications (e.g., Elias and Gingold 2006;Matsue and Takayasu 2020a) is the Poincaré type, which is easy to understand from geometric viewpoints and widely applied in many fields of mathematics. However, the Poincaré-type compactification has an unavoidable defect, the presence of radicals in the definition. Radicals generally lose the smoothness of desingularized vector fields on the horizon. In other words, desingularized vector fields under the Poincaré-type compactification are C 0 but not C 1 in general around the horizon. Therefore, typical "linear stability analysis" in the theory of dynamical systems does not always make sense on the horizon. Nevertheless, it should be noted that there is an exception where the Poincaré-type compactifications can be applied without losing the smoothness of resulting vector fields, which is the case if f is quasi-homogeneous (not only in the asymptotic sense), or the residual term f − f α,k has sufficiently low degree. In this case, the associated desingularized vector field is also smooth, and hence, no obstruction of C 1 smoothness on the horizon arises. Details are discussed in Matsue (2018).
Although the degree of polynomials significantly increases when we apply the parabolic-type compactifications, we do not worry about the lack of smoothness of desingularized vector fields. Indeed, parabolic-type transformations of the present type originally transforms rational functions into rational ones, unlike the Poincarétype ones (cf. Gingold 2004). We thus do not suffer from obstructions to consider dynamics at infinity when we apply parabolic-type compactifications.

The Other Choice?
The geometrically simplest compactification would be the one-point compactifications such as embedding of R n into S n , which is known as the Bendixson's compactification. One can use the Bendixson's compactification to map the infinity to a bounded point, where the corresponding dynamics possess the high degeneracy in general (e.g., Hell 2010). In order to avoid the degeneracy at infinity, we have to apply an additional desingularization (blowing up) of the infinity. The Poincaré-type and the parabolictype compactifications avoid such extra tasks for obtaining desingularized dynamics at infinity.

Preliminary 2: Parameterization Method
In this section, we introduce the theory of the parameterization method (Cabré et al. 2003a(Cabré et al. , b, 2005 to compute rigorous charts of local stable and unstable manifolds of fixed points of ODEs of the formẋ = g(x), where g is a desingularized vector field. We begin by making some assumptions, which will be sufficient for the purpose of the present paper. A1. Assume g : R n → R n is a polynomial vector field with a steady statex ∈ R n (i.e., g(x) = 0).
A2. Assume that the eigenvalues of the Jacobian matrix Dg(x) are real, nonzero and distinct. (Hence, the Jacobian matrix Dg(x) is diagonalizable over the real and x is hyperbolic.) Denote by λ 1 , . . . , λ m < 0 the stable eigenvalues of Dg(x) with ξ 1 , . . . , ξ m ∈ R n some associated stable eigenvectors. From now on, we focus on the computation of a local stable manifold, which we denote by W s loc (x), and note that dim W s loc (x) = m ≤ n. The computation of the unstable manifold is similar (e.g., see Breden et al. 2016). The idea of the computational approach is to represent the chart of the local stable manifold using a Taylor series representation P : B m → R n of the form where B m ⊂ R m is a domain (usually chosen to be a ball) on which the Taylor series converges, and where This requires making an extra assumption, which involves the notion of a resonance.
We are ready to state our third hypothesis.
Construct the following real-valued matrices: an m × m diagonal matrix with the diagonal entries made up of the stable eigenvalues and an n × m matrix whose columns are the associated eigenvectors Using the basis defined by the stable eigenvectors, the linearized equation forẋ = g(x) restricted to the stable subspace takes the forṁ The associated flow is given by e t . As indicated above our goal is to construct an analytic function P : B m → R n such that P(B m ) = W s loc (x). To obtain constraints, The dynamics in the parameter space is generated by exponentiating the matrix of stable eigenvalues . The dynamics in phase space is generated by the flow ϕ associated with the vector field g. The diagram commutes in the sense that applying first the chart map P, and then, nonlinear flow ϕ(t, ·) is required to be the same as applying the linear dynamics e t , and then, the chart map P. The result is that the dynamics on the local stable manifold are described by the stable linear dynamics (Color figure online) so that we can solve for P, we begin by insisting that P be a conjugacy between the flow ϕ ofẋ = g(x) restricted to W s loc (x) and the flow e t of the linear equation. The most obvious restriction is that P must map fixed points to fixed points, and hence, To obtain the conjugacy, we assume that for all θ ∈ B m . The geometric meaning of this conjugacy is illustrated in Fig. 2 because the entries of are negative.
Note that any function P(θ ) satisfying Eq. (3.4) is one-to-one on B m . To see this observe that P is tangent to the stable eigenspace at the origin as D P(0) = A 0 . Moreover, recall that A 0 is of full rank as its columns are linearly independent. By the implicit function theorem, P is of rank m, and hence one-to-one, in some neighborhood U ⊂ B m of 0. Now suppose that θ 1 , θ 2 ∈ B m and that P(θ 1 ) = P(θ 2 ). Then for any t ∈ R, ϕ (t, P(θ 1 )) = ϕ (t, P(θ 2 )) by the uniqueness of the initial value problem. Choose T > 0 large enough so that e T θ 1 , e T θ 2 ∈ U . By the conjugacy relation, we have that P e t θ 1 = P e t θ 2 , and because the arguments are in U , the local immersion gives that e T θ 1 = e T θ 2 . But e T is an isomorphism and we have θ 1 = θ 2 . We therefore conclude from the above discussion that P(B m ) = W s loc (x). The utility of (3.4) is limited by the appearance of the flow ϕ in the equation. In practice, the flow is only known implicitly, that is, it is determined by solving the differential equation. The following lemma establishes a more practical infinitesimal version of (3.4).

Then P(θ ) satisfies the conjugacy relationship (3.4) if and only if P is a solution of the partial differential equation (PDE)
Proof Let P : B m → R n be a smooth function with P(0) =x and D P(0) = A 0 .
(⇐ ) Suppose that P(θ ) solves the partial differential equation (3.6) in B m . Choose a fixed θ ∈ B m and fix t > 0. Define the function γ : [0, t] → R n by (3.7) Then, γ (0) = P(θ ) and where we pass from the first to the second equality by the chain rule, from the second to the third equality by the invariance Eq. (3.6) and the fact that e t θ ∈ B m when t > 0, and from the third to the fourth equation by the definition of γ . Hence, γ is the solution of the initial value problem Therefore, by definition ϕ(t, γ (0)) = γ (t), and it follows from (3.7) and (3.8) that ϕ(t, P(θ )) = P(e t θ).
( ⇒) Suppose that P satisfies the conjugacy relationship (3.4) for all θ ∈ B m . Fix θ ∈ B m and differentiate both sides with respect to t in order to obtain g(ϕ(t, P(θ ))) = D P(e t θ) e t θ.
Taking the limit as t → 0 gives that P(θ ) is a solution of (3.6).
As a consequence of Lemma 3.2, it should now be clear that computing a local m-dimensional stable manifold is equivalent to find a solution P : B m → R n of the PDE (3.6). As mentioned earlier, the idea is to use a Taylor series representation of the form (3.1). Note that since g : R n → R n is a polynomial vector field, the power series expansion of g(P(θ )) involves Cauchy products. Denote the Taylor expansion of g(P(θ )) as where we abuse slightly the notation and used the same notation g(a) to denote the vector field g where the monomial terms in the variables x 1 , . . . , x n are replaced by Cauchy products in the variables a 1 , . . . , a n .
Formally plugging the Taylor expansion (3.1) in the PDE (3.6) results in where α · λ def = α 1 λ 1 + · · · + α m λ m and a α = ((a 1 ) α , . . . , (a n ) α ) ∈ R n . The first-order constraints (3.5) imply that where e j is the jth vector of the canonical basis of R n . In other words, the Taylor coefficients a α for |α| ∈ {0, 1} are fixed and do not need to be solved for.
Denote the Banach space Moreover, denoting the Banach spacẽ and X def = (˜ 1 ) n , we get that F : X → X . Slightly generalized Banach spaces can be also applied, which will be mentioned in Sect. 5.
Denote by B m . . , m} the unit polydisc in C m . We have the following result.
Theorem 3.4 Assume that Assumptions A1, A2 and A3 are satisfied. If there exists a ∈ X such that F(ã) = 0 with F given in (3.9), then the corresponding Taylor expansion P : B m 1 → R n given by provides a parameterization of a local stable manifold ofx, that is, P(B m 1 ) = W s loc (x). Proof Assume thatã ∈ X solves F(ã) = 0. Then by construction, the function P(θ ) given in (3.13) converges absolutely and uniformly on B m 1 as for each j ∈ {1, . . . , n} since ã 1 < ∞. By construction, the function P : B m 1 → R n given in (3.13) satisfies the first-order constraints (3.5) and the PDE (3.6). By Lemma 3.2, P satisfies the conjugacy relationship (3.4). Finally, we conclude that P : B m 1 → R n provides a parameterization of a local stable manifold ofx, that is, The strategy to compute a parameterization of W s loc (x) is now clear. Fix the lengths of the eigenvectors ξ 1 , . . . , ξ m such that we can computeã ∈ X such that F(ã) = 0. This is achieved with a Newton-Kantorovich type argument, which we now state.
the closed ball of radius r > 0 centered at a given b ∈ X and B(X 1 , X 2 ) the space of bounded linear operators between two Banach spaces X 1 and X 2 .
Theorem 3.5 (A Newton-Kantorovich type theorem) Let X and X be Banach spaces, where · B(X ) denotes the operator norm. Define the radii polynomial by If there exists r 0 > 0 such that p(r 0 ) < 0, then there exists a uniqueã ∈ B r 0 (ā) such that F(ã) = 0.
The strategy of Theorem 3.5 requires obtainingā (a numerical approximation), the operator A † ∈ B(X , X ) (an approximation of the Fréchet derivative D F(ā)) and the operator A ∈ B(X , X ) (an approximate inverse of D F(ā)).
To compute the numerical approximationā, we first consider a finite-dimensional projection of the map F : X → X . Fixing a dimensional Taylor projection number N , denote by X (N ) the finite-dimensional space . . , N }} the number of multiindices α with order between 2 and N . Given a vector b = (b ) | |≥0 ∈ 1 , consider the projection We generalize that projection to get N : Given a ∈ X , we denote Moreover, we define the natural inclusion ι N : Assume that a numerical approximationā (N ) = ā (N ) ) ≈ 0 has been computed (e.g., using Newton's method). Given j = 1, . . . , n, denoteā j = ι Nā (N ) j ∈ 1 and denoteā = (ā 1 , . . . ,ā n ), and for the sake of simplicity of the presentation, we use the same notationā to denoteā ∈ X and a (N ) ∈ X (N ) . Denote by D F (N ) (ā) the Jacobian of F (N ) atā, and let us write it as The next step is to construct the linear operator A † (an approximate derivative of the derivative D F(ā)) and the linear operator A (an approximate inverse of D F(ā)). Let This allows defining the linear operator A as Having obtained an approximate solutionā and the linear operators A † and A, the next step is to construct the bounds Y 0 , Z 0 , Z 1 and Z 2 (r ) satisfying (3.14), (3.15), (3.16) and (3.17), respectively.

The Y 0 Bound
Denote by d the highest order nonlinear term of the vector field f . Then sinceā consists of Taylor coefficients of order N , then (F(ā)) α = 0 for all |α| > d N .
which is a collection of finite sums that can be evaluated with interval arithmetic. We conclude that (3.22)

The Z 0 Bound
We look for a bound of the form Recalling the definitions of A and A † given in (3.21) and (3.20) We remark that (B i, j ) n 1 ,n 2 = 0 for any i, j ∈ {1, . . . , n}, whenever n 1 > N or n 2 > N . Hence, we can compute the norms B i, j B( 1 ) using the following standard result.
where each norm B i, j B( 1 ) can be computed using formula (3.23) with vanishing tail terms.

The Z 2 Bound
For a fixed r * > 0, set which satisfies (by the mean value inequality in Banach spaces) Evaluating the bound (3.26) is straightforward with interval arithmetic and the easily computed formulas of the second derivatives of each component f i of the vector field f .

Rigorous Enclosure of the Points on W s loc (x)
Assume that assumptions A1, A2 and A3 are satisfied for a fixed pointx. Let λ 1 , . . . , λ m < 0 be the corresponding non-resonant real (stable) eigenvalues and ξ 1 , . . . , ξ m ∈ R n be some associated stable eigenvectors.
Consider N ≥ 2 the order of the Taylor approximation, and as before, assume that a numerical approximationā (N ) (N ) ) ≈ 0 has been computed. Denote by (3.27) Using a computer program in MATLAB using the interval arithmetic package INT-LAB, we can compute rigorously the bounds Y 0 , Z 0 , Z 1 and Z 2 satisfying (3.22), (3.24), (3.25) and (3.26), respectively. Define the radii polynomial p(r ) defined in (3.18), and assume the existence of r 0 > 0 such that p(r 0 ) < 0. From Theorem 3.5, there exists a uniqueã ∈ B r 0 (ā) such that F(ã) = 0. By Theorem 3.4, the corresponding Taylor expansion P : B m 1 → R n given by (3.13) provides a parameterization of a local stable manifold ofx, that is, . From the computer-assisted proof, we immediately obtain a rigorous upper bound for the C 0 error bound between the approximate parameterization (3.27) and the true parameterization. More explicitly, for a fixed j = 1, . . . , n Using that estimate, given a point z ∈ B m 1 in parameter space, one may evaluate rigorously the corresponding value P(z) ∈ W s loc (x) on the local stable manifold using the following enclosure where P (N ) j (z) can be computed with interval arithmetic using the formula (3.27).

Saddle-Type Blow-Up Solutions: Basic Methodology for Validations and Extensions
In this section, we provide a methodology for validating saddle-type blow-up solutions with computer-assisted proofs. A remarkable feature obtained from theorems mentioned in Sect. 2 is that stability of equilibria on the horizon does not matter for characterizing blow-up solutions. Therefore, we can characterize blow-up solutions whose blow-up direction is characterized by unstable equilibria 9 in the same way as stable ones. When we emphasize structure of equilibria on the horizon, we shall call them as follows.
Definition 4.1 We say that a blow-up solution is sink-type (resp. saddle-type) if it is transformed into a trajectory on W s loc (x * ; g) with a sink (resp. saddle) equilibrium x * on the horizon for the associated desingularized vector field g introduced in Sect. 2.
There are many studies of blow-up solutions through analytic arguments (e.g., Fila and Matano 2002;Herrero and Velázquez 1997;Winkler 2013) or numerical simulations (e.g., Anada et al. 2017;Cho 2016;Cho et al. 2007;Zhou and Saito 2017), many of which would be sink-type through related numerical simulations and computer-assisted proofs (e.g., Matsue 2018; Matsue and Takayasu 2020a, b). Whereas, saddle-type blow-up solutions are quite difficult to calculate and to understand the role in global dynamics, because generic small perturbations of initial points [for (2.1)] break the structure. Even in the context of dynamics at infinity (i.e., without concerning the blow-up nature of solutions), there are very limited studies for characterizing trajectories asymptotic to the horizon themselves and their global nature, except special cases such as planar dynamical systems (e.g., Dumortier and Herssens 1999;Dumortier et al. 2006). On the other hand, saddle-type blow-up solutions themselves can exist in various types of differential equations, many of which do not concern with its sensitivity under perturbations of initial points, but are interested only in their existence and/or persistence of blow-up structure under perturbation of initial points is mentioned implicitly (cf. Harada 2016Harada , 2017Nouaili and Zaag 2015 for complex-valued PDEs).
Here we will see that our validation methodology provide not only a systematic way to capture saddle-type blow-up solutions but also distributions of a collection of blowup solutions in the phase space, both of which are with mathematical rigor. Moreover, as seen in preceding works, methodologies with computer-assisted proofs provide explicit enclosures of computation objects. This property enables us to visualize the distribution of solution profiles and blow-up times depending on initial points of blowup solutions. As a by-product of the application of parameterization method reviewed in Sect. 3, we obtain an explicit formula of t max as a function of initial points and its smoothness.

Basic Methodology
First we discuss a basic methodology for validating (locally defined) blow-up solutions and their extension. The fundamental steps consist of the following: 1. Validation of local stable manifolds of equilibria on the horizon for desingularized vector fields; 2. Extension of validated stable manifolds via rigorous integration of desingularized vector fields. These steps are shown to provide a collection of divergent solutions of (2.1), according to Theorems 2.4, 2.8, and 2.12 except the evaluation of t max . When an equilibrium p on the horizon is stable, the validation procedures reported in Matsue and Takayasu (2020a, b); Takayasu et al. (2017) allow (a) studying the local stable manifold of p by means of locally defined Lyapunov functions and (b) computing rigorous enclosure of solutions converging to p, hence yielding a rigorous bound of the blow-up time. Although the same strategy or similar topological arguments such as covering relations (e.g., Zgliczyński 2009) can work effectively, we apply the parameterization method to validating local stable manifolds of equilibria on the horizon for desingularized vector fields here instead.
An important merit of the parameterization method is that only the stable information of equilibria can be treated through the whole computations involving invariant manifolds, no matter how unstable equilibria or general invariant sets are. In other words, if we can compute stable eigenvectors at the equilibria and a topological conjugacy P with high accuracy, we obtain the local stable manifold without containing intrinsic unstable information of equilibria. 10 Moreover, we obtain the embedding of the parameterized invariant manifolds, and hence, the distribution of locally defined invariant manifolds in the whole phase space can be captured at the same time. This distribution greatly helps us with investigating the behavior of trajectories far from invariant manifolds. Universality of such features in the parameterization is shown in many preceding works (e.g., Barker et al. 2020;Breden et al. 2016;Gonzalez and Mireles-James 2017;Mireles-James 2018;van den Berg et al. 2011) for obtaining global nature of dynamical systems.
After validating the locally parameterized stable manifolds, these manifolds can be extended through time integrations of the time-reversed desingularized vector fields. In particular, we obtain globalized stable manifolds whose preimages under compactifications are (candidates of) families of saddle-type blow-up solutions. 11 Globalization of invariant manifolds enables us to investigate global nature of dynamical systems, including blow-up solutions in the present study, while the methodology itself is standard and essentially identical with the one used in preceding works (e.g., Matsue and Takayasu 2020a, b;Takayasu et al. 2017).

Remark 4.2 (Rigorous integrators of ODEs)
Many methods for rigorously integrating solution trajectories of vector fields have been proposed over the last thirty years. 10 In topological arguments such as local Lyapunov functions and covering relations, topological information of both stable and unstable directions around equilibria are necessary to validate locally defined invariant manifolds, which cause a big difference of treatments between stable and unstable invariant sets. 11 Needless to say, the proposing methodology can be also applied to sink-type blow-up solutions.
The most famous achievement is the resolution of Smale's 14th problem by Tucker (2002). We refer to Berz and Makino (1998); Bünger (2020); Immler (2018); Kashiwagi and Oishi (1994); Lessard and Reinhardt (2014); Lohner (1987); Zgliczynski (2002) for different methods for rigorous integration of ODEs. These methods are based on fixed-point arguments, which is equivalent to show the existence of solution trajectories, and several techniques of interval arithmetic. For the sake of forward time integration, we use a C++ Library for rigorous integration of ODEs, which is named the kv library (Kashiwagi). This integrator is based on an interval representation of the solutions' Taylor series and the Affine arithmetic (Rump and Kashiwagi 2015), which is a technique for preventing the so-called wrapping effect in interval analysis.
The remaining issue is the finiteness of t max and its explicit enclosure to assure that our validated trajectories indeed correspond to blow-up solutions for the original system. The blow-up time t max generally depends on initial points of solutions. We now introduce an explicit estimate methodology for obtaining blow-up times.
The blow-up time t max is defined by the improper integral as τ → ∞ in (2.6), (2.15) or (2.19), where τ is the corresponding timescale. The basic approach to enclose t max is to divide the integral into two parts: for someτ > τ 0 , where h is a functional representing integrands for characterizing t max depending on solutions of the desingularized vector field g. The key point of the successive treatments is to enclose t max,2 (τ ) by using the asymptotic information of trajectories, say the fact that trajectories of our interest are located on a local stable manifold W s loc ( p; g) of a saddle equilibrium p. Once the manifold W s loc ( p; g) is constructed through the parameterization P, the functional h is expressed by means of a (nonlinear) combination of P, and the enclosure of t max,2 (τ ) is also computed through P itself or its enclosure. As for t max,1 (τ ), we directly enclose the integral through the enclosed trajectories via ODE integrations.
When we extend the local stable manifold, then integrate the vector field g in the reverse time direction and evaluate t max,1 (τ ) by The functional h is given as follows, depending on the choice of compactifications:

Remark 4.3
The absence of constant terms in the integrand of t max is the most essential property to show that t max < ∞ in the preceding work Matsue (2018) when equilibria on the horizon are hyperbolic, where the Hartman-Grobman-type argument is applied to extracting the exponentially decaying property of the integrand. This property is essentially independent of the choice of compactifications associated with appropriately chosen timescale desingularizations. The present argument explicitly extracts this property to verify t max < ∞ by means of the parameterization method.
Summarizing the above arguments, our methodology for validating (saddle-type) blow-up solutions consists of the following.
1. Validate the local stable manifold of an equilibrium on the horizon for desingularized vector fields via the parameterization method. 2. Extend the validated stable manifold via (backward) integration, which is done by considering the time-reversed desingularized vector fields. 3. Compute a rigorous enclosure of the blow-up time t max through the decomposition of the form (4.1) as well as direct integrations through trajectories and parameterizations.
In the subsequent sections, applicability of the present methodology is shown. In particular, we aim at showing the following features, respectively: • Section 5 shows an application of directional compactifications for validating saddle-type blow-up profiles and computing the validated curve t max as a function of initial points. • Section 6 shows an application to higher-dimensional systems. In the present study, we consider an artificial three-dimensional system. The Poincaré-type compactification is applied to an asymptotically homogeneous vector field. This example shows the global phase portrait involving multiple saddle-type blow-up solutions. • Section 7 shows a characteristic nature of saddle-type blow-up solutions with bounded global solutions which separate the whole phase space into four sets, one of which is the set of points such that solutions through them determine time global solutions for both time directions and the others are the sets of points such that solutions through them are blow-up solutions in positive and/or negative time directions. Dependence of t max as a function of initial points including saddletype blow-up solutions is also addressed. The parabolic-type compactification is applied to an asymptotically quasi-homogeneous vector field.

Smooth Dependence of t max on Initial Points
Explicit expressions of t max shown in Sect. 4.1 indicate that t max depends continuously, possibly smoothly, on initial points within stable manifolds of hyperbolic equilibria (for desingularized vector fields) on the horizon, which is just a consequence of standard calculus. One of benefits of applying the parameterization method reviewed in Sect. 3 is that t max can be treated as a locally analytic function on initial points of solutions. Here we discuss the dependence of t max on initial points in more details. Consider the desingularized vector field g associated with the directional (resp. Poincaré type and parabolic type) compactification with the associated timescale desingularization. Let p * ∈ E be a hyperbolic equilibrium for g. First note that typical choices of timescale desingularizations h, namely the integrand of t max mentioned in (4.1), satisfy the following properties (so that trajectories for g is orbitally equivalent to the original dynamical system (cf. Matsue 2019)): • It vanishes at p * .
• It is positive along W s loc ( p * ; g). • It is smooth, in particular analytic, except k/2c / ∈ N in the case of the Poincaré-type compactifications.
Assume that all assumptions in Theorem 3.5 for F given in (3.9) as well as (A1), (A2) and (A3) associated with the hyperbolic equilibrium p * for g are satisfied, in which case the local stable manifold W s loc ( p * ; g) is parameterized by an analytic function P defined on the m-dimensional unit polydisc B m 1 so that W s loc ( p * ; g) = P(B m 1 ). For typical asymptotically quasi-homogeneous fields, the function h can be chosen as a polynomial or a rational function whose denominator is polynomial and positive on W s loc ( p * ; g). From these observations, we obtain the following proposition.
is an analytic function on B m 1 satisfying U (0) = 0. Proof Because h and P are analytic, then so is h • P, and hence, the integrand of U is written by the convergent series h • P(e τ θ) = |α|≥0 c α e τ θ α .
Denoting α · λ = m i=1 α i λ i , the assumption h • P(0) = 0 implies that c 0 = 0 and which converges uniformly on B m 1 . Indeed, letting σ gap def = min j=1,...,m |λ j | > 0, then (c j ) α θ α holds for j = 1, . . . , n, uniformly in B m 1 . This proposition provides a fundamental feature of blow-up times. Combining with the choice of functionals h mentioned in Sect. 4.1, the compositions of h and the parameterization P are given as follows for each compactification: and x(τ ) = P(e τ θ).
The function U (θ ) in (4.2) equals to t max = t max (θ ) under corresponding compactifications and timescale desingularizations. Note that the function h corresponds to the timescale transformation factor. Different choice of timescale desingularizations provides different h and, consequently, different determination of U (θ ) = t max (θ ).
The inverse T −1 of compactifications away from the horizon can be described by analytic functions because it is defined by the n-tuples of composite functions of radicals and rational functions whose singularities in the sense of the loss of regularity and convergence of infinite series are located on the horizon. The analytic dependence of t max on bounded initial points for (2.1) near blow-up is therefore inherited by restricting our attention to stable manifolds of equilibria on the horizon for desingularized vector fields.
Theorem 4.5 (Analyticity of blow-up times) Let t max be given by the directional (resp. the Poincaré type with k/2c ∈ N, or the parabolic type) compactification given by (2.6) (resp. (2.15) and (2.19)). Let W s loc ( p * ; g) be a local stable manifold of a hyperbolic saddle p * on the horizon for the desingularized vector field g = g d (resp. g = g q P and g = g para ) given by the parameterization P satisfying all requirements presented in Proposition 4.4. Let y 0 ∈ R n be a point such that the solution y(t) to (2.1) with y(t 0 ) = y 0 is mapped into the trajectory {(T (y))(τ )} τ ≥τ 0 included in W s loc ( p * ; g) through T = T d (resp. T = T q P and T = T para ) and the corresponding timescale desingularization.
Then the blow-up time t max is real analytic at y 0 in int T −1 (D) T −1 (W s loc ( p * ; g)), where D is given by (2.3) (resp. (2.11) for Poincaré and parabolic types). Moreover, t max = t max (y 0 ) converges to 0 as y 0 goes to infinity along the solution y(t).
Proof Let y 0 ∈ int T −1 (D) T −1 (W s loc ( p * ; g)) be arbitrary. Then there is a unique point θ 0 ∈ B m 1 such that y 0 = T −1 (P(θ 0 )). We shall write θ 0 = P −1 (T (y 0 )), where the expression of P −1 reflects the one-to-one property of P on B m 1 . Because T is analytic in D (Remark 2.13), then so is P −1 • T in int T −1 (D) T −1 (W s loc ( p * ; g)). 12 Then the blow-up time t max = t max (y 0 ) at y 0 is written by (T (y 0 ))))dτ ≡ t max (y 0 ), 12 Analyticity of P −1 follows from that of P by assumption, linear isomorphism property of D P and the inverse function theorem for analytic functions. See, e.g., Dieudonné (1960) for the latter argument.
where h is a function mentioned just after the proof of Proposition 4.4. The integrand is analytic in int T −1 (D) T −1 (W s loc ( p * ; g)), according to the same argument as the proof of Proposition 4.4. Therefore, t max = t max (y 0 ) is analytic at y 0 ∈ int T −1 (D) T −1 (W s loc ( p * ; g)). Our assumption for the solution y(t) implies that the property of y(t) going to infinity as t → t max corresponds to (T (y))(τ ) → p * as τ → ∞. Moreover, for any y 0 ∈ int T −1 (D) T −1 (W s loc ( p * ; g)), we can choose a pointỹ 0 ∈ R n such that T (ỹ 0 ) ∈ W s loc ( p * ; g) and that y 0 = y(t) = y(t;ỹ 0 , t 0 ) for some t > t 0 , where t is uniquely determined byỹ 0 . The last assertion is equivalent to T (y 0 ) = (T (ỹ 0 ))(τ ) forτ > τ 0 uniquely determined by t and the timescale desingularization. Fix the pointỹ 0 . Consider the decomposition (4.1) of t max = t max (ỹ 0 ) with the integrand S h (τ ) ≡ h • P(e τ (P −1 (T (y 0 )))). Notice that the second term in the right-hand side is the contribution of y 0 to the determination of t max . As a result, we have As mentioned, the convergence of T (y 0 ) to p * corresponds toτ → ∞. Because S h is analytic in B m 1 , the integral t max (y 0 ) goes to 0 asτ → ∞. This implies the final statement in the theorem.
Theorem 4.5 indicates that t max depends analytically on initial points on T −1 (W s loc ( p * ; g)), provided that the non-resonance condition holds for eigenvalues of Dg( p * ). Furthermore, a computer-assisted proof for the existence of P as discussed in Sect. 3 provides the explicit region where the analyticity of t max as a function of initial points of trajectories is guaranteed. Extending W s loc ( p * ; g) through the flow and using the smooth dependence of the flow on initial points, we can extend t max as a smooth function of the initial points whose smoothness depends on that for the flow, as long as W s loc ( p * ; g) is smoothly continued. In particular, t max can be analytically continued if the vector field g is analytic. Note that the analyticity or even continuity of t max is not guaranteed as a function of y in R n because the expression of t max as an analytic function only makes sense on W s loc ( p * ; g). The different choice of p * induces a different expression of P and hence of t max .
In the end of this section, we shall derive a detailed implementation of t max for directional compactifications, namely an estimate of the integral ∞ τ s(τ ) k dτ given in (2.6). The corresponding calculations of t max for Poincaré type with k/2c ∈ N and parabolic-type compactifications are achieved in the same manner. Let p * ∈ E be a hyperbolic equilibrium for the desingularized vector field g = g d . Assume that the parameterization method around p * works and the local stable manifold W s loc ( p * ; g) is obtained through the (m-dimensional) stable polydisk B m 1 and the parameterization P. For simplicity, stable eigenvalues {λ i } m i=1 of the linearized matrix of the desingularized vector field at p * are assumed to be simple and real. In particular, λ i < 0 for i = 1, . . . , m. Recalling (3.3), write Then the solution (s(τ ),x(τ )) ∈ W s loc ( p * ; g) is written by 3) is the parameterization of W s loc ( p * ; g). The rightmost integral in (2.6) can be calculated as follows, once we obtain a concrete form of P: Here we observe that (a 1 ) 0 = 0, since P(0) = p * is the equilibrium for the desingularized vector field (2.7) and P 1 (0) = 0 from our choice of compactifications. Using the previous notation and the above fact, the above integral is formally written as follows:

Denote the Cauchy product over multi-indices by
(4.5) In particular, the denominator α · λ is strictly negative for all possible α, and the analyticity of P implies that the above infinite sum is convergent uniformly in B m 1 . The final formula (4.5) implies that we can calculate the rigorous value of t max near blow-up, once we obtain the parameterization of the local stable manifold W s loc ( p * ; g) and fix the point θ ∈ B m , namely P(θ ) ∈ W s loc ( p * ; g). As seen below, the similar expressions of t max to (4.5) can be obtained for Poincaré-type and parabolic-type compactifications.
Remark 4.6 (Special case) If n s = 1, the explicit expression (4.5) admits the simpler form: Indeed, α becomes a single index l and where we have used the fact that (a 1 ) 0 = 0 and that the Cauchy product (a k 1 ) l = ( k times a 1 * · · · * a 1 ) l , with l < k contains at least one (a 1 ) 0 .

Remark 4.7 (Integrands and smoothness of t max )
The concrete procedure to compute the integral (4.5) or its upper bound depends on problems, namely the choice of compactifications and timescale desingularizations.
• Our first example (Sect. 5) applies a directional compactification, while the timescale desingularization has the different form from (2.5) so that the resulting desingularized vector field is polynomial. Instead, t max requires integrations of rational-type functions. Nevertheless, the essence of the above argument, namely the absence of constant terms in the integrand of t max , can be applied to verifying that t max < ∞. Analyticity of the integrand follows from that for both the numerator and the denominator with additional boundedness property of the denominator. Detailed derivation of t max or its upper bound is shown in subsequent sections. • In the case of Poincaré-type compactifications, analyticity of t max is not guaranteed when k/2c / ∈ N, because the function h(x) = x k/2c is not analytic at x = 0. This failure comes from the "mismatch" of properties of vector fields in the sense that the order k + 1 and the type α, consequently the natural number c, determining an appropriate Poincaré-type compactifications are determined by the asymptotic quasi-homogeneity of vector fields. We then need further estimates for calculating t max in such a case. The difficulty originated from this issue can be overcome by choosing the parabolic-type compactifications.

Remark 4.8 (Lyapunov functions versus parameterizations for expressing t max )
In the preceding studies (e.g., Matsue and Takayasu 2020a, b;Takayasu et al. 2017), t max in all examples there are enclosed by means of Lyapunov functions. Local Lyapunov functions only provide upper bounds of t max , because they do not trace concrete trajectories on stable manifolds, but values of functionals on trajectories, implying that smoothness arguments for t max as a function of initial points cannot be derived. Instead, simple inequalities by means of Lyapunov functions provide upper bounds of t max even in the case of Poincaré-type compactifications with k/2c / ∈ N, as demonstrated in Takayasu et al. (2017). Moreover, non-resonance condition (A3) is not required for estimations.
On the other hand, we can trace trajectories on stable manifolds by means of parameterizations, indicating that t max is "exactly" calculated through the integration of given functions depending on solutions. In particular, we can explicitly discuss properties of t max as a functions of initial points. In compensation for these precise information, however, we have to take care of analytic information of dynamical systems to ensure smoothness or analyticity of functions of interests, such as non-resonance condition (A3) for analyticity of P providing the conjugacy to linearizations, matching of integers k and c for Poincaré-type compactifications mentioned in Remark 4.7.

Example 1: Validation and Visualization of Globally Extended Saddle-Type Blow-Ups
In what follows, we show several applications of our proposed methodology not only to show its applicability but also to reveal several remarkable features of saddle-type blow-up solutions. The first problem is concerned with saddle-type blow-up solutions for the following system: and ρ 2 > ρ 1 are positive constants. Moreover, Points (β L , v L ) and (β R , v R ) are given in advance.
The system (5.1) is the reduced problem of (5.4) satisfying viscosity profile criterion, namely the traveling wave problem with respect to the frame coordinate ζ = x − ct with the boundary condition where c is the speed of traveling waves. Saddle-type blow-up solutions for (5.1) are considered as components of singular shock wave solutions 13 to (5.4).

A Local One-Dimensional Stable Manifold of p 2 in (5.7)
Consider the system of desingularized ODEṡ .
We focus on the one-dimensional stable manifold of the steady state x (2) with stable eigenvalue λ def = − ρ 2 2 (ρ 2 − ρ 1 ) < 0 and corresponding stable eigenvector Our goal is to produce an analytic function P : (−ν, ν) → R 2 with ν = 1 that parameterizes W s loc (x (2) ). Note that the parameter ν can generalize the setting in Sect. 3, and corresponds to that appeared in Sect. 5.1.3. The Taylor series representation has the form P(θ ) = ∞ n=0 a n θ n where a n = (a 1 ) n (a 2 ) n .

A Computer-Assisted Proof
Fixing N = 300, we computed the bounds Y 0 , Z 0 , Z 1 and Z 2 as presented in Sects. 3.1, 3.2, 3.3 and 3.4, respectively. Then, we applied Theorem 3.5 to proving the existence ofã ∈ B r (ā) such that F 1 (ã) = F 2 (ã) = 0 with F 1 and F 2 given in (5.10) and (5.11), respectively. More explicitly, we got that ã −ā X ≤ r = 4.2 × 10 −13 . The Taylor series representation of the parameterization of the local stable manifold has the form and denote by n θ n whereā n = (ā 1 ) n (ā 2 ) n the numerical approximation of the local stable manifold. Then,

Computing the Blow-Up Time
Given a point (x 1 (0), s(0)) ∈ W s loc ( p 2 ) (with p 2 = (2, 0)), the blow-up time is given by (5.12) Given that (x 1 (0), s(0)) = (P 1 (θ ), P 2 (θ )) for a given θ ∈ (−ν, ν), we get from (3.4) that ϕ (t, P(θ )) = P e λt θ for all t ≥ 0. Hence, the solution (x 1 (t), s(t)) with the initial point (x 1 (0), s(0)) = (P 1 (θ ), P 2 (θ )) is given by (x 1 (t), s(t)) = P e λt θ . Rescaling the time interval η ∈ [0, ∞] to u ∈ [θ, 0] leads (via the change of coordinates u = e λη θ ) to (5.13) Now, note that  Note that the plotted local stable manifold is defined for the desingularized vector field (5.8), which itself makes sense for both positive and negative x 2 . On the other hand, this makes sense only in {x 2 > 0} as the corresponding object to the original vector field (5.1), while the horizon {x 2 = 0} corresponds to the infinity in the original (β, v)-phase space Assume now that we have (again using rigorous numerics) obtained which is in essence computable (that is, we can provide a numerical approximation together with rigorous error bounds). In Fig. 3, we present a rigorous numerical computation (with rigorous bounds) of the value of t max as a function of θ , that is, as a function of the initial points P(θ ) on W s loc ( p 2 ). The rigorous error bound is obtained by computing rigorously the Taylor coefficients of r n in the expansion (5.14). We present how to do that next.

Rigorous Computation of the Coefficients r n
Throughã = (ã 1 ,ã 2 ) with ã −ā X ≤ r = 4.2 × 10 −13 the power series P i (u) = n≥0 (ã i ) n u n are determined. The goal in this section is to compute rigorously the coefficients r n of R(u) = n≥0 r n u n such that [P 1 (u)] 2 R(u) = Q(u). This amounts to solve the Taylor coefficients equation Using Newton's method, assume that we computedr such that ψ(r ) ≈ 0. Denote by Dψ (N ) (r ) the Jacobian of ψ (N ) atr . The next step is to construct the linear operator A † (an approximate derivative of the derivative Dψ(r )) and the linear operator A (an approximate inverse of Dψ(r )). Let A † be defined as This allows defining the linear operator A whose action on an element h ∈ 1 ν is where, given ν ≥ 1, In this expression, X = ( 1 ν ) 2 is applied and the norm is given by a X := max ( a 1 ν , a 2 ν ) for a = (a 1 , a 2 ) ∈ X . We have used a generalized setting of Banach spaces towards further applications. Having obtained an approximate solution r and the linear operators A † and A, the next step is to construct the bounds Y 0 , Z 0 , Z 1 and Z 2 (r ) satisfying (3.14), (3.15), (3.16) and (3.17), respectively. Note that since problem (5.15) is linear, then Z 2 = 0. The bound Y 0 We look for a bound such that Aψ(r ) ν ≤ Y 0 . Expand where we used that The bound Z 0 It is the same computation as the one presented in Sect. 3.2. The bound Z 1 Given h ∈ 1 ν , denote Define β k = (ã 2 1 ) k for k > 0 and β 0 = 0. Hence, We therefore set |(ā 2 1 ) n |ν n + 2 ā 1 ν r 0 + r 2 0 .
Assume that using the radii polynomial approach of Theorem 3.5, we prove the existencer ∈ B r min (r ) such that ψ(r ) = 0. Hence, given θ ∈ (−ν, ν), t max given in (5.14) can be controlled which can be evaluated rigorously with interval arithmetic.

Remark 5.3
The above estimate directly shows the analyticity of t max on θ , which is implicitly guaranteed by analyticity of the parameterization P and the uniform boundedness of the denominator x 1 (η) = P 1 (u) away from 0 on W s loc ( p 2 ). See Fig. 4 about the latter fact. Finally, we have applied X = ( 1 ν ) 2 with ν = 1 in the present validations. Different choice of ν can be also valid.

Extension of the Stable Manifold of p 2 in (5.7) and Blow-Up Time Validations
Once we validate the local stable manifold of a saddle equilibrium, we can extend the manifold integrating (5.7) in the backward time direction, which is achieved by standard rigorous integrator of ODEs. Recall that we rewrite the system of differential equations (5.7) as in (5.8), that is, We integrate (5.7) backward in time. Taking ξ def = −η, we integrate (5.16) from 0 to ξ 0 with the initial point (x 1 (0), x 2 (0)) = p 0 = P(θ )| θ=−1 , which is on the local stable manifold W s loc ( p 2 ). The rigorous integrator we have used is mentioned in Remark 4.2. Furthermore, we rigorously compute the passing time in the original timescale using the following formula: where x 1 (ξ ) and x 2 (ξ ) denote the solution of (5.16).
In the present example, (5.7) is integrated with the initial point at the boundary of locally validated stable manifold, which is the boundary of the red curve in Fig. 4 with x 2 > 0, in the backward time direction and compute an enclosure of the evolution time in the original timescale: Extended stable manifold W s ( p 2 ) for (5.7) and corresponding blow-up times. The blue curve is the validated stable manifold W s ( p 2 ), while the black curve is the projection onto the (x 1 , x 2 )-plane. Numbers near points along the curve correspond to those shown in Table 1 where the rigorous enclosures of blow-up times are shown (Color figure online) x 1 x 2 Blow-up time  ( p 2 )). Note that the point in the figure where the corresponding blow-up time tends to infinity is the source equilibrium for (5.7), which corresponds to the bounded source for (5.1). Rigorous enclosures of t max on several sample points are shown in Table 1. Finally, we can reconstruct the true blow-up profile of the validated saddle-type blow-up solution through the directional compactification T d , which is drawn in Fig. 6. Note that this profile cannot be computed in the direct way since small perturbations of initial points violate the profile. 14 Remark 5.4 The integrand of t max has a different form from typical integrands shown in Sect. 2. Indeed, the integrand of (5.12) is a rational function consisting of two analytic functions. Nevertheless, the function x 1 (η) determining the denominator attains the value around 2 with sufficiently small error bounds so that the function 1/x 1 (η) 2 is analytic at x 1 (0), which is justified through the parameterization P, provided the trajectory {x 1 (η), s(η)} η∈[0,∞) is located on the interior of W s loc ( p 2 ). In particular, Proposition 4.4 and Theorem 4.5 can be still applied to showing that t max defined by  (0), v(0)) (5.12) depends analytically on initial points. Note that arguments in Sect. 5.1.3 directly confirm the analyticity of t max .

Example 2: Application to Higher-Dimensional Systems
The second example is the following (artificial) system in R 3 : ⎧ ⎪ ⎨ ⎪ ⎩ y 1 = y 1 (y 2 1 − 1), y 2 = y 2 1 y 2 + y 2 1 y 3 , y 3 = y 2 1 y 3 + δ −1 cy 2 1 y 3 − y 2 (y 2 − ay 1 )(y 1 − y 2 ) + wy 3 1 . (6.1) The present system is asymptotically homogeneous of order 3, namely asymptotically quasi-homogeneous of type α = (1, 1, 1). We thus apply the Poincaré-type compactification 15 to obtain the associated desingularized vector field as written by (2.16). In the present case, k = 2, n = 3, α j = β j = c = 1 for j = 1, . . . , n, and hence, derived by (2.12), is applied to determining (2.16). The concrete form iṡ The direct calculation of the Jacobian matrix of (6.3) is quite lengthy. Assuming that the Jacobian matrix off with respect to x is calculated, the Jacobian matrix of g with respect to x is calculated as follows: where δ i j is the Kronecker's delta. In the present case, the Jacobian matrix off is We observe that there are (at least) three equilibria on the horizon { p(x) 2 ≡ 3 i=1 x 2 i = 1}, one of which, denoted by p 0 , has a one-dimensional stable manifold and two of which, denoted by p 1 and p 2 , have two-dimensional stable manifolds. In the present study, we fix the following parameters: (a, c, δ, w) = (0.3, 0.7, 9.0, 0.02).

Fig. 7
Rigorously computed local stable manifolds for hyperbolic equilibria for (6.3). The C 0 rigorous error bound for the manifold around p 1 (left) is P − P (N ) ∞ ≤ r = 8.2 × 10 −9 with N = 50, while it is P − P (N ) ∞ ≤ r = 9.8 × 10 −10 with N = 60 for the manifold around p 2 (right) and P − P (N ) ∞ ≤ r = 9.8 × 10 −13 with N = 160 for the manifold around p 0 (center). The black dots are equilibria on the horizon; denoting p 1 , p 0 and p 2 from the left to the right We have computed the concrete position and associated eigenvalues, which are approximately given as follows: On the other hand, (6.2) possesses a source in a bounded region, namely The parameterization method applied to three equilibria on the horizon; p 0 , p 1 and p 2 , for (6.3) provides local stable manifolds with rigorous error enclosures. Distributions of these local stable manifolds are drawn in Fig. 7.

Blow-Up Time Computation
Since the compactification is homogeneous (namely α = (1, . . . , 1) for defining compactifications) and k = 2 in the present example, the maximal existence time t max is according to (2.15). Let P be a parameterization around x * ∈ E whose image of B n s determines the local stable manifold W s loc (x * ) of x * such that P(0) = x * . P is assumed to have a polynomial expression (cf. (4.3)) (α 1 , . . . , α n s ) ∈ Z n s ≥0 denotes the multi-index and θ α = θ α 1 1 · · · θ α ns n s . Assuming that the solution trajectory x(τ ) is on W s loc (x * ), the parameterization argument indicates that For a while, we further assume that Q = I , λ i ∈ R for i = 1, . . . , n s and k = 2. Then where (a * b) α denotes the discrete convolution over the multi-index α ∈ Z n s ≥0 given in (4.4) and θ α 0 = ((θ 1 ) 0 ) α 1 · · · ((θ n s ) 0 ) α ns . Here we use the fact where the denominator is strictly negative for all possible α and the analyticity of P ensures the convergence of the above series. Finally, we have the following expression of t max : Remark that the above expression makes sense only if by definition of the Poincaré compactification. With an explicit expression or enclosure of P(θ ), the quantity (6.5) or its enclosure is rigorously calculated for each θ 0 ∈ B n s . The above procedure is applied with n = 3 and n s = 1 or 2 in the present problem. If n s = 1, the expression (6.5) can be simplified by considering the single index l ≥ 1 instead of the multi-index α to obtain In practice, the computation of the Taylor coefficients a 1 , . . . , a n comes from a successful application of the Newton-Kantorovich type theorem (Theorem 3.5) applied to F : X → X given in (3.9) with X = ( 1 ) 2 = ( 1 1 ) 2 , namely ν = 1 (cf. Sect. 5.1.3). More precisely, denote byā 1 , . . . ,ā n the numerical approximations (of order N ) and r 0 > 0 such that the true coefficients satisfy a −ā X = max j=1,...,n a j −ā j 1 ≤ r 0 .
Denote b = a −ā and note that Denote the spectral gap of the stable eigenvalues by and note that σ gap = min |α|>0 |α · λ|. Hence, for all θ 0 ∈ B n s 1 , Similarly, we can show that then a rigorous enclosure of t max is given by the computable formula

Distribution of t max Near Blow-Up
In the present example, saddle equilibria p 1 and p 2 on the horizon both have two-dimensional stable manifolds. Once the parameterization method is applied to validating these invariant manifolds, the blow-up time t max defined by (6.4) is obtained as a function of the parameter θ determining local stable manifolds. In particular, we can validate distributions of t max on local stable manifolds. Figure 8 shows the distributions of t max . Because the vector field (6.3) itself can be defined outside D, namely in { x > 1} also, t max can attain negative values. Nevertheless, from the viewpoint that (6.3) is obtained from (6.1) through the compactification, only the positive values make sense as the blow-up time of solutions to (6.1). Now we pay attention to the following facts, which follow from fundamental arguments of compactifications (cf. Matsue 2018): • The horizon E is a codimension one invariant submanifold of R 3 .
• The integrand determining t max (e.g., (6.4)) is identically zero on E. Fig. 8 indeed reflect the above nature. For example, one-dimensional submanifold of two-dimensional stable manifolds of p 1 and p 2 are located on the horizon where t max is identically zero. Our computations further indicate that the region {t max > 0} is included in { x < 1}. Looking at the region {t max > 0}, like Fig. 3 in the previous example, we can discuss the distribution of blow-up times.

Results in
From our present observations, we have an interesting result about the distribution of blow-up times. In the present example, eigenvalues determining stable submanifolds on the horizon have smaller moduli than the transverse direction. In other words, the leading (stable) eigendirections are directed tangent to the horizon (red curves in Fig. 8) in both manifolds. Asymptotic behavior of trajectories around equilibria is therefore essentially determined by the exponential decay behavior in the direction parallel to the horizon. On the other hand, level sets of t max are distributed so that they are foliated parallel to the horizon, equivalently the level set t max = 0, in both cases. These observations may look strange from the viewpoint of the asymptotic behavior around (hyperbolic) equilibria. Indeed, dynamics around hyperbolic equilibria of interest are essentially governed by leading eigendirection, implying that the behavior along the leading eigendirection should mainly contribute to estimate t max . However, the integrand in (6.4) is almost zero near the horizon. More precisely, according to the proof of the blow-up criterion theorem (Theorem 2.8 whose proof is found in Matsue 2018), the integrand as a function of τ decays exponentially fast near the horizon. 16 Therefore, asymptotic behavior of solution trajectories near the horizon does little contributions to t max . As a consequence, blow-up time is essentially foliated parallel to the horizon, no matter where the leading eigendirection is distributed. This is a reason why the level set of t max is distributed parallel to the horizon.

Extension of Blow-Up Solutions
As demonstrated in Sect. 5, we can extend local stable manifolds globally by rigorous integration of (6.3) in backward time direction. In the present case, we have a (bounded) source equilibrium p b and we have succeeded in validating connecting orbits between three equilibria on the horizon and p b . The validated global stable manifolds are drawn in Fig. 9. These stable manifolds separate the asymptotic behavior of solution trajectories outside the manifolds, although we omit the detailed description of phase portraits because it is hard to clearly visualize.
Note that the present validation of connecting orbits is done by the method typically used in the similar works (e.g., Matsue and Takayasu 2020a). In particular, solutions approaching to trapping regions of equilibria are validated for the existence of globalin-time existence of solutions. In the present work, trapping regions of sink equilibria are validated by means of local Lyapunov functions (cf. Matsue and Takayasu 2020a), while the parameterization for sink equilibria can be also applied to constructing trapping regions.

Example 3: Presence of Separatrix Involving Blow-Ups
The final example is (7.1) The present vector field originally comes from the Keyfitz-Kranzer model (Kranzer and Keyfitz 1990) demonstrating a non-trivial example of system of conservation laws including singular shock waves. See Kranzer and Keyfitz (1990) or references therein Fig. 9 Rigorously computed trajectories on global stable manifolds of hyperbolic equilibria for (6.3). Local stable manifolds for (6.3) colored by pink and red are validated by the parameterization method (Fig. 7). The green dot denotes the (bounded) source equilibrium p b (Color figure  online) for details. A brief introduction of the model is also shown in Matsue (2018) (see also Sect. 8.1.1). Our purpose here is to validate blow-up solutions for (7.1) as well as bounded heteroclinic connections among bounded equilibria toward the global phase portrait. The present study unravels a significant characteristic of saddle-type blow-up solutions.
Firstly, a direct calculation yields the following.
Note that (7.1) is not quasi-homogeneous. On the other hand, the system (7.1) possesses the symmetry if (u(t), v(t)) is a solution to (7.1), then so is (−u(−t), v(−t)). This property is used to understand the global phase portrait of (7.1) including infinity.
To study the dynamics at infinity, we introduce the quasi-parabolic compactification of type (1, 2) given by Then the corresponding desingularized vector field g is given by the following: Fortunately, we know that all equilibria (including the origin) are hyperbolic, and hence, we do not need additional desingularization. Detailed information of our targeting equilibria is as follows: • The origin p 0 = (x 1 , x 2 ) = (0, 0), which is saddle. −0.7328506362011802, 0.537070054 9804747), which is sink.
• Equilibria on the horizon p ± ∞ = (x 1 , x 2 ) ≈ (±0.989136995894977, 0.2067585 57005180). The point p + ∞ is sink, while p − ∞ is source. Sample (non-rigorous) numerical computations indicate that there is a chain of global trajectories connecting p 0 and p + b , and p + b and p + ∞,s , respectively. The numerically computed global phase portrait including the horizon is shown in Fig. 10. The figure indicates that the whole phase space is separated into two subdomains by a heteroclinic chain among equilibria, including those on the horizon.

Remark 7.2
Here we have chosen the parabolic-type compactification in the present argument for the following reasons. First, our objective here is the global phase portrait for (7.1), which is insufficient to study only one local chart, namely directional compactifications. The change of coordinates by numerics (both in rigorous and non-rigorous sense) requires unnecessary and difficult tasks. Second, Poincarétype compactifications are inappropriate to study (7.1) including dynamics at infinity, because (7.1) is quasi-homogeneous only in the asymptotic sense, and the application to Poincaré-type compactifications to such a system cause the loss of regularity of the desingularized vector field on the horizon, as mentioned in Sect. 2.3.4.
One of our main goals here is to construct the chain, mainly connecting orbits among { p + ∞,s , p + b , p 0 }. Like in the previous examples, the local stable manifold W s loc ( p + ∞,s ) of the saddle p + ∞,s on the horizon can be validated by the parameterization method. Validated local stable manifolds of p + ∞,s as well as p 0 are shown in Fig. 11. These are validated through the parameterization method in the same way as Sects. 5 and 6. We omit the detailed implementation of the method applied to the present problem because the basic idea is identical, while we need lengthy calculations of terms we should enclose.
We then extend the manifold inside D ≡ {p(x) < 1} by the rigorous integration of (7.3). According to numerical simulations (Fig. 10), W s loc ( p + ∞,s ) is connected to the source p + b . Rigorous integration of (7.3) in backward time direction provide the computer-assisted validation of the connecting orbit from p + ∞,s to p + b by constructing a trapping region of p + b in backward time, which is a standard techniques for validating global-in-time trajectories and applied in, e.g., Matsue and Takayasu (2020a). On the other hand, we have another bounded equilibrium; the origin p 0 . Eigenvalue validation indicates that p 0 is a saddle, and the global trajectory connecting the source p + b and the origin p 0 is also validated by extending the local stable manifold W s loc ( p 0 ) of p 0 via the parameterization and the rigorous integration of (7.3) in backward time direction. By symmetry, we obtain the chain of connecting orbits among the points { p ± ∞,s , p ± b , p 0 }. Note that all these points are validated with rigorous errors through the parameterization method. Also note that the connecting orbit between p ± ∞,s exists through the fact that the horizon E is invariant and there are no equilibria between them (cf. Matsue 2018).
As a consequence, an invariant closed curve consisting of connecting orbits among equilibria { p ± ∞,s , p ± b , p 0 } is constructed, as indicated in Fig. 10, with computerassisted proof. The well-known Jordan's closed curve theorem indicates that the invariant closed curve decomposes the phase space D into two regions. 17 In the sequel, we study the nature of solutions through points on these separated regions from the viewpoint of blow-up behavior.

Blow-Up Time Computation
The maximal existence time of the solution (y 1 (t), y 2 (t)) for the original vector field is given as follows (see (2.19) and Matsue and Takayasu 2020a): and p + ∞,s from the left, respectively. Red curves are validated local stable manifolds with rigorous error bounds P − P (N ) ∞ ≤ r = 5.171 × 10 −14 for p 0 and P − P (N ) ∞ ≤ r = 1.381 × 10 −10 for p + ∞,s , respectively. The black curve denotes the horizon E. Although validated local stable manifolds are characterized for (7.3) which themselves make sense outside the horizon also, they make sense inside the horizon as the corresponding objects to the original vector field (7.1). In both validations, the approximation order N is chosen as N = 100. Because p + b is source, it does not admit a non-trivial stable manifold (Color figure online) Let x * = (x * ,1 , x * ,2 ) ∈ E be a saddle equilibrium. Note that x 4 * ,1 + x 2 * ,2 = 1 by definition of the present parabolic-type compactification.
As in the previous case, let P be a parameterization whose image of B n s determines the local stable manifold of x * such that P(0) = x * . P is assumed to have a power series expression (4.3) satisfying a 0 = x * . Assume that the trajectory {(x 1 (τ ), x 2 (τ ))} is included in W s loc (x * ) for the desingularized vector field. In the present case, n = 2, and we consider only the case n s = 1. Calculations below are slightly simplified by introducing u = e λτ θ 0 , where λ be the stable eigenvalue at x * . Indeed, we have Using that we have the following exact formula for t max : (7.4)

Chain of Connecting Orbits as Separatrix
In what follows, we discuss a global nature of saddle-type blow-up solutions in dynamical systems. In Fig. 10, we numerically observe that the compactified phase space is separated into four domains, one of which consists of points whose trajectories tend to the origin as τ → ±∞, while another consists of points whose trajectories tend to equilibria on the horizon as either both τ → ±∞, or only τ → −∞ or τ → +∞. Namely, the latter sets consist of initial points which solutions through these points blow up in finite times in the original coordinate. A significant importance of this observation is that these four domains are divided by sequences of trajectories including ones inducing blow-up solutions. In particular, saddle-type blow-up solutions themselves or bounded global-in-time trajectories connecting blow-up solutions can locally divide initial points into the above domains. As demonstrated in Sect. 6.3 and mentioned previously, connecting orbits between equilibria can be validated through the parameterization, extension of local (un)stable manifolds and construction of trapping regions (namely, local stable manifolds of sink equilibria). In two-dimensional systems like (7.1), the detailed nature of global dynamics can be easily considered by studying asymptotic behavior of solutions through neighborhoods of connecting orbits. Moreover, our validated connecting orbits involve blow-up solutions, and the characteristic value t max is associated with all points on validated connecting orbits and solutions close to them. Here we study connecting orbits involving hyperbolic saddles on the horizon and global-in-time solutions for the desingularized vector field, and the corresponding characteristics in the original vector field, yielding significantly different nature of asymptotic behavior. In particular, we investigate the following issues: • Dependence of blow-up characterizations on magnitude of initial points.
• Continuous dependence of t max on initial points.
(Local) stable manifolds of saddle equilibria locally separate neighborhoods of the equilibria, as well as those asymptotic behavior, unlike sink and source equilibria. The first issue is then equivalent to a non-trivial question here is whether such a separation around the horizon can significantly change the asymptotic behavior of solutions for the original vector field. Now we have a hyperbolic saddle on the horizon p + ∞,s , a bounded source p + b and the origin p 0 as a hyperbolic saddle. As shown in Figs. 10 and 11, local stable manifolds of p + ∞,s and p + b are validated through the parameterization method and extended through the integration of (7.3) like connecting orbits in Fig. 9. Let C sep be the union of validated connecting orbits: (7.5)

Dependence of Blow-Up Characterizations on Magnitude of Initial Points
First we consider the following issue.

Problem 7.3 Does the blow-up behavior depend on magnitudes of initial points?
In arguments of blow-up criteria, magnitudes (equivalently, norms) or values of several functionals of initial points are typically concerned for determining whether or not the corresponding solutions blow up. In many cases, there are mathematical arguments showing that initial points whose norms or associated functionals are sufficiently large induce finite-time blow-up. On the other hand, there are also several mathematical results of blow-up behavior which do not mention the magnitude of initial points. The aim of the present issue here is to reveal a qualitative characterization of asymptotic behavior around saddle-type blow-up solutions, which partially gives an answer to the above question. Now we choose two pairs of initial points. One pair is located close to p + ∞,s , while another pair is located close to the origin. In both pairs, two initial points are located at the opposite side to each other across C sep . More precisely, the former pair is chosen close to (x 1 , x 2 ) = (0.83, 0.53), while the latter pair is chosen close to (x 1 , x 2 ) = (0.32, 0.32). The corresponding points in the original coordinate are approximately (u, v) = (3.39444993, 4.69501202) and (u, v) = (0.36072017, 0.40662201), (7.6) respectively. Details are drawn in Fig. 12. The methodology shown in Sect. 4 is applied to validating global-in-time trajectories for (7.3) through each point, showing that the asymptotic behavior of trajectories are completely separated for both pairs of initial points. More precisely, • across saddle-type blow-up solutions, the asymptotic behavior of solutions as those for (7.1) significantly change, one of which attains t max = ∞, while another attains t max < ∞.
Moreover, we also observe that • such a nature can be observed even near the origin, where another connecting orbit between p 0 and p + b locally separates the phase space and is connected to the saddle-type blow-up solution generated by p + ∞,s . See Figs. 12 and 13. From the above observation, we can say that the magnitude of initial points is not always essential to determine the blow-up behavior. In other words, the chain C sep plays a role as the separatrix dividing global-in-time solutions and blow-up solutions. The key point is that the chain C sep including the saddle on the horizon locally separates the phase space and that there are sinks p + ∞ ∈ E and p − b ∈ D inducing global-in-time solutions for (7.3) approaching to them. The significant change of solutions in the original vector field is then responsible for the existence of saddletype blow-up solutions, in particular C sep , sinks on the horizon and another sinks on the other side of C sep . Nevertheless, saddle-type blow-up solutions themselves plays a role as the trigger of the above nature. Finally note that the present observation can be applied to other dynamical systems like (6.1), where the global extension of stable manifolds characterizing saddle-type blow-up solutions is validated in Sect. 6.3 (cf. Fig. 9).

Continuous Dependence of t max on Initial Points
Next we investigate the continuous dependence of t max on initial points across C sep given in (7.5). Here we consider a line segment which is transverse to C sep . See Fig. 14. The segment is chosen so that C sep and are orthogonal to each other at the boundary p 0,s of W s loc ( p + ∞,s ) validated by the parameterization method (cf. Fig. 11). The boundary p 0,s of W s loc ( p + ∞,s ) in D is then uniquely determined as the intersection C sep ∩ ≡ {p 0,s }. Our problem here is then stated as follows.
Problem 7.4 Does the blow-up time vary continuously on ? If not, study whether t max is discontinuous only in each side of C sep on , or discontinuous in both sides of C sep .
Indeed, the concrete dependence of t max cannot be unraveled unless explicit formulae (or both lower and upper bounds) for t max as functions of initial points are obtained. Our present methodology enables us to unravel this hidden nature in a reasonable way.
To study the above problem, the following steps are operated. The point p 0,s decomposes the line segment into two pieces, denoted by l and r consisting of points on in the left side (red in Fig. 14) and the right side (blue in Fig. 14) of p 0,s , respectively. Our validations, rigorous integrations of (7.3) in forward time direction, show that all sample points on r converge to p + ∞ as τ → ∞, which correspond to a family of sink-type blow-up solutions. Their validated blow-up times as well as the blow-up time of the solution through p 0,s are shown in Fig. 15 with their rigorous error bounds. Looking at Fig. 15, the corresponding blow-up times increase as sectional points on r become close to W s loc ( p ∞,s ). On the other hand, all points on l converge to the sink equilibrium p − b (Fig. 10). Because the preimage of p − b under the compactification is bounded, the corresponding solution in the original timescale exists for all t ≥ 0. This fact is easily confirmed by showing that t max ( p) = ∞ for p ∈ l . These observations show that t max is discontinuous as a function of points on at p 0,s from l . Next we discuss the continuity of t max at p 0,s on { p 0,s } ∪ r . Our validations show that t max ( p 0,s ) ∈ 3.109637008 441572 391221 , which is much higher than t max = t max ( p) through p ∈ r , according to Fig. 15. However, t max = t max ( p) drastically increases as p ∈ r approaches to p 0,s . At the point p ∈ with | p 0,s − p| = 1.0 × 10 −13 , validation of blow-up solutions did Recall that trajectories through points in the red region correspond to global-in-time solutions for (7.1), while trajectories through points in the blue region correspond to blow-up solutions for (7.1). A line is chosen so that it is transverse to C sep and is divided into two segments l (purple line) and r (black line) across C sep and it is orthogonal to C sep at p 0,s mentioned below. The intersection point { p 0,s } ≡ ∩ C sep is denoted by the green star (Color figure online) (a) (b) Fig. 15 Blow-up times of solutions with initial points on . We have totally chosen 10,000 points on r for validating t max . a Relationship of points on r and the blow-up times of solutions through those points. Horizontal: distance from p 0,s on . Vertical: blow-up time t max of the corresponding solution. The value 0 on the horizontal axis corresponds to p 0,s . The blow-up time t max = t max ( p) looks discontinuous at p = p 0,s . The red point denotes t max ( p 0,s ), while green points denote t max = t max ( p) at p ∈ r \{ p 0,s }. All plotted blow-up times here except t max ( p 0,s ) have rigorous error bounds less than 3.7235×10 −5 , while the rigorous error bound of t max ( p 0,s ) is 2.5175 × 10 −11 . b Enlarged view of the graph (a) for points within the distance ≤ 1.0 × 10 −8 from p 0,s . All plotted blow-up times here except t max ( p 0,s ) have rigorous error bounds less than 1.8288 × 10 −2 . As p ∈ r approaches to p 0,s , t max significantly increases. In the present study we do not have validations for t max associated with points p ∈ r within the distance ≤ 1.0 × 10 −13 (Color figure online) not succeed. As long as we have validated, we cannot conclude the discontinuity of t max at p 0,s in both sides. Nevertheless, we can still conclude that t max behaves in a singular manner around p 0,s where the trajectory approaches to different invariant sets as τ → ∞.
Remark 7.5 Rigorous enclosures of t max on W s ( p + ∞ ), namely sink-type blow-up solutions, are validated by local Lyapunov functions and rigorous integrations of (7.3), which are exactly machineries applied in Matsue and Takayasu (2020a), and hence, the detailed validation methodology is omitted. The difference of orders of (the worst) rigorous error bounds of t max on and off W s ( p ∞,s ) shown in Fig. 15 comes from that of the methodology for validating rigorous bounds of t max . Nevertheless, there is no significant influence on the qualitative tendency of t max in the present study. Remark 7.7 (Behavior of t max : a numerical experiment) We have numerically calculated the behavior of t max as a function of distance to the stable manifold in Fig. 15b. Let x be the distance of a point p from p 0,s in r and t max (x) be the corresponding blow-up time. As far as we have calculated, we could not match t max (x) by functions of the form x a , e ax , c(ln x) a and C x a (ln x) b for constants a, b, c. It is needless to say that this asymptotic form can be different for smaller x and a different choice of .

Short Summary of Our Observations
Our observations here are summarized as follows.
• Blow-up characterizations such as the asymptotic behavior and blow-up times do not always depend continuously on initial points in the presence of saddle-type blow-up solutions. • The blow-up time t max varies in a singular manner near the chain of connecting orbits involving saddle-type blow-ups, like C sep .
Note that these features cannot be unraveled only from local information around invariant objects, because local invariant manifolds themselves do not determine the asymptotic behavior of solutions through all points around the manifolds. In other words, global information of solutions are necessary to investigate this issue. It should be also noted that the above nature is observed not only by the presence of invariant sets like C sep , but also by the presence of another invariant sets like p + ∞ and p − b , at least one of which is included in the horizon E. This consequence strongly supports the importance of investigations of global dynamical structure to unravel the significantly different asymptotic behavior of solutions for the original vector field. Computerassisted proofs provide a systematic and mathematical rigorous way to investigate such global information of solutions. Moreover, the presence of saddle-type blow-up solutions provides an easy prediction of the existence of the above nature.

Concluding Remarks
In this paper, we have shown several characteristics of blow-up solutions for autonomous ODEs which are unstable under perturbations of initial points, referred to as saddle-type blow-up solutions, with the computer-assisted proofs of their existence and analytic characterization of blow-up times. Combining compactifications, timescale desingularizations of vector fields, parameterization of invariant manifolds and their extensions via ODE integrations with computer-assisted proofs, blow-up solutions and their extensions are validated systematically, no matter how stable equilibria on the horizon characterizing these blow-up solutions are. It should be noted that, as seen in all examples, our methodology does not require a priori information about the existence of blow-up solutions. This is a big advantage so that the present methodology can be applied to various dynamical systems and blow-up problems under mild assumptions.
Characteristics we have unraveled in the present paper are just examples of intrinsic natures which saddle-type blow-up solutions induce. But it is not an easy task to predict the presence of such features theoretically, because these are observed as the composite of multiple structures. For example, distribution of t max can be investigated by the combination of an analytic expression of t max and explicit distribution of local stable manifolds of equilibria on the horizon for desingularized vector fields. As for the separatrix nature among global-in-time solutions and blow-up solutions, it cannot be characterized without concrete distribution of global-in-time solutions, sink-type and saddle-type blow-up solutions. Computer-assisted proofs, on the other hand, connect features of explicitly validated objects to extract global nature as the composite of local characteristics, like the above features. These computation techniques efficiently work to gain insights into blow-up solutions.
We end this paper by leaving comments about topics involving saddle-type blowup solutions, which can relate to the present study toward further insights into global nature of blow-up solutions, dynamics at infinity and general finite-time singularities.

Remarks on Saddle-Type Blow-Up Solutions in Science and Engineering
Saddle-type blow-up solutions can arise in scientific and engineering studies. We review several preceding studies to assert the importance of saddle-type blow-up solutions and believe that our present methodology will contribute to unravel the dynamical nature of finite-time singularities involving saddle-type blow-up solutions in the following kinds of problems.

Singular Shock Waves
In the Riemann problem of the systems of conservation laws U t + f (U ) x = 0 (8.1) for some smooth f : R n → R n , namely the initial value problem of (8.1) with shock waves are characterized by locally integrable (weak) solutions with discontinuities with the constraints called jump conditions or the Rankine-Hugoniot conditions. With the assumption of viscous shock criterion, the Riemann problem is reduced to find connecting orbits of the traveling wave ODE associated with (8.1) connecting U L and U R . In the 1980s and 1990s, shock waves with a singular nature on the front were observed for a simple system of conservation laws, which are referred to as delta shocks or singular shocks. Roughly speaking, singular shocks are characterized by shocks with Dirac's delta singularity on the shock front (see, e.g., Keyfitz et al. 2003;Kranzer and Keyfitz 1990; Sever 2007 for precise discussions of delta shocks and singular shocks). A typical feature of singular shocks with the presence of the deltalike singularity is that several constraints in jump conditions are violated, 18 which is referred to as the presence of the Rankine-Hugoniot deficit of a shock measuring the magnitude of singularity on the shock front. From the viewpoint of dynamical systems, there is a characterization of singular shocks (e.g., Schaeffer et al. 1993), showing that singular shocks can consist of a collection of blow-up solutions and "invariant sets at infinity". In several concrete problems such as the Keyfitz-Kranzer model (Kranzer and Keyfitz 1990) and the two-phase model (Keyfitz et al. 2003), the geometric singular perturbation theory plays a key role in characterizing singular shocks as a singular perturbation of blow-up connections for the traveling wave problems associated with the original conservation laws with the regularization keeping the self-similarity of waves (well known as Dafermos regularization). Preceding studies with blowing up (desingularization) of singularities and the geometric singular perturbation theory indicate that singular shocks are characterized by trajectories approaching to normally hyperbolic invariant manifolds, corresponding to the infinity for appropriately transformed dynamical systems (Hsu 2016;Schecter 2004). 19 We believe that saddle-type blow-up solutions can play key roles in characterizing such singular nature both qualitatively and quantitatively (e.g., Rankine-Hugoniot deficits).

Suspension Bridge
The equation of the following form is well studied as a model expressing scientific and engineering phenomena: where k ∈ R is a parameter and f is a locally Lipschitzian. This equation arises in the dynamical phase space analogy of a nonlinearly supported elastic structure (Hunt et al. 1989) and a model characterizing pattern formations in physical, chemical and biological systems (Bonheure and Sanchez 2006). See also, e.g., Peletier and Troy (2012). In Berchio et al. (2011), a possible finite-time blow-up for the solution of (8.2) is discussed with a mild assumption f ∈ Lip loc (R), f (t)t > 0 for every t ∈ R\{0}.
A fundamental result involving blow-up is that the existence of a blow-up solution w(t) for (8.2) as t → t max < ∞ implies that lim inf namely a blow-up with oscillation. Moreover, the existence of the above oscillatory blow-up for (8.2) with a specific nonlinearity f is proved. There are several reports about the relationship between the system (8.2) to traveling waves for the model equation of a suspension bridge proposed by Lazer-McKenna (Lazer and McKenna 1990). According to many preceding works and historical sources, one of the most interesting behaviors for suspension bridges (including the Tacoma Narrow Bridge where was collapsed in November 1940) is the following: Large vertical oscillations can rapidly change, almost instantaneously, to a torsional oscillation (quotation from Gazzola and Pavani 2011).
Preceding works involving this catastrophic phenomenon discuss the mechanism of torsional oscillations in detail, 20 one of which is considered to be the oscillatory blowup behavior mentioned above. It should be noted that there is another direction to the origin of such torsional oscillations. In Arioli and Gazzola (2015), it is explained that internal resonances can trigger the torsional instability. Later successive works (e.g., Gazzola and Pavani 2013) have reported the qualitative nature of the above blow-up such as infinitely many change of signs before blow-up, vanishing intervals of oscillations via several quantitative estimates. In order to obtain the nature, several growth conditions of f (but generalized under these conditions unlike Gazzola and Pavani 2011), restrictions to k and an inequality for derivatives of solution w at an initial time are assumed. It should be noted that norms of initial points are not essential to characterize the above behavior. See Gazzola and Pavani (2013) for details. Recently, the first author and collaborators (D'Ambrosio et al. 2015) have characterized the above blow-up nature for particular nonlinearity f in (8.2) by constructing a concrete asymptotic form of blow-up profiles and validating a periodic solution with computer-assisted proofs. In D' Ambrosio et al. (2015), it is also validated that the periodic solution for an auxiliary equation is unstable, which indicates that the corresponding blow-up solution is unstable under perturbations of initial points. It is thus expected that the blow-up nature which is unstable under perturbations of initial points plays a key role in describing rich and interesting, sometimes catastrophic, scientific and engineering nature. Remark 8.1 In Matsue (2018), it is proved that blow-up behavior with wide oscillations like (8.3) can be characterized by periodic orbits at infinity, which is referred to as a periodic blow-up. More precisely, global trajectories on the stable manifold of a hyperbolic periodic orbit on the horizon for the desingularized vector field correspond to blow-up solutions with oscillations whose asymptotic behavior, such as the blow-up rate and the oscillatory nature, are uniquely determined by the order of the original vector field and the periodic orbit on the horizon. The fundamental machinery for this characterization is the same as that shown in Sect. 2. Arguments in the present paper will also contribute to reveal universal mechanisms of this kind of blow-up solutions which are saddle-type both quantitatively and qualitatively, and their validations.

More Comments
We leave several comments about the link to blow-up behavior arising in the suspension bridge problem. As noted, it is proved in D' Ambrosio et al. (2015) with the computer assistance that there is an unstable hyperbolic periodic orbit = {w(t)} expressing an asymptotic behavior of blow-up behavior for (8.2) with specific k and f . It is then conjectured in D' Ambrosio et al. (2015) that, for the appropriately transformed dynamics from the problem of the form (8.2), the boundary of the basin of attraction of the origin coincides with W s ( ). A consequence of the conjecture is the existence of a three-dimensional manifold which "separates" the phase space and for which solutions with initial points taken on one side of the manifold blow-up in finite time while on the other side, solutions converge to the origin. In the present paper, we have focused on unstable, in particular saddle-type, blow-up solutions which can extract the above nature. We have revealed here that saddle-type blow-up solutions, even with the simpler asymptotic behavior than (D'Ambrosio et al. 2015), can separate the phase space so that initial points on one side determine global-in-time solutions, while those on the other side induce blow-up solutions. We have mainly investigated asymptotic behavior of solutions near a chain of connecting orbits for desingularized vector fields including saddles on the horizon, like C sep given in (7.5) for (7.3), and shown that C sep triggers the above significantly different asymptotic behavior among solutions. In particular, C sep have played a role as a separatrix among solutions for the original vector field. We believe that such invariant objects can characterize the "boundary" of the basin of attraction mentioned in D' Ambrosio et al. (2015).
Note that the above object is characterized only for stationary blow-up (Theorems 2.4, 2.8 and 2.12) so far. On the other hand, a computer-assisted proof of the existence of (un)stable manifolds of hyperbolic periodic orbits is already established in, e.g., (Castelli et al. 2018), and the treatment of blow-up solutions involving periodic orbits at infinity is also established in Matsue (2018Matsue ( , 2019. In other words, the same machinery as shown in Sect. 2 can be applied. Going back to the suspension bridge problem, combination of preceding works with the arguments in the present paper can contribute to unravel the nature of blow-up behavior in (8.2) only with a few mild assumptions.