The phenomenon of one-point blow-up is a typical example in the numerical analysis of blow-up problems

25  Download (0)

Full text




Abstract.We study finite difference schemes for axisymmetric blow-up solutions of a nonlinear heat equation in higher spatial dimensions. The phenomenology of blow-up in higher-dimensional space is much more complex than that in one space dimension. To obtain a more complete picture for such phenomena from computational results, it is useful to know the technical details of the numerical schemes for higher spatial dimensions. Since first-order differentiation appears in the differential equation, we pay special attention to it. A sufficient condition for stability is derived. In addition to the convergence of the numerical blow-up time, certain blow-up behaviors, such as blow-up sets and blow-up in theLp-norm, are taken into consideration. It is sometimes experienced that a certain property of solutions of a partial differential equation may be lost by a faithfully constructed convergent numerical scheme.

The phenomenon of one-point blow-up is a typical example in the numerical analysis of blow-up problems. We prove that our scheme can preserve such a property. It is also remarkable that theLp-norm(1p <∞)of the solution of the nonlinear heat equation may blow up simultaneously with theL-norm or remains bounded in[0, T), whereTdenotes the blow-up time of theL-norm. We propose a systematic way to compute numerical evidence of theLp-norm blow-up. The computational results are also analyzed. Moreover, we prove an abstract theorem which shows the relationship between the numericalLp-norm blow-up and the exactLp-norm blow-up. Numerical examples for higher-dimensional blow-up solutions are presented and discussed.

Key words.blow-up, finite difference method, nonlinear heat equation,Lp-norm blow-up AMS subject classifications.65M06, 65M12

1. Introduction. LetΩbe the unit ball inRN with the origin as its center. We consider inΩthe axisymmetric solutions of a nonlinear heat equation of Fujita type, i.e.,

(1.1) ut=4u+u1+α,

where4 denotes the Laplace operator andα >0 is a given constant. We consider only Dirichlet boundary conditions in the present paper. Since the solution is assumed to be axisymmetric, our target is the following initial and boundary value problem:


r ur+u1+α (0≤r <1,0< t), (1.2)

ur= 0 (r= 0),


u= 0 (r= 1),


u(0, r) =u0(r) (0≤r≤1).


Here, we suppose that the initial datau0(r)are continuous in[0,1],u0(r)≥0everywhere, andu0(1) = 0. It was proved (see [5,14]) that the solutions of (1.2)–(1.5) may become infinite in finite time. This phenomenon is called blow-up, and the time at which the solution becomes infinite is called the blow-up time.

Friedman & McLeod [14] and Mueller & Weissler [25] proved that if the initial datau0

are monotonically decreasing inr, then the solution blows up only atr= 0. Later, Chen [5]

proved a similar proposition for a more general boundary condition with a milder assumption.

Numerical experiments tell us that a blow-upin generaloccurs at a single point, no matter

Received January 13, 2020. Accepted July 7, 2020. Published online on August 25, 2020. Recommended by Bill Layton.

