On the numerical solutions for a parabolic system with blow-up

1.   Introduction

We consider axisymmetric solutions for the semilinear parabolic system

where $\triangle$ denotes the Laplace operator, $p, q > 1$ and $\Omega\subseteq\mathbb{R}^N$ is the unit ball with the origin as its center. Since we consider the radially symmetric solutions, letting $r = |x|$, (1.1) reads as

Here, the initial data $u_0(r), v_0(r)$ are assumed to be nonnegative. It was proved (see ^[20,21]) that the solutions of (1.2) may become unbounded in a finite time $T$. This phenomenon is known as blow-up and the finite time $T$ is called the blow-up time.

Friedman & Giga ^[19] proved, for $N = 1$ and $p = q$, that the solutions of (1.1) blow up at a single point if the initial data are of bell-shaped, that is, symmetric with a single peak. Later, Souplet ^[32] proved that blow-up only occurs at the origin for the initial-boundary value problem (1.2) if the initial data $u_0(r), v_0(r)\geq 0$ is decreasing in $(0, 1)$ and

In fact, (1.3) is known to hold if one assumes that either

or

See for instance ^[8,16,32]. For more results concerning (1.3), one may also consult ^[5,18]. Moreover, among these results, it was proved in ^[16] that the solution of (1.1) for any bounded domain satisfies

namely,

for some positive constants $C_i\; (i = 1, 2, 3, 4)$ if we assume

for certain $a\in(0, 1)$. In ^[32], (1.4) was also shown to be true for the solution of (1.2) under the assumptions $u_r, v_r\leq 0$ and $u_t, v_t\geq 0$.

Although the solution is known to become unbounded in a finite time $T$ in the sense of $L^\infty$-norm, that is,

it is also interesting to investigate whether the solution blows up in other measures. Quittner and Souplet ^[29] proved, for the solutions of (1.1), that

Later, Souplet ^[32] extended (1.7) for the solutions of (1.2) to the case of $\sigma_1 = \frac{N(pq-1)}{2(p+1)}$ and $\sigma_2 = \frac{N(pq-1)}{2(q+1)}$ under the assumptions $u_r, v_r\leq 0$ and $u_t, v_t\geq 0$.

In this paper, we focus on the numerical aspects for the problem (1.2). The first numerical study for blow-up problems can be traced back to the paper ^[26] by Nakagawa. He considered the semilinear heat equation

in the case of $p = 2$ and space dimension $N = 1$ and proposed a numerical scheme with adaptive temporal increments to compute an approximate blow-up time. His idea was generalized in ^{[6,9,13,14,22]}, in which the asymptotic behaviors for the numerical solutions of (1.8) in space dimension one were also analyzed. For the cases of space dimensions $N\geq 2$, numerical approximation was first studied by Chen ^[7]. He again used Nakagawa's adaptive strategy and considered a finite difference scheme for the axisymmetric solutions of (1.8) to reproduce the phenomenon of finite-time blow-up. The numerical blow-up sets were also classified. Later, Groisman ^[22] proposed a fairly general numerical scheme for (1.1) whose temporal grid size is also defined adaptively and derived the blow-up rate and the blow-up set for his numerical solutions. Recently, the author proposed schemes with both adaptive and uniform time meshes for the computation of blow-up solutions of (1.8) and showed that our schemes can faithfully reproduce several blow-up behaviors, including the blow-up set, blow-up rate and blow-up of $L^\sigma$-norm. See ^[11,15] for the detail.

Besides the finite difference approximation, it is worth mentioning that ^[27,28] used the finite element method (FEM) to compute the numerical solution and an approximate blow-up time for the solution of (1.8) with convergence proofs. In fact, FEM is a strong tool which can deal with problems in general domains and is also often used to numerically resolve the concentration of the singularity of PDEs. See for instance ^[4,31,33].

It is known that adaptive temporal increments can reproduce the phenomenon of finite-time blow-up very well. See for instance ^{[1,6,7,9,13,22,26]}. Such a strategy, however, can not faithfully reproduce several blow-up behaviors such as the blow-up rate, blow-up curve and so on. We refer the readers to ^{[11,12,15,22]} for the details. We thus use uniform temporal increments and the algorithm proposed in ^[10] for the computation of blow-up solutions for (1.2).

The rest of the paper is organized as follows. As a first step to numerically analyze the solutions of (1.2), we first consider a finite difference analogue for the ODE system

for the computation of approximate blow-up times in Section 2. Then we propose a finite difference scheme to compute blow-up solutions of (1.2) in Section 3. Convergence for the numerical solution and the numerical blow-up time are also derived in the same section. In Section 4, we show how our numerical solutions reproduce various blow-up behaviors such as the blow-up sets, the blow-up rates and $L^\sigma$-norm blow-ups for both $u$ and $v$. Finally, the paper is ended up with a conclusion.

2.   An ODE problem

