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

FINITE DIFFERENCE SCHEMES FOR AN AXISYMMETRIC NONLINEAR HEAT EQUATION WITH BLOW-UP^∗

CHIEN-HONG CHO^†ANDHISASHI OKAMOTO^‡

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 theL^p-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 theL^p-norm(1≤p <∞)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 theL^p-norm blow-up. The computational results are also analyzed. Moreover, we prove an abstract theorem which shows the relationship between the numericalL^p-norm blow-up and the exactL^p-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,L^p-norm blow-up AMS subject classifications.65M06, 65M12

1. Introduction. LetΩbe the unit ball inR^N 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+u^1+α,

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:

u_t=u_rr+N−1

r u_r+u^1+α (0≤r <1,0< t), (1.2)

ur= 0 (r= 0),

(1.3)

u= 0 (r= 1),

(1.4)

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

(1.5)

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. (chcho20@ccu.edu.tw).

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

391

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 theL^p-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,·)kL^p=∞ if p > N α 2 , (1.6)

lim sup

t→T

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

Weissler [30] showed that (1.6) can be extended to the critical casep = ^{N α}₂ 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 nonlinearityu^1+α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 L^p-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

(2.1)







u_t=u_xx+u² (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:











U_jⁿ⁺¹−U_jⁿ

∆tn

=U_j+1ⁿ −2U_jⁿ+U_j−1ⁿ

(∆x)² + (U_jⁿ)² (j= 1, . . . , J−1, n≥0), U_j⁰=u₀(x_j) (j= 1, . . . , J−1),

U₀ⁿ =U_Jⁿ = 0 (n≥0),

whose temporal grid sizes are defined adaptively by

(2.2) ∆t_n=τ·min

1, c

kUⁿk2

.

Here,J ∈N,∆x= ¹_J is the spatial grid size,xj =j∆x(0 ≤j ≤J)are the spatial grid points,τ andcare prescribed parameters, andλ= _(∆x)^τ ₂ ≤ ¹₂ is fixed. Moreover,t₀ = 0 andtn = t_n−1+ ∆t_n−1 (n ≥ 1)denote the temporal grid points, and U_jⁿ denotes the approximation of the exact solutionuof (2.1) at(tn, xj). The discreteL^p-norm is defined by

kUⁿk_p =









 max

j=1,...,J−1|U_jⁿ| ifp=∞,

J−1

P

j=1

∆x|U_jⁿ|^p

!¹_p

if1≤p <∞.

The temporal grid size∆tn in (2.2) is defined at each time step according to the discrete L²-norm of the numerical solution. As a matter of fact, the discreteL²-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 whileUⁿis not so large, and an adaptive mesh begins to be in effect onceUⁿbecomes large. The author in [26] defined the numerical blow-up time asT(τ,∆x) =P∞

n=0∆tn. Then he proved that the numerical solutionU_jⁿconverges 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 theL^p-norms (1≤p <∞) blow up simultaneously with theL^∞-norm. In fact, the unboundedness of the discreteL^∞-norm is equivalent to the unboundedness of other discreteL^p-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 theL^p-norm.

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

u⁰(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, G⁰, G⁰⁰>0andR∞

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

We now consider the forward Euler scheme

(2.3) Uⁿ⁺¹−Uⁿ

∆t =G(Uⁿ), U⁰=u(0).

Note that the numerical solutions can be computed for any timet_n=n∆teven though the exact solutionushows a blown up in finite timeT_O. 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,Uⁿis monotonically increasing innand

Uⁿ→ ∞ as n→ ∞.

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

∆t·H U^n∆t−1

<1 and ∆t·H(U^n∆t)≥1.

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

THEOREM2.1.It holds that

− Z ∞

H⁻¹(∆t¹ ) ds

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

H⁻¹(∆t¹ ) ds

G(s)+ ∆t 1 + logG H⁻¹ _∆t¹ G(U⁰)

! , whereH⁻¹denotes the inverse function ofH. In particular, we have lim

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

(2.4) ∆t·log

G

H⁻¹

1

∆t

→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,

U^n∆t ≥H⁻¹ 1

∆t

→ ∞ as ∆t→0,

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

T_O(∆t)→T_O 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) =u^l(l >1), then (2.4) can be replaced by

(2.5) ∆t·log

H⁻¹

1

∆t

→0 as ∆t→0.

In this case, it is easy to see thatH(s) =s^q(∀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) U_jⁿ⁺¹−U_jⁿ

∆t_n = U_j+1ⁿ⁺¹−2U_jⁿ⁺¹+U_j−1ⁿ⁺¹

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

U_j+1ⁿ⁺¹−U_j−1ⁿ⁺¹

2∆r + U_j^n1+α , whereU_jⁿdenotes an approximation foru(t_n, j∆r). Here,∆ris the spatial grid size, and the temporal increment∆t_nis defined adaptively by

∆t_n=τ·min

1, c kUⁿk^α_p

(1≤p≤ ∞),

whereτandcare prescribed parameters andkUⁿkpdenotes the discreteL^p-norm:

(3.2) kUⁿkp =











J−1

P

j=0

∆r((j+ 1)∆r)^N−1|U_jⁿ|^p

!¹_p

(1≤p <∞), max

j=0,...,J|U_jⁿ| (p=∞).

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

(3.3) lim

r→0

ur(t, r) r = lim

r→0

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

r =u_rr(t,0).

Accordingly, one may assume the following equation:

(3.4) U₀ⁿ⁺¹−U₀ⁿ

∆t_n =N−2U₀ⁿ⁺¹+ 2U₁ⁿ⁺¹

(∆r)² + (U₀ⁿ)^1+α.

This is nothing but an approximation of axisymmetric solutions of u_t = 4u+u^1+α 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):

U_jⁿ⁺¹−U_jⁿ

∆tn

=NU_j+1ⁿ⁺¹−2U_jⁿ⁺¹+U_j−1ⁿ⁺¹

(∆r)² + (U_jⁿ)^1+α for1≤j < N0:=

N+ 1 2

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

2

denotes the largest integer≤ ^N₂⁺¹. He proved that his scheme is stable for allN ≥2and converges in the orderO (∆r)²

. 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

1,_N−1

2

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

(3.5) U₀ⁿ⁺¹−U₀ⁿ

∆t_n =N−2U₀ⁿ+ 2U₁ⁿ

(∆r)² + (U₀ⁿ)^1+α and

U_jⁿ⁺¹−U_jⁿ

∆tn

= U_j+1ⁿ −2U_jⁿ+U_j−1ⁿ

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

U_j+1ⁿ −U_j−1ⁿ

2∆r + U_j^n1+α , J0+ 1≤j≤J−1.

(3.6)

For1≤j ≤J₀, we introduce a new parameterσand approximateu_rin the following way:

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)

∆r

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

2r u_rr(t, r) +O((∆r)²).

Note that, for smallr,

ur(t, r)

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

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

N−1 r u_r(t, r)

= N−1 r

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

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

∆r

+O(r²) +O (∆r)² .

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,

U_jⁿ⁺¹−U_jⁿ

∆tn

=U_j−1ⁿ −2U_jⁿ+U_j+1ⁿ (∆r)²

+ 2j

2j+ 1−2σ N−1

rj

σU_jⁿ−U_j−1ⁿ

∆r + (1−σ)U_j+1ⁿ −U_jⁿ

∆r

+ U_j^n1+α . (3.7)

Note that the truncation error for the scheme (3.7) is of order(∆r)²for 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 setU_j⁰=u0(rj)≥0, for0≤j≤J, and

(3.8) U_Jⁿ = 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

kUⁿk^γp

(γ >0),

or be given uniformly by

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

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

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

(3.11) U₀ⁿ=U₁ⁿ or U₀ⁿ= (4U₁ⁿ−U₂ⁿ)/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,U₀ⁿ=U₁ⁿwould 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) sinceU₁ⁿ/U₀ⁿ≥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)^τ ₂ andλ_n=_(∆r)^∆tn₂ (∀n≥0). Note that now (3.5)–(3.7) can be rewritten as (3.12) U₀ⁿ⁺¹= (1−2N λn)U₀ⁿ+ 2N λnU₁ⁿ+ ∆tn(U₀ⁿ)^1+α

for1≤j ≤J0, U_jⁿ⁺¹=

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

U_jⁿ +λn

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

U_j−1ⁿ +λn

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