Department of Mathematics and Advanced Institute of Manufacturing with High-tech Innovations, National Chung Cheng University, Chia-Yi, 621, Taiwan, R.O.C. (

Department of Mathematics, Gakushuin University, Tokyo, 171–8588, Japan.



where the blow-up occurs. But proving this statement is rather difficult. In this regard, [7] is worthy of notice since it proves that the blow-up set consists of finite points. However, this result seems to be proved only in one dimension, namely, only whenN= 1.

It is also remarkable that the solutions of (1.2)–(1.5) may become infinite in a finite time Tin the sense that

ku(t,·)kL → ∞ as t→T,

while theLp-norms(1≤p <∞)may either become infinite ast→T or remain bounded in [0, T). In fact, Friedman & McLeod [14] showed that

lim inf

t→T ku(t,·)kLp=∞ if p > N α 2 , (1.6)

lim sup


ku(t,·)kLp<∞ if p < N α 2 . (1.7)

Weissler [30] showed that (1.6) can be extended to the critical casep = N α2 under some restrictions on the initial datau0(x).

We would like to examine similar problems in a finite difference scheme. Specifically, we would like to construct a finite difference scheme that reproduces the blow-up characteristics described above. It is sometimes experienced that a certain property of solutions of a partial differential equation may be lost by a faithfully constructed convergent numerical solution.

Therefore, we consider it important to derive a numerical scheme that preserves such a property.

We assume thatN ≥2. In the case ofN= 1, we proposed in [8,11] a set of fairly general schemes for computing blow-up solutions. The purpose of the present paper is to compute blow-up solutions in higher dimensions. In two respects this extension is not straightforward.

First, the phenomenology of blow-up in higher dimensions is considerably more complex than that in one dimension. In fact, if we consider a higher-dimensional space, then many curious or bizarre phenomena are known to exist. One of them is the recovery of smoothness after a blow-up, which is prohibited in low dimensions, but is known to occur if the spatial dimension is large enough. See [13,22,23,24]. So far we know no successful numerical treatment for reproducing such phenomena. To reproduce them, we believe that we should know the technical details of the numerical schemes for higher dimensions, which were not treated in [8,9,10,11,12]. This is one of our motivations for the present paper. Secondly, in computing axisymmetric solutions, we must deal with a term involving first-order derivatives.

We need special care for the term in a finite difference scheme.

The number of blow-up points of numerical solutions of (1.2)–(1.5) was studied by Chen [6]. He assumed monotonicity of the initial datau0(r), in which case one-point blow-up is known to occur in (1.2). Nevertheless, he proved that the solution of a finite difference scheme may blow up at more than one point. A similar phenomenon was also proved for a numerical scheme of (1.1) with zero Dirichlet boundary conditions by Groisman [17]. Both of their results suggest that it is difficult for a scheme with a uniform spatial mesh to completely reflect the blow-up phenomena. Besides the number of blow-up points for the numerical solutions, it is worth mentioning that Groisman also proved that the blow-up rate for his numerical solution is completely determined byα, while there exist solutions that blow up in different rates for the same nonlinearityu1+αin the continuous problem [19,21]. We will visit this problem in a brief discussion in Section6.

The rest of this paper is organized as follows. In Section2, we briefly review two numerical algorithms for the computation of blow-up solutions. In Section3, we first recall a pioneering work by Chen [6], in which a finite difference scheme was considered for the computation of blow-up solutions of (1.2)–(1.5). Then we propose a class of finite difference


schemes and derive a sufficient condition for stability. The convergence of the numerical solution is also proved in the same section. In Section4, we focus on the behavior of the numerical solution, including the numerical blow-up time, the numerical blow-up set, and the Lp-norm blow-up. In Section5, we are concerned with the phenomenon of the recovery of smoothness after a blow-up in higher space dimensions. Namely, we try to compute what is called the minimal solution. Finally, the paper ends with concluding remarks.

After having finished the present paper, we learned about the paper [27] by Nakanishi and Saito, where the authors considered blow-up problems of the same kind with the finite element method. In [27], they showed convergence of their scheme and the numerical blow-up time.

An analysis of the asymptotic behavior of their numerical solution such as in Section4below, however, seems to be out of their scope.

2. Numerical methods for blow-up problems. In this section, we review two algo- rithms for the computation of blow-up solutions: adaptive and uniform time meshes.

2.1. Adaptive time mesh. Nakagawa [26] considered the one-dimensional semilinear heat equation




ut=uxx+u2 (0< x <1, t >0), u(0, x) =u0(x)≥0 (0≤x≤1), u(t,0) =u(t,1) = 0 (t >0),

and proposed the following numerical scheme with a uniform spatial mesh:







=Uj+1n −2Ujn+Uj−1n

(∆x)2 + (Ujn)2 (j= 1, . . . , J−1, n≥0), Uj0=u0(xj) (j= 1, . . . , J−1),

U0n =UJn = 0 (n≥0),

whose temporal grid sizes are defined adaptively by

(2.2) ∆tn=τ·min

1, c



Here,J ∈N,∆x= 1J is the spatial grid size,xj =j∆x(0 ≤j ≤J)are the spatial grid points,τ andcare prescribed parameters, andλ= (∆x)τ 212 is fixed. Moreover,t0 = 0 andtn = tn−1+ ∆tn−1 (n ≥ 1)denote the temporal grid points, and Ujn denotes the approximation of the exact solutionuof (2.1) at(tn, xj). The discreteLp-norm is defined by

kUnkp =



 max

j=1,...,J−1|Ujn| ifp=∞,






if1≤p <∞.

The temporal grid size∆tn in (2.2) is defined at each time step according to the discrete L2-norm of the numerical solution. As a matter of fact, the discreteL2-norm can be replaced by a certain measure of the numerical solution determined a priori; see for instance [11,17].

A mesh defined in this sense will be called here and hereafter an adaptive mesh. The meaning of (2.2) is that we choose∆tnuniformly whileUnis not so large, and an adaptive mesh begins to be in effect onceUnbecomes large. The author in [26] defined the numerical blow-up time asT(τ,∆x) =P

n=0∆tn. Then he proved that the numerical solutionUjnconverges to the


exact solutionuand that the numerical blow-up timeT(τ,∆x)is finite and converges to the exact blow-up time ofuas∆x→0.

Nakagawa’s method seems to be the first successful attempt to the computation of the blow-up time at which theL-norm of the solution becomes infinite. Later, his idea was generalized in [1,4,9,11,17,29], in which the asymptotic behavior of the numerical solution was also analyzed. One of our purposes in the present paper is to extend these results to higher dimensions.

2.2. Uniform time mesh. An adaptive temporal grid strategy seems to give a good approximation for blow-up problems, however, we cannot tell from the computational results whether or not theLp-norms (1≤p <∞) blow up simultaneously with theL-norm. In fact, the unboundedness of the discreteL-norm is equivalent to the unboundedness of other discreteLp-norms if time-independent spatial grids are applied. To reproduce phenomena such as (1.6), (1.7), a different approach for the computation of blow-up solutions is necessary.

For this purpose, we use the idea proposed in [8] and consider finite difference schemes whose temporal grid size is given uniformly(∆tn = ∆t for all n≥0)for the computation of a blow-up in theLp-norm.

For convenience, we illustrate the idea given in [8] by the following nonlinear ODE problem

u0(t) =G(u), u(0)>0, whose solution blows up atTO≡R

u(0) ds

G(s). Here, we assume for simplicity that the nonlinear termG(u)satisfiesG, G0, G00>0andR

u(0) ds G(s)<∞.

We now consider the forward Euler scheme

(2.3) Un+1−Un

∆t =G(Un), U0=u(0).

Note that the numerical solutions can be computed for any timetn=n∆teven though the exact solutionushows a blown up in finite timeTO. This implies that our computation should be stopped at a certain finite step. To determine such a step, we impose a nonnegative, strictly increasing functionH satisfying lim

s→∞H(s) =∞. One can easily derive by (2.3) and the monotonicity ofGthat for any given∆t,Unis monotonically increasing innand

Un→ ∞ as n→ ∞.

Accordingly there should exist a stepn∆t∈Ndepending on∆tsuch that

∆t·H Un∆t−1

<1 and ∆t·H(Un∆t)≥1.

We then defineTO(∆t) = ∆t·n∆t, which we call the numerical blow-up time.

THEOREM2.1.It holds that

− Z

H−1(∆t1 ) ds

G(s) ≤TO(∆t)−TO≤ − Z

H−1(∆t1 ) ds

G(s)+ ∆t 1 + logG H−1 ∆t1 G(U0)

! , whereH−1denotes the inverse function ofH. In particular, we have lim

∆t→0TO(∆t) =TOif the functionHsatisfies

(2.4) ∆t·log





→0 as ∆t→0.


For the proof, see [8].

REMARK 2.2. Although the computation is stopped at a finite step n∆t where the numerical solution is still finite, the finite-time blow-up is reproduced in the sense that the numerical solution atTO(∆t)becomes unbounded as∆ttends to0,

Un∆t ≥H−1 1


→ ∞ as ∆t→0,

and that the numerical blow-up timeTO(∆t)converges to the exact blow-up timeTOas∆t tends to0,

TO(∆t)→TO as ∆t→0.

Although we have explained the idea via an ODE for the sake of convenience, it can be applied to PDEs, too, as will be shown in Theorem4.13below.

REMARK2.3. We do not use an adaptive spatial mesh in this paper, not because we are uninterested in it. The reason is that our primary goal is to propose a numerical scheme that is supported by a convergence proof. Adaptive spatial meshes are very effective in numerically resolving the concentration of singularities of a solution of PDEs. There may be many references to such schemes, but we cite here only [28]. However, for those schemes it is difficult to prove their convergence mathematically, as far as we know.

REMARK2.4. IfG(u)is a power ofu, sayG(u) =ul(l >1), then (2.4) can be replaced by

(2.5) ∆t·log




→0 as ∆t→0.

In this case, it is easy to see thatH(s) =sq(∀q >0)satisfies (2.5) and thus can be used to compute an approximate blow-up time.

3. Finite difference scheme. We begin this section by recalling the pioneering work by Chen [6], where he considered a scheme for (1.2)–(1.5) in which the nonlinear term is evaluated explicitly and the diffusion term and the first-order derivative term are discretized implicitly by central differences. His scheme is

(3.1) Ujn+1−Ujn

∆tn = Uj+1n+1−2Ujn+1+Uj−1n+1

(∆r)2 +N−1 j∆r


2∆r + Ujn1+α , whereUjndenotes an approximation foru(tn, j∆r). Here,∆ris the spatial grid size, and the temporal increment∆tnis defined adaptively by


1, c kUnkαp

(1≤p≤ ∞),

whereτandcare prescribed parameters andkUnkpdenotes the discreteLp-norm:

(3.2) kUnkp =








∆r((j+ 1)∆r)N−1|Ujn|p


(1≤p <∞), max

j=0,...,J|Ujn| (p=∞).


To discretize the boundary conditionur(t,0) = 0, he proceeded as follows: Note that (1.3) implies that

(3.3) lim


ur(t, r) r = lim


ur(t, r)−ur(t,0)

r =urr(t,0).

Accordingly, one may assume the following equation:

(3.4) U0n+1−U0n

∆tn =N−2U0n+1+ 2U1n+1

(∆r)2 + (U0n)1+α.

This is nothing but an approximation of axisymmetric solutions of ut = 4u+u1+α by the central difference atr= 0. However, the scheme (3.1), together with (3.4), is unstable ifN > 3, as is proved in [6]. In order to overcome this difficulty, Chen [6] proposed an intermediate scheme between (3.1) and (3.4):




(∆r)2 + (Ujn)1+α for1≤j < N0:=

N+ 1 2

, and (3.1) forN0≤j < J. Here,N+1


denotes the largest integer≤ N2+1. He proved that his scheme is stable for allN ≥2and converges in the orderO (∆r)2

. He also analyzed the numerical blow-up set in the case thatu0(r)is nonnegative and decreasing inr∈[0,1].

He showed that his numerical solution blows up at not only the maximum point but also the point neighboring to it forα≤1, while the theoretical result tells us that the solution blows up only atr= 0, namely, the maximum point. A similar result was also mentioned by Groisman in his very interesting work [17] on (1.1). He discretized the partial differential equation by the so-called method of lines to obtain a system of ordinary differential equations. His spatial discretization may be quite general and grids which are not equidistant may be included. But once the spatial grids are set, they are fixed throughout the time evolution.

Both Chen’s and Groisman’s papers proved, under certain assumptions, that the maxima of their numerical solutions are notisolatedblow-up points in the case ofα≤1. Consequently, the question arises as to whether or not we can construct a convergent scheme that faithfully reproduces the one-point blow-up phenomenon for (1.2)–(1.5).

In view of these observations in the previous studies, we propose in this paper a new scheme by which we can reproduce some of the blow-up properties more faithfully than Chen’s scheme.

We divide[0,1]intoJ equal subintervals. Letrj =j∆r (0 ≤j ≤J)be the spatial grid points, where∆r= 1/J. SetJ0= max



] . Forj= 0andj ≥J0+ 1, we consider an explicit version of Chen’s scheme, that is,