In this section, we consider the ODE system (1.9) and its finite difference analogue. It is not difficult to see that the solutions $u$ and $v$ blow up simultaneously in a finite time. In fact, multiplying respectively the first equation and the second equation of (1.9) by $u^q$ and $v^p$, one has

which implies

It then follows

where

This implies that the solutions $u$ and $v$ become unbounded at the same time if there is a blow-up. Substituting (2.1) into (1.9), we have

where

Since $pq > 1$, it is not difficult to derive

where

Now we consider the following finite difference scheme for (1.9):

Here, we use $\left(\cdot\right)_t$ to denote the forward difference operator. $\Delta t$ is the grid size, $t_n = n\Delta t\; (n\geq 0)$ are the grid points and $U^n, V^n$ denote the approximations for $u(t_n), \; v(t_n)$ respectively. Since $U^0 = u_0, \; V^0 = v_0 > 0$, one has

Moreover, multiplying $(U^{n+1})^q$ and $(V^n)^p$ to the first and the second equations of (2.5) respectively, we have

Theorem 2.1. Let $U^n, V^n$ be the solution of (2.5). Let $T_0 < T_O$ be given arbitrarily. Then

Since the proof can be carried out in a standard way, we thus omit it.

To define the numerical blow-up time, we need the following lemma.

Lemma 2.1. It holds that

In particular, one has

and

Proof. By (2.6), one has for all $n\geq 0$

which implies (2.7).

Now (2.8) and (2.9) are direct results of (2.7) and (2.5).

Lemma 2.2. It holds that

and

Proof. First, one has, by (2.5),

which implies

On the other hand, one has again by (2.5)

which implies

Observe that $G_1(U^0) = (v_0)^p > 0$. By (2.8) and the monotonicity of $G_1$ and $\{U^n\}$, one has

which, together with (2.5), also implies

Now we use the algorithm proposed in ^[10] to define the numerical blow-up time. We import a strictly increasing function $H:(0, \infty)\rightarrow (0, \infty)$, which is to be determined, satisfying $\lim\limits_{s\rightarrow \infty}H(s) = \infty$. Given any $\Delta t > 0$, since $U^n, V^n\rightarrow \infty$ as $n\rightarrow \infty$, there exist $n^u_{\Delta t}, \; n^v_{\Delta t}\in\mathbb{N}$ such that

and

We define $T^u_O(\Delta t) = \Delta t\cdot n^u_{\Delta t}$ and $T^v_O(\Delta t) = \Delta t\cdot n^v_{\Delta t}$ and call them the numerical blow-up time of $u$ and $v$ respectively.

Theorem 2.2. Assume that $H$ satisfies

Then $\lim\limits_{\Delta t\rightarrow 0}T^u_O(\Delta t) = \lim\limits_{\Delta t\rightarrow 0}T^v_O(\Delta t) = T_O$.

Proof. By virtue of (2.8) and (2.11), the proof is similar to [Theorem 3.2, ^[10]]. We outline the proof for $T^u_O(\Delta t)$ for the readers' convenience. The proof for $T^v_O(\Delta t)$ can be done in exactly the same way.

We complete the proof by showing that

First, we consider the finite difference scheme

It is easy to show, by (2.8), that $U^n\geq Y^n\; (\forall n\geq 0)$. Let $\bar{n}^u_{\Delta t}\in\mathbb{N}$ be the positive integer such that

Then one has $T^u_O(\Delta t) = \Delta t\cdot n^u_{\Delta t}\leq \Delta t\cdot\bar{n}^u_{\Delta t}$, which, together with [Theorem 2.1, ^[10]], gives

We thus have

Here use has been made of (2.14).

On the other hand, we assume that

Then there exists $\{\tau_k\}$ such that $\tau_k\rightarrow 0$ as $k\rightarrow \infty$ and

Let $\{U^n(\tau_k)\}$ be the solutions corresponding to $\tau_k$. Then

while the exact solution $u(t)$ remains bounded in $\left[0, \frac{T_O+\underline{T^u_O}}{2}\right]$, which is a contradiction to Theorem 2.1. Therefore,

and we are done.

Remark 2.1. Observe that $G_1(s)\sim s^\frac{p(q+1)}{p+1}$ and $G_2(s)\sim s^\frac{q(p+1)}{q+1}$. Therefore, (2.14) can be replaced by

It is easy to see that $H(s) = s^\gamma\; (\gamma > 0)$ satisfies (2.16) and thus can be used to compute an approximate blow-up time. In addition, by (2.15), one has

where $C$ denotes a generic constant. The error bound

can also be derived in a similar way.

Although we can only prove the error bounds (2.17) and (2.18) currently, our numerical results seem to suggest that

and

hold true. In fact, our computational results suggest that

and