U_j+1ⁿ + ∆tn(U_jⁿ)^1+α

(3.13)

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

U_jⁿ⁺¹= (1−2λn)U_jⁿ+λn

1−N−1 2j

U_j−1ⁿ +λn

1 +N−1 2j

U_j+1ⁿ + ∆tn(U_jⁿ)^1+α.

(3.14)

LEMMA3.3.Suppose thatu0(r)is nonnegative. Letλ≤ _2N¹ be fixed andσ≤ _2N³ . Let {U_jⁿ}be a solution of (3.5)–(3.8). Then,U_jⁿ ≥0, for alln≥0and0≤j≤J.

Proof. The proof is carried out by induction. Assume thatU_jⁿ ≥0for all0 ≤j ≤J. Forj = 0andJ0+ 1 ≤ j ≤ J −1,U_jⁿ⁺¹ ≥ 0follows directly from (3.12), (3.14), and the assumptionλ≤ _2N¹ . In the case1≤j ≤J₀, 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

aⁿ_j = 1−2λ_n−2(N−1)(1−2σ) 2j+ 1−2σ λ_n, bⁿ_j =λn

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

, cⁿ_j =λn

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

. Sinceσ≤_2N³ and1≤j, it is easy to see thatcⁿ_j ≥0for alln≥0andj. Observe also that bⁿ_j >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 bⁿ_j. Now it remains to show the nonnegativity of aⁿ_j. For

1

2≤σ≤ _2N³ , we have

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

whence we obtainaⁿ_j ≥1−2λn≥0. In the case ofσ < ¹₂, note that 2(N−1)(1−2σ)

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

aⁿ_j >1−2λ_n−2(N−1)λ_n = 1−2N λ_n≥0.

Thus, we have proved the nonnegativity ofaⁿ_j,bⁿ_j,cⁿ_j, and we are done.

LEMMA 3.4. Suppose thatu₀(r)is nonnegative and monotonically decreasing. Let λ≤ _3N¹ andσ ≤ _2N³ . Let{U_jⁿ}be a solution of (3.5)–(3.8). Then,U_jⁿ is monotonically decreasing inj. Namely,U_jⁿ ≥U_j+1ⁿ ≥0, for alln≥0and0≤j≤J−1.

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

U₀ⁿ⁺¹−U₁ⁿ⁺¹= (1−2N λn)U₀ⁿ+ 2N λnU₁ⁿ−(aⁿ₁U₁ⁿ+bⁿ₁U₀ⁿ+cⁿ₁U₂ⁿ) + ∆t_n

(U₀ⁿ)^1+α−(U₁ⁿ)^1+α

≥(1−2N λn−bⁿ₁)U₀ⁿ−(aⁿ₁+cⁿ₁−2N λn)U₁ⁿ. Since

1−2N λn−bⁿ₁ = 1−2N λn−3−2N σ 3−2σ λn

and

aⁿ₁+cⁿ₁−2N λn= 1−λn+2(N−1)σ

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

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

(3.15) U₀ⁿ⁺¹−U₁ⁿ⁺¹≥

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

(U₀ⁿ−U₁ⁿ).

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

U₀ⁿ⁺¹−U₁ⁿ⁺¹≥(1−3N λn)(U₀ⁿ−U₁ⁿ)≥0.

For1≤j ≤J0−1,

U_jⁿ⁺¹−U_j+1ⁿ⁺¹=aⁿ_jU_jⁿ+bⁿ_jU_j−1ⁿ +cⁿ_jU_j+1ⁿ −(aⁿ_j+1U_j+1ⁿ +bⁿ_j+1U_jⁿ+cⁿ_j+1U_j+2ⁿ ) + ∆tn

(U_jⁿ)^1+α−(U_j+1ⁿ )^1+α

≥(aⁿ_j +bⁿ_j)U_jⁿ+cⁿ_jU_j+1ⁿ −bⁿ_j+1U_jⁿ−(aⁿ_j+1+cⁿ_j+1)U_j+1ⁿ

= (aⁿ_j +bⁿ_j −bⁿ_j+1)U_jⁿ−(aⁿ_j+1+cⁿ_j+1−cⁿ_j)U_j+1ⁿ

=

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

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

(U_jⁿ−U_j+1ⁿ ).

