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 α}_{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 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^{2} (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}^{n+1}−U_{j}^{n}

∆tn

=U_{j+1}^{n} −2U_{j}^{n}+U_{j−1}^{n}

(∆x)^{2} + (U_{j}^{n})^{2} (j= 1, . . . , J−1, n≥0),
U_{j}^{0}=u_{0}(x_{j}) (j= 1, . . . , J−1),

U_{0}^{n} =U_{J}^{n} = 0 (n≥0),

whose temporal grid sizes are defined adaptively by

(2.2) ∆t_{n}=τ·min

1, c

kU^{n}k2

.

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

kU^{n}k_{p} =

max

j=1,...,J−1|U_{j}^{n}| ifp=∞,

J−1

P

j=1

∆x|U_{j}^{n}|^{p}

!^{1}_{p}

if1≤p <∞.

The temporal grid size∆tn in (2.2) is defined at each time step according to the discrete
L^{2}-norm of the numerical solution. As a matter of fact, the discreteL^{2}-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^{n}is not so large, and an adaptive mesh begins
to be in effect onceU^{n}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}^{n}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^{0}(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^{0}, G^{00}>0andR∞

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

We now consider the forward Euler scheme

(2.3) U^{n+1}−U^{n}

∆t =G(U^{n}), U^{0}=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^{n}is monotonically increasing innand

U^{n}→ ∞ 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^{−1}(∆t^{1} )
ds

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

H^{−1}(∆t^{1} )
ds

G(s)+ ∆t 1 + logG H^{−1} _{∆t}^{1}
G(U^{0})

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

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

(2.4) ∆t·log

G

H^{−1}

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}
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}

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}^{n+1}−U_{j}^{n}

∆t_{n} = U_{j+1}^{n+1}−2U_{j}^{n+1}+U_{j−1}^{n+1}

(∆r)^{2} +N−1
j∆r

U_{j+1}^{n+1}−U_{j−1}^{n+1}

2∆r + U_{j}^{n}^{1+α}
,
whereU_{j}^{n}denotes an approximation foru(t_{n}, j∆r). Here,∆ris the spatial grid size, and the
temporal increment∆t_{n}is defined adaptively by

∆t_{n}=τ·min

1, c
kU^{n}k^{α}_{p}

(1≤p≤ ∞),

whereτandcare prescribed parameters andkU^{n}kpdenotes the discreteL^{p}-norm:

(3.2) kU^{n}kp =

J−1

P

j=0

∆r((j+ 1)∆r)^{N}^{−1}|U_{j}^{n}|^{p}

!^{1}_{p}

(1≤p <∞), max

j=0,...,J|U_{j}^{n}| (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_{0}^{n+1}−U_{0}^{n}

∆t_{n} =N−2U_{0}^{n+1}+ 2U_{1}^{n+1}

(∆r)^{2} + (U_{0}^{n})^{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}^{n+1}−U_{j}^{n}

∆tn

=NU_{j+1}^{n+1}−2U_{j}^{n+1}+U_{j−1}^{n+1}

(∆r)^{2} + (U_{j}^{n})^{1+α} for1≤j < N0:=

N+ 1 2

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

2

denotes the largest integer≤ ^{N}_{2}^{+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

1,_{N−1}

2

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

(3.5) U_{0}^{n+1}−U_{0}^{n}

∆t_{n} =N−2U_{0}^{n}+ 2U_{1}^{n}

(∆r)^{2} + (U_{0}^{n})^{1+α}
and

U_{j}^{n+1}−U_{j}^{n}

∆tn

= U_{j+1}^{n} −2U_{j}^{n}+U_{j−1}^{n}

(∆r)^{2} +N−1
j∆r

U_{j+1}^{n} −U_{j−1}^{n}

2∆r + U_{j}^{n}^{1+α}
,
J0+ 1≤j≤J−1.

(3.6)

For1≤j ≤J_{0}, we introduce a new parameterσand approximateu_{r}in 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)^{2}).

Note that, for smallr,

ur(t, r)

r −urr(t, r) =O r^{2}
.
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^{2}) +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,

U_{j}^{n+1}−U_{j}^{n}

∆tn

=U_{j−1}^{n} −2U_{j}^{n}+U_{j+1}^{n}
(∆r)^{2}

+ 2j

2j+ 1−2σ N−1

rj

σU_{j}^{n}−U_{j−1}^{n}

∆r + (1−σ)U_{j+1}^{n} −U_{j}^{n}

∆r

+ U_{j}^{n}^{1+α}
.
(3.7)

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