In the following example, we set $p = 3, \; q = 2$. The initial data is given by $v_0 = 3$, $u_0 = \left(\frac{q+1}{p+1}v_0^{p+1}\right)^\frac{1}{q+1}$ so that $C_0$ in (2.2) vanishes. Then the blow-up time $T_O$ can be computed explicitly by $T_O = \int_3^\infty\frac{ds}{\left(\frac{3}{4}s^4\right)^\frac{2}{3}}\approx 0.1164774$. Figures 1–3 illustrate the convergence and the convergence order when different choices of $H$ are used to compute the numerical blow-up times of $u$ and $v$. In addition, we remark that we use different $H$ to determine the numerical blow-up times for $u$ and $v$ in each numerical experiment. This is necessary when we want to compute also the asymptotic behaviors of the numerical solutions. The reason will become clearer in Section 4.

3.   Finite difference solutions for (1.2)

In this section, we consider a finite difference scheme to the parabolic system (1.2). To discretize (1.2), we borrow the idea proposed in ^[7] for the axisymmetric solution of the semilinear heat equation (1.8).

We discretize (1.2) as follows. Let $J\in\mathbb{N}$. $\Delta r = 1/J$ is the spatial grid size and $r_j = j\Delta r, \; (j = 0, \cdots, J)$ are the spatial grid points. $\Delta t > 0$ is the temporal grid size and $t_n = n\Delta t\; (n\geq 0)$ are the temporal grid points. Let $N_0 = \left\lfloor\frac{N-1}{2}\right\rfloor$, where $\left\lfloor\frac{N-1}{2}\right\rfloor$ is the largest integer among those that are smaller or equal to $\frac{N-1}{2}$. $U^n_j, V^n_j$ are the approximations of $u, v$ at $(t_n, r_j)$ respectively. At the origin, we discretize as follows due to the boundary condition $u_r(t, 0) = v_r(t, 0) = 0$:

For $j = 1, \cdots, N_0$, we consider the following discretization:

while, for $j = N_0+1, \cdots, J-1$, we consider

The boundary conditions are given by

and the initial data are defined by $U^0_j = u_0(r_j)$ and $V^0_j = v_0(r_j)$ for $j = 0, \cdots, J$. Note that (3.1)(3.2) are nothing but approximations for axisymmetric solutions of

by the central difference at $r = 0$. (3.3)(3.4) are proposed for stability in the cases of $N > 3$. One may consult ^[7] for the detail.

Remark 3.1. We remark that ^[15] also proposed a finite difference scheme for axisymmetric solutions of (1.8), which can reproduce the phenomenon of single-point blow-up if adaptive temporal increments are applied. The discretization, however, is more complicated than Chen's idea. Since we do not pursue reproducing the phenomenon of single-point blow-up by schemes with adaptive temporal meshes in the current paper, we use an explicit version of Chen's idea for our analysis for simplicity. Nevertheless, one can also use the recipe proposed in ^[15] to discretize (1.2) and investigate numerically the blow-up behaviors. All the results given below can be derived parallel.

Remark 3.2. Although we do not consider adaptive spatial increments in this paper, an adaptive spatial mesh indeed enhances numerical resolution of singularities. See for instance ^[2,23,24,30] and the references therein. It is, however, very difficult to prove convergence and analyze asymptotic behaviors mathematically for those schemes. Our scheme is simple but is able to reproduce asymptotic blow-up profile to some extent with rigorous convergence proofs.

Lemma 3.1. Let $U^n_j, \; V^n_j$ be the solution of (3.1)–(3.7). Let $\lambda\equiv\frac{\Delta t}{(\Delta r)^2}\leq\frac{1}{2N}$ be fixed. Assume that $u_0(r), v_0(r)\geq 0$ for all $r\in[0, 1]$. Then $U^n_j, V^n_j\geq 0$ for all $j = 0, \cdots, J$ and $n\geq 0$.

Proof. The proof is completed by induction. Assume that $U^n_j, V^n_j\geq 0$ for all $j = 0, \cdots, J$. Then it follows directly from (3.1)(3.3) and the assumption $\lambda\leq\frac{1}{2N}$ that $U^{n+1}_j\geq 0$, for $j = 0, \cdots, N_0$, which, together with (3.2)(3.4), subsequently imply the nonnegativity of $V^{n+1}_j\; (j = 1, \cdots, N_0)$. For $j = N_0+1, \cdots, J-1$, one has by (3.5)

Since $j\geq N_0+1 > \frac{N-1}{2}$, the coefficients on the right-hand side of the above equation are nonnegative, which implies $U^{n+1}_j\geq 0$. Then one has by (3.6)

Thus, we are done.

Theorem 3.1. Let $T$ be the blow-up time of (1.2) and $T_0 < T$ be given. Let $\lambda\leq\frac{1}{2N}$ be fixed. Then

and

Since the proof is a standard one, we omit it. We refer the readers to ^[7,15] for the convergence proofs of a single equation.

Lemma 3.2. Let $U^n_j, \; V^n_j$ be the solution of (3.1)–(3.7). Let $\lambda\leq\frac{1}{3N}$ be fixed. Assume that $u_0(r), v_0(r)\geq 0$ and $u_0(r), v_0(r)$ are monotonically decreasing in $(0, 1)$. Then $U^n_j\geq U^n_{j+1}\geq 0$ and $V^n_j\geq V^n_{j+1}\geq 0$ for all $j = 0, \cdots, J-1$ and $n\geq 0$.

Proof. We give the detail for the case of $N > 3$ and the cases of $N = 2, 3$ can be carried out similarly. Assume that $U^n_j\geq U^n_{j+1}\geq 0$ and $V^n_j\geq V^n_{j+1}\geq 0$ for all $j = 0, \cdots, J-1$. We first show $U^{n+1}_j\geq U^{n+1}_{j+1}\; (j = 0, \cdots, J-1)$. To this end, it suffices to show $U^{n+1}_0\geq U^{n+1}_1, \; U^{n+1}_{N_0}\geq U^{n+1}_{N_0+1}$. For other $j'$s, the arguments given in Lemma 3.1 can be applied to derive $U^{n+1}_j-U^{n+1}_{j+1}\geq 0$.

By (3.1)(3.3), one has

and

from which follow

For $j = N_0$, one has by (3.3)(3.5)

and

which imply

Now $V^{n+1}_j\geq V^{n+1}_{j+1} > 0$ can be proved in a similar way and the proof is completed by induction.

From now on, we assume that $u_0, v_0\geq 0$ and $u_0, v_0$ are monotone decreasing in $(0, 1)$. To illustrate how blow-up is reproduced numerically, we assume in addition that the initial data satisfy a discrete analogue of (1.6)

$\rm(A)$ There exists $a\in(0, 1)$ such that $\Delta_d U^0_j+(1-a)(V^0_j)^p\geq 0$ and $\Delta_d V^0_j+(1-a)(U^1_j)^q\geq 0$ for all $j = 0, \cdots, J-1$.

Here, the operator $\Delta_d$ is defined by

Remark 3.3. Observe that the assumption (A) implies

and

This can be regarded as a discrete analogue for

which is a sufficient condition of the finite-time blow-up for the solutions of (1.2). See for instance ^[32]. We remark that Deng ^[16] also derived the blow-up rate for the solution of (1.2) under the assumption of the initial data (1.6). Although (1.6) is more restrictive than (3.9) because of the positive parameter $a$, it simplifies the analysis. We thus consider those initial data satisfying (A) in the following discussion.

Define, for $0\leq j\leq J-1$ and $n\geq 0$,

Note that, by (3.1)–(3.6), $W^n_j$ and $Z^n_j$ can also be written as

Since $U^n_J = V^n_J = 0\; (\forall n\geq 0)$, we set $W^n_J = Z^n_J = 0\; (\forall n\geq 0)$.

Lemma 3.3. Let $U^n_j, \; V^n_j$ be the solutions of (3.1)-(3.7). Let $\lambda\leq\frac{1}{3N}$. Suppose that there exists $a\in(0, 1)$ such that (A) holds. Then, for all $j = 0, \cdots, J-1$ and $n\geq 0$,

That is, $W^n_j, Z^n_j\geq 0$ for all $j = 0, \cdots, J-1$ and $n\geq 0$. In particular, one has

where $\|U^n\|_\infty = \max\limits_{j = 0, \cdots, J}|U^n_j|$ and $\|V^n\|_\infty = \max\limits_{j = 0, \cdots, J}|V^n_j|$.

Proof. Assume that $W^n_j, Z^n_j\geq 0\; (j = 0, \cdots., J-1)$ holds for certain $n$. Observe that

and

imply

Since $Z^n_j\geq 0\; (\forall j)$, one has by Lemma 3.1 $V^{n+1}_j\geq V^n_j$, which implies

On the other hand, one has by Lemma 3.2

Consequently,

This, together with the assumption $\lambda\leq\frac{1}{3N} < \frac{1}{2N}$, implies the nonnegativity of $W^{n+1}_j\; (\forall 0\leq j\leq J-1)$.

Similarly, one can easily derive

By virtue of the fact $W^{n+1}_j\geq 0$ and Lemma 3.1, 3.2, we have

from which follow

Again, the stability condition $\lambda\leq\frac{1}{3N}$ yields $Z^{n+1}_j\geq 0$. Now the proof of (3.10) is completed by induction.

(3.11) is a direct result of Lemma 3.2 and (3.10).

By Lemma 3.3, it is easy to see

and

Let $H:(0, \infty)\longmapsto(0, \infty)$ be a strictly increasing function satisfying $\lim\limits_{s\rightarrow \infty}H(s) = \infty$. Then for any given $\Delta t > 0$, there exist positive integers $n^u_{\Delta t}, n^v_{\Delta t}\in\mathbb{N}$ such that

and