Since, forσ≥0,

1−2λ_n−2(N−1)(1−σ)

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,

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

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

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

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 haveU_jⁿ⁺¹≥U_j+1ⁿ⁺¹for1≤j≤J0−1. Forj=J0, note that U_Jⁿ⁺¹

0 −U_Jⁿ⁺¹

0+1= (aⁿ_J0U_Jⁿ₀+bⁿ_J0U_Jⁿ₀₋₁+cⁿ_J0U_Jⁿ₀₊₁)

−

"

(1−2λn)U_Jⁿ₀₊₁+λn

1− N−1 2(J0+ 1)

U_Jⁿ₀

+λ_n

1 + N−1 2(J0+ 1)

U_Jⁿ

0+2

#

+ ∆tn

(U_Jⁿ₀)^1+α−(U_Jⁿ₀₊₁)^1+α

≥

aⁿ_J0+bⁿ_J0−λn

1− N−1 2(J₀+ 1)

U_Jⁿ₀

−

(1−2λn) +λn

1 + N−1 2(J0+ 1)

−cⁿ_J0

U_Jⁿ₀₊₁

=

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

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

(U_Jⁿ₀−U_Jⁿ₀₊₁).

NowU_Jⁿ⁺¹

0 ≥U_Jⁿ⁺¹

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σ= _2N³ , instead ofλ≤ _3N¹ , it suffices to assumeλ≤ _2N¹ 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λ≤ _2N¹ andσ≤ _2N³ . Let{U_jⁿ}be a solution of (3.5)–(3.8)and u(t, x)∈C^2,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 |U_jⁿ−u(tn, rj)|< C(∆r)², fortn ≤T0.

Proof.Under the stability conditionsλ≤_2N¹ andσ≤ _2N³ , the proof is essentially well known. We give a sketch of it for the sake of convenience. Seteⁿ_j =U_jⁿ−u(tn, rj). Note that the truncation errors for (3.5)–(3.7) are of order(∆r)², whence we have by (3.12)–(3.14),

eⁿ⁺¹₀ = (1−2N λn)eⁿ₀+ 2N λneⁿ₁ + ∆tn

(U₀ⁿ)^1+α−(u(tn,0))^1+α+O(∆r)² , eⁿ⁺¹_j =

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

eⁿ_j +λn

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

eⁿ_j−1 +λn

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

eⁿ_j+1 + ∆tn

(U_jⁿ)^1+α−(u(tn, rj))^1+α+O(∆r)² , for1≤j ≤J0, and forJ0+ 1≤j≤J−1,

eⁿ⁺¹_j = (1−2λ_n)eⁿ_j +λ_n

1−N−1 2j

eⁿ_j−1+λ_n

1 +N−1 2j

eⁿ_j+1 + ∆t_n

(U_jⁿ)^1+α−(u(t_n, r_j))^1+α+O(∆r)² . LetEⁿ = max

0≤j≤J|eⁿ_j|. Since all the coefficients are nonnegative under the stability assumptions λ≤_2N¹ andσ≤_2N³ , it follows that

eⁿ⁺¹_j

≤Eⁿ+ ∆tn(1 +α) (u(tn, rj) +|Eⁿ|)^αEⁿ+ ∆tnO (∆r)² . Hence,

Eⁿ⁺¹≤(1 + ¯C∆t_n)Eⁿ+ ∆t_nO (∆r)² ,

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

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

Then it is used once more to conclude thatEⁿ=O (∆r)²

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 defineW_jⁿby

(4.1) W_jⁿ =δ²U_jⁿ+ (1−a)(U_jⁿ)^1+α (∀j= 0, . . . , J−1), where

δ²U₀ⁿ= 2NU₁ⁿ−U₀ⁿ (∆r)² , δ²U_jⁿ= U_j−1ⁿ −2U_jⁿ+U_jⁿ⁺¹

(∆r)²

+ 2j

2j+ 1−2σ N−1

rj

σU_jⁿ−U_j−1ⁿ

∆r + (1−σ)U_jⁿ⁺¹−U_jⁿ

∆r

!

(1≤j≤J0) and

δ²U_jⁿ= U_j−1ⁿ −2U_jⁿ+U_j+1ⁿ

(∆r)² +N−1 r_j

