Processing math: 69%
Research article Special Issues

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

  • 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

    Related Papers:

    [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+N1rur+vp,r(0,1),t>0vt=vrr+N1rvr+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)(Tt)p+1pq1u(t,)L<andsupt(0,T)(Tt)q+1pq1v(t,)L<. (1.3)

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

    u0+vp00,v0+uq00

    or

    max{p+1pq1,q+1pq1}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(Tt)p+1pq1andv(t,)L(Tt)q+1pq1 (1.4)

    namely,

    C1<maxx¯Ωu(t,x)(Tt)p+1pq1<C2andC3<maxx¯Ωv(t,x)(Tt)q+1pq1<C4, (1.5)

    for some positive constants Ci(i=1,2,3,4) if we assume

    u0+(1a)vp00andv0+(1a)uq00, (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,vr0 and ut,vt0.

    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,)LastT,

    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 suptTu(t,)Lσ1=,if σ1>N(pq1)2(p+1)lim suptTv(t,)Lσ2=,if σ2>N(pq1)2(q+1). (1.7)

    Later, Souplet [32] extended (1.7) for the solutions of (1.2) to the case of σ1=N(pq1)2(p+1) and σ2=N(pq1)2(q+1) under the assumptions ur,vr0 and ut,vt0.

    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 N2, 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+11p+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

    C01q+1uq+101p+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)astTO<,

    where

    TOu0dsG1(s)=v0dsG2(s). (2.4)

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

    {(Un)tUn+1UnΔt=(Vn)p,U0=u0(Vn)tVn+1VnΔ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(n0) 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(n0).

    Moreover, multiplying (Un+1)q and (Vn)p to the first and the second equations of (2.5) respectively, we have

    (Un+1Un)(Un+1)q=Δt(Vn)p(Un+1)q=(Vn+1Vn)(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]{|Unu(tn)|,|Vnv(tn)|}0asΔt0.

    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+11p+1(Vn)p+1C0(n0). (2.7)

    In particular, one has

    (Un)tG1(Un), (2.8)

    and

    (Vn)tG2(Vn+1). (2.9)

    Proof. By (2.6), one has for all n0

    [1q+1(Un+1)q+11p+1(Vn+1)p+1][1q+1(Un)q+11p+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+1Un(Un+1)q1q+1[(Un+1)q+1(Un)q+1]}+{Vn+1(Vn)p(Vn)p+11p+1[(Vn+1)p+1(Vn)p+1]}=qq+1(Un+1)q(Un+1Un)+1q+1Un[(Un+1)q(Un)q]1p+1Vn+1[(Vn+1)p(Vn)p]+pp+1(Vn)p(Vn+1Vn)qq+1(Un+1)q(Un+1Un)+qq+1(Un+1)q1Un(Un+1Un)pp+1(Vn)p1Vn+1(Vn+1Vn)+pp+1(Vn)p(Vn+1Vn)=qq+1(Un+1)q1(Un+1Un)2pp+1(Vn)p1(Vn+1Vn)20,

    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)tG1(Un+1), (2.10)

    and

    (Vn)tG2(Vn). (2.11)

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

    [(Un)t]p+1p[(Un1)t]p+1p=(Vn)p+1(Vn1)p+1(p+1)(Vn)p(VnVn1)=(p+1)(Un+1Un)(Un)q,

    which implies

    [(Un)t]p+1p[(U0)t]p+1p+(p+1)nk=0(Uk+1Uk)(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[(Vn1)t]q+1q=(Un+1)q+1(Un)q+1(q+1)(Un)q(Un+1Un)=(q+1)(VnVn1)(Vn)p,

    which implies

    [(Vn)t]q+1q[(V0)t]q+1q+(q+1)nk=1(VkVk1)(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+1Un+ΔtG1(Un)Un+ΔtG1(U0)U0+(n+1)ΔtG1(U0)asn,

    which, together with (2.5), also implies

    Vn+1=Vn+Δt(Un+1)qasn.

    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 limsH(s)=. Given any Δt>0, since Un,Vn as n, there exist nuΔt,nvΔtN such that

    ΔtH(UnuΔt1)<1andΔtH(UnuΔt)1. (2.12)

    and

    ΔtH(UnvΔt1)<1andΔtH(UnvΔt)1. (2.13)

    We define TuO(Δt)=ΔtnuΔt and TvO(Δt)=ΔtnvΔt and call them the numerical blow-up time of u and v respectively.

    Theorem 2.2. Assume that H satisfies

    ΔtlnGi(H1(1Δt))0asΔt0(i=1,2). (2.14)

    Then limΔt0TuO(Δt)=limΔt0TvO(Δ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Δt0TuO(Δt)TOandlim infΔt0TuO(Δt)TO.

    First, we consider the finite difference scheme

    Yn+1YnΔt=G1(Yn),Y0=U0.

    It is easy to show, by (2.8), that UnYn(n0). Let ˉnuΔtN be the positive integer such that

    ΔtH(YˉnuΔt1)<1andΔtH(YˉnuΔt)1.

    Then one has TuO(Δt)=ΔtnuΔtΔtˉnuΔt, which, together with [Theorem 2.1, [10]], gives

    TuO(Δt)ΔtˉnuΔtU0dsG1(s)H1(1Δt)dsG1(s)+Δt(1+lnG1(H1(1Δt))G1(U0)). (2.15)

    We thus have

    lim supΔt0TuO(Δt)U0dsG1(s)=TO.

    Here use has been made of (2.14).

    On the other hand, we assume that

    lim infΔt0TuO(Δt)TuO_<TO.

    Then there exists {τk} such that τk0 as k and

    TuO(τk)=τknuτk<TO+TuO_2<TO.

    Let {Un(τk)} be the solutions corresponding to τk. Then

    Unuτk(τk)H1(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Δt0TuO(Δ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

    ΔtlnH1(1Δt)0asΔt0. (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)TOC(Δt+(Δt)pq1γ(p+1)+Δtln|Δt|), (2.17)

    where C denotes a generic constant. The error bound

    TvO(Δt)TOC(Δt+(Δt)pq1γ(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)pq1γ(p+1)+Δtln|Δt|)

    and

    |TvO(Δt)TO|C(Δt+(Δt)pq1γ(q+1)+Δtln|Δt|),

    hold true. In fact, our computational results suggest that

    |TuO(Δt)TO|={O(Δtln|Δt|), if γpq1p+1O((Δt)pq1γ(p+1)) if γ>pq1p+1,

    and

    |TvO(Δt)TO|={O(Δtln|Δt|), if γpq1q+1O((Δt)pq1γ(q+1)) if γ>pq1q+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)230.1164774. Figures 13 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.

    Figure 1.  We choose Hu(s)=spq1p+1 and Hv(s)=spq1q+1 for our computation. (Left) Δt vs. TuO(Δt) and TvO(Δt). (Right) The solid lines are the log-log plot of the errors |TuO(Δt)TO| and |TvO(Δt)TO|, while the dashed line is the log-log plot of the function y=0.3x|logx|.
    Figure 2.  We choose Hu(s)=spq3p+1 and Hv(s)=spq3q+1 for our computation. (Left) Δt vs. TuO(Δt) and TvO(Δt). (Right) The solid lines are the log-log plot of the errors |TuO(Δt)TO| and |TvO(Δt)TO|, while the dashed line is the log-log plot of the function y=0.3x|logx|.
    Figure 3.  We choose Hu(s)=spq+1p+1 and Hv(s)=spq+1q+1 for our computation. (Left) Δt vs. TuO(Δt) and TvO(Δt). (Right) The solid lines are the log-log plot of the errors |TuO(Δt)TO| and |TvO(Δt)TO|, while the dashed line is a straight line with slope pq1pq+1.

    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 JN. Δ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(n0) are the temporal grid points. Let N0=N12, where N12 is the largest integer among those that are smaller or equal to N12. 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+10Un0Δt=N2Un12Un0(Δr)2+(Vn0)p, (3.1)
    Vn+10Vn0Δt=N2Vn12Vn0(Δr)2+(Un+10)q. (3.2)

    For j=1,,N0, we consider the following discretization:

    Un+1jUnjΔt=NUnj+12Unj+Unj1(Δr)2+(Vnj)p, (3.3)
    Vn+1jVnjΔt=NVnj+12Vnj+Vnj1(Δr)2+(Un+1j)q, (3.4)

    while, for j=N0+1,,J1, we consider

    Un+1jUnjΔt=Unj+12Unj+Unj1(Δr)2+N1rjUnj+1Unj12Δr+(Vnj)p, (3.5)
    Vn+1jVnjΔt=Vnj+12Vnj+Vnj1(Δr)2+N1rjVnj+1Vnj12Δ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)212N be fixed. Assume that u0(r),v0(r)0 for all r[0,1]. Then Unj,Vnj0 for all j=0,,J and n0.

    Proof. The proof is completed by induction. Assume that Unj,Vnj0 for all j=0,,J. Then it follows directly from (3.1)(3.3) and the assumption λ12N that Un+1j0, 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,,J1, one has by (3.5)

    Un+1j=(12λ)Unj+λ(1+N12j)Unj+1+λ(1N12j)Unj1+Δt(Vnj)p.

    Since jN0+1>N12, the coefficients on the right-hand side of the above equation are nonnegative, which implies Un+1j0. Then one has by (3.6)

    Vn+1j=(12λ)Vnj+λ(1+N12j)Vnj+1+λ(1N12j)Vnj1+Δt(Un+1j)q0.

    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

    max0jJ,tn[0,T0]|u(tn,rj)Unj|0asΔt0,

    and

    max0jJ,tn[0,T0]|v(tn,rj)Vnj|0asΔt0,

    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 UnjUnj+10 and VnjVnj+10 for all j=0,,J1 and n0.

    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 UnjUnj+10 and VnjVnj+10 for all j=0,,J1. We first show Un+1jUn+1j+1(j=0,,J1). To this end, it suffices to show Un+10Un+11,Un+1N0Un+1N0+1. For other js, the arguments given in Lemma 3.1 can be applied to derive Un+1jUn+1j+10.

    By (3.1)(3.3), one has

    Un+10=(12Nλ)Un0+2NλUn1+Δt(Vn0)p,

    and

    Un+11=(12Nλ)Un1+Nλ(Un0+Un2)+Δt(Vn1)p(1Nλ)Un1+NλUn0+Δt(Vn1)p,

    from which follow

    Un+10Un+11(13Nλ)(Un0Un1)+Δt[(Vn0)p(Vn1)p]0.

    For j=N0, one has by (3.3)(3.5)

    Un+1N0=(12Nλ)UnN0+Nλ(UnN01+UnN0+1)+Δt(VnN0)p(1Nλ)UnN0+NλUnN0+1+Δt(VnN0)p,

    and

    Un+1N0+1=(12λ)UnN0+1+λ(1N12j)UnN0+λ(1+N12j)UnN0+2+Δt(VnN0+1)p(1λ)UnN0+1+λUnN0+Δt(VnN0+1)p,

    which imply

    Un+1N0Un+1N0+1(1(N+1)λ)(UnN0UnN0+1)+Δt[(VnN0)p(VnN0+1)p]0.

    Now Vn+1jVn+1j+1>0 can be proved in a similar way and the proof is completed by induction.

    From now on, we assume that u0,v00 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+(1a)(V0j)p0 and ΔdV0j+(1a)(U1j)q0 for all j=0,,J1.

    Here, the operator Δd is defined by

    ΔdYj={N2Yj+2Yj+1(Δr)2,j=0NYj+12Yj+Yj1(Δr)2,1jN0Yj+12Yj+Yj1(Δr)2+N1jΔrYj+1Yj12Δr,N0+1jJ1 (3.8)

    Remark 3.3. Observe that the assumption (A) implies

    (U0j)t=ΔdU0j+(V0j)pa(V0j)p0,

    and

    (V0j)t=ΔdV0j+(U1j)qa(U1j)q0.

    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 0jJ1 and n0,

    Wnj=dUnj+(1a)(Vnj)pandZnj=dVnj+(1a)(Un+1j)q.

    Note that, by (3.1)–(3.6), Wnj and Znj can also be written as

    Wnj=(Unj)ta(Vnj)pandZnj=(Vnj)ta(Un+1j)q.

    Since UnJ=VnJ=0(n0), we set WnJ=ZnJ=0(n0).

    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,,J1 and n0,

    (Unj)ta(Vnj)pand(Vnj)ta(Un+1j)q. (3.10)

    That is, Wnj,Znj0 for all j=0,,J1 and n0. In particular, one has

    {Un+1UnΔt=Un+10Un0Δta(Vn0)pVn+1VnΔt=Vn+10Vn0Δta(Un+10)q, (3.11)

    where Un=maxj=0,,J|Unj| and Vn=maxj=0,,J|Vnj|.

    Proof. Assume that Wnj,Znj0(j=0,.,J1) holds for certain n. Observe that

    (Wnj)t=(dUnj)t+(1a)[(Vnj)p]t=d[(Unj)t]+(1a)[(Vnj)p]t,

    and

    dWnj=d[(Unj)t]ad[(Vnj)p],

    imply

    (Wnj)tdWnj=(1a)[(Vnj)p]t+ad[(Vnj)p].

    Since Znj0(j), one has by Lemma 3.1 Vn+1jVnj, which implies

    [(Vnj)p]tp(Vnj)p1(Vnj)t.

    On the other hand, one has by Lemma 3.2

    Δd[(Vnj)p]p(Vnj)p1ΔdVnj.

    Consequently,

    (Wnj)tdWnjp(Vnj)p1[(1a)(Vnj)t+aΔdVnj]=p(Vnj)p1[(Vnj)ta(Un+1j)q]=p(Vnj)p1Znj0.

    This, together with the assumption λ13N<12N, implies the nonnegativity of Wn+1j(0jJ1).

    Similarly, one can easily derive

    (Znj)tdZnj=(1a)[(Un+1j)q]t+ad[(Un+1j)q].

    By virtue of the fact Wn+1j0 and Lemma 3.1, 3.2, we have

    [(Un+1j)q]tq(Un+1j)q1(Un+1j)tandΔd[(Un+1j)q]q(Un+1j)q1ΔdUn+1j,

    from which follow

    (Znj)tdZnj=q(Un+1j)q1[(1a)(Un+1j)t+aΔdUn+1j]=q(Un+1j)q1[(Un+1j)ta(Vn+1j)q]=q(Un+1j)q1Wn+1j0.

    Again, the stability condition λ13N yields Zn+1j0. 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+10U0+(n+1)aΔt(V00)pasn,

    and

    Vn+1=Vn+10V0+(n+1)aΔt(U00)qasn.

    Let H:(0,)(0,) be a strictly increasing function satisfying limsH(s)=. Then for any given Δt>0, there exist positive integers nuΔt,nvΔtN such that

    ΔtH(UnuΔt1)<1andΔtH(UnuΔt)1, (3.12)

    and

    ΔtH(VnvΔt1)<1andΔtH(VnvΔt)1. (3.13)

    We define the numerical blow-up time for the solutions u and v of (1.2) by Tu(Δt)=ΔtnuΔt and Tv(Δt)=ΔtnvΔ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Δt0Tu(Δt)=limΔt0Tv(Δt)=T.

    Proof. We write down the prove of limΔt0Tu(Δt)=T. Convergence of Tv(Δt) can be shown in parallel.

    We complete the proof by showing

    Tu_lim infΔt0Tu(Δt)Tand¯Tulim supΔt0Tu(Δt)T.

    Assume that Tu_<T. Then there exists {(τi,hi)} satisfying 0<λ=τih2i13N, τi,hi0 as i and that

    Tu(τi)=τinuτi<Tu_+T2<T.

    Since, by (3.12), τiH(Unuτi)1, we have

    UnuτiH1(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<λ=τih2i13N, τi,hi0 as i and that

    Tu(τi)=τinuτ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

    Δtn0<TandUn0M.

    For simplicity, we denote the solutions of (3.1)–(3.7) corresponding to the grid sizes Δt,Δr by UnjUnj(Δt,Δr),VnjVnj(Δt,Δr). We remark that n0<nuΔt since

    Δtn0<T<¯Tu+T2<Tu(Δt)=ΔtnuΔt.

    Now we consider the finite difference system

    {An+1AnΔt=a(Bn)p,An0=MBn+1BnΔt=a(An+1)q,Bn0=Vn0

    By (3.11), it is easy to show

    UnAnandVnBn(nn0).

    Let R(Δt) be the numerical blow-up time for An. Namely, R(Δt)=Δt(mΔtn0), where mΔtN is determined by

    ΔtH(AmΔt1)<1andΔtH(AmΔt)1.

    Since H is strictly increasing, one has H(Un)H(An)(nn0), which implies

    Δt(nuΔtn0)Δt(mΔtn0)=R(Δt).

    Note that we have, by Theorem 2.2,

    R(Δt)MdsaG1(s)asΔt0.

    Since we can choose sufficiently large M and sufficiently small Δt such that

    MdsaG1(s)<¯TuT4and|R(Δt)MdsaG1(s)|<¯TuT4,

    it then follows

    Tu(Δt)=Δtn0+Δt(nuΔtn0)<T+R(Δt)<T+¯TuT2=¯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)=spq1p+1andHv(s)=spq1q+1 (3.14)
    Figure 4.  Δt vs. UnuΔt and VnvΔt. The numerical solution UnuΔt and VnvΔt become unbounded as Δt0.
    Figure 5.  The numerical blow-up time Tu(Δt) and Tv(Δt).

    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.8pq1p+1andHv(s)=s0.8pq1q+1 (3.15)
    Figure 6.  The numerical blow-up time Tu(Δt) and Tv(Δt).

    and

    Hu(s)=s1.2pq1p+1andHv(s)=s1.2pq1q+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 N2, 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 N2 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 Δt0. 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Δrx0<(jΔr+1)Δr. Then we compute Unj and Vnj at (Tu(Δt),rjΔr) and (Tv(Δt),rjΔr) respectively.

    ● Let Δt,Δr0 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,Δr0. See Figure 7. Similar results can be observed for x0=0.001 and x0=0.0001.

    Figure 7.  x0=0.01. (Left) The behavior of UnjΔr at Tu(Δt) as Δt0. (Right) The behavior of VnjΔr at Tv(Δt) as Δt0.

    We also compute the behaviors of Unj and Vnj at (Tu(Δt),rjΔr) and (Tv(Δt),rjΔr) respectively as Δt0 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 Δt0. 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 Δt0, they are always bounded by the ones computed with (3.14) or (3.15).

    Figure 8.  x0=0.01. (Left) The behavior of U at t=Tu(Δt) as Δt0. (Right) The behavior of V at t=Tv(Δt) as Δt0.

    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Δt0UnuΔtjΔrM<.

    Then

    lim suptTu(t,x0)<.

    (b) Assume that

    lim supΔt0VnvΔtjΔr<.

    Then

    lim suptTv(t,x0)<.

    Proof. We outline the proof for u. Assume on the contrary that lim suptTu(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εΔt1)<tεΔtnεΔ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],0jJ|Unju(tn,rj)|0asΔt0,

    it follows that, for sufficiently small Δt and Δr,

    |UnϵΔtjΔru(tϵ,x0)|<M2,

    which implies

    UnϵΔtjΔr>u(tϵ,x0)M2>2MM2=3M2.

    On the other hand, since limΔt0Tu(Δt)=T, one has, for sufficiently small Δt,

    ΔtnεΔt<T+tε2<Tu(Δt)=ΔtnuΔt.

    Now (3.10) yields

    UnuΔtjΔr>UnϵΔtjΔr>3M2,

    for all sufficiently small Δt. This contradicts the assumption lim supΔt0UnuΔ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 Δt0, it suggests the boundedness of the exact solution u(t,x0),v(t,x0) as tT. 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 γpq1p+1. Then there exist constants c1,c2>0 such that

    c1(Tu(Δt)tn)(Un0)pq1p+1c2(n=0,,nuΔt1). (4.2)

    (b) Let γpq1q+1. Then there exist constants c3,c4>0 such that

    c3(Tv(Δt)tn)(Vn0)pq1q+1c4(n=0,,nvΔt1). (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 γpq1p+1. There exist C,c>0 such that

    c(Un0)p(q+1)p+1(Un0)tC(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 Vk0Vk10(k1),

    (Vk0)p+1(Vk10)p+1(p+1)(Vk0)p(Vk0Vk10),

    which, together with (4.5), implies

    (Vk0)p+1(Vk10)p+1+(p+1)Uk+10Uk0aΔtΔt(Uk0)q(V00)p+1+p+1akl=1(Ul+10Ul0)(Ul0)q(V00)p+1+p+1aUk+10U10xqdx=(V00)p+1+p+1a(q+1)[(Uk+10)q+1(U10)q+1].

    We thus have

    (Un0)t(Vn0)p(Vn10+Δ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(0knuΔt1),

    and that

    qγqpq1p+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γ}pC(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,
    Figure 9.  (Left) The discrete L^3 -norm of U at t = T^u(\Delta t) as \Delta t\rightarrow 0 . (Right) The discrete L^4 -norm of U at t = T^u(\Delta t) as \Delta t\rightarrow 0 .
    Figure 10.  (Left) The discrete L^4 -norm of V at t = T^v(\Delta t) as \Delta t\rightarrow 0 . (Right) The discrete L^5 -norm of V at t = T^v(\Delta t) as \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.
  • This article has been cited by:

    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
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3171) PDF downloads(179) Cited by(2)

Figures and Tables

Figures(10)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog