
We study the finite difference approximation for axisymmetric solutions of a parabolic system with blow-up. A scheme with adaptive temporal increments is commonly used to compute an approximate blow-up time. There are, however, some limitations to reproduce the blow-up behaviors for such schemes. We thus use an algorithm, in which uniform temporal grids are used, for the computation of the blow-up time and blow-up behaviors. In addition to the convergence of the numerical blow-up time, we also study various blow-up behaviors numerically, including the blow-up set, blow-up rate and blow-up in Lσ-norm. Moreover, the relation between blow-up of the exact solution and that of the numerical solution is also analyzed and discussed.
Citation: Chien-Hong Cho, Ying-Jung Lu. On the numerical solutions for a parabolic system with blow-up[J]. AIMS Mathematics, 2021, 6(11): 11749-11777. doi: 10.3934/math.2021683
[1] | Jia Li, Changchun Bi . Existence and blowup of solutions for non-divergence polytropic variation-inequality in corn option trading. AIMS Mathematics, 2023, 8(7): 16748-16756. doi: 10.3934/math.2023856 |
[2] | Huafei Di, Yadong Shang, Jiali Yu . Blow-up analysis of a nonlinear pseudo-parabolic equation with memory term. AIMS Mathematics, 2020, 5(4): 3408-3422. doi: 10.3934/math.2020220 |
[3] | Yaxin Zhao, Xiulan Wu . Asymptotic behavior and blow-up of solutions for a nonlocal parabolic equation with a special diffusion process. AIMS Mathematics, 2024, 9(8): 22883-22909. doi: 10.3934/math.20241113 |
[4] | Zhiqiang Li . The finite time blow-up for Caputo-Hadamard fractional diffusion equation involving nonlinear memory. AIMS Mathematics, 2022, 7(7): 12913-12934. doi: 10.3934/math.2022715 |
[5] | Hao Zhang, Zejian Cui, Xiangyu Xiao . Decay estimate and blow-up for a fourth order parabolic equation modeling epitaxial thin film growth. AIMS Mathematics, 2023, 8(5): 11297-11311. doi: 10.3934/math.2023572 |
[6] | Sen Ming, Xiaodong Wang, Xiongmei Fan, Xiao Wu . Blow-up of solutions for coupled wave equations with damping terms and derivative nonlinearities. AIMS Mathematics, 2024, 9(10): 26854-26876. doi: 10.3934/math.20241307 |
[7] | Zhanwei Gou, Jincheng Shi . Blow-up phenomena and global existence for nonlinear parabolic problems under nonlinear boundary conditions. AIMS Mathematics, 2023, 8(5): 11822-11836. doi: 10.3934/math.2023598 |
[8] | W. Y. Chan . Blow-up for degenerate nonlinear parabolic problem. AIMS Mathematics, 2019, 4(5): 1488-1498. doi: 10.3934/math.2019.5.1488 |
[9] | Wai Yuen Chan . Simultaneous blow-up of the solution for a singular parabolic system with concentrated sources. AIMS Mathematics, 2024, 9(3): 6951-6963. doi: 10.3934/math.2024339 |
[10] | Mengyang Liang, Zhong Bo Fang, Su-Cheol Yi . Blow-up analysis for a reaction-diffusion equation with gradient absorption terms. AIMS Mathematics, 2021, 6(12): 13774-13796. doi: 10.3934/math.2021800 |
We study the finite difference approximation for axisymmetric solutions of a parabolic system with blow-up. A scheme with adaptive temporal increments is commonly used to compute an approximate blow-up time. There are, however, some limitations to reproduce the blow-up behaviors for such schemes. We thus use an algorithm, in which uniform temporal grids are used, for the computation of the blow-up time and blow-up behaviors. In addition to the convergence of the numerical blow-up time, we also study various blow-up behaviors numerically, including the blow-up set, blow-up rate and blow-up in Lσ-norm. Moreover, the relation between blow-up of the exact solution and that of the numerical solution is also analyzed and discussed.
We consider axisymmetric solutions for the semilinear parabolic system
{ut=△u+vp,x∈Ω,t>0vt=△v+uq,x∈Ω,t>0u(0,x)=u0(x),v(0,x)=v0(x),x∈Ωu(t,x)=v(t,x)=0,x∈∂Ω,t>0, | (1.1) |
where △ denotes the Laplace operator, p,q>1 and Ω⊆RN is the unit ball with the origin as its center. Since we consider the radially symmetric solutions, letting r=|x|, (1.1) reads as
{ut=urr+N−1rur+vp,r∈(0,1),t>0vt=vrr+N−1rvr+uq,r∈(0,1),t>0u(0,r)=u0(r),v(0,r)=v0(r),r∈(0,1)ur(t,0)=vr(t,0)=0,t>0u(t,1)=v(t,1)=0,t>0. | (1.2) |
Here, the initial data u0(r),v0(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 u0(r),v0(r)≥0 is decreasing in (0,1) and
supt∈(0,T)(T−t)p+1pq−1‖u(t,⋅)‖L∞<∞andsupt∈(0,T)(T−t)q+1pq−1‖v(t,⋅)‖L∞<∞. | (1.3) |
In fact, (1.3) is known to hold if one assumes that either
△u0+vp0≥0,△v0+uq0≥0 |
or
max{p+1pq−1,q+1pq−1}≥N2. |
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
‖u(t,⋅)‖L∞∼(T−t)−p+1pq−1and‖v(t,⋅)‖L∞∼(T−t)−q+1pq−1 | (1.4) |
namely,
C1<maxx∈¯Ωu(t,x)(T−t)p+1pq−1<C2andC3<maxx∈¯Ωv(t,x)(T−t)q+1pq−1<C4, | (1.5) |
for some positive constants Ci(i=1,2,3,4) if we assume
△u0+(1−a)vp0≥0and△v0+(1−a)uq0≥0, | (1.6) |
for certain a∈(0,1). In [32], (1.4) was also shown to be true for the solution of (1.2) under the assumptions ur,vr≤0 and ut,vt≥0.
Although the solution is known to become unbounded in a finite time T in the sense of L∞-norm, that is,
‖u(t,⋅)‖L∞,‖v(t,⋅)‖L∞→∞ast→T, |
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
{lim supt→T‖u(t,⋅)‖Lσ1=∞,if σ1>N(pq−1)2(p+1)lim supt→T‖v(t,⋅)‖Lσ2=∞,if σ2>N(pq−1)2(q+1). | (1.7) |
Later, Souplet [32] extended (1.7) for the solutions of (1.2) to the case of σ1=N(pq−1)2(p+1) and σ2=N(pq−1)2(q+1) under the assumptions ur,vr≤0 and ut,vt≥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
ut=△u+up(p>1) | (1.8) |
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≥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σ-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
{u′(t)=vp(t)v′(t)=uq(t)u(0)=u0>0,v(0)=v0>0,(pq>1), | (1.9) |
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σ-norm blow-ups for both u and v. Finally, the paper is ended up with a conclusion.
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 uq and vp, one has
uq(t)u′(t)=uq(t)vp(t)=vp(t)v′(t), |
which implies
(1q+1uq+1−1p+1vp+1)′(t)=0. |
It then follows
uq+1(t)=q+1p+1vp+1(t)+(q+1)C0andvp+1(t)=p+1q+1uq+1(t)−(p+1)C0, | (2.1) |
where
C0≡1q+1uq+10−1p+1vp+10. | (2.2) |
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
u′(t)=G1(u(t))andv′(t)=G2(v(t)), | (2.3) |
where
G1(s)=[p+1q+1sq+1−(p+1)C0]pp+1andG2(s)=[q+1p+1sp+1+(q+1)C0]qq+1. |
Since pq>1, it is not difficult to derive
u(t),v(t)→∞ast→TO<∞, |
where
TO≡∫∞u0dsG1(s)=∫∞v0dsG2(s). | (2.4) |
Now we consider the following finite difference scheme for (1.9):
{(Un)t≡Un+1−UnΔt=(Vn)p,U0=u0(Vn)t≡Vn+1−VnΔt=(Un+1)q,V0=v0. | (2.5) |
Here, we use (⋅)t to denote the forward difference operator. Δt is the grid size, tn=nΔt(n≥0) are the grid points and Un,Vn denote the approximations for u(tn),v(tn) respectively. Since U0=u0,V0=v0>0, one has
Un+1>Un>0,andVn+1>Vn>0(∀n≥0). |
Moreover, multiplying (Un+1)q and (Vn)p to the first and the second equations of (2.5) respectively, we have
(Un+1−Un)(Un+1)q=Δt(Vn)p(Un+1)q=(Vn+1−Vn)(Vn)p. | (2.6) |
Theorem 2.1. Let Un,Vn be the solution of (2.5). Let T0<TO be given arbitrarily. Then
maxtn∈[0,T0]{|Un−u(tn)|,|Vn−v(tn)|}→0asΔt→0. |
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
1q+1(Un)q+1−1p+1(Vn)p+1≤C0(∀n≥0). | (2.7) |
In particular, one has
(Un)t≥G1(Un), | (2.8) |
and
(Vn)t≤G2(Vn+1). | (2.9) |
Proof. By (2.6), one has for all n≥0
[1q+1(Un+1)q+1−1p+1(Vn+1)p+1]−[1q+1(Un)q+1−1p+1(Vn)p+1]=1q+1[(Un+1)q+1−(Un)q+1]−1p+1[(Vn+1)p+1−(Vn)p+1]=−{(Un+1)q+1−Un(Un+1)q−1q+1[(Un+1)q+1−(Un)q+1]}+{Vn+1(Vn)p−(Vn)p+1−1p+1[(Vn+1)p+1−(Vn)p+1]}=−qq+1(Un+1)q(Un+1−Un)+1q+1Un[(Un+1)q−(Un)q]−1p+1Vn+1[(Vn+1)p−(Vn)p]+pp+1(Vn)p(Vn+1−Vn)≤−qq+1(Un+1)q(Un+1−Un)+qq+1(Un+1)q−1Un(Un+1−Un)−pp+1(Vn)p−1Vn+1(Vn+1−Vn)+pp+1(Vn)p(Vn+1−Vn)=−qq+1(Un+1)q−1(Un+1−Un)2−pp+1(Vn)p−1(Vn+1−Vn)2≤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
(Un)t≤G1(Un+1), | (2.10) |
and
(Vn)t≥G2(Vn). | (2.11) |
Proof. First, one has, by (2.5),
[(Un)t]p+1p−[(Un−1)t]p+1p=(Vn)p+1−(Vn−1)p+1≤(p+1)(Vn)p(Vn−Vn−1)=(p+1)(Un+1−Un)(Un)q, |
which implies
[(Un)t]p+1p≤[(U0)t]p+1p+(p+1)n∑k=0(Uk+1−Uk)(Uk)q≤(V0)p+1+(p+1)∫Un+1U0sqds=(V0)p+1+p+1q+1[(Un+1)q+1−(U0)q+1]=G1(Un+1)p+1p. |
On the other hand, one has again by (2.5)
[(Vn)t]q+1q−[(Vn−1)t]q+1q=(Un+1)q+1−(Un)q+1≥(q+1)(Un)q(Un+1−Un)=(q+1)(Vn−Vn−1)(Vn)p, |
which implies
[(Vn)t]q+1q≥[(V0)t]q+1q+(q+1)n∑k=1(Vk−Vk−1)(Vk)p≥(U1)q+1+(q+1)∫VnV0spds≥(U0)q+1+q+1p+1[(Vn)p+1−(V0)p+1]=G2(Vn)q+1q. |
Observe that G1(U0)=(v0)p>0. By (2.8) and the monotonicity of G1 and {Un}, one has
Un+1≥Un+ΔtG1(Un)≥Un+ΔtG1(U0)≥U0+(n+1)ΔtG1(U0)→∞asn→∞, |
which, together with (2.5), also implies
Vn+1=Vn+Δt(Un+1)q→∞asn→∞. |
Now we use the algorithm proposed in [10] to define the numerical blow-up time. We import a strictly increasing function H:(0,∞)→(0,∞), which is to be determined, satisfying lims→∞H(s)=∞. Given any Δt>0, since Un,Vn→∞ as n→∞, there exist nuΔt,nvΔt∈N such that
Δt⋅H(UnuΔt−1)<1andΔt⋅H(UnuΔt)≥1. | (2.12) |
and
Δt⋅H(UnvΔt−1)<1andΔt⋅H(UnvΔt)≥1. | (2.13) |
We define TuO(Δt)=Δt⋅nuΔt and TvO(Δt)=Δt⋅nvΔt and call them the numerical blow-up time of u and v respectively.
Theorem 2.2. Assume that H satisfies
Δt⋅lnGi(H−1(1Δt))→0asΔt→0(i=1,2). | (2.14) |
Then limΔt→0TuO(Δt)=limΔt→0TvO(Δt)=TO.
Proof. By virtue of (2.8) and (2.11), the proof is similar to [Theorem 3.2, [10]]. We outline the proof for TuO(Δt) for the readers' convenience. The proof for TvO(Δt) can be done in exactly the same way.
We complete the proof by showing that
lim supΔt→0TuO(Δt)≤TOandlim infΔt→0TuO(Δt)≥TO. |
First, we consider the finite difference scheme
Yn+1−YnΔt=G1(Yn),Y0=U0. |
It is easy to show, by (2.8), that Un≥Yn(∀n≥0). Let ˉnuΔt∈N be the positive integer such that
Δt⋅H(YˉnuΔt−1)<1andΔt⋅H(YˉnuΔt)≥1. |
Then one has TuO(Δt)=Δt⋅nuΔt≤Δt⋅ˉnuΔt, which, together with [Theorem 2.1, [10]], gives
TuO(Δt)≤Δt⋅ˉnuΔt≤∫∞U0dsG1(s)−∫∞H−1(1Δt)dsG1(s)+Δt(1+lnG1(H−1(1Δt))G1(U0)). | (2.15) |
We thus have
lim supΔt→0TuO(Δt)≤∫∞U0dsG1(s)=TO. |
Here use has been made of (2.14).
On the other hand, we assume that
lim infΔt→0TuO(Δt)≡TuO_<TO. |
Then there exists {τk} such that τk→0 as k→∞ and
TuO(τk)=τk⋅nuτk<TO+TuO_2<TO. |
Let {Un(τk)} be the solutions corresponding to τk. Then
Unuτk(τk)≥H−1(1τk)→∞ask→∞, |
while the exact solution u(t) remains bounded in [0,TO+TuO_2], which is a contradiction to Theorem 2.1. Therefore,
lim infΔt→0TuO(Δt)≥TO, |
and we are done.
Remark 2.1. Observe that G1(s)∼sp(q+1)p+1 and G2(s)∼sq(p+1)q+1. Therefore, (2.14) can be replaced by
Δt⋅lnH−1(1Δt)→0asΔt→0. | (2.16) |
It is easy to see that H(s)=sγ(γ>0) satisfies (2.16) and thus can be used to compute an approximate blow-up time. In addition, by (2.15), one has
TuO(Δt)−TO≤C(Δt+(Δt)pq−1γ(p+1)+Δtln|Δt|), | (2.17) |
where C denotes a generic constant. The error bound
TvO(Δt)−TO≤C(Δt+(Δt)pq−1γ(q+1)+Δtln|Δt|), | (2.18) |
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
|TuO(Δt)−TO|≤C(Δt+(Δt)pq−1γ(p+1)+Δtln|Δt|) |
and
|TvO(Δt)−TO|≤C(Δt+(Δt)pq−1γ(q+1)+Δtln|Δt|), |
hold true. In fact, our computational results suggest that
|TuO(Δt)−TO|={O(Δtln|Δt|), if γ≤pq−1p+1O((Δt)pq−1γ(p+1)) if γ>pq−1p+1, |
and
|TvO(Δt)−TO|={O(Δtln|Δt|), if γ≤pq−1q+1O((Δt)pq−1γ(q+1)) if γ>pq−1q+1, |
In the following example, we set p=3,q=2. The initial data is given by v0=3, u0=(q+1p+1vp+10)1q+1 so that C0 in (2.2) vanishes. Then the blow-up time TO can be computed explicitly by TO=∫∞3ds(34s4)23≈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.
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∈N. Δr=1/J is the spatial grid size and rj=jΔr,(j=0,⋯,J) are the spatial grid points. Δt>0 is the temporal grid size and tn=nΔt(n≥0) are the temporal grid points. Let N0=⌊N−12⌋, where ⌊N−12⌋ is the largest integer among those that are smaller or equal to N−12. Unj,Vnj are the approximations of u,v at (tn,rj) respectively. At the origin, we discretize as follows due to the boundary condition ur(t,0)=vr(t,0)=0:
Un+10−Un0Δt=N2Un1−2Un0(Δr)2+(Vn0)p, | (3.1) |
Vn+10−Vn0Δt=N2Vn1−2Vn0(Δr)2+(Un+10)q. | (3.2) |
For j=1,⋯,N0, we consider the following discretization:
Un+1j−UnjΔt=NUnj+1−2Unj+Unj−1(Δr)2+(Vnj)p, | (3.3) |
Vn+1j−VnjΔt=NVnj+1−2Vnj+Vnj−1(Δr)2+(Un+1j)q, | (3.4) |
while, for j=N0+1,⋯,J−1, we consider
Un+1j−UnjΔt=Unj+1−2Unj+Unj−1(Δr)2+N−1rjUnj+1−Unj−12Δr+(Vnj)p, | (3.5) |
Vn+1j−VnjΔt=Vnj+1−2Vnj+Vnj−1(Δr)2+N−1rjVnj+1−Vnj−12Δr+(Un+1j)q. | (3.6) |
The boundary conditions are given by
UnJ=VnJ=0, | (3.7) |
and the initial data are defined by U0j=u0(rj) and V0j=v0(rj) for j=0,⋯,J. Note that (3.1)(3.2) are nothing but approximations for axisymmetric solutions of
ut=△u+vpandvt=△v+uq |
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 Unj,Vnj be the solution of (3.1)–(3.7). Let λ≡Δt(Δr)2≤12N be fixed. Assume that u0(r),v0(r)≥0 for all r∈[0,1]. Then Unj,Vnj≥0 for all j=0,⋯,J and n≥0.
Proof. The proof is completed by induction. Assume that Unj,Vnj≥0 for all j=0,⋯,J. Then it follows directly from (3.1)(3.3) and the assumption λ≤12N that Un+1j≥0, for j=0,⋯,N0, which, together with (3.2)(3.4), subsequently imply the nonnegativity of Vn+1j(j=1,⋯,N0). For j=N0+1,⋯,J−1, one has by (3.5)
Un+1j=(1−2λ)Unj+λ(1+N−12j)Unj+1+λ(1−N−12j)Unj−1+Δt(Vnj)p. |
Since j≥N0+1>N−12, the coefficients on the right-hand side of the above equation are nonnegative, which implies Un+1j≥0. Then one has by (3.6)
Vn+1j=(1−2λ)Vnj+λ(1+N−12j)Vnj+1+λ(1−N−12j)Vnj−1+Δt(Un+1j)q≥0. |
Thus, we are done.
Theorem 3.1. Let T be the blow-up time of (1.2) and T0<T be given. Let λ≤12N be fixed. Then
max0≤j≤J,tn∈[0,T0]|u(tn,rj)−Unj|→0asΔt→0, |
and
max0≤j≤J,tn∈[0,T0]|v(tn,rj)−Vnj|→0asΔt→0, |
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 Unj,Vnj be the solution of (3.1)–(3.7). Let λ≤13N be fixed. Assume that u0(r),v0(r)≥0 and u0(r),v0(r) are monotonically decreasing in (0,1). Then Unj≥Unj+1≥0 and Vnj≥Vnj+1≥0 for all j=0,⋯,J−1 and n≥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 Unj≥Unj+1≥0 and Vnj≥Vnj+1≥0 for all j=0,⋯,J−1. We first show Un+1j≥Un+1j+1(j=0,⋯,J−1). To this end, it suffices to show Un+10≥Un+11,Un+1N0≥Un+1N0+1. For other j′s, the arguments given in Lemma 3.1 can be applied to derive Un+1j−Un+1j+1≥0.
By (3.1)(3.3), one has
Un+10=(1−2Nλ)Un0+2NλUn1+Δt(Vn0)p, |
and
Un+11=(1−2Nλ)Un1+Nλ(Un0+Un2)+Δt(Vn1)p≤(1−Nλ)Un1+NλUn0+Δt(Vn1)p, |
from which follow
Un+10−Un+11≥(1−3Nλ)(Un0−Un1)+Δt[(Vn0)p−(Vn1)p]≥0. |
For j=N0, one has by (3.3)(3.5)
Un+1N0=(1−2Nλ)UnN0+Nλ(UnN0−1+UnN0+1)+Δt(VnN0)p≥(1−Nλ)UnN0+NλUnN0+1+Δt(VnN0)p, |
and
Un+1N0+1=(1−2λ)UnN0+1+λ(1−N−12j)UnN0+λ(1+N−12j)UnN0+2+Δt(VnN0+1)p≤(1−λ)UnN0+1+λUnN0+Δt(VnN0+1)p, |
which imply
Un+1N0−Un+1N0+1≥(1−(N+1)λ)(UnN0−UnN0+1)+Δt[(VnN0)p−(VnN0+1)p]≥0. |
Now Vn+1j≥Vn+1j+1>0 can be proved in a similar way and the proof is completed by induction.
From now on, we assume that u0,v0≥0 and u0,v0 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)
(A) There exists a∈(0,1) such that ΔdU0j+(1−a)(V0j)p≥0 and ΔdV0j+(1−a)(U1j)q≥0 for all j=0,⋯,J−1.
Here, the operator Δd is defined by
ΔdYj={N−2Yj+2Yj+1(Δr)2,j=0NYj+1−2Yj+Yj−1(Δr)2,1≤j≤N0Yj+1−2Yj+Yj−1(Δr)2+N−1jΔrYj+1−Yj−12Δr,N0+1≤j≤J−1 | (3.8) |
Remark 3.3. Observe that the assumption (A) implies
(U0j)t=ΔdU0j+(V0j)p≥a(V0j)p≥0, |
and
(V0j)t=ΔdV0j+(U1j)q≥a(U1j)q≥0. |
This can be regarded as a discrete analogue for
ut(0,r),vt(0,r)≥0(0<r<1), | (3.9) |
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≤j≤J−1 and n≥0,
Wnj=△dUnj+(1−a)(Vnj)pandZnj=△dVnj+(1−a)(Un+1j)q. |
Note that, by (3.1)–(3.6), Wnj and Znj can also be written as
Wnj=(Unj)t−a(Vnj)pandZnj=(Vnj)t−a(Un+1j)q. |
Since UnJ=VnJ=0(∀n≥0), we set WnJ=ZnJ=0(∀n≥0).
Lemma 3.3. Let Unj,Vnj be the solutions of (3.1)-(3.7). Let λ≤13N. Suppose that there exists a∈(0,1) such that (A) holds. Then, for all j=0,⋯,J−1 and n≥0,
(Unj)t≥a(Vnj)pand(Vnj)t≥a(Un+1j)q. | (3.10) |
That is, Wnj,Znj≥0 for all j=0,⋯,J−1 and n≥0. In particular, one has
{‖Un+1‖∞−‖Un‖∞Δt=Un+10−Un0Δt≥a(Vn0)p‖Vn+1‖∞−‖Vn‖∞Δt=Vn+10−Vn0Δt≥a(Un+10)q, | (3.11) |
where ‖Un‖∞=maxj=0,⋯,J|Unj| and ‖Vn‖∞=maxj=0,⋯,J|Vnj|.
Proof. Assume that Wnj,Znj≥0(j=0,⋯.,J−1) holds for certain n. Observe that
(Wnj)t=(△dUnj)t+(1−a)[(Vnj)p]t=△d[(Unj)t]+(1−a)[(Vnj)p]t, |
and
△dWnj=△d[(Unj)t]−a△d[(Vnj)p], |
imply
(Wnj)t−△dWnj=(1−a)[(Vnj)p]t+a△d[(Vnj)p]. |
Since Znj≥0(∀j), one has by Lemma 3.1 Vn+1j≥Vnj, which implies
[(Vnj)p]t≥p(Vnj)p−1(Vnj)t. |
On the other hand, one has by Lemma 3.2
Δd[(Vnj)p]≥p(Vnj)p−1ΔdVnj. |
Consequently,
(Wnj)t−△dWnj≥p(Vnj)p−1[(1−a)(Vnj)t+aΔdVnj]=p(Vnj)p−1[(Vnj)t−a(Un+1j)q]=p(Vnj)p−1Znj≥0. |
This, together with the assumption λ≤13N<12N, implies the nonnegativity of Wn+1j(∀0≤j≤J−1).
Similarly, one can easily derive
(Znj)t−△dZnj=(1−a)[(Un+1j)q]t+a△d[(Un+1j)q]. |
By virtue of the fact Wn+1j≥0 and Lemma 3.1, 3.2, we have
[(Un+1j)q]t≥q(Un+1j)q−1(Un+1j)tandΔd[(Un+1j)q]≥q(Un+1j)q−1ΔdUn+1j, |
from which follow
(Znj)t−△dZnj=q(Un+1j)q−1[(1−a)(Un+1j)t+aΔdUn+1j]=q(Un+1j)q−1[(Un+1j)t−a(Vn+1j)q]=q(Un+1j)q−1Wn+1j≥0. |
Again, the stability condition λ≤13N yields Zn+1j≥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
‖Un+1‖∞=Un+10≥U0+(n+1)aΔt(V00)p→∞asn→∞, |
and
‖Vn+1‖∞=Vn+10≥V0+(n+1)aΔt(U00)q→∞asn→∞. |
Let H:(0,∞)⟼(0,∞) be a strictly increasing function satisfying lims→∞H(s)=∞. Then for any given Δt>0, there exist positive integers nuΔt,nvΔt∈N such that
Δt⋅H(‖UnuΔt−1‖∞)<1andΔt⋅H(‖UnuΔt‖∞)≥1, | (3.12) |
and
Δt⋅H(‖VnvΔt−1‖∞)<1andΔt⋅H(‖VnvΔt‖∞)≥1. | (3.13) |
We define the numerical blow-up time for the solutions u and v of (1.2) by Tu(Δt)=Δt⋅nuΔt and Tv(Δt)=Δt⋅nvΔt.
Theorem 3.2. Let T be the blow-up time for the solution of (1.2) and let Unj,Vnj be the solutions of (3.1)–(3.7). Let λ≤13N be fixed. Assume that the assumption (A) holds true and that the function H satisfies (2.16). Then we have limΔt→0Tu(Δt)=limΔt→0Tv(Δt)=T.
Proof. We write down the prove of limΔt→0Tu(Δt)=T. Convergence of Tv(Δt) can be shown in parallel.
We complete the proof by showing
Tu_≡lim infΔt→0Tu(Δt)≥Tand¯Tu≡lim supΔt→0Tu(Δt)≤T. |
Assume that Tu_<T. Then there exists {(τi,hi)} satisfying 0<λ=τih2i≤13N, τi,hi→0 as i→∞ and that
Tu(τi)=τi⋅nuτi<Tu_+T2<T. |
Since, by (3.12), τi⋅H(‖Unuτi‖∞)≥1, we have
‖Unuτi‖∞≥H−1(1τi)→∞asi→∞, |
while the exact solution u(t,x) remains bounded in [0,Tu_+T2]. This is a contradiction to Theorem 3.1. We thus have Tu_≥T.
On the other hand, we assume that ¯Tu>T. Then there exists {(τi,hi)} satisfying 0<λ=τih2i≤13N, τi,hi→0 as i→∞ and that
Tu(τi)=τi⋅nuτi>¯Tu+T2>T. |
Given any M>0, Theorem 3.1 guarantees the existence of sufficiently small (Δt,Δr)∈{(τi,hi)} and a positive integer n0, depending on Δt and Δr, such that
Δt⋅n0<Tand‖Un0‖∞≥M. |
For simplicity, we denote the solutions of (3.1)–(3.7) corresponding to the grid sizes Δt,Δr by Unj≡Unj(Δt,Δr),Vnj≡Vnj(Δt,Δr). We remark that n0<nuΔt since
Δt⋅n0<T<¯Tu+T2<Tu(Δt)=Δt⋅nuΔt. |
Now we consider the finite difference system
{An+1−AnΔt=a(Bn)p,An0=MBn+1−BnΔt=a(An+1)q,Bn0=‖Vn0‖∞ |
By (3.11), it is easy to show
‖Un‖∞≥Anand‖Vn‖∞≥Bn(∀n≥n0). |
Let R(Δt) be the numerical blow-up time for An. Namely, R(Δt)=Δt⋅(mΔt−n0), where mΔt∈N is determined by
Δt⋅H(AmΔt−1)<1andΔt⋅H(AmΔt)≥1. |
Since H is strictly increasing, one has H(‖Un‖∞)≥H(An)(∀n≥n0), which implies
Δt⋅(nuΔt−n0)≤Δt⋅(mΔt−n0)=R(Δt). |
Note that we have, by Theorem 2.2,
R(Δt)→∫∞MdsaG1(s)asΔt→0. |
Since we can choose sufficiently large M and sufficiently small Δt such that
∫∞MdsaG1(s)<¯Tu−T4and|R(Δt)−∫∞MdsaG1(s)|<¯Tu−T4, |
it then follows
Tu(Δt)=Δt⋅n0+Δt(nuΔt−n0)<T+R(Δt)<T+¯Tu−T2=¯Tu+T2<Tu(Δt), |
which is a contradiction. Thus, we are done.
In the following example, we set N=5,p=3,q=2,λ=1/(3N). The initial data are given by u0(r)=150(1+cos(πr)),v0(r)=100(1+cos(πr)). In Figures 4 and 5, we choose
Hu(s)=spq−1p+1andHv(s)=spq−1q+1 | (3.14) |
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∈[0,1],
ut(t,r)∼G1(u(t,r))andvt(t,r)∼G2(v(t,r)). |
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
Hu(s)=s0.8pq−1p+1andHv(s)=s0.8pq−1q+1 | (3.15) |
and
Hu(s)=s1.2pq−1p+1andHv(s)=s1.2pq−1q+1. | (3.16) |
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.
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 u0,v0 are nonnegative, monotonically decreasing and satisfy assumption (A).
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
ut=uxx+up(p>1). | (4.1) |
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≥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∈(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≥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 x0∈(0,1], the numerical solution remains bounded at x0 at the numerical blow-up time as Δt→0. We proceed our computation as follows:
● Given Δt,Δr>0 satisfying λ≤13N, we first compute the numerical blow-up time Tu(Δt) and Tv(Δt).
● Let jΔr satisfy jΔrΔr≤x0<(jΔr+1)Δr. Then we compute Unj and Vnj at (Tu(Δt),rjΔr) and (Tv(Δt),rjΔr) respectively.
● Let Δt,Δr→0 to see whether or not the numerical solution remains bounded at rjΔr.
Here, we choose (3.14) for the computation of Tu(Δt) and Tv(Δt) respectively. In the following example, we set N=5, (p,q)=(3,2), and u0(r)=150(1+cos(πr)),v0(r)=100(1+cos(πr)). Let x0=0.01. Our computational results show that the numerical solution remains bounded at r=x0 as Δt,Δr→0. See Figure 7. Similar results can be observed for x0=0.001 and x0=0.0001.
We also compute the behaviors of Unj and Vnj at (Tu(Δt),rjΔr) and (Tv(Δt),rjΔr) respectively as Δt→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 x0=0.01 decreases as Δt→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 Δt→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. (a) Assume that
lim supΔt→0UnuΔtjΔr≡M<∞. |
Then
lim supt→Tu(t,x0)<∞. |
(b) Assume that
lim supΔt→0VnvΔtjΔr<∞. |
Then
lim supt→Tv(t,x0)<∞. |
Proof. We outline the proof for u. Assume on the contrary that lim supt→Tu(t,x0)=∞. Then given ε>0, there exists tε∈(T−ε,T) such that u(tε,x0)>2M. Observe that
0<T−ε<tε<T+tε2<T. |
For small Δt>0, let nεΔt be the positive integer satisfying that
T−ε<Δt⋅(nεΔt−1)<tε≤Δt⋅nεΔt<T+tε2. |
Since u remains smooth on [0,T+tε2]×[0,1] and, by Theorem 3.1, we also have
maxtn∈[0,T+tε2],0≤j≤J|Unj−u(tn,rj)|→0asΔt→0, |
it follows that, for sufficiently small Δt and Δr,
|UnϵΔtjΔr−u(tϵ,x0)|<M2, |
which implies
UnϵΔtjΔr>u(tϵ,x0)−M2>2M−M2=3M2. |
On the other hand, since limΔt→0Tu(Δt)=T, one has, for sufficiently small Δt,
Δt⋅nεΔt<T+tε2<Tu(Δt)=Δt⋅nuΔt. |
Now (3.10) yields
UnuΔtjΔr>UnϵΔtjΔr>3M2, |
for all sufficiently small Δt. This contradicts the assumption lim supΔt→0UnuΔtjΔr=M.
Thanks to this theorem, if one can know the boundedness of the numerical solution at certain spatial point x0 at the numerical blow-up time as Δt→0, it suggests the boundedness of the exact solution u(t,x0),v(t,x0) as t→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).
In this subsection, we show that our numerical solutions satisfy a discrete analogue of (1.5). Put H(s)=sγ(γ>0).
Theorem 4.2. (a) Let γ≥pq−1p+1. Then there exist constants c1,c2>0 such that
c1≤(Tu(Δt)−tn)(Un0)pq−1p+1≤c2(∀n=0,⋯,nuΔt−1). | (4.2) |
(b) Let γ≥pq−1q+1. Then there exist constants c3,c4>0 such that
c3≤(Tv(Δt)−tn)(Vn0)pq−1q+1≤c4(∀n=0,⋯,nvΔt−1). | (4.3) |
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 γ≥pq−1p+1. There exist C,c>0 such that
c(Un0)p(q+1)p+1≤(Un0)t≤C(Un0)p(q+1)p+1. | (4.4) |
Proof. Note that, by (3.11) and (3.1)(3.2),
{a(Vn0)p<(Un0)t≤(Vn0)pa(Un+10)q<(Vn0)t≤(Un+10)q. | (4.5) |
Since Vk0≥Vk−10(∀k≥1),
(Vk0)p+1−(Vk−10)p+1≤(p+1)(Vk0)p(Vk0−Vk−10), |
which, together with (4.5), implies
(Vk0)p+1≤(Vk−10)p+1+(p+1)Uk+10−Uk0aΔtΔt(Uk0)q≤(V00)p+1+p+1ak∑l=1(Ul+10−Ul0)(Ul0)q≤(V00)p+1+p+1a∫Uk+10U10xqdx=(V00)p+1+p+1a(q+1)[(Uk+10)q+1−(U10)q+1]. |
We thus have
(Un0)t≤(Vn0)p≤(Vn−10+Δt(Un0)q)p≤{((V00)p+1+p+1a(q+1)[(Un0)q+1−(U10)q+1])1p+1+Δt(Un0)q}p. |
Note that (3.12) implies
Δt⋅(Uk0)γ<1(∀0≤k≤nuΔt−1), |
and that
q−γ≤q−pq−1p+1=q+1p+1. |
Since {Un0} is increasing in n and Un0→∞ as n→∞, it then follows
(Un0)t≤{((V00)p+1+p+1a(q+1)[(Un0)q+1−(U10)q+1])1p+1+(Un0)q−γ}p≤C(Un0)p(q+1)p+1, |
for some positive constant C .
On the other hand, note that
\left(V^k_0\right)^{p+1}-\left(V^{k-1}_0\right)^{p+1}\geq(p+1)\left(V^{k-1}_0\right)^{p}\left(V^k_0-V^{k-1}_0\right), |
which, together with (4.5), implies
\begin{eqnarray*} \left(V^k_0\right)^{p+1}& \geq &\left(V^{k-1}_0\right)^{p+1}+(p+1)\frac{U^{k}_0-U^{k-1}_0}{\Delta t}a\Delta t\left(U^k_0\right)^q\\ & \geq &\left(V^{0}_0\right)^{p+1}+a(p+1)\sum\limits_{l = 0}^{k-1}\left(U^{l+1}_0-U^l_0\right)\left(U^{l+1}_0\right)^q\\ & \geq &\left(V^{0}_0\right)^{p+1}+a(p+1)\int_{U^0_0}^{U^{k}_0}x^qdx\\ & = &\left(V^{0}_0\right)^{p+1}+\frac{a(p+1)}{q+1}\left[\left(U^{k}_0\right)^{q+1}-\left(U^0_0\right)^{q+1}\right]. \end{eqnarray*} |
It then follows that there exists c > 0 such that
\begin{equation*} \label{blowuprategeq} \left(U^n_0\right)_t\geq \left\{\left(V^{0}_0\right)^{p+1}+\frac{a(p+1)}{q+1}\left[\left(U^{n}_0\right)^{q+1}-\left(U^0_0\right)^{q+1}\right]\right\}^\frac{p}{p+1}\geq c\left(U^{n}_0\right)^\frac{p(q+1)}{p+1}. \end{equation*} |
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
f_k(z) = z+C\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}z^\frac{p(q+1)}{p+1}-z_k, \quad\mathit{\text{for}} \ z\in[0, 1]. |
Then \{z_k\} is decreasing. Moreover, we have
\begin{equation} \left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-k}_0\geq z_k\quad(\forall k\geq 1). \end{equation} | (4.6) |
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
U^{n^u_{\Delta t}-k}_0\leq U^{n^u_{\Delta t}-(k+1)}_0+C\Delta t\left(U^{n^u_{\Delta t}-(k+1)}_0\right)^\frac{p(q+1)}{p+1}. |
Then it follows
\begin{eqnarray*} \label{ykineq} z_k& \leq &\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0+C\left(\Delta t\right)^{1+\frac{1}{\gamma}}\left(U^{n^u_{\Delta t}-(k+1)}_0\right)^\frac{p(q+1)}{p+1}\\ & = &\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0+C\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}\left(\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0\right)^\frac{p(q+1)}{p+1}, \end{eqnarray*} |
or equivalently,
f_k\left(\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0\right)\geq 0. |
By the monotonicity of f_k and the fact f_k(z_{k+1}) = 0 , we thus have
\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0\geq z_{k+1}. |
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
g_k(z) = z+c\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}z^\frac{p(q+1)}{p+1}-y_k, \quad\mathit{\text{for}} \ z\in[0, 1). |
Then \{y_k\} is decreasing. Moreover, we have
\begin{equation} \left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-k}_0\leq y_k\quad(\forall k\geq 1). \end{equation} | (4.7) |
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
U^{n^u_{\Delta t}-k}_0\geq U^{n^u_{\Delta t}-(k+1)}_0+c\Delta t\left(U^{n^u_{\Delta t}-(k+1)}_0\right)^\frac{p(q+1)}{p+1}, |
from which follows
y_k \geq \left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0+c\left(\Delta t\right)^{1+\frac{1}{\gamma}}\left(U^{n^u_{\Delta t}-(k+1)}_0\right)^\frac{p(q+1)}{p+1}\\ = \left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0+c\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}\left(\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0\right)^\frac{p(q+1)}{p+1}, |
or equivalently,
g_k\left(\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0\right)\leq 0. |
By the monotonicity of f_k and the fact f_k(y_{k+1}) = 0 , we thus have
\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-(k+1)}_0\leq y_{k+1}, |
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
\frac{\left(z_{n-1}-z_n\right)z_{n}^{-\frac{p(q+1)}{p+1}}}{\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}} = C, |
and
\frac{\left(y_{n-1}-y_n\right)y_{n}^{-\frac{p(q+1)}{p+1}}}{\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}} = c, |
which imply
n\sim\frac{1}{\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}}\sum\limits_{k = 1}^n\left(z_{n-1}-z_n\right)z_{n}^{-\frac{p(q+1)}{p+1}}\sim\frac{1}{\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}}\int_{z_0}^{z_n}s^{-\frac{p(q+1)}{p+1}}ds, |
and
n\sim\frac{1}{\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}}\sum\limits_{k = 2}^n\left(y_{n-1}-y_n\right)y_{n}^{-\frac{p(q+1)}{p+1}}\sim\frac{1}{\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}}\int_{y_1}^{y_n}s^{-\frac{p(q+1)}{p+1}}ds. |
It is not difficult to derive that there exist c_1, c_2 > 0 such that
c_1\leq n\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}z_n^\frac{pq-1}{p+1}\leq n\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}y_n^\frac{pq-1}{p+1}\leq c_2. |
Let k = n^u_{\Delta t}-n . Then one has
(T^u(\Delta t)-t_n)\left(U^n_0\right)^\frac{pq-1}{p+1} = k\left(\Delta t\right)\left(U^n_0\right)^\frac{pq-1}{p+1} = k\left(\Delta t\right)^{1-\frac{pq-1}{\gamma(p+1)}}\left[\left(\Delta t\right)^\frac{1}{\gamma}U^{n^u_{\Delta t}-k}_0\right]^\frac{pq-1}{p+1}, |
which, together with (4.6)(4.7) and the above inequality, implies (4.2).
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
\begin{equation} \|U^n\|_\sigma = \left\{ \begin{array}{ll} \max\limits_{j = 0, \cdots, J}|U^n_j| & , \;\text{if }p = \infty\\ \left(\sum\limits_{j = 0}^{J-1}r_{j+1}^{N-1}\Delta r|U^n_j|^\sigma\right)^\frac{1}{\sigma} & , \;\text{if }1\leq p < \infty \end{array} \right.. \end{equation} | (4.8) |
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
\|U^{n^u_{\Delta t}}\|_4, \;\|V^{n^v_{\Delta t}}\|_5\rightarrow \infty\quad\text{as}\quad \Delta t\rightarrow 0, |
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} ,
\|U^{n^u_{\Delta t}}\|_\sigma\rightarrow \infty\quad\mathit{\text{as}}\quad \Delta t\rightarrow 0. |
\rm(b) Let H^v(s) = s^{\gamma_2}\; (\gamma_2 > 0) . Then it holds that, for all \sigma > \frac{N \gamma_2}{2} ,
\|V^{n^v_{\Delta t}}\|_\sigma\rightarrow \infty\quad\mathit{\text{as}}\quad \Delta t\rightarrow 0. |
Proof. By (4.8) and (3.12) we have
\|U^{n^u_{\Delta t}}\|_\sigma\geq\left(\Delta r\right)^\frac{N}{\sigma}\|U^{n^u_{\Delta t}}\|_\infty\geq\lambda^{-\frac{N}{2\sigma}}\left(\Delta t\right)^{\frac{N}{2\sigma}-\frac{1}{\gamma_1}}\rightarrow \infty\quad\text{as}\quad\Delta t\rightarrow 0, |
provided that \sigma > \frac{N\gamma_1}{2} . Similarly, we have by (3.13)
\|V^{n^v_{\Delta t}}\|_\sigma\geq\left(\Delta r\right)^\frac{N}{\sigma}\|V^{n^v_{\Delta t}}\|_\infty\geq\lambda^{-\frac{N}{2\sigma}}\left(\Delta t\right)^{\frac{N}{2\sigma}-\frac{1}{\gamma_2}}\rightarrow \infty\quad\text{as}\quad\Delta t\rightarrow 0, |
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
\limsup\limits_{\Delta t\rightarrow 0}\max\limits_{0\leq n\leq n^u_{\Delta t}}\|U^n\|_\sigma < \infty. |
Then
\lim\limits_{t\rightarrow T}\|u(t, \cdot)\|_{L^\sigma} < \infty. |
\rm(b) Suppose that
\limsup\limits_{\Delta t\rightarrow 0}\max\limits_{0\leq n\leq n^v_{\Delta t}}\|V^n\|_\sigma < \infty. |
Then
\lim\limits_{t\rightarrow T}\|v(t, \cdot)\|_{L^\sigma} < \infty. |
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.
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
T^u(\Delta t), \;T^v(\Delta t)\rightarrow T\quad\text{as}\quad\Delta t\rightarrow 0. |
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.
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.
The authors declare that they have no conflict of interest.
[1] |
J. Abia, J. C. López-Marcos, J. Martínez, The Euler method in the numerical integration of reaction-diffusion problems with blow-up, Appl. Numer. Math., 38 (2001), 287–313. doi: 10.1016/S0168-9274(01)00035-6
![]() |
[2] |
C. J. Budd, W. Huang, R. D. Russel, Moving mesh methods for problems with blow-up, J. SIAM J. Sci. Comput., 17 (1996), 305–327. doi: 10.1137/S1064827594272025
![]() |
[3] |
C. J. Budd, O. Koch, L. Taghizadeh, E. B. Weinmüller, Asymptotic properties of the space-time adaptive numerical solution of a nonlinear heat equation, Calcolo, 55 (2018), 43. doi: 10.1007/s10092-018-0286-z
![]() |
[4] |
M. Burger, J. A. Carrillo, M. T. Wolfram, A mixed finite element method for nonlinear diffusion equations, Kinet. Relat. Models, 3 (2010), 59–83. doi: 10.3934/krm.2010.3.59
![]() |
[5] |
G. Caristi, E. Mitidieri, Blow-up estimates of positive solutions of a parabolic system, J. Diff. Eqn., 113 (1994), 265–271. doi: 10.1006/jdeq.1994.1124
![]() |
[6] | Y. G. Chen, Asymptotic behaviours of blowing-up solutions for finite difference analogue of u_t = u_xx+u^{1+\alpha}, J. Fac. Sci. Univ. Tokyo Sect. IA, 33 (1986), 541–574. |
[7] |
Y. G. Chen, Blow-up solutions to a finite difference analogue of u_t = u_xx+u^{1+\alpha} in N-dimensional ball, Hokkaido Math. J., 21 (1992), 447–474. doi: 10.14492/hokmj/1381413684
![]() |
[8] | M. Chlebík, M. Fila, From critical exponents to blow-up rates for parabolic problems, Rend. Mat. Appl., 19 (1999), 449–470. |
[9] |
C. H. Cho, On the finite difference approximation for blow-up solutions of the porous medium equation with a source, Appl. Numer. Math., 65 (2013), 1–26. doi: 10.1016/j.apnum.2012.11.001
![]() |
[10] |
C. H. Cho, On the computation of the numerical blow-up time, Japan J. Indust. Appl. Math., 30 (2013), 331–349. doi: 10.1007/s13160-013-0101-9
![]() |
[11] |
C. H. Cho, A numerical algorithm for blow-up problems revisited, Numer. Algor., 75 (2017), 675–697. doi: 10.1007/s11075-016-0216-6
![]() |
[12] |
C. H. Cho, On the computation for blow-up solutions of the nonlinear wave equation, Numerische Mathematik, 138 (2018), 537–556. doi: 10.1007/s00211-017-0919-1
![]() |
[13] |
C. H. Cho, S. Hamada, H. Okamoto, On the finite difference approximation for a parabolic blow-up problem, Japan J. Indust. Appl. Math., 24 (2007), 105–134. doi: 10.1007/BF03167510
![]() |
[14] | C. H. Cho, H. Okamoto, Further remarks on asymptotic behavior of the numerical solutions of parabolic blow-up problems, Methods Appl. Anal., 14 (2007), 213–226. |
[15] |
C. H. Cho, H. Okamoto, Finite difference schemes for an axisymmetric nonlinear heat equation with blow-up, Electronic Trans. Numer. Anal., 52 (2020), 391–415. doi: 10.1553/etna_vol52s391
![]() |
[16] |
K. Deng, Blow-up rates for parabolic systems, Z. Angew. Math. Phys., 47 (1996), 132–143. doi: 10.1007/BF00917578
![]() |
[17] |
A. Friedman, B. McLeod, Blow-up of positive solutions of semilinear heat equations, Indiana Uni. Math. J., 34 (1985), 425–447. doi: 10.1512/iumj.1985.34.34025
![]() |
[18] |
M. Fila, P. Souplet, The blow-up rate for semilinear parabolic problems on general domains, NoDEA Nonlinear Differ. Equ. Appl., 8 (2001), 473–480. doi: 10.1007/PL00001459
![]() |
[19] | A. Friedman, Y. Giga, A single point blow-up for solutions of semilinear parabolic systems, J. Fac. Sci. Univ. Tokyo Sect. IA Math., 34 (1987), 65–79. |
[20] | V. A. Galaktionov, S. P. Kurdyumov, A. A. Samarskii, A parabolic system of quasilinear equations I, Differential Equations, 19 (1983), 2123–2140. |
[21] | V. A. Galaktionov, S. P. Kurdyumov, A. A. Samarskii, A parabolic system of quasilinear equations II, Differential Equations, 21 (1985), 1544–1559. |
[22] |
P. Groisman, Totally discrete explicit and semi-implicit Euler methods for a blow-up problem in several space dimensions, Computing, 76 (2006), 325–352. doi: 10.1007/s00607-005-0136-0
![]() |
[23] | W. Huang, R. D. Russel, Adaptive moving mesh methods, Springer, New York, 2010. |
[24] |
W. Huang, J. Ma, R. D. Russel, A study of moving mesh PDE methods for numerical simulation of blowup in reaction diffusion equations, J. Comput. Phys., 227 (2008), 6532–6552. doi: 10.1016/j.jcp.2008.03.024
![]() |
[25] | Y. J. Lu, On a finite difference scheme for blow-up solutions of a semilinear parabolic system, Master's Thesis in National Chung Cheng University, 2018. |
[26] |
T. Nakagawa, Blowing up of a finite difference solution to u_t = u_xx+u^2, Appl. Math. Optim., 2 (1975), 337–350. doi: 10.1007/BF01448176
![]() |
[27] | T. Nakagawa, T. Ushijima, Numerical analysis of the semi-linear equation of blow-up type, Publications mathématiques et informatique de Rennes, S5 (1976), 1–24. |
[28] |
T. Nakanishi, N. Saito, Finite element method for radially symmetric solution of a multidimensional semilinear heat equation, Japan J. Indust. Appl. Math., 37 (2020), 165–191. doi: 10.1007/s13160-019-00393-z
![]() |
[29] |
P. Quittner, P. Souplet, Global existence from single-component Lp estimates in a semilinear reaction-diffusion system, Proc. Amer. Math. Soc., 130 (2002), 2719–2724. doi: 10.1090/S0002-9939-02-06453-5
![]() |
[30] |
W. Ren, X. P. Wang, An iterative grid redistribution method for singular problems in multiple dimensions, J. Comput. Phys., 159 (2000), 246–273. doi: 10.1006/jcph.2000.6435
![]() |
[31] |
N. Saito, Error analysis of a conservative finite-element approximation for the Keller-Segel system of chemotaxis, Commun. Pure Appl. Anal., 11 (2012), 339–364. doi: 10.3934/cpaa.2012.11.339
![]() |
[32] | P. Souplet, Single-point blow-up for a semilinear parabolic system, J. Eur. Math. Soc., 11 (2009), 169–188. |
[33] | G. Viglialoro, On the blow-up time of a parabolic system with damping terms, Comptes Rendus de L'Academie Bulgare des Sciences, 67 (2014), 1223–1232. |
[34] | F. Weissler, Single point blowup of semilinear initial value problems, J. Diff. Eqn., 55 (1984), 202–224. |
1. | Hui Guo, Xueting Liang, Yang Yang, Provable convergence of blow-up time of numerical approximations for a class of convection-diffusion equations, 2022, 466, 00219991, 111421, 10.1016/j.jcp.2022.111421 | |
2. | J.I. Ramos, Carmen María García López, Blow-up phenomena in a one-dimensional, bidirectional model equation for small-amplitude waves, 2025, 0961-5539, 10.1108/HFF-07-2024-0552 |