We define the numerical blow-up time for the solutions $u$ and $v$ of (1.2) by $T^u(\Delta t) = \Delta t\cdot n^u_{\Delta t}$ and $T^v(\Delta t) = \Delta t\cdot n^v_{\Delta t}$.

Theorem 3.2. Let $T$ be the blow-up time for the solution of (1.2) and let $U^n_j, \; V^n_j$ be the solutions of (3.1)–(3.7). Let $\lambda\leq\frac{1}{3N}$ be fixed. Assume that the assumption (A) holds true and that the function $H$ satisfies (2.16). Then we have $\lim\limits_{\Delta t\rightarrow 0}T^u(\Delta t) = \lim\limits_{\Delta t\rightarrow 0}T^v(\Delta t) = T.$

Proof. We write down the prove of $\lim\limits_{\Delta t\rightarrow 0}T^u(\Delta t) = T$. Convergence of $T^v(\Delta t)$ can be shown in parallel.

We complete the proof by showing

Assume that $\underline{T^u} < T$. Then there exists $\{(\tau_i, h_i)\}$ satisfying $0 < \lambda = \frac{\tau_i}{h_i^2}\leq \frac{1}{3N}$, $\tau_i, h_i\rightarrow 0$ as $i\rightarrow \infty$ and that

Since, by (3.12), $\tau_i\cdot H\left(\|U^{n^u_{\tau_i}}\|_\infty\right)\geq 1$, we have

while the exact solution $u(t, x)$ remains bounded in $\left[0, \frac{\underline{T^u}+T}{2}\right]$. This is a contradiction to Theorem 3.1. We thus have $\underline{T^u}\geq T$.

On the other hand, we assume that $\overline{T^u} > T$. Then there exists $\{(\tau_i, h_i)\}$ satisfying $0 < \lambda = \frac{\tau_i}{h_i^2}\leq \frac{1}{3N}$, $\tau_i, h_i\rightarrow 0$ as $i\rightarrow \infty$ and that

Given any $M > 0$, Theorem 3.1 guarantees the existence of sufficiently small $(\Delta t, \Delta r)\in\{(\tau_i, h_i)\}$ and a positive integer $n_0$, depending on $\Delta t$ and $\Delta r$, such that

For simplicity, we denote the solutions of (3.1)–(3.7) corresponding to the grid sizes $\Delta t, \Delta r$ by $U^n_j\equiv U^{n}_j(\Delta t, \Delta r), V^n_j\equiv V^{n}_j(\Delta t, \Delta r)$. We remark that $n_0 < n^u_{\Delta t}$ since

Now we consider the finite difference system

By (3.11), it is easy to show

Let $R(\Delta t)$ be the numerical blow-up time for $A^n$. Namely, $R(\Delta t) = \Delta t\cdot(m_{\Delta t}-n_0)$, where $m_{\Delta t}\in\mathbb{N}$ is determined by

Since $H$ is strictly increasing, one has $H\left(\|U^n\|_\infty\right)\geq H(A^n)\; (\forall n\geq n_0)$, which implies

Note that we have, by Theorem 2.2,

Since we can choose sufficiently large $M$ and sufficiently small $\Delta t$ such that

it then follows

which is a contradiction. Thus, we are done.

In the following example, we set $N = 5, \; p = 3, \; q = 2, \; \lambda = 1/(3N)$. The initial data are given by $u_0(r) = 150(1+\cos(\pi r)), \; v_0(r) = 100(1+\cos(\pi r))$. In Figures 4 and 5, we choose

for the computation of the numerical blow-up time.

Note that, under the assumption (3.9), the solution of (1.2) satisfies, for all $r\in[0, 1]$,

Consequently, (3.14) are considered to be suitable choices for the computation of blow-up solutions $u$ and $v$ of (1.2) respectively. See the discussions in ^[10,11]. Our computational results also support this viewpoint. In fact, we also compute the numerical blow-up times with different choices of $H$. Note that it takes more steps for the numerical solutions to satisfy (3.12) or (3.13) if a smaller $H$ is chosen. That is, a smaller $H$ results in a larger numerical blow-up time. Figure 6 shows the computational results of the numerical blow-up times for

and

We remark that both (3.15) and (3.16) satisfy (2.16). In the case of (3.15), the numerical blow-up time converges from above. This suggests that the computation is overcalculated and should be stopped at an earlier step. A larger $H$, compared with (3.15), is considered to be better. On the other hand, the convergence is from below for the choice (3.16), which implies insufficiency of our computation for an approximate blow-up time. That is, a smaller $H$ is better.

4.   Numerical blow-up behaviors

We now pay our attention to the blow-up behaviors for the numerical solutions of (3.1)–(3.7). In this section, we always assume that the initial data $u_0, v_0$ are nonnegative, monotonically decreasing and satisfy assumption (A).

4.1.   Blow-up set

Numerical blow-up sets were first discussed by Chen ^[6]. He proposed a finite difference scheme whose temporal increment is defined adaptively for the one-dimensional semilinear heat equation

