June 29, 2004 15:47 WSPC/INSTRUCTION FILE bgr

FULLY DISCRETE ADAPTIVE METHODS FOR A BLOW-UP PROBLEM

CRISTINA BRÄNDLE

Depto. de Matemáticas, U. Autónoma de Madrid, 28049 Madrid, Spain.

                           cristina.brandle@uam.es

PABLO GROISMAN

Inst. de Cálculo, FCEyN., Universidad de Buenos Aires, (1428) Buenos Aires, Argentina.

                              pgroisma@dm.uba.ar

JULIO D. ROSSI

Dpto. de Matemática, FCEyN, Universidad de Buenos Aires, (1428) Buenos Aires, Argentina.

                               jrossi@dm.uba.ar

Received (Day Month Year) Revised (Day Month Year) Communicated by (xxxxxxxxxx)

We present adaptive procedures in space and time for the numerical study of positive solutions to the following problem,

\{ \begin{matrix}u_{t}(x, t) = (u^{m})_{xx}(x, t) & (x, t) ∈ (0, 1) × [0, T), \\ (u^{m})_{x}(0, t) = 0 & t ∈ [0, T), \\ (u^{m})_{x}(1, t) = u^{p}(1, t) & t ∈ [0, T), \\ u(x, 0) = u_{0}(x) & x ∈ (0, 1), \end{matrix}

with p, m > 0. We describe how to perform adaptive methods in order to reproduce the exact asymptotic behavior (the blow-up rate and the blow-up set) of the continuous problem.

Keywords: numerical blow-up; nonlinear diffusion; nonlinear boundary conditions; adaptive schemes.

AMS Subject Classification: 35B40, 35K20, 65M50 2 Brändle-Groisman-Rossi

1. Introduction

We are interested in developing fully discrete adaptive numerical approximations for the following problem,

\{ \begin{matrix}u_{t}(x, t) = (u^{m})_{xx}(x, t) & (x, t) ∈ (0, 1) × [0, T), \\ (u^{m})_{x}(0, t) = 0 & t ∈ [0, T), \\ (u^{m})_{x}(1, t) = u^{p}(1, t) & t ∈ [0, T), \\ u(x, 0) = u_{0}(x) & x ∈ (0, 1). \end{matrix}  (1.1)

In order to have a regular solution, we assume that u_{0} is positive and compatible with the boundary condition. Parabolic problems with nonlinear boundary conditions like (1.1) appear in several branches of applied mathematics. They have been used to model, for example, chemical reactions, heat transfer or population dynamics. See ^{23} and the references therein.

The constant T denotes the maximal time of existence of the solution. If T is infinite we say that the solution is global. However, it is known that, for certain choices of the parameters p and m, the solution may become unbounded in finite time, which means that this solution is defined for t ∈ [0, T), with

lim sup_{t↗T} ∥u(·, t)∥_{L^{∞}} = ∞.

This phenomena is known as blow-up, ^{16}.

Let us summarize the known results concerning solutions of (1.1).

If m > 1 the equation is known as porous medium equation. It is proved in 16 that for every positive initial data the solution blows up if and only if p > 1. Moreover the blow-up rate for increasing in space initial data is given by,

• ∥u(·, t)∥_{L^{∞}} ∼ (T − t)^{− \frac{1}{p−1}} if 1 < p ≤ m,
• ∥u(·, t)∥_{L^{∞}} ∼ (T − t)^{− \frac{1}{2p−m−1}} if p > m.

Here and throughout the paper f(t) ∼ g(t) means that there exist two positive constants c, C, such that cg(t) ≤ f(t) ≤ Cg(t).

For this kind of problems it is important to study the set of points where the solution becomes unbounded. To this end the blow-up set is defined as

B(u) = {x : there exist x_{n} → x, t_{n} ↗ T, such that u(x_{n}, t_{n}) → ∞},

In ^{10,?,?} it is proved that the blow-up set for solutions to (1.1) is given by

• B(u) = [0, 1] if 1 < p ≤ m (global blow-up),

• B(u) = {1} if p > m (single point blow-up).

In the case 0 < m ≤ 1 (known as fast diffusion equation if 0 < m < 1 or heat equation for m = 1) the existence of blow-up solutions depends on m and p. In fact

Adaptive fully discrete methods 3

if p ≤ \frac{m+1}{2} every solution is global (see ^{16}) and every positive solution blows up if p > \frac{m+1}{2}. In this case the blow-up rate is given by

• ∥u(·, t)∥_{L^{∞}} ∼ (T − t)^{− \frac{1}{2p−m−1}} ,

and the blow-up set is

• B(u) = {1} (single point blow-up).

It is remarkable that the blow-up condition depends strongly on the diffusivity, through the exponent m. Indeed, the condition changes across the critical value m = 1, and, if m > 1 the properties of the solution change also across the value p = m.

Numerical approximations for (1.1) have been already studied in ^{13}, where the authors deal with a semidiscrete numerical approximation in a fixed mesh of size h, keeping the time variable t continuos. They compare the behavior of the continuous problem with the discrete one and observed that significant differences appear. For example, they proved that positive solutions of the numerical problem blow up if and only if p > 1. Hence, in the case 0 < m < 1, if \frac{m+1}{2} < p ≤ 1 the continuos solutions blow up in finite time, while the numerical approximations are globally defined. Moreover, it is shown that the blow-up rate for the numerical approximation, U, is given by

∥U(t)∥_{∞} ∼ (T_{h} − t)^{− \frac{1}{p−1}} .

Therefore the blow-up rate does not coincide with the expected for the continuous problem when p > m. In addition, it is shown that the numerical blow-up set is [1 − Lh, 1]. The constant L depends only on p and m, hence the numerical blow-up sets concentrates in a neighborhood of B(u) as h → 0.

In view of these facts, we conclude that usual methods with fixed grid are not well suited in order to reproduce the asymptotic behavior of problem (1.1): the blow-up cases do not coincide and even if both, the discrete and the continuous solutions blow up, the blow-up rates can be different.

The aim of this paper is to combine ideas from ^{2} and ^{14} to obtain fully discrete adaptive schemes that reproduce the blow-up cases, blow-up rates and blow-up sets. We deal with two different improvements of the previous works simultaneously. On the one hand we extend the adaptive in space ideas of ^{14} to the case of nonlinear diffusion. On the other hand we use the ideas from^{2} to develop a method to integrate the time variable adequately.

We use two different methods in order to define the mesh refinement: the first one adds points near the boundary when the numerical solution becomes large. The second method moves the last point near the boundary when the solution becomes large. This procedure is inspired in the moving mesh algorithms developed in ^{6,7,8,18} combined with the time step adaptive procedures of ^{2}. In our case we

4 Brändle-Groisman-Rossi

take advantage of the a priori knowledge of the spatial location of the singularity at x = 1 and, instead of moving the whole mesh, we concentrate only the last point near the boundary, leaving the rest of the mesh fixed. This allows us to use a unified approach to analyze rigorously both schemes simultaneously. One advantage of this moving procedure is keeping the size of the problem to be solved constant in time.

Since we are dealing with totally discrete schemes we also need to determine how to integrate in time. First of all we present an explicit Euler method with adaptive time steps. The implicit Euler scheme is briefly mentioned. We also describe how to discretize using second order Runge-Kutta scheme. These discretizations may be combined with any of the two mentioned mesh refinements (adding points or moving point) to get new numerical methods. Higher order schemes can also be considered for the time discretization.

For previous work on numerical approximations of solutions with blow-up we refer to ^{1,2,3,5,8,11,12,20,21,22} the survey ^{4} and references therein. In ^{2,5,8,14}, adaptive numerical methods are described for problems with linear diffusion, that is m = 1.

To state our main result we need to introduce some notation. Let u_{h,λ} be the numerical solution of (1.1). We use h to denote the size of the spatial mesh and λ is the parameter used in the time discretization. We will say that u_{h,λ} blows up if there exists a finite time T_{h,λ} such that

lim_{t^{j}↗T_{h,λ}} ∥u_{h,λ}(·, t^{j})∥_{∞} = +∞.

In this case, we will call T_{h,λ} the numerical blow-up time. The main theorem of the paper states that, for an explicit Euler scheme, the numerical solution, u_{h,λ}, reproduces the properties of the continuous blowing up solution, u.

Theorem 1.1. Let u be an increasing in x, smooth solution of (1.1), that blows up at time T, and u_{h,λ} the numerical approximation. Then

(i) for every τ > 0, the numerical solution u_{h,λ} converges to the continuous one uniformly in [0, 1]×[0, T −τ]. In fact, there exists a constant C = C(τ) such that

∥u − u_{h,λ}∥_{L^{∞}([0,1]×[0,T −τ])} ≤ C(h^{2} + λ).

(ii) The numerical solution u_{h,λ} blows up if and only if p > 1 when m > 1 and

if and only if p > (m + 1)/2 when m ≤ 1. (iii) The numerical blow-up rate is given by

lim_{t^{j}↗T_{h,λ}}(T_{h,λ} − t^{j})^{\frac{1}{(2p−m−1)}} ∥u_{h,λ}(·, t)∥_{∞} = Γ_{1}, if p > m,

lim_{t^{j}↗T_{h,λ}}(T_{h,λ} − t^{j})^{\frac{1}{(p−1)}} ∥u_{h,λ}(·, t)∥_{∞} = Γ_{2}, if p ≤ m.

Adaptive fully discrete methods 5

(iv) The numerical blow-up time, T_{h,λ}, converges to the continuous one when h and λ tend to 0,

lim_{(h,λ)→(0,0)} |T_{h,λ} − T | = 0.

Moreover, when m = 1, there exist α > 0 and C > 0 such that

|T_{h,λ} − T | ≤ C(h^{2} + λ)^{α}.

(v) Finally, the numerical blow-up set is given by

B(u_{h,λ}) = \{ \begin{matrix}\{1\} & if p > m, \\ [0, 1] & if p ≤ m. \end{matrix}

If higher order methods are considered the conclusions of the theorem still holds, but the order of convergence increases.

Organization of the paper: In Section 2 we describe the numerical adaptive procedures. Section 3 is devoted to the proof of the main result, Theorem 1.1. To begin the analysis we prove that numerical approximations converge uniformly in sets of the form [0, 1] × [0, T − τ]. Next we prove that the scheme reproduces the blow-up rate and set. We also prove that the numerical blow-up time converges to the continuous one.

2. Adaptive numerical schemes

To obtain a fully discrete adaptive scheme we proceed in three steps:

(1) space discretization;

(2) time discretization with adaptive time steps;

(3) mesh refinement as time evolves.

2.1. Space discretization.

For the space discretization we propose piecewise linear finite elements with mass lumping. Consider a partition (that can be non-uniform), x_{1}, ..., x_{N+1} of [0, 1] of size h (h = max(x_{i} − x_{i−1})) and its associated standard piecewise linear finite element space V_{h}. Let {φ_{j}}_{1≤j≤N} be the usual Lagrange basis of V_{h}. We denote by u_{h}(t) the semidiscrete approximation obtained by the finite element method with mass lumping (that in this case coincides with a finite differences discretization) and U(t) = (u_{1}(t), …, u_{N+1}(t)) the values of the numerical approximation at the nodes x_{k}, at time t. Writing,

u_{h}(t) = ∑_{j=1}^{N} u_{j}(t)φ_{j},

6 Brändle-Groisman-Rossi

a simple computation shows that U(t) satisfies the following system of ordinary differential equations (see ^{9}):

\{ \begin{matrix}MU′(t) = −AU^{m}(t) + BU^{p}(t), \\ U(0) = u_{0}^{I}. \end{matrix}  (2.1)

Here M is the mass matrix obtained with lumping, A is the stiffness matrix and u_{0}^{I} is the Lagrange interpolation of the initial datum, u_{0}.

2.2. Time discretization.

Several methods can be used to discretize the time variable. Nevertheless, we will only focus on explicit and implicit adaptive Euler schemes and a second order explicit Runge-Kutta method. Higher order Runge-Kutta methods can be considered as well. Euler methods. If we apply an explicit Euler method to (2.1), we get

\{ \begin{matrix}MU^{j+1} = MU^{j} − τ_{j}A(U^{j})^{m} + τ_{j}B(U^{j})^{p}, \\ U^{0} = u_{0}^{I}, \end{matrix}

where U^{j} = (u_{1}^{j}, …, u_{N+1}^{j}) is the discretization of U(t) at time t^{j}.

The time steps are selected, following ^{2}, according to the rule

\begin{matrix}τ_{j}(u_{N+1}^{j})^{2p−m} = λ, & if p > m, \\ τ_{j}(u_{N+1}^{j})^{p} = λ, & if p ≤ m, \end{matrix}  (2.2)

where λ is the parameter of the method that has to be selected small in order to verify the usual stability condition of explicit schemes. Remark that the time step τ_{j} becomes smaller as the solution increases.

In a similar way we can consider an implicit method by evaluating the discrete laplacian (A) at time t^{j+1}. That is

\{ \begin{matrix}MU^{j+1} = MU^{j} − τ_{j}A(U^{j+1})^{m} + τ_{j}B(U^{j})^{p}, \\ U^{0} = u_{0}^{I}. \end{matrix}

Note that the scheme is not totally implicit since the nonlinear boundary condition u^{p} is evaluated at time t^{j} while the discrete laplacian is evaluated a time t^{j+1}. This discretization has the advantage that, if m = 1, the explicit evaluation of (U^{j})^{p} avoids the problem of solving a nonlinear system in each step. For m ≠ 1 one can consider either B(U^{j})^{p} or B(U^{j+1})^{p}, since the nonlinearity is originated by the diffusion term. The time steps τ_{j} are chosen in the same way as before, but in this case no stability restrictions have to be imposed.

However, we will show that if U^{j} is large enough, the required stability restriction for the explicit method is always verified. Therefore any implicit method can be

Adaptive fully discrete methods 7

used until the stability condition is satisfied, and from that moment on the explicit scheme is recommended. Runge-Kutta methods. To produce higher order methods we proceed following ideas from^{2}. We introduce a new time variable s = s(t) for the semidiscrete problem (2.1), such that Û(s) = U(t(s)). We select s in such a way that Û(s) is globally defined. Hence any method with fixed time step can be applied to Û(s). The time rescaling is given by

\{ \begin{matrix}\frac{dt}{ds} = \frac{1}{u_{N+1}^{η}(t)}, \\ t(0) = 0. \end{matrix}  (2.3)

The parameter η will be chosen below in terms of p and m. Observe that since U is positive, t(s) defines a one to one function. In terms of (Û(s), s), the system (2.1) takes the form

\{ \begin{matrix}MÛ′(s) = F (Û), \\ Û(0) = u_{0}^{I}, \end{matrix}  (2.4)

where F is given by

F (Û) = \frac{1}{û_{N+1}^{η}(s)} (−AÛ^{m}(s) + BÛ^{p}(s)) . (2.5)

We can approximate (2.3)–(2.4) by any method with fixed step λ. This allows us to define adaptive in time procedures for (2.1). For example, if we apply an explicit Euler method we obtain

\begin{matrix}\{ & \begin{matrix}\frac{t^{j+1} − t^{j}}{λ} = \frac{1}{(u_{N+1}^{j})^{η}} , \\ t(0) = 0, \end{matrix} \\ \{ & \begin{matrix}\frac{MÛ^{j+1} − MÛ^{j}}{λ} = \frac{1}{(û_{N+1}^{j})^{η}} (−A(Û^{j})^{m} + B(Û^{j})^{p}) , \\ Û(0) = u_{0}^{I}. \end{matrix} \end{matrix}

Calling τ_{j} = t^{j+1} − t^{j}, these equations can be rewritten in terms of U^{j} = Û^{j} and τ_{j} as

\{ \begin{matrix}MU^{j+1} − MU^{j} = τ_{j} (−A(U^{j})^{m} + B(U^{j})^{p}) , \\ τ_{j} = \frac{λ}{(u_{N+1}^{j})^{η}} . \end{matrix}

Observe that this coincides with the adaptive explicit Euler method previously described. So, as before we take η as

η = \{ \begin{matrix}2p − m & if p > m, \\ p & if p ≤ m. \end{matrix}

8 Brändle-Groisman-Rossi

If we apply the half step second order R-K method to (2.3)–(2.4) we obtain

\begin{matrix}\{ & \begin{matrix}\frac{t^{j+1} − t^{j}}{λ} = \frac{1}{(û_{N+1}^{j} + \frac{λ}{2}F (Û^{j})_{N+1})^{η}} , \\ t(0) = 0, \end{matrix} \\ \{ & \begin{matrix}\frac{MÛ^{j+1} − MÛ^{j}}{λ} = F (Û^{j} + \frac{λ}{2}F (Û^{j})), \\ Û(0) = u_{0}^{I}. \end{matrix} \end{matrix}

The vector function F is defined as in (2.5). This system can be rewritten in the original variables as

\{ \begin{matrix}MU^{j+1} − MU^{j} = τ_{j} (−A(U^{j} + \frac{λ}{2}F (U^{j}))^{m} + B((U^{j} + \frac{λ}{2}F (U^{j}))^{p}) , \\ τ_{j} = \frac{λ}{(u_{N+1}^{j} + \frac{λ}{2}F (U^{j})_{N+1})^{η}} . \end{matrix}

Using these ideas any Runge-Kutta method can be used to approximate (2.4). This provides higher order adaptive in time procedures.

2.3. Mesh refinement.

Once the method to integrate the time variable is selected, it remains to show how to perform the spatial mesh refinement. To this end let us pay special attention to the numerical solution at the point x_{N+1} = 1, where the continuous solution develops the singularity.

The numerical mesh refinement that we describe is based on the scale invariance of the problem in the half line with the nonlinear boundary condition placed at x = 0, −(u^{m})_{x}(0, t) = u^{p}(0, t). We mean the following: there exists a self-similar blow-up solution in the half-line of the form u_{S}(x, t) = (T −t)^{−α}ψ(ξ), ξ = x(T −t)^{−β}, where α = 1/(2p − m − 1), β = (p − m)/(2p − m − 1). This solution u_{S}(x, t) gives the behavior near the blow-up time T for solutions of (1.1). Our numerical schemes use this fact to add or move points of the mesh near x = 1 trying to reproduce the scaling invariance in the half-line. Adding points. For the sake of simplicity we will consider the explicit Euler method for time discretization to describe the spatial refinement. We will describe how to proceed just in the case p > m since for this choice of parameters fixed mesh methods do not reproduce the asymptotic behavior. If p ≤ m, as will be shown, the adaptive procedure described here does not refine the mesh and hence the method is a time discretization of the system analyzed in ^{13}.

Writing the equation satisfied by u_{N+1}^{j} explicitly we obtain the following,

u_{N+1}^{j+1} − u_{N+1}^{j} = τ_{j} (\frac{2}{h_{N}^{2}} ((u_{N}^{j})^{m} − (u_{N+1}^{j})^{m}) + \frac{2}{h_{N}} (u_{N+1}^{j})^{p}) ,

Adaptive fully discrete methods 9

where h_{N} = 1 − x_{N}. As we want the numerical blow-up rate to be

u_{N+1}^{j} ∼ (T_{h,λ} − t^{j})^{− \frac{1}{(2p−m−1)}}

we impose the last node u_{N+1}^{j} to satisfy

c_{1}(u_{N+1}^{j})^{2p−m} ≤ \frac{u_{N+1}^{j+1} − u_{N+1}^{j}}{τ_{j}} ≤ c_{2}(u_{N+1}^{j})^{2p−m}. (2.6)

A straightforward calculation shows that this implies the behavior described above. Recalling (2.2), we get

c_{1}λ ≤ u_{N+1}^{j+1} − u_{N+1}^{j} ≤ c_{2}λ. (2.7)

On the other hand, imposing (2.6) is equivalent to,

c_{1}(u_{N+1}^{j})^{2p−m} ≤ \frac{2}{h_{N}^{2}} ((u_{N}^{j})^{m} − (u_{N+1}^{j})^{m}) + \frac{2}{h_{N}} (u_{N+1}^{j})^{p} ≤ c_{2}(u_{N+1}^{j})^{2p−m}. (2.8)

As the scheme with fixed mesh blows up at the last node, so we have that

R(j; h_{N}) = \frac{\frac{2}{h_{N}^{2}} ((u_{N}^{j})^{m} − (u_{N+1}^{j})^{m}) + \frac{2}{h_{N}} (u_{N+1}^{j})^{p}}{(u_{N+1}^{j})^{2p−m}} → 0,

as j increases. Let t^{j_{1}} be the first time such that R(j; h_{N}) ≤ c_{1} (hereafter we will use t_{i} to denote t^{j_{i}}). Namely t_{1} is the first time where R(j_{1}; h_{N}) does not verify the lower bound c_{1}. At that moment, we add a point z between x_{N} and x_{N+1} = 1 to the mesh and give the value to u_{h,λ}(z, t_{1}) such that the slope of the line between (z, (u_{h,λ})^{m}(z, t_{1})) and (1, (u_{h,λ})^{m}(1, t_{1})) is the same as the slope between (x_{N}, (u_{h,λ})^{m}(x_{N}, t_{1})) and (1, (u_{h,λ})^{m}(1, t_{1})). Hence we have a new value for the length of the last interval, [z, 1],

h_{N,1} = 1 − z < h_{N} = 1 − x_{N}.

In other words, we add the new node z at 1 − h_{N,1} in such a way that,

\frac{1}{h_{N}} ((u_{N}^{j_{1}})^{m} − (u_{N+1}^{j_{1}})^{m}) = \frac{1}{h_{N,1}} ((u_{z}^{j_{1}})^{m} − (u_{N+1}^{j_{1}})^{m}).

10 Brändle-Groisman-Rossi

Therefore, we obtain that R(j_{1}; h_{N,1}) satisfies

R(j_{1}; h_{N,1}) = \frac{\frac{2}{h_{N,1}^{2}} ((u_{z}^{j_{1}})^{m} − (u_{N+1}^{j_{1}})^{m}) + \frac{2}{h_{N,1}} (u_{N+1}^{j_{1}})^{p}}{(u_{N+1}^{j_{1}})^{2p−m}}

= \frac{\frac{2}{h_{N,1}} ((u_{z}^{j_{1}})^{m} − (u_{N+1}^{j_{1}})^{m}) + 2(u_{N+1}^{j_{1}})^{p}}{h_{N,1} (u_{N+1}^{j_{1}})^{2p−m}}

= \frac{\frac{2}{h_{N}} ((u_{N}^{j_{1}})^{m} − (u_{N+1}^{j_{1}})^{m}) + 2(u_{N+1}^{j_{1}})^{p}}{h_{N,1} (u_{N+1}^{j_{1}})^{2p−m}}

= \frac{h_{N}}{h_{N,1}} R(j_{1}; h_{N}) > c_{1}.

With this new mesh x_{1}, ..., x_{N}, z, x_{N+1} we have that

R(j_{1}; h_{N,1}) > c_{1},

and we can continue the method with initial time step τ_{j} in the new mesh. This gives a solution u_{h,λ} that verifies (2.8) in a time interval [t_{1}, t_{2}] where at time t_{2}, the quotient R(j_{2}; h_{N,1}) fails again the lower bound c_{1}. At that time we have to add another point in the last interval. As before this increases R(j_{2}; h_{N,2}) and we can continue in a new mesh, made up by the old one plus the point that we have added near the boundary x = 1. This procedure generates an increasing sequence of times t_{i} (the times at which we refine the spatial mesh), an increasing sequence of added points accumulating at x = 1 and a numerical solution u_{h,λ}(x, t) defined for (x, t) ∈ [0, 1] × [0, T_{h,λ}).

Let us point out that, when using this procedure, the system verified by the numerical solution changes whenever we add a new point to the mesh. Indeed, U^{j} is a solution of

\{ \begin{matrix}M_{i}U^{j+1} = M_{i}U^{j} − τ_{j}A_{i}(U^{j})^{m} + τ_{j}B_{i}(U^{j})^{p}, \\ U^{0} = u_{0}^{I}. \end{matrix}  (2.9)

The matrices A_{i}, B_{i} and M_{i} are modified each time the mesh is refined. In fact, as we add a new point to the mesh, an equation is added to the system, and hence we need to add a new row (which takes into account the added point) to the three matrices. The matrix A_{i} is modified in the elements that correspond to the new point. The matrix B_{i} has zeros in all of its rows but the last one, where it collects the boundary condition. Nevertheless, the matrix A_{i} preserves the properties of the original matrix A, i.e.

a_{ii} > 0,  a_{ij} ≤ 0, if j ≠ i  ∑_{k=1}^{N(j)+1} a_{kj} ≥ 0,

Adaptive fully discrete methods 11

where N(j) is the dimension of U^{j}. For simplicity, we will keep the original nomenclature of the nodes, i.e

U^{j} = \overbrace{(u_{1}^{j}, u_{2}^{j}, …, u_{N}^{j}, u_{N+1}^{j})}^{N(j)+1} = (u_{1}^{j}, u_{2}^{j}, …, u_{N(j)}^{j}, u_{N(j)+1}^{j}),

and we will assume that the mesh is fixed and uniform, whenever we can do it. If this is not the case, we will mention it explicitly. Under these conditions, we obtain the following discretization for the original problem (1.1)

\{ \begin{matrix}\frac{u_{1}^{j+1} − u_{1}^{j}}{τ_{j}} & = \frac{2((u_{2}^{j})^{m} − (u_{1}^{j})^{m})}{h^{2}}, \\ \frac{u_{i}^{j+1} − u_{i}^{j}}{τ_{j}} & = \frac{(u_{i−1}^{j})^{m} − 2(u_{i}^{j})^{m} + (u_{i+1}^{j})^{m}}{h^{2}} , 2 ≤ i ≤ N, \\ \frac{u_{N+1}^{j+1} − u_{N+1}^{j}}{τ_{j}} & = \frac{2((u_{N}^{j})^{m} − (u_{N+1}^{j})^{m})}{h^{2}} + \frac{2(u_{N+1}^{j})^{p}}{h} , \end{matrix}  (2.10)

and the initial datum u_{i}^{0} = u_{0}(x_{i}).

It remains to be more concrete on the election of the sequence h_{i} and the constants c_{1}, c_{2}. If we take,

c_{1} = Γ_{1}^{−(2p−m−1)}(2p − m − 1)^{−1},

and h_{i} and c_{2} such that

\frac{h_{i}}{h_{i+1}} → 1, and c_{2}(t_{i}) = \frac{h_{i}}{h_{i+1}} R(j, h_{i}), (2.11)

we get

c_{1} ≤ R(j, h_{i+1}) ≤ \frac{h_{i}}{h_{i+1}} R(j, h_{i}) ∼ c_{1}.

In fact, if h_{i} is chosen such that

h_{i+1}(u_{N+1}^{j})^{p−m} = \frac{2}{c_{1}} − \frac{A}{u_{N+1}^{j}} , (2.12)

12 Brändle-Groisman-Rossi

we have

c_{2} = \frac{h_{i}}{h_{i+1}} R(j_{i}, h_{i}) = \frac{\frac{2}{h_{i+1}} ((u_{z}^{j_{i}})^{m} − (u_{N+1}^{j_{i}})^{m}) + 2(u_{N+1}^{j_{i}})^{p}}{h_{i+1} (u_{N+1}^{j_{i}})^{2p−m}}

= \frac{\frac{2}{h_{i+1}} ((u_{z}^{j_{i}})^{m} − (u_{N+1}^{j_{i}})^{m}) + 2(u_{N+1}^{j_{i}})^{p}}{(u_{N+1}^{j_{i}})^{p} (\frac{2}{c_{1}} − _{\frac{A}{u_{N+1}^{j_{i}}}} )}

= \frac{\frac{2}{h_{i+1}} ((u_{z}^{j_{i}})^{m} − (u_{N+1}^{j_{i}})^{m})}{(u_{N+1}^{j_{i}})^{p} (\frac{2}{c_{1}} − _{\frac{A}{u_{N+1}^{j_{i}}}} )} + c_{1} (\frac{1}{1 − \frac{c_{1}A}{2u_{N+1}^{i}}} )

≤ c_{1} (\frac{1}{1 − \frac{c_{1}A}{2u_{N+1}^{i}}} ) → c_{1}.

Observe that

\frac{τ_{j}}{h_{i+1}} ≤ Cλ\frac{(u_{N+1}^{j})^{p−m}}{(u_{N+1}^{j})^{2p−m}} = \frac{Cλ}{u_{N+1}^{j}} .

Hence usual stability conditions for explicit schemes are verified if u_{N+1}^{j} is large. In addition we will show that this selection of h_{i} makes

h_{i} = 1 − x_{N} ∼ C(T_{h,λ} − t^{j})^{β},

hence the method reproduces both the asymptotic behavior and the spatial structure of the continuous solution near the blow up time.

Observe that if p ≤ m, R(j, h_{N}) is bounded from below and the method keeps the mesh fixed. Moving point. An alternative procedure to refine the mesh near the singularity is to move a point towards x_{n+1} = 1. This method is inspired in the moving mesh methods described in ^{8}.

We begin with a mesh composed by two types of nodes, a uniform mesh of size h = 1/N and a node placed between x_{N} = (N − 1)h and x_{N+1} = 1, that we are going to move when appropriate. Let us call 0 = x_{1} < .... < x_{N} < x_{N+1} = 1 the fixed mesh and w the moving node, x_{N} < w < x_{N+1}.

As before, let us consider an explicit Euler procedure for the time discretization and let us use the mesh composed by the x_{i} and the w together for the space discretization. If we denote by z the numerical solution at the moving node w, we

Adaptive fully discrete methods 13

[figure: x_{N−1}; x_{N}; w; h; h − h_{i}; h_{i}; x_{N+1}; Moving point method]

arrive to the following system

\{ \begin{matrix}\frac{u_{1}^{j+1} − u_{1}^{j}}{τ_{j}} & = \frac{2 ((u_{2}^{j})^{m} − (u_{1}^{j})^{m})}{h^{2}}, \\ \frac{u_{i}^{j+1} − u_{i}^{j}}{τ_{j}} & = \frac{(u_{i−1}^{j})^{m} − 2(u_{i}^{j})^{m} + (u_{i+1}^{j})^{m}}{h^{2}} , 2 ≤ i ≤ N − 1 \\ \frac{u_{N}^{j+1} − u_{N}^{j}}{τ_{j}} & = \frac{2 ((u_{N−1}^{j})^{m}(h − h_{i}) − (u_{N}^{j})^{m}(2h − h_{i}) + (z^{j})^{m}h)}{h(h − h_{i})(2h − h_{i})} \\ \frac{z^{j+1} − z^{j}}{τ_{j}} & = \frac{2 ((u_{N}^{j})^{m}(h_{i}) − (z^{j})^{m}h + (u_{N+1}^{j})^{m}(h − h_{i}))}{(h − h_{i})hh_{i}} \\ \frac{u_{N+1}^{j+1} − u_{N+1}^{j}}{τ_{j}} & = \frac{2 ((z^{j})^{m} − (u_{N+1}^{j})^{m})}{h_{i}^{2}} + \frac{2(u_{N+1}^{j})^{p}}{h_{i}} , \end{matrix}

and the initial data u_{i}^{0} = u_{0}(x_{i}).

The moving point method uses the ideas developed for the adding points method to move the point when R(j, h_{N}) ≤ c_{1}. When the constant c_{1} is reached the node w is moved towards x_{N+1} = 1 choosing the new length h_{i} = 1 − w in the same way as in the adding points method. This new value makes R(j, h_{N,1}) > c_{1} and the method continues until the constant c_{1} is reached again.

Let us remark that the criteria that we use to modify the mesh is the same in the adding points method and in the moving point method. This allows us to make a unified approach in the course of the proofs contained in the following section.

3. Proof of the main results

We may now proceed with the proof of Theorem 1.1. We will do it with full details for the adding points scheme with an explicit Euler method for the time variable and briefly point out the main differences with the other schemes described in Section 2. We state first some preliminary results.

Lemma 3.1. If U^{0} is increasing, then the numerical solution U^{j} is increasing for every j. Hence it satisfies max_{i} u_{i}^{j} = u_{N+1}^{j}, for every j.

14 Brändle-Groisman-Rossi

Proof. As U^{0} is increasing then u_{i+1}^{0} > u_{i}^{0}, let us see that this holds for every j > 0. Assume not, then there exists a first positive j_{0} and an index i such that u_{i+1}^{j_{0}} ≤ u_{i}^{j_{0}}. From the equations satisfied by u_{i+1}^{j} and u_{i}^{j} we have that at j_{0} it holds u_{i+1}^{j_{0}+1} − u_{i+1}^{j_{0}} − (u_{i}^{j_{0}+1} − u_{i}^{j_{0}}) > 0, a contradiction that proves the result.

From now on we will assume that the discrete initial data U^{0} is increasing, u_{i}^{0} ≤ u_{i+1}^{0} for all 1 ≤ i ≤ N + 1.

Definition 3.2. We will say that Z^{j} is a supersolution (resp. a subsolution) of (2.10) if verifies

\{ \begin{matrix}MZ^{j+1} ≥ MZ^{j} − τ_{j}A(Z^{j})^{m} + τ_{j}B(Z^{j})^{p}, (resp. ≤) \\ Z^{0} ≥ z_{0}^{I}, (resp. ≤). \end{matrix}

Lemma 3.3. Let \overline{U}^{j} be a supersolution of (2.10) and \underline{U}^{j} a subsolution, such that \underline{U}^{0} < \overline{U}^{0}. If the time step verifies

τ_{j} < \frac{Ch_{i}^{2}}{∥\overline{U}^{j}∥_{∞}^{m}} , (3.1)

then \underline{U}^{j} < \overline{U}^{j}.

Proof. Let us define Z^{j} = \overline{U}^{j} −\underline{U}^{j}. By an approximation argument we can assume that we have strict inequalities in (2.10). Then Z^{j} is a solution of the following system

\{ \begin{matrix}\frac{MZ^{j+1} − MZ^{j}}{τ_{j}} > −A((\overline{U}^{j})^{m} − (\underline{U}^{j})^{m}) + B((\overline{U}^{j})^{p} − (\underline{U}^{j})^{p}), \\ Z^{0} > 0. \end{matrix}

If the Lemma is false, then there exist a first time t^{j+1} and a node x_{i}, such that z_{i}^{j+1} ≤ 0 and z_{i}^{j} > 0. Then we have,

z_{i}^{j+1} > z_{i}^{j} − \frac{τ_{j}}{m_{i}} ( ∑_{k=1}^{N+1} a_{ki}((\overline{u}_{i}^{j})^{m} − (\underline{u}_{i}^{j})^{m}) + b_{i}((\overline{u}_{i}^{j})^{p} − (\underline{u}_{i}^{j})^{p}))

= z_{i}^{j} (1 − \frac{Cτ_{j}}{m_{i}} a_{ii} ((\overline{u}_{i}^{j})^{m} − (\underline{u}_{i}^{j})^{m}))

+\frac{τ_{j}}{m_{i}} (∑_{k≠i} a_{ki} ((\overline{u}_{i}^{j})^{m} − (\underline{u}_{i}^{j})^{m}) z_{k}^{j} + b_{i}((\overline{u}_{i}^{j})^{p} − (\underline{u}_{i}^{j})^{p}))

> z_{i}^{j} (1 − \frac{Cτ_{j}}{m_{i}} a_{ii}∥\overline{U}^{j}∥_{∞}^{m}) ≥ z_{i}^{j} (1 − \frac{Cτ_{j}}{h_{i}^{2}} ∥\overline{U}^{j}∥_{∞}^{m}) ≥ 0,

a contradiction.

Adaptive fully discrete methods 15

Corollary 3.4. If \overline{U}^{j} is a supersolution of (2.10), such that

∥\overline{U}^{j}∥_{∞} ≤ Cu_{N+1}^{j},

then Lemma 3.3 holds for all j, if p > m. If p ≤ m, the lemma holds as long as (3.1) is valid.

3.1. Convergence of the numerical schemes.

We are now ready to establish an uniform convergence result: for any τ > 0 we want that u_{h,λ} → u (when h, λ → 0) uniformly in [0, 1]×[0, T −τ]. This is a natural requirement since on such an interval the continuous solution is regular.

Lemma 3.5. Let u(x, t) ∈ C^{4,2}([0, 1]×[0, T −τ]) be a positive solution of (1.1) and let u_{h,λ} its discrete approximation obtained by any of the adaptive schemes described in Section 2. Then there exists a constant C depending on τ such that,

∥u − u_{h,λ}∥_{L^{∞}([0,1]×[0,T −τ])} ≤ C(h^{2} + λ).

Proof. In the course of the proof we will consider an uniform mesh. We are able to do this since in the time interval [0, T − τ] both, the numerical and the continuous solution, are bounded if h is small enough.

If we rewrite the system (2.10) in terms of Z = (U^{j})^{m}, we obtain

\{ \begin{matrix}\frac{(z_{1}^{j+1})^{1/m} − (z_{1}^{j})^{1/m}}{τ_{j}} & = \frac{2(z_{2}^{j} − z_{1}^{j})}{h^{2}}, \\ \frac{(z_{i}^{j+1})^{1/m} − (z_{i}^{j})^{1/m}}{τ_{j}} & = \frac{(z_{i−1}^{j} − 2z_{i}^{j} + z_{i+1}^{j})}{h^{2}},  2 ≤ i ≤ N, \\ \frac{(z_{N+1}^{j+1})^{1/m} − (z_{N+1}^{j})^{1/m}}{τ_{j}} & = \frac{2(z_{N}^{j} − z_{N+1}^{j})}{h^{2}} + \frac{2(z_{N+1}^{j})^{p/m}}{h} , \end{matrix}

and the initial datum z_{i}^{0} = (u_{0}(x_{i}))^{m}. Let v_{i}^{j} = u^{m}(x_{i}, t^{j}), with u the solution of the continuous problem (1.1). The error function

e_{i}^{j} = z_{i}^{j} − v_{i}^{j}

16 Brändle-Groisman-Rossi

satisfies, for 2 ≤ i ≤ N,

\frac{1}{m}ξ^{(1−m)/m} \frac{e_{i}^{j+1} − e_{i}^{j}}{τ_{j}}

= \frac{(e_{i−1}^{j} − 2e_{i}^{j} + e_{i+1}^{j})}{h^{2}} − \frac{1 − m}{m^{2}}ζ^{(2−m)/m} \frac{v_{i}^{j+1} − v_{i}^{j}}{τ_{j}} |ξ − η|

+C_{1}(h^{2} + λ)

≤ \frac{(e_{i−1}^{j} − 2e_{i}^{j} + e_{i+1}^{j})}{h^{2}} − \frac{1 − m}{m^{2}}ζ^{(2−m)/m} \frac{v_{i}^{j+1} − v_{i}^{j}}{τ_{j}} (v_{i}^{j+1} − v_{i}^{j})|e_{i}^{j+1}|

+C_{1}(h^{2} + λ),

where ξ is an intermediate value between z_{i}^{j+1} and z_{i}^{j}, η is between v_{i}^{j+1} and v_{i}^{j} and ζ is an intermediate value between ξ and η. Taking into account that there exist constants, c and C, such that c ≤ z_{i}^{j} ≤ C for every j ∈ [0, j_{0}] we have,

\frac{e_{i}^{j+1} − e_{i}^{j}}{τ_{j}} ≤ \frac{C_{1}(e_{i−1}^{j} − 2e_{i}^{j} + e_{i+1}^{j})}{h^{2}} + C_{2}|e_{i}^{j+1}| + C_{3}(h^{2} + λ). (3.2)

Making analogous calculations for the first and the last nodes, we get

\frac{e_{1}^{j+1} − e_{1}^{j}}{τ_{j}} ≤ \frac{C_{1}(e_{2}^{j} − e_{1}^{j})}{h^{2}} + C_{2}|e_{1}^{j+1}| + C_{3}(h^{2} + λ), (3.3)

and

\begin{aligned}\frac{e_{N+1}^{j+1} − e_{N+1}^{j}}{τ_{j}} & ≤ \frac{C_{1}(e_{N}^{j} − e_{N+1}^{j})}{h^{2}} + C_{2}|e_{N+1}^{j+1}| \\ & +\frac{C_{3}}{h} κ^{(p−m)/m}|e_{N+1}^{j}| + C_{4}(h^{2} + λ), \end{aligned}  (3.4)

where κ is an intermediate value between z_{N+1}^{j} and v_{N+1}^{j}, and hence κ is also bounded.

Now we use a comparison argument. The error e_{i}^{j} is a subsolution of a discretization of the problem

\{ \begin{matrix}e_{t} = e_{xx} + e + C(h^{2} + λ) & (x, t) ∈ (0, 1) × [0, T − τ), \\ e_{x}(0, t) = 0 & t ∈ [0, T − τ), \\ e_{x}(1, t) = Ce(1, t) & t ∈ [0, T − τ), \\ e(x, 0) = C(h^{2} + λ) & x ∈ (0, 1). \end{matrix}  (3.5)

Since solutions of problem (3.5) are bounded by

|e(x, t)| ≤ C exp{CT }(h^{2} + λ),

Adaptive fully discrete methods 17

we can construct a supersolution \overline{e}_{i}^{j} of (3.2)–(3.4), the system verified by e_{i}^{j}, by taking \overline{e}_{i}^{j} = e(x_{i}, t^{j}). Therefore we obtain a similar bound for the error function e_{i}^{j}. Arguing similarly with −e_{i}^{j} we obtain

|e_{i}^{j}| ≤ C exp{CT }(h^{2} + λ), t^{j} ∈ (0, T − τ],

and this proves the convergence of the method.

3.2. Numerical blow-up.

We prove now that solutions of the numerical problem (2.10) blow up in the same range of exponents as the solutions of the continuous problem do.

Let us recall, for notational purposes

t^{j} = ∑_{k=0}^{j} τ_{k}.

Observe that since, for each j, the maximum of the numerical solution is attained at the last node, in order to have a numerical solution that blows up in finite time T_{h,λ}, we have to see that the series converge and that u_{N+1}^{j} → ∞ as j increases.

Lemma 3.6. The numerical solution u_{h,λ} blows up if and only if p > 1, when m > 1 and p > (m + 1)/2, when m ≤ 1.

Proof. Since U^{0} is increasing, by Lemma 3.1 we have that U^{j} reaches its maximum at i = N + 1. On the other hand, from (2.7) we obtain u_{N+1}^{j} ∼ λj, and hence

u_{N+1}^{j} → ∞, as j → ∞. (3.6)

We still have to show conditions which imply that the time T_{h,λ} at which (3.6) occurs is finite. Indeed, for p > m, we get

T_{h,λ} = ∑_{k=1}^{∞} τ_{k} = ∑_{k=1}^{∞}(u_{N+1}^{j})^{−(2p−m)}λ ∼ C ∑_{k=1}^{∞}(λj)^{−(2p−m)}λ.

Observe that this series converges or diverges according to 2p−m > 1 or 2p−m ≤ 1. Analogously, for p ≤ m, we get

T_{h,λ} = ∑_{k=1}^{∞} τ_{k} = ∑_{k=1}^{∞}(u_{N+1}^{j})^{−p}λ ∼ C ∑_{k=1}^{∞}(λj)^{−p}λ.

18 Brändle-Groisman-Rossi

3.3. Blow-up rate.

This section deals with the blow-up rate of the numerical solutions. We show that they coincide with the one for the continuous problem.

Lemma 3.7. If p > m (for m > 1) or p > (m + 1)/2 (for m ≤ 1), then every positive solution that blows up at time T_{h,λ} verifies

∥U^{j}∥_{∞} = u_{N+1}^{j} ∼ (T_{h,λ} − t^{j})^{\frac{−1}{2p−m−1}} .

Proof. From Lemma 3.6 we know that u_{h,λ} blows up in finite time T_{h,λ} = ∑_{k=1}^{∞} τ_{k}. In order to prove the blow-up rate we observe that for any k > j, we have

u_{N+1}^{k} ≥ u_{N+1}^{j} + c_{1}λ(k − j).

Therefore,

(T_{h,λ} − t^{j})

= ∑_{k=j+1}^{∞} τ_{k} = ∑_{k=j+1}^{∞} \frac{λ}{(u_{N+1}^{k})^{2p−m}} ≤ ∑_{k=j+1}^{∞} \frac{λ}{(u_{N+1}^{j} + c_{1}λ(k − j))^{2p−m}}

= ∑_{l=1}^{∞} \frac{λ}{(u_{N+1}^{j} + c_{1}λl)^{2p−m}} ≤ ∫_{0}^{∞} \frac{λ}{(u_{N+1}^{j} + c_{1}λs)^{2p−m}} ds

= C ∫_{u_{N+1}^{j}}^{∞} \frac{1}{z^{2p−m}} dz = \frac{C}{(u_{N+1}^{j})^{2p−m−1}} .

Thus,

u_{N+1}^{j} ≤ C(T_{h,λ} − t^{j})^{\frac{−1}{2p−m−1}} .

On the other hand, since u_{N+1}^{k} ≤ u_{N+1}^{j} +c_{2}λ(k −j), using the same ideas as before, we get

u_{N+1}^{j} ≥ C(T_{h,λ} − t^{j})^{\frac{−1}{2p−m−1}} ,

and we conclude

u_{N+1}^{j} ∼ (T_{h,λ} − t^{j})^{\frac{−1}{2p−m−1}} ,

as expected from the description of the method. Moreover, if we take into account that c_{1} and c_{2} were chosen as

c_{1} = Γ_{1}^{−(2p−m−1)}(2p − m − 1)^{−1}, c_{2}(t^{j}) ∼ c_{1},

we have

lim_{j→∞} u_{N+1}^{j}(T_{h,λ} − t^{j})^{\frac{1}{2p−m−1}} → c_{1}.

Adaptive fully discrete methods 19

Observe that from the blow-up rate we have

h_{i} ∼ (u_{N+1}^{j})^{−p+m} ∼ C(T_{h,λ} − t^{j})^{\frac{p−m}{2p−m−1}} = C(T_{h,λ} − t^{j})^{β},

and, as mentioned before, the selfsimilar behavior is reproduced.

Lemma 3.8. If p ≤ m, then every positive solution that blows up at time T_{h,λ} verifies

∥U^{j}∥_{∞} = u_{N+1}^{j} ∼ (T_{h,λ} − t^{j})^{− \frac{1}{p−1}} .

Proof. It is similar to the previous one, using τ_{j}(u_{N+1}^{j})^{p} = λ.

3.4. Convergence of the numerical blow-up times.

We gather now the proof of convergence of the numerical blow-up times, T_{h,λ}, to the continuous one, T as (h, λ) → (0, 0). When m = 1 we state a stronger result, that gives a bound for |T_{h,λ} − T | in terms of h and λ.

Lemma 3.9. If the numerical solution u_{h,λ} converges uniformly to the continuous one, u, in sets of the form [0, 1] × [0, T − τ] and verifies the blow-up rates described in lemmas 3.7 and 3.8, then the numerical blow-up times converge to the continuous one:

lim_{(h,λ)→(0,0)} T_{h,λ} = T.

Observe that the schemes we are considering verify the assumptions of this lemma. In fact, in Lemma 3.5 we proved the uniform convergence and from Lemmas 3.7 and 3.8 we get the blow-up rate. Recall that in the proof of convergence we used a comparison result that needs a stability restriction on λ when the method is explicit.

Proof. From the blow-up rate we have that,

T_{h,λ} − t^{j} ≤ \frac{C}{(u_{N+1}^{j})^{γ}} ,

where γ = 2p − m − 1 or γ = p − 1 according to p > m or p ≤ m. Given ε > 0, we can choose K large enough in order to have

\frac{C}{K^{γ}} < \frac{ε}{2}.

20 Brändle-Groisman-Rossi

Let τ < ε/4 such that the continuous solution verifies u(1, T − 2τ) > K. On the other hand, from the uniform convergence in [0, 1] × [0, T − τ], for every h, λ small enough, there exists j_{0} such that T − 2τ < t^{j_{0}} < T − τ and u_{N+1}^{j_{0}} > K. Hence,

\begin{aligned}|T_{h,λ} − T | & ≤ & |T_{h,λ} − t^{j_{0}}| + |t^{j_{0}} − T | ≤ \frac{C}{(u_{N+1}^{j_{0}})^{γ}} + 2τ \\ & ≤ & \frac{C}{K^{γ}} + 2τ ≤ ε. \end{aligned}

This ends the proof.

Lemma 3.10. If m = 1, then there exist α > 0 and C > 0 such that

|T_{h,λ} − T | ≤ C(h^{2} + λ)^{α}.

Proof. The idea of the proof is as follows: first we get bounds for the first j_{0}, such that the error verifies

\begin{matrix}∥u_{h,λ}(·, t^{j}) − u(·, t^{j})∥_{L^{∞}([0,1])} < 1, for all j ∈ [0, j_{0}], \\ ∥u_{h,λ}(·, t^{j_{0}+1}) − u(·, t^{j_{0}+1})∥_{L^{∞}([0,1])} ≥ 1, \end{matrix}  (3.7)

then we use this bound to prove the convergence result.

To get a bound for j_{0} let us look at our scheme (2.10). From the description of the numerical scheme we have

\frac{u_{N+1}^{j+1} − u_{N+1}^{j}}{τ_{j}} = \frac{2(u_{N}^{j} − u_{N+1}^{j})}{h^{2}} + \frac{2(u_{N+1}^{j})^{p}}{h} . (3.8)

If we call v_{i}^{j} = u(x_{i}, t^{j}) we get that v satisfies

\frac{v_{N+1}^{j+1} − v_{N+1}^{j}}{τ_{j}} = \frac{2(v_{N}^{j} − v_{N+1}^{j})}{h^{2}} + \frac{2(v_{N+1}^{j})^{p}}{h} + C(h^{2} + λ)(T − t)^{θ}, (3.9)

where θ depends only on p. Let us remark that C(T − t)^{θ} is a bound for the fourth spatial derivatives of u(x, t), a regular solution of (1.1).

Subtracting (3.8) from (3.9) we have for e_{N+1}^{j} = u_{N+1}^{j} − v_{N+1}^{j},

\frac{e_{N+1}^{j+1} − e_{N+1}^{j}}{τ_{j}} ; ≤ \frac{C_{1}(e_{N}^{j} − e_{N+1}^{j})}{h^{2}} + \frac{C_{2}}{h} ξ^{p−1}|e_{N+1}^{j}| + C_{3}(h^{2} + λ)(T − t)^{θ},  (3.10)

where ξ is an intermediate value between z_{N+1}^{j} and v_{N+1}^{j}. Using that j ∈ [0, j_{0}), the bound (3.7) and the blow-up rate, we get

ξ^{p−1} ≤ C(T − t)^{− \frac{1}{2}} .

Adaptive fully discrete methods 21

Hence (3.10) gives

\frac{e_{N+1}^{j+1} − e_{N+1}^{j}}{τ_{j}} ≤ \frac{C_{1}(e_{N}^{j} − e_{N+1}^{j})}{h^{2}} + C_{2} \frac{|e_{N+1}^{j}|}{h(T − t)^{\frac{1}{2}}} + C_{3}(h^{2} + λ)(T − t)^{θ}.

This inequality is a discretization of the problem

\{ \begin{matrix}e_{t} = e_{xx} + C(T − t)^{θ}(h^{2} + λ) & (x, t) ∈ (0, 1) × [0, T), \\ e_{x}(0, t) ≤ 0 & t ∈ [0, T), \\ e_{x}(1, t) ≤ C(T − t)^{− \frac{1}{2}} e(1, t) & t ∈ [0, T), \\ e(x, 0) ≤ C(h^{2} + λ) & x ∈ (0, 1). \end{matrix}  (3.11)

Below, we will construct a supersolution for problem (3.11), such that

\overline{e}(x, t) ≤ C(h^{2} + λ)(T − t)^{−η}, (3.12)

and with the fourth spatial derivatives positives in [0, 1] and the first two time derivatives positives in [0, T).

We take e_{i}^{j} = \overline{e}(x_{i}, t^{j}), which is a supersolution. Applying a comparison argument we obtain e_{i}^{j} ≤ C(T − t^{j})^{−η}(h^{2} + λ). Similarly, applying the same ideas to ẽ_{i}^{j} = v_{i}^{j} − u_{i}^{j}, we get e_{i}^{j} ≥ −C(T − t^{j})^{−η}(h^{2} + λ). Therefore

|e_{i}^{j}| ≤ C(T − t^{j})^{−η}(h^{2} + λ).

We get a bound for j_{0} as follows: since

1 ≤ |e_{i}^{j_{0}+1}| ≤ C(T − t^{j_{0}+1})^{−η}(h^{2} + λ),

we get

(T − t^{j_{0}+1}) ≤ C(h^{2} + λ)^{\frac{1}{η}} .

With this bound for (T − t^{j_{0}+1}) we obtain the desired result. Indeed,

|T_{h,λ} − T | ≤ |T_{h,λ} − t^{j_{0}+1}| + |T − t^{j_{0}+1}| ≤ |T_{h,λ} − t^{j_{0}+1}| + C(h^{2} + λ)^{\frac{1}{η}} .

On the other hand

|T_{h,λ} − t^{j_{0}+1}| ≤ |T_{h,λ} − t^{j_{0}}| + |t^{j_{0}} − t^{j_{0}+1}| ≤ C(u_{N+1}^{j_{0}})^{−2(p−1)} + λ
≤ C(v_{N+1}^{j_{0}} + 1)^{−2(p−1)} + λ ≤ C(T − t^{j_{0}}) + λ
≤ C(h^{2} + λ)^{\frac{1}{η}} + λ ≤ C(h^{2} + λ)^{\frac{1}{η}} .

It remains to be more concrete in the construction of the supersolution to (3.11). This kind of construction was previously done in ^{14}. We include it here for the sake of completeness.

22 Brändle-Groisman-Rossi

Let us consider the following problem,

\{ \begin{matrix}E_{t} = E_{xx} + C(T − t)^{θ}(h^{2} + λ) & (x, t) ∈ (0, 1) × [0, T), \\ E_{x}(0, t) = 0 & t ∈ [0, T), \\ E_{x}(1, t) = C(T − t)^{−1/2}E(1, t) & t ∈ [0, T), \\ E(x, 0) = C(h^{2} + λ) & x ∈ (0, 1). \end{matrix}

In order to build a supersolution, \overline{E}, to this problem such that

\overline{E}(x, t) ≤ C(h^{2} + λ)(T − t)^{−η}

and with the first four spatial derivatives positives in [0, 1], we look for a supersolution of the form

\overline{E}(x, t) = C(h^{2} + λ)(T − t)^{θ}a(x, t),

with a(x, t) a solution of

\{ \begin{matrix}a_{t} = a_{xx} & (x, t) ∈ (0, 1) × [0, T), \\ a_{x}(0, t) = 0 & t ∈ [0, T), \\ a_{x}(1, t) = C(T − t)^{−1/2}a(1, t) & t ∈ [0, T), \\ a(x, 0) = a_{0}(x) & x ∈ (0, 1). \end{matrix}  (3.13)

As we want \overline{E} to have their first four spatial derivatives positives, we impose that a_{0}(x) is a smooth compatible initial datum with the first four spatial derivatives positives. The positivity of the derivatives is preserved for every t ∈ [0, T). Now let us see that there exists r such that

a(x, t) ≤ \frac{C}{(T − t)^{r}} . (3.14)

To this end we want to construct a supersolution, v(x, t), to (3.13) such that (3.14) holds. This can be easily done by the following procedure: take v(x, t) a solution of (1.1) with the boundary condition given by v_{x}(1, t) = v^{q}(1, t), with q small, and initial datum v_{0} such that v(x, t) blows up exactly at time T, see ^{24} for a proof of the fact that for every time T there exist an initial datum such that v(x, t) blows up at time T. From the known blow-up rate for solutions of (1.1), see ^{15,19}, we get that v(x, t) verifies

v(1, t) ∼ \frac{L}{(T − t)^{\frac{1}{2(q−1)}}} ,

with L^{q−1} → +∞ as q ↘ 1, see ^{15} for an explicit formula for L(q). Let us fix q such that L^{q−1}(q) > C. This choice leads to a supersolution of (3.13), since

v_{x}(1, t) = v^{q}(1, t) = v^{q−1}v(1, t) ≥ \frac{L^{q−1}}{(T − t)^{\frac{1}{2}}} v(1, t).

This ends the proof.

Adaptive fully discrete methods 23

3.5. The blow-up set.

Now we turn our interest to the blow-up set of the numerical solution. We want to look at the set of points, x, such that u_{h,λ}(x, t) → +∞ as t^{j} ↗ T_{h,λ}. Throughout this section we assume that the parameters h and λ are fixed.

We begin our analysis with the case p > m. In ^{2} and ^{13} it is proved that for a fixed mesh the numerical blow-up set is given by B(U) = [1 − Lh, 1], where L depends only on p and m. Notice that as L does not depend on h, the blow-up set converges to the continuous blow-up set B(u) = {1}, as h → 0.

We remark that the adding points method, adds at least L points near x = 1 when the solution blows up, therefore the parameter h_{i} goes to zero and formally the blow-up set for the method must be B(u_{h,λ}) = {1}. In order to prove this, we first study the propagation of blow-up.

Lemma 3.11. Assume p > m. If u_{k+1}^{j} ≤ j^{α_{k+1}}, 0 < α_{k+1} ≤ 1, then

u_{k}^{j} ≤ j^{α_{k}}, α_{k} = 1 − (2p − m − mα_{k+1}).

Proof. Recall that from Lemma 3.1 and Lemma 3.7 we have that

max_{i} u_{i}^{j} = u_{N+1}^{j} ∼ C(T_{h,λ} − t)^{− \frac{1}{2p−m−1}} ∼ j.

Therefore,

u_{k}^{j+1} − u_{k}^{j} = \frac{τ_{j}}{h^{2}} ((u_{k−1}^{j})^{m} − 2(u_{k}^{j})^{m} + (u_{k+1}^{j})^{m})

≤ \frac{Cτ_{j}}{h^{2}} (u_{k+1}^{j})^{m} ≤ \frac{Cλ}{(u_{N+1}^{j})^{2p−m}} j^{m(α_{k+1})}

≤ \frac{Cλ}{j^{2p−m}} j^{mα_{k+1}} = \frac{Cλ}{j^{2p−m−mα_{k+1}}} .

This bound implies that

u_{k}^{j+1} − u_{k}^{0} = ∑_{l=0}^{j} (u_{k}^{l+1} − u_{k}^{l}) ≤ ∑_{l=0}^{j} \frac{Cλ}{l^{2p−m−mα_{k+1}}} ,

which means that

u_{k}^{j+1} ≤ u_{k}^{0} + Cj^{−(2p−m−mα_{k+1})+1}

.

In view of this, u_{k}^{j} is bounded if α_{k} = 1 − (2p − m − mα_{k+1}) ≤ 0. Now, given a blow-up node x_{k} at most u_{k−1}^{j}, …, u_{k−L}^{j} blow up and (at least) u_{1}^{j}, …, u_{k−L−1}^{j}

24 Brändle-Groisman-Rossi

are bounded. Indeed, since u_{k}^{j} ≤ cj we can apply the previous lemma with α_{k} ≤ 1, hence

α_{k} ≤ 1
α_{k−1} ≤ 1 − (2p − 2m)
α_{k−2} ≤ 1 − (2p − 2m) − m(2p − 2m)
⋮
α_{k−l} ≤ 1 − (2p − 2m) ∑_{i=0}^{l} m^{i}.

Let L be the first integer such that (2p−2m) ∑_{i=0}^{L+1} m^{i} ≥ 1, then u_{k−L−1}^{j} is bounded.

Lemma 3.12. Assume that u_{h,λ} is given by the adding points method. If p > m then B(u_{h,λ}) = {1}.

Proof. Let x̄ < 1, we claim that the numerical solution is bounded in [0, x̄], i.e there exists C = C(x̄) > 0 such that u_{h,λ}(x, t) ≤ C, for all x ∈ [0, x̄], t ∈ [0, T_{h,λ}). Since we adapt near x = 1 collapsing (at least) L points in {1} as the solution gets large, we can consider a node x_{k+1} > x̄ such that, as j → ∞, there exits L+1 nodes between x̄ and x_{k+1}. Hence, from the previous computations, for t^{j} close to T_{h,λ}, we get that for all nodes x_{i} < x̄, u_{i}^{j} is bounded. This implies that B(u_{h,λ}) = {1}.

For the moving point method the result is different:

Lemma 3.13. Assume that u_{h,λ} is given by the moving point method. If p > m then B(u_{h,λ}) ⊆ [1 − Lh, 1], where L is given as above.

Proof. Let w be the moving node and z the numerical solution at this node. Since z^{j} ≤ u_{N+1}^{j} ≤ Cj, we have

\begin{aligned}u_{N}^{j+1} − u_{N}^{j} & = \frac{2τ_{j} ((u_{N−1}^{j})^{m}(h − h_{i}) − (u_{N}^{j})^{m}(2h − h_{i}) + (z^{j})^{m}h)}{(h − h_{i})h(2h − h_{i})} \\ & ≤ \frac{2Cj^{m−2p}(z^{j})^{m}h}{(h − h_{i})h(2h − h_{i})} ≤ Cj^{2m−2p}. \end{aligned}

Hence, we get that

u_{N}^{j+1} ≤ u_{N}^{j} + Cj^{2m−2p+1}.

This means that the node u_{N}^{j} is bounded if 2m − 2p + 1 ≤ 0. Otherwise, we apply Lemma 3.11 and get that at most the nodes w, x_{N}, …, x_{N−L+1} blow up. Therefore, we obtain B(u_{h,λ}) ⊆ [1 − Lh, 1].

Adaptive fully discrete methods 25

We still have to prove that when p ≤ m the blow-up set is given by [0, 1], i.e. there is global blow-up. Recall that in this range of parameters, the methods adapt just the time variable, letting the spacial mesh fixed.

Inspired in the self-similar variables introduced in ^{17} we define (Y^{j}) as

y_{i}^{j} = u_{i}^{j}j^{−1},  1 ≤ i ≤ N + 1.  (3.15)

In the sequel of the proof we will use Δy_{i}^{j+1} to denote

\frac{y_{i}^{j+1} − y_{i}^{j}}{τ_{j}/j^{−(p−1)}} ,

This can be thought as τ_{j}/j^{−(p−1)} to be the time step in the new variables. With this notation Y^{j} verifies

\{ \begin{matrix}Δy_{1}^{j+1} = & −\frac{2j^{m−(p−1)} ((y_{2}^{j})^{m} − (y_{1}^{j})^{m})}{(j + 1)h^{2}} − \frac{j^{−(p−1)}u_{i}^{j}}{τ_{j}j(j + 1)}, \\ Δy_{i}^{j+1} = & −\frac{j^{m−(p−1)} ((y_{i−1}^{j})^{m} − 2(y_{i}^{j})^{m} + (y_{i+1}^{j})^{m})}{(j + 1)h^{2}} \\ & −\frac{j^{−(p−1)}u_{i}^{j}}{τ_{j}j(j + 1)},  2 ≤ i ≤ N, \\ Δy_{N+1}^{j+1} = & −\frac{2j^{m−(p−1)} ((y_{N}^{j})^{m} − (y_{N+1}^{j})^{m})}{(j + 1)h^{2}} \\ & +\frac{2j^{p−(p−1)}(y_{N+1}^{j})^{p}}{(j + 1)h} − \frac{j^{−(p−1)}u_{i}^{j}}{τ_{j}j(j + 1)}, \end{matrix}  (3.16)

and the initial datum becomes

y_{i}^{0} = j^{−1}u_{0}(x_{i}), 1 ≤ i ≤ N + 1.

From the blow-up rate proved in Lemma 3.8, we have that the vector Y^{j} is bounded, i.e

Y^{j} ≤ C, for all j.

Moreover, there exist positive constants, c and C such that the last node verifies

c ≤ y_{N+1}^{j} ≤ C, for all j.

Lemma 3.14. If 1 < p < m then

u_{k}^{j} ∼ \frac{C}{(T_{h,λ} − t)^{\frac{1}{p−1}}} ∼ Cj, 1 ≤ k ≤ N,

therefore B(u_{h,λ}) = [0, 1]. This phenomena is known as uniform global blow-up.

26 Brändle-Groisman-Rossi

Proof. From the equation for the last node in (3.16), the expression for τ_{j} and the asymptotic behavior we get

y_{N+1}^{j+1} − y_{N+1}^{j} ≥ \frac{a_{j}}{j^{1+p−m}} − \frac{c}{j} , (3.17)

where a_{j} = (y_{N+1}^{j})^{m} − (y_{N}^{j})^{m} and c is a positive constant. If we write y_{N+1}^{j+1} as

y_{N+1}^{j+1} = ∑_{l=0}^{j} (y_{N+1}^{l+1} − y_{N+1}^{l}) + y_{N+1}^{0},

from (3.17), we get

y_{N+1}^{j+1} ≥ ∑_{l=0}^{j} (\frac{a_{l}}{l^{1+p−m}} − \frac{c}{l}) + y_{N+1}^{0}.

On the other hand, since y_{N+1}^{j} is bounded, we have

(\frac{a_{j}}{j^{1+p−m}} − \frac{c}{j} ) j → 0,

which implies that

\frac{a_{j}}{j^{p−m}} → c.

Hence, since p < m, we get

0 ≤ (y_{N+1}^{j})^{m} − (y_{N}^{j})^{m} = a_{j} → 0, as j → ∞.

This means that y_{N}^{j} ≥ c. Applying the same argument with y_{k}^{j} we obtain that y_{k−1}^{j} ≥ c for all 2 ≤ k ≤ N, which in terms of the original variables leads to

u_{k}^{j} ≥ c(T_{h,λ} − t^{j})^{− \frac{1}{p−1}} for 1 ≤ k ≤ N.

The reverse inequality is obtained trivially from the fact that u_{k}^{j} ≤ u_{N+1}^{j} ≤ c(T_{h,λ} − t^{j})^{−1/(p−1)}. We obtain uniform global blow-up.

Lemma 3.15. If p = m then

u_{k}^{j} ∼ \frac{C}{(T_{h,λ} − t)^{\frac{1}{p−1}}} ∼ Cj, 1 ≤ k ≤ N,

therefore B(u_{h,λ}) = [0, 1].

Proof. As we have already seen u_{N+1}^{j} ≥ cj. Let us prove that from this estimate we get that all nodes behave as the last one, i.e. u_{k}^{j} ≥ cj, for all 1 ≤ k ≤ N, and hence the blow-up is global.

First of all observe that U^{j} is a supersolution of

M(U^{j+1} − U^{j}) ≥ −τ_{j}A(U^{j})^{m} + τ_{j}Bcj^{m}.

Adaptive fully discrete methods 27

Let us look for a subsolution of this problem of the form V^{j} = jZ, where Z = (z_{1}, …, z_{N+1}), and does not depend on j. Then Z must verify

CMZ ≤ −AZ^{m} + Bc,

which can be written explicitly as

\{ \begin{matrix}Cz_{1} & ≤ \frac{2(z_{2}^{m} − z_{1}^{m})}{h^{2}}, \\ Cz_{i} & ≤ \frac{z_{i−1}^{m} − 2z_{i}^{m} + z_{i+1}^{m}}{h^{2}} , 2 ≤ i ≤ N \\ Cz_{N+1} & ≤ \frac{2(z_{N}^{m} − z_{N+1}^{m})}{h^{2}} + c. \end{matrix}  (3.18)

In order to find Z > 0 verifying the above system, we proceed as follows: we fix z_{1} > 0 which will be determined later and compute z_{2}, …, z_{N+1} in terms of z_{1}. We choose z_{2} in such a way that the first inequality in (3.18) is verified with an equality. Iterating this procedure we obtain z_{3}, …, z_{N+1} using the first N equations. It remains to check that the last inequality holds. Observe that z_{N+1} is a continuous positive function of z_{1} and that z_{N+1} = 0 when z_{1} = 0. Therefore the last inequality is verified, if z_{1} is small enough. Hence, as the spatial mesh is fixed, we may use a standard comparison argument, which is proved in a similar way as Lemma 3.3, and get v_{i}^{j} = z_{i}j ≤ u_{i}^{j}, 1 ≤ i ≤ N.

Remark that, as in the case p < m, the upper bound u_{i}^{j} ≤ cj follows easily from the behaviour of the last node. Therefore, we obtain uniform complete blow-up.

References

1. L. M. Abia, J. C. Lopez-Marcos, J. Martinez. Blow-up for semidiscretizations of reaction diffusion equations. Appl. Numer. Math. Vol. 20, (1996), 145-156.

2. G. Acosta, R. Duran and J. D. Rossi. An adaptive time step procedure for a parabolic problem with blow-up. Computing. Vol. 68 (4), (2002), 343-373.

3. G. Acosta, J. Fernandez Bonder, P. Groisman and J.D. Rossi. Numerical approximation of a parabolic problem with a nonlinear boundary condition in several space dimensions. Discr. Cont. Dyn. Sys. B, Vol. 2 (2), (2002), 279-294.

4. C. Bandle and H. Brunner. Blow-up in diffusion equations: a survey. J. Comp. Appl. Math. Vol. 97, (1998), 3-22.

5. M. Berger and R. V. Kohn. A rescaling algorithm for the numerical calculation of blowing up solutions. Comm. Pure Appl. Math. Vol. 41, (1988), 841-863.

6. C. Budd and G. Collins. An invariant moving mesh scheme for the nonlinear diffusion equation. Appl. Numer. Math. Vol. 26, (1998), 23-39.

7. C. Budd, S. Chen and R. D. Russell. New self-similar solutions of the nonlinear Schrodinger equation with moving mesh computations. J. Comp. Phys. Vol. 152, (1999), 756-789.

8. C. J. Budd, W. Huang and R. D. Russell. Moving mesh methods for problems with blow-up. SIAM Jour. Sci. Comput. Vol. 17(2), (1996), 305-327.

9. P. Ciarlet. The finite element method for elliptic problems. North Holland,

28 Brändle-Groisman-Rossi

10. C. Cortázar, M. Elgueta and O. Venegas. On the Blow-up set for u_{t} = (u^{m})_{xx}, m > 1, with nonlinear boundary conditions. Preprint.

11. R. G. Duran, J. I. Etcheverry and J. D. Rossi. Numerical approximation of a parabolic problem with a nonlinear boundary condition. Discr. and Cont. Dyn. Sys. Vol. 4 (3), (1998), 497-506.

12. J. Fernández Bonder and J. D. Rossi. Blow-up vs. spurious steady solutions. Proc. Amer. Math. Soc. Vol. 129 (1), (2001), 139-144.

13. R. Ferreira, P. Groisman and J. D. Rossi. Numerical blow-up for a nonlinear problem with a nonlinear boundary condition. Math. Models Methods Appl. Sci. M3AS. Vol. 12 (4), (2002), 461-484.

14. R. Ferreira, P. Groisman and J.D. Rossi. Adaptive numerical schemes for a parabolic problem with blow-up. IMA J. Numer. Anal. Vol. 23 (3), (2003), 439– 463.

15. M. Fila and P. Quittner. The blowup rate for the heat equation with a nonlinear boundary condition. Math. Methods Appl. Sci., Vol. 14, (1991), 197–205.

16. J. Filo. Diffusivity versus absorption through the boundary. J. Differ. Equations, Vol. 99, (1992), n. 2, 281–305.

17. Y. Giga and R. V. Kohn. Nondegeneracy of blow up for semilinear heat equations. Comm. Pure Appl. Math. Vol. 42, (1989), 845-884.

18. W. Huang, Y. Ren and R. D. Russell. Moving mesh partial differential equations (MMPDEs) based on the equidistribution principle. SIAM J. Numer. Anal. Vol 31, (1994), 709-730.

19. B. Hu and H. M. Yin. The profile near blowup time for solution of the heat equation with a nonlinear boundary condition. Trans. Amer. Math. Soc., Vol. 346(1), (1994), 117–135.

20. M. N. Le Roux. Semidiscretizations in time of nonlinear parabolic equations with blow-up of the solutions. SIAM J. Numer. Anal. Vol. 31, (1994), 170-195.

21. T. Nakagawa. Blowing up of a finite difference solution to u_{t} = u_{xx} + u^{2}. Appl. Math. Optim., 2, (1976), 337-350.

22. T. Nakagawa and T. Ushijima. Finite element analysis of the semilinear heat equation of blow-up type. Topics in Numerical Analysis III (J.J.H Miller, ed.), Academic Press, London, (1977), 275-291.

23. C. V. Pao, “Nonlinear parabolic and elliptic equations”. Plenum Press, New York, 1992.

24. J. D. Rossi. The blow-up rate for a system of heat equations with nontrivial coupling at the boundary. Math. Methods Appl. Sci., vol 20, (1997), 1-11.