(3.5) U0n+1−U0n

∆tn =N−2U0n+ 2U1n

(∆r)2 + (U0n)1+α and



= Uj+1n −2Ujn+Uj−1n

(∆r)2 +N−1 j∆r

Uj+1n −Uj−1n

2∆r + Ujn1+α , J0+ 1≤j≤J−1.


For1≤j ≤J0, we introduce a new parameterσand approximateurin the following way:


r ur(t, r) =N−1 r

σu(t, r)−u(t, r−∆r)

∆r + (1−σ)u(t, r+ ∆r)−u(t, r)



2r urr(t, r) +O((∆r)2).


Note that, for smallr,

ur(t, r)

r −urr(t, r) =O r2 . Thus, we have the following approximation

1 +(1−2σ)∆r 2r

N−1 r ur(t, r)

= N−1 r

σu(t, r)−u(t, r−∆r)

∆r + (1−σ)u(t, r+ ∆r)−u(t, r)


+O(r2) +O (∆r)2 .

Here, use has been made of the fact thaturrr(t,0) = 0.σis a parameter to be chosen. This observation leads us to consider the following scheme: For1≤j≤J0,



=Uj−1n −2Ujn+Uj+1n (∆r)2

+ 2j

2j+ 1−2σ N−1



∆r + (1−σ)Uj+1n −Ujn