It is well-known that the solution of (4.1) blows up at a single point if bell-shaped initial data is considered. See for instance ^[17,34]. However, Chen showed that the blow-up points for his numerical solution may blow up at more than one point. Later, more complete results, including the cases of the space dimension $N\geq 2$, were derived in ^{[7,9,13,14,22]}. All these results considered numerical schemes with adaptive temporal increments, a strategy proposed by Nakagawa ^[26], and suggested that single-point blow-up can not be faithfully reproduced. As a matter of fact, the numerical blow-up sets contain more than one point if $p\in (1, 2]$. In this regard, it is worth mentioning that the authors ^[15] proposed a finite difference scheme, which again used Nakagawa's adaptive algorithm, for radially symmetric solutions of (1.8) in space dimension $N\geq 2$ and reproduced the single-point blow-up phenomenon successfully.

In this section, we would like to discuss the numerical blow-up sets in a different way. Instead of adaptive temporal increments, we use uniform ones for the computation of blow-up solutions. See ^[11] for the application to (4.1).

It was proved in ^[32] that the solution of (1.2) blows up only at the origin under our assumptions on the initial data. To see how single-point blow-up is reproduced, we show numerically that, for arbitrarily given $x_0\in(0, 1]$, the numerical solution remains bounded at $x_0$ at the numerical blow-up time as $\Delta t\rightarrow 0$. We proceed our computation as follows:

● Given $\Delta t, \Delta r > 0$ satisfying $\lambda\leq \frac{1}{3N}$, we first compute the numerical blow-up time $T^u(\Delta t)$ and $T^v(\Delta t)$.

● Let $j_{\Delta r}$ satisfy $j_{\Delta r}\Delta r\leq x_0 < (j_{\Delta r}+1)\Delta r$. Then we compute $U^n_j$ and $V^n_j$ at $(T^u(\Delta t), r_{j_{\Delta r}})$ and $(T^v(\Delta t), r_{j_{\Delta r}})$ respectively.

● Let $\Delta t, \Delta r\rightarrow 0$ to see whether or not the numerical solution remains bounded at $r_{j_{\Delta r}}$.

Here, we choose (3.14) for the computation of $T^u(\Delta t)$ and $T^v(\Delta t)$ respectively. In the following example, we set $N = 5$, $(p, q) = (3, 2)$, and $u_0(r) = 150(1+\cos(\pi r)), \; v_0(r) = 100(1+\cos(\pi r))$. Let $x_0 = 0.01$. Our computational results show that the numerical solution remains bounded at $r = x_0$ as $\Delta t, \Delta r\rightarrow 0$. See Figure 7. Similar results can be observed for $x_0 = 0.001$ and $x_0 = 0.0001$.

We also compute the behaviors of $U^n_j$ and $V^n_j$ at $(T^u(\Delta t), r_{j_{\Delta r}})$ and $(T^v(\Delta t), r_{j_{\Delta r}})$ respectively as $\Delta t\rightarrow 0$ with (3.15) and (3.16). As Figure 8 shows, the boundedness can be observed numerically from the computational results with (3.14) and (3.15). Especially in the case of (3.15), the numerical solution at $x_0 = 0.01$ decreases as $\Delta t\rightarrow 0$. As for the case (3.16), note that, by (3.10), $U$ and $V$ are monotonically increasing in $n$ at each spatial grid points so that the numerical solution will be smaller if a large $H$ is used for the computation. Although the numerical solution computed with (3.16) seems to be increasing as $\Delta t\rightarrow 0$, they are always bounded by the ones computed with (3.14) or (3.15).

Concerning the relation between the boundedness of the numerical solution (3.1)-(3.7) and the exact solution (1.2), we have the following theorem.

Theorem 4.1. $\rm (a)$ Assume that

Then

$\rm (b)$ Assume that

Then

Proof. We outline the proof for $u$. Assume on the contrary that $\limsup\limits_{t\rightarrow T}u(t, x_0) = \infty.$ Then given $\varepsilon > 0$, there exists $t_\varepsilon\in(T-\varepsilon, T)$ such that $u(t_\varepsilon, x_0) > 2M$. Observe that

For small $\Delta t > 0$, let $n^\varepsilon_{\Delta t}$ be the positive integer satisfying that

Since $u$ remains smooth on $\left[0, \frac{T+t_\varepsilon}{2}\right]\times[0, 1]$ and, by Theorem 3.1, we also have

it follows that, for sufficiently small $\Delta t$ and $\Delta r$,

which implies

On the other hand, since $\lim\limits_{\Delta t\rightarrow 0}T^u(\Delta t) = T$, one has, for sufficiently small $\Delta t$,

Now (3.10) yields

for all sufficiently small $\Delta t$. This contradicts the assumption $\limsup\limits_{\Delta t\rightarrow 0}U^{n^u_{\Delta t}}_{j_{\Delta r}} = M.$