U_j+1ⁿ −U_j−1ⁿ

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) W_j⁰≥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),W_jⁿcan also be written as

(4.2) W_jⁿ= U_jⁿ⁺¹−U_jⁿ

∆t_n −a(U_jⁿ)^1+α.

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

REMARK4.2. Since (H) and (4.2) imply that U_j¹−U_j⁰

∆t0

≥a(U_j⁰)^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{U_jⁿ}be the solution of (3.5)–(3.8). Letλ≤ _3N¹ andσ≤_2N³ . Assume that there existsa∈(0,1)such that(H)holds. ThenW_jⁿ ≥0, for allj= 0, . . . , J −1and n≥0. In particular, we have

(4.4) kUⁿ⁺¹k_∞− kUⁿk_∞

∆t_n = U₀ⁿ⁺¹−U₀ⁿ

∆t_n ≥a(U₀ⁿ)^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 thatW_jⁿ ≥0for allj = 0, . . . , J −1.

By (4.2) and Lemma3.3, we have

U_jⁿ⁺¹−U_jⁿ

∆tn

≥a(U_jⁿ)^1+α≥0, which implies thatU_jⁿ⁺¹≥U_jⁿfor allj. SetV_jⁿ :=^U

n+1 j −U_jⁿ

∆t_n . We then obtain by (4.1), W_jⁿ⁺¹−W_jⁿ

∆tn

= δ²U_jⁿ⁺¹−δ²U_jⁿ

∆tn

+ (1−a)(U_jⁿ⁺¹)^1+α−(U_jⁿ)^1+α

∆tn

=δ²V_jⁿ+ (1−a)(U_jⁿ⁺¹)^1+α−(U_jⁿ)^1+α

∆t_n ,

and, by (4.2),δ²W_jⁿ=δ²V_jⁿ−aδ²

(U_jⁿ)^1+α

. Therefore, we have W_jⁿ⁺¹−W_jⁿ

∆tn

−δ²W_jⁿ= (1−a)(U_jⁿ⁺¹)^1+α−(U_jⁿ)^1+α

∆tn

+aδ²

(U_jⁿ)^1+α

for allj = 0, . . . , J −1. SinceU_jⁿ⁺¹≥U_jⁿand0≤U_j+1ⁿ ≤U_jⁿ (∀j= 0, . . . , J−1), we obtain the following four inequalities:

(U_jⁿ⁺¹)^1+α−(U_jⁿ)^1+α

∆t_n ≥(1 +α)(U_jⁿ)^αV_jⁿ (0≤ ∀j≤J−1), δ²

(U₀ⁿ)^1+α

≥(1 +α)(U₀ⁿ)^αδ²U₀ⁿ, δ²

(U_jⁿ)^1+α

= 1

(∆r)²

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

(U_j+1ⁿ )^1+α−(U_jⁿ)^1+α

−

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

(U_jⁿ)^1+α−(U_j−1ⁿ )^1+α

≥(1 +α)(U_jⁿ)^αδ²U_jⁿ (1≤ ∀j≤J₀), δ²

(U_jⁿ)^1+α

= 1

(∆r)²

1 +N−1 2j

(U_j+1ⁿ )^1+α−(U_jⁿ)^1+α

−

1−N−1 2j

(U_jⁿ)^1+α−(U_j−1ⁿ )^1+α

≥(1 +α)(U_jⁿ)^αδ²U_jⁿ (J0+ 1≤ ∀j ≤J−1).

These inequalities imply W_jⁿ⁺¹−W_jⁿ

∆tn

−δ²W_jⁿ ≥(1−a)(1 +α)(U_jⁿ)^αV_jⁿ+a(1 +α)(U_jⁿ)^αδ²U_jⁿ

= (1 +α)(U_jⁿ)^αW_jⁿ≥0

for allj = 0, . . . , J −1. Then the stability assumptionsλ≤ _3N¹ andσ≤ _2N³ prove that W_jⁿ⁺¹ ≥ 0. Now (4.4) is a direct consequence of Lemma3.4and the fact thatW₀ⁿ ≥ 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) =

∞

P

n=0