+ Ujn1+α . (3.7)

Note that the truncation error for the scheme (3.7) is of order(∆r)2for allσ. Ifσ = 1/2, then (3.7) is nothing but (3.6). We would like to know which value ofσis suited to approximate the blow-up solutions.

We now setUj0=u0(rj)≥0, for0≤j≤J, and

(3.8) UJn = 0 (∀n≥0).

Here∆tn>0is the temporal grid size,t0= 0, andtn =tn−1+ ∆tnare the temporal grid points. The temporal increment can either be defined adaptively by

(3.9) ∆tn=τ·min

1, c


(γ >0),

or be given uniformly by

(3.10) ∆tn=τ (∀n≥0).

In these equations,τ >0is a prescribed small parameter, andkUnkpis defined as in (3.2).

REMARK 3.1. Forj = 0, besides (3.5), it is also natural to consider the following candidates

(3.11) U0n=U1n or U0n= (4U1n−U2n)/3.

In any of the cases of (3.5) and (3.11), smoothness is implicitly assumed. Since smoothness is lost at the blow-up point, it is not clear which is the best. We again note that the theoretical results tell us that the solution blows up at and only atr= 0ifu0(r)is nonnegative and de- creasing inr∈[0,1]. In view of this,U0n=U1nwould not be suitable for our purpose, since it automatically implies blow-ups at two points. We therefore exclude this from our consideration.


A similar situation also occurs for the second choice in (3.11) sinceU1n/U0n≥3/4 (∀n≥0) if nonnegative solutions are considered. Accordingly, we use (3.5) in what follows.

REMARK3.2. Our scheme is not a mere explicit version of Chen’s scheme except in the case thatσ = 1/2andN = 2. Nor is it a generalization of the scheme of Chen. It is true that our scheme is more complicated than Chen’s. However, our scheme has the merit that, in addition to stability and convergence, more blow-up properties are preserved. See the discussions in Section4.

Now we investigate the stability for the scheme (3.5)–(3.8). For the sake of convenience, we setλ=(∆r)τ 2 andλn=(∆r)∆tn2 (∀n≥0). Note that now (3.5)–(3.7) can be rewritten as (3.12) U0n+1= (1−2N λn)U0n+ 2N λnU1n+ ∆tn(U0n)1+α

for1≤j ≤J0, Ujn+1=