Thanks to this theorem, if one can know the boundedness of the numerical solution at certain spatial point $x_0$ at the numerical blow-up time as $\Delta t\rightarrow 0$, it suggests the boundedness of the exact solution $u(t, x_0), \; v(t, x_0)$ as $t\rightarrow T$. From this point of view, a smaller $H$ seems to provide numerical results which can suggest numerically the boundedness of the numerical solution more strongly. (See Figure 8). However, it should be noted that such a choice may lead to unnecessary calculation when computing the numerical blow-up time. (See Figure 6).

4.2.   Blow-up rate

In this subsection, we show that our numerical solutions satisfy a discrete analogue of (1.5). Put $H(s) = s^\gamma\; (\gamma > 0)$.

Theorem 4.2. $\rm(a)$ Let $\gamma\geq \frac{pq-1}{p+1}$. Then there exist constants $c_1, c_2 > 0$ such that

$\rm(b)$ Let $\gamma\geq \frac{pq-1}{q+1}$. Then there exist constants $c_3, c_4 > 0$ such that

We sketch out the proof for (4.2). (4.3) can be carried out in a similar way.

To show the validity of (4.2), we need the following lemmas.

Lemma 4.1. Let $\gamma\geq \frac{pq-1}{p+1}$. There exist $C, c > 0$ such that

Proof. Note that, by (3.11) and (3.1)(3.2),

Since $V^{k}_0\geq V^{k-1}_0\; (\forall k\geq 1)$,

which, together with (4.5), implies

We thus have

Note that (3.12) implies

and that

Since $\{U^n_0\}$ is increasing in $n$ and $U^n_0\rightarrow \infty$ as $n\rightarrow \infty$, it then follows

for some positive constant $C$.

On the other hand, note that

which, together with (4.5), implies

It then follows that there exists $c > 0$ such that

Lemma 4.2. Assume that $\gamma\geq\frac{pq-1}{p+1}$. Let $z_0 = 1$ and $z_{k+1}\; (k\geq 0)$ be the positive root of

Then $\{z_k\}$ is decreasing. Moreover, we have

Proof. Since $f_k$ is monotonically increasing and satisfies $f_k(0) = -z_k < 0$, $f_{k}(z_{k+1}) = 0$ and $f_k(z_k) > 0$, one has $0 < z_{k+1} < z_k$.

We now prove (4.6). Assume that (4.6) holds for certain $k$. By virtue of (4.4), one has

Then it follows

or equivalently,

By the monotonicity of $f_k$ and the fact $f_k(z_{k+1}) = 0$, we thus have

Lemma 4.3. Assume that $\gamma\geq\frac{pq-1}{p+1}$. Let $y_1 = 1$ and $y_{k+1}\; (k\geq 1)$ is the positive root of

Then $\{y_k\}$ is decreasing. Moreover, we have

Proof. Since $g_k$ is monotonically increasing and satisfies $g_k(0) = -y_k < 0$, $g_{k}(y_{k+1}) = 0$ and $g_k(y_k) > 0$, the monotonicity of $\{y_k\}$ follows.

Now we assume that (4.7) holds for certain $k$. By virtue of (4.4), one has

from which follows

or equivalently,

By the monotonicity of $f_k$ and the fact $f_k(y_{k+1}) = 0$, we thus have

and the proof is completed by induction.

Proof of (4.2).

Note first that $\{z_k\}, \{y_k\}$ are decreasing sequences. By the definition of $\{z_k\}$ and $\{y_k\}$, we have

and

which imply

and

It is not difficult to derive that there exist $c_1, c_2 > 0$ such that

Let $k = n^u_{\Delta t}-n$. Then one has

which, together with (4.6)(4.7) and the above inequality, implies (4.2).

4.3.   $L^\sigma$-norm blow-ups

For blow-up solutions of (1.2), the $L^\infty$-norm becomes unbounded in a finite time $T$. However, whether $L^\sigma$-norms $(1\leq \sigma < \infty)$ also blow up in the finite time $T$ or not depends on the space dimension $N$ and the parameters $p, q$ of the nonlinear terms. In fact, (1.7) gives that $L^\sigma$-norm blows up simultaneously with the $L^\infty$-norm if $\sigma$ is large enough. It should be noted that, however, whether the $L^\sigma$-norm becomes unbounded as $t\rightarrow T$ or remains bounded in $[0, T)$ seems to remain open for small $\sigma$. In this subsection, we would like to numerically examine this problem. To this end, we first define discrete $L^\sigma$-norms by

We compute as follows. As in the computation of blow-up set, we use (3.14) to compute numerical blow-up times $T^u(\Delta t)$ and $T^v(\Delta t)$ respectively. Then we compute the discrete $L^\sigma$-norm for $U$ and $V$ at $T^u(\Delta t)$ and $T^v(\Delta t)$ respectively and let $\Delta t\rightarrow 0$ to see whether the discrete $L^\sigma$-norms tend to infinity or remain bounded.