∆tn. THEOREM4.4. LetT denote the blow-up time for the solution of (1.2)–(1.5), and let {U_jⁿ}be the solution of (3.5)–(3.8)in which∆tnis given by(3.9). Letλ≤ _3N¹ andσ≤_2N³ . Assume that assumption(H)holds and that0< γ <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.4thatU₀ⁿ =kUⁿk_∞. Then the results in [5,14,25], namely blow-up occurs only at r= 0, suggest thatU_jⁿremain 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λ≤ _3N¹ andσ≤ _2N³ , then a blow-up solution{U_jⁿ}satisfies

n→∞lim

kUⁿk∞

kUⁿ⁺¹k_∞ = 1 1 +cτ.

Proof. Note that we are concerned with only thosen’s for whichU₀ⁿ = kUⁿk∞ is sufficiently large, whence ∆tn = cτ(U₀ⁿ)^−α and λn = cλ(U₀ⁿ)^−α → 0 asn → ∞.

Therefore, we have by Lemma3.4, kUⁿ⁺¹k_∞

kUⁿk_∞ = U₀ⁿ⁺¹

U₀ⁿ = 1−2N λn+ 2N λn

U₁ⁿ

U₀ⁿ +cτ →1 +cτ as n→ ∞.

Note here that0< U₁ⁿ/U₀ⁿ≤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λ≤ _3N¹ andσ≤ _2N³ , thenU₁ⁿ/U₀ⁿtends to0asn→ ∞.

Proof. Note again that we are concerned with only thosen’s whenU₀ⁿ =kUⁿk_∞ is sufficiently large, whence∆tn=cτ(U₀ⁿ)^−α. 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

U₁ⁿ⁺¹

U₀ⁿ⁺¹ ={1−(σ1+σ2)λn}U₁ⁿ+λnσ1U₀ⁿ+λnσ2U₂ⁿ+ ∆tn(U₁ⁿ)^1+α (1−2N λ_n)U₀ⁿ+ 2N λ_nU₁ⁿ+ ∆t_n(U₀ⁿ)^1+α . Setan=U₁ⁿ/U₀ⁿ. Then we have by monotonicity

an+1=

{1−(σ1+σ2)λn}an+λnσ1+λnσ2 U₂ⁿ

U₀ⁿ +a^1+α_n ∆tn(U₀ⁿ)^α 1−2N λ_n+ 2N λ_na_n+ ∆t_n(U₀ⁿ)^α

≤a_n−λ_nσ₁a_n+λ_nσ₁+a^1+α_n ∆t_n(U₀ⁿ)^α 1−2N λn+ 2N λnan+ ∆tn(U₀ⁿ)^α . (4.5)

Assume that0≤a:= lim inf

n→∞ an <lim sup

n→∞

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

a_nk < κ and a_nk+1≥κ.

By (4.5), we have

κ≤lim sup

k→∞

an_k+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) a_n+1

an

≤ a_n+λ_nσ₁(1−a_n) + (a_n)^1+α∆t_n(U₀ⁿ)^α an−2N λnan(1−an) +an∆tn(U₀ⁿ)^α =: C_n

Dn

. 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)⁻²σ1+ 2N an(∆r)⁻²−an−a^1+α_n 1−a_n (U₀ⁿ)^α

≤(1−an)∆tn[(∆r)⁻²(2N+σ1)−αan(U₀ⁿ)^α].

(4.7)

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

(4.8) C_n−D_n≤(1−a_n)∆t_n[(∆r)⁻²(2N+σ₁)−a_n(U₀ⁿ)^α].

Ifliman =a= 1, thenan(U₀ⁿ)^α→ ∞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 a_n+1

an

=

( (1 +cτ)⁻¹ σ= _2N³ , (1 +cτ)^−min{α,1} σ < _2N³ .

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

a_n+1≥ (1−σ₂λ_n)a_n+cτ a^1+α_n 1−2N λn+ 2N λnan+cτ, which implies

an+1

a_n ≥ 1−σ2λn+cτ a^α_n 1−2N λ_n+ 2N λ_na_n+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σ < _2N³ 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{U_jⁿ}be a blow-up solution of (3.5)–(3.8), and letλ ≤ _2N¹ and σ= _2N³ . Then, for allα >0, we haveU₀ⁿ → ∞asn→ ∞, whileU_jⁿremain bounded for all1≤j≤J.