(3.8) U_{J}^{n} = 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^{n}k^{γ}p

(γ >0),

or be given uniformly by

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

In these equations,τ >0is a prescribed small parameter, andkU^{n}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_{0}^{n}=U_{1}^{n} or U_{0}^{n}= (4U_{1}^{n}−U_{2}^{n})/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_{0}^{n}=U_{1}^{n}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_{1}^{n}/U_{0}^{n}≥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)}^{∆t}^{n}_{2} (∀n≥0). Note that now (3.5)–(3.7) can be rewritten as
(3.12) U_{0}^{n+1}= (1−2N λn)U_{0}^{n}+ 2N λnU_{1}^{n}+ ∆tn(U_{0}^{n})^{1+α}

for1≤j ≤J0,
U_{j}^{n+1}=

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

U_{j}^{n}
+λn

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

U_{j−1}^{n} +λn

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

U_{j+1}^{n}
+ ∆tn(U_{j}^{n})^{1+α}

(3.13)

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

U_{j}^{n+1}= (1−2λn)U_{j}^{n}+λn

1−N−1 2j

U_{j−1}^{n} +λn

1 +N−1 2j

U_{j+1}^{n}
+ ∆tn(U_{j}^{n})^{1+α}.

(3.14)

LEMMA3.3.Suppose thatu0(r)is nonnegative. Letλ≤ _{2N}^{1} be fixed andσ≤ _{2N}^{3} . Let
{U_{j}^{n}}be a solution of (3.5)–(3.8). Then,U_{j}^{n} ≥0, for alln≥0and0≤j≤J.

Proof. The proof is carried out by induction. Assume thatU_{j}^{n} ≥0for all0 ≤j ≤J.
Forj = 0andJ0+ 1 ≤ j ≤ J −1,U_{j}^{n+1} ≥ 0follows directly from (3.12), (3.14), and
the assumptionλ≤ _{2N}^{1} . In the case1≤j ≤J_{0}, 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^{n}_{j} = 1−2λ_{n}−2(N−1)(1−2σ)
2j+ 1−2σ λ_{n},
b^{n}_{j} =λn

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

, c^{n}_{j} =λn

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

.
Sinceσ≤_{2N}^{3} and1≤j, it is easy to see thatc^{n}_{j} ≥0for alln≥0andj. Observe also that
b^{n}_{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^{n}_{j}. Now it remains to show the nonnegativity of a^{n}_{j}. For

1

2≤σ≤ _{2N}^{3} , we have

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

whence we obtaina^{n}_{j} ≥1−2λn≥0. In the case ofσ < ^{1}_{2}, note that
2(N−1)(1−2σ)

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

a^{n}_{j} >1−2λ_{n}−2(N−1)λ_{n} = 1−2N λ_{n}≥0.

Thus, we have proved the nonnegativity ofa^{n}_{j},b^{n}_{j},c^{n}_{j}, and we are done.

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

Proof.We prove this by induction. Assume thatU_{j}^{n} ≥U_{j+1}^{n} ≥0for a certainnand all
0≤j≤J−1. Note that forj= 0,

U_{0}^{n+1}−U_{1}^{n+1}= (1−2N λn)U_{0}^{n}+ 2N λnU_{1}^{n}−(a^{n}_{1}U_{1}^{n}+b^{n}_{1}U_{0}^{n}+c^{n}_{1}U_{2}^{n})
+ ∆t_{n}

(U_{0}^{n})^{1+α}−(U_{1}^{n})^{1+α}

≥(1−2N λn−b^{n}_{1})U_{0}^{n}−(a^{n}_{1}+c^{n}_{1}−2N λn)U_{1}^{n}.
Since

1−2N λn−b^{n}_{1} = 1−2N λn−3−2N σ
3−2σ λn

and

a^{n}_{1}+c^{n}_{1}−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_{0}^{n+1}−U_{1}^{n+1}≥

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

(U_{0}^{n}−U_{1}^{n}).

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

U_{0}^{n+1}−U_{1}^{n+1}≥(1−3N λn)(U_{0}^{n}−U_{1}^{n})≥0.

For1≤j ≤J0−1,

U_{j}^{n+1}−U_{j+1}^{n+1}=a^{n}_{j}U_{j}^{n}+b^{n}_{j}U_{j−1}^{n} +c^{n}_{j}U_{j+1}^{n} −(a^{n}_{j+1}U_{j+1}^{n} +b^{n}_{j+1}U_{j}^{n}+c^{n}_{j+1}U_{j+2}^{n} )
+ ∆tn

(U_{j}^{n})^{1+α}−(U_{j+1}^{n} )^{1+α}

≥(a^{n}_{j} +b^{n}_{j})U_{j}^{n}+c^{n}_{j}U_{j+1}^{n} −b^{n}_{j+1}U_{j}^{n}−(a^{n}_{j+1}+c^{n}_{j+1})U_{j+1}^{n}

= (a^{n}_{j} +b^{n}_{j} −b^{n}_{j+1})U_{j}^{n}−(a^{n}_{j+1}+c^{n}_{j+1}−c^{n}_{j})U_{j+1}^{n}

=

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

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

(U_{j}^{n}−U_{j+1}^{n} ).

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}^{n+1}≥U_{j+1}^{n+1}for1≤j≤J0−1. Forj=J0, note that
U_{J}^{n+1}