1−2λn−2(N−1)(1−2σ) 2j+ 1−2σ λn


1− 2σ(N−1) 2j+ 1−2σ


1 + 2(1−σ)(N−1) 2j+ 1−2σ

Uj+1n + ∆tn(Ujn)1+α


and, forJ0+ 1≤j≤J−1,

Ujn+1= (1−2λn)Ujnn

1−N−1 2j


1 +N−1 2j

Uj+1n + ∆tn(Ujn)1+α.


LEMMA3.3.Suppose thatu0(r)is nonnegative. Letλ≤ 2N1 be fixed andσ≤ 2N3 . Let {Ujn}be a solution of (3.5)–(3.8). Then,Ujn ≥0, for alln≥0and0≤j≤J.

Proof. The proof is carried out by induction. Assume thatUjn ≥0for all0 ≤j ≤J. Forj = 0andJ0+ 1 ≤ j ≤ J −1,Ujn+1 ≥ 0follows directly from (3.12), (3.14), and the assumptionλ≤ 2N1 . In the case1≤j ≤J0, it suffices to show that the coefficients at the right-hand side of (3.13) are nonnegative. For convenience, we set the coefficients at the right-hand side of (3.13) to

anj = 1−2λn−2(N−1)(1−2σ) 2j+ 1−2σ λn, bnjn

1− 2σ(N−1) 2j+ 1−2σ

, cnjn

1 + 2(1−σ)(N−1) 2j+ 1−2σ

. Sinceσ≤2N3 and1≤j, it is easy to see thatcnj ≥0for alln≥0andj. Observe also that bnj >0is obvious ifσ <0and that forσ≥0, we have

1− 2σ(N−1)

2j+ 1−2σ ≥1−2σ(N−1)

3−2σ = 3−2N σ 3−2σ ≥0,

which gives nonnegativity of bnj. Now it remains to show the nonnegativity of anj. For


2≤σ≤ 2N3 , we have

2(N−1)(1−2σ) 2j+ 1−2σ ≤0,


whence we obtainanj ≥1−2λn≥0. In the case ofσ < 12, note that 2(N−1)(1−2σ)

2j+ 1−2σ <2(N−1), which implies

anj >1−2λn−2(N−1)λn = 1−2N λn≥0.

Thus, we have proved the nonnegativity ofanj,bnj,cnj, and we are done.

LEMMA 3.4. Suppose thatu0(r)is nonnegative and monotonically decreasing. Let λ≤ 3N1 andσ ≤ 2N3 . Let{Ujn}be a solution of (3.5)–(3.8). Then,Ujn is monotonically decreasing inj. Namely,Ujn ≥Uj+1n ≥0, for alln≥0and0≤j≤J−1.

Proof.We prove this by induction. Assume thatUjn ≥Uj+1n ≥0for a certainnand all 0≤j≤J−1. Note that forj= 0,

U0n+1−U1n+1= (1−2N λn)U0n+ 2N λnU1n−(an1U1n+bn1U0n+cn1U2n) + ∆tn


≥(1−2N λn−bn1)U0n−(an1+cn1−2N λn)U1n. Since

1−2N λn−bn1 = 1−2N λn−3−2N σ 3−2σ λn


an1+cn1−2N λn= 1−λn+2(N−1)σ

3−2σ λn−2N λn= 1−3−2N σ

3−2σ λn−2N λn, we obtain

(3.15) U0n+1−U1n+1

1−2N λn−3−2N σ 3−2σ λn


The lower bound of the quantity in parenthesis as a function ofσis easily verified, and we have

U0n+1−U1n+1≥(1−3N λn)(U0n−U1n)≥0.

For1≤j ≤J0−1,

Ujn+1−Uj+1n+1=anjUjn+bnjUj−1n +cnjUj+1n −(anj+1Uj+1n +bnj+1Ujn+cnj+1Uj+2n ) + ∆tn

(Ujn)1+α−(Uj+1n )1+α

≥(anj +bnj)Ujn+cnjUj+1n −bnj+1Ujn−(anj+1+cnj+1)Uj+1n

= (anj +bnj −bnj+1)Ujn−(anj+1+cnj+1−cnj)Uj+1n



2j+ 1−2σ λn+ 2(N−1)σ 2(j+ 1) + 1−2σλn

(Ujn−Uj+1n ).

Since, forσ≥0,


2j+ 1−2σ λn+ 2(N−1)σ 2(j+ 1) + 1−2σλn

≥1−2λn− 2(N−1)(1−σ)

2j−1 + 2(1−σ)λn≥1−2λn−(N−1)λn≥1−2N λn≥0,


and forσ <0,


2j+ 1−2σ λn+ 2(N−1)σ 2(j+ 1) + 1−2σλn


2j+ 1−2σ λn+ 2(N−1)σ 2j+ 1−2σλn

= 1−2λn−2(N−1)(1−2σ) 2j+ 1−2σ λn

≥1−2λn−2(N−1)λn= 1−2N λn≥0,