We use the same example as was used in Section 3 and 4.1 to illustrate our computational results. We set $(p, q) = (3, 2)$, $N = 5$, and $u_0(r) = 150(1+\cos(\pi r)), \; v_0(r) = 100(1+\cos(\pi r))$. As a consequence, (1.7) suggest that $\|u(t, \cdot)\|_{L^\sigma}$ becomes unbounded at the blow-up time if $\sigma > \frac{25}{8}$ and $\|v(t, \cdot)\|_{L^\sigma}$ blows up at the blow-up time if $\sigma > \frac{25}{6}$. As shown in Figures 9 and 10, the computational results suggest that

and $\|U^{n^u_{\Delta t}}\|_3, \; \|V^{n^v_{\Delta t}}\|_4$ remain bounded as $\Delta t\rightarrow 0$. This suggests that $L^\sigma$-norms for both $u$ and $v$ might perhaps remain bounded for small $\sigma$.

As a matter of fact, we have the following theorem:

Theorem 4.3. $\rm(a)$ Let $H^u(s) = s^{\gamma_1}\; (\gamma_1 > 0)$. Then it holds that, for all $\sigma > \frac{N \gamma_1}{2}$,

$\rm(b)$ Let $H^v(s) = s^{\gamma_2}\; (\gamma_2 > 0)$. Then it holds that, for all $\sigma > \frac{N \gamma_2}{2}$,

Proof. By (4.8) and (3.12) we have

provided that $\sigma > \frac{N\gamma_1}{2}$. Similarly, we have by (3.13)

if $\sigma > \frac{N\gamma_2}{2}$.

Remark 4.1. If we choose $\gamma_1 = \frac{pq-1}{p+1}$ and $\gamma_2 = \frac{pq-1}{q+1}$, then Theorem 4.3 is consistent with (1.7). This suggests that (3.14) might perhaps be the most appropriate choice for the computation to determine the $L^\sigma$-norm blow-up for the solutions of (1.2).

Moreover, we have the following theorem.

Theorem 4.4. Let $1\leq \sigma < \infty$.

$\rm(a)$ Suppose that

Then

$\rm(b)$ Suppose that

Then

By virtue of Theorem 3.1 and 3.2, Theorem 4.4 can actually be proved in exactly the same way as was given in ^[15]. We thus omit it.

5.   Conclusions

We compute the blow-up time and several blow-up behaviors for (1.2) by a finite difference scheme. To obtain a better approximation for the blow-up phenomenon, we compute the numerical blow-up time $T^u(\Delta t)$ and $T^v(\Delta t)$ for $u$ and $v$ and investigate various asymptotic behaviors for the numerical solution $U^n_j$ and $V^n_j$ at $T^u(\Delta t)$ and $T^v(\Delta t)$ respectively. This is necessary in order to observe numerically the blow-up features for $u$ and $v$ near the blow-up time. Although two numerical blow-up times for (1.2) are computed, we remark that

Here, $T$ is the blow-up time of (1.2).

We use functions $H^u$ and $H^v$ to determine the numerical blow-up times for the solutions $u$ and $v$ of (1.2). It is natural to ask which choice is the best. For the computation of the numerical blow-up time, by Theorem 3.2, the function $H(s) = s^\gamma\; (\forall \gamma > 0)$ can be used to determine both $T^u(\Delta t)$ and $T^v(\Delta t)$ with rigorous convergence proofs. Our numerical results, however, suggest that the choice (3.14) gives a better approximation. If we want to determine whether the numerical solution at a certain point $x_0$ remains bounded at the numerical blow-up time as $\Delta t\rightarrow 0$, it seems that it will be easier for us to observe the boundedness of the numerical solutions if $\gamma$ is not greater than that given in (3.14). See Section 4.1. On the other hand, in the analysis of the blow-up rate for the numerical solutions, we need the assumption that $\gamma$ is not smaller than that given in (3.14). See Theorem 4.2. For the computation of $L^\sigma$-norm blow-up, our analysis again suggests that the numerical results coincide with that of the continuous problem if (3.14) is applied. All these results suggest (3.14) might be a suitable choice for the computation of blow-up solutions of (1.2).

Finally, we remark that, although we can prove the convergence of the numerical blow-up time for (1.2) (Theorem 3.2), we have no idea about the convergence order. This still remains open even in the case of a single equation (1.8). We thus leave this to future study.

Acknowledgments

The present paper is an outgrowth of the second author's master's thesis ^[25] under the guidance of the first author in National Chung Cheng University. The first author is supported by the grant MOST 109-2115-M-194 -001 and was partially supported by the Advanced Institute of Manufacturing with High-tech Innovations (AIM-HI) from The Featured Areas Research Center Program within the framework of the Higher Education Sprout Project by the Ministry of Education in Taiwan.

Conflict of interest

The authors declare that they have no conflict of interest.