0 −U_{J}^{n+1}

0+1= (a^{n}_{J}_{0}U_{J}^{n}_{0}+b^{n}_{J}_{0}U_{J}^{n}_{0}_{−1}+c^{n}_{J}_{0}U_{J}^{n}_{0}_{+1})

−

"

(1−2λn)U_{J}^{n}_{0}_{+1}+λn

1− N−1 2(J0+ 1)

U_{J}^{n}_{0}

+λ_{n}

1 + N−1 2(J0+ 1)

U_{J}^{n}

0+2

#

+ ∆tn

(U_{J}^{n}_{0})^{1+α}−(U_{J}^{n}_{0}_{+1})^{1+α}

≥

a^{n}_{J}_{0}+b^{n}_{J}_{0}−λn

1− N−1
2(J_{0}+ 1)

U_{J}^{n}_{0}

−

(1−2λn) +λn

1 + N−1 2(J0+ 1)

−c^{n}_{J}_{0}

U_{J}^{n}_{0}_{+1}

=

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

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

(U_{J}^{n}_{0}−U_{J}^{n}_{0}_{+1}).

NowU_{J}^{n+1}

0 ≥U_{J}^{n+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σ= _{2N}^{3} , instead ofλ≤ _{3N}^{1} , it suffices to assumeλ≤ _{2N}^{1} 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}^{1} andσ≤ _{2N}^{3} . Let{U_{j}^{n}}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}^{n}−u(tn, rj)|< C(∆r)^{2},
fortn ≤T0.

Proof.Under the stability conditionsλ≤_{2N}^{1} andσ≤ _{2N}^{3} , the proof is essentially well
known. We give a sketch of it for the sake of convenience. Sete^{n}_{j} =U_{j}^{n}−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),

e^{n+1}_{0} = (1−2N λn)e^{n}_{0}+ 2N λne^{n}_{1} + ∆tn

(U_{0}^{n})^{1+α}−(u(tn,0))^{1+α}+O(∆r)^{2}
,
e^{n+1}_{j} =

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

e^{n}_{j} +λn

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

e^{n}_{j−1}
+λn

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

e^{n}_{j+1}
+ ∆tn

(U_{j}^{n})^{1+α}−(u(tn, rj))^{1+α}+O(∆r)^{2}
,
for1≤j ≤J0, and forJ0+ 1≤j≤J−1,

e^{n+1}_{j} = (1−2λ_{n})e^{n}_{j} +λ_{n}

1−N−1 2j

e^{n}_{j−1}+λ_{n}

1 +N−1 2j

e^{n}_{j+1}
+ ∆t_{n}

(U_{j}^{n})^{1+α}−(u(t_{n}, r_{j}))^{1+α}+O(∆r)^{2}
.
LetE^{n} = max

0≤j≤J|e^{n}_{j}|. Since all the coefficients are nonnegative under the stability assumptions
λ≤_{2N}^{1} andσ≤_{2N}^{3} , it follows that

e^{n+1}_{j}

≤E^{n}+ ∆tn(1 +α) (u(tn, rj) +|E^{n}|)^{α}E^{n}+ ∆tnO (∆r)^{2}
.
Hence,

E^{n+1}≤(1 + ¯C∆t_{n})E^{n}+ ∆t_{n}O (∆r)^{2}
,

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

This inequality is then used to show thatE^{n}≤1for allnif the initial error is small enough.

Then it is used once more to conclude thatE^{n}=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 defineW_{j}^{n}by

(4.1) W_{j}^{n} =δ^{2}U_{j}^{n}+ (1−a)(U_{j}^{n})^{1+α} (∀j= 0, . . . , J−1),
where