we haveUjn+1≥Uj+1n+1for1≤j≤J0−1. Forj=J0, note that UJn+1

0 −UJn+1

0+1= (anJ0UJn0+bnJ0UJn0−1+cnJ0UJn0+1)



1− N−1 2(J0+ 1)



1 + N−1 2(J0+ 1)




+ ∆tn



1− N−1 2(J0+ 1)


(1−2λn) +λn

1 + N−1 2(J0+ 1)





2J0+ 1−2σ λn+ N−1 2(J0+ 1)λn



0 ≥UJn+1

0+1follows from 1−2λn−2(N−1)(1−σ)

2J0+ 1−2σ λn+ N−1 2(J0+ 1)λn

≥1−2λn− 2(N−1)(1−σ)

2J0−1 + 2(1−σ)λn≥1−2λn−(N−1)λn≥1−2N λn ≥0.

Forj ≥J0+ 1, the monotonicity can easily be proved by a similar argument given in [6].

Now the proof of monotonicity is completed by induction.

COROLLARY 3.5. Forσ= 2N3 , instead ofλ≤ 3N1 , it suffices to assumeλ≤ 2N1 to preserve the monotonicity of the numerical solution by virtue of (3.15). If0≤σ <3/(2N), then it is enough to assume thatλ≤1/(2N+ 1).

We now prove convergence of the approximate solution when the exact solution is smooth.

The restriction forλin the preceding lemma and corollary is not used here. It will be used in later sections when we consider the monotonicity of the solutions.

THEOREM3.6. Letλ≤ 2N1 andσ≤ 2N3 . Let{Ujn}be a solution of (3.5)–(3.8)and u(t, x)∈C2,4([0, T)×[0,1])be a solution of (1.2)–(1.5). Then for anyT0< T, there exists a positive constantC, depending only onT0and the initial data, such that

0≤j≤Jmax |Ujn−u(tn, rj)|< C(∆r)2, fortn ≤T0.


Proof.Under the stability conditionsλ≤2N1 andσ≤ 2N3 , the proof is essentially well known. We give a sketch of it for the sake of convenience. Setenj =Ujn−u(tn, rj). Note that the truncation errors for (3.5)–(3.7) are of order(∆r)2, whence we have by (3.12)–(3.14),

en+10 = (1−2N λn)en0+ 2N λnen1 + ∆tn

(U0n)1+α−(u(tn,0))1+α+O(∆r)2 , en+1j =

1−2λn−2(N−1)(1−2σ) 2j+ 1−2σ λn


1− 2σ(N−1) 2j+ 1−2σ


1 +2(1−σ)(N−1) 2j+ 1−2σ

enj+1 + ∆tn

(Ujn)1+α−(u(tn, rj))1+α+O(∆r)2 , for1≤j ≤J0, and forJ0+ 1≤j≤J−1,

en+1j = (1−2λn)enjn

1−N−1 2j


1 +N−1 2j

enj+1 + ∆tn

(Ujn)1+α−(u(tn, rj))1+α+O(∆r)2 . LetEn = max

0≤j≤J|enj|. Since all the coefficients are nonnegative under the stability assumptions λ≤2N1 andσ≤2N3 , it follows that


≤En+ ∆tn(1 +α) (u(tn, rj) +|En|)αEn+ ∆tnO (∆r)2 . Hence,

En+1≤(1 + ¯C∆tn)En+ ∆tnO (∆r)2 ,

as long asEn≤1, where the constantC¯depends only onuand its derivatives on[0, T0]×[0,1].

This inequality is then used to show thatEn≤1for allnif the initial error is small enough.

Then it is used once more to conclude thatEn=O (∆r)2

for alltn≤T0.

4. Behavior of the numerical solution. In this section, we assume for simplicity that the initial data are nonnegative and monotonically decreasing inr. We will give some sufficient conditions by which the numerical solution reproduces the blow-up characteristics.

Leta∈(0,1)be given, and defineWjnby

(4.1) Wjn2Ujn+ (1−a)(Ujn)1+α (∀j= 0, . . . , J−1), where

δ2U0n= 2NU1n−U0n (∆r)2 , δ2Ujn= Uj−1n −2Ujn+Ujn+1


+ 2j

2j+ 1−2σ N−1



∆r + (1−σ)Ujn+1−Ujn



(1≤j≤J0) and

δ2Ujn= Uj−1n −2Ujn+Uj+1n

(∆r)2 +N−1 rj

Uj+1n −Uj−1n

2∆r (J0+ 1≤j≤J−1).


In the following discussion, in addition to the nonnegativity and the monotonicity, we assume that the initial data satisfy

(H) Wj0≥0 (0≤ ∀j≤J−1) for somea∈(0,1).

REMARK 4.1. The assumption (H) is borrowed from [1]. It certainly imposes some restriction on the initial data, but it simplifies the analysis. Note that, by (3.5)–(3.7),Wjncan also be written as

(4.2) Wjn= Ujn+1−Ujn

∆tn −a(Ujn)1+α.

