Numerical validation of blow-up solutions with quasi-homogeneous compactifications

Abstract

We provide a numerical validation method of blow-up solutions for finite dimensional vector fields admitting asymptotic quasi-homogeneity at infinity. Our methodology is based on quasi-homogeneous compactifications containing quasi-parabolic-type and directional-type compactifications. Divergent solutions including blow-up solutions then correspond to global trajectories of associated vector fields with appropriate time-variable transformation tending to equilibria on invariant manifolds representing infinity. We combine standard methodology of rigorous numerical integration of differential equations with Lyapunov function validations around equilibria corresponding to divergent directions, which yields rigorous upper and lower bounds of blow-up time as well as rigorous profile enclosures of blow-up solutions.

Appendices

A Schur decompositions

In this section we review Schur decompositions of squared matrices.

Proposition A.1

(Schur decomposition, e.g., [11]). Let $A\in M_n(\mathbb {C})$: complex $n\times n$ matrix. Then there exists a unitary matrix $Q\in U(n)$ such that

$$\begin{aligned} Q^H A Q = T \equiv D + N, \end{aligned}$$

where $Q^H$ is the Hermitian transpose of Q, $D = \mathrm{diag}(\lambda _1,\ldots , \lambda _n)$ and $N\in M_n(\mathbb {C})$ is strictly upper triangular. Furthermore, Q can be chosen so that the eigenvalues $\lambda _i$ appear in any order along the diagonal. We shall call Ta Schur normal form ofA.

When we treat all computations in real floating number or interval arithmetic, the real version of Schur decompositions can be applied.

Proposition A.2

(Real Schur decomposition, e.g., [11]). Let $A\in M_n(\mathbb {R})$: real $n\times n$ matrix. Then there exists an orthogonal matrix $Q\in O(n)$ such that

$$\begin{aligned} Q^T A Q = T\equiv \begin{pmatrix} R_{11} &{}\quad R_{12} &{}\quad \ldots &{}\quad R_{1m} \\ 0 &{}\quad R_{22} &{}\quad \ldots &{}\quad R_{2m} \\ \vdots &{}\quad \vdots &{}\quad \ddots &{}\quad \vdots \\ 0 &{}\quad 0 &{}\quad \ldots &{}\quad R_{mm} \end{pmatrix}, \end{aligned}$$

where each $R_{ii}$ is either a $1\times 1$ or a $2\times 2$ matrix having complex conjugate eigenvalues. We shall call Ta real Schur normal form ofA.

A merit of Schur decompositions is that we can apply it to arbitrary square matrices. In particular, change of coordinates via Schur decompositions can be realized no matter what the multiplicities of any eigenvalues are.

B concrete calculations of an upper bound of $t_{\max }$ with quasi-parabolic compactifications

In this section, we consider the rigorous validation of the maximal existence time

$$\begin{aligned} t_{para,\max } = \int _0^\infty \kappa _{para}(y)^{-k}\left( 1-\frac{2c-1}{2c}\kappa _{para}(y)^{-1}\right) d\tau _{para} \end{aligned}$$

of solution trajectories with quasi-parabolic compactifications and computer assistance. First of all, we compute the following integral representing the time of integration of computed trajectory for desingularized vector fields int-timescale in advance:

$$\begin{aligned} t_{para,N}=\int _0^{\tau _{para,N}}\left( 1-p(x(\tau _{para}))^{2c}\right) ^{k}\left( 1-\frac{2c-1}{2c}\left( 1-p(x(\tau _{para}))^{2c}\right) \right) d\tau _{para}. \end{aligned}$$

As mentioned in Sect. 4.3, the estimate of $|1-p(x)^{2c}|$ is essential to computation of an upper bound $C_{n,\alpha ,N}(L)$. At first, we derive the estimate with the type $\alpha =(1,2)$ and $c=2$ as an example. Let $x_*= (x^*_1, x^*_2) \in \mathcal {E}$ and assume that a Lyapunov function L(x) is validated in a vicinity of $x_*$. Then

$$\begin{aligned} x_1^4+x_2^2 =&(x_1-x_1^*+x_1^*)^4+(x_2-x_2^*+x_2^*)^2\\ =&(x_1-x_1^*)^4+4(x_1-x_1^*)^3x_1^*+6(x_1-x_1^*)^2(x_1^*)^2+4(x_1-x_1^*)(x_1^*)^3+(x_1^*)^4\\&+(x_2-x_2^*)^2+2(x_2-x_2^*)x_2^*+(x_2^*)^2. \end{aligned}$$