δ^{2}U_{0}^{n}= 2NU_{1}^{n}−U_{0}^{n}
(∆r)^{2} ,
δ^{2}U_{j}^{n}= U_{j−1}^{n} −2U_{j}^{n}+U_{j}^{n+1}

(∆r)^{2}

+ 2j

2j+ 1−2σ N−1

rj

σU_{j}^{n}−U_{j−1}^{n}

∆r + (1−σ)U_{j}^{n+1}−U_{j}^{n}

∆r

!

(1≤j≤J0) and

δ^{2}U_{j}^{n}= U_{j−1}^{n} −2U_{j}^{n}+U_{j+1}^{n}

(∆r)^{2} +N−1
r_{j}

U_{j+1}^{n} −U_{j−1}^{n}

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 (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}^{n}can
also be written as

(4.2) W_{j}^{n}= U_{j}^{n+1}−U_{j}^{n}

∆t_{n} −a(U_{j}^{n})^{1+α}.

In view of (4.2) and the Dirichlet boundary condition (3.8), we thus setW_{J}^{n}= 0 (n≥0)for
the sake of convenience.

REMARK4.2. Since (H) and (4.2) imply that
U_{j}^{1}−U_{j}^{0}

∆t0

≥a(U_{j}^{0})^{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}^{n}}be the solution of (3.5)–(3.8). Letλ≤ _{3N}^{1} andσ≤_{2N}^{3} . Assume
that there existsa∈(0,1)such that(H)holds. ThenW_{j}^{n} ≥0, for allj= 0, . . . , J −1and
n≥0. In particular, we have

(4.4) kU^{n+1}k_{∞}− kU^{n}k_{∞}

∆t_{n} = U_{0}^{n+1}−U_{0}^{n}

∆t_{n} ≥a(U_{0}^{n})^{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}^{n} ≥0for allj = 0, . . . , J −1.

By (4.2) and Lemma3.3, we have

U_{j}^{n+1}−U_{j}^{n}

∆tn

≥a(U_{j}^{n})^{1+α}≥0,
which implies thatU_{j}^{n+1}≥U_{j}^{n}for allj. SetV_{j}^{n} :=^{U}

n+1
j −U_{j}^{n}

∆t_{n} . We then obtain by (4.1),
W_{j}^{n+1}−W_{j}^{n}

∆tn

= δ^{2}U_{j}^{n+1}−δ^{2}U_{j}^{n}

∆tn

+ (1−a)(U_{j}^{n+1})^{1+α}−(U_{j}^{n})^{1+α}

∆tn

=δ^{2}V_{j}^{n}+ (1−a)(U_{j}^{n+1})^{1+α}−(U_{j}^{n})^{1+α}

∆t_{n} ,

and, by (4.2),δ^{2}W_{j}^{n}=δ^{2}V_{j}^{n}−aδ^{2}

(U_{j}^{n})^{1+α}

. Therefore, we have
W_{j}^{n+1}−W_{j}^{n}

∆tn

−δ^{2}W_{j}^{n}= (1−a)(U_{j}^{n+1})^{1+α}−(U_{j}^{n})^{1+α}

∆tn

+aδ^{2}

(U_{j}^{n})^{1+α}

for allj = 0, . . . , J −1. SinceU_{j}^{n+1}≥U_{j}^{n}and0≤U_{j+1}^{n} ≤U_{j}^{n} (∀j= 0, . . . , J−1), we
obtain the following four inequalities:

(U_{j}^{n+1})^{1+α}−(U_{j}^{n})^{1+α}

∆t_{n} ≥(1 +α)(U_{j}^{n})^{α}V_{j}^{n} (0≤ ∀j≤J−1),
δ^{2}

(U_{0}^{n})^{1+α}

≥(1 +α)(U_{0}^{n})^{α}δ^{2}U_{0}^{n},
δ^{2}

(U_{j}^{n})^{1+α}

= 1

(∆r)^{2}

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

(U_{j+1}^{n} )^{1+α}−(U_{j}^{n})^{1+α}

−

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

(U_{j}^{n})^{1+α}−(U_{j−1}^{n} )^{1+α}

≥(1 +α)(U_{j}^{n})^{α}δ^{2}U_{j}^{n} (1≤ ∀j≤J_{0}),
δ^{2}

(U_{j}^{n})^{1+α}

= 1

(∆r)^{2}

1 +N−1 2j

(U_{j+1}^{n} )^{1+α}−(U_{j}^{n})^{1+α}

−

1−N−1 2j

(U_{j}^{n})^{1+α}−(U_{j−1}^{n} )^{1+α}