In view of (4.2) and the Dirichlet boundary condition (3.8), we thus setWJn= 0 (n≥0)for the sake of convenience.

REMARK4.2. Since (H) and (4.2) imply that Uj1−Uj0


≥a(Uj0)1+α≥0, it can be regarded as a discrete analogue of

(4.3) ut(0, r)≥0 (0< r <1),

which is a sufficient condition of the finite-time blow-up for the solution of (1.2)–(1.5).

See [14] for the details. Since we assume thata >0, assumption (H) is stronger than (4.3).

However, sincea >0can be chosen arbitrarily small, the restriction does not appear so serious.

Furthermore, it plays an important role when the blow-up rate of the discrete solution is taken into consideration. As a matter of fact, if assumption (H) is satisfied, then the blow-up of the discrete solution is always of Type I. See [10,18] and also the discussion in Section6. We however admit that assumption (H) may be a technical one.

LEMMA4.3.Let{Ujn}be the solution of (3.5)–(3.8). Letλ≤ 3N1 andσ≤2N3 . Assume that there existsa∈(0,1)such that(H)holds. ThenWjn ≥0, for allj= 0, . . . , J −1and n≥0. In particular, we have

(4.4) kUn+1k− kUnk

∆tn = U0n+1−U0n

∆tn ≥a(U0n)1+α (∀n≥0).

Proof.We follow the recipe used in [1], in which a one-dimensional discrete semilinear heat equation was considered. We write down the details for the readers’ convenience.

The proof is carried out by induction. Assume thatWjn ≥0for allj = 0, . . . , J −1.

By (4.2) and Lemma3.3, we have



≥a(Ujn)1+α≥0, which implies thatUjn+1≥Ujnfor allj. SetVjn :=U

n+1 j −Ujn

∆tn . We then obtain by (4.1), Wjn+1−Wjn


= δ2Ujn+1−δ2Ujn


+ (1−a)(Ujn+1)1+α−(Ujn)1+α


2Vjn+ (1−a)(Ujn+1)1+α−(Ujn)1+α

∆tn ,


and, by (4.2),δ2Wjn2Vjn−aδ2


. Therefore, we have Wjn+1−Wjn


−δ2Wjn= (1−a)(Ujn+1)1+α−(Ujn)1+α




for allj = 0, . . . , J −1. SinceUjn+1≥Ujnand0≤Uj+1n ≤Ujn (∀j= 0, . . . , J−1), we obtain the following four inequalities:


∆tn ≥(1 +α)(Ujn)αVjn (0≤ ∀j≤J−1), δ2


≥(1 +α)(U0n)αδ2U0n, δ2


= 1


1 +(1−σ)(N−1) 2j+ 1−2σ

(Uj+1n )1+α−(Ujn)1+α

1− σ(N−1) 2j+ 1−2σ

(Ujn)1+α−(Uj−1n )1+α

≥(1 +α)(Ujn)αδ2Ujn (1≤ ∀j≤J0), δ2


= 1


1 +N−1 2j

(Uj+1n )1+α−(Ujn)1+α

1−N−1 2j

(Ujn)1+α−(Uj−1n )1+α

≥(1 +α)(Ujn)αδ2Ujn (J0+ 1≤ ∀j ≤J−1).

These inequalities imply Wjn+1−Wjn


−δ2Wjn ≥(1−a)(1 +α)(Ujn)αVjn+a(1 +α)(Ujn)αδ2Ujn

= (1 +α)(Ujn)αWjn≥0

for allj = 0, . . . , J −1. Then the stability assumptionsλ≤ 3N1 andσ≤ 2N3 prove that Wjn+1 ≥ 0. Now (4.4) is a direct consequence of Lemma3.4and the fact thatW0n ≥ 0 (∀n≥0).

In the following two sections, we consider (3.5)–(3.8) with the temporal increments (3.9) or (3.10). In addition to the convergence of the numerical blow-up time, we show that each scheme can numerically reproduce a certain blow-up behavior well. For the adaptive one (3.9), we classify the exact numerical blow-up set. On the other hand, we use the uniform one to compute numerical evidence for (1.6) and (1.7).

4.1. Adaptive temporal increment(3.9). In this section, we consider (3.5)–(3.8) in which∆tnis given by (3.9). We define the numerical blow-up time byT(τ,∆r) =



∆tn. THEOREM4.4. LetT denote the blow-up time for the solution of (1.2)–(1.5), and let {Ujn}be the solution of (3.5)–(3.8)in which∆tnis given by(3.9). Letλ≤ 3N1 andσ≤2N3 . Assume that assumption(H)holds and that0< γ <1 +α+τ−1. Then we have

T(τ,∆r)<∞ and lim

∆r→0T(τ,∆r) =T.

By virtue of Lemma4.3, the proof can be carried out in almost the same way as in [11]. For details, we refer the readers to [11].