Now $p(x_*)=1$ holds since $x_*\in \mathcal {E}$. Thus we have

$$\begin{aligned} \left| 1-p(x)^{2c}\right| =&\Big |(x_1-x_1^*)^4+4(x_1-x_1^*)^3x_1^*+6(x_1-x_1^*)^2(x_1^*)^2+4(x_1-x_1^*)(x_1^*)^3\\&+(x_2-x_2^*)^2+2(x_2-x_2^*)x_2^*\Big |\\ =&\left| \left[ 4(x_1^*)^3~2x_2^*\right] \left[ \begin{array}{l} x_1-x_1^*\\ x_2-x_2^*\end{array}\right] +\left[ 6(x_1^*)^2~1\right] \left[ \begin{array}{l} (x_1-x_1^*)^2\\ (x_2-x_2^*)^2 \end{array}\right] \right. \\&\left. +\left[ 4x_1^*~0\right] \left[ \begin{array}{l} (x_1-x_1^*)^3\\ (x_2-x_2^*)^3 \end{array}\right] +\left[ 1~0\right] \left[ \begin{array}{l} (x_1-x_1^*)^4\\ (x_2-x_2^*)^4 \end{array}\right] \right| \\ \le&\left\| \left[ \begin{array}{l} 4(x_1^*)^3\\ 2x_2^*\end{array}\right] \right\| \left\| x-x_*\right\| \\&+\max \left\{ 6(x_1^*)^2,1\right\} \left\| x-x_*\right\| ^2+\left| 4x_1^*\right| \left\| x-x_*\right\| ^3+\left\| x-x_*\right\| ^4. \end{aligned}$$

By $\Vert x-x_*\Vert \le \left( c_1L\right) ^{1/2}$ followed by the value of Lyapunov function L(x), we obtain

$$\begin{aligned} \left| 1-p(x)^{2c}\right|&\le \left\{ 16(x_1^*)^6+4(x_2^*)^2\right\} ^{1/2}\left( c_1L\right) ^{1/2}\\&\quad +\max \left\{ 6(x_1^*)^2,1\right\} c_1L+\left| 4x_1^*\right| \left( c_1L\right) ^{3/2}+\left( c_1L\right) ^{2}\\&=:C_{n,\alpha ,N}(L). \end{aligned}$$

Finally we obtain an upper bound of $t_{para,\max }$ as follows:

$$\begin{aligned} t_{para,\max }&=t_{para,N}+\int _{\tau _{para,N}}^\infty \left( 1-p(x(\tau _{para}))^{2c}\right) ^{k}\\&\quad \left( 1-\frac{2c-1}{2c}\left( 1-p(x(\tau _{para}))^{2c}\right) \right) d\tau _{para}\\&=t_{para,N}+\int _{\tau _{para,N}}^\infty \left( 1-p(x(\tau _{para}))^{2c}\right) ^{k}\\&\quad \left( \frac{1}{2c}+\frac{2c-1}{2c}p(x(\tau _{para}))^{2c}\right) d\tau _{para}\\&\le t_{para,N}+\int _{\tau _{para,N}}^\infty \left| 1-p(x(\tau _{para}))^{2c}\right| ^{k}d\tau _{para}\\&\le t_{para,N} +\frac{1}{c_{\tilde{N}}c_1} \int _0^{L(x(\tau _{para,N}))} \frac{C_{n,\alpha ,N}(L)^{k}}{L}dL, \end{aligned}$$

where we have used the estimate $\frac{dL}{d\tau _{para}}\le -c_{\tilde{N}}c_1L$ along the trajectory $\{x(\tau _{para})\}$, which follows from the inequality of Lyapunov functions. The positive constants $c_{\tilde{N}}, c_1$ are shown in [21].

Next we show an estimate of $|1-p(x)^{2c}|$ with compactifications of general type $\alpha =(\alpha _1,\dots ,\alpha _n)$. As in the previous case, let $x_*= (x^*_1, \ldots , x^*_n) \in \mathcal {E}$ and assume that a Lyapunov function L(x) is validated in a vicinity of $x_*$. Then