≥(1 +α)(U_{j}^{n})^{α}δ^{2}U_{j}^{n} (J0+ 1≤ ∀j ≤J−1).

These inequalities imply
W_{j}^{n+1}−W_{j}^{n}

∆tn

−δ^{2}W_{j}^{n} ≥(1−a)(1 +α)(U_{j}^{n})^{α}V_{j}^{n}+a(1 +α)(U_{j}^{n})^{α}δ^{2}U_{j}^{n}

= (1 +α)(U_{j}^{n})^{α}W_{j}^{n}≥0

for allj = 0, . . . , J −1. Then the stability assumptionsλ≤ _{3N}^{1} andσ≤ _{2N}^{3} prove that
W_{j}^{n+1} ≥ 0. Now (4.4) is a direct consequence of Lemma3.4and the fact thatW_{0}^{n} ≥ 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}^{n}}be the solution of (3.5)–(3.8)in which∆tnis given by(3.9). Letλ≤ _{3N}^{1} andσ≤_{2N}^{3} .
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.4thatU_{0}^{n} =kU^{n}k_{∞}. Then the results in [5,14,25], namely blow-up occurs only at
r= 0, suggest thatU_{j}^{n}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}^{1} andσ≤ _{2N}^{3} , then a blow-up solution{U_{j}^{n}}satisfies

n→∞lim

kU^{n}k∞

kU^{n+1}k_{∞} = 1
1 +cτ.

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

Therefore, we have by Lemma3.4,
kU^{n+1}k_{∞}

kU^{n}k_{∞} = U_{0}^{n+1}

U_{0}^{n} = 1−2N λn+ 2N λn

U_{1}^{n}

U_{0}^{n} +cτ →1 +cτ as n→ ∞.

Note here that0< U_{1}^{n}/U_{0}^{n}≤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}^{1} andσ≤ _{2N}^{3} , thenU_{1}^{n}/U_{0}^{n}tends to0asn→ ∞.

Proof. Note again that we are concerned with only thosen’s whenU_{0}^{n} =kU^{n}k_{∞} is
sufficiently large, whence∆tn=cτ(U_{0}^{n})^{−α}. 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_{1}^{n+1}

U_{0}^{n+1} ={1−(σ1+σ2)λn}U_{1}^{n}+λnσ1U_{0}^{n}+λnσ2U_{2}^{n}+ ∆tn(U_{1}^{n})^{1+α}
(1−2N λ_{n})U_{0}^{n}+ 2N λ_{n}U_{1}^{n}+ ∆t_{n}(U_{0}^{n})^{1+α} .
Setan=U_{1}^{n}/U_{0}^{n}. Then we have by monotonicity

an+1=

{1−(σ1+σ2)λn}an+λnσ1+λnσ2
U_{2}^{n}

U_{0}^{n} +a^{1+α}_{n} ∆tn(U_{0}^{n})^{α}
1−2N λ_{n}+ 2N λ_{n}a_{n}+ ∆t_{n}(U_{0}^{n})^{α}

≤a_{n}−λ_{n}σ_{1}a_{n}+λ_{n}σ_{1}+a^{1+α}_{n} ∆t_{n}(U_{0}^{n})^{α}
1−2N λn+ 2N λnan+ ∆tn(U_{0}^{n})^{α} .
(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_{n}_{k} < κ and a_{n}_{k}_{+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}(1−a_{n}) + (a_{n})^{1+α}∆t_{n}(U_{0}^{n})^{α}
an−2N λnan(1−an) +an∆tn(U_{0}^{n})^{α} =: 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)^{−2}σ1+ 2N an(∆r)^{−2}−an−a^{1+α}_{n}
1−a_{n} (U_{0}^{n})^{α}

≤(1−an)∆tn[(∆r)^{−2}(2N+σ1)−αan(U_{0}^{n})^{α}].

(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)^{−2}(2N+σ_{1})−a_{n}(U_{0}^{n})^{α}].

Ifliman =a= 1, thenan(U_{0}^{n})^{α}→ ∞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τ)^{−1} σ= _{2N}^{3} ,
(1 +cτ)^{−}^{min{α,1}} σ < _{2N}^{3} .

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

a_{n+1}≥ (1−σ_{2}λ_{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 λ_{n}a_{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}^{3} 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}^{n}}be a blow-up solution of (3.5)–(3.8), and letλ ≤ _{2N}^{1} and
σ= _{2N}^{3} . Then, for allα >0, we haveU_{0}^{n} → ∞asn→ ∞, whileU_{j}^{n}remain bounded for
all1≤j≤J.