Since we assume that the initial datau0(r)are monotonically decreasing, it follows from Lemma3.4thatU0n =kUnk. Then the results in [5,14,25], namely blow-up occurs only at r= 0, suggest thatUjnremain bounded for all1≤j. We will examine if this is actually the case. The reader will see that the following argument goes parallel to those in [4,11]. However, because of the existence ofσ, it will be of some use if we explicitly state the theorems and their proofs.

For simplicity, we putγ=αandp=∞in (3.9) in the following discussion. For other choices ofγandp, the arguments given in [6] and [9,12] for the one-dimensional nonlinear heat equations can be applied if suitably modified.

LEMMA4.5.Ifλ≤ 3N1 andσ≤ 2N3 , then a blow-up solution{Ujn}satisfies



kUn+1k = 1 1 +cτ.

Proof. Note that we are concerned with only thosen’s for whichU0n = kUnk is sufficiently large, whence ∆tn = cτ(U0n)−α and λn = cλ(U0n)−α → 0 asn → ∞.

Therefore, we have by Lemma3.4, kUn+1k

kUnk = U0n+1

U0n = 1−2N λn+ 2N λn


U0n +cτ →1 +cτ as n→ ∞.

Note here that0< U1n/U0n≤1.

The contents of the next lemma is weaker than what we actually wish to obtain. But it subsequently plays an important role.

LEMMA4.6.Ifλ≤ 3N1 andσ≤ 2N3 , thenU1n/U0ntends to0asn→ ∞.

Proof. Note again that we are concerned with only thosen’s whenU0n =kUnk is sufficiently large, whence∆tn=cτ(U0n)−α. In order to save notational spaces, we set

σ1= 1−2σ(N−1)

3−2σ and σ2= 1 +2(1−σ)(N−1) 3−2σ . Then, by (3.12), (3.13), we have


U0n+1 ={1−(σ12n}U1nnσ1U0nnσ2U2n+ ∆tn(U1n)1+α (1−2N λn)U0n+ 2N λnU1n+ ∆tn(U0n)1+α . Setan=U1n/U0n. Then we have by monotonicity


{1−(σ12n}annσ1nσ2 U2n

U0n +a1+αn ∆tn(U0n)α 1−2N λn+ 2N λnan+ ∆tn(U0n)α

≤an−λnσ1annσ1+a1+αn ∆tn(U0n)α 1−2N λn+ 2N λnan+ ∆tn(U0n)α . (4.5)

Assume that0≤a:= lim inf

n→∞ an <lim sup


an =:a≤1. Then, for any givenκ∈(a, a), one may choose a subsequence{ank}such that

ank < κ and ank+1≥κ.

By (4.5), we have

κ≤lim sup


ank+1≤ 1 +cτ κα 1 +cτ κ < κ,


which is a contradiction. This proves the existence ofa= lim

n→∞an ∈[0,1]. Again by (4.5) we have

a≤ 1 +cτ aα 1 +cτ a,

which, together with the fact thata∈[0,1], impliesa= 0or1. Note that we have by (4.5),

(4.6) an+1


≤ annσ1(1−an) + (an)1+α∆tn(U0n)α an−2N λnan(1−an) +an∆tn(U0n)α =: Cn


. We now use the following inequality, which is proved in an elementary way:

(1−x)≤ 1

α(1−xα) (0< α <1, 0< x <1).

Then for0< α <1we have Cn−Dn = (1−an)∆tn

(∆r)−2σ1+ 2N an(∆r)−2−an−a1+αn 1−an (U0n)α



On the other hand, if1≤α, we have1−an ≤1−(an)α, whence

(4.8) Cn−Dn≤(1−an)∆tn[(∆r)−2(2N+σ1)−an(U0n)α].

Ifliman =a= 1, thenan(U0n)α→ ∞asn→ ∞. Then it follows in either case of (4.7) or (4.8) that{an}is decreasing for largen, which contradictsa= 1. We therefore havea= 0.

COROLLARY4.7.It holds that

n→∞lim an+1



( (1 +cτ)−1 σ= 2N3 , (1 +cτ)min{α,1} σ < 2N3 .

Proof. We prove the statement in the case thatσ= 2N3 . Note that we have by (3.12), (3.13), and Lemmas3.3,3.4,

an+1≥ (1−σ2λn)an+cτ a1+αn 1−2N λn+ 2N λnan+cτ, which implies


an ≥ 1−σ2λn+cτ aαn 1−2N λn+ 2N λnan+cτ.

Therefore we havelim infan+1/an ≥ 1/(1 +cτ). On the other hand, the upper bound lim supan+1/an≤1/(1 +cτ)is easily derived from (4.6).

The proof in the case ofσ < 2N3 is rather complicated. However, it can be carried out in a way similar to that in Corollary 4.2 of [11] or Theorem 1 of [4], and we may omit it.

We now state and prove the asymptotic behavior in two cases separately:

THEOREM 4.8. Let{Ujn}be a blow-up solution of (3.5)–(3.8), and letλ ≤ 2N1 and σ= 2N3 . Then, for allα >0, we haveU0n → ∞asn→ ∞, whileUjnremain bounded for all1≤j≤J.




Related subjects :