$$\begin{aligned} \left| 1-p(x)^{2c}\right|&=\left| 1-\sum _{i=1}^nx_i^{2\beta _i}\right| \\&=\left| \sum _{i=1}^n\left( x_i^*\right) ^{2\beta _i}-\sum _{i=1}^n\left( x_i-x_i^*+x_i^*\right) ^{2\beta _i}\right| \\&=\left| \sum _{i=1}^n\sum _{j=1}^{2\beta _i}\left( {\begin{array}{c}2\beta _i\\ j\end{array}}\right) \left( x_i-x_i^*\right) ^{j}\left( x_i^*\right) ^{2\beta _i-j}\right| \\&=\left| \sum _{j=1}^{\max \left\{ 2\beta _i\right\} }v_j^T\left[ \begin{array}{c} (x_1-x_1^*)^j\\ (x_2-x_2^*)^j\\ \vdots \\ (x_n-x_n^*)^j \end{array}\right] \right| , \end{aligned}$$

where $v_j\in \mathbb {R}^n$ is the vector given by

$$\begin{aligned} \left( v_j\right) _i=\left\{ \begin{array}{ll} \left( {\begin{array}{c}2\beta _i\\ j\end{array}}\right) \left( x_i^*\right) ^{2\beta _i-j}&{}(j\le 2\beta _i),\\ 0&{}(j>2\beta _i). \end{array}\right. \end{aligned}$$

Thus we have

$$\begin{aligned} \left| 1-p(x)^{2c}\right|&=\left| \sum _{j=1}^{\max \left\{ 2\beta _i\right\} }v_j^T\left[ \begin{array}{c} (x_1-x_1^*)^j\\ (x_2-x_2^*)^j\\ \vdots \\ (x_n-x_n^*)^j \end{array}\right] \right| \\&\le \Vert v_1\Vert \Vert x-x_*\Vert +\sum _{j=2}^{\max \left\{ 2\beta _i\right\} }\Vert v_j\Vert _\infty \Vert x-x_*\Vert ^j\\&\le \Vert v_1\Vert \left( c_1L\right) ^{1/2}+\sum _{j=2}^{\max \left\{ 2\beta _i\right\} }\Vert v_j\Vert _\infty \left( c_1L\right) ^{j/2}=:C_{n,\alpha ,N}(L), \end{aligned}$$

where we have used $\Vert x-x_*\Vert \le \left( c_1L\right) ^{1/2}$.

If $k=1$, which is the case shown in Sect. 5.3, then an upper bound estimate of $t_{\max }$ is realized as follows, for example:

$$\begin{aligned} t_{para,\max }&\le t_{para,N} +\frac{1}{c_{\tilde{N}}c_1} \int _0^{L(x(\tau _{para,N}))} \frac{C_{n,\alpha ,N}(L)}{L}dL\\&= t_{para,N} +\frac{1}{c_{\tilde{N}}} \int _0^{L(x(\tau _{para,N}))} \left\{ \Vert v_1\Vert \left( c_1L\right) ^{-1/2}+\sum _{j=2}^{\max \left\{ 2\beta _i\right\} }\Vert v_j\Vert _\infty \left( c_1L\right) ^{j/2-1}\right\} dL\\&= t_{para,N} +\frac{1}{c_{\tilde{N}}} \left\{ 2\Vert v_1\Vert c_1^{-1/2}L(x(\tau _{para,N}))^{1/2}+\sum _{j=2}^{\max \left\{ 2\beta _i\right\} }\frac{2}{j}\Vert v_j\Vert _\infty c_1^{j/2-1}L(x(\tau _{para,N}))^{j/2}\right\} . \end{aligned}$$

C Proofs of statements

1.1 C.1 Proof of Proposition 2.4

Now compute the Jacobian matrix of T for verifying its bijectivity. Direct computations yield

$$\begin{aligned} \frac{\partial x_i}{\partial y_j} = \kappa ^{-\alpha _i}\left( \delta _{ij} - \kappa ^{-1} \alpha _i y_i \frac{\partial \kappa }{\partial y_j}\right) \end{aligned}$$

with the matrix form

$$\begin{aligned} J&= \left( \frac{\partial x_i}{\partial y_j}\right) _{i,j = 1,\ldots , n} = A_{\alpha } \left( I_n - \kappa ^{-1} y_\alpha (\nabla \kappa )^T\right) ,\\ A_\alpha&= \mathrm{diag}(\kappa ^{-\alpha _1}, \ldots , \kappa ^{-\alpha _n}),\quad y_\alpha = (\alpha _1 y_1,\ldots , \alpha _n y_n)^T. \end{aligned}$$

We following arguments in [8], for any (column) vectors $y,z\in \mathbb {R}^n$, to have

$$\begin{aligned} (I_n + \beta yz^T)(I_n + \delta yz^T)&= I_n + (\beta + \delta )yz^T + \beta \delta yz^T yz^T\\&= I_n + (\beta + \delta + \beta \delta \langle z,y\rangle )yz^T, \end{aligned}$$

so $I_n+\delta yz^T = (I_n +\beta yz^T )^{-1}$ if $\delta = -\beta / (1 + \beta \langle z,y\rangle )$.

In this case, we choose $\beta = -\kappa ^{-1}, y = y_\alpha , z = \nabla \kappa$ and have

$$\begin{aligned} \left( \frac{\partial y_j}{\partial x_i}\right) = \left( \frac{\partial x_i}{\partial y_j}\right) ^{-1} = \left( I_n - \frac{1}{\kappa - \langle y_\alpha , \nabla \kappa \rangle } y_\alpha (\nabla \kappa )^T\right) A_{\alpha }^{-1} \end{aligned}$$

By (A3) we have $\kappa > \langle y_\alpha , \nabla \kappa \rangle$, which indicates that the transformation T as well as $T^{-1}$ are well-defined and $C^1$ locally bijective including $y=0$. On the other hand, the map T maps any one-dimensional curve $y = (r^{\alpha _1}v_1, \ldots , r^{\alpha _n}v_n)$, $0\le r < \infty$, with some fixed direction $v\in \mathbb {R}^n$, into itself (cf. [16]). For continuous mappings from $\mathbb {R}$ to $\mathbb {R}$, local bijectivity implies global bijectivity. Consequently, (A3) guarantees also the global bijectivity of Tinto$T(\mathbb {R}^n)\subset \mathcal {D}$. Finally we prove that T is onto. First let $y\in \mathbb {R}^n\setminus \{0\}$. Then the correspondence

$$\begin{aligned} \iota : y_i \mapsto \frac{y_i}{p(y)^{\alpha _i}} \end{aligned}$$

maps $\mathbb {R}^n\setminus \{0\}$ onto the set $\{p(y) = 1\}$. If $p(y) = 1$ then $\kappa (y)$ attains a constant $\kappa _1 > 1$ from (A0). From $\kappa (y)^{2c}p(x)^{2c} = p(y)^{2c}$, we know that the compactification T maps the set $\{p(y) = 1\}$ onto the set $\{x\in \mathcal {D} \mid p(x) = 1/\kappa _1\}$. In particular, the set $\mathbb {R}^n\setminus \{0\}$ is mapped onto $\{x\in \mathcal {D} \mid p(x) = 1/\kappa _1\}$ via the map $T\circ \iota$. Therefore the surjectivity of T is reduced to that on the ray $C_y = \{(r^{\alpha _1}y_1, \ldots , r^{\alpha _n}y_n)\mid 0\le r < \infty \}$ for each $(y_1, \ldots y_n)\in \mathbb {R}^n$ with $p(y) = 1$. By definition $T(C_y)$ is

$$\begin{aligned} T(C_y) = \left\{ (x_1, \ldots , x_n) = \left( \frac{r^{\alpha _1}}{q(r)^{\alpha _1}}y_1, \ldots , \frac{r^{\alpha _n}}{q(r)^{\alpha _n}}y_n \right) \mid 0\le r < \infty \right\} . \end{aligned}$$

From $\kappa (y)^{2c}p(x)^{2c} = p(y)^{2c}$, we have $p(x) = r/q(r)$ on $T(C_y)$. From (A0) and (A1), for any value $c_x\in (0,1)$, we can choose the value $r \in (0,1)$ so that $r/q(r) = c_x$. Correspondingly we can define y from x on $T(C_y)$. Obviously, $T(\mathbf{0})=\mathbf{0}$ and hence $T:\mathbb {R}^n\rightarrow \mathcal {D}$ is onto and the proof is completed.

Note that the condition (A2) is not actually used in the present argument.

1.2 C.2 Proof of Lemma 3.1

Suppose that $y_*$ is an equilibrium of (2.1), i.e., $f(y_*) = 0$. Then the right-hand side of (3.1) obviously vanishes at the corresponding $x_*$.

Conversely, suppose that the right-hand side of (3.1) vanishes at a point $x\in \mathcal {D}, p(x) < 1$: namely,

$$\begin{aligned} f(\kappa x) - \kappa (y)^{-1}\langle \nabla \kappa , f(\kappa x) \rangle y_\alpha = 0. \end{aligned}$$

Multiplying $\nabla \kappa$, we have

$$\begin{aligned} \langle \nabla \kappa , f(\kappa x) \rangle \left( 1- \kappa (y)^{-1}\langle \nabla \kappa , y_\alpha \rangle \right) = 0. \end{aligned}$$

Due to (A3), we have $|\kappa (y)^{-1}\langle \nabla \kappa , y_\alpha \rangle | < 1$ and hence $\langle \nabla \kappa , f(\kappa x) \rangle = 0$. Thus we have $f(y) = f(\kappa x) = 0$ by the assumption.

1.3 C.3 Proof of Lemma 3.2

By admissibility (A1)-(A2), we have

$$\begin{aligned} (\nabla \kappa (y))_i \sim \frac{1}{\alpha _i}\frac{y_i^{2\beta _i-1}}{\kappa (y)^{2c-1}} = \frac{1}{\alpha _i}\frac{\kappa ^{\alpha _i(2\beta _i-1)}x_i^{2\beta _i-1}}{\kappa ^{2c-1}} = \frac{1}{\alpha _i}\frac{x_i^{2\beta _i-1}}{\kappa ^{\alpha _i-1}} \quad \text { as }\quad p(y)\rightarrow \infty , \end{aligned}$$

where we used the condition $\alpha _j \beta _j \equiv c$ for all j from (2.2). Therefore the vector field (3.3) near infinity becomes

$$\begin{aligned} x_i'&\sim \kappa ^k \tilde{f}_i(x) - \alpha _i x_i \sum _{j=1}^n \frac{1}{\alpha _j}\frac{x_j^{2\beta _j-1}}{\kappa ^{\alpha _j-1}} \kappa ^{k+\alpha _j - 1}\tilde{f}_j(x)\nonumber \\&= \kappa ^k \left\{ \tilde{f}_i(x) - \alpha _i x_i \sum _{j=1}^n \frac{x_j^{2\beta _j-1}}{\alpha _j} \tilde{f}_j(x)\right\} \quad \text { as }\quad \kappa \rightarrow \infty . \end{aligned}$$

(C.1)

Since $\tilde{f}_i$ is O(1) as $\kappa \rightarrow \infty$, then right-hand side of (C.1) is $O(\kappa ^k)$ as $\kappa \rightarrow \infty$.

1.4 C.4 Proof of Theorem 3.7

The property $b = \sup \{t\mid y(t) \text { is a solution of (2.1)}\}$ corresponds to the property that

$$\begin{aligned} \sup \{\tau \mid x(\tau ) = T(y(t)) \text { is a solution of } (3.5) \text { in the time variable } \tau \} = \infty . \end{aligned}$$

Indeed, if not, then $\tau \rightarrow \tau _0 < \infty$ and $\lim _{\tau \rightarrow \tau _0-0}x(\tau ) = x_*$ as $t\rightarrow b-0$. The condition $x(\tau ) = x_*$ is the regular initial condition of (3.5). The vector field (3.5) with the new initial point $x(\tau ) = x_*$ thus has a locally unique solution $x(\tau )$ in a neighborhood of $\tau _0$, which contradicts the maximality of b. Therefore we know that $\tau \rightarrow +\infty$ as $t\rightarrow b-0$. Since $\lim _{\tau \rightarrow \infty }x(\tau ) = x_*$, then $x_*$ is an equilibrium of (3.5) on $\partial {\mathcal {D}}$. The similar arguments show that $t\rightarrow a+0$ corresponds to $\tau \rightarrow -\infty$ and that the same consequence holds true.

1.5 C.5 Proof of Proposition 3.9

Each $\tilde{f}_j(x)$ given by (3.2) with $\kappa = \kappa _{para}$ is $C^1$ on $\overline{\mathcal {D}}$, since all terms of $\tilde{f}_j$ are multiples of powers of $(1-p(x)^{2c})$ and smooth asymptotically quasi-homogeneous terms in $f_j(y)$. Consequently, we know that the right-hand side of (3.7) is $C^1$ on $\overline{\mathcal